---
title: "Tracking intratumoral heterogeneity in glioblastoma via regularized classification of single-cell RNA-Seq data"
author: "Marta B. Lopes and Susana Vinga"
date: "January, 2020"
output:
  html_document:
    number_sections: yes
    toc: yes
  pdf_document:
    number_sections: yes"N"
    toc: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# GBM dataset

```{r}
# Data available at 'http://gbmseq.org/' 
GBM_metadata <- read.table('GBM_metadata.csv', sep = "")
summary(GBM_metadata)

# tumor cells (# 2343)
summary(GBM_metadata[which(GBM_metadata$Location == "Tumor"),2])

tumor_id <- rownames(GBM_metadata[which(GBM_metadata$Location == "Tumor"),]) 
length(tumor_id) # 2343

# 1091 cells in the 3 neoplastic clusters; only 1029 are neoplastic cells from the tumor core; 62 are peripheral neoplastic cells that migrated from the tumor core to the surrounding peritumural space, i.e., infiltrating cells. 

tumor_core_c1_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "1"),])
tumor_core_c4_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "4"),])
tumor_core_c11_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "11"),])
tumor_core_id <- c(tumor_core_c1_id,tumor_core_c4_id,tumor_core_c11_id)
length(tumor_core_id) # 1091

tumor_core_neoplastic_c1_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "1" & GBM_metadata$Location == "Tumor"),])
tumor_core_neoplastic_c4_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "4" & GBM_metadata$Location == "Tumor"),])
tumor_core_neoplastic_c11_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "11" & GBM_metadata$Location == "Tumor"),])
tumor_core_neoplastic_id <- c(tumor_core_neoplastic_c1_id,tumor_core_neoplastic_c4_id,tumor_core_neoplastic_c11_id)
length(tumor_core_neoplastic_id) # 1029

GBM_tumor_core_neoplastic <- GBM_metadata[tumor_core_neoplastic_id,]
dim(GBM_tumor_core_neoplastic)

tumor_periphery_neoplastic_c1_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "1" & GBM_metadata$Location == "Periphery"),])
tumor_periphery_neoplastic_c4_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "4" & GBM_metadata$Location == "Periphery"),])
tumor_periphery_neoplastic_c11_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "11" & GBM_metadata$Location == "Periphery"),])
tumor_periphery_neoplastic_id <- c(tumor_periphery_neoplastic_c1_id,tumor_periphery_neoplastic_c4_id,tumor_periphery_neoplastic_c11_id)
length(tumor_periphery_neoplastic_id) # 62

GBM_tumor_periphery_neoplastic <- GBM_metadata[tumor_periphery_neoplastic_id,]
dim(GBM_tumor_periphery_neoplastic)
  
tumor_distant_c1_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "1" & GBM_metadata$Location == "Distant"),])
tumor_distant_c4_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "4" & GBM_metadata$Location == "Distant"),])
tumor_distant_c11_id <- rownames(GBM_metadata[which(GBM_metadata$Cluster_2d == "11" & GBM_metadata$Location == "Distant"),])
tumor_distant_id <- c(tumor_distant_c1_id,tumor_distant_c4_id,tumor_distant_c11_id)
length(tumor_distant_id) # 0

GBM_tumor_distant <- GBM_metadata[tumor_distant_id,]
dim(GBM_tumor_distant)

distant_id <- rownames(GBM_metadata[which(GBM_metadata$Location == "Distant"),]) 
length(distant_id) # 57 distant
GBM_distant <- GBM_metadata[distant_id,]

periphery_id <- rownames(GBM_metadata[which(GBM_metadata$Location == "Periphery"),]) 
length(periphery_id) # 1189 periphery
# from which # 62 are in the neoplastic clusters
length(which(tumor_periphery_neoplastic_id %in% periphery_id))
GBM_periphery <- GBM_metadata[periphery_id,]
GBM_periphery_nonneoplastic <- GBM_periphery[-which(rownames(GBM_periphery) %in% tumor_periphery_neoplastic_id),]
dim(GBM_periphery_nonneoplastic)

# the cell types per cluster
summary(GBM_metadata[which(GBM_metadata$Cluster_2d == "1"),2])
summary(GBM_metadata[which(GBM_metadata$Cluster_2d == "4"),2])
summary(GBM_metadata[which(GBM_metadata$Cluster_2d == "11"),2])

# the cell location origin per cluster
summary(GBM_metadata[which(GBM_metadata$Cluster_2d == "1"),3])
summary(GBM_metadata[which(GBM_metadata$Cluster_2d == "4"),3])
summary(GBM_metadata[which(GBM_metadata$Cluster_2d == "11"),3])

# peritumoral cells (# 1246)
summary(GBM_metadata[which(GBM_metadata$Location == "Distant"),2])
summary(GBM_metadata[which(GBM_metadata$Location == "Periphery"),2])

summary(GBM_periphery_nonneoplastic)

summary(GBM_metadata[which(GBM_metadata$Location=="Distant"),2])
 
summary(GBM_metadata[which(GBM_metadata$Location=="Periphery"),2])
 
summary(GBM_metadata[which(GBM_metadata$Location=="Tumor"),2])
 
sum(summary(GBM_metadata[which(GBM_metadata$Location=="Distant"),2]))

sum(summary(GBM_metadata[which(GBM_metadata$Location=="Periphery"),2]))

sum(summary(GBM_metadata[which(GBM_metadata$Location=="Tumor"),2]))
                 
table(GBM_metadata[which(GBM_metadata$Cluster_2d == "1"),c(2,3)])
table(GBM_metadata[which(GBM_metadata$Cluster_2d == "4"),c(2,3)])
table(GBM_metadata[which(GBM_metadata$Cluster_2d == "11"),c(2,3)])

apply(rbind(table(GBM_metadata[which(GBM_metadata$Cluster_2d == "1" & GBM_metadata$Location == "Tumor"),2]),table(GBM_metadata[which(GBM_metadata$Cluster_2d == "4" & GBM_metadata$Location == "Tumor"),2]), table(GBM_metadata[which(GBM_metadata$Cluster_2d == "11" & GBM_metadata$Location == "Tumor"),2])),2,sum)

apply(rbind(table(GBM_metadata[which(GBM_metadata$Cluster_2d == "1" & GBM_metadata$Location == "Periphery"),2]),table(GBM_metadata[which(GBM_metadata$Cluster_2d == "4" & GBM_metadata$Location == "Periphery"),2]), table(GBM_metadata[which(GBM_metadata$Cluster_2d == "11" & GBM_metadata$Location == "Periphery"),2])),2,sum)

apply(rbind(table(GBM_metadata[which(GBM_metadata$Cluster_2d == "1" & GBM_metadata$Location == "Distant"),2]),table(GBM_metadata[which(GBM_metadata$Cluster_2d == "4" & GBM_metadata$Location == "Distant"),2]), table(GBM_metadata[which(GBM_metadata$Cluster_2d == "11" & GBM_metadata$Location == "Distant"),2])),2,sum)

GBM_tumor_core_neoplastic_astrocytes_id <- rownames(GBM_tumor_core_neoplastic[which(GBM_tumor_core_neoplastic$Selection == "Astrocytes(HEPACAM)"),])
length(GBM_tumor_core_neoplastic_astrocytes_id)

GBM_tumor_periphery_neoplastic_astrocytes_id <- rownames(GBM_tumor_periphery_neoplastic[which(GBM_tumor_periphery_neoplastic$Selection == "Astrocytes(HEPACAM)"),])
length(GBM_tumor_periphery_neoplastic_astrocytes_id)

GBM_tumor_periphery_astrocytes_id <- rownames(GBM_periphery_nonneoplastic[which(GBM_periphery_nonneoplastic$Selection == "Astrocytes(HEPACAM)"),])
length(GBM_tumor_periphery_astrocytes_id)

# Gene counts

GBM_normalized_gene_counts <- t(read.table('GBM_normalized_gene_counts.csv', sep = ""))
dim(GBM_normalized_gene_counts)

rownames(GBM_normalized_gene_counts) <- sub("X", "", rownames(GBM_normalized_gene_counts))

GBM_normalized_gene_counts_tsne_cell_type <- GBM_normalized_gene_counts[c(GBM_tumor_core_neoplastic_astrocytes_id,GBM_tumor_periphery_neoplastic_astrocytes_id,GBM_tumor_periphery_astrocytes_id),]

GBM_cell_type <- as.factor(c(rep("T_core",length(GBM_tumor_core_neoplastic_astrocytes_id)),rep("T_periphery",length(GBM_tumor_periphery_neoplastic_astrocytes_id)),rep("N_periphery",length(GBM_tumor_periphery_astrocytes_id))))

   # T-SNE plot

library("Rtsne")
tsne_color <- c(rep("T-core",length(GBM_tumor_core_neoplastic_astrocytes_id)),rep("T-periphery",length(GBM_tumor_periphery_neoplastic_astrocytes_id)),rep("N-core",length(GBM_tumor_periphery_astrocytes_id)))

set.seed(39)
tsne <- Rtsne(GBM_normalized_gene_counts_tsne_cell_type, perplexity=30,verbose=FALSE, max_iter = 500)

tsne_GBM_data <- as.data.frame(tsne$Y)

tsne_GBM_tumor_core_neoplastic_astrocytes_id = tsne_GBM_data[1:length(GBM_tumor_core_neoplastic_astrocytes_id),]  
tsne_GBM_tumor_periphery_neoplastic_astrocytes_id = tsne_GBM_data[(length(GBM_tumor_core_neoplastic_astrocytes_id)+1):(length(GBM_tumor_core_neoplastic_astrocytes_id)+length(GBM_tumor_periphery_neoplastic_astrocytes_id)),] 
tsne_GBM_tumor_periphery_astrocytes_id = tsne_GBM_data[(length(GBM_tumor_core_neoplastic_astrocytes_id)+length(GBM_tumor_periphery_neoplastic_astrocytes_id)+1):(length(GBM_tumor_core_neoplastic_astrocytes_id)+length(GBM_tumor_periphery_neoplastic_astrocytes_id)+length(GBM_tumor_periphery_astrocytes_id)),] 

library("ggplot2")
ggplot(tsne_GBM_tumor_core_neoplastic_astrocytes_id, aes(x=V1, y=V2)) +  
  geom_point(color = "#666666",size=3) +
  xlab("") + ylab("") +
  geom_point(colour = "#CC6699", size = 1.2) +
  theme(axis.text.x=element_blank(),legend.title = element_text(size=12, color="white"),legend.text =    element_text(size=14),axis.text.y=element_blank()) +
  geom_point(data=tsne_GBM_tumor_periphery_neoplastic_astrocytes_id,aes(x=V1, y=V2), col="#FF9966",size=3) + 
  geom_point(data=tsne_GBM_tumor_periphery_neoplastic_astrocytes_id,aes(x=V1, y=V2),colour = "#CC6699", size = 1.2) +
  geom_point(data=tsne_GBM_tumor_periphery_astrocytes_id,aes(x=V1, y=V2), col="#FF9966",size=4) +
  geom_point(data=tsne_GBM_tumor_periphery_astrocytes_id,aes(x=V1, y=V2),colour = "white", size = 1.2) +
  theme_minimal()
```

