ImmunogenomicLandscape-BC / Git / [8e0848] /Fig7_Univariate_Coxph

Models:
WandaB/
ImmunogenomicLandscape-BC
Downloads: 1
[8e0848]: / Fig7_Univariate_Coxph_survival.R
History
Download this file
914 lines (681 with data), 49.4 kB

GIT_HOME="/research/users/ppolonen/git_home/ImmunogenomicLandscape-BloodCancers/"
source(file.path(GIT_HOME, "common_scripts/visualisation/plotting_functions.R"))
source(file.path(GIT_HOME, "common_scripts/statistics/functions_statistics.R"))
source(file.path(GIT_HOME, "common_scripts/statistics/statistics_wrappers.R"))
source(file.path(GIT_HOME, "common_scripts/featurematrix/compute.pairwise.R"))
source(file.path(GIT_HOME, "common_scripts/featurematrix/functions_generate_fm.R"))

library(RColorBrewer)
library(survival)
library(data.table)
library(ggplot2)
library(ComplexHeatmap)
library(survMisc)
library(survminer)
library(plyr)

setwd("/research/groups/sysgen/PROJECTS/HEMAP_IMMUNOLOGY/petri_work/HEMAP_IMMUNOLOGY/Published_data_figures")

#******************************************** Hemap *********************************************

annot = get(load("Hemap_immunology_Annotations.Rdata"))
annot=annot[!is.na(annot$OS_Time),]

# survival time and status
TIME=annot$OS_Time
STATUS=as.numeric(annot$OS_Status)
TIME2=annot$PFS_Time
STATUS2=as.numeric(annot$PFS_Status)

#************************************** Process gene expression data ************************************
profile=data.matrix(get(load("mixtureM_profile.Rdata")))
profile[profile==-1] = 0
profile2=profile[,colnames(profile)%in%annot$GSM.identifier..sample.]
data=t(get(load("data9544_with_gene_symbols.RData")))
data=data[,colnames(data)%in%annot$GSM.identifier..sample.]

# take only high expressed into account, for CGA score
profile2[data<5]=0

#********************************************************************************************************


#********************************** Make necessary gene lists ***************************************

# other lists from each analysis
cga = read.delim("t.antigen_df.txt", stringsAsFactors=F, header=T)[,1]
co.stim=fread("costim_ligands_final.txt", data.table = F)

hlagenes=c("B2M", "HLA-A", "HLA-B", "HLA-C", "HLA-DMA", "HLA-DMB", "HLA-DPA1", "HLA-DPB1", "HLA-DRA", "HLA-DRB1", "CIITA")
immunoscore=c("HLAI", "HLAII", "CytolyticScore", "Freq.CGA")

# for AML:
MDS=get(load("MDS_genesets.Rdata"))

# significant per disease
microenv=get(load("Hemap_cytolytic_correlated_genes_TableS2_onlysignif.Rdata"))

#****************************************************************************************************
# data for survival
df=data.frame("time"=annot$OS_Time, "status"=annot$OS_Status, "HLAI"=annot$HLAIScore, "HLAII"=annot$HLAIIScore, "CytolyticScore"=annot$CytolyticScore, "Freq.CGA"=colSums(profile2[rownames(profile2)%in%cga,]), stringsAsFactors = F)

#********************************* just AML *********************************
filterv = annot$subclasses%in%"AML"&annot$CELLS_SORTED==0
logicalv=get.logical(annovector = list(annot$GSE.identifier..experiment.), filterv = filterv)
names(logicalv)=unique(paste(names(logicalv), annot$subclasses[filterv]))

filterv = annot$subclasses%in%"AML"&annot$CELLS_SORTED==0
logicalv2=get.logical(annovector = list(annot$subclasses), filterv = filterv)
names(logicalv2)="Hemap_AML"
logicalv=append(logicalv2, logicalv)
logicalv=logicalv[lapply(logicalv, sum)>2]

# HLA, cytolytic
DATA=data.frame("HLAI"=scale(annot$HLAIScore), "HLAII"=scale(annot$HLAIIScore), "CytolyticScore"=scale(annot$CytolyticScore))
aml_res=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
aml_res_multivar=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = F, pretty=F))

DATA_clin=data.frame("Age"=annot$AGE, "Gender"=annot$GENDER)
aml_res_clin=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA_clin,TIME,STATUS, univariate = T, pretty=F))
aml_res_multivar=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA_clin,TIME,STATUS, univariate = F, pretty=F))

DATA_hla=data.frame(t(data[hlagenes,]))
aml_res_HLA=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA_hla,TIME,STATUS, univariate = T, pretty=F))

# costim
DATA=data.frame(scale(t(data[rownames(data)%in%c(co.stim[,1]),])))
aml_res_costim=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# MDS:
DATA=data.frame(scale(t(data[rownames(data)%in%c(MDS$MDS_signature_all_filt),])))
aml_res_MDS=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# microenvironment
ME=microenv[microenv$disease%in%"AML",]
DATA=data.frame(scale(t(data[rownames(data)%in%ME$gene,])))
aml_res_microenvironment=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
aml_res_microenvironment[aml_res_microenvironment$Adj.P<0.1,]

# subtype:
load("Hemap_AML_subtypes.Rdata")
coordinates.subtype=coordinates.subtype[!is.na(coordinates.subtype$subtype),]

samples=lapply(unique(coordinates.subtype$subtype), function(type)coordinates.subtype$ID[coordinates.subtype$subtype%in%type])
DATA=data.frame(do.call(cbind, lapply(samples, function(id)annot$GSM.identifier..sample.%in%id))*1)
colnames(DATA)=gsub("-", ".", unique(coordinates.subtype$subtype))

aml_res_subtype=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# make table S6, adjusted p-value set here to correct for number of comparisons in total:
tableS7=rbind(aml_res, aml_res_subtype, aml_res_costim, aml_res_MDS, aml_res_microenvironment, aml_res_clin)
tableS7=tableS7[tableS7$Name%in%c("Hemap_AML"),]
tableS7$Adj.P=p.adjust(tableS7$P, method="BH")

tableS7[,2]=prettyNum(tableS7[,2])
tableS7[,3]=prettyNum(tableS7[,3])
tableS7[,4]=prettyNum(tableS7[,4])
tableS7[,5]=prettyNum(tableS7[,5])
tableS7[,6]=prettyNum(tableS7[,6])
tableS7[,8]=prettyNum(tableS7[,8])

# annotate these genes, needed later:
tableS7[,1]=gsub("\\.", "-", tableS7[,1])

genelist_signif=data.frame(tableS7[,1], type="Clinical", stringsAsFactors = F)
genelist_signif$type[genelist_signif[,1]%in%c("HLAI", "HLAII", "CytolyticScore", "Freq.CGA")]="ImmunoScores"
genelist_signif$type[genelist_signif[,1]%in%unique(coordinates.subtype$subtype)]="Subtype"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho > 0)",1]]="Stromal/cancer gene (Rho > 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho < 0)",1]]="Stromal/cancer gene (Rho < 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"CTL/NK gene",1]]="CTL/NK gene"
genelist_signif$type[genelist_signif[,1]%in%aml_res_MDS[,1]]="MDS-signature gene"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Inhibitory", co.stim$`Immune checkpoint function`),1]]="Inhibitory ligand"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Stimulatory", co.stim$`Immune checkpoint function`),1]]="Stimulatory ligand"
tableS7$Type=genelist_signif$type
tableS7=tableS7[order(tableS7$Type),]

