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

Models:
WandaB/
ImmunogenomicLandscape-BC
Downloads: 1
[8e0848]: / Fig5_DE_analysis_costim.R
History
Download this file
486 lines (373 with data), 26.1 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/pathway_analysis/functions.GSEA.R"))
source(file.path(GIT_HOME, "common_scripts/scRNA/functions.scRNA.analysis.R"))
source(file.path(GIT_HOME, "common_scripts/statistics/useful_functions.R"))

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

# analyze FIMM and Galen AML and find costim associations to certain clusters.
co.stim=data.table::fread("costim_ligands_final.txt", data.table = F)[,c(1,3,5)]
co.stim=rbind(co.stim, c("FGL1", "LAG3", "Inhibitory"))
co.stimR.f=unlist(strsplit(co.stim$`Receptor gene`, ", "))

co.stim.all=data.table::fread("costim_ligands_final.txt", data.table = F)
ligand.name=gsub(" \\(\\)", "", paste0(co.stim.all$Gene, " (", co.stim.all$`Common name`, ")"))
Receptor.name=gsub(" \\(\\)", "", paste0(co.stim.all$`Receptor gene`, " (", co.stim.all$`Receptor common name`, ")"))
Receptor.name[co.stim.all$`Receptor gene`==co.stim.all$`Receptor common name`]=co.stim.all$`Receptor gene`[co.stim.all$`Receptor gene`==co.stim.all$`Receptor common name`]

rec=strsplit(co.stim$`Receptor gene`, ", ")
co.stimR.func=unique(do.call(rbind, lapply(seq(rec), function(i)cbind(rec[[i]], co.stim$`Immune checkpoint function`[i]))))

files=list.files(path = ".", "subtypes.Rdata")
names(files)=gsub("_subtypes.Rdata", "", files)

# plot subtypes
wrapper.de.analysis=function(i, files){
  load(files[i])
  library(ggplot2)
  
  if(!is.null(dim(coordinates.subtype))){
    subtype=factor(coordinates.subtype$subtype)
    
    coordinates.subtype=coordinates.subtype[order(subtype),]
    subtype=subtype[order(subtype)]
    
    SIZE=ifelse(dim(gexp)[2]>500, 0.5, 2)
    
    # # plot colorv for all:
    # colorv = data.frame("subtype"=c("PML-RARA", "CBFB-MYH11", "Progenitor-like", "MDS-like", "CEBPA", "Monocyte-like","Monocyte-like-MLL", "RUNX1-RUNX1T1"), "color"=c("#aff558", "#dfe4c3", "#d3c684", "#fac75d", "#8a9979", "#d4bcc5", "brown", "#edb2c6"))
    # 
    # colorv=colorv[match(unique(subtype), colorv[,1]),2]
    # 
    # plot.scatter(x=coordinates.subtype$x, y = coordinates.subtype$y, group = subtype, namev = subtype, main = names(files)[i],colorv =as.character(colorv), rasterize = F, width = 70*2, height = 74*2, SIZE = SIZE)
    
    # plot colorv for all:
    # filt=as.character(subtype)%in%c("Ph", "KMT2A", "TCF3-PBX1","Hyperdiploid","ETV6-RUNX1", "PAX5 P80R", "PAX5alt","DUX4", "ZNF384", "MEF2D", "Other")
    # 
    # coordinates.subtype=coordinates.subtype[filt,]
    # subtype=subtype[filt]
      
    plot.scatter(x=coordinates.subtype$x, y = coordinates.subtype$y, group = subtype, namev = subtype, main = names(files)[i], rasterize = F, width = 70*2, height = 74*2, SIZE = SIZE)
    
  }else{
    return(NULL)
  }
}
plot.list=lapply(seq(files), wrapper.de.analysis, files)
plot.list=plot.list[!sapply(plot.list, is.null)]
# 74mm * 99mm per panelwidth = 77, height = 74
figure <-multipanelfigure:: multi_panel_figure(width = 170*2, height = ceiling(length(plot.list)/2)*75*2+175, rows = ceiling(length(plot.list)/2)+1, columns = 2, panel_label_type = "none", unit = "mm", row_spacing = unit(1, "mm"), column_spacing = unit(4, "mm"))