# Classification of GBM cells

## Regularized classification

### Discrimination between neoplastic astrocytes from the periphery and tumor core

```{r}
# Building the dataset
xdata <- GBM_normalized_gene_counts[c(GBM_tumor_core_neoplastic_astrocytes_id,GBM_tumor_periphery_neoplastic_astrocytes_id),]
dim(xdata)
ydata <- c(rep(1,length(GBM_tumor_core_neoplastic_astrocytes_id)),rep(0,length(GBM_tumor_periphery_neoplastic_astrocytes_id))) 
length(ydata)

xdata_sd <- sapply(seq(ncol(xdata)), function(ix) {sd(xdata[,ix])})
xdata <- xdata[,xdata_sd != 0]
dim(xdata)

# cross-validation for finding optimum alpha for a random given sample
my_alpha <- seq(0.1,0.9,0.1)

# Optimizing alpha and lambda for logistic regression with Elastic net regularization
nvar_selected_EN <- matrix(0,1,length(my_alpha))
pred_EN_cv <- matrix(0,dim(xdata)[1],length(my_alpha))
MSE_EN_cv <- matrix(0,1,length(my_alpha))

# assigning samples to folds, to be used in cross-validation when tunning alpha
set.seed(2010)
my_foldid <- sample(1:10,size=length(ydata),replace=TRUE)

library(glmnet)

for (j in 1:length(my_alpha)){

    # Logistic model fitting with 10-fold cross-validation for glmnet:
    fit_EN_cv <- cv.glmnet(as.matrix(xdata),as.factor(ydata),family="binomial",foldid=my_foldid,alpha=my_alpha[j],type.measure="mse")
    var_selected_EN <- which(fit_EN_cv$glmnet.fit$beta[,which(fit_EN_cv$cvm == min(fit_EN_cv$cvm))] != 0)
    nvar_selected_EN[j] <- length(var_selected_EN)

    # Predictions obtained by model i:
    pred_EN_cv[,j] <- predict(fit_EN_cv,as.matrix(xdata),s=fit_EN_cv$lambda[which(fit_EN_cv$cvm == min(fit_EN_cv$cvm))],type="response")

    # Mean squared error of prediction (MSE):
    MSE_EN_cv[j] <- mean((ydata-pred_EN_cv[,j])^2)
}

MSE_EN_cv
nvar_selected_EN
my_alpha[which(MSE_EN_cv == min(MSE_EN_cv))]


# Performing 1000 train and test boostrap samples
times_boot <- 1000

nvar_selected <- vector()
miscl_train <- vector()
miscl_test <- vector()
mse_train <- vector()
mse_test <- vector()
var_selected <- vector()
var_selected_idx <- vector()
pr_auc_train <- vector()
pr_auc_test <- vector()

new_test_id_matrix <- matrix(0,round(sum(ydata)*0.25)+round((length(ydata)-sum(ydata))*0.25),times_boot) 

# building random selected training and test sets (25 % assigned for testing)
set.seed(2019)
for (i in 1:times_boot){
new_test_id_matrix[,i] <- c(sample(1:sum(ydata),round(sum(ydata)*0.25),replace=FALSE),sample((sum(ydata)+1):length(ydata),round((length(ydata)-sum(ydata))*0.25),replace=FALSE))
}

for (i in 1:times_boot){

new_xtrain <- xdata[-new_test_id_matrix[,i],]
new_ytrain <- ydata[-new_test_id_matrix[,i]]

new_xtest <- xdata[new_test_id_matrix[,i],]
new_ytest <- ydata[new_test_id_matrix[,i]]

# Classification by sparse logistic regression

fit_EN_cv <- cv.glmnet(new_xtrain, as.factor(new_ytrain), family="binomial", nfolds=10, alpha=my_alpha[which(MSE_EN_cv == min(MSE_EN_cv))], type.measure="mse")

var_selected_EN_idx <- which(fit_EN_cv$glmnet.fit$beta[,which(fit_EN_cv$cvm == min(fit_EN_cv$cvm))] != 0)

nvar_selected[i] <- length(var_selected_EN_idx)

var_selected_idx <- c(var_selected_idx,var_selected_EN_idx)

var_selected <- c(var_selected,names(var_selected_idx))

# misclassifications, MSE, ROC and PR

# train set
pred_EN_train <- predict(fit_EN_cv,new_xtrain,type="response")
table(new_ytrain,round(pred_EN_train))
miscl_train[i] <- length(which(new_ytrain !=round(pred_EN_train)))

mse_train[i] <- mean((new_ytrain - pred_EN_train)^2)

library("PRROC")
pr_auc_train[i] <- pr.curve(scores.class0=as.vector(round(pred_EN_train)), weights.class0=new_ytrain)$auc.integral

# test set
pred_EN_test <- predict(fit_EN_cv,new_xtest,type="response")
table(new_ytest,round(pred_EN_test))
miscl_test[i] <- length(which(new_ytest !=round(pred_EN_test)))

mse_test[i] <- mean((new_ytest - pred_EN_test)^2)

pr_auc_test[i] <- pr.curve(scores.class0=as.vector(round(pred_EN_test)), weights.class0=new_ytest)$auc.integral

rm(new_xtrain,new_xtest,new_ytrain,new_ytest)
}


# median number of variables selected
median(nvar_selected)

# median number of misclassifications in the train set
median(round(miscl_train))

# median number of misclassifications in the test set
median(round(miscl_test))

# median squared error in the train set
median(mse_train)

# median squared error in the test set
median(mse_test)

# median PR AUC in the train set
median(pr_auc_train)

# median PR AUC in the test set
median(pr_auc_test)

# variables always selected by penalized logistic regression in the 1000 bootstrap samples
var_selected_alw_model1 <- as.numeric(names(which(table(var_selected_idx)>999)))
var_selected_alw_code_model1 <- colnames(xdata[,var_selected_alw_model1])
sort(var_selected_alw_code_model1)

# variables selected in more than 75% of the bootstrap runs by penalized logistic regression
var_selected_75_model1 <- as.numeric(names(which(table(var_selected_idx)>750)))
var_selected_75_code_model1 <- colnames(xdata[,var_selected_75_model1])
sort(var_selected_75_code_model1)
```

