GIT_HOME="/research/users/ppolonen/git_home/ImmunogenomicLandscape-BloodCancers/"
source(file.path(GIT_HOME, "common_scripts/featurematrix/functions_generate_fm.R"))
source(file.path(GIT_HOME, "common_scripts/visualisation/plotting_functions.R"))
source(file.path(GIT_HOME, "common_scripts/pathway_analysis/functions.GSEA.R"))
library(GSVA)
library(parallel)
# FigureS2 A-B related analysis
setwd("/research/groups/sysgen/PROJECTS/HEMAP_IMMUNOLOGY/petri_work/HEMAP_IMMUNOLOGY/Published_data_figures")
# annotations
annot=get(load("Hemap_immunology_Annotations.Rdata"))
load("/research/groups/sysgen/PROJECTS/HEMAP_IMMUNOLOGY/data/data9544_with_gene_symbols.RData")
data = t(data[rownames(data)%in%annot$GSM.identifier..sample.,])
subclass2=get.logical(annovector = list(annot$subclasses), filterv = annot$CELLS_SORTED==0&annot$Sample.type%in%c("Prolif", "Cancer"))
tbly2=get.logical(annovector = list(annot$tbLY), filterv = annot$CELLS_SORTED==0&annot$Sample.type%in%c("Prolif", "Cancer"))
tbly2=tbly2[!names(tbly2)%in%"Lymphoma_BCL"]
list_cancers=lapply(c(subclass2, tbly2), list)
list_cancers=unlist(list_cancers, recursive=F)
l=c("AML", "BCL_DLBCL", "CLL", "CML", "Lymphoma_BCL_CHL", "Lymphoma_BCL_MCL", "Lymphoma_BCL_FL", "pre-B-ALL", "Prolif_Myeloproliferative_MDS", "T-ALL")
logicalvectors=list_cancers[names(list_cancers)%in%l]
#*********************** run GSEA using cytolyticScore as ranking metric on each of the disease: ***********************
GENESETS=file.path(getwd(), "HALLMARKS.gmt")
WD=file.path(getwd(), "GSEA")
OUTDIR=file.path(getwd(), "GSEA")
names(logicalvectors)=gsub("-", "_", names(logicalvectors))
res=parallel::mclapply(seq(logicalvectors), function(i){
filt=logicalvectors[[i]]
command=run.GSEA(data=data[,filt], cls.vector = as.numeric(annot$CytolyticScore[filt]), datatype = "N", GENESETS = GENESETS, dataname = paste0(names(logicalvectors)[i]), clsname = "CytolyticScore", WD=WD, OUTDIR=OUTDIR)
try(system(command))
}, mc.cores=5)
#***********************************************************************************************************************
#************** run GSEA using cytolyticScore as ranking metric on TNK cells vs other normal ******************
cls=annot$immunoNormals%in%c("CD8+Tcell", "NKCell")
filt=!annot$immunoNormals%in%c("")
command=run.GSEA(data=data[,filt], cls.vector = as.numeric(cls[filt]), datatype = "B", GENESETS = GENESETS, dataname = "CD8_NK_pathways", clsname = "CD8_NK", WD=WD, OUTDIR=OUTDIR)
try(system(command))
#***********************************************************************************************************************
# combine results to same matrix:
datp=list.files(path = OUTDIR, pattern = "gsea_report_for_1_*.*.xls|gsea_report_for_feat_pos.*..xls", recursive = T, full.names = T)
datn=list.files(path = OUTDIR, pattern = "gsea_report_for_0_*.*.xls|gsea_report_for_feat_neg.*..xls", recursive = T, full.names = T)
# filter to contain only significant
pwnames=scan(pipe("cut -f1 HALLMARKS.gmt"), "genesets")
FUN=function(f){
dat=try(read.delim(f, stringsAsFactors=F, header=T), silent = T)
dat$FDR.q.val[dat$FDR.q.val==0]=0.0001
# variable.1, variable.2, features (y name), id (x name)
d=data.frame("features"=gsub("-MSIGDB_HALLMARKS", "", dat$NAME), "variable.1"=dat$NES, "variable.2"=-log10(as.numeric(dat$FDR.q.val)), "id"=gsub("_CytolyticScore_.*.|_pathways_.*.|Prolif_Myeloproliferative_|Lymphoma_BCL_|^/BCL_|/","", gsub(eval(OUTDIR), "", f)))
}
pwp=mclapply(datp, FUN, mc.cores=10)
pwp2=mclapply(datn, FUN, mc.cores=10)
dat=do.call(rbind, pwp)
dat2=do.call(rbind, pwp2)
dat_comb=rbind(dat, dat2)
dat_comb$id=factor(gsub("_", "-", dat_comb$id), levels=c("DLBCL","CHL", "FL", "MCL", "CLL", "CML", "AML", "MDS", "pre-B-ALL", "T-ALL", "CD8-NK"))
pdf("FigureS2A.pdf", width = 5.5, height = 6.5)
plot.DotPlot.df(data.plot = dat_comb, cols = c("blue", "white", "red"), fontsize = 8, name.variable.1 = "NES", name.variable.2 = "-log10 FDR", dot.scale = 3)
dev.off()
# plot correlation plot between protein and Pathway
fm=get(load("TCGA_DLBCL_FM_DUFVA.Rdata"))
# run GSVA for certain pathways and add to matrix:
gexp=data.matrix(fm[grepl("N:GEXP", rownames(fm)),])
rownames(gexp)=gsub(":.*.", "", gsub("N:GEXP:", "",rownames(gexp)))
# gene sets:
GENESETS="Combined_pathway_signatures_2017_filtered_robust.gmt"
# get GSVA visualization for the pathways
# Geneset list
Onc.pathways=read.delim(GENESETS, stringsAsFactors = FALSE, header=F, col.names = paste("V",1:max(count.fields(GENESETS, sep = '\t'), na.rm = T)), fill = TRUE)
# Make list
listA=mclapply(1:length(Onc.pathways[,1]), function(i){A=as.character(Onc.pathways[i,3:length(Onc.pathways),])
B=A[!A==""&!A=="NA"]}, mc.cores=6)
names(listA) <- toupper(Onc.pathways[,1])
viz_scores=gsva(expr = gexp, gset.idx.list = listA, parallel.sz=8, tau=0.25)
rownames(viz_scores)=paste0("N:GSVA:", rownames(viz_scores))
fm.add=rbind(fm, viz_scores)
fm.add=fm.add[!(!grepl("^B:GNAB:*.*chr*.*y_n_somatic", rownames(fm.add))&grepl("^B:GNAB:", rownames(fm.add))),!is.na(fm.add["N:RPPA:PARP1:chr1:226548392:226595780:-:PARP_cleaved",])]
annotdf=data.frame("CytolyticScore"=as.numeric(fm["N:SAMP:CytolyticScore",]), "CASP7"=as.numeric(fm["N:RPPA:CASP7:chr10:115438942:115490662:+:Caspase-7_cleavedD198",]), "IRF1"=as.numeric(fm["N:RPPA:IRF1:chr5:131817301:131826490:-:IRF-1",]))
rownames(annotdf)=colnames(fm)
annotdf=annotdf[!is.na(annotdf$CASP7),]
pdf("FigureS2B.pdf", height = 2.5, width = 2.5)
ggscatter(annotdf, x = "CytolyticScore", y = "CASP7",
size = 1.5,
add = "reg.line", # Add regression line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE # Add confidence interval
) +
stat_cor(method = "pearson", label.sep = "\n") +
xlab("Cytolytic Score") +
ylab("CASP7")
ggscatter(annotdf, x = "CytolyticScore", y = "IRF1",
size = 1.5,
add = "reg.line", # Add regression line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE # Add confidence interval
) +
stat_cor(method = "pearson", label.sep = "\n") +
xlab("Cytolytic Score") +
ylab("IRF1")
dev.off()
Datasets
Models
