Evaluation_Survival_Analysis_Inter_Dataset.R

#########################################################################
###      Evaluation of Survial Analysis Models (Training - Test)      ###
#########################################################################


setwd('~/bigdata/PCa/')
source('script/PCa_Functions.R')
source('~/bigdata/LABDATA/BLUPHAT/script/Helper_Functions.R')


# BRB-ArrayTools: https://brb.nci.nih.gov/BRB-ArrayTools/download.html
# randomForestSRC: https://kogalur.github.io/randomForestSRC/theory.html
# superpc: http://statweb.stanford.edu/~tibs/superpc/tutorial.html
# plsRcox: https://fbertran.github.io/plsRcox/
# BhGLM: https://github.com/nyiuab/BhGLM

#library(remotes)
#install_github("nyiuab/BhGLM", force=T, build_vignettes=T)
#BiocManager::install('survcomp')

library(Biobase)
library(readxl)
library(survminer)
library(survcomp)

library(survival) # coxph
library(glmnet) # coxnet
library(superpc) # pca
#library(survivalsvm) # svm
#library(randomForestSRC)
#library(plsRcox)

#library(BhGLM)

# Elastic Net/Lasso-Cox
# survivalSVM
# randomForestSRC
# plsRcox
# Superpc

# Bayesian Losso
# Neural Network Cox
# Dr. Yi


# TCGA-PRAD     == Cell 2015                == GDC         == RNAseq: CPM              == Time to BCR
# CPC-GENE      == Nature_2017              == GSE107299   == HuGene 2.0 ST & HTA 2.0  == Time to BCR
# MSKCC         ==_Cancer_Cell_2010         == GSE21034    == HuEx 1.0 ST              == Time to BCR
# DKFZ          == Cancer Cell 2018         == cBioPortal  == RNAseq: RPKM             == Time to BCR
# GSE54460      == Cancer Research 2014     == GSE54460    == RNAseq: CPM              == Time to BCR
# CamCap        == Nature Genetics 2016     == GSE70768    == Illumina HumanHT-12 V4.0 == Time to BCR
# Stockholm     == Nature Genetics 2016     == GSE70769    == Illumina HumanHT-12 V4.0 == Time to BCR
# DESNT         == Eur Urol Focus 2017      == GSE94767    == HuEx 1.0 ST              == Time to BCR
# E-MTAB-6128   == Annuals of Oncology 2018 == E-MTAB-6128 == HuGene 2.0 ST            == Time to BCR
# Belfast       == Journal of Clinical Oncology 2017  == GSE116918 == ADXPCv1a520642   == Time to BCR/Metastasis

datasets <- c('TCGA_PRAD','GSE107299','GSE21034','DKFZ2018','GSE54460','GSE70768','GSE70769',
              'GSE94767','E_MTAB_6128','GSE116918_BCR') # , 'GSE116918_Metastasis'

# Published signatures
signatures <- c('Agell','Bibikova','Bismar','Decipher','Ding','Glinsky','Irshad',
                'Jennifer','Jia','Kamoun','Long','Luca','Mo','Nakagawa','Olmos',
                'Oncotype','Penney','Planche','Prolaris','Ramaswamy','Ramos_Montoya',
                'Ross_Adams','Ross_Robert','Sharma','Talantov','Varambally','Wu','Yang',
                'Yu')

training.set <- datasets[1]
training.set
eSet <- readRDS(paste0('data/Database/Primary/', training.set, '_eSet.RDS'))
exprData <- exprs(eSet)
phenoData <- pData(eSet)
#View(phenoData)
table(phenoData$sample_type)

keep <- which(phenoData$sample_type=='Primary' | phenoData$sample_type=='Tumor' | phenoData$sample_type=='Tumour' )
exprData <- exprData[,keep]
phenoData <- phenoData[keep,]

filter <- which(duplicated(phenoData$patient_id, incomparables = NA))

if (length(filter)>0) {
  exprData <- exprData[,-filter]
  phenoData <- phenoData[-filter,]
}