### Discrimination between periphery and neoplastic periphery astrocytes

```{r}
# Building the dataset
xdata <- GBM_normalized_gene_counts[c(GBM_tumor_periphery_astrocytes_id,GBM_tumor_periphery_neoplastic_astrocytes_id),]
dim(xdata)
ydata <- c(rep(1,length(GBM_tumor_periphery_astrocytes_id)),rep(0,length(GBM_tumor_periphery_neoplastic_astrocytes_id))) 
length(ydata)

xdata_sd <- sapply(seq(ncol(xdata)), function(ix) {sd(xdata[,ix])})
xdata <- xdata[,xdata_sd != 0]
dim(xdata)

# cross-validation for finding optimum alpha
my_alpha <- seq(0.1,0.9,0.1)

# Optimizing alpha and lambda for logistic regression with Elastic net regularization
nvar_selected_EN <- matrix(0,1,length(my_alpha))
pred_EN_cv <- matrix(0,dim(xdata)[1],length(my_alpha))
MSE_EN_cv <- matrix(0,1,length(my_alpha))

# assigning samples to folds, to be used in cross-validation when tunning alpha
set.seed(2013)
my_foldid <- sample(1:10,size=length(ydata),replace=TRUE)

library(glmnet)

for (j in 1:length(my_alpha)){

    # Logistic model fitting with 10-fold cross-validation for glmnet:
    fit_EN_cv <- cv.glmnet(as.matrix(xdata),as.factor(ydata),family="binomial",foldid=my_foldid,alpha=my_alpha[j],type.measure="mse")
    var_selected_EN <- which(fit_EN_cv$glmnet.fit$beta[,which(fit_EN_cv$cvm == min(fit_EN_cv$cvm))] != 0)
    nvar_selected_EN[j] <- length(var_selected_EN)

    # Predictions obtained by model i:
    pred_EN_cv[,j] <- predict(fit_EN_cv,as.matrix(xdata),s=fit_EN_cv$lambda[which(fit_EN_cv$cvm == min(fit_EN_cv$cvm))],type="response")

    # Mean squared error of prediction (MSE):
    MSE_EN_cv[j] <- mean((ydata-pred_EN_cv[,j])^2)
}

MSE_EN_cv
nvar_selected_EN
my_alpha[which(MSE_EN_cv == min(MSE_EN_cv))]


# Performing 1000 train and test boostrap samples
times_boot <- 1000

nvar_selected <- vector()
miscl_train <- vector()
miscl_test <- vector()
mse_train <- vector()
mse_test <- vector()
var_selected <- vector()
var_selected_idx <- vector()
pr_auc_train <- vector()
pr_auc_test <- vector()

new_test_id_matrix <- matrix(0,round(sum(ydata)*0.25)+round((length(ydata)-sum(ydata))*0.25),times_boot) 

# building random selected training and test sets (25 % assigned for testing)
set.seed(2019)
for (i in 1:times_boot){
new_test_id_matrix[,i] <- c(sample(1:sum(ydata),round(sum(ydata)*0.25),replace=FALSE),sample((sum(ydata)+1):length(ydata),round((length(ydata)-sum(ydata))*0.25),replace=FALSE))
}

for (i in 1:times_boot){

new_xtrain <- xdata[-new_test_id_matrix[,i],]
new_ytrain <- ydata[-new_test_id_matrix[,i]]

new_xtest <- xdata[new_test_id_matrix[,i],]
new_ytest <- ydata[new_test_id_matrix[,i]]

# Classification by sparse logistic regression

fit_EN_cv <- cv.glmnet(new_xtrain, as.factor(new_ytrain), family="binomial", nfolds=10, alpha=my_alpha[which(MSE_EN_cv == min(MSE_EN_cv))], type.measure="mse")

var_selected_EN_idx <- which(fit_EN_cv$glmnet.fit$beta[,which(fit_EN_cv$cvm == min(fit_EN_cv$cvm))] != 0)

nvar_selected[i] <- length(var_selected_EN_idx)

var_selected_idx <- c(var_selected_idx,var_selected_EN_idx)

var_selected <- c(var_selected,names(var_selected_idx))

# misclassifications, MSE, ROC and PR

# train set
pred_EN_train <- predict(fit_EN_cv,new_xtrain,type="response")
table(new_ytrain,round(pred_EN_train))
miscl_train[i] <- length(which(new_ytrain !=round(pred_EN_train)))

mse_train[i] <- mean((new_ytrain - pred_EN_train)^2)


pr_auc_train[i] <- pr.curve(scores.class0=as.vector(round(pred_EN_train)), weights.class0=new_ytrain)$auc.integral

# test set
pred_EN_test <- predict(fit_EN_cv,new_xtest,type="response")
table(new_ytest,round(pred_EN_test))
miscl_test[i] <- length(which(new_ytest !=round(pred_EN_test)))

mse_test[i] <- mean((new_ytest - pred_EN_test)^2)

pr_auc_test[i] <- pr.curve(scores.class0=as.vector(round(pred_EN_test)), weights.class0=new_ytest)$auc.integral

rm(new_xtrain,new_xtest,new_ytrain,new_ytest)
}

# median number of variables selected
median(nvar_selected)

# median number of misclassifications in the train set
median(round(miscl_train))

# median number of misclassifications in the test set
median(round(miscl_test))

# median squared error in the train set
median(mse_train)

# median squared error in the test set
median(mse_test)

# median PR AUC in the train set
median(pr_auc_train)

# median PR AUC in the test set
median(pr_auc_test)

# variables always selected by penalized logistic regression in the 1000 bootstrap samples
var_selected_alw_model2 <- as.numeric(names(which(table(var_selected_idx)>999)))
var_selected_alw_code_model2 <- colnames(xdata[,var_selected_alw_model2])
sort(var_selected_alw_code_model2)

# variables selected in more than 75% of the bootstrap runs by penalized logistic regression
var_selected_75_model2 <- as.numeric(names(which(table(var_selected_idx)>750)))
var_selected_75_code_model2 <- colnames(xdata[,var_selected_75_model2])
sort(var_selected_75_code_model2)
```