for(i in seq(plot.list)){
  figure <- multipanelfigure::fill_panel(figure,  plot.list[[i]])
}

name="FigS5A_FigS3D_Fig6E_FigS6F_subtypes"
multipanelfigure::save_multi_panel_figure(figure, filename=paste0(name, "_scatterplot.pdf"), limitsize = FALSE)

files=files[!files%in%"PecanALL_subtypes.Rdata"]

#************************************** Differential expression ****************************************

wrapper.de.analysis=function(i, files, genelist){
  load(files[i])
  
  if(!is.null(dim(coordinates.subtype))){
    subtype=coordinates.subtype$subtype
  }else{
    subtype=coordinates.subtype
  }
  
  # make lv of the subtype:
  lv=get.logical(list(subtype))
  
  res=wrapper.wilcoxtest(genelist[genelist%in%rownames(gexp)], data = gexp, logicalVectors = lv, ALTERNATIVE = "two.sided", adj.method = "BH", CORES = 8, prettynum = F)
  
  res$Name=names(files)[i]
  # res=res[res$FDR<0.05,]
  return(res)
}

genelist=co.stim[,1]

res=lapply(seq(files), wrapper.de.analysis, files, genelist)

genelist=co.stimR.f

resR=lapply(seq(files), wrapper.de.analysis, files, genelist)

wrapper.immunoscore.analysis=function(i, files, genesets){
  load(files[i])
  
  if(!is.null(dim(coordinates.subtype))){
    subtype=coordinates.subtype$subtype
  }else{
    subtype=coordinates.subtype
  }
  
  # make lv of the subtype:
  lv=get.logical(list(subtype))
  
  # scores
  gm.objects=do.call(rbind, lapply(seq(genesets), function(i){
    dat3=2^gexp[rownames(gexp)%in%genesets[[i]],,drop=F]+0.01
    gm=log2(t(apply(dat3, 2, gm_mean))) # done to normalized values
    rownames(gm)=names(genesets)[i]
    return(gm)
  }))
  
  
  res=wrapper.wilcoxtest(rownames(gm.objects), data = gm.objects, logicalVectors = lv, ALTERNATIVE = "two.sided", adj.method = "BH", CORES = 8, prettynum = F)
  
  res$Name=names(files)[i]
  # res=res[res$FDR<0.05,]
  return(res)
}
# compute cytolytic and HLA scores
genesets=list(HLAIScore=c("B2M", "HLA-A", "HLA-B","HLA-C"), HLAIIScore=c("HLA-DMA","HLA-DMB","HLA-DPA1","HLA-DPB1","HLA-DRA","HLA-DRB1"), CytolyticScore=c("GZMA", "PRF1", "GNLY", "GZMH", "GZMM"))

res.immunoscores=lapply(seq(files), wrapper.immunoscore.analysis, files, genesets)

# diseases to test:
AML=grep("AML", files)
ALL=grep("ALL", files)
MM=grep("MM", files)
DLBCL=grep("DLBCL", files)

# all results
AML.res=do.call(rbind, res[AML])
ALL.res=do.call(rbind, res[ALL])
ALL.res$Gene[ALL.res$Gene=="VSIR"]="C10orf54"
MM.res=do.call(rbind, res[MM])
DLBCL.res=do.call(rbind, res[DLBCL])

# all results
AML.resR=do.call(rbind, resR[AML])
ALL.resR=do.call(rbind, resR[ALL])
MM.resR=do.call(rbind, resR[MM])
DLBCL.resR=do.call(rbind, resR[DLBCL])

# all results
AML.res.immunoscores=do.call(rbind, res.immunoscores[AML])
ALL.res.immunoscores=do.call(rbind, res.immunoscores[ALL])
MM.res.immunoscores=do.call(rbind, res.immunoscores[MM])
DLBCL.res.immunoscores=do.call(rbind, res.immunoscores[DLBCL])