data.table::fwrite(tableS7[,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Hemap_AML.tsv", sep="\t")
data.table::fwrite(tableS7[tableS7$Adj.P<0.2,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Hemap_AML_signif.tsv", sep="\t")

# save cohort and survival
gexp=data.frame(scale(t(data)))
immunoscore=data.frame("HLAI"=scale(annot$HLAIScore), "HLAII"=scale(annot$HLAIIScore), "CytolyticScore"=scale(annot$CytolyticScore))

samples=lapply(unique(coordinates.subtype$subtype), function(type)coordinates.subtype$ID[coordinates.subtype$subtype%in%type])
subtypes=do.call(cbind, lapply(samples, function(id)annot$GSM.identifier..sample.%in%id))*1 
colnames(subtypes)=gsub("-", ".", unique(coordinates.subtype$subtype))
samp=data.frame(subtypes, "Age"=annot$AGE, "Gender"=annot$GENDER)

save(list = c("gexp","immunoscore", "TIME", "STATUS", "logicalv", "samp"), file="Hemap_AML_survival_data.Rdata")

#************************************* just MM *************************************
filterv = annot$subclasses%in%"Cancer_Myeloma"&!is.na(annot$MM_ISS)
logicalv=get.logical(annovector = list(annot$GSE.identifier..experiment.), filterv = filterv)
names(logicalv)=unique(paste(names(logicalv), annot$subclasses[filterv]))

# just MM
filterv = annot$subclasses%in%"Cancer_Myeloma"&!is.na(annot$MM_ISS)
logicalv2=get.logical(annovector = list(annot$subclasses), filterv = filterv)
names(logicalv2)="MM_all"
logicalv=append(logicalv2, logicalv)
names(logicalv)=c("Hemap_MM", "GSE19784_Hemap_MM", "GSE16716,GSE24080_Hemap_MM")

data.test=data.frame("time"=TIME[logicalv[[1]]], "status"=STATUS[logicalv[[1]]])

ggsurvplot(survfit(Surv(time, status) ~ 1, data = data.test), 
           xlab = "months", 
           ylab = "Overall survival probability")


DATA=data.frame("HLAI"=scale(annot$HLAIScore), "HLAII"=scale(annot$HLAIIScore), "ISS3"=(annot$MM_ISS==3)*1, "ISS1"=(annot$MM_ISS==1)*1, "Freq.CGA"=as.numeric(df$Freq.CGA),"Age"=annot$AGE, "Gender"=annot$GENDER, check.names = F)

mm_res=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
# mm_res_multivar=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA[,3:5],TIME,STATUS, univariate = F, pretty=F))

# costim
DATA=data.frame(scale(t(data[rownames(data)%in%c(co.stim[,1]),])))
mm_res_costim=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

DATA=data.frame(scale(t(data[rownames(data)%in%cga,])))
mm_res_antigen=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# subtype:
load("Hemap_MM_subtypes.Rdata")

samples=lapply(unique(coordinates.subtype$subtype), function(type)coordinates.subtype$ID[coordinates.subtype$subtype%in%type])
DATA=data.frame(do.call(cbind, lapply(samples, function(id)annot$GSM.identifier..sample.%in%id))*1)
colnames(DATA)=gsub("-", "_", unique(coordinates.subtype$subtype))

mm_res_subtype=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# make table S6, adjusted p-value set here to correct for number of comparisons in total:
tableS7=rbind(mm_res, mm_res_costim, mm_res_antigen, mm_res_subtype)
tableS7=tableS7[tableS7$Name%in%c("GSE16716,GSE24080_Hemap_MM"),]
tableS7$Adj.P=p.adjust(tableS7$P, method="BH")


tableS7[,2]=prettyNum(tableS7[,2])
tableS7[,3]=prettyNum(tableS7[,3])
tableS7[,4]=prettyNum(tableS7[,4])
tableS7[,5]=prettyNum(tableS7[,5])
tableS7[,6]=prettyNum(tableS7[,6])
tableS7[,8]=prettyNum(tableS7[,8])

# annotate these genes, needed later:
genelist_signif=data.frame(tableS7[,1], type="Clinical", stringsAsFactors = F)
genelist_signif$type[genelist_signif[,1]%in%c("HLAI", "HLAII", "CytolyticScore", "Freq.CGA")]="ImmunoScores"
genelist_signif$type[genelist_signif[,1]%in%gsub("-", "_", unique(coordinates.subtype$subtype))]="Subtype"
genelist_signif$type[genelist_signif[,1]%in%cga]="CGA"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Inhibitory", co.stim$`Immune checkpoint function`),1]]="Inhibitory ligand"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Stimulatory", co.stim$`Immune checkpoint function`),1]]="Stimulatory ligand"
tableS7$Type=genelist_signif$type
tableS7=tableS7[order(tableS7$Type),]
tableS7[,1]=gsub("\\.", "-", tableS7[,1])

data.table::fwrite(tableS7[,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Hemap_MM_GSE24080.tsv", sep="\t")
data.table::fwrite(tableS7[tableS7$Adj.P<0.2,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Hemap_MM_GSE24080_signif.tsv", sep="\t")


# make table S6, adjusted p-value set here to correct for number of comparisons in total:
tableS7=rbind(mm_res, mm_res_costim, mm_res_antigen)
tableS7=tableS7[tableS7$Name%in%c("GSE19784_Hemap_MM"),]
tableS7$Adj.P=p.adjust(tableS7$P, method="BH")

tableS7[,2]=prettyNum(tableS7[,2])
tableS7[,3]=prettyNum(tableS7[,3])
tableS7[,4]=prettyNum(tableS7[,4])
tableS7[,5]=prettyNum(tableS7[,5])
tableS7[,6]=prettyNum(tableS7[,6])
tableS7[,8]=prettyNum(tableS7[,8])

# annotate these genes, needed later:
genelist_signif=data.frame(tableS7[,1], type="Clinical", stringsAsFactors = F)
genelist_signif$type[genelist_signif[,1]%in%c("HLAI", "HLAII", "CytolyticScore", "Freq.CGA")]="ImmunoScores"
genelist_signif$type[genelist_signif[,1]%in%gsub("-", "_", unique(coordinates.subtype$subtype))]="Subtype"
genelist_signif$type[genelist_signif[,1]%in%cga]="CGA"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Inhibitory", co.stim$`Immune checkpoint function`),1]]="Inhibitory ligand"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Stimulatory", co.stim$`Immune checkpoint function`),1]]="Stimulatory ligand"
tableS7$Type=genelist_signif$type
tableS7=tableS7[order(tableS7$Type),]
tableS7[,1]=gsub("\\.", "-", tableS7[,1])

data.table::fwrite(tableS7[,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Hemap_MM_GSE19784.tsv", sep="\t")
data.table::fwrite(tableS7[tableS7$Adj.P<0.2,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Hemap_MM_GSE19784_signif.tsv", sep="\t")

# save cohort and survival
gexp=data.frame(scale(t(data)))
immunoscore=data.frame("HLAI"=scale(annot$HLAIScore), "HLAII"=scale(annot$HLAIIScore), "Freq.CGA"=as.numeric(df$Freq.CGA), check.names = F)

samples=lapply(unique(coordinates.subtype$subtype), function(type)coordinates.subtype$ID[coordinates.subtype$subtype%in%type])
subtype=do.call(cbind, lapply(samples, function(id)annot$GSM.identifier..sample.%in%id))*1
colnames(subtype)=c(gsub("-", ".", unique(coordinates.subtype$subtype)), colnames(subtype))
samp=data.frame(subtype, "Age"=annot$AGE, "Gender"=annot$GENDER, "ISS3"=(annot$MM_ISS==3)*1, "ISS1"=(annot$MM_ISS==1)*1)

save(list = c("gexp","immunoscore", "TIME", "STATUS", "logicalv", "samp"), file="Hemap_MM_survival_data.Rdata")

#***************************** just DLBCL ************************************
# remove genes not in chapyu dataset:
gexp=get(load("/research/groups/sysgen/PROJECTS/HEMAP_IMMUNOLOGY/data/GSE98588/GSE98588_symbol_rma_normalized.Rdata"))

ME=microenv[microenv$disease%in%"DLBCL"&microenv$category=="Stromal/cancer gene (Rho > 0)"&microenv$Rho>0.4,]
microenv_dlbcl=ME[,1]
microenv_dlbcl=microenv_dlbcl[microenv_dlbcl%in%rownames(gexp)]

filterv = annot$subclasses%in%"BCL_DLBCL"&!is.na(annot$dlbcl_ipi)
logicalv=get.logical(annovector = list(annot$GSE.identifier..experiment.), filterv = filterv)
names(logicalv)=unique(paste(names(logicalv), annot$subclasses[filterv]))

filterv = annot$subclasses%in%"BCL_DLBCL"
logicalv2=get.logical(annovector = list(annot$subclasses), filterv = filterv)
names(logicalv2)="BCL_DLBCL_all"

logicalv=append(logicalv2, logicalv)
logicalv=logicalv[lapply(logicalv, sum)>2]

RCHOP=list(annot$Chemotherapy_RCHOP==1&!is.na(STATUS)&!TIME==0)
CHOP=list(annot$Chemotherapy_CHOP==1&!is.na(STATUS)&!TIME==0)
names(RCHOP)="Hemap_DLBCL_RCHOP"
names(CHOP)="Hemap_DLBCL_CHOP"

logicalv <- append(logicalv, append(RCHOP, CHOP))

# HLA, cytolytic
DATA=data.frame("HLAI"=scale(annot$HLAIScore), "HLAII"=scale(annot$HLAIIScore), "CytolyticScore"=scale(annot$CytolyticScore), "IPI_0to1"=(annot$dlbcl_ipi%in%c(0,1))*1, "IPI_4to5"=(annot$dlbcl_ipi%in%c(4,5))*1, "Freq.CGA"=as.numeric(df$Freq.CGA),  "ABC"=(grepl("ABC", annot$tbLY))*1, "GCB"=(grepl("GCB", annot$tbLY))*1,check.names = F)

dlbcl_res=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
# dlbcl_res_multivar=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA[,4:6],TIME,STATUS, univariate = F, pretty=F))

# costim
DATA=data.frame(scale(t(data[rownames(data)%in%c(co.stim[,1]),])))
dlbcl_res_costim=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# microenv
DATA=data.frame(scale(t(data[rownames(data)%in%c(microenv_dlbcl),])))
dlbcl_res_microenv=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

DATA=data.frame(scale(t(data[rownames(data)%in%cga,])))
dlbcl_res_antigen=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# make table S6, adjusted p-value set here to correct for number of comparisons in total:
tableS7=rbind(dlbcl_res, dlbcl_res_costim, dlbcl_res_microenv, dlbcl_res_antigen)
tableS7=tableS7[tableS7$Name%in%c("Hemap_DLBCL_RCHOP"),]
tableS7$Adj.P=p.adjust(tableS7$P, method="BH")

tableS7[,2]=prettyNum(tableS7[,2])
tableS7[,3]=prettyNum(tableS7[,3])
tableS7[,4]=prettyNum(tableS7[,4])
tableS7[,5]=prettyNum(tableS7[,5])
tableS7[,6]=prettyNum(tableS7[,6])
tableS7[,8]=prettyNum(tableS7[,8])

# annotate these genes, needed later:
genelist_signif=data.frame(tableS7[,1], type="Clinical", stringsAsFactors = F)
genelist_signif$type[genelist_signif[,1]%in%c("HLAI", "HLAII", "CytolyticScore", "Freq.CGA")]="ImmunoScores"
genelist_signif$type[genelist_signif[,1]%in%cga]="CGA"
genelist_signif$type[tableS7$Feature%in%c("ABC", "GCB")]="Subtype"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho > 0)",1]]="Stromal/cancer gene (Rho > 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho < 0)",1]]="Stromal/cancer gene (Rho < 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"CTL/NK gene",1]]="CTL/NK gene"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Inhibitory", co.stim$`Immune checkpoint function`),1]]="Inhibitory ligand"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Stimulatory", co.stim$`Immune checkpoint function`),1]]="Stimulatory ligand"

tableS7$Type=genelist_signif$type
tableS7=tableS7[order(tableS7$Type),]
tableS7[,1]=gsub("\\.", "-", tableS7[,1])

data.table::fwrite(tableS7[,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Hemap_DLBCL_RCHOP.tsv", sep="\t")
data.table::fwrite(tableS7[tableS7$Adj.P<0.2,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Hemap_DLBCL_RCHOP_signif.tsv", sep="\t")


# make table S6, adjusted p-value set here to correct for number of comparisons in total:
tableS7=rbind(dlbcl_res, dlbcl_res_costim, dlbcl_res_microenv, dlbcl_res_antigen) 
tableS7=tableS7[tableS7$Name%in%c("Hemap_DLBCL_CHOP"),]
tableS7$Adj.P=p.adjust(tableS7$P, method="BH")

tableS7[,2]=prettyNum(tableS7[,2])
tableS7[,3]=prettyNum(tableS7[,3])
tableS7[,4]=prettyNum(tableS7[,4])
tableS7[,5]=prettyNum(tableS7[,5])
tableS7[,6]=prettyNum(tableS7[,6])
tableS7[,8]=prettyNum(tableS7[,8])

# annotate these genes, needed later:
genelist_signif=data.frame(tableS7[,1], type="Clinical", stringsAsFactors = F)
genelist_signif$type[genelist_signif[,1]%in%c("HLAI", "HLAII", "CytolyticScore", "Freq.CGA")]="ImmunoScores"
genelist_signif$type[genelist_signif[,1]%in%cga]="CGA"
genelist_signif$type[tableS7$Feature%in%c("ABC", "GCB")]="Subtype"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho > 0)",1]]="Stromal/cancer gene (Rho > 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho < 0)",1]]="Stromal/cancer gene (Rho < 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"CTL/NK gene",1]]="CTL/NK gene"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Inhibitory", co.stim$`Immune checkpoint function`),1]]="Inhibitory ligand"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Stimulatory", co.stim$`Immune checkpoint function`),1]]="Stimulatory ligand"

