library(Seurat) library(openxlsx) library(readxl) library(tidyr) library("RColorBrewer") library(dplyr) library(reshape2) library(fgsea) library(matrixStats) library(fmsb) cluster_dir <- "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_NotaroTAD_A019/Analysis_MM" wdir <- paste0(cluster_dir, "/Result3_MNOT") ddir <- paste0(cluster_dir, "/results3/Liver_Tumor_APCs/") pdir <- paste0(cluster_dir, "/Result3_MNOT/Plots") # Custom functions source(file = "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_NotaroTAD_A018/Analysis_MM/last_scripts_MNOT/TAD018/plot_stellare_function.R") # Load RDS obs1 <- readRDS(paste0(ddir, "Liver_Tumor_APCs_final.rds")) ######################## #####Annotations####### ######################## obs1@meta.data$Tissue <- obs1@meta.data$PROP.Source obs1@meta.data$RNA_Group <- obs1@meta.data$PROP.Condition obs1@meta.data$RNA_Group[obs1@meta.data$RNA_Group == "TA_only"] <- "TA33" obs1@meta.data$RNA_Group[obs1@meta.data$RNA_Group == "TA_combo"] <- "TA33.Combo" obs1@meta.data$RNA_Group <- factor(x = obs1@meta.data$RNA_Group, levels = c("TA33","TA33.Combo")) # Tissue Sample assign obs1$Tube <- paste0(obs1$hash.ID,obs1$RNA_Group) obs1@meta.data$Orig<-paste0(obs1@meta.data$Tissue, ".", obs1@meta.data$RNA_Group) table(obs1@meta.data$Tissue,obs1@meta.data$RNA_Group) table(obs1@meta.data$Orig) obs1@meta.data$allData <- "all" #############Analysis APCs################### features1 = c("H2-Q7", "H2-T22", "H2-Q4", "H2-K1", "H2-T23", "H2-M3", "H2-D1","Tap1","Tap2", #MHCI "H2-Eb1", "H2-Ab1", "H2-Aa", "H2-DMb1", "H2-DMa", "Ciita", #MHCII "Ifi44", "Oasl1","Irf7","Isg15", "Isg20","Oas1a","Oas1g", #IFNa "Ccl5", "Upp1", "Slamf7", "Cxcl9", "Gbp4","Cd74", #IFNg "Vav2","Tgfb1","Il10","Ccl24","Mmp8","Tmem176b","Trem2","Fn1") #protumor #####Figure 6A############## pdf(paste0(pdir, "/liver.apc.19.LandT_Fig6A.pdf"), width=14, height=3) plot.a <- PlotStellare(obs1, features1=features1, valueG="allData", i="all", group.by0="Orig", controGroup="Liver.TA33") print(plot.a) par(mfrow=c(1,2)) dev.off() #######################GSEA####################################################################### #######################Figure S6B - gsea on Liver and tumor APCs################################## # Load signature #miDB_sig4.MLS <- readRDS("C:/Users/notaro.marco/Dropbox (HSR Global)/90-756142658_RNAseq_Squadrito/05-DGE-Tumor/miDB_sig4.MLS.rds") miDB_sig4.MLS <- readRDS("/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_NotaroTAD_A018/reference/miDB_sig4.MLS.rds") sub.obj<-subset(obs1, subset= NewlablesM == "KCs"| NewlablesM == "Macrophages"| NewlablesM == "Monocytes") sub.obj <- SetIdent(sub.obj, value = "Orig") obs1.dge <- FindMarkers(sub.obj, ident.1 = "Liver.TA33.Combo", ident.2 = "Liver.TA33", min.pct=0.01, logfc.threshold=0, min.cells.feature=1, min.cells.group=1) obs2.dge <- FindMarkers(sub.obj, ident.1 = "Tumor.TA33.Combo", ident.2 = "Tumor.TA33", min.pct=0.01, logfc.threshold=0, min.cells.feature=1, min.cells.group=1) obs1.ranks <- obs1.dge$avg_log2FC obs2.ranks <- obs2.dge$avg_log2FC names(obs2.ranks) <- rownames(obs2.dge) names(obs1.ranks) <- rownames(obs1.dge) res.combo.liver<-fgsea(pathways = miDB_sig4.MLS, stats = obs1.ranks, minSize=0, maxSize=500, eps=0) # nPermSimple=50000 res.combo.tumor<-fgsea(pathways = miDB_sig4.MLS, stats = obs2.ranks, minSize=0, maxSize=500, eps=0) # PermSimple=50000 # Select columns of interest res.combo.liver2<-res.combo.liver[,c(1,3,6,8)] res.combo.tumor2<-res.combo.tumor[,c(1,3,6,8)] # Merge in a single file and save xlsx merge1 <- merge(x=res.combo.liver2, y=res.combo.tumor2, by="pathway", suffixes = c("-liver", "-tumor")) write.xlsx(merge1, paste0(pdir, "/GSEA_LandT_KC-mono-macs.xlsx"), asTable=T, rowNames=T) ##################GSEA Heatmap of selected####################################### #rearrange data to plot heatmap res.combo.liver$origin<-"Liver" res.combo.tumor$origin<-"Tumor" # res.OVA$origin<- factor(x = res.OVA$origin, levels = c("OVA.Combo","OVA.Il12","OVA.IFNa","liOVA")) merge.results<-rbind(res.combo.liver,res.combo.tumor) #saveRDS(merge.results, "GSEA.liverandtumormacs.nperm50000.rds") #create variable with star for statistic merge.results$stats <- ifelse(merge.results$padj < 0.0005, "***", ifelse(merge.results$padj < 0.005, "**", ifelse(merge.results$padj < 0.05, "*", ""))) # Manually selected pathways a1<-c("IFNa_RO", "HALLMARK_INTERFERON_GAMMA_RESPONSE", "GOBP_DEFENSE_RESPONSE_TO_VIRUS", "GOBP_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION", "GOBP_POSITIVE_REGULATION_OF_CELL_KILLING", "GOCC_MHC_PROTEIN_COMPLEX", "GOBP_POSITIVE_REGULATION_OF_T_CELL_MEDIATED_CYTOTOXICITY", "GOCC_MHC_CLASS_II_PROTEIN_COMPLEX", "GOBP_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_ENDOGENOUS_ANTIGEN", "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", "GOBP_POSITIVE_REGULATION_OF_EPITHELIAL_CELL_PROLIFERATION", "GOBP_POSITIVE_REGULATION_OF_ENDOTHELIAL_CELL_PROLIFERATION", "PGE2_RO", "HALLMARK_ANGIOGENESIS", "IL4_RO" ) # Subset selected pathways merge.results1 <-merge.results[merge.results$pathway %in% a1,] # Rename pathways merge.results1[merge.results1 == "IFNa_RO"] <- "IFNa response (Cilenti et all)" merge.results1[merge.results1 == "HALLMARK_INTERFERON_GAMMA_RESPONSE"] <- "Hallmark_IFNg response" merge.results1[merge.results1 == "GOBP_DEFENSE_RESPONSE_TO_VIRUS"] <- "Antiviral response" merge.results1[merge.results1 == "GOBP_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION"] <- "Positive regulation of cell adhesion" merge.results1[merge.results1 == "GOBP_POSITIVE_REGULATION_OF_CELL_KILLING"] <- "Positive regulation of cell killing" merge.results1[merge.results1 == "GOCC_MHC_PROTEIN_COMPLEX"] <- "MHC protein complex" merge.results1[merge.results1 == "GOBP_POSITIVE_REGULATION_OF_T_CELL_MEDIATED_CYTOTOXICITY"] <- "Positive regulation of T cell cytotoxicity" merge.results1[merge.results1 == "GOCC_MHC_CLASS_II_PROTEIN_COMPLEX"] <- "MHC-II protein complex" merge.results1[merge.results1 == "GOBP_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_ENDOGENOUS_ANTIGEN"] <- "Antigen processing and presentation" merge.results1[merge.results1 == "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION"] <- "Hallmark EMT" merge.results1[merge.results1 == "GOBP_POSITIVE_REGULATION_OF_EPITHELIAL_CELL_PROLIFERATION"] <- "Positive regulation of epithelial cell proliferation" merge.results1[merge.results1 == "GOBP_POSITIVE_REGULATION_OF_ENDOTHELIAL_CELL_PROLIFERATION"] <- "Positive regulation of endothelial cell proliferation" merge.results1[merge.results1 == "PGE2_RO"] <- "PGE2 signature from Cilenti et al" merge.results1[merge.results1 == "HALLMARK_ANGIOGENESIS"] <- "Angiogenesis" merge.results1[merge.results1 == "IL4_RO"] <- "IL-4 signature from Cilenti et al" # Order in the way we want to see merge.results1$pathway<- factor(merge.results1$pathway, levels = c("IFNa response (Cilenti et all)", "Hallmark_IFNg response", "Antiviral response", "Positive regulation of cell adhesion", "Positive regulation of cell killing", "Positive regulation of T cell cytotoxicity", "Antigen processing and presentation", "MHC protein complex", "MHC-II protein complex", "Hallmark EMT", "Positive regulation of epithelial cell proliferation", "Positive regulation of endothelial cell proliferation", "Angiogenesis", "IL-4 signature from Cilenti et al", "PGE2 signature from Cilenti et al")) # Heatmap y<-ggplot(merge.results1, aes (x=origin, y=pathway, fill=NES)) + geom_tile() + scale_fill_gradientn(limits = c(-4,5.2), colours=c("Blue","White","Red"))+ geom_text(aes(label = stats), vjust=0.7,angle=0) + labs(y='',x='')+ theme_classic()+ theme(axis.text.y = element_text(color = "Black"), axis.text.x = element_text(angle=90, vjust=0.5, hjust=1, color="Black"), axis.line.x=element_line(color="Black"), axis.line.y=element_line(color="Black"))+ ggtitle("GSEA_MacroKCandmono vs OVA_Liver") ggsave(filename = paste0(pdir, "/Liver_and_Tumor_MacroKCsMonovsOVA_FigS6B.pdf"), plot=y, width=5, height=3) # Unbias heatmap merge.results <- merge.results[order(merge.results$padj, decreasing = F)] merge.results_temp <-merge.results[merge.results$padj < 0.0001,] filter_in <- head(unique(merge.results_temp$pathway), 100) merge.results2 <-merge.results[merge.results$pathway %in% filter_in,] merge.results2$origin<- factor(x = merge.results2$origin, levels = c("Liver","Tumor")) y<-ggplot(merge.results2, aes (x=origin, y=pathway, fill=NES)) + geom_tile() + scale_fill_distiller(palette = "RdBu") + geom_text(aes(label = stats), vjust = 0.7,angle = 0) + labs(y='',x='')+ theme_classic()+ theme(axis.text.y = element_text(color = "Black"), text = element_text(size = 15), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, color = "Black"), axis.line.x=element_line(color="Black"), axis.line.y=element_line(color="Black")) ggtitle("GSEA_MacroKCandmono vs OVA_Liver") ggsave(filename = paste0(pdir, "/Liver_and_Tumor_MacroKCsMonovsOVA_unbias.pdf"), plot=y, width=15, height=30) # limitsize = FALSE ############################################################################## ##################Figure S6D - RADAR PLOT###################################### ############################################################################## plotRadar <- function(i){ ploTData<- rbind(max1,min1,get(i)) id <- sub(".db","",i) radarchart(ploTData, title = id, pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty=1, cglwd=0.8) } #function to plot radarplot ####Create label to cluster dc and granulocytes together obs1@meta.data$Shortlables<-obs1@meta.data$NewlablesM #DCs<-c("CD8 cDC1","moDCs", "Ccr7 DCs","pDCs","pre DCs") #obs1@meta.data$Shortlables[obs1@meta.data$NewlablesM%in%DCs] <- "DCs" #granu<-c("Neutrophils","Basophils") #obs1@meta.data$Shortlables[obs1@meta.data$NewlablesM%in%granu] <- "Granulocytes" obs1$Shortlables <- paste0(obs1$Shortlables,".",obs1$Tissue) # AddModuleScore MHCI_genes<- list(c("H2-Q7", "H2-T22", "H2-Q4", "H2-K1", "H2-T23", "H2-M3", "H2-D1","Tap1","Tap2")) MHCII_genes<- list( c("H2-Eb1", "H2-Ab1", "H2-Aa", "H2-DMb1", "H2-DMa", "Ciita")) Protumoral_genes<- list(c("Vav2","Tgfb1","Il10","Ccl24","Mmp8","Tmem176b","Trem2","Fn1") ) IFNa_genes<- list(c("Ifi44", "Oasl1","Irf7","Isg15", "Isg20","Oas1a","Oas1g") ) IFNg_genes<-list(c("Ccl5", "Upp1", "Slamf7", "Cxcl9", "Gbp4","Cd74")) obs1<- AddModuleScore(obs1, features = MHCI_genes, name = "MHCI_genes") obs1<- AddModuleScore(obs1, features = MHCII_genes, name = "MHCII_genes") obs1<- AddModuleScore(obs1, features = Protumoral_genes, name = "Protumoral_genes") obs1<- AddModuleScore(obs1, features = IFNa_genes, name = "IFNa_genes") obs1<- AddModuleScore(obs1, features = IFNg_genes, name = "IFNg_genes") # Normalize data #we needed to find a way to normalize the module scores. #we used a quantile normalization method. #M1 <- obs1@meta.data[,c(37,42, 43:47)] M1 <- obs1@meta.data[,c("Shortlables", "RNA_Group", "MHCI_genes1","MHCII_genes1","Protumoral_genes1","IFNa_genes1","IFNg_genes1")] M1[,3:7]<-limma::normalizeQuantiles(M1[,3:7]) Groups <-unique(obs1$RNA_Group) CellType<-unique(M1$Shortlables) #with Median #TAMs #Groups for (i in CellType){ M2<-M1[M1$Shortlables%in%i,] 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) M4<-rbind(M3a,M3b) rownames(M4)<-Groups colnames(M4)<-colnames(M2[,3:7]) x<-as.data.frame(M4) assign(paste0(i,".db"),x)} lsCells <- grep(".db",ls(),value=T) MinMax<-lapply(lsCells,get) max1 <- rep(do.call("max", MinMax),5) min1 <- rep(do.call("min", MinMax),5) colors_border=c(rgb(0,0.8,0,0.9), rgb(0,0,0,0.9)) colors_in=c(rgb(0,0.3,0,0.), rgb(0,0,0.0,0.1)) pdf(paste0(pdir, "/Radar_plot_liver_and_tumor_APCs_FigS6D_new2.pdf"), width=15, height=12) par(mfrow=c(4,5)) for (i in 1:length(unique(obs1$Shortlables))){ plotRadar(lsCells[i]) } dev.off() #######################Plot UMAP################ liver.apcs <- readRDS("/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_NotaroTAD_A019/Analysis_MM/results3/Liver_APCs/Liver_APCs_final.rds") tumor.apcs <- readRDS("/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_NotaroTAD_A019/Analysis_MM/results3/Tumor_APCs/Tumor_APCs_final.rds") liver.apcs@meta.data$Tissue <- liver.apcs@meta.data$PROP.Source liver.apcs@meta.data$RNA_Group <- liver.apcs@meta.data$PROP.Condition liver.apcs@meta.data$RNA_Group[liver.apcs@meta.data$RNA_Group == "TA_only"] <- "TA33" liver.apcs@meta.data$RNA_Group[liver.apcs@meta.data$RNA_Group == "TA_combo"] <- "TA33.Combo" liver.apcs@meta.data$RNA_Group <- factor(x = liver.apcs@meta.data$RNA_Group, levels = c("TA33","TA33.Combo")) liver.apcs$Tube <- paste0(liver.apcs$hash.ID,liver.apcs$RNA_Group) tumor.apcs@meta.data$Tissue <- tumor.apcs@meta.data$PROP.Source tumor.apcs@meta.data$RNA_Group <- tumor.apcs@meta.data$PROP.Condition tumor.apcs@meta.data$RNA_Group[tumor.apcs@meta.data$RNA_Group == "TA_only"] <- "TA33" tumor.apcs@meta.data$RNA_Group[tumor.apcs@meta.data$RNA_Group == "TA_combo"] <- "TA33.Combo" tumor.apcs@meta.data$RNA_Group <- factor(x = tumor.apcs@meta.data$RNA_Group, levels = c("TA33","TA33.Combo")) tumor.apcs$Tube <- paste0(tumor.apcs$hash.ID,tumor.apcs$RNA_Group) # Save plot pdf(paste0(pdir, "/liver_apcs_umap_FigS11A-1.pdf"), width=7, height=6) # Figure S11A-1 DimPlot(liver.apcs, reduction = "umap.harmony.orig.ident",group.by = "NewlablesM_clu", pt.size = 0.5, label = T, label.size = 3, repel = T) dev.off() pdf(paste0(pdir, "/tumor_apcs_umap_FigS11A-2.pdf"), width=7, height=6) # Figure S11A-2 DimPlot(tumor.apcs, reduction = "umap.harmony.orig.ident",group.by = "NewlablesM_clu", pt.size = 0.5, label = T, label.size = 3, repel = T) dev.off() # Export DimPlot as excel var_list <- c("UMAPh_1", "UMAPh_2","Tissue","orig.ident","Sample","RNA_Group","Tube","RNA_snn_h.orig.ident_res.1.2","NewlablesM_clu") liver.apcs_df <- FetchData(liver.apcs, vars = var_list) liver.apcs_df$cell_ID <- rownames(liver.apcs_df) write.xlsx(liver.apcs_df, file=paste0(pdir, "/liver_apcs_umap_FigS11A-1.xlsx"), overwrite=T) tumor.apcs_df <- FetchData(tumor.apcs, vars = var_list) tumor.apcs_df$cell_ID <- rownames(tumor.apcs_df) write.xlsx(tumor.apcs_df, file=paste0(pdir, "/tumor_apcs_umap_FigS11A-2.xlsx"), overwrite=T) #################################ENDS here##############################