get.summary.df=function(gene, res.df, stars=T){
  a=sapply(unique(res.df$Group1), function(subtype){
    mean(res.df$FC[res.df$Gene%in%gene&res.df$Group1%in%subtype])
  })
  b=sapply(unique(res.df$Group1), function(subtype){
    survcomp::combine.test(res.df$FDR[res.df$Gene%in%gene&res.df$Group1%in%subtype], method = "z.transform")
  })
  
  pval=-log10(b)
  if(stars){
    pval[b>0.05]=0
    pval[b<0.05]=1
    pval[b<0.01]=2
    pval[b<0.001]=3
    pval[b<1e-5]=4
    pval[b<1e-16]=6
  }
  
  data.frame("variable.1"=log2(a), "variable.2"=pval, "features"=gene, "id"=unique(res.df$Group1))
}

genes.costim=c('CD80', 'CD86', 'CD274', 'PDCD1LG2', 'CD276','C10orf54', 'ICOSLG', 'TNFRSF14', 'PVR','PVRL2', 'PVRL3', 'LGALS9', 'CD48', 'CD58', 'SLAMF6', 'SLAMF7', 'CD84', 'LY9', 'TNFSF4', 'TNFSF8', 'TNFSF9', 'TNFSF18', 'TNFSF15', 'CD70', 'BTN1A1', 'BTN2A2', 'BTN3A1', 'BTNL2', 'BTNL8', 'CD72', 'CD200',  'MICA', 'MICB', 'ULBP1', 'ULBP2', 'ULBP3','RAET1E', 'CLEC2B', 'CLEC2D', 'VTCN1', 'HHLA2', 'IDO1', 'IDO2', 'TDO2', 'NT5E', 'ENTPD1', 'ARG1')

cat(ligand.name[match(genes.costim, co.stim$Gene)], sep="\n")
cat(Receptor.name[match(genes.costim, co.stim$Gene)], sep="\n")
cat(co.stim$`Immune checkpoint function`[match(genes.costim, co.stim$Gene)], sep="\n")

# summarize all the values costim:
df.AML=do.call(rbind, lapply(genes.costim, get.summary.df, AML.res))
df.ALL=do.call(rbind, lapply(genes.costim, get.summary.df, ALL.res))
df.MM=do.call(rbind, lapply(genes.costim, get.summary.df, MM.res))
df.DLBCL=do.call(rbind, lapply(genes.costim, get.summary.df, DLBCL.res))

# all together:
all.subtypes=c(c("MDS-like", "Progenitor-like", "Monocyte-like", "Monocyte-like-MLL", "CEBPA", "RUNX1-RUNX1T1", "CBFB-MYH11", "PML-RARA"),c("Ph", "KMT2A", "TCF3-PBX1","Hyperdiploid","ETV6-RUNX1", "PAX5 P80R", "PAX5alt","DUX4", "ZNF384", "MEF2D", "Other"),  c("ABC", "GCB"), c("Hyperdiploid_gain1q", "Hyperdiploid_gain11q", "CCND1_Ig", "MAF_Ig", "WHSC1_FGFR3_Ig", "TRAF3_Aberrated"))
all.df=rbind(df.AML, df.ALL, df.MM, df.DLBCL)
all.df=all.df[all.df$id%in%all.subtypes,]
all.df=all.df[order(match(as.character(all.df$id), all.subtypes)),]

# summarize all the values, costim R:
df.AML.R=do.call(rbind, lapply(unique(co.stimR.f), get.summary.df, AML.resR))
df.ALL.R=do.call(rbind, lapply(unique(co.stimR.f), get.summary.df, ALL.resR))
df.MM.R=do.call(rbind, lapply(unique(co.stimR.f), get.summary.df, MM.resR))
df.DLBCL.R=do.call(rbind, lapply(unique(co.stimR.f), get.summary.df, DLBCL.resR))

# summarize all the values, immunoscores:
df.AML.immunoscores=do.call(rbind, lapply(rev(unique(AML.res.immunoscores[,1])), get.summary.df, AML.res.immunoscores))
df.ALL.immunoscores=do.call(rbind, lapply(rev(unique(AML.res.immunoscores[,1])), get.summary.df, ALL.res.immunoscores))
df.MM.immunoscores=do.call(rbind, lapply(rev(unique(AML.res.immunoscores[,1])), get.summary.df, MM.res.immunoscores))
df.DLBCL.immunoscores=do.call(rbind, lapply(rev(unique(AML.res.immunoscores[,1])), get.summary.df, DLBCL.res.immunoscores))