tableS7$Type=genelist_signif$type
tableS7=tableS7[order(tableS7$Type),]
tableS7[,1]=gsub("\\.", "-", tableS7[,1])

data.table::fwrite(tableS7[,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Hemap_DLBCL_CHOP.tsv", sep="\t")
data.table::fwrite(tableS7[tableS7$Adj.P<0.2,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Hemap_DLBCL_CHOP_signif.tsv", sep="\t")

# save cohort and survival
gexp=data.frame(scale(t(data)))
immunoscore=data.frame("HLAI"=scale(annot$HLAIScore), "HLAII"=scale(annot$HLAIIScore), "Freq.CGA"=as.numeric(df$Freq.CGA), check.names = F)

samp=data.frame("Age"=annot$AGE, "Gender"=annot$GENDER, "IPI_0to1"=(annot$dlbcl_ipi%in%c(0,1))*1, "IPI_4to5"=(annot$dlbcl_ipi%in%c(4,5))*1,"ABC"=(grepl("ABC", annot$tbLY))*1, "GCB"=(grepl("GCB", annot$tbLY))*1)

save(list = c("gexp","immunoscore","samp", "TIME", "STATUS", "logicalv"), file="Hemap_DLBCL_survival_data.Rdata")


#******************************************** TCGA AML *****************************************************
fm_org=get(load("TCGA_AML_FM_DUFVA.Rdata"))
fm=fm_org[,!is.na(fm_org["N:SAMP:CytolyticScore",])]

risks=c("N:CLIN:Age:::::",
        "C:CLIN:acute_myeloid_leukemia_calgb_cytogenetics_risk_category:::::" ,
        "C:CLIN:FISH_test_component:::::",
        "B:GNAB:NPM1:chr5:170814708:170837888:+:y_n_somatic",
        "B:GNAB:FLT3:chr13:28577411:28674729:-:y_n_somatic",
        "B:GNAB:CEBPA:chr19:33790840:33793430:-:y_n_somatic",
        "B:GNAB:TP53:chr17:7565097:7590863:-:y_n_somatic")

df=t(fm[rownames(fm)%in%risks,])
colnames(df)=do.call(rbind, strsplit(colnames(df), ":"))[,3]

data=data.matrix(fm[grepl("GEXP", rownames(fm)),])
rownames(data)=make.unique(do.call(rbind, strsplit(rownames(data), ":"))[,3])


OS=as.numeric(fm["N:CLIN:OS.months..3.31.12:::::",])
TIME=OS
STATUS=as.numeric(fm["C:CLIN:vital_status_TCGA_paper:::::",]=="DECEASED")
PFS=as.numeric(fm["N:CLIN:EFS.months....4.30.13:::::",])

DATA=data.frame("HLAI"=scale(as.numeric(fm["N:SAMP:HLAIScore",])), "HLAII"=scale(as.numeric(fm["N:SAMP:HLAIIScore",])), "CytolyticScore"=scale(as.numeric(fm["N:SAMP:CytolyticScore",])))

logicalv=list("TCGA_AML"=rep(T, dim(DATA)[1]))

TCGA_AML_res=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
# TCGA_AML_multivar=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = F, pretty=F))

# clinical
DATA_clin=data.frame("Age"=scale(as.numeric(fm["N:CLIN:Age:::::",])), "Blast.percentage"=scale(as.numeric(fm["N:CLIN:X.BM.Blast:::::",])), "Gender"=as.character(fm["C:CLIN:Sex:::::",]))
TCGA_AML_clin=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA_clin,TIME,STATUS, univariate = T, pretty=F))

# costim
DATA=data.frame(scale(t(data[rownames(data)%in%c(co.stim[,1]),])))
TCGA_AML_costim=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

DATA_hla=data.frame(scale(t(data[rownames(data)%in%hlagenes,])))
TCGA_AML_HLA=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA_hla,TIME,STATUS, univariate = T, pretty=F))

