Kidney renal clear cell carcinoma (KIRC).Rmd

---
title: "Systems and Precision medicine - Kidney renal clear cell carcinoma (KIRC); "
author: "André Simão (68816), Maria Inês Gomes (68828), Mariana Vasques (68829)"
date: "2024-05-23"
output:
  pdf_document:
    number_sections: true
---

# Kidney renal clear cell carcinoma (KIRC)

```{r setup, include=FALSE}
KIRC<- load("C:/Users/inesv/OneDrive - FCT NOVA/Nano Project 1 - 3/KIRC_geneExpression_allGenes.RData")
library(dplyr)
library(ggplot2)
library(glmnet)
library(pROC)
library(gmodels)
library(caret)
library(e1071)
```

```{r}
dim(clinicalInfo)
```

Inspect how many normal and tumor samples we have:
```{r, echo=FALSE}
# Extract sample IDs
sample_ids <- colnames(readCounts)

tumor_samples <- grep("0[1-9]$", sample_ids, value = TRUE)
normal_samples <- grep("1[0-9]$", sample_ids, value = TRUE)
control_samples <- grep("2[0-9]$", sample_ids, value = TRUE)


cat("Normal samples:\n")
cat(length(normal_samples), "\n")
cat("Tumor samples:\n")
cat(length(tumor_samples), "\n")
cat("Control samples:\n")
cat(length(control_samples), "\n")
```


```{r}
clinical_df <- data.frame(clinicalInfo)
genes_df <- data.frame(readCounts)
# head(clinical_df)
# head(genes_df)
```


```{r}
colnames(genes_df) <- gsub("\\.\\d+$", "", colnames(genes_df))
```

```{r}
dim(genes_df)
dim(clinical_df)
```

RNA-seq data from 164 patients measured over 17979 genes

```{r}
clinical_df %>% 
  count(clinical_df$pathology_M_stage, sort = TRUE)
```


```{r}
patients_to_remove <- rownames(clinical_df)[clinical_df$pathology_M_stage == "mx" | is.na(clinical_df$pathology_M_stage)]
clinical_filtered <- clinical_df[!rownames(clinical_df) %in% patients_to_remove, ]
pattern_to_remove <- paste(patients_to_remove, collapse = "|")
columns_to_remove <- grep(pattern_to_remove, colnames(genes_df), value = TRUE)

genes_filtered <- genes_df[, !colnames(genes_df) %in% columns_to_remove]

genes_filtered <- t(genes_filtered)
```


```{r}
clinical_filtered %>% 
  count(clinical_filtered$pathology_M_stage, sort = TRUE)
```

```{r}
# Get the common patient IDs
common_patient_ids <- intersect(rownames(genes_filtered), rownames(clinical_filtered))

# Subset genes_filtered and clinical_filtered based on the common patient IDs
genes_filtered <- genes_filtered[common_patient_ids, , drop = FALSE]
clinical_filtered <- clinical_filtered[common_patient_ids, , drop = FALSE]
```

```{r}
dim(clinical_filtered)
dim(genes_filtered)
```
```{r}
Xdata = genes_filtered
Ydata <- clinical_filtered$pathology_M_stage

dim(Xdata)
```
Now we have RNA-seq data from 155 patients measured over 17979 genes


```{r}
Xdata[1:5,1:8]
Ydata[1:8]
length(Ydata)
summary(as.factor(Ydata))
```
## Data preprocessing

**Data filtering** by removing samples with standard deviation zero.
First, it calculates the standard deviation of each gene expression column. Then remove the columns (genes) with zero variance.
```{r}
Xdata_sd<-sapply(seq(ncol(Xdata)), function(ix){sd(Xdata[,ix])})
Xdata<-Xdata[,Xdata_sd!=0]
dim(Xdata)
```
Now we have RNA-seq data from 155 patients measured over 17661 genes

**Data normalization** by computing the z-score to standardize the gene expression data to have mean 0 and standard deviation 1.
```{r}
Xdata_sc<-scale(Xdata)
```

# Exploratory data analysis
## Principal component analysis (PCA)