all.df.immunoscores=rbind(df.AML.immunoscores, df.ALL.immunoscores, df.MM.immunoscores, df.DLBCL.immunoscores)
all.df.immunoscores=all.df.immunoscores[all.df.immunoscores$id%in%all.subtypes,]
all.df.immunoscores=all.df.immunoscores[order(match(as.character(all.df.immunoscores$id), all.subtypes)),]

all.df$features=factor(all.df$features, levels=unique(all.df$features))
all.df$id=factor(all.df$id, levels=unique(all.df$id))
all.df.immunoscores$features=factor(all.df.immunoscores$features, levels=unique(all.df.immunoscores$features))
all.df.immunoscores$id=factor(all.df.immunoscores$id, levels=unique(all.df.immunoscores$id))

plot.dat=rbind(all.df.immunoscores, all.df)

# Final:
pdf("FigS5B_dotplot.pdf", width = 6, height = 6)
plot.DotPlot.df(data.plot = plot.dat, name.variable.1 = "Fold-Change (log2)", name.variable.2 = "FDR (-log10)", cols = c("blue", "white","red"), col.min = -2, col.max = 2, scale.min = 1, scale.max = 6, dot.scale = 2.5, number.legend.points = 6, fontsize = 7)
dev.off()


df.AML=df.AML[order(-df.AML$variable.2, df.AML$id),]
df.AML$features=factor(df.AML$features, levels = unique(as.character(df.AML$features))[order(df.AML$variable.1)])
df.AML$id=factor(df.AML$id, levels = unique(as.character(df.AML$id)))

# find for each subtype genes that are upregulated/downregulated:
find=ALL.res$FDR<0.05&abs(log2(ALL.res$FC))>0.15
ALL.res=ALL.res[find,]

# find for each subtype genes that are upregulated/downregulated:
find=MM.res$FDR<0.05&abs(log2(MM.res$FC))>0.15
MM.res=MM.res[find,]

# find for each subtype genes that are upregulated/downregulated:
find=DLBCL.res$FDR<0.05&abs(log2(DLBCL.res$FC))>0.15
DLBCL.res=DLBCL.res[find,]

# write tables for supplement:
write.table(AML.res, file = "TableS5_AML_costim_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")
write.table(ALL.res, file = "TableS5_ALL_costim_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")
write.table(MM.res, file = "TableS5_MM_costim_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")
write.table(DLBCL.res, file = "TableS5_DLBCL_costim_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")


# diseases to test:
AML=grep("AML", files)
ALL=grep("ALL", files)
MM=grep("MM", files)
DLBCL=grep("DLBCL", files)

# find for each subtype genes that are upregulated/downregulated:
find=AML.resR$FDR<0.05&abs(log2(AML.resR$FC))>0.15
AML.resR=AML.resR[find,]

# find for each subtype genes that are upregulated/downregulated:
find=ALL.resR$FDR<0.05&abs(log2(ALL.resR$FC))>0.15
ALL.resR=ALL.resR[find,]

# find for each subtype genes that are upregulated/downregulated:
find=MM.resR$FDR<0.05&abs(log2(MM.resR$FC))>0.15
MM.resR=MM.resR[find,]

# find for each subtype genes that are upregulated/downregulated:
find=DLBCL.resR$FDR<0.05&abs(log2(DLBCL.resR$FC))>0.15
DLBCL.resR=DLBCL.resR[find,]

# write tables for supplement:
# write.table(AML.resR, file = "AML_costimR_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")
# write.table(ALL.resR, file = "ALL_costimR_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")
# write.table(MM.resR, file = "MM_costimR_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")
# write.table(DLBCL.resR, file = "DLBCL_costimR_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")



# similar analysis to immunoscores:
# find for each subtype genes that are upregulated/downregulated:
find=AML.res.immunoscores$FDR<0.05&abs(log2(AML.res.immunoscores$FC))>0.15
AML.res.immunoscores=AML.res.immunoscores[find,]

