#folder = "Add the source code folder here" folder = "C:/Users/Shennor/Dropbox/research_documents/Postdoc/Shai/multipleAlignment/Figures/New/Revision/V2/sourceCode" source(file.path(folder, "enrichmentTest.R")) library(TimeAx) library(ggplot2) dataPlotter = function(currMetaData, currTraj, currMeta, additionalCol = NULL, box = T, yValues = NULL, xValues = NULL, regLine = F, DoF = 4){ x = currMetaData[,currMeta] indsToUse = which(!is.na(x) & !is.na(currTraj)) metaDataLevels = x[indsToUse] currTrajLevels = currTraj[indsToUse] currCorValue = cor(metaDataLevels, currTrajLevels) if(length(unique(metaDataLevels))<6 & box){ if(length(unique(metaDataLevels)) == 2 & min(table(metaDataLevels))>3){ currCorValue = t.test(currTrajLevels[metaDataLevels == unique(metaDataLevels)[1]], currTrajLevels[metaDataLevels == unique(metaDataLevels)[2]])$p.value } ggplot(data = NULL, aes(x = factor(metaDataLevels), y = currTrajLevels))+geom_boxplot()+ggtitle(paste(currMeta, currCorValue,sep = " _ "))+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) }else{ if(is.null(xValues)){ xValues = c(min(currTrajLevels),max(currTrajLevels)) } if(is.null(yValues)){ yValues = c(min(metaDataLevels),max(metaDataLevels)) } regLineFunc = NULL if(regLine){ regLineFunc = geom_smooth(method = "lm", formula = y ~ poly(x, DoF), se = F) } if(!is.null(additionalCol)){ additionalCol = factor(additionalCol[indsToUse]) ggplot(data = NULL, aes(x = currTrajLevels, y = metaDataLevels))+geom_point(size = 3, aes(color = additionalCol))+ggtitle(paste(currMeta, currCorValue,sep = " _ "))+xlim(xValues)+ylim(yValues)+regLineFunc+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) }else{ ggplot(data = NULL, aes(x = currTrajLevels, y = metaDataLevels))+geom_point(size = 3)+ggtitle(paste(currMeta, currCorValue,sep = " _ "))+xlim(xValues)+ylim(yValues)+regLineFunc+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) } } } proteinCodingGenes = as.character(as.matrix(read.table(file.path(folder,"proteinCodingGenes.txt"),sep = "\t", row.names = NULL, col.names = F))) #### Load and norm data #### # Main cohort # UBCFullData = readRDS(file.path(folder,"UBCFullDataRaw.rds")) annotations = UBCFullData$annotations GEData = UBCFullData$GEData annotationDataFull = read.table(file.path(folder, "annotationTableFull.txt"),sep = "\t", row.names = 1, header = T, check.names = F)[colnames(GEData),] annotationDataFull = annotationDataFull[colnames(GEData),] sampleNames = annotations$Patient_ID trajectory = annotations$Tumor_Number batchInfo = read.table(file.path(folder,"batchInformation.txt"),sep="\t", row.names = 1, header = T, check.names = F) batchInfo = batchInfo[colnames(GEData),] GEDataCombat = sva::ComBat(GEData, batchInfo) trajectoryFixed = unlist(lapply(unique(sampleNames),function(currName){ 1:length(which(currName == sampleNames))})) trajectoryTime = as.numeric(annotations$Months_from_Primary)/12 GEDataToUse = GEDataCombat trainSamples = names(which(table(sampleNames)>3)) testSamples = unique(sampleNames)[!(unique(sampleNames) %in% trainSamples)] trainIndexes = which(sampleNames %in% trainSamples) testIndexes = which(sampleNames %in% testSamples) sampleNamesTrain = sampleNames[trainIndexes] sampleNamesTest = sampleNames[testSamples] GEDataTrain = GEDataToUse[,trainIndexes] GEDataTest = GEDataToUse[,testSamples] trajectoryFixedTrain = trajectoryFixed[trainIndexes] trajectoryFixedTest = trajectoryFixed[testSamples] trajectoryTimeTrain = trajectoryTime[trainIndexes] trajectoryTimeTest = trajectoryTime[testSamples] annotationsTrain = annotations[trainIndexes,] annotationsTest = annotations[testSamples,] metaDataNumTrain = apply(annotationsTrain,2,function(x){ xFixed = as.numeric(x) if(length(which(is.na(xFixed)))>30){ xFixed = as.numeric(as.factor(x)) } xFixed[x=="N/A"] = NA xFixed }) metaDataNumTest = apply(annotationsTest,2,function(x){ xFixed = as.numeric(x) if(length(which(is.na(xFixed)))>30){ xFixed = as.numeric(as.factor(x)) } xFixed[x=="N/A"] = NA xFixed }) # Microarray cohort # dataGSE83586Full = readRDS(file.path(folder,"UBCFullDataRawGSE83586.rds")) dataGSE83586 = dataGSE83586Full$GEData annotationsGSE83586 = dataGSE83586Full$annotations batchInfoGSE83586 = read.table(file.path(folder,"GSE83586_batches.txt"),sep="\t", row.names = 1, header = T, check.names = F) batchInfoGSE83586 = batchInfoGSE83586[colnames(dataGSE83586),] dataGSE83586 = dataGSE83586[intersect(row.names(dataGSE83586),row.names(GEDataToUse)),] dataGSE83586 = sva::ComBat(dataGSE83586, batchInfoGSE83586) # TCGA cohort # dataFullTCGA = readRDS(file.path(folder,"TCGAData.rds")) dataTCGA = dataFullTCGA$data colnames(dataTCGA) = sapply(colnames(dataTCGA),function(x){unlist(strsplit(x,"\\."))[1]}) dataTCGA = dataTCGA[which(!is.na(rowSums(dataTCGA))),] annotationsTCGA = dataFullTCGA$clinical dataTCGANorm = as.matrix(log(dataTCGA+1)) dataTCGANorm = dataTCGANorm[which(!is.nan(rowSums(dataTCGANorm))),] indexesTCGAOneSample = match(names(which(table(annotationsTCGA$case_submitter_id)==1)),annotationsTCGA$case_submitter_id) dataTCGANorm = dataTCGANorm[,indexesTCGAOneSample] annotationsTCGA = annotationsTCGA[indexesTCGAOneSample,] #### Calculate model and pseudotime #### proteinCodingGenes = intersect(proteinCodingGenes, row.names(GEDataTrain)) multiAlignment = modelCreation(GEDataTrain[proteinCodingGenes,], sampleNamesTrain, numOfTopFeatures = 100) seedGenes = multiAlignment$seed metPoint = 0.7 trajToUseProfileTrain = predictByConsensus(multiAlignment,GEDataTrain)$predictions trajToUseProfile = predictByConsensus(multiAlignment,GEDataTest)$predictions annotationsToUseProfile = annotationsTest trajToUseProfileGSE83586 = predictByConsensus(multiAlignment,dataGSE83586)$predictions trajToUseProfileTCGA = predictByConsensus(multiAlignment,dataTCGANorm)$predictions ##### Figure 2 ##### #### 2E ##### trajectoryTimeTrainForVis = trajectoryTimeTrain trajectoryTimeTrainForVis[which.max(trajectoryTimeTrainForVis)] = 8 chosenPatients = c("20","26","30","40","57","65","67","72") tmpData = as.data.frame(cbind(trajectoryTimeTrainForVis,trajToUseProfileTrain,sampleNamesTrain)) colnames(tmpData) = c("Pre","Post","Sample") tmpData = tmpData[tmpData$Sample %in% chosenPatients,] ggplot(data = tmpData, aes(x = Pre,y = Post, colour = factor(Sample)))+geom_line(size = 1)+ylim(0,1)+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### 2F #### regPre = summary(lm(t(GEDataTrain)~poly(trajectoryTimeTrain,2))) genePValuesPre = sapply(regPre, function(x){ fstat <- x$fstatistic pf(fstat[1], fstat[2], fstat[3], lower.tail=FALSE) #x[[4]][2,4] }) regAligned = summary(lm(t(GEDataTrain)~poly(trajToUseProfileTrain,2))) genePValuesAligned = sapply(regAligned, function(x){ fstat <- x$fstatistic pf(fstat[1], fstat[2], fstat[3], lower.tail=FALSE) #x[[4]][2,4] }) geneRegPreLog = -log(genePValuesPre,base=10) geneRegAlignedLog = -log(genePValuesAligned,base=10) newCutOff = -log(0.01, base=10) geneColourProVsAligned = rep(0, length(geneRegPreLog)) geneColourProVsAligned[geneRegPreLog>newCutOff & geneRegAlignedLognewCutOff] = 2 geneColourProVsAligned[geneRegPreLog>newCutOff & geneRegAlignedLog>newCutOff] = 3 ggplot(data = NULL, aes(x = geneRegPreLog, y = geneRegAlignedLog, colour = factor(geneColourProVsAligned)))+geom_point()+ geom_hline(yintercept = newCutOff, linetype = "dashed")+geom_vline(xintercept = newCutOff, linetype = "dashed")+ scale_color_manual(values = c("0" = "lightgray","1" = "darkgray","2" = "purple", "3" = "orange"))+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### S3A #### factorToPlot = "CCL2" #factorToPlot = "IFITM2" #factorToPlot = "SGPL1" dataPlotter(t(GEDataTrain), trajToUseProfileTrain,factorToPlot,regLine = T, DoF=3) dataPlotter(t(GEDataTrain), trajectoryTimeTrain,factorToPlot,regLine = T, DoF=3) #### S3B #### geneSets = readRDS(file.path(folder, "geneSets.rds")) corGeneTrajTime = cor(t(GEDataTrain), trajectoryTimeTrain)[,1] corGeneTrajPseudo = cor(t(GEDataTrain), trajToUseProfileTrain)[,1] enrichmentTime = runEnrichment(sort(corGeneTrajTime),geneSets) enrichmentPseudo = runEnrichment(sort(corGeneTrajPseudo),geneSets) minPoint = min(enrichmentTime,enrichmentPseudo) maxPoint = max(enrichmentTime,enrichmentPseudo) ggplot(data = NULL, aes(x = enrichmentTime, y = enrichmentPseudo))+geom_point()+ xlim(minPoint,maxPoint)+ylim(minPoint,maxPoint)+geom_hline(yintercept = 0, linetype = "dashed")+geom_vline(xintercept = 0, linetype = "dashed")+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) ##### Figure 3 ##### #### 3A #### stagesFull = annotations$Stage stagesFull[stagesFull == "N/A" | stagesFull == "cis" | stagesFull == "Met" | stagesFull == "T1M1"] = "Unknown" stagesFull[stagesFull == "Ta"] = "T0" stagesFull[grep("T3",stagesFull)] = "T3" stages = stagesFull[trainIndexes] stageCont = as.numeric(gsub("T","",stages[stages!="Unknown"])) tStage = summary(lm(stageCont~trajToUseProfileTrain[stages!="Unknown"]))[[4]][2,4] ggplot(data = NULL, aes(x = factor(stages[stages!="Unknown"]), y = trajToUseProfileTrain[stages!="Unknown"]))+geom_boxplot()+ggtitle(paste("Main",tStage,sep=" _ "))+geom_jitter(height = 0, width = 0.1)+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### S5A #### tTime = summary(lm(stageCont~trajectoryTimeTrain[stages!="Unknown"]))[[4]][2,4] ggplot(data = NULL, aes(x = factor(stages[stages!="Unknown"]), y = trajectoryTimeTrain[stages!="Unknown"]))+geom_boxplot()+ggtitle(tTime)+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### S5B #### stagesTest = annotationsToUseProfile$Stage stagesTest[stagesTest == "N/A" | stagesTest == "cis" | stagesTest == "Met" | stagesTest == "T1M1"] = "Unknown" stagesTest[stagesTest == "Ta"] = "T0" stagesTest[grep("T3",stagesTest)] = "T3" stageContTest = as.numeric(gsub("T","",stagesTest[stagesTest!="Unknown"])) tStageTest = summary(lm(stageContTest~trajToUseProfile[stagesTest!="Unknown"]))[[4]][2,4] ggplot(data = NULL, aes(x = factor(stagesTest[stagesTest!="Unknown"]), y = trajToUseProfile[stagesTest!="Unknown"]))+geom_boxplot()+ggtitle(paste("Main",tStageTest,sep=" _ "))+geom_jitter(height = 0, width = 0.1)+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) stagesTCGA = annotationsTCGA$ajcc_pathologic_t stagesTCGA[stagesTCGA == "--"] = "Unknown" stagesTCGA[stagesTCGA == ""] = "Unknown" stagesTCGA[stagesTCGA == "TX"] = "Unknown" stagesTCGA[grep("T1",stagesTCGA)] = "T1" stagesTCGA[grep("T2",stagesTCGA)] = "T2" stagesTCGA[grep("T3",stagesTCGA)] = "T3" stagesTCGA[grep("T4",stagesTCGA)] = "T4" stageTCGACont = as.numeric(gsub("T","",stagesTCGA[stagesTCGA!="Unknown"])) tStageTCGA = summary(lm(stageTCGACont~trajToUseProfileTCGA[stagesTCGA!="Unknown"]))[[4]][2,4] ggplot(data = NULL, aes(x = factor(stagesTCGA[stagesTCGA!="Unknown"]), y = trajToUseProfileTCGA[stagesTCGA!="Unknown"]))+geom_boxplot()+geom_jitter(height = 0, width = 0.1)+ggtitle(paste("TCGA",tStageTCGA,sep=" _ "))+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) stagesGSE83586 = annotationsGSE83586$Pathologically_Reviewed_TURB_Sample_Stage stagesGSE83586[stagesGSE83586 == "pTa" | stagesGSE83586 == "pTis"] = "Stage0" stagesGSE83586[stagesGSE83586 == "pT1"] = "Stage1" stagesGSE83586[stagesGSE83586 == "pT2"] = "Stage2" stagesGSE83586[stagesGSE83586 == "pT3"] = "Stage3" stagesGSE83586[stagesGSE83586 == "pT4"] = "Stage4" stagesGSE83586[stagesGSE83586 == "pTx"] = "Unknown" stageGSE83586Cont = as.numeric(gsub("Stage","",stagesGSE83586[stagesGSE83586!="Unknown"])) tStageGSE83586 = summary(lm(stageGSE83586Cont~trajToUseProfileGSE83586[stagesGSE83586!="Unknown"]))[[4]][2,4] ggplot(data = NULL, aes(x = factor(stagesGSE83586[stagesGSE83586!="Unknown"]), y = trajToUseProfileGSE83586[stagesGSE83586!="Unknown"]))+geom_boxplot()+geom_jitter(height = 0, width = 0.1)+ggtitle(paste("GSE83586",tStageGSE83586,sep=" _ "))+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### S5C #### clinicalFullTCGA = read.table(file.path(folder,"clinicalFull.tsv"),sep="\t",row.names = 1,header = T,check.names = F) clinicalFullTCGA = clinicalFullTCGA[intersect(annotationsTCGA$case_id, row.names(clinicalFullTCGA)),] clinicalTCGAIndexes = which(!is.na(match(annotationsTCGA$case_submitter_id,clinicalFullTCGA$case_submitter_id))) trajForTCGAGroups = trajToUseProfileTCGA[clinicalTCGAIndexes] cancerTypeTCGA = clinicalFullTCGA$primary_diagnosis cancerTypeTCGA[grep("transitional",tolower(cancerTypeTCGA),invert = T)] = "Unknown" tTypeTCGA = t.test(trajForTCGAGroups[cancerTypeTCGA=="Papillary transitional cell carcinoma"], trajForTCGAGroups[cancerTypeTCGA=="Transitional cell carcinoma"],alternative = "less")$p.value ggplot(data = NULL, aes(x = factor(cancerTypeTCGA[cancerTypeTCGA!="Unknown"]), y = trajForTCGAGroups[cancerTypeTCGA!="Unknown"]))+geom_boxplot()+geom_jitter(height = 0, width = 0.1)+ggtitle(paste("TCGA",tTypeTCGA,sep=" _ "))+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### S5D #### cancerPriorTCGA = clinicalFullTCGA$prior_malignancy tPriorTCGA = t.test(trajForTCGAGroups[cancerPriorTCGA=="no"],trajForTCGAGroups[cancerPriorTCGA=="yes"],alternative = "less")$p.value ggplot(data = NULL, aes(x = factor(cancerPriorTCGA), y = trajForTCGAGroups))+geom_boxplot()+geom_jitter(height = 0, width = 0.1)+ggtitle(paste("TCGA",tPriorTCGA,sep=" _ "))+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### 3B #### factorToPlot = "ESTIMATE_purity" dataPlotter(metaDataNumTrain, trajToUseProfileTrain,factorToPlot,regLine = T, DoF=3) dataPlotter(metaDataNumTrain, trajectoryTimeTrain,factorToPlot,regLine = T, DoF=3) #### S5E-F #### ageTCGA = annotationsTCGA$age_at_index beforeAfterPointTCGA = rep(-1, length(trajToUseProfileTCGA)) beforeAfterPointTCGA[trajToUseProfileTCGA>metPoint] = 1 agePTCGA = t.test(ageTCGA[beforeAfterPointTCGA==-1],ageTCGA[beforeAfterPointTCGA==1])$p.value ggplot(data = NULL, aes(x = factor(beforeAfterPointTCGA), y = ageTCGA))+geom_boxplot()+geom_jitter(height = 0, width = 0.1)+ggtitle(agePTCGA)+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) sexTableTCGA = table(annotationsTCGA$gender,beforeAfterPointTCGA) fisher.test(sexTableTCGA) beforeAfterPoint = rep(-1, length(trajToUseProfileTrain)) beforeAfterPoint[trajToUseProfileTrain>metPoint] = 1 sexTableTrain = table(annotationsTrain$Sex,beforeAfterPoint) fisher.test(sexTableTrain) fullSexTable = rbind(t(t(sexTableTCGA)/colSums(sexTableTCGA)),t(t(sexTableTrain)/colSums(sexTableTrain))) row.names(fullSexTable) = c("F-TCGA","M-TCGA","F-Long","M-Long") sexTableDF = as.data.frame(reshape2::melt(fullSexTable)) colnames(sexTableDF) = c("General","PrePost","Ratio") sexTableDF$Sex = sapply(as.character(sexTableDF$General),function(x){unlist(strsplit(x,"-"))[1]}) sexTableDF$Data = sapply(as.character(sexTableDF$General),function(x){unlist(strsplit(x,"-"))[2]}) ggplot(data = sexTableDF, aes(x = interaction(PrePost,Data), y = Ratio, fill = factor(Sex)))+geom_bar(stat="identity")+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### 3C #### source(file.path(folder,"CIBERSORT_LibLinear.R")) LM22 = read.table(file.path(folder, "LM22Matrix.txt"),sep="\t", header = T, row.names = 1, check.names = F) LM22 = LM22[intersect(row.names(LM22),row.names(GEDataTrain)),] singleCellData = readRDS(file.path(folder, "GSM4307111_SingleCellData.rds")) singleCellInfo = read.table(file.path(folder, "GSM4307111_CellInfo.txt"),sep="\t", header = T, row.names = 1, check.names = F) cellTypes = singleCellInfo[colnames(singleCellData),"cell_type"] selectedCells = grep("Unknown",cellTypes,invert = T) singleCellData = singleCellData[,selectedCells] cellTypes = cellTypes[selectedCells] geneMarkers = unique(unlist(lapply(unique(cellTypes), function(currCellType){ currCells = singleCellData[,cellTypes==currCellType] otherCells = singleCellData[,cellTypes!=currCellType] names(sort(rowMeans(currCells)-rowMeans(otherCells),decreasing = T)[1:20]) }))) signatureMatrix = do.call(cbind,lapply(unique(cellTypes), function(currCellType){ rowMeans(singleCellData[geneMarkers,cellTypes==currCellType]) })) row.names(signatureMatrix) = geneMarkers colnames(signatureMatrix) = unique(cellTypes) geneMarkersMain = intersect(geneMarkers, row.names(GEDataTrain)) cibersortRes = CIBERSORT(signatureMatrix[geneMarkersMain,],GEDataTrain[geneMarkersMain,]) cibersortResLM22 = CIBERSORT(LM22,GEDataTrain[geneMarkersMain,]) cibersortResLM22 = cbind(cibersortResLM22,cibersortResLM22[,"T cells CD4 memory resting"]+cibersortResLM22[,"T cells CD4 memory activated"]) colnames(cibersortResLM22)[dim(cibersortResLM22)[2]] = "T cells CD4 memory" ciberCorWithTraj = cor(cibersortRes,trajToUseProfileTrain)[,1] ciberCorWithTrajLM22 = cor(cibersortResLM22,trajToUseProfileTrain)[,1] cellTypeOrder = names(sort(ciberCorWithTraj)) cellTypeOrderLM22 = names(sort(ciberCorWithTrajLM22)) cellTypeDeconPlots = lapply(cellTypeOrder, function(currCellType){ currDeconLevels = cibersortRes[,currCellType] trajDisc = rep(1, length(currDeconLevels)) trajDisc[trajToUseProfileTrainmetPoint) stageGroups = as.numeric(gsub("T","",stagesTCGA[clinicalTCGAIndexes])) classGroups = gsub("-Inf","",TCGAClassifier[clinicalTCGAIndexes,]$LundTax_subclasses) TCGAAgeAtIndex = clinicalFullTCGA$age_at_index TCGAGender = clinicalFullTCGA$gender survDF = data.frame(time = survTime, vital = survGroup, pseudo = trajForTCGAGroups, inflection = pseudoGroups, stage = stageGroups, age = TCGAAgeAtIndex, sex = TCGAGender) survDF = survDF[survTime>-Inf,] survFit <- survfit(Surv(time, vital) ~ inflection, data = survDF) surv_pvalue(survFit) ggsurvplot(survFit, data = survDF, risk.table = F) #### 3E #### TCGAClassifier = read.table(file.path(folder,"LundTax_and TCGA_subtypes_TCGA_data.txt"), sep="\t", row.names=1, header = T, check.names = F) TCGAClassifier = TCGAClassifier[annotationsTCGA$case_submitter_id,] uroASamplesTCGA = which(TCGAClassifier$LundTax_subclasses=="Basal_squamous") trajUroASamplesTCGA = trajToUseProfileTCGA[uroASamplesTCGA] annotationsTCGAUroA = annotationsTCGA[uroASamplesTCGA,] highUro = as.numeric(trajUroASamplesTCGA>metPoint) survTable = table(annotationsTCGAUroA$vital_status,highUro) survTableP = fisher.test(survTable)$p.value survDFFull = as.data.frame(reshape2::melt(t(t(survTable)/colSums(survTable)))) ggplot(data = survDFFull, aes(x = factor(highUro), y = value, fill = factor(Var1)))+geom_col()+ggtitle(survTableP)+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### S5I-J #### classForSurvival = "Urothelial-like" #classForSurvival = "Basal_squamous" specificIndexUro = which(TCGAClassifier[clinicalTCGAIndexes,]$LundTax_classes==classForSurvival) survFitSpecific <- survfit(Surv(time, vital) ~ inflection, data = survDF[specificIndexUro,]) surv_pvalue(survFitSpecific) ggsurvplot(survFitSpecific, data = survDF[specificIndexUro,], risk.table = F) ##### Figure 4 ##### #### S6A #### tumorSubType = annotationsTrain$Subtype_LundTax_RNA patientsNoBasal = unique(sampleNamesTrain)[sapply(unique(sampleNamesTrain),function(currPatient){ length(which(c("Ba/Sq") %in% tumorSubType[sampleNamesTrain == currPatient]))==0 })] GEDataTrainNoBasal = GEDataTrain[,sampleNamesTrain %in% patientsNoBasal] sampleNamesTrainNoBasal = sampleNamesTrain[sampleNamesTrain %in% patientsNoBasal] multiAlignmentNoBasal = modelCreation(GEDataTrainNoBasal[proteinCodingGenes,], sampleNamesTrainNoBasal, numOfTopFeatures = 100) trajToUseProfileTrainNoBasal = predictByConsensus(multiAlignmentNoBasal,GEDataTrain)$predictions noBasalOverlap = cor(trajToUseProfileTrain,trajToUseProfileTrainNoBasal) ggplot(data = NULL, aes(x = trajToUseProfileTrain,y = trajToUseProfileTrainNoBasal))+geom_point()+ggtitle(noBasalOverlap)+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### S6B #### regAlignedNoBasal = summary(lm(t(GEDataTrain)~trajToUseProfileTrainNoBasal)) genePValuesAlignedNoBasal = sapply(regAlignedNoBasal, function(x){ x[[4]][2,4] }) geneRegAlignedNoBasalLog = -log(genePValuesAlignedNoBasal,base=10) geneColourNoBasalVsAligned = rep(0, length(geneRegAlignedNoBasalLog)) geneColourNoBasalVsAligned[geneRegAlignedNoBasalLog>newCutOff & geneRegAlignedLognewCutOff] = 2 geneColourNoBasalVsAligned[geneRegAlignedNoBasalLog>newCutOff & geneRegAlignedLog>newCutOff] = 3 ggplot(data = NULL, aes(x = geneRegAlignedNoBasalLog, y = geneRegAlignedLog, colour = factor(geneColourNoBasalVsAligned)))+geom_point()+ geom_hline(yintercept = newCutOff, linetype = "dashed")+geom_vline(xintercept = newCutOff, linetype = "dashed")+ scale_color_manual(values = c("0" = "lightgray","1" = "darkgray","2" = "purple", "3" = "orange"))+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### 4A #### lundTrain = annotationsTrain$Subtype_LundTax_RNA tLund =aov(trajToUseProfileTrain~lundTrain) summary(tLund)[[1]][["Pr(>F)"]][1] ggplot(data = NULL, aes(x = reorder(factor(lundTrain),trajToUseProfileTrain), y = trajToUseProfileTrain))+geom_boxplot()+geom_jitter(height = 0, width = 0.1)+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### 4B #### tumorSubType = annotationsTrain$Subtype_LundTax_RNA inflectionVector = rep("Before",length(tumorSubType)) inflectionVector[trajToUseProfileTrain>metPoint] = "After" inflectionTypeTable = table(tumorSubType,inflectionVector) inflectionTypeTable = t(t(inflectionTypeTable)/colSums(inflectionTypeTable)) inflectionDF = as.data.frame(reshape2::melt(inflectionTypeTable)) ggplot(data = inflectionDF, aes(x = reorder(factor(tumorSubType),value,na.rm = TRUE), y = value, fill = factor(inflectionVector)))+ geom_bar(stat="identity", color="black", position=position_dodge())+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### 4C #### tumorSubTypeTrain = gsub(" ", "", annotationsTrain$Prim_Rec_Prog) indexesForFigure4 = which(trajToUseProfileTrain(i-0.1) & survDFForBins$pseudo0.25) survTableUroEarlyLate = table(survDFUroPre$vital,survDFUroPre$progSetClass) survDFEarlyLate = as.data.frame(reshape2::melt(t(t(survTableUroEarlyLate)/colSums(survTableUroEarlyLate)))) survTablePEarlyLate = fisher.test(survTableUroEarlyLate)$p.value ggplot(data = survDFEarlyLate, aes(x = factor(Var2), y = value, fill = factor(Var1)))+geom_col()+ggtitle(survTablePEarlyLate)+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### S6G #### uroASamples = which(annotationsTrain$Subtype_LundTax_RNA=="UroA-Prog" & trajToUseProfileTrain0.25])$p.value }) uroAGenePAllQ = p.adjust(uroAGenePAll, method = "fdr") #signGeneUroA = names(uroAGenePAllQ[uroAGenePAll<0.05]) signGeneUroA = names(uroAGenePAllQ[uroAGenePAllQ<0.05]) signGeneUroA = signGeneUroA[!(signGeneUroA %in% seedGenes)] signGenesCorMatrix = cor(t(GEDataUroAllGenes[signGeneUroA,])) createGeneHeatmap = function(coExp, prvStats = NULL){ library(WGCNA) if(is.null(prvStats)){ dissTOM=1-coExp hierADJ=hclust(as.dist(dissTOM), method="average") currOrder = hierADJ$order }else{ genesToCheck = prvStats$genes currOrder = prvStats$order coExp = coExp[genesToCheck,genesToCheck] } collectGarbage() pd.m <- reshape2::melt(coExp[currOrder,currOrder], id.vars = "Gene1", variable.name = "Gene2") currPlot = ggplot( pd.m, aes(Var1, Var2) ) + geom_tile(aes(fill = value)) + scale_fill_gradient2(low="blue", mid = 'white',midpoint = 0, high="red")+ theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())+ theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) list(plot = currPlot, order = currOrder) } signGenesHeatmap = createGeneHeatmap(signGenesCorMatrix) #### S6H-I #### firstGroup = 1:1055 secondGroup = (max(firstGroup)+1):dim(signGenesCorMatrix)[1] smallGroup = row.names(signGenesCorMatrix)[signGenesHeatmap$order][firstGroup] bigGroup = row.names(signGenesCorMatrix)[signGenesHeatmap$order][secondGroup] #geneTotalScore = colMeans(GEDataTrain[bigGroup,]) geneTotalScore = colMeans(GEDataTrain[smallGroup,]) ggplot(data = NULL, aes(x = trajToUseProfileTrain[uroASamples], y = geneTotalScore[uroASamples]))+geom_point(size = 4)+geom_smooth()+geom_vline(xintercept = metPoint,linetype="dashed")+ theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) #### 5B #### lcRNATypes = c("^RNU","^MIR","^RN7SK","^RN7SL","^RNA5S") pseudoEnrichment = function(groups,groupForEnrich, PW = F, groupCutOff = 0){ sapply(groups,function(currGroup){ if(PW){ groupGenes = currGroup }else{ groupGenes = row.names(GEDataTrain)[grep(currGroup,row.names(GEDataTrain))] } groupSelected = length(intersect(groupForEnrich,groupGenes)) if(groupSelected0.25) geneDataForUroDF = as.data.frame(reshape2::melt(geneDataForUroNorm[kinetochoreGenes,])) uroGenesPs = apply(geneDataForUroNorm[kinetochoreGenes,],1,function(x){ t.test(x[colnames(geneDataForUroNorm)==0],x[colnames(geneDataForUroNorm)==1])$p.value }) uroGenesP = paste(formatC(uroGenesPs, digits = 4, format = "f"), collapse = " _ ") ggplot(data = geneDataForUroDF, aes(x = factor(Var1), y = value, colour = factor(Var2)))+geom_boxplot()+geom_hline(yintercept = 0, linetype = "dashed")+ggtitle(uroGenesP)+ theme_bw() + theme(title = element_text(size=9), panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))