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

Models:
WandaB/
ImmunogenomicLandscape-BC
Downloads: 1
[8e0848]: / Statistical_analysis_Szabo_TCell_analysis.R
History
Download this file
170 lines (133 with data), 7.5 kB

GIT_HOME="/research/users/ppolonen/git_home/ImmunogenomicLandscape-BloodCancers/"
source(file.path(GIT_HOME, "common_scripts/scRNA/functions.scRNA.analysis.R"))
source(file.path(GIT_HOME, "common_scripts/visualisation/plotting_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")

#****************************** plot Szabo scores and immunoscores, singleR etc ******************************

a=fread("Szabo_Tcell_markers.txt", data.table = F)
rownames(a)=make.unique(a[,1])
dat=data.matrix(a[,-1])
dat[is.na(dat)]=0
# plot markers in

# compare NK cells
geneset=data.frame("genes"=a[,1], "name"=NA, stringsAsFactors = F)
geneset$name[1:70]="Treg"
geneset$name[71:140]="CD4 Tn / Tcm resting"
geneset$name[141:210]="CD4 / CD8 resting"
geneset$name[211:280]="IFN response"
geneset$name[281:350]="Proliferation"
geneset$name[351:420]="CD8 Cytotoxic"
geneset$name[421:490]="CD8 Cytokine"

geneset2=data.frame("genes"=a[,1], "name"=NA, stringsAsFactors = F)
geneset2$name[1:70]="Resting"
geneset2$name[71:140]="Resting"
geneset2$name[141:210]="Resting"
geneset2$name[211:280]="Activated"
geneset2$name[281:350]="Activated"
geneset2$name[351:420]="Activated"
geneset2$name[421:490]="Activated"

geneset.l=lapply(unique(geneset[,2]), function(g)geneset[geneset[,2]%in%g,1])
names(geneset.l)=unique(geneset[,2])

geneset.l2=lapply(unique(geneset2[,2]), function(g)geneset2[geneset2[,2]%in%g,1])
names(geneset.l2)=unique(geneset2[,2])

geneset.exhaustion=c("PDCD1", "CTLA4", "LAYN", "LAG3", "HAVCR2", "CD244", "CD160")
geneset.cytokine=c("CCL3", "CCL4", "XCL1", "XCL2", "IFNG", "CSF2", "IL10","TNFRSF9")
geneset.cytotoxic=c("CCL5", "GZMK", "GNLY", "EOMES", "ZEB2", "ZNF683", "KLRG1", "NKG7")

geneset.l3=list("Exhaustion"=geneset.exhaustion, "Cytokine"=geneset.cytokine, "Cytotoxic"=geneset.cytotoxic)

# add gm based score:
add.scores=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"))
add.scores=append(append(geneset.l, append(geneset.l2, add.scores)), geneset.l3)



load("Szabo_Tcell_BM_scRNA.Rdata")
szabo=scmat
szabo <- NormalizeData(szabo, assay = "RNA", normalization.method = "LogNormalize", scale.factor = 10000)
fm.szabo <- data.matrix(GetAssayData(szabo, assay = "RNA", slot="data"))

# same in FIMM data:
load("FIMM_AML_HCA_T_scRNA.Rdata")

scmat=FindClusters(scmat, resolution = 0.5)

# make sure that using logNormalized data and scaled to 10000:
scmat <- NormalizeData(scmat, assay = "RNA", normalization.method = "LogNormalize", scale.factor = 10000)
fm.f <- data.matrix(GetAssayData(scmat, assay = "RNA", slot="data"))

sample=gsub("_.*.", "", colnames(scmat))
sample[grepl("Manton", sample)]="HCA BM"
sample[!grepl("3667|5249|5897", sample)&!grepl("HCA", sample)]="AML other"
scmat[["sample"]]=sample

# take cluster assignment for AML cells:
scmat2=FindNeighbors(scmat[,grepl("3667|5897|5249", sample)], dims = 1:25, reduction = "mnnCorrect")
scmat2=FindClusters(scmat2, resolution = 0.5, algorithm = 1)

AML_cluster_Tcell=Idents(scmat2)
save(AML_cluster_Tcell, file="Szabo_HCA_AML_cluster_cell.Rdata")

pdf("Figure_S3M_FIMM_AML_reclustered_Tcell.pdf", width = 6, height = 5)
DimPlot(scmat, pt.size = 0.1, cells.highlight = colnames(scmat2)[Idents(scmat2)%in%0], cols.highlight = "#006f00") # cytotoxic
DimPlot(scmat, pt.size = 0.1, cells.highlight = colnames(scmat2)[Idents(scmat2)%in%3], cols.highlight = "#006f00") # Cytokine
DimPlot(scmat2, label = T) # all clusters
dev.off()

geneset.l=lapply(unique(geneset[,2]), function(g){
  genes=geneset[geneset[,2]%in%g,1]
  
  #take only significant:
  genes=genes[genes%in%markers.all$gene[grepl("CD|Treg", markers.all$cluster)]]
})
names(geneset.l)=paste(unique(geneset[,2]), "filtered")

# add.scores=append(add.scores, geneset.l)

# add gm based score:
gm.objects=do.call(rbind, lapply(seq(add.scores), function(i){
  dat3=fm.f[rownames(fm.f)%in%add.scores[[i]],]
  gm=log2(t(apply(dat3, 2, gm_mean))) # done to normalized values
  rownames(gm)=names(add.scores)[i]
  return(gm)
}))

# also add to seurat object:
for(i in seq(add.scores)){
  scmat[[names(add.scores)[i]]] <- gm.objects[i,]
}

# add gm based score:
gm.objects=do.call(rbind, lapply(seq(add.scores), function(i){
  dat3=fm.szabo[rownames(fm.szabo)%in%add.scores[[i]],]
  gm=log2(t(apply(dat3, 2, gm_mean))) # done to normalized values
  rownames(gm)=names(add.scores)[i]
  return(gm)
}))

# also add to seurat object:
for(i in seq(add.scores)){
  szabo[[names(add.scores)[i]]] <- gm.objects[i,]
}

# pdf("FigureS3M_Szabo_markers_FIMM_AML_HCA_Tcells_Scores.pdf", width=16, height=10)
# FeaturePlot(scmat, features = names(add.scores), cols=c("grey75", "red"), min.cutoff = 0.1)
# FeaturePlot(szabo, features = names(add.scores), cols=c("grey75", "red"), min.cutoff = 0.1)
# dev.off()

pdf("Figure3M_Szabo_markers_FIMM_AML_HCA_Tcells_Scores_scaled.pdf", width=16, height=14)
p1=FeaturePlot(scmat, features = names(add.scores), combine = F, min.cutoff = 0.2, max.cutoff = 1)
fix.sc <- scale_color_gradientn( colours = c('grey75', 'red'),  limits = c(0.2, 1))
p2 <- lapply(p1, function (x) x + fix.sc)
CombinePlots(p2)

p1=FeaturePlot(szabo, features = names(add.scores), combine = F, min.cutoff = 0.2, max.cutoff = 1)
fix.sc <- scale_color_gradientn( colours = c('grey75', 'red'),  limits = c(0.2, 1))
p2 <- lapply(p1, function (x) x + fix.sc)
CombinePlots(p2)
dev.off()

pdf("Figure3I_FigureS3M_Szabo_markers_FIMM_AML_HCA_Tcells_Scores_blend.pdf", width=20, height=5)
FeaturePlot(scmat, features = names(add.scores)[c(14,15)], blend.threshold = 0, blend = T,cols = c("grey74", "darkblue", "red"))
FeaturePlot(scmat, features = names(add.scores)[c(8,9)], blend = T, blend.threshold = 0.6, cols = c("grey74", "darkblue", "red"))
dev.off()

colors.group=data.table::fread("colors_lineage.txt", data.table = F, header = F)
colors.group=colors.group[grepl("CD|Treg", colors.group$V1),]
plot.multi.scatter.scRNA(data.frame("Sample"="T cells HCA AML", "Data"="FIMM_AML_HCA_T_scRNA.Rdata", "Clusters"="", "Type"="", stringsAsFactors = F), colors.group = colors.group, seurat.feature = "SingleR.label", name="FigureS3M_singleR", cores=9, SIZE = 0.1, rasterize=F, width = 77*4, height = 74*4, text.size = 10*4)

# plot MDS-like samples to map:
scmat[["MDSlike"]]=ifelse(grepl("5897|3667|5249", scmat[["batch"]][,1]), "MDS-like", "other")
scmat[["MDSlike"]][grepl("BM", scmat[["batch"]][,1]),1]="normal BM"

lv = scmat[["MDSlike"]][,1]%in%c("MDS-like", "other")
# plot.scatter.seurat(scmat = scmat, colors.group = data.frame("V1"=c("normal BM", "other", "MDS-like"), "V2"=c("grey75", "orange", "#21a366"), stringsAsFactors = F), seurat.feature = "MDSlike", name="T cells (AML, HCA)",lv = lv, cores=1, add.density = T, add.proportions = T, SIZE = 1, rasterize=F, width = 64*4, height = 74*4, text.size=10*4)

xy=data.frame(Embeddings(scmat, "umap"))

pdf("Figure3I_Density_MDSlike_Tcell.pdf", width = 5, height = 4)
ggplot(xy[lv,], aes(x=UMAP_1, y=UMAP_2) ) +
stat_density2d(aes(fill = ..density..), contour = F, geom = 'tile') +
  viridis::scale_fill_viridis()
ggplot(xy[scmat[["MDSlike"]]=="normal BM",], aes(x=UMAP_1, y=UMAP_2) ) +
  stat_density2d(aes(fill = ..density..), contour = F, geom = 'tile') +
  viridis::scale_fill_viridis()
dev.off()