# MDS
DATA=data.frame(scale(t(data[rownames(data)%in%c(MDS$MDS_signature_all_filt),])))
TCGA_AML_MDS=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# mikroenvironment
ME=microenv[microenv$disease%in%"AML",]
DATA=data.frame(scale(t(data[rownames(data)%in%ME$gene,])))
TCGA_res_microenvironment=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
TCGA_res_microenvironment[TCGA_res_microenvironment$Adj.P<0.1,]

# subtype:
load("TCGA_AML_subtypes.Rdata")

DATA_subtypes=data.frame(do.call(cbind, get.logical(list(coordinates.subtype$subtype)))*1)

TCGA_aml_res_subtype=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA_subtypes,TIME,STATUS, univariate = T, pretty=F))

# make table S6, adjusted p-value set here to correct for number of comparisons in total:
tableS7=rbind(TCGA_AML_res, TCGA_AML_costim, TCGA_AML_MDS, TCGA_res_microenvironment, TCGA_aml_res_subtype, TCGA_AML_clin)
tableS7$Adj.P=p.adjust(tableS7$P, method="BH")

tableS7[,2]=prettyNum(tableS7[,2])
tableS7[,3]=prettyNum(tableS7[,3])
tableS7[,4]=prettyNum(tableS7[,4])
tableS7[,5]=prettyNum(tableS7[,5])
tableS7[,6]=prettyNum(tableS7[,6])
tableS7[,8]=prettyNum(tableS7[,8])

# annotate these genes, needed later:
tableS7[,1]=gsub("\\.", "-", tableS7[,1])

genelist_signif=data.frame(tableS7[,1], type="Clinical", stringsAsFactors = F)
genelist_signif$type[genelist_signif[,1]%in%c("HLAI", "HLAII", "CytolyticScore", "Freq.CGA")]="ImmunoScores"
genelist_signif$type[genelist_signif[,1]%in%unique(coordinates.subtype$subtype)]="Subtype"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho > 0)",1]]="Stromal/cancer gene (Rho > 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho < 0)",1]]="Stromal/cancer gene (Rho < 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"CTL/NK gene",1]]="CTL/NK gene"
genelist_signif$type[genelist_signif[,1]%in%aml_res_MDS[,1]]="MDS-signature gene"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Inhibitory", co.stim$`Immune checkpoint function`),1]]="Inhibitory ligand"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Stimulatory", co.stim$`Immune checkpoint function`),1]]="Stimulatory ligand"

tableS7$Type=genelist_signif$type
tableS7=tableS7[order(tableS7$Type),]