First, it will perform a PCA on the normalized data (scale = TRUE).
PCA transforms the data into a new coordinate system where the greatest variances by any projection of the data come to lie on the first few principal components.
Then, it will extract the PCA results, using the function: metastasis_pca$x, that extracts the PC scores from the PCA object, which are the transformed coordinates of the original data in the new PC space.
Conversion: The scores are converted into a data frame for easier manipulation and plotting.
Adding Class Labels: The original class labels (Ydata) are added to the scores data frame. This allows us to color-code the points in plots according to their tumor type.

**variance_explained**: Squaring the standard deviations gives the variance explained by each principal component.
**prop_variance_explained**: Proportion of the total variance explained by each principal component.
**cumulative_prop_variance**: Cumulative proportion of variance explained, summing up the proportions progressively to show how much total variance is explained as we add more principal components.
```{r}
metastasis_pca <- prcomp(Xdata, scale = TRUE) 

# Scores of the principal components
scores <- as.data.frame(metastasis_pca$x)
scores$class <- Ydata

# Variance explained by each component
variance_explained <- metastasis_pca$sdev^2
prop_variance_explained <- variance_explained / sum(variance_explained)
cumulative_prop_variance <- cumsum(prop_variance_explained)
```

Plot PCA scores plot:

```{r}
pc1_variance <- round(prop_variance_explained[1] * 100, 2)
pc2_variance <- round(prop_variance_explained[2] * 100, 2)

# Create the plot with updated axis labels
ggplot(scores, aes(x = PC1, y = PC2, color = class)) +
  geom_point() +
  labs(
    title = "PCA Scores plot",
    x = paste("PC1 (", pc1_variance, "%)", sep = ""),
    y = paste("PC2 (", pc2_variance, "%)", sep = "")
  ) +
  theme_minimal()
```

For the **PCA: Cumulative Variance Explained** a data frame pca_summary is created with two columns:
* PC: The principal component numbers.
* Cumulative_Prop_Variance: The cumulative proportion of variance explained by the principal components.
```{r}
# Plot % of variance explained
pca_summary <- data.frame(
  PC = 1:length(cumulative_prop_variance),
  Cumulative_Prop_Variance = cumulative_prop_variance
)

ggplot(pca_summary, aes(x = PC, y = Cumulative_Prop_Variance)) +
  geom_line(color = "red") +
  geom_point(color = "gray") +
  labs(title = "PCA: Cumulative Variance Explained",
       x = "Principal Component (PC)",
       y = "Cumulative Proportion of Variance Explained") +
  theme_minimal()
```
The **PCA** reduces the dimensionality of the data by transforming it into principal components that capture the most variance. The **Scores** represent the data in the new principal component space. The **Variance Explained** helps understand how much of the total data variability is captured by each principal component. Finally, plotting the **cumulative variance explained** helps visualize how many principal components are needed to capture most of the variance in the data, which is essential for understanding the underlying structure and for dimensionality reduction.


# Regularized logistic regression
## Model training

**Partition the data into training and test sets**
A vector test_ID is created containing the indices of the samples that will be included in the test set. First, creates a sequence from 1 to the number of rows in Xdata_sc (number of samples = 91). Then, it calculates 25% of the total number of samples to determine the size of the test set (22.75 =~ 23). replace=FALSE, ensures that each sample is chosen only once.

Then, it converts the class labels from categorical strings to numeric values: m1 = 1 and m0 = 0.

Training Set:

Xdata_train: Selects all rows from Xdata_sc except those indexed by test_ID and converts the result to a matrix.
Ydata_train: Selects all elements from Ydata except those indexed by test_ID and converts the result to a factor.
Test Set:

Xdata_test: Selects the rows indexed by test_ID from Xdata_sc and converts the result to a matrix.
Ydata_test: Selects the elements indexed by test_ID from Ydata and converts the result to a factor.
```{r}
set.seed(11) # for reproducibility
test_ID <- sample(1:dim(Xdata_sc)[1], round(dim(Xdata_sc)[1]*0.25), replace=FALSE)

Ydata[Ydata=="m1"] <- 1
Ydata[Ydata=="m0"] <- 0

# train set
Xdata_train <- as.matrix(Xdata_sc[-test_ID,])
Ydata_train <- as.factor(Ydata[-test_ID])
# test set
Xdata_test <- as.matrix(Xdata_sc[test_ID,])
Ydata_test <- as.factor(Ydata[test_ID])
```