# find for each subtype genes that are upregulated/downregulated:
find=ALL.res.immunoscores$FDR<0.05&abs(log2(ALL.res.immunoscores$FC))>0.15
ALL.res.immunoscores=ALL.res.immunoscores[find,]

# find for each subtype genes that are upregulated/downregulated:
find=MM.res.immunoscores$FDR<0.05&abs(log2(MM.res.immunoscores$FC))>0.15
MM.res.immunoscores=MM.res.immunoscores[find,]

# find for each subtype genes that are upregulated/downregulated:
find=DLBCL.res.immunoscores$FDR<0.05&abs(log2(DLBCL.res.immunoscores$FC))>0.15
DLBCL.res.immunoscores=DLBCL.res.immunoscores[find,]

# write tables for supplement:
# write.table(AML.res.immunoscores, file = "AML_immunoscores_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")
# write.table(ALL.res.immunoscores, file = "ALL_immunoscores_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")
# write.table(MM.res.immunoscores, file = "MM_immunoscores_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")
# write.table(DLBCL.res.immunoscores, file = "DLBCL_immunoscores_DEgenes.txt", col.names = T, row.names = F, quote = F, sep="\t")

extract.res=function(de.genes, geneannot){
  subtypes=de.genes$Group1
  
  de.genes.up=de.genes$FC>1
  de.genes.down=de.genes$FC<1
  
  res=lapply(unique(subtypes), function(type){
    
    counts.up=table(de.genes$Gene[subtypes%in%type&de.genes.up])
    counts.down=table(de.genes$Gene[subtypes%in%type&de.genes.down])
    
    res.up=list()
    res.down=list()
    
    # up
    if(sum(counts.up>1)>0){
      significant=names(counts.up)[counts.up>1]
      funct=geneannot[match(significant, geneannot[,1]),2]
      
      res.up=data.frame(significant, funct, direction="upregulated", stringsAsFactors = F)
      res.up=res.up[order(res.up$funct),]
    }
    
    # down
    if(sum(counts.down>1)>0){
      significant=names(counts.down)[counts.down>1]
      funct=geneannot[match(significant, geneannot[,1]),2]
      
      res.down=data.frame(significant, funct, direction="downregulated", stringsAsFactors = F)
      res.down=res.down[order(res.down$funct),]
    }
    
    res=rbind(res.up, res.down)   
  })
  
  names(res)=unique(subtypes)
  
  return(res)
}

# all results:
data.AML=extract.res(AML.res[AML.res$FDR<0.05&abs(log2(AML.res$FC))>0.15,], co.stim[,c(1,3)])
# sapply(names(data.AML),function (x) write.table(data[[x]], file=paste(x, "_AML_costim.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))

data.ALL=extract.res(ALL.res[ALL.res$FDR<0.05&abs(log2(ALL.res$FC))>0.15,], co.stim[,c(1,3)])
# sapply(names(data.ALL)[!sapply(data.ALL, length)==0],function (x) write.table(data[[x]], file=paste(x, "_ALL_costim.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))

data.MM=extract.res(MM.res, co.stim[,c(1,3)])
# sapply(names(data.MM),function (x) write.table(data[[x]], file=paste(x, "_MM_costim.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))

data.DLBCL=extract.res(DLBCL.res[DLBCL.res$FDR<0.05&abs(log2(DLBCL.res$FC))>0.15,], co.stim[,c(1,3)])
# sapply(names(data.DLBCL),function (x) write.table(data[[x]], file=paste(x, "_DLBCL_costim.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))

data.AML.R=extract.res(AML.resR, co.stimR.func)
# sapply(names(data.AML.R)[names(data.AML.R)%in%"MDS-like"],function (x) write.table(data[[x]], file=paste(x, "_AML_costimR.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))

data.ALL.R=extract.res(ALL.resR[ALL.resR$FDR<0.05&abs(log2(ALL.resR$FC))>0.15,], co.stimR.func)
# sapply(names(data)[-15],function (x) write.table(data[[x]], file=paste(x, "_ALL_costimR.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))

data=extract.res(MM.resR, co.stimR.func)
# sapply(names(data),function (x) write.table(data[[x]], file=paste(x, "_MM_costimR.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))