## Correlation-based regularized classification

### Discrimination between neoplastic periphery + core astrocytes and periphery astrocytes

```{r,warning=FALSE,message=FALSE}
# correlation matrices
tumor_core_neoplastic_astrocytes_aux <- GBM_normalized_gene_counts[GBM_tumor_core_neoplastic_astrocytes_id,sapply(seq(ncol(GBM_normalized_gene_counts)), function(ix) {sd(GBM_normalized_gene_counts[GBM_tumor_core_neoplastic_astrocytes_id,ix])}) != 0]
dim(tumor_core_neoplastic_astrocytes_aux)
 
tumor_periphery_neoplastic_astrocytes_aux <- GBM_normalized_gene_counts[GBM_tumor_periphery_neoplastic_astrocytes_id,sapply(seq(ncol(GBM_normalized_gene_counts)), function(ix) {sd(GBM_normalized_gene_counts[GBM_tumor_periphery_neoplastic_astrocytes_id,ix])}) != 0]
dim(tumor_periphery_neoplastic_astrocytes_aux)
 
tumor_core_neoplastic_astrocytes_less <- tumor_core_neoplastic_astrocytes_aux[,which(colnames(tumor_core_neoplastic_astrocytes_aux) %in% colnames(tumor_periphery_neoplastic_astrocytes_aux))]
dim(tumor_core_neoplastic_astrocytes_less)

tumor_periphery_neoplastic_astrocytes_less <- tumor_periphery_neoplastic_astrocytes_aux[,which(colnames(tumor_periphery_neoplastic_astrocytes_aux) %in% colnames(tumor_core_neoplastic_astrocytes_aux))]
dim(tumor_periphery_neoplastic_astrocytes_less)

tumor_periphery_astrocytes_aux <- GBM_normalized_gene_counts[GBM_tumor_periphery_astrocytes_id,]
dim(tumor_periphery_astrocytes_aux)

tumor_periphery_astrocytes_less <- tumor_periphery_astrocytes_aux[,which(colnames(tumor_periphery_astrocytes_aux) %in% colnames(tumor_periphery_neoplastic_astrocytes_less))]
dim(tumor_periphery_astrocytes_less)

xdata <- rbind(tumor_core_neoplastic_astrocytes_less,tumor_periphery_neoplastic_astrocytes_less,tumor_periphery_astrocytes_less)
dim(xdata)

ydata <- c(rep(1,dim(tumor_core_neoplastic_astrocytes_less)[1]),rep(1,dim(tumor_periphery_neoplastic_astrocytes_less)[1]),rep(0,dim(tumor_periphery_astrocytes_less)[1])) 
length(ydata)

library("propagate")
library("lsa")
tumor_core_neoplastic_astrocytes_cor <- bigcor(tumor_core_neoplastic_astrocytes_less, y = NULL, fun = "cor", size = 2000, verbose=TRUE)
tumor_core_neoplastic_astrocytes_cor <- as.data.frame(as.ffdf(tumor_core_neoplastic_astrocytes_cor))

tumor_periphery_neoplastic_astrocytes_cor <- bigcor(tumor_periphery_neoplastic_astrocytes_less, y = NULL, fun = "cor", size = 2000, verbose=TRUE)
tumor_periphery_neoplastic_astrocytes_cor <- as.data.frame(as.ffdf(tumor_periphery_neoplastic_astrocytes_cor))

# Building a weight vector that penalizes genes with larger distances between the correlation matrices of neoplastic periphery astrocytes and neoplastic core astrocytes
weights <- vector()

for (i in 1:dim(xdata)[2]){
weights[i] <- acos(cosine(tumor_core_neoplastic_astrocytes_cor[,i],tumor_periphery_neoplastic_astrocytes_cor[,i]))/pi
}

# normalized weights (between 0 and 1)
weights_norm <- weights / max(weights)
weights <- weights_norm
```