if (training.set=='GSE116918_Metastasis') {
  phenoData$time_to_bcr <- phenoData$time_to_metastasis
  phenoData$bcr_status <- phenoData$metastasis_status
}

filter <- which(is.na(phenoData$bcr_status) | is.na(phenoData$time_to_bcr))

if (length(filter)>0) {
  exprData <- exprData[,-filter]
  phenoData <- phenoData[-filter,]
}

##### Smples used for classification
### ========================================= ###
# phenoData$bcr_status_5yr <- getEventFun(n=5, time.to.event=phenoData$time_to_bcr, status=phenoData$bcr_status)
# 
# table(phenoData$bcr_status_5yr)
# 
# filter <- which(is.na(phenoData$bcr_status_5yr))
# phenoData <- phenoData[-filter,]
# exprData <- exprData[,-filter]
### ========================================= ###

training.pheno <- phenoData
dim(training.pheno)
training.pheno$time_to_bcr[training.pheno$time_to_bcr<=0] <- 0.01

training.geno.all <- as.matrix(t(exprData))
training.geno.all <- scale(training.geno.all)

test.sets <- datasets[which(!datasets %in% training.set)]
test.sets

for (test.set in test.sets) {
  message (test.set)
  
  eSet <- readRDS(paste0('data/Database/Primary/', test.set, '_eSet.RDS'))
  exprData <- exprs(eSet)
  phenoData <- pData(eSet)
  #View(phenoData)
  table(phenoData$sample_type)
  
  keep <- which(phenoData$sample_type=='Primary' | phenoData$sample_type=='Tumor' | phenoData$sample_type=='Tumour' )
  exprData <- exprData[,keep]
  phenoData <- phenoData[keep,]
  
  filter <- which(duplicated(phenoData$patient_id, incomparables = NA))
  
  if (length(filter)>0) {
    exprData <- exprData[,-filter]
    phenoData <- phenoData[-filter,]
  }
  
  
  if (test.set=='GSE116918_Metastasis') {
    phenoData$time_to_bcr <- phenoData$time_to_metastasis
    phenoData$bcr_status <- phenoData$metastasis_status
  }
  
  filter <- which(is.na(phenoData$bcr_status) | is.na(phenoData$time_to_bcr))
  
  if (length(filter)>0) {
    exprData <- exprData[,-filter]
    phenoData <- phenoData[-filter,]
  }
  
  ##### Smples used for classification
  ### ========================================= ###
  # phenoData$bcr_status_5yr <- getEventFun(n=5, time.to.event=phenoData$time_to_bcr, status=phenoData$bcr_status)
  # 
  # table(phenoData$bcr_status_5yr)
  # 
  # filter <- which(is.na(phenoData$bcr_status_5yr))
  # phenoData <- phenoData[-filter,]
  # exprData <- exprData[,-filter]
  ### ========================================= ###
  
  
  test.pheno <- phenoData
  dim(test.pheno)
  test.pheno$time_to_bcr[test.pheno$time_to_bcr<=0] <- 0.01
  
  test.geno.all <- as.matrix(t(exprData))
  test.geno.all <- scale(test.geno.all)
  
  
  for (signature.name in signatures) {
    message(signature.name)
    
    signature <- read_xlsx(path = 'data/Classifiers.xlsx', sheet=signature.name)
    signature.genes <- signature$Ensembl
    
    genes <- Reduce(intersect, list(signature.genes, colnames(training.geno.all), colnames(test.geno.all)))
    #genes <- Reduce(intersect, list(genes, colnames(training.geno.all), colnames(test.geno.all)))
    training.geno <- training.geno.all[,genes]
    test.geno <- test.geno.all[,genes]
    
    
    ################################ CoxPH

    model <- 'CoxPH'
    print (model)

    training.bcr.time <- training.pheno$time_to_bcr
    training.bcr.status <-training.pheno$bcr_status

    test.bcr.time <- test.pheno$time_to_bcr
    test.bcr.status <-  test.pheno$bcr_status

    training.surv.data <- data.frame(training.geno, bcr.time=training.bcr.time, bcr.status=training.bcr.status)

    multi.var <- paste0(genes, collapse = '+')

    fml <- as.formula(paste0('Surv(bcr.time, bcr.status) ~ ', multi.var))

    coxtest <- coxph(fml, data = training.surv.data)
    summcph <- summary(coxtest)
    coeffs <- summcph$coefficients[,1]



    training.risk.score <- predict(coxtest, data.frame(training.geno), type="lp", se.fit=FALSE) # type='risk'
    training.risk.score

    training.risk.threshold <- median(training.risk.score, na.rm = T)
    training.risk.threshold

    test.risk.score <- predict(coxtest, data.frame(test.geno), type="lp", se.fit=FALSE)
    test.risk.score

    test.risk.threshold <- median(test.risk.score, na.rm = T)
    test.risk.threshold

    training.risk.group <- test.risk.score > training.risk.threshold
    training.risk.group

    test.risk.group <- test.risk.score > test.risk.threshold
    test.risk.group

    riskScore <- test.risk.score
    trainingRiskGroup <- training.risk.group
    testRiskGroup <- test.risk.group

    riskThreshold <- median(riskScore, na.rm=T)
    riskGroup <- riskScore > riskThreshold

    res <- data.frame(risk.score=riskScore, risk.group=riskGroup, training.risk.group=trainingRiskGroup,
                      test.risk.group=testRiskGroup, bcr.time=test.pheno$time_to_bcr, bcr.status=test.pheno$bcr_status)
    res

    write.csv(res, file=paste0('report/Survival/TrainTestScale/Signature/Survival_',
                               model, '_', signature.name, '_', training.set, '_', test.set, '_TrainTest_Signature.csv'),
              quote = F)

    ##### Smples used for classification
    # ### ========================================= ###
    # write.csv(res, file=paste0('report/SurvivalC/TrainTestScale/Signature/Survival_',
    #                            model, '_', signature.name, '_', training.set, '_', test.set, '_TrainTest_Signature.csv'),
    #           quote = F)
    # ### ========================================= ###
    
    
    
    ################################ CoxNet
    
    ### LASSO by default, alpha = 1
    
    alpha <- 1
    model <- 'CoxNet'
    
    #alpha <- 0
    #model <- paste0('CoxNetAlpha',alpha)


    training.bcr.time <- training.pheno$time_to_bcr
    training.bcr.status <-training.pheno$bcr_status

    test.bcr.time <- test.pheno$time_to_bcr
    test.bcr.status <-  test.pheno$bcr_status

    training.surv.data <- data.frame(training.geno, bcr.time=training.bcr.time, bcr.status=training.bcr.status)

    set.seed(777)
    cv.fit <- cv.glmnet(training.geno, Surv(training.bcr.time, training.bcr.status), family="cox", maxit = 1000,
                        alpha=alpha)

    coeffs <- coef(cv.fit, s = cv.fit$lambda.min)
    coeffs <- as.numeric(coeffs)

    message(paste0(sum(coeffs!=0),'/',length(genes)))

    training.risk.score <- predict(cv.fit, s=cv.fit$lambda.min, newx = training.geno, type="link") # response: reltive risk; link: linear prediction
    training.risk.score
    training.risk.score <- training.risk.score[,1]

    training.risk.threshold <- median(training.risk.score, na.rm = T)
    training.risk.threshold

    test.risk.score <- predict(cv.fit, s=cv.fit$lambda.min, newx = test.geno, type="link") # response: reltive risk; link: linear prediction
    test.risk.score
    test.risk.score <- test.risk.score[,1]

    test.risk.threshold <- median(test.risk.score, na.rm = T)
    test.risk.threshold

    training.risk.group <- test.risk.score > training.risk.threshold
    training.risk.group

    test.risk.group <- test.risk.score > test.risk.threshold
    test.risk.group

    riskScore <- test.risk.score
    trainingRiskGroup <- training.risk.group
    testRiskGroup <- test.risk.group

    riskThreshold <- median(riskScore, na.rm=T)
    riskGroup <- riskScore > riskThreshold

    res <- data.frame(risk.score=riskScore, risk.group=riskGroup, training.risk.group=trainingRiskGroup,
                      test.risk.group=testRiskGroup, bcr.time=test.pheno$time_to_bcr, bcr.status=test.pheno$bcr_status)
    res

    write.csv(res, file=paste0('report/Survival/TrainTestScale/Signature/Survival_',
                              model, '_', signature.name, '_', training.set, '_', test.set, '_TrainTest_Signature.csv'),
             quote = F)

    ##### Smples used for classification
    # ### ========================================= ###
    # write.csv(res, file=paste0('report/SurvivalC/TrainTestScale/Signature/Survival_',
    #                            model, '_', signature.name, '_', training.set, '_', test.set, '_TrainTest_Signature.csv'),
    #           quote = F)
    # ### ========================================= ###
    
    
    
    ##################################### superpc
    
    threshold <- 0.3
    model <- 'SuperPC'
    print (model)
    
    # threshold <- 0 # 0.1
    # model <- paste0('SuperPC',threshold)

    training.bcr.time <- training.pheno$time_to_bcr
    training.bcr.status <-training.pheno$bcr_status

    test.bcr.time <- test.pheno$time_to_bcr
    test.bcr.status <-  test.pheno$bcr_status

    #training.surv.data <- data.frame(training.geno, bcr.time=training.bcr.time, bcr.status=training.bcr.status)

    training.data.pca <-list(x=t(training.geno), y=training.bcr.time, censoring.status= training.bcr.status, featurenames=genes)
    test.data.pca <- list(x=t(test.geno), y=test.bcr.time, censoring.status= test.bcr.status, featurenames=genes)

    pca.fit <- superpc.train(training.data.pca, type="survival")
    
    message(paste0(sum(abs(pca.fit$feature.scores)>threshold),'/',length(genes)))

    training.risk.score <- superpc.predict(pca.fit, training.data.pca, training.data.pca,
                                           threshold=threshold, n.components=1, prediction.type="continuous")

    training.risk.score <- training.risk.score$v.pred[,1]
    names(training.risk.score) <- rownames(training.geno)

    training.risk.threshold <- median(training.risk.score, na.rm = T)

    test.risk.score <- superpc.predict(pca.fit, training.data.pca, test.data.pca,
                                       threshold=threshold, n.components=1, prediction.type="continuous")

    test.risk.score <- test.risk.score$v.pred[,1]
    names(test.risk.score) <- rownames(test.geno)

    test.risk.threshold <- median(test.risk.score, na.rm = T)
    test.risk.threshold

    training.risk.group <- test.risk.score > training.risk.threshold
    training.risk.group

    test.risk.group <- test.risk.score > test.risk.threshold
    test.risk.group

    riskScore <- test.risk.score
    trainingRiskGroup <- training.risk.group
    testRiskGroup <- test.risk.group

    riskThreshold <- median(riskScore, na.rm=T)
    riskGroup <- riskScore > riskThreshold

    res <- data.frame(risk.score=riskScore, risk.group=riskGroup, training.risk.group=trainingRiskGroup,
                      test.risk.group=testRiskGroup, bcr.time=test.pheno$time_to_bcr, bcr.status=test.pheno$bcr_status)
    res

    write.csv(res, file=paste0('report/Survival/TrainTestScale/Signature/Survival_',
                              model, '_', signature.name, '_', training.set, '_', test.set, '_TrainTest_Signature.csv'),
             quote = F)

    ##### Smples used for classification
    # ### ========================================= ###
    # write.csv(res, file=paste0('report/SurvivalC/TrainTestScale/Signature/Survival_',
    #                            model, '_', signature.name, '_', training.set, '_', test.set, '_TrainTest_Signature.csv'),
    #           quote = F)
    # ### ========================================= ###
    
  }
}