**Building the sparse logistic regression model (lambda optimized by cross-validation)**
Fitting the Model using glmnet and the training data matrix (Xdata_train), and the training response variable (class labels) (Ydata_train).

First, we specify that we are performing logistic regression (binomial) and using 10-fold cross-validation to find the optimal regularization parameter (lambda).
Alpha specifies the type of regularization:
alpha=1 indicates Lasso regression (L1 regularization), which tends to produce sparse models by setting some coefficients to zero.
Finally, uses the Area Under the ROC Curve (AUC) as the performance measure during cross-validation.
```{r}
set.seed(10) # for reproducibility
# alpha parameter gives the balance between ...
metastasis_fit <- cv.glmnet(Xdata_train, Ydata_train, family="binomial", nfolds=10, alpha=1, type.measure="auc")

metastasis_coef <- metastasis_fit$glmnet.fit$beta[,which(metastasis_fit$cvm == max(metastasis_fit$cvm))]

genes_selected <- which(metastasis_coef != 0)
length(genes_selected)
```

```{r}
genes_selected
```


```{r}
data <- data.frame(
  x = 1:length(genes_selected),
  y = metastasis_coef[genes_selected],
  labels = names(metastasis_coef[genes_selected])
)

ggplot(data, aes(x = x, y = y, label = labels)) +
  geom_point(color = "dodgerblue1") +
  geom_text(hjust = 0, vjust = 0) +
  labs(title = "Selected variable's coefficients",
       x = "Variable",
       y = "Coefficients") +
  theme_minimal()
```

# Model predictive performance
## Predicting for the training set
```{r}
pred_train <- predict(metastasis_fit, Xdata_train, type="class", type.measure = "auc", s = "lambda.min")
# Confusion matrix for the train set
conf_matrix_train <- table(Ydata_train,pred_train)
conf_matrix_train
# 127 classified good LGG-od
# Calculate AUC for the train set
roc_obj <- roc(as.numeric(as.character(Ydata_train)), as.numeric(pred_train))
auc(roc_obj)

# Predicting for a test set
pred_test <- predict(metastasis_fit,Xdata_test,type="class")
# Confusion matrix for the test set
conf_matrix_test <- table(Ydata_test,pred_test)
conf_matrix_test
# Calculate AUC value for the test set
roc_obj <- roc(as.numeric(as.character(Ydata_test)), as.numeric(pred_test))
auc(roc_obj)
```

```{r, echo=FALSE}
# Create ROC curves for both training and test sets
roc_train <- roc(as.numeric(as.character(Ydata_train)), as.numeric(pred_train))
roc_test <- roc(as.numeric(as.character(Ydata_test)), as.numeric(pred_test))

# Plot ROC curves for both training and test sets
plot(roc_train, main = "ROC Curves - Training vs. Test Set", col = "lightskyblue1", lwd = 2, lty = 1,
     xlab = "False Positive Rate", ylab = "True Positive Rate")
lines(roc_test, col = "tomato1", lwd = 2, lty = 1)
legend("bottomright", legend = c(paste("Training AUC =", round(auc(roc_train), 2)),
                                 paste("Test AUC =", round(auc(roc_test), 2))),
       col = c("lightskyblue1", "tomato1"), lty = 1, cex = 0.8)
```
With statistics:

```{r}
# Calculate performance metrics for the train set
confusionMatrix(data = factor(pred_train), reference = factor(Ydata_train))
# Calculate performance metrics for the test set
confusionMatrix(data = factor(pred_test), reference = factor(Ydata_test))
```
More details conf matrix:
```{r}
#Computes the crosstable calculations
CrossTable(Ydata_test,pred_test)
```
Predicted
  TN |  FP | 
-----|-----|
  FN |  TP | 
  
Calculate the performance metrics like accuracy, precision, recall, F1
```{r}
recall <- (15/(15+6))
```