Comparing model predictive performing either using network information or not, based on resampling (1000 random train and test samples)

```{r,warning=FALSE,message=FALSE}

# cross-validation for finding optimum alpha
my_alpha <- seq(0.1,0.9,0.1)

# Optimizing alpha and lambda for logistic regression with Elastic net regularization
nvar_selected_EN <- matrix(0,1,length(my_alpha))
pred_EN_cv <- matrix(0,dim(xdata)[1],length(my_alpha))
MSE_EN_cv <- matrix(0,1,length(my_alpha))

# assigning samples to folds, to be used in cross-validation when tunning alpha
set.seed(2013)
my_foldid <- sample(1:10,size=length(ydata),replace=TRUE)

library(glmnet)

for (j in 1:length(my_alpha)){

    # Logistic model fitting with 10-fold cross-validation for glmnet:
    fit_EN_cv <- cv.glmnet(as.matrix(xdata),as.factor(ydata),family="binomial",foldid=my_foldid,alpha=my_alpha[j],type.measure="mse")
    var_selected_EN <- which(fit_EN_cv$glmnet.fit$beta[,which(fit_EN_cv$cvm == min(fit_EN_cv$cvm))] != 0)
    nvar_selected_EN[j] <- length(var_selected_EN)

    # Predictions obtained by model i:
    pred_EN_cv[,j] <- predict(fit_EN_cv,as.matrix(xdata),s=fit_EN_cv$lambda[which(fit_EN_cv$cvm == min(fit_EN_cv$cvm))],type="response")

    # Mean squared error of prediction (MSE):
    MSE_EN_cv[j] <- mean((ydata-pred_EN_cv[,j])^2)
}

MSE_EN_cv
nvar_selected_EN
my_alpha[which(MSE_EN_cv == min(MSE_EN_cv))]


library(glmnet)
library(PRROC)

# Performing 100 train and test boostrap samples
times_boot <- 1000

nvar_selected <- matrix(0,2,times_boot)
miscl_train <- matrix(0,4,times_boot)
miscl_test <- matrix(0,4,times_boot)
mse_train <- matrix(0,4,times_boot)
mse_test <- matrix(0,4,times_boot)
var_selected <- vector()
var_selected_netw <- vector()
var_selected_idx <- vector()
var_selected_netw_idx <- vector()
pr_auc_train <- matrix(0,4,times_boot)
pr_auc_test <- matrix(0,4,times_boot)

new_test_id_matrix <- matrix(0,round(sum(ydata)*0.25)+round((length(ydata)-sum(ydata))*0.25),times_boot) 

# building random selected training and test sets (25 % assigned for testing)
set.seed(2019)
for (i in 1:times_boot){
new_test_id_matrix[,i] <- c(sample(1:sum(ydata),round(sum(ydata)*0.25),replace=FALSE),sample((sum(ydata)+1):length(ydata),round((length(ydata)-sum(ydata))*0.25),replace=FALSE))
}

for (i in 1:times_boot){

new_xtrain <- xdata[-new_test_id_matrix[,i],]
new_ytrain <- ydata[-new_test_id_matrix[,i]]

new_xtest <- xdata[new_test_id_matrix[,i],]
new_ytest <- ydata[new_test_id_matrix[,i]]

# assigning samples to folds
set.seed(40)
my_foldid <- sample(1:10,size=dim(xdata[-new_test_id_matrix[,i],])[1]
,replace=TRUE)

# Classification by sparse logistic regression

## with the elastic net (EN) penalty

fit_EN_cv <- cv.glmnet(new_xtrain, as.factor(new_ytrain), family="binomial", nfolds=10, alpha=0.9, foldid=my_foldid, type.measure="mse")

var_selected_EN_idx <- which(fit_EN_cv$glmnet.fit$beta[,which(fit_EN_cv$cvm == min(fit_EN_cv$cvm))] != 0)

nvar_selected[1,i] <- length(var_selected_EN_idx)

var_selected_idx <- c(var_selected_idx,var_selected_EN_idx)

var_selected <- c(var_selected,names(var_selected_idx))

# misclassifications, MSE, ROC and PR

# train set
pred_EN_train <- predict(fit_EN_cv,new_xtrain,type="response")
table(new_ytrain,round(pred_EN_train))
miscl_train[1,i] <- length(which(new_ytrain !=round(pred_EN_train)))

mse_train[1,i] <- mean((new_ytrain - pred_EN_train)^2)

pr_auc_train[1,i] <- pr.curve(scores.class0=as.vector(round(pred_EN_train)), weights.class0=new_ytrain)$auc.integral

# test set
pred_EN_test <- predict(fit_EN_cv,new_xtest,type="response")
table(new_ytest,round(pred_EN_test))
miscl_test[1,i] <- length(which(new_ytest !=round(pred_EN_test)))

mse_test[1,i] <- mean((new_ytest - pred_EN_test)^2)

pr_auc_test[1,i] <- pr.curve(scores.class0=as.vector(round(pred_EN_test)), weights.class0=new_ytest)$auc.integral

## with network-based information

fit_ENnetw_cv <- cv.glmnet(new_xtrain, as.factor(new_ytrain), family="binomial", nfolds=10, alpha=0.9, foldid=my_foldid, penalty.factor=weights,type.measure="mse")

var_selected_ENnetw_idx <- which(fit_ENnetw_cv$glmnet.fit$beta[,which(fit_ENnetw_cv$cvm == min(fit_ENnetw_cv$cvm))] != 0)

nvar_selected[2,i] <- length(var_selected_ENnetw_idx)

var_selected_netw_idx <- c(var_selected_netw_idx,var_selected_ENnetw_idx)

var_selected_netw <- c(var_selected_netw,names(var_selected_idx))

# misclassifications, MSE, ROC and PR 

# train set
pred_ENnetw_train <- predict(fit_ENnetw_cv,new_xtrain,type="response")
table(new_ytrain,round(pred_ENnetw_train))
miscl_train[2,i] <- length(which(new_ytrain !=round(pred_ENnetw_train)))

mse_train[2,i] <- mean((new_ytrain - pred_ENnetw_train)^2)

pr_auc_train[2,i] <- pr.curve(scores.class0=as.vector(round(pred_ENnetw_train)), weights.class0=new_ytrain)$auc.integral

# test set
pred_ENnetw_test <- predict(fit_ENnetw_cv,new_xtest,type="response")
table(new_ytest,round(pred_ENnetw_test))
miscl_test[2,i] <- length(which(new_ytest !=round(pred_ENnetw_test)))

mse_test[2,i] <- mean((new_ytest - pred_ENnetw_test)^2)

pr_auc_test[2,i] <- pr.curve(scores.class0=as.vector(round(pred_ENnetw_test)), weights.class0=new_ytest)$auc.integral


# Naive Bayes Classifier

## based on the variables selected by EN

library(naivebayes)
fit_NB_EN <- naive_bayes(new_xtrain[,var_selected_EN_idx],as.factor(new_ytrain))

# misclassifications, MSE, ROC and PR

# train set

pred_NB_train <- predict(fit_NB_EN, new_xtrain[,var_selected_EN_idx], type='prob')

miscl_train[3,i] <- length(which(new_ytrain != round(pred_NB_train[,2])))

mse_train[3,i] <- mean((new_ytrain - pred_NB_train[,2])^2)

pr_auc_train[3,i] <- pr.curve(scores.class0=as.vector(round(pred_NB_train[,2])), weights.class0=new_ytrain)$auc.integral

# test set

pred_NB_test <- predict(fit_NB_EN, new_xtest[,var_selected_EN_idx], type='prob')

miscl_test[3,i] <- length(which(new_ytest != round(pred_NB_test[,2])))

mse_test[3,i] <- mean((new_ytest - pred_NB_test[,2])^2)

pr_auc_test[3,i] <- pr.curve(scores.class0=as.vector(round(pred_NB_test[,2])), weights.class0=new_ytest)$auc.integral


## based on the variables selected by twiner

fit_NB_twiner <- naive_bayes(new_xtrain[,var_selected_ENnetw_idx],as.factor(new_ytrain))

# misclassifications, MSE, ROC and PR

# train set

pred_NB_train <- predict(fit_NB_twiner, new_xtrain[,var_selected_ENnetw_idx], type='prob')

miscl_train[4,i] <- length(which(new_ytrain != round(pred_NB_train[,2])))

mse_train[4,i] <- mean((new_ytrain - pred_NB_train[,2])^2)

pr_auc_train[4,i] <- pr.curve(scores.class0=as.vector(round(pred_NB_train[,2])), weights.class0=new_ytrain)$auc.integral

# test set

pred_NB_test <- predict(fit_NB_twiner, new_xtest[,var_selected_ENnetw_idx], type='prob')

miscl_test[4,i] <- length(which(new_ytest != round(pred_NB_test[,2])))

mse_test[4,i] <- mean((new_ytest - pred_NB_test[,2])^2)

pr_auc_test[4,i] <- pr.curve(scores.class0=as.vector(round(pred_NB_test[,2])), weights.class0=new_ytest)$auc.integral


rm(new_xtrain,new_xtest,new_ytrain,new_ytest)
}
```