data.table::fwrite(tableS7[,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_TCGA_AML.tsv", sep="\t")
data.table::fwrite(tableS7[tableS7$Adj.P<0.2,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_TCGA_AML_signif.tsv", sep="\t")

# save cohort and survival
gexp=data.frame(scale(t(data)))
immunoscore=data.frame("HLAI"=scale(as.numeric(fm["N:SAMP:HLAIScore",])), "HLAII"=scale(as.numeric(fm["N:SAMP:HLAIIScore",])), "CytolyticScore"=scale(as.numeric(fm["N:SAMP:CytolyticScore",])))

samp=cbind(DATA_clin, DATA_subtypes)

save(list = c("gexp","immunoscore", "samp", "TIME", "STATUS", "logicalv"), file="TCGA_AML_survival_data.Rdata")

#********************************************************* Compass MM *********************************************************
fm=get(load("MM_COMPASS_FM.Rdata"))
annot=get(load("MM_COMPASS_ANNOT.Rdata"))

annot=annot[match(colnames(fm), rownames(annot)),]

data=fm[grepl("N:GEXP:", rownames(fm)),]
rownames(data)=gsub("N:GEXP:", "", rownames(data))

data.mut=fm[grepl("B:GNAB:", rownames(fm)),]
rownames(data.mut)=gsub("B:GNAB:", "", rownames(data.mut))
data.mut.filt=data.mut[!rowSums(data.mut, na.rm = T)<10,]

TIME=as.numeric(fm["N:CLIN:OS",])
STATUS=as.numeric(fm["B:CLIN:STATUS",])

STATUS[TIME>1825&!is.na(STATUS)&STATUS==1]=0 # change to 4year survival, 5 year sharp drop
TIME[TIME>1825]=1825

# compute antigen scores
t.df = read.delim("t.antigen_df.txt", stringsAsFactors=F, header=T)
t.df=t.df[order(t.df[,3]),]
genelist=t.df[,1]

data.test=data.frame("time"=TIME[logicalv[[1]]]*0.0328767, "status"=STATUS[logicalv[[1]]])

ggsurvplot(survfit(Surv(time, status) ~ 1, data = data.test), 
           xlab = "months", 
           ylab = "Overall survival probability")

expressed_testis_num=as.numeric(fm["N:SAMP:nCGA",])

l.regulon.gene=regulon.feats(fm, co.stim[,1])

hlagenes=c("B2M", "HLA-A", "HLA-B", "HLA-C", "HLA-DMA", "HLA-DMB", "HLA-DPA1", "HLA-DPB1", "HLA-DRA", "HLA-DRB1", "CIITA")

DATAcostim=scale(data.frame(t(gexp[rownames(gexp)%in%co.stim[,1],])))
DATAhla=scale(data.frame(t(gexp[rownames(gexp)%in%hlagenes,])))
DATAcga=scale(data.frame(t(gexp[rownames(gexp)%in%t.df[,1],])))

DATA=data.frame("HLAI"=scale(as.numeric(fm["N:SAMP:HLAIScore",])), "HLAII"=scale(as.numeric(fm["N:SAMP:HLAIIScore",])), "ISS1"=as.numeric(fm["B:CLIN:R_ISS_1",]),"ISS2"=as.numeric(fm["B:CLIN:R_ISS_2",]),"ISS3"=as.numeric(fm["B:CLIN:R_ISS_3",]), "Freq.CGA"=as.numeric(expressed_testis_num), stringsAsFactors = F)

logicalv=list(!is.na(expressed_testis_num))

names(logicalv)="CoMMPass"

# plot overall survival
r=lapply(seq(logicalv), function(i){
  data.test=data.frame("time"=TIME[logicalv[[i]]], "status"=STATUS[logicalv[[i]]])
  
  ggsurvplot(survfit(Surv(time, status) ~ 1, data = data.test), 
             xlab = "months", 
             ylab = "Overall survival probability",title=(names(logicalv)[i])
  )
})

names(r)=names(logicalv)

pdf("CoMMpass_MM_cohorts.pdf", width =5, height = ceiling(length(r)/2)*2.75)
plots.together=arrange_ggsurvplots(r, print = TRUE, ncol = 2, nrow = ceiling(length(r)/2))
dev.off()

commpass_res=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
commpass_res_multivar=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = F, pretty=F))

# costim, no mut
DATA=data.frame(scale(t(data[rownames(data)%in%c(co.stim[,1]),])))
commpass_costim=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# mutations
mut=t(data.mut.filt[rownames(data.mut.filt)%in%c(t.df[,1]),])
colnames(mut)=paste0("MUT:", colnames(mut))
DATA=data.frame(scale(t(data[rownames(data)%in%c(t.df[,1]),])), mut)

commpass_antig=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# subtype:
load("CoMMpass_MM_subtypes.Rdata")
coordinates.subtype=coordinates.subtype[match(colnames(fm), coordinates.subtype$ID),]

coordinates.subtype$subtype[coordinates.subtype$cluster=="CGA_Prolif"&!is.na(coordinates.subtype$cluster)]="CGA_Prolif"
DATA_subtypes=data.frame(do.call(cbind, get.logical(list(coordinates.subtype$subtype)))*1)
commpass_subtype=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA_subtypes,TIME,STATUS, univariate = T, pretty=F))

# make table S6, adjusted p-value set here to correct for number of comparisons in total:
tableS7=rbind(commpass_res, commpass_costim, commpass_antig, commpass_subtype)
tableS7$Adj.P=p.adjust(tableS7$P, method="BH")

# annotate these genes, needed later:
tableS7[,2]=prettyNum(tableS7[,2])
tableS7[,3]=prettyNum(tableS7[,3])
tableS7[,4]=prettyNum(tableS7[,4])
tableS7[,5]=prettyNum(tableS7[,5])
tableS7[,6]=prettyNum(tableS7[,6])
tableS7[,8]=prettyNum(tableS7[,8])

# annotate these genes, needed later:
genelist_signif=data.frame(tableS7[,1], type="Clinical", stringsAsFactors = F)
genelist_signif$type[genelist_signif[,1]%in%c("HLAI", "HLAII", "CytolyticScore", "Freq.CGA")]="ImmunoScores"
genelist_signif$type[genelist_signif[,1]%in%cga]="CGA"
genelist_signif$type[genelist_signif[,1]%in%commpass_subtype$Feature]="Subtype"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Inhibitory", co.stim$`Immune checkpoint function`),1]]="Inhibitory ligand"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Stimulatory", co.stim$`Immune checkpoint function`),1]]="Stimulatory ligand"

tableS7$Type=genelist_signif$type
tableS7=tableS7[order(tableS7$Type),]
tableS7[,1]=gsub("\\.", "-", tableS7[,1])

