Upload New File

scripts/5_Analysis_TAMS_harmony.R

+library(Seurat)
+library(harmony)
+library(ggplot2)
+library(openxlsx)
+library(RColorBrewer)
+library(fmsb)
+library(matrixStats)
+
+wdir <- "squadrito_livertumor2022_scrnaseq"
+out_dir <- paste(wdir, "results", sep = "/")
+# Full dir
+full_dir <- paste(out_dir, "Full", sep = "/")
+full_obj <- readRDS(paste(full_dir, "Full_final_DBR_labeled_Und_rem.rds", sep = "/"))
+Idents(full_obj) <- "Newlabels_detailed"
+tams_obj_min <- subset(full_obj, idents = c("TAMs", "IFNa TAMs", "KCs"))
+
+#Normalize data
+tams_obj_min <- NormalizeData(tams_obj_min)
+
+# Find most variable genes (Feature selection) using SB code
+tams_obj_min <- FindVariableFeatures(tams_obj_min, selection.method = "vst", nfeatures = 3000, verbose = T)
+
+# Cell-Cycle scoring
+s.genes <- c("Mcm5","Pcna","Tyms","Fen1","Mcm2","Mcm4","Rrm1","Ung","Gins2","Mcm6","Cdca7","Dtl","Prim1","Uhrf1","E2f8",
+             "Cenpu","Hells","Rfc2","Rpa2","Nasp","Rad51ap1","Gmnn","Wdr76","Slbp","Ccne2","Ubr7","Pold3","Msh2","Atad2",
+             "Rad51","Rrm2","Cdc45","Cdc6","Exo1","Tipin","Dscc1","Blm","Casp8ap2","Usp1","Clspn","Pola1","Chaf1b","Brip1")
+g2m.genes <- c("Hmgb2","Cdk1","Nusap1","Ube2c","Birc5","Tpx2","Top2a","Ndc80","Cks2","Nuf2","Cks1b","Mki67","Tmpo","Cenpf","Tacc3",
+               "Pimreg","Smc4","Ccnb2","Ckap2l","Ckap2","Aurkb","Bub1","Kif11","Anp32e","Tubb4b","Gtse1","Kif20b","Hjurp","Cdca3","Notch1",
+               "Cdc20","Ttk","Cdc25c","Kif2c","Rangap1","Ncapd2","Dlgap5","Cdca2","Cdca8","Ect2","Kif23","Hmmr","Aurka","Psrc1","Anln",
+               "Lbr","Ckap5","Cenpe","Ctcf","Nek2","G2e3","Gas2l3","Cbx5","Cenpa")
+tams_obj_min <- CellCycleScoring(object = tams_obj_min, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
+tams_obj_min@meta.data$CC.Difference <- tams_obj_min@meta.data$S.Score - tams_obj_min@meta.data$G2M.Score
+
+# Scale Data, SCTransform was not used, but it is an alternative, TAKES 2 hours!!
+reg_vars <- c("CC.Difference", "percent.mt", "nCount_RNA")
+tams_obj_min <- ScaleData(object = tams_obj_min, vars.to.regress = reg_vars,  display.progress = T, features = rownames(tams_obj_min))
+
+##Run PCA
+tams_obj_min <- RunPCA(tams_obj_min, npcs = 50)
+#Determine number of PCs
+ElbowPlot(tams_obj_min, ndims = 50)
+
+integration.var <- "orig.ident"
+prefix <- "RNA"
+# Harmony
+cat("Harmony: tSNE & UMAP", "\n")
+tams_obj_min <- RunHarmony(object = tams_obj_min, group.by.vars = integration.var, plot_convergence = FALSE,
+                           reduction.save = paste0("harmony.",integration.var))
+tams_obj_min <- RunUMAP(object = tams_obj_min, seed.use = 123, reduction = paste0("harmony.",integration.var), dims = 1:25,
+                        reduction.name = paste0("umap.harmony.",integration.var), reduction.key = "UMAPh_", return.model = TRUE)
+
+tams_obj_min <- RunTSNE(object = tams_obj_min, do.fast = T, seed.use = 123, reduction = paste0("harmony.",integration.var), perplexity = 30,
+                        dims = 1:25, reduction.name = paste0("tsne.harmony.",integration.var), reduction.key = "tSNEh_", return.model=TRUE)
+tams_obj_min <- FindNeighbors(object = tams_obj_min, reduction = paste0("harmony.",integration.var), dims = 1:25, force.recalc = T,
+                              graph.name = c(paste0(prefix,"_nn_h.",integration.var),paste0(prefix,"_snn_h.",integration.var)))
+
+tams_obj_min <- FindClusters(object = tams_obj_min, algorithm = 1, resolution = 0.6,
+                             graph.name = paste0(prefix,"_snn_h.",integration.var))
+
+liver_moccart <- "DEGs_Macrophages.xlsx"
+liver_signatures <- list()
+for (sn in grep("Mouse", grep("DEGs", openxlsx::getSheetNames(liver_moccart), value = T), value = T)) {
+    sname <- gsub(")" , "", gsub("[(&.]", "_", gsub(" ", "", gsub(" DEGs", "", gsub("Mouse ", "", sn)))))
+    t <- read.xlsx(xlsxFile = liver_moccart, sheet = sn)
+    t <- head(x = t, n = 50)
+    liver_signatures[[sname]] <- t$X1
+    tams_obj_min <- AddModuleScore(object = tams_obj_min, features = list(liver_signatures[[sname]]), name = sname)
+}
+
+Idents(tams_obj_min) <- "RNA_snn_h.orig.ident_res.0.6"
+marks <- FindAllMarkers(object = tams_obj_min, logfc.threshold = 0, min.pct = 0.1, return.thresh = 1)
+write.xlsx(x = list("ClusterMarkers" = marks),
+           file = paste(out_dir, "TAMS_ClusterMarkersFull.xlsx", sep = "/"))
+
+tams_obj_min@misc[["markers_h"]] <- marks
+
+saveRDS(object = tams_obj_min, file = paste(out_dir, "TAMS_final.rds", sep = "/"))
+
+
+APC_obj<-readRDS(paste(out_dir, "TAMS_final.rds", sep = "/"))
+APC_obj@meta.data$TAM_clusters<- APC_obj@meta.data$seurat_clusters
+APC_obj@meta.data$TAM_clusters[APC_obj@meta.data$seurat_clusters%in%c(6,7,8)] <- 6
+
+quartz(width=15,height=7)
+DimPlot(APC_obj, reduction = "umap.harmony.orig.ident", group.by = "TAM_clusters", label = T, split.by = "RNA_Group",pt.size = 1.2)
+
+TamGenes<- c("Mmp8", "Chil3", "Fn1", "Anxa1","Sell",
+             "Ace","Adgre4","Rgs2","Sat1",
+             "Cxcl10","Ifi205","Cd40","Marcksl1",
+             "Ceacam1","Ear2","Fabp4","Klf4",
+             "Ccl24","Arg1","Dab2","Mif","Cd274",
+             "Mrc1","C1qa","C1qb","Clec4a2","Apoe")
+
+
+DotPlot(APC_obj,features = TamGenes,group.by = "TAM_clusters",
+        dot.scale=14)+ theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
+
+
+#Differential expressed genes in macrophages, NOt Used in the paper
+All.Markers <- NULL
+for (i in 0:5){
+    APC_obj<-SetIdent(APC_obj, value = "TAM_clusters")
+    obj.i<-subset(APC_obj,idents = i)
+    obj.i<-SetIdent(obj.i, value = "RNA_Group")
+    Markers.i<- FindAllMarkers(obj.i,min.pct = 0.3,only.pos = T)
+    Markers.i$obj <- paste("Cluster",i,sep="_")
+    All.Markers <- rbind(All.Markers,Markers.i)
+}
+All.Markers$dpct <- All.Markers$pct.1-All.Markers$pct.2
+All.Markers<- All.Markers[rev(order(All.Markers$dpct)),]
+
+for ( i in unique(All.Markers$obj)){
+    a <- unique(All.Markers[All.Markers$obj==i & All.Markers$cluster=="CTRL","gene"])[1:9]
+    b <- unique(All.Markers[All.Markers$obj==i & All.Markers$cluster=="IFN_PR","gene"])[1:30]
+    b <-b[!a%in%b][1:9]
+    c <- unique(All.Markers[All.Markers$obj==i & All.Markers$cluster=="IFN_R","gene"])[1:30]
+    c <-c[!c%in%c(a,b)][1:9]
+    selected.i <- c(a,b,c)
+    assign(value =selected.i,x=paste("Markers",i,sep=".") )
+}
+Markers.Cluster_5
+
+obj.i<-subset(APC_obj,idents = 5)
+#quartz()
+DotPlot(obj.i,features =Markers.Cluster_5,group.by = "RNA_Group",
+        dot.scale=15)+ theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
+
+obj.i<-subset(APC_obj,idents = 2)
+quartz()
+DotPlot(obj.i,features =Markers.Cluster_2,group.by = "RNA_Group",
+        dot.scale=15)+ theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
+
+
+#Other try with genes we like
+obj.i<-subset(APC_obj,idents = 1)
+MHCI_genes<- c("H2-T22", "H2-T23", "H2-D1","B2m", "Tap1", "Tap2", "Tapbp", "Psmb8", "Psmb9", "Cd86")
+MHCII_genes<- c("H2-Aa","H2-Ab1", "H2-Eb1","H2-DMb1","H2-Oa", "Ciita", "Cd74", "Cd40")
+IL10_genes<- c("Tgfb1", "Cebpb", "Il4ra", "Socs3", "Ccl24")
+CTRL_genes<- c("Mmp8", "Trem2", "Fn1")
+IFNa_genes<-c("Stat1", "Socs1", "Nlrc5")#
+quartz()
+DotPlot(obj.i,features = c(MHCI_genes,MHCII_genes,IL10_genes, CTRL_genes, IFNa_genes),group.by = "RNA_Group",
+        dot.scale=15)+ theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
+
+######
+##Figure Spider plot.. spider plot FIGURE S5 of revision
+#Add module score with selected features
+MHCI_genes<- c("H2-T22", "H2-T23", "H2-D1","B2m", "Tap1", "Tap2", "Tapbp", "Psmb8", "Psmb9", "Cd86")
+MHCI_genes<-list(MHCI_genes)
+MHCII_genes<- c("H2-Aa","H2-Ab1", "H2-Eb1","H2-DMb1","H2-Oa", "Ciita", "Cd74", "Cd40")
+MHCII_genes<-list(MHCII_genes)
+IL10_genes<- c("Tgfb1", "Cebpb", "Il4ra", "Socs3", "Ccl24")
+IL10_genes<-list(IL10_genes)
+CTRL_genes<- c("Mmp8", "Trem2", "Fn1")
+CTRL_genes<-list(CTRL_genes)
+IFNa_genes<-c("Stat1", "Socs1", "Nlrc5", "Irf7", "Ifit1", "Oas1a")
+IFNa_genes<-list(IFNa_genes)
+
+APC_obj<- AddModuleScore(APC_obj, features = MHCI_genes, name = "MHCI_genes")
+APC_obj<- AddModuleScore(APC_obj, features = MHCII_genes, name = "MHCII_genes")
+APC_obj<- AddModuleScore(APC_obj, features = IL10_genes, name = "IL10_genes")
+APC_obj<- AddModuleScore(APC_obj, features = CTRL_genes, name = "CTRL_genes")
+APC_obj<- AddModuleScore(APC_obj, features = IFNa_genes, name = "IFNa_genes")
+
+#<Normalize data
+#we needed to find a way to normalize the module scores. 
+#we used a quantile normalization method.
+M1<-APC_obj@meta.data[,c(4, 25:30)]
+M1[,3:7]<-limma::normalizeQuantiles(M1[,3:7])
+Groups <-unique(APC_obj$RNA_Group)
+
+
+#with Median
+#Cluster0
+M2<-M1[M1$TAM_clusters%in%0,]
+M3a<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[1],3:7]),na.rm=T)
+M3b<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[2],3:7]),na.rm=T)
+M3c<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[3],3:7]),na.rm=T)
+M4<-rbind(M3a,M3b,M3c)
+rownames(M4)<-Groups
+colnames(M4)<-colnames(M2[,3:7])
+C0.db<-as.data.frame(M4)
+
+#Cluster1
+M2<-M1[M1$TAM_clusters%in%1,]
+M3a<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[1],3:7]),na.rm=T)
+M3b<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[2],3:7]),na.rm=T)
+M3c<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[3],3:7]),na.rm=T)
+M4<-rbind(M3a,M3b,M3c)
+rownames(M4)<-Groups
+colnames(M4)<-colnames(M2[,3:7])
+C1.db<-as.data.frame(M4)
+
+#Cluster2
+M2<-M1[M1$TAM_clusters%in%2,]
+M3a<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[1],3:7]),na.rm=T)
+M3b<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[2],3:7]),na.rm=T)
+M3c<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[3],3:7]),na.rm=T)
+M4<-rbind(M3a,M3b,M3c)
+rownames(M4)<-Groups
+colnames(M4)<-colnames(M2[,3:7])
+C2.db<-as.data.frame(M4)
+
+#Cluster3
+M2<-M1[M1$TAM_clusters%in%3,]
+M3a<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[1],3:7]),na.rm=T)
+M3b<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[2],3:7]),na.rm=T)
+M3c<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[3],3:7]),na.rm=T)
+M4<-rbind(M3a,M3b,M3c)
+rownames(M4)<-Groups
+colnames(M4)<-colnames(M2[,3:7])
+C3.db<-as.data.frame(M4)
+
+#Cluster4
+M2<-M1[M1$TAM_clusters%in%4,]
+M3a<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[1],3:7]),na.rm=T)
+M3b<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[2],3:7]),na.rm=T)
+M3c<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[3],3:7]),na.rm=T)
+M4<-rbind(M3a,M3b,M3c)
+rownames(M4)<-Groups
+colnames(M4)<-colnames(M2[,3:7])
+C4.db<-as.data.frame(M4)
+
+#Cluster5
+M2<-M1[M1$TAM_clusters%in%5,]
+M3a<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[1],3:7]),na.rm=T)
+M3b<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[2],3:7]),na.rm=T)
+M3c<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[3],3:7]),na.rm=T)
+M4<-rbind(M3a,M3b,M3c)
+rownames(M4)<-Groups
+colnames(M4)<-colnames(M2[,3:7])
+C5.db<-as.data.frame(M4)
+
+
+MinMax<-rbind(C0.db, C1.db, C2.db, C3.db, C4.db, C5.db)
+C0.db<-rbind(rep(max(MinMax),5),rep(min(MinMax),5),C0.db)
+C1.db<-rbind(rep(max(MinMax),5),rep(min(MinMax),5),C1.db)
+C2.db<-rbind(rep(max(MinMax),5),rep(min(MinMax),5),C2.db)
+C3.db<-rbind(rep(max(MinMax),5),rep(min(MinMax),5),C3.db)
+C4.db<-rbind(rep(max(MinMax),5),rep(min(MinMax),5),C4.db)
+C5.db<-rbind(rep(max(MinMax),5),rep(min(MinMax),5),C5.db)
+
+colors_border=c( rgb(1,0,0,0.9), rgb(0,0.5,0,0.9) , rgb(0,0,1,0.9) )
+colors_in=c( rgb(1,0,0,0.1), rgb(0,0.5,0,0.1) , rgb(0,0,1,0.1) )
+
+quartz()
+par(mfrow=c(2,3))
+radarchart(C0.db, title = "Cluster 0", pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty = 1,cglwd=0.8)
+radarchart(C1.db, title = "Cluster 1", pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty = 1,cglwd=0.8)
+radarchart(C2.db, title = "Cluster 2", pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty = 1,cglwd=0.8)
+radarchart(C3.db, title = "Cluster 3", pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty = 1,cglwd=0.8)
+radarchart(C4.db, title = "Cluster 4", pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty = 1,cglwd=0.8)
+radarchart(C5.db, title = "Cluster 5", pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty = 1,cglwd=0.8)
+