Summary numbers for each modeling strategy

```{r}
# median number of variables selected
apply(nvar_selected,1,median)

# median number of misclassifications in the train set
round(apply(round(miscl_train),1,median))

# median number of misclassifications in the test set
round(apply(round(miscl_test),1,median))

# median squared error in the train set
apply(mse_train,1,median)

# median squared error in the test set
apply(mse_test,1,median)

# median PR AUC in the train set
apply(pr_auc_train,1,median)

# median PR AUC in the test set
apply(pr_auc_test,1,median)



# violin plots fro MSE and AUC
MSE_violin_data <- data.frame('Model'=c(rep("EN train",dim(mse_train)[2]),rep("Twiner train",dim(mse_train)[2]),rep("EN test",dim(mse_train)[2]),rep("Twiner test",dim(mse_train)[2])),'MSE'=c(mse_train[1,],mse_train[2,],mse_test[1,],mse_test[2,]))
ggplot(data=MSE_violin_data, aes(x=Model, y=MSE)) +
geom_violin(aes(fill=Model, color=Model)) +
theme_minimal() +
scale_colour_manual(values=c("lightsteelblue","royalblue","rosybrown1","indianred3"))+
scale_fill_manual(values=c("lightsteelblue","royalblue","rosybrown1","indianred3")) +
theme(axis.text=element_text(size=14),axis.title=element_text(size=14),legend.position='none')

pr_auc_violin_data <- data.frame('Model'=c(rep("EN train",dim(pr_auc_train)[2]),rep("Twiner train",dim(pr_auc_train)[2]),rep("EN test",dim(pr_auc_train)[2]),rep("Twiner test",dim(pr_auc_train)[2])),'AUC'=c(pr_auc_train[1,],pr_auc_train[2,],pr_auc_test[1,],pr_auc_test[2,]))
ggplot(data=pr_auc_violin_data, aes(x=Model, y=AUC)) +
geom_violin(aes(fill=Model, color=Model)) +
theme_minimal() +
scale_colour_manual(values=c("lightsteelblue","royalblue","rosybrown1","indianred3"))+
scale_fill_manual(values=c("lightsteelblue","royalblue","rosybrown1","indianred3")) +
theme(axis.text=element_text(size=14),axis.title=element_text(size=14),legend.position='none')
```

Genes selected