data.DLBCL.R=extract.res(DLBCL.resR, co.stimR.func)
# sapply(names(data.DLBCL.R),function (x) write.table(data[[x]], file=paste(x, "_DLBCL_costimR.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))

data=extract.res(AML.res.immunoscores, co.stim)
# sapply(names(data),function (x) write.table(data[[x]], file=paste(x, "_AML_immunoscores.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))

data=extract.res(ALL.res.immunoscores, co.stim)
# sapply(names(data),function (x) write.table(data[[x]], file=paste(x, "_ALL_immunoscores.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))

data=extract.res(MM.res.immunoscores, co.stim)
# sapply(names(data),function (x) write.table(data[[x]], file=paste(x, "_MM_immunoscores.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))

data=extract.res(DLBCL.res.immunoscores, co.stim)
# sapply(names(data),function (x) write.table(data[[x]], file=paste(x, "_DLBCL_immunoscores.txt", sep=""), quote = F, sep="\t", col.names = T, row.names = F))


inhibitory=c("CytolyticScore", "HLAIIScore","HLAIScore", co.stim[grepl("Inhibitory|unknown", co.stim$`Immune checkpoint function`),1], co.stimR.func[grepl("Inhibitory|unknown", co.stimR.func[,2]),1])
stimulatory=c("CytolyticScore", "HLAIIScore","HLAIScore", co.stim[grepl("Stimulatory", co.stim$`Immune checkpoint function`),1], co.stimR.func[grepl("Stimulatory", co.stimR.func[,2]),1])
stimulatory=stimulatory[!stimulatory%in%"CTLA4"]

# AML.interesting:
subtype.order=c("MDS-like", "Progenitor-like", "Monocyte-like", "Monocyte-like-MLL", "CEBPA", "RUNX1-RUNX1T1", "CBFB-MYH11", "PML-RARA")
df.AML.pick=rbind(df.AML,df.AML.R, df.AML.immunoscores)
df.AML.pick$features=as.character(df.AML.pick$features)
df.AML.pick$id=as.character(df.AML.pick$id)
interesting.feats=unique(c("CytolyticScore", "HLAIIScore","PVRL3", "SLAMF7","HHLA2","CD274", "TNFRSF14", "ARG1","PVR","CD86", "C10orf54", "CD276","ENTPD1","NT5E", "CD48", "ICOSLG","TNFSF8", "LY9",  "BTN3A1","CD200",  "ULBP1", "LGALS9", "CD70", "MICB", "TNFSF9", "ULBP1", "ULBP3","CLEC2D","TNFSF9","CD84","LY9","CD58","CD72","LGALS9", "CD86", "CD48", "CLEC2B", "LY9", "TNFSF15", "TNFSF8","CD276", "PVR", "BTNL8", "CD58",
                    "LAG3", "SLAMF7","TIGIT","TMIGD2", "PDCD1", "CD160"))

# take all
interesting.feats.L=unlist(sapply(data.AML[subtype.order], function(d)d$significant[d$direction=="upregulated"]))
interesting.feats.R=unlist(sapply(data.AML.R[names(data.AML.R)%in%"MDS-like"], function(d)d$significant[d$direction=="upregulated"]))

interesting.feats=c("CytolyticScore", "HLAIIScore", interesting.feats.L, interesting.feats.R)

df.AML.pick=df.AML.pick[df.AML.pick$features%in%interesting.feats,]
df.AML.pick=df.AML.pick[order(match(df.AML.pick$id, subtype.order)),]
df.AML.pick=df.AML.pick[order(match(df.AML.pick$features, interesting.feats)),]
df.AML.pick$features=factor(df.AML.pick$features, levels=unique(df.AML.pick$features))
df.AML.pick$id=factor(df.AML.pick$id, levels=unique(df.AML.pick$id))

