GIT_HOME="/research/users/ppolonen/git_home/"

source(file.path(GIT_HOME, "common_scripts/featurematrix/functions_generate_fm.R"))

source(file.path(GIT_HOME, "common_scripts/visualisation/plotting_functions.R"))

library(limma)

library(edgeR)

library(parallel)

library(gridExtra)

library(ComplexHeatmap)

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

# all CGAss:

t.df = read.delim("t.antigen_df.txt", stringsAsFactors=F, header=T)

# CCLE gexp data:

ccle=get(load("CCLE_RNASEQ_SYMBOL_COUNTS.Rdata"))

dge=DGEList(ccle)

# keep.exprs <- filterByExpr(dge) # these remove some of the CGA that are rarely expressed

# dge <- dge[keep.exprs,, keep.lib.sizes=FALSE]

dge <- calcNormFactors(dge, method="TMM")

ccle=voom(dge, plot=T)$E

# CCLE RPKM

ccle2=get(load("CCLE_RNASEQ_SYMBOL_RPKM.Rdata"))

# CCLE meth data:

ccle.meth=get(load("CCLE_combined_methylation.Rdata"))

# CCLE annot:

annot=read.delim("CCLE_sample_info_file_2012-10-18.txt", header=T, stringsAsFactors = F)

find=intersect(intersect(colnames(ccle.meth), intersect(annot$CCLE.name, colnames(ccle))), colnames(ccle2))

ccle.meth=ccle.meth[,match(find, colnames(ccle.meth))]

ccle=ccle[,match(find, colnames(ccle))]

annot=annot[match(find, annot$CCLE.name),]

ccle2=ccle2[,match(find ,colnames(ccle2))]

# compute HLA-scores:

dat_a3=2^ccle[rownames(ccle)%in%c("HLA-DMA",

                                  "HLA-DMB",

                                  "HLA-DPB1",

                                  "HLA-DRA",

                                  "HLA-DRB1"),]+1

gm3=log2(t(apply(dat_a3, 2, gm_mean)))

rownames(gm3)="HLAII_SCORE"

dat2=rbind(ccle, gm3)

annotv=colnames(d)[colnames(d)%in%colnames(ccle.meth)[ccle.meth["CIITA@chr16_10971271",]>0.2]]

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"plasma_cell_myeloma"]="MM"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"mantle_cell_lymphoma"]="MCL"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"diffuse_large_B_cell_lymphoma"]="DLBCL"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"chronic_lymphocytic_leukaemia-small_lymphocytic_lymphoma"]="CLL"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"blast_phase_chronic_myeloid_leukaemia"]="CML"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"anaplastic_large_cell_lymphoma"]="ALCL"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"adult_T_cell_lymphoma-leukaemia"]="TCL"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"acute_myeloid_leukaemia"]="AML"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"acute_lymphoblastic_T_cell_leukaemia"]="T-ALL"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"acute_lymphoblastic_B_cell_leukaemia"]="pre-B-ALL"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"Hodgkin_lymphoma"]="CHL"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%"Burkitt_lymphoma"]="BL"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%c("B_cell_lymphoma_unspecified")]="BCL, unspecified"

annot$Hist.Subtype_hem[annot$Hist.Subtype1%in%c("peripheral_T_cell_lymphoma_unspecified")]="TCL"

d=dat2

groups=factor(annot$Hist.Subtype_hem[grep("lymphoma|leukaemia|myeloma", annot$Hist.Subtype1)])

logicalVectors=lapply(unique(annot$Hist.Subtype_hem), function(g){annot$Hist.Subtype_hem%in%g})

names(logicalVectors)=unique(annot$Hist.Subtype_hem)

logicalVectors=logicalVectors[!is.na(names(logicalVectors))]

logicalVectors=logicalVectors[-6] # remove unspesified

a=names(sort(table(t.df[,3]), decreasing = T))

t.df=do.call(rbind, lapply(a, function(n){

  b=t.df[t.df[,3]%in%n,]

  b=b[order(b[,2], decreasing = T),]

}))

# data:

genelist=unique(t.df[,1])

ccle_antg=ccle[match(genelist,rownames(ccle)),]

antigens_expresssed=colSums(ccle2[match(genelist,rownames(ccle2)),]>0.5, na.rm = T)