```{r}
# variables always selected by penalized logistic regression in the 100 bootstrap samples
var_selected_alw_model3 <- as.numeric(names(which(table(var_selected_idx)>990)))
var_selected_alw_model3_code <- colnames(xdata[,var_selected_alw_model3])
sort(var_selected_alw_model3_code)

# variables selected in more than 75% of the bootstrap runs by penalized logistic regression
var_selected_75_model3 <- as.numeric(names(which(table(var_selected_idx)>750)))
var_selected_75_model3_code <- colnames(xdata[,var_selected_75_model3])
sort(var_selected_75_model3_code)

# variables always selected by network-based penalized logistic regression in the 100 bootstrap samples
var_selected_netw_alw_model3 <- as.numeric(names(which(table(var_selected_netw_idx)>999)))
var_selected_netw_alw_model3_code <- colnames(xdata[,var_selected_netw_alw_model3])
sort(var_selected_netw_alw_model3_code)

# variables selected in more than 75% of the bootstrap runs by network-based penalized logistic regression
var_selected_netw_75_model3 <- as.numeric(names(which(table(var_selected_netw_idx)>750)))
var_selected_netw_75_model3_code <- colnames(xdata[,var_selected_netw_75_model3])
sort(var_selected_netw_75_model3_code)

# variables selected between penalized logistic regression and network-based penalized logistic regression
common_var_selected_netw_75_model3 <- var_selected_netw_75_model3[which(var_selected_netw_75_model3 %in% var_selected_75_model3)]
common_var_selected_netw_75_model3_code <- var_selected_netw_75_model3_code[which(var_selected_netw_75_model3_code %in% var_selected_75_model3_code)]
sort(common_var_selected_netw_75_model3_code)

# distinct variables selected by penalized logistic regression
distinct_var_selected_75_model3 <- var_selected_75_model3[-which(var_selected_75_model3 %in% var_selected_netw_75_model3)]
distinct_var_selected_75_model3_code <- var_selected_75_model3_code[-which(var_selected_75_model3_code %in% var_selected_netw_75_model3_code)]
sort(distinct_var_selected_75_model3_code)

# distinct variables selected by network-based penalized logistic regression
distinct_var_selected_netw_75_model3 <- var_selected_netw_75_model3[-which(var_selected_netw_75_model3 %in% var_selected_75_model3)]
distinct_var_selected_netw_75_model3_code <- var_selected_netw_75_model3_code[-which(var_selected_netw_75_model3_code %in% var_selected_75_model3_code)]
sort(distinct_var_selected_netw_75_model3_code)

# plot weights

weights_color <- c(rep("EN",length(distinct_var_selected_75_model3)),rep("twiner",length(distinct_var_selected_netw_75_model3)),rep("common",length(common_var_selected_netw_75_model3)))

weights_plot <- data.frame(Gene = c(distinct_var_selected_75_model3_code,distinct_var_selected_netw_75_model3_code,common_var_selected_netw_75_model3_code), Weight = weights[c(distinct_var_selected_75_model3,distinct_var_selected_netw_75_model3,common_var_selected_netw_75_model3)], label = factor(weights_color))

ggplot(weights_plot, aes(Gene, Weight)) + geom_point(aes(color = label),size=1.8) + 
scale_colour_manual(values=c("gray55","royalblue4","firebrick3"))+
theme_minimal() +
theme(axis.text=element_text(size=12),axis.title=element_text(size=14),legend.title = element_text(size=12, color="white"),legend.text = element_text(size=16),axis.text.x = element_text(angle = 90,vjust=0.5, hjust=1)) 
  

# correlation networks

# Gene networks
qgraph_matrix <- xdata

qgraph_vars <- c(distinct_var_selected_netw_75_model3,common_var_selected_netw_75_model3)

# tumor core neoplastic astrocytes
qgraph_cor_tumor_core <- cor(qgraph_matrix[1:length(GBM_tumor_core_neoplastic_astrocytes_id),qgraph_vars], method="pearson")

# tumor periphery neoplastic astrocytes
qgraph_cor_tumor_periphery <- cor(qgraph_matrix[(length(GBM_tumor_core_neoplastic_astrocytes_id)+1):(length(GBM_tumor_core_neoplastic_astrocytes_id)+length(GBM_tumor_periphery_neoplastic_astrocytes_id)),qgraph_vars], method="pearson")

# tumor periphery normal astrocytes
qgraph_cor_normal_periphery <- cor(qgraph_matrix[(length(GBM_tumor_core_neoplastic_astrocytes_id)+length(GBM_tumor_periphery_neoplastic_astrocytes_id)+1):(length(GBM_tumor_core_neoplastic_astrocytes_id)+length(GBM_tumor_periphery_neoplastic_astrocytes_id)+length(GBM_tumor_periphery_astrocytes_id)),qgraph_vars], method="pearson")

qgraph_label <- c(distinct_var_selected_netw_75_model3_code,common_var_selected_netw_75_model3_code)

qgraph_groups <- list(netEN = c(1:length(distinct_var_selected_netw_75_model3_code)), common = (length(distinct_var_selected_netw_75_model3_code)+1):(length(distinct_var_selected_netw_75_model3_code)+length(common_var_selected_netw_75_model3_code)))


library(qgraph)

qgraph(qgraph_cor_tumor_core, layout = "circle", minimum = 0.2, vsize = 1, groups = qgraph_groups, legend=FALSE,borders = FALSE, node.width=4, color=c("#FF0000","#999999"),labels = qgraph_label,label.cex=2, theme="TeamFortress",title="T-core",title.cex=2)

qgraph(qgraph_cor_tumor_periphery, layout = "circle", minimum = 0.2, vsize = 1.1, groups = qgraph_groups, legend=FALSE,borders = FALSE, node.width=4, color=c("#FF0000","#999999"),labels = qgraph_label,label.cex=2,theme="TeamFortress",title="T-periphery",title.cex=2)

qgraph(qgraph_cor_normal_periphery, layout = "circle", minimum = 0.2, vsize = 1.1, groups = qgraph_groups, legend=FALSE,borders = FALSE, node.width=4, color=c("#FF0000","#999999"),labels = qgraph_label, label.cex=2,theme="TeamFortress",title="N-periphery",title.cex=2)
```

# Survival analysis