# ALL.interesting:
subtype.order=c("Ph", "KMT2A", "TCF3-PBX1","Hyperdiploid","ETV6-RUNX1", "PAX5 P80R", "PAX5alt","DUX4", "ZNF384", "MEF2D", "Other")
df.ALL.pick=rbind(df.ALL,df.ALL.R, df.ALL.immunoscores)
df.ALL.pick=df.ALL.pick[df.ALL.pick$id%in%subtype.order,]
df.ALL.pick$features=as.character(df.ALL.pick$features)
df.ALL.pick$id=as.character(df.ALL.pick$id)
interesting.feats=unique(c("CytolyticScore", "HLAIIScore","HLAIScore", 
                           "CD200", "NT5E", "TNFRSF14", "BTN3A1", "CLEC2B", "TNFSF4", "TNFSF8",
                           "CLEC2D", "CD58", "CD72","TNFSF8", "TNFSF9",
                           "CLEC2D","CD48", "CD58", "CD72","CD84",
                           "CD200", "CD86", "ICOSLG", "MICA",
                           "BTN2A2", "CD200", "NT5E", "BTN3A1", "TNFSF4",
                           "TNFRSF14", "CD48", "LY9", "TNFSF4",
                           "NT5E", "TNFRSF14", "CLEC2D", "CD48", "CD70", "CLEC2B", "LY9",
                           "BTN2A2", "LGALS9", "TNFRSF14", "CD84", "LY9", "TNFSF4", "TNFSF9",
                           "ENTPD1", "LGALS9", "TNFRSF14", "CD84", "CLEC2B", "TNFSF9",
                           "NT5E", "CLEC2D", "CD48", "CD72", "CD84",
                           "CD274", "CD80"
                           ))

interesting.feats.L=unlist(sapply(data.ALL[subtype.order], function(d)d$significant[d$direction=="upregulated"]))
interesting.feats.R=unlist(sapply(data.ALL.R[subtype.order], function(d)d$significant[d$direction=="upregulated"]))
interesting.feats=c("CytolyticScore", "HLAIIScore","HLAIScore", interesting.feats.L, interesting.feats.R)

df.ALL.pick=df.ALL.pick[df.ALL.pick$features%in%interesting.feats,]
df.ALL.pick=df.ALL.pick[order(match(df.ALL.pick$id, subtype.order)),]
df.ALL.pick=df.ALL.pick[order(match(df.ALL.pick$features, interesting.feats)),]
df.ALL.pick$features=factor(df.ALL.pick$features, levels=unique(df.ALL.pick$features))
df.ALL.pick$id=factor(df.ALL.pick$id, levels=unique(df.ALL.pick$id))

# MM.interesting:
subtype.order=c("Hyperdiploid_gain1q", "Hyperdiploid_gain11q", "MAF_Ig", "WHSC1_FGFR3_Ig", "TRAF3_Aberrated", "CCND1_Ig")
df.MM.pick=rbind(df.MM,df.MM.R, df.MM.immunoscores)
df.MM.pick$features=as.character(df.MM.pick$features)
df.MM.pick$id=as.character(df.MM.pick$id)
interesting.feats=unique(c("HLAIIScore","HLAIScore", "MICA", "MICB", 
                           "CD274","LGALS9", "CD200",  "CD48", "CD86", "ICOSLG", "PVRL3"))
interesting.feats=unlist(sapply(data.MM[subtype.order], function(d)d$significant[d$direction=="upregulated"]))
interesting.feats=c("HLAIIScore","HLAIScore", interesting.feats)

df.MM.pick=df.MM.pick[df.MM.pick$features%in%interesting.feats,]
df.MM.pick=df.MM.pick[order(match(df.MM.pick$id, subtype.order)),]
df.MM.pick=df.MM.pick[order(match(df.MM.pick$features, interesting.feats)),]
df.MM.pick$features=factor(df.MM.pick$features, levels=unique(df.MM.pick$features))
df.MM.pick$id=factor(df.MM.pick$id, levels=unique(df.MM.pick$id))

# DLBCL.interesting:
df.DLBCL.pick=rbind(df.DLBCL,df.DLBCL.R, df.DLBCL.immunoscores)
df.DLBCL.pick$features=as.character(df.DLBCL.pick$features)
df.DLBCL.pick$id=as.character(df.DLBCL.pick$id)
interesting.feats=unique(c("CytolyticScore", "HLAIIScore","HLAIScore","CD274","TNFRSF14", "ENTPD1","SLAMF7",  "ICOSLG", "LY9","CD86",
                           "PDCD1", "BTLA","CD160",  "SLAMF7", "LAG3"))