data.table::fwrite(tableS7[,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_CoMMpass.tsv", sep="\t")
data.table::fwrite(tableS7[tableS7$Adj.P<0.2,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_CoMMpass_signif.tsv", sep="\t")

# save cohort and survival
gexp=data.frame(scale(t(data)))
immunoscore=data.frame("HLAI"=scale(as.numeric(fm["N:SAMP:HLAIScore",])), "HLAII"=scale(as.numeric(fm["N:SAMP:HLAIIScore",])), "Freq.CGA"=as.numeric(expressed_testis_num), stringsAsFactors = F)

samp=cbind(data.frame("ISS1"=as.numeric(fm["B:CLIN:R_ISS_1",]),"ISS2"=as.numeric(fm["B:CLIN:R_ISS_2",]),"ISS3"=as.numeric(fm["B:CLIN:R_ISS_3",]), stringsAsFactors = F), DATA_subtypes)

save(list = c("gexp","immunoscore", "samp", "TIME", "STATUS", "logicalv"), file="CoMMpass_survival_data.Rdata")

#********************************************************* Chapyu DLBCL *********************************************************

fm=get(load("GSE98588_fm.Rdata"))
annot=get(load("GSE98588_annot.Rdata"))

TIME=as.numeric(fm["N:CLIN:OS",])
PFS=as.numeric(fm["N:CLIN:PFS",])
STATUS=as.numeric(fm["B:CLIN:OS_STAT",])
STATUS2=as.numeric(fm["B:CLIN:PFS_STAT",])

# compute antigen scores
genelist=t.df[,1]

data=get(load("GSE98588_symbol_rma_normalized.Rdata"))
colnames(data)=gsub("_DLBCL", "", colnames(data))
gexp=data

load("GSE98588_DLBCL_mixtureM_profile.Rdata")
profile[profile==-1] = 0
profile[data.matrix(data)<5]=0

expressed_testis_num=as.numeric(colSums(profile[rownames(profile)%in%genelist,]))

cnv_annot=fread("41591_2018_16_MOESM8_ESM_CNV_ANNOT.txt", data.table=F)

l.regulon.gene=regulon.feats(fm, c(co.stim[,1], hlagenes), cnv_annot)

ME=microenv[microenv$disease%in%"DLBCL"&microenv$category=="Stromal/cancer gene (Rho > 0)"&microenv$Rho>0.4,]

# all costim-cytolytic-CGA correlated mutations:
costim_feats=data.matrix(fm[rownames(fm)%in%unlist(l.regulon.gene),])
rownames(costim_feats)=sapply(rownames(costim_feats), function(n){
  if(!grepl("CNVR", n))return(n)
  a=names(l.regulon.gene)[sapply(l.regulon.gene, function(a)any(a%in%n))]
  if(length(a))paste0(n,"@", paste(a, collapse=","))
})

DATAmut=data.frame(t(costim_feats[grepl("GNAB|CNVR", rownames(costim_feats)),]))

DATA=data.frame("HLAI"=scale(as.numeric(fm["N:SAMP:HLAIScore",])), "HLAII"=scale(as.numeric(fm["N:SAMP:HLAIIScore",])), "CytolyticScore"=scale(as.numeric(fm["N:SAMP:CytolyticScore",])), "Freq.CGA"=as.numeric(expressed_testis_num), stringsAsFactors = F)

logicalv=list(rep(T, dim(DATA)[1]))
names(logicalv)="GSE98588_DLBCL"

GSE98588_DLBCL_res=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
# GSE98588_DLBCL_multivar=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA[,4:6],TIME,STATUS, univariate = F, pretty=F))

# Clinical
DATA_clin=data.frame("IPI_0to1"=annot$IPI%in%c(0,1)*1,"IPI_4to5"=annot$IPI%in%c(4,5)*1, "ABC"=(annot$COO_byGEP=="ABC")*1 ,"GCB"=(annot$COO_byGEP=="GCB")*1,  stringsAsFactors = F)
GSE98588_DLBCL_clin=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA_clin,TIME,STATUS, univariate = T, pretty=F))

# immuno-editing:
GSE98588_DLBCL_immunoediting=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATAmut,TIME,STATUS, univariate = T, pretty=F))
GSE98588_DLBCL_immunoediting$Feature=gsub("B.GNAB.", "MUT:", GSE98588_DLBCL_immunoediting$Feature)
GSE98588_DLBCL_immunoediting$Feature=gsub("HLA.", "HLA-", GSE98588_DLBCL_immunoediting$Feature)
GSE98588_DLBCL_immunoediting$Feature=gsub("B.CNVR.", "", GSE98588_DLBCL_immunoediting$Feature)

# costim
DATA=data.frame(scale(t(data[rownames(data)%in%c(co.stim[,1]),])), DATAmut, "costim.mut"=(rowSums(DATAmut)))
GSE98588_DLBCL_costim=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# monocyte
DATA=data.frame(scale(t(data[rownames(data)%in%c(ME[,1]),])))
GSE98588_DLBCL_microenvironment=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

DATA=data.frame(scale(t(data[rownames(data)%in%c(t.df[,1]),])))
GSE98588_DLBCL_antigen=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# make table S6, adjusted p-value set here to correct for number of comparisons in total:
tableS7=rbind(GSE98588_DLBCL_res,GSE98588_DLBCL_clin, GSE98588_DLBCL_costim, GSE98588_DLBCL_microenvironment, GSE98588_DLBCL_antigen, GSE98588_DLBCL_immunoediting) 
tableS7$Adj.P=p.adjust(tableS7$P, method="BH")


tableS7[,2]=prettyNum(tableS7[,2])
tableS7[,3]=prettyNum(tableS7[,3])
tableS7[,4]=prettyNum(tableS7[,4])
tableS7[,5]=prettyNum(tableS7[,5])
tableS7[,6]=prettyNum(tableS7[,6])
tableS7[,8]=prettyNum(tableS7[,8])

# annotate these genes, needed later:
genelist_signif=data.frame(tableS7[,1], type="Clinical", stringsAsFactors = F)
genelist_signif$type[genelist_signif[,1]%in%c("HLAI", "HLAII", "CytolyticScore", "Freq.CGA")]="ImmunoScores"
genelist_signif$type[genelist_signif[,1]%in%colnames(DATAmut)]="Immune Editing mutation"
genelist_signif$type[genelist_signif[,1]%in%cga]="CGA"
genelist_signif$type[tableS7$Feature%in%c("ABC", "GCB")]="Subtype"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho > 0)",1]]="Stromal/cancer gene (Rho > 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho < 0)",1]]="Stromal/cancer gene (Rho < 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"CTL/NK gene",1]]="CTL/NK gene"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Inhibitory", co.stim$`Immune checkpoint function`),1]]="Inhibitory ligand"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Stimulatory", co.stim$`Immune checkpoint function`),1]]="Stimulatory ligand"

tableS7$Type=genelist_signif$type
tableS7=tableS7[order(tableS7$Type),]
tableS7[,1]=gsub("\\.", "-", tableS7[,1])

data.table::fwrite(tableS7[,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_GSE98588_DLBCL.tsv", sep="\t")
data.table::fwrite(tableS7[tableS7$Adj.P<0.2,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_GSE98588_DLBCL_signif.tsv", sep="\t")

# save cohort and survival
gexp=data.frame(scale(t(data)))
immunoscore=data.frame("HLAI"=scale(as.numeric(fm["N:SAMP:HLAIScore",])), "HLAII"=scale(as.numeric(fm["N:SAMP:HLAIIScore",])), "CytolyticScore"=scale(as.numeric(fm["N:SAMP:CytolyticScore",])), "Freq.CGA"=as.numeric(expressed_testis_num), stringsAsFactors = F)

samp=DATA_clin

save(list = c("gexp","immunoscore", "samp", "TIME", "STATUS", "logicalv"), file="GSE98588_DLBCL_survival_data.Rdata")

#********************************************************* Reddy DLBCL *********************************************************
fm=get(load("REDDY_DLBCL_fm.Rdata"))
annot=get(load("REDDY_DLBCL_annot.Rdata"))

data=data.matrix(fm[grepl("GEXP", rownames(fm)),])
rownames(data)=gsub("N:GEXP:", "", rownames(data))

TIME=as.numeric(annot$Overall.Survival.years)
STATUS=(as.numeric(annot$Censored)==0)*1 # 0 indicates no censoring, meaning that the death was observed, whereas a 1 indicates that the patient was alive

ME=microenv[microenv$disease%in%"DLBCL"&microenv$category=="Stromal/cancer gene (Rho > 0)"&microenv$Rho>0.4,]

l.regulon.gene=regulon.feats(fm, c(co.stim[,1], hlagenes))

# all costim-cytolytic-CGA correlated mutations:
costim_feats=data.matrix(fm[rownames(fm)%in%unlist(l.regulon.gene),])
rownames(costim_feats)=sapply(rownames(costim_feats), function(n){
  if(!grepl("CNVR", n))return(n)
  a=names(l.regulon.gene)[sapply(l.regulon.gene, function(a)any(a%in%n))]
  # if(length(a))paste0(n,"@", paste(a, collapse=","))
})

DATAmut=data.frame(t(costim_feats[grepl("GNAB|CNVR", rownames(costim_feats)),]))

DATA=data.frame("HLAII"=scale(as.numeric(fm["N:SAMP:HLAIIScore",])), "CytolyticScore"=scale(as.numeric(fm["N:SAMP:CytolyticScore",])), stringsAsFactors = F)

logicalv=list("Reddy_DLBCL"=!is.na(data[1,]), "Reddy_DLBCL_ABC"=!is.na(data[1,])&annot$ABC.GCB..RNAseq.=="ABC", "Reddy_DLBCL_GCB"=!is.na(data[1,])&annot$ABC.GCB..RNAseq.=="GCB")

Reddy_DLBCL_res=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
Reddy_DLBCL_multivar=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = F, pretty=F))

# Clinical
DATA_clin=data.frame("IPI_0to1"=annot$IPI%in%c(0,1)*1,"IPI_4to5"=annot$IPI%in%c(4,5)*1, "ABC"=(annot$ABC.GCB..RNAseq.=="ABC")*1 ,"GCB"=(annot$ABC.GCB..RNAseq.=="GCB")*1,  stringsAsFactors = F)
Reddy_DLBCL_clin=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA_clin,TIME,STATUS, univariate = T, pretty=F))

