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"))

library(Matrix)
library(Seurat)
library(data.table)
library(ComplexHeatmap)
library(circlize)
library(parallel)
library(ggplot2)

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.stimR=data.table::fread("costim_ligands_final.txt", data.table = F)
co.stimR.f=c(unlist(strsplit(co.stimR$`Receptor gene`, ", ")), "LAG3")
co.stim=rbind(co.stim, c("FGL1", "LAG3", "Inhibitory"))

load("AML_Galen_scRNA.Rdata")
galen=scmat

load("HCA_scRNA.Rdata")
HCA=scmat

load("FIMM_AML_scRNA.Rdata")

# plot from each dataset all significant genes:
subtype.order=c("MDS-like", "Progenitor-like", "Monocyte-like", "Monocyte-like-MLL", "CEBPA", "RUNX1-RUNX1T1", "CBFB-MYH11", "PML-RARA")

fimm=scmat[,scmat[["SingleR.label"]][,1]%in%c("HSC", "MPP", "GMP", "CMP", "MEP", "Monocytes", "Erythrocytes")]
galen.sub=galen[,galen[["SingleR.label"]][,1]%in%c("HSC", "MPP", "GMP", "CMP", "MEP", "Monocytes", "Erythrocytes")]
hca=HCA[,HCA[["SingleR.label"]][,1]%in%c("HSC", "MPP", "GMP", "CMP", "MEP", "Monocytes", "Erythrocytes")]

Idents(fimm)=fimm[["SingleR.label"]][,1]
Idents(galen.sub)=factor(galen.sub[["SingleR.label"]][,1], levels=c("HSC", "MPP", "GMP", "CMP", "MEP", "Monocytes", "Erythrocytes"))
Idents(hca)=hca[["SingleR.label"]][,1]

# do DE gene tests:
DE.FIMM=FindAllMarkers(object = fimm, features = genelist[genelist%in%rownames(scmat)], only.pos = T, logfc.threshold = 0.15)
DE.galen=FindAllMarkers(object = galen.sub, features = genelist[genelist%in%rownames(galen)], only.pos = T, logfc.threshold = 0.15)
DE.HCA=FindAllMarkers(object = hca, features = genelist[genelist%in%rownames(HCA)], only.pos = T, logfc.threshold = 0.15)

all=rbind(DE.FIMM, DE.galen, DE.HCA)
all=all[order(match(all$cluster, c("HSC", "MPP", "GMP", "CMP", "MEP", "Monocytes", "Erythrocytes"))),]

a=table(all$cluster, all$gene)
genelist.filt=c(colnames(a[,apply(a, 2, function(v)any(v>1))]), "C10orf54")

# these were selected for the figure:
genelist.filt=c("CD34", "CLEC2B", "TNFRSF14", "CD84", "VSIR", "C10orf54", "CD68", "CD48","CD86","ENTPD1") 

pdf("FigS5C_VISTA_FIMM_dotplot.pdf", width = 4.5, height = 2.5)
plot.DotPlot(fimm, features = genelist.filt[genelist.filt%in%rownames(scmat)], cols=c("grey75", "red"), scale.max = 50, scale.min = 10, dot.scale = 4)
plot.DotPlot(galen.sub, features = genelist.filt[genelist.filt%in%rownames(galen)], cols=c("grey75", "red"), scale.max = 50, scale.min = 10, dot.scale = 4)
plot.DotPlot(hca, features = genelist.filt[genelist.filt%in%rownames(HCA)], cols=c("grey75", "red"), scale.max = 50, scale.min = 10, dot.scale = 4)
dev.off()