# Regularized logistic regression (DEGs)
## Select DEGs
```{r}
setwd("C:/Users/inesv/OneDrive - FCT NOVA/Nano Project 1 - 3")
library(readr)
DEG <- read_csv("DEG.csv")
DEG <- column_to_rownames(DEG, var = "ID")

Xdata <- genes_filtered
Ydata <- clinical_filtered$pathology_M_stage
dim(Xdata)
length(Ydata)
summary(as.factor(Ydata))


DEGs_genes <- rownames(DEG)
Xdata_degs <- Xdata[, colnames(Xdata) 
                    %in% DEGs_genes]
dim(Xdata_degs)
```


```{r}
Xdata_degs_sc<-scale(Xdata_degs)
dim(Xdata_degs_sc)
```

## Model training

```{r}
set.seed(9) # for reproducibility
test_ID <- sample(1:dim(Xdata_degs_sc)[1], round(dim(Xdata_degs_sc)[1]*0.25), replace=FALSE)

Ydata[Ydata=="m1"] <- 1
Ydata[Ydata=="m0"] <- 0

# train set
Xdata_train <- as.matrix(Xdata_degs_sc[-test_ID,])
Ydata_train <- as.factor(Ydata[-test_ID])
# test set
Xdata_test <- as.matrix(Xdata_degs_sc[test_ID,])
Ydata_test <- as.factor(Ydata[test_ID])
```

```{r}
set.seed(0) # for reproducibility
# alpha parameter gives the balance between ...
metastasis_fit <- cv.glmnet(Xdata_train, Ydata_train, family="binomial", nfolds=10, alpha=1, type.measure="auc")

metastasis_coef <- metastasis_fit$glmnet.fit$beta[,which(metastasis_fit$cvm == max(metastasis_fit$cvm))]

genes_selected <- which(metastasis_coef != 0)
length(genes_selected)
```




```{r}
data <- data.frame(
  x = 1:length(genes_selected),
  y = metastasis_coef[genes_selected],
  labels = names(metastasis_coef[genes_selected])
)

ggplot(data, aes(x = x, y = y, label = labels)) +
  geom_point(color = "dodgerblue1") +
  geom_text(hjust = 0, vjust = 0) +
  labs(title = "Selected variable's coefficients",
       x = "Variable",
       y = "Coefficients") +
  theme_minimal()
```


```{r}
genes_selected
```

## Model predictive performance

```{r}
pred_train <- predict(metastasis_fit, Xdata_train, type="class", type.measure = "auc", s = "lambda.min")
# Confusion matrix for the train set
conf_matrix_train <- table(Ydata_train,pred_train)
conf_matrix_train
# 127 classified good LGG-od
# Calculate AUC for the train set
roc_obj <- roc(as.numeric(as.character(Ydata_train)), as.numeric(pred_train))
auc(roc_obj)

# Predicting for a test set
pred_test <- predict(metastasis_fit,Xdata_test,type="class")
# Confusion matrix for the test set
conf_matrix_test <- table(Ydata_test,pred_test)
conf_matrix_test
# Calculate AUC value for the test set
roc_obj <- roc(as.numeric(as.character(Ydata_test)), as.numeric(pred_test))
auc(roc_obj)
```


```{r, echo=FALSE}
# Create ROC curves for both training and test sets
roc_train <- roc(as.numeric(as.character(Ydata_train)), as.numeric(pred_train))
roc_test <- roc(as.numeric(as.character(Ydata_test)), as.numeric(pred_test))

# Plot ROC curves for both training and test sets
plot(roc_train, main = "ROC Curves - Training vs. Test Set", col = "lightskyblue1", lwd = 2, lty = 1,
     xlab = "False Positive Rate", ylab = "True Positive Rate")
lines(roc_test, col = "tomato1", lwd = 2, lty = 1)
legend("bottomright", legend = c(paste("Training AUC =", round(auc(roc_train), 2)),
                                 paste("Test AUC =", round(auc(roc_test), 2))),
       col = c("lightskyblue1", "tomato1"), lty = 1, cex = 0.8)
```


```{r}
# Calculate performance metrics for the train set
confusionMatrix(data = factor(pred_train), reference = factor(Ydata_train))
# Calculate performance metrics for the test set
confusionMatrix(data = factor(pred_test), reference = factor(Ydata_test))
```