# immuno-editing:
Reddy_DLBCL_immunoediting=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATAmut,TIME,STATUS, univariate = T, pretty=F))
Reddy_DLBCL_immunoediting$Feature=gsub("B.GNAB.", "MUT:", Reddy_DLBCL_immunoediting$Feature)
Reddy_DLBCL_immunoediting$Feature=gsub("HLA.", "HLA-", Reddy_DLBCL_immunoediting$Feature)
Reddy_DLBCL_immunoediting$Feature=gsub("N.CNVR.", "", Reddy_DLBCL_immunoediting$Feature)

# costim
DATA=data.frame(scale(t(data[rownames(data)%in%c(co.stim[,1]),])))
Reddy_DLBCL_costim=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# microenvironment
DATA=data.frame(scale(t(data[rownames(data)%in%c(ME[,1]),])))
Reddy_DLBCL_microenvironment=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# make table S6, adjusted p-value set here to correct for number of comparisons in total:
tableS7=rbind(Reddy_DLBCL_res,Reddy_DLBCL_clin[!is.na(Reddy_DLBCL_clin$`exp(coef)`),], Reddy_DLBCL_immunoediting, Reddy_DLBCL_costim, Reddy_DLBCL_microenvironment) 
tableS7=tableS7[tableS7$Name%in%"Reddy_DLBCL",]
tableS7$Adj.P=p.adjust(tableS7$P, method="BH")


tableS7[,2]=prettyNum(tableS7[,2])
tableS7[,3]=prettyNum(tableS7[,3])
tableS7[,4]=prettyNum(tableS7[,4])
tableS7[,5]=prettyNum(tableS7[,5])
tableS7[,6]=prettyNum(tableS7[,6])
tableS7[,8]=prettyNum(tableS7[,8])

# annotate these genes, needed later:
genelist_signif=data.frame(tableS7[,1], type="Clinical", stringsAsFactors = F)
genelist_signif$type[genelist_signif[,1]%in%c("HLAI", "HLAII", "CytolyticScore", "Freq.CGA")]="ImmunoScores"
genelist_signif$type[genelist_signif[,1]%in%Reddy_DLBCL_immunoediting$Feature]="Immune Editing mutation"
genelist_signif$type[genelist_signif[,1]%in%cga]="CGA"
genelist_signif$type[tableS7$Feature%in%c("ABC", "GCB")]="Subtype"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho > 0)",1]]="Stromal/cancer gene (Rho > 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho < 0)",1]]="Stromal/cancer gene (Rho < 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"CTL/NK gene",1]]="CTL/NK gene"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Inhibitory", co.stim$`Immune checkpoint function`),1]]="Inhibitory ligand"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Stimulatory", co.stim$`Immune checkpoint function`),1]]="Stimulatory ligand"

tableS7$Type=genelist_signif$type
tableS7=tableS7[order(tableS7$Type),]
tableS7[,1]=gsub("\\.", "-", tableS7[,1])

data.table::fwrite(tableS7[,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Reddy_DLBCL.tsv", sep="\t")
data.table::fwrite(tableS7[tableS7$Adj.P<0.2,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_Reddy_DLBCL_signif.tsv", sep="\t")

# save cohort and survival
gexp=data.frame(scale(t(data)))
immunoscore=data.frame("HLAII"=scale(as.numeric(fm["N:SAMP:HLAIIScore",])), "CytolyticScore"=scale(as.numeric(fm["N:SAMP:CytolyticScore",])), stringsAsFactors = F)

samp=DATA_clin

save(list = c("gexp","immunoscore", "samp", "TIME", "STATUS", "logicalv"), file="Reddy_DLBCL_survival_data.Rdata")


#************************************ beatAML:
load("BeatAML_fm.Rdata")
annot=get(load("BeatAML_fm_annot.Rdata"))

gexp=fm[grepl("GEXP", rownames(fm)),]
rownames(gexp)=gsub("N:GEXP:", "", rownames(gexp))

annot$vitalStatus[annot$vitalStatus=="Unknown"]=NA
TIME=as.numeric(annot$overallSurvival)
STATUS=as.numeric(annot$vitalStatus=="Dead")

STATUS[TIME>1825&!is.na(STATUS)&STATUS==1]=0 # change to 5year survival
TIME[TIME>1825]=1825
TIME[TIME==0]=0.1

filtv=annot$specimenType%in%"Bone Marrow Aspirate"&!(is.na(TIME)|is.na(STATUS)|is.na(annot$TCGA_coord))&annot$causeOfDeath%in%c("Alive", "Dead-Disease", "Dead-Unknown")
filtv=annot$specimenType%in%"Bone Marrow Aspirate"&!(is.na(annot$TCGA_coord))&annot$causeOfDeath%in%c("Alive", "Dead-Disease", "Dead-Unknown")
filtv2=annot$specimenType%in%"Bone Marrow Aspirate"&!(is.na(annot$TCGA_coord))&annot$TCGA_coord%in%c("CMP-like","MDS-like","Monocyte-like")&annot$causeOfDeath%in%c("Alive", "Dead-Disease", "Dead-Unknown")
filtv3=annot$specimenType%in%"Bone Marrow Aspirate"&!(is.na(annot$TCGA_coord))&annot$TCGA_coord%in%c("MDS-like")&annot$causeOfDeath%in%c("Alive", "Dead-Disease", "Dead-Unknown")
filtv4=annot$specimenType%in%"Bone Marrow Aspirate"&!(is.na(annot$TCGA_coord))&annot$TCGA_coord%in%c("CMP-like")&annot$causeOfDeath%in%c("Alive", "Dead-Disease", "Dead-Unknown")
filtv5=annot$specimenType%in%"Bone Marrow Aspirate"&!(is.na(annot$TCGA_coord))&annot$TCGA_coord%in%c("Monocyte-like")&annot$causeOfDeath%in%c("Alive", "Dead-Disease", "Dead-Unknown")

logicalv=list("BeatAML"=!(is.na(TIME)|is.na(STATUS)|is.na(annot$TCGA_coord))&annot$causeOfDeath%in%c("Alive", "Dead-Disease", "Dead-Unknown"), "beatAML_BMonly"=filtv, "normal_karyotype"=filtv2, "MDS-like"=filtv3,  "CMP-like"=filtv4,  "Monocyte-like"=filtv5)

# plot overall survival
r=lapply(seq(logicalv), function(i){
  ggsurvplot(survfit(Surv(time, status) ~ 1, data = data.frame("time"=TIME[logicalv[[i]]], "status"=STATUS[logicalv[[i]]])), 
             xlab = "months", 
             ylab = "Overall survival probability",title=(names(logicalv)[i])
  )
})

names(r)=names(logicalv)

pdf("BeatAML_cohorts.pdf", width =5, height = ceiling(length(r)/2)*2.75)
plots.together=arrange_ggsurvplots(r, print = TRUE, ncol = 2, nrow = ceiling(length(r)/2))
dev.off()

data.test=data.frame("time"=TIME[logicalv[[1]]]*0.0328767, "status"=STATUS[logicalv[[1]]])

ggsurvplot(survfit(Surv(time, status) ~ 1, data = data.test), 
           xlab = "months", 
           ylab = "Overall survival probability")

# HLA, cytolytic
DATA=data.frame(scale(t(fm[grepl("N:SAMP:.*.Score", rownames(fm)),])), check.names = F)
colnames(DATA)=c("HLAI", "HLAII", "CytolyticScore")
beatAML=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))

# Clinical:
DATA_clin=data.frame(t(fm[grepl("ELN2017", rownames(fm)),]),"Age"=as.numeric(fm[grepl("ageAtDiagnosis", rownames(fm)),,drop=F]), t(fm[grepl("BM_Transplant", rownames(fm)),,drop=F]),t(fm[grepl("B:CLIN:priorMDS_TRUE", rownames(fm)),,drop=F]), t(fm[grepl("in_PB|in_BM|B:CLIN:is|finalFusion_Complex", rownames(fm)),,drop=F]))
DATA_clin$Age[is.na(DATA_clin$Age)]=median(DATA_clin$Age, na.rm = T) #one observation --> set to median

colnames(DATA_clin)=gsub("..CLIN.", "", colnames(DATA_clin))
beatAML_clin=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA_clin,TIME,STATUS, univariate = T, pretty=F))