# take all
interesting.feats.L=unlist(sapply(data.DLBCL[c("ABC", "GCB")], function(d)d$significant[d$direction=="upregulated"]))
interesting.feats.R=unlist(sapply(data.DLBCL.R[c("ABC", "GCB")], function(d)d$significant[d$direction=="upregulated"]))

interesting.feats=c("CytolyticScore", "HLAIIScore", interesting.feats.L, interesting.feats.R)

df.DLBCL.pick=df.DLBCL.pick[df.DLBCL.pick$features%in%interesting.feats&df.DLBCL.pick$id%in%c("ABC", "GCB"),]
df.DLBCL.pick=df.DLBCL.pick[order(match(df.DLBCL.pick$id, c("ABC", "GCB"))),]
df.DLBCL.pick=df.DLBCL.pick[order(match(df.DLBCL.pick$features, interesting.feats)),]
df.DLBCL.pick$features=factor(df.DLBCL.pick$features, levels=unique(df.DLBCL.pick$features))
df.DLBCL.pick$id=factor(df.DLBCL.pick$id, levels=unique(df.DLBCL.pick$id))

# Most interesting shown in main figure, rest in FigS5B. These were significant in at least 2 data sets. Also checked expression in scRNA for various targets.
pdf("Fig5B_AML.pdf", height = 3.5, width = 4.25)
plot.DotPlot.df(data.plot = df.AML.pick[df.AML.pick$features%in%c("CytolyticScore", "HLAIIScore", "PVRL3", "SLAMF7", "HHLA2", "CD274", "TNFRSF14", "ARG1", "PVR", "CD86", "C10orf54", "CD276", "ENTPD1", "NT5E", "CLEC2B", "CD84", "CD48", "LAG3", "TIGIT", "SLAMF7", "TMIGD2", "PDCD1", "CD160", "CD2", "KLRF1"),], name.variable.1 = "Fold-Change (log2)", name.variable.2 = "FDR (-log10)", cols = c("blue", "white","red"), col.min = -2, col.max = 2, scale.min = 1, scale.max = 6, dot.scale = 2.5, number.legend.points = 6, fontsize = 8)
dev.off()

pdf("Fig5B_ALL.pdf", height = 1.5, width = 4.25)
plot.DotPlot.df(data.plot = df.ALL.pick[df.ALL.pick$features%in%c("CD274", "C10orf54", "NT5E"),], name.variable.1 = "Fold-Change (log2)", name.variable.2 = "FDR (-log10)", cols = c("blue", "white","red"), col.min = -2, col.max = 2, scale.min = 1, scale.max = 6, dot.scale = 2.5, number.legend.points = 6, fontsize = 8)
dev.off()

pdf("Fig5B_MM.pdf", height = 2.25, width = 3.8)
plot.DotPlot.df(data.plot = df.MM.pick[df.MM.pick$features%in%c("CytolyticScore", "HLAIIScore","MICA", "MICB","CD274","LGALS9", "CD48", "CD86", "ICOSLG", "PVRL3"),], name.variable.1 = "Fold-Change (log2)", name.variable.2 = "FDR (-log10)", cols = c("blue", "white","red"), col.min = -2, col.max = 2, scale.min = 1, scale.max = 6, dot.scale = 2.5, number.legend.points = 6, fontsize = 8)
dev.off()

pdf("Fig5B_DLBCL.pdf", height = 3.5, width = 3.5)
# all interesting:
plot.DotPlot.df(data.plot = df.DLBCL.pick, name.variable.1 = "Fold-Change (log2)", name.variable.2 = "FDR (-log10)", cols = c("blue", "white","red"), col.min = -2, col.max = 2, scale.min = 1, scale.max = 6, dot.scale = 3, number.legend.points = 6)
plot.DotPlot.df(data.plot = df.DLBCL.pick, name.variable.1 = "Fold-Change (log2)", name.variable.2 = "FDR (-log10)", cols = c("blue", "white","red"), dot.scale = 3)
dev.off()