disease=annot$Hist.Subtype_hem

go.thr=gsub("@.*", "", rownames(ccle.meth))

ccle.meth=ccle.meth[go.thr%in%genelist,]

go.thr=go.thr[go.thr%in%genelist]

# load("/research/groups/sysgen/PROJECTS/HEMAP_IMMUNOLOGY/data/CCLE_METHYLATION/probelists.Rdata")

highest_anticor=unlist(lapply(unique(go.thr), function(g){

  if(!g%in%go.thr)return(NULL)

  if(!g%in%rownames(ccle))return(NULL)

  m=ccle.meth[go.thr%in%g,grep("lymphoma|leukaemia|myeloma", annot$Hist.Subtype1)]

  cc=as.numeric(ccle[rownames(ccle)%in%g,grep("lymphoma|leukaemia|myeloma", annot$Hist.Subtype1)])

  if(all(rowSums(is.na(m))>140))return(NULL)

  # remove probes that are not found in most datasets:

  m=m[!rowSums(!is.na(m))<20,]

  if(dim(m)[1]==0)return(NULL)

  res=cor(cc,t(m), use="pairwise.complete.obs", method = "spearman")

  # rm poor correlations:

  res=res[res<(-0.2)]

  if(!length(res))return(NULL)

  anticor=rownames(m[which.min(res),])

}))

# data frame

df_anno=data.frame("Disease"=disease, "Number_Antigens_expressed"=as.integer(antigens_expresssed))

# sort based on antigens expressed and disease

sample_order=order(match(annot$Hist.Subtype_hem, c("T-ALL", "pre-B-ALL", "AML", "CML", "CLL","MM","BL","CHL","DLBCL","MCL","BCL, unspecified","TCL","ALCL")), antigens_expresssed)

df_anno=df_anno[sample_order,]

ccle_antg=ccle_antg[,sample_order[!is.na(df_anno$Disease)]]

df_anno=df_anno[!is.na(df_anno$Disease),]

ccle_antg_meth=ccle.meth[,match(colnames(ccle_antg), colnames(ccle.meth))]

cname=colnames(ccle_antg)

colnames(ccle_antg)=gsub("_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE","", colnames(ccle_antg))

colnames(ccle_antg_meth)=gsub("_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE","", colnames(ccle_antg_meth))

ccle_antg_meth=ccle_antg_meth[rownames(ccle_antg_meth)%in%highest_anticor,]

ccle_antg_meth=ccle_antg_meth[match(genelist, gsub("@.*", "", rownames(ccle_antg_meth))),]

ccle_antg_meth2=ccle_antg_meth

ccle_antg_meth2[ccle_antg_meth2<=0.2]=0 # this works best, tested

ccle_antg_meth2[ccle_antg_meth2>0.2]=1

ccle_antg_meth2=ccle_antg_meth2[!rowSums(is.na(ccle_antg_meth2))>100,]

sum_hypometh=colSums(ccle_antg_meth2==0, na.rm = T)

cor(df_anno$Number_Antigens_expressed, sum_hypometh, use="complete.obs")

df_anno=data.frame(df_anno, "Number_Antigens_methylated"=sum_hypometh)

ha = HeatmapAnnotation(df=df_anno[,1,drop=F], Number.Antigens = anno_barplot(df_anno$Number_Antigens_expressed, axis = T, gp = gpar(fill = "indianred")), Number.Me.Antigens = anno_barplot(sum_hypometh, axis = T, gp = gpar(fill = "darkgoldenrod")), height = unit(4, "inch"))

ccle_antg=t(scale(t(ccle_antg)))

ccle_antg[ccle_antg<(-4)]=-4

ccle_antg[ccle_antg>(4)]=4

Heatmap(ccle_antg_meth, top_annotation = ha, name = "TestisAntigens", cluster_columns = F, cluster_rows = F)

pdf("Fig6D_CCLE_ComplexHeatmap_testis_antigens.pdf", height = 15, width = 30)

Heatmap(ccle_antg, top_annotation = ha, name = "TestisAntigens", cluster_columns = F, cluster_rows = F, col=colorRamp2(c(-4,-2, 0,2, 4), c("#063061","#579ec9", "white", "#d45d4a", "#67001e")))

dev.off()