# stroma
ME=microenv[microenv$disease=="AML",]
DATA=data.frame(t(gexp[rownames(gexp)%in%ME$gene,]))
beatAML_microenv=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
beatAML_microenv[beatAML_microenv$P<0.01,]

# MDS
DATA=data.frame(t(gexp[rownames(gexp)%in%MDS$MDS_signature_all_filt,]))
beatAML_MDS=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
beatAML_MDS[beatAML_MDS$P<0.05,]

# costim
DATA=data.frame(t(gexp[rownames(gexp)%in%co.stim[,1],]))
beatAML_costim=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA,TIME,STATUS, univariate = T, pretty=F))
beatAML_costim[beatAML_costim$P<0.05,]

# subtype:
load("BeatAML_subtypes.Rdata")
coordinates.subtype$subtype[coordinates.subtype$subtype%in%"Progenitor-like"]="CMP-like"
coordinates.subtype=coordinates.subtype[match(colnames(fm), coordinates.subtype$ID),]

DATA_subtypes=data.frame(do.call(cbind, get.logical(list(coordinates.subtype$subtype)))*1)

beataml_res_subtype=do.call(rbind, lapply(seq(logicalv), fun.get.cox, logicalv, DATA_subtypes,TIME,STATUS, univariate = T, pretty=F))

fun.kapplanMeier(TIME[logicalv[[1]]], STATUS[logicalv[[1]]],GROUPS=coordinates.subtype$subtype[logicalv[[1]]], MONTHS=T, PVAL=1, INDIVIDUAL_GROUPS=F,LWD = 1, NAME = "Prognostic Index - validation")

# make table S6, adjusted p-value set here to correct for number of comparisons in total:
tableS7=rbind(beatAML, beataml_res_subtype, beatAML_costim, beatAML_MDS, beatAML_microenv)
tableS7=tableS7[tableS7$Name%in%"BeatAML",]
tableS7$Adj.P=p.adjust(tableS7$P, method="BH")

tableS7[,2]=prettyNum(tableS7[,2])
tableS7[,3]=prettyNum(tableS7[,3])
tableS7[,4]=prettyNum(tableS7[,4])
tableS7[,5]=prettyNum(tableS7[,5])
tableS7[,6]=prettyNum(tableS7[,6])
tableS7[,8]=prettyNum(tableS7[,8])

# annotate these genes, needed later:
genelist_signif=data.frame(tableS7[,1], type="", stringsAsFactors = F)
genelist_signif$type[genelist_signif[,1]%in%c("HLAI", "HLAII", "CytolyticScore", "Freq.CGA")]="ImmunoScores"
genelist_signif$type[genelist_signif[,1]%in%beataml_res_subtype$Feature]="Subtype"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho > 0)",1]]="Stromal/cancer gene (Rho > 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"Stromal/cancer gene (Rho < 0)",1]]="Stromal/cancer gene (Rho < 0)"
genelist_signif$type[genelist_signif[,1]%in%ME[ME$category%in%"CTL/NK gene",1]]="CTL/NK gene"
genelist_signif$type[genelist_signif[,1]%in%MDS$MDS_signature_all_filt]="MDS-signature gene"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Inhibitory", co.stim$`Immune checkpoint function`),1]]="Inhibitory ligand"
genelist_signif$type[genelist_signif[,1]%in%co.stim[grepl("Stimulatory", co.stim$`Immune checkpoint function`),1]]="Stimulatory ligand"

tableS7$Type=genelist_signif$type
tableS7=tableS7[order(tableS7$Type),]
tableS7[,1]=gsub("\\.", "-", tableS7[,1])

data.table::fwrite(tableS7[,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_beatAML.tsv", sep="\t")
data.table::fwrite(tableS7[tableS7$Adj.P<0.2,c(1,11,10,2,3,4,5,6,8,9)], "tableS7_beatAML_signif.tsv", sep="\t")

# save cohort and survival
gexp=data.frame(scale(t(gexp)))
immunoscore=data.frame(scale(t(fm[grepl("N:SAMP:.*.Score", rownames(fm)),])), check.names = F)
colnames(immunoscore)=c("HLAI", "HLAII", "CytolyticScore")

samp=cbind(DATA_clin, DATA_subtypes)

save(list = c("gexp","immunoscore", "samp", "TIME", "STATUS", "logicalv"), file="BeatAML_survival_data.Rdata")

#******************************** make excel table
files=list.files(pattern = "tableS7")
files=files[!grepl("_CHOP", files)]

univariate.results=lapply(files[grepl("signif", files)], fread, data.table=F)
names(univariate.results)=gsub("tableS7_|_signif.tsv", "", files[grepl("signif", files)])
univariate.results=univariate.results[c(grep("DLBCL", names(univariate.results)), grep("MM", names(univariate.results)), grep("AML", names(univariate.results)))]

require(openxlsx)

write.xlsx(univariate.results, file = "TableS7_cox.xlsx", )