```{r}
# GBM RNA-seq data from the TCGA

library(SummarizedExperiment)
library(TCGAbiolinks)

query_exp <- GDCquery(project = "TCGA-GBM", 
                      legacy = TRUE,
                      data.category = "Gene expression",
                      data.type = "Gene expression quantification",
                      platform = "Illumina HiSeq", 
                      file.type = "normalized_results",
                      experimental.strategy = "RNA-Seq",
                      sample.type = "Primary Tumor")

GDCdownload(query_exp)

GBM_exp <- GDCprepare(query = query_exp, save = TRUE, save.filename = "GBMExp.rda")

# searching for possible outliers using the TCGAanalyze_Preprocessing function, which performs an Array Array Intensity correlation AAIC
data_prep <- TCGAanalyze_Preprocessing(object = GBM_exp, cor.cut = 0.6)                      
# using the TCGAanalyze_Normalization function for normalizing mRNA transcripts and miRNA
data_norm <- TCGAanalyze_Normalization(tabDF = data_prep,
                                      geneInfo = geneInfo,
                                      method = "gcContent")                

# using the TCGAanalyze_Filtering function to remove features
data_filt <- t(TCGAanalyze_Filtering(tabDF = data_norm,
                                    method = "quantile", 
                                    qnt.cut =  0.25))   

# information and clinical status and survival time / days to last follow-up
time_status <- cbind(GBM_exp$days_to_death,GBM_exp$vital_status)
time_status[which(time_status[,2]=="Alive"),1] <- GBM_exp$days_to_last_follow_up[which(time_status[,2]=="Alive")]

rownames(time_status) <- GBM_exp$barcode
data_filt <- data_filt[-which(is.na(time_status[,1])),]

time_status <- time_status[-which(is.na(time_status[,1])),]
time_status[which(time_status[,2] == "Dead"),2] <- 1
time_status[which(time_status[,2] == "Alive"),2] <- 0

time_status <- data.frame(time_status)
time_status[,1] <- as.numeric(as.character(time_status[,1]))
time_status[,2] <- as.numeric(as.character(time_status[,2]))

colnames(time_status) <- c("time","status")

data_filt <- data_filt[which(rownames(data_filt) %in% rownames(time_status)),]
time_status <- time_status[rownames(data_filt),]

library(survival)
library(survminer)

my_fold <- sample(1:10,size=dim(data_filt)[1],replace=TRUE)

# survival, neoplastic core v. neoplastic periphery
data_filt_EN <- data_filt[,which(colnames(data_filt) %in% var_selected_75_code_model1)]

fit_cox_EN_cv <- cv.glmnet(data_filt_EN, Surv(as.numeric(time_status[,1]),as.numeric(time_status[,2])), family="cox", alpha=0, foldid=my_fold)
fit_cox_EN_coefficients <- as.numeric(coef(fit_cox_EN_cv, s=fit_cox_EN_cv$lambda.min))
length(colnames(data_filt_EN[,which(fit_cox_EN_coefficients != 0)]))

fit_cox_EN_risk <- data_filt_EN %*% fit_cox_EN_coefficients
fit_cox_EN_risk_group <- fit_cox_EN_risk
fit_cox_EN_risk_group[which(fit_cox_EN_risk > median(fit_cox_EN_risk))] <- 2
fit_cox_EN_risk_group[which(fit_cox_EN_risk <= median(fit_cox_EN_risk))] <- 1

KM_mat <- cbind(time_status,fit_cox_EN_risk_group)
fit_KM <- survfit(Surv(as.numeric(time),as.numeric(status)) ~ fit_cox_EN_risk_group, data=as.data.frame(KM_mat), type="kaplan-meier", conf.type="plain")

survdiff(Surv(as.numeric(time),as.numeric(status)) ~ fit_cox_EN_risk_group, data=as.data.frame(KM_mat))

ggsurvplot(fit_KM,title="T-core vs. T-periphery",font.main = c(14, "bold", "black"),legend.labs = c("low","high"), legend.title = "Risk group",pval=TRUE,palette = c("#999999","#993366"), xlab= "Time (days)", pval.coord=c(1800,0.25), pval.size=4.5)


# survival, neoplastic periphery vs. normal periphery
data_filt_EN <- data_filt[,which(colnames(data_filt) %in% var_selected_75_code_model2)]

fit_cox_EN_cv <- cv.glmnet(data_filt_EN, Surv(as.numeric(time_status[,1]),as.numeric(time_status[,2])), family="cox", alpha=0, foldid=my_fold)
fit_cox_EN_coefficients <- as.numeric(coef(fit_cox_EN_cv, s=fit_cox_EN_cv$lambda.min))
length(colnames(data_filt_EN[,which(fit_cox_EN_coefficients != 0)]))

fit_cox_EN_risk <- data_filt_EN %*% fit_cox_EN_coefficients
fit_cox_EN_risk_group <- fit_cox_EN_risk
fit_cox_EN_risk_group[which(fit_cox_EN_risk > median(fit_cox_EN_risk))] <- 2
fit_cox_EN_risk_group[which(fit_cox_EN_risk <= median(fit_cox_EN_risk))] <- 1

KM_mat <- cbind(time_status,fit_cox_EN_risk_group)
fit_KM <- survfit(Surv(as.numeric(time),as.numeric(status)) ~ fit_cox_EN_risk_group, data=as.data.frame(KM_mat), type="kaplan-meier", conf.type="plain")

survdiff(Surv(as.numeric(time),as.numeric(status)) ~ fit_cox_EN_risk_group, data=as.data.frame(KM_mat))

ggsurvplot(fit_KM,title="T-periphery vs. N-periphery",font.main = c(14, "bold", "black"),legend.labs = c("low","high"), legend.title = "Risk group",pval=TRUE,palette = c("#999999","#993366"), xlab= "Time (days)", pval.coord=c(1800,0.25), pval.size=4.5)


# Survival based on twiner variables, neoplastic core/periphery vs. normal periphery
data_filt_twiner <- data_filt[,which(colnames(data_filt) %in% var_selected_netw_75_model3_code)]

fit_cox_twiner_cv <- cv.glmnet(data_filt_twiner, Surv(time_status$time,time_status$status), family="cox", alpha=0, foldid=my_fold)
fit_cox_twiner_coefficients <- as.numeric(coef(fit_cox_twiner_cv, s=fit_cox_twiner_cv$lambda.min))

fit_cox_twiner_risk <- data_filt_twiner %*% fit_cox_twiner_coefficients
fit_cox_twiner_risk_group <- fit_cox_twiner_risk
fit_cox_twiner_risk_group[which(fit_cox_twiner_risk > median(fit_cox_twiner_risk))] <- 2
fit_cox_twiner_risk_group[which(fit_cox_twiner_risk <= median(fit_cox_twiner_risk))] <- 1

KM_mat <- cbind(time_status,fit_cox_twiner_risk_group)
fit_KM <- survfit(Surv(time,status) ~ fit_cox_twiner_risk_group, data=KM_mat, type="kaplan-meier", conf.type="plain")

survdiff(Surv(time,status) ~ fit_cox_twiner_risk_group, data=KM_mat)

ggsurvplot(fit_KM,title="T-core/periphery vs. N-periphery",font.main = c(14, "bold", "black"),legend.labs = c("low","high"), legend.title = "Risk group",pval=TRUE,palette = c("#999999","#993366"), xlab= "Time (days)", pval.coord=c(1800,0.25), pval.size=4.5)

```

