library(Seurat) library(ggplot2) library(ComplexHeatmap) library(dplyr) library(circlize) library(RColorBrewer) library(gridExtra) library(cowplot) library(stringr) library(scales) library(clusterProfiler) library(ggrepel) library(openxlsx) # Functions plotVolcano <- function(dge_res, adj.pval.thr, logfc.thr, hgenes = NULL) { data <- data.frame(gene = row.names(dge_res), pval = -log10(dge_res$p_val_adj), lfc = dge_res$avg_logFC) data <- na.omit(data) data <- mutate(data, color = case_when(data$lfc > logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Overexpressed", data$lfc < -logfc.thr & data$pval > -log10(adj.pval.thr) ~ "Underexpressed", abs(data$lfc) < logfc.thr | data$pval < -log10(adj.pval.thr) ~ "NonSignificant")) data <- data[order(data$pval, decreasing = TRUE),] if (is.null(hgenes)) { highl <- head(subset(data, color != "NonSignificant"), 20) } else { highl <- subset(data, gene %in% hgenes) } vol <- ggplot(data, aes(x = lfc, y = pval, colour = color)) + theme_bw(base_size = 20) + theme(legend.position = "none", panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + xlab("log2 Fold Change") + ylab(expression(-log[10]("adjusted p-value"))) + geom_hline(yintercept = -log10(adj.pval.thr), colour = "darkgray") + geom_point(size = 2, alpha = 0.8, na.rm = T) + scale_color_manual(name = "Expression", values = c(Overexpressed = "#CD4F39", Underexpressed = "#008B00", NonSignificant = "darkgray")) + geom_label_repel(data = subset(highl, lfc > 0), aes(label = gene), size = 8, show.legend = FALSE, direction = "y", nudge_x = 1, nudge_y = 40, hjust = 0) + geom_label_repel(data = subset(highl, lfc <= 0), aes(label = gene), size = 8, show.legend = FALSE, direction = "y", nudge_x = -.5) return(vol) } doSCHeatmapMean <- function(fname, clustering, ntop, num_color, high_genes, clu_lab, duplicated_labs = F) { col_scales <- list(colorRampPalette(rev(brewer.pal(n = 9, name ="RdBu")))(100), colorRamp2(c(-4, -3, -2, -1, 0, 1, 2, 3, 4), rev(brewer.pal(n = 9, name ="RdBu"))), colorRamp2(c(-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2), rev(brewer.pal(n = 9, name ="RdBu"))), colorRamp2(c(-1.5, -1, -0.5, 0, 0.5, 1, 1.5), rev(brewer.pal(n = 7, name ="RdBu"))), colorRamp2(c(-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1), rev(brewer.pal(n = 9, name ="RdBu")))) if (num_color > length(col_scales)) { col_sc <- col_scales[[length(col_scales)]] } else { col_sc <- col_scales[[num_color]] } obj <- readRDS(fname) top_genes <- obj@misc$markers[[clustering]] %>% group_by(cluster) %>% top_n(ntop, avg_logFC) if (nrow(obj@assays$RNA@scale.data) == 0) { hm_data <- as.matrix(obj@assays$RNA@data[top_genes$gene,]) #print(hm_data[1:3, 1:5]) hm_data <- t(scale(t(hm_data))) #print(hm_data[1:3, 1:5]) } else { hm_data <- obj@assays$RNA@scale.data[top_genes$gene,] } mm <- merge(obj@meta.data, t(hm_data), by = 0) rownames(mm) <- mm$Row.names mm$Row.names <- NULL mm_filt <- mm[, c(clustering, top_genes$gene)] mm_matrix <- mm_filt %>% group_by_at(clustering) %>% summarise_all(list(mean)) %>% data.frame() #print(mm_matrix) c19 <- hue_pal()(nrow(mm_matrix)) if (length(clu_lab) > 0) { rownames(mm_matrix) <- clu_lab names(c19) <- clu_lab } else { rownames(mm_matrix) <- mm_matrix[[clustering]] names(c19) <- mm_matrix[[clustering]] } mm_matrix[[clustering]] <- NULL hm_matrix <- as.matrix(mm_matrix) if (duplicated_labs) { colnames(hm_matrix) <- str_split_fixed(colnames(hm_matrix), pattern = "[.]", n = 2)[,1] # !!! } gg <- intersect(colnames(hm_matrix), high_genes) gg_index <- which(colnames(hm_matrix) %in% gg) labels <- colnames(hm_matrix)[gg_index] ha_tre <- columnAnnotation(link = anno_mark(at = gg_index, labels = labels, labels_rot = 45, padding = 0.75, labels_gp = gpar(fontsize = 14), link_width = unit(15, "mm"), side = "bottom")) if (length(clu_lab) > 0) { rsplit <- factor(clu_lab, levels = clu_lab) } else { rsplit <- as.numeric(levels(obj@meta.data[[clustering]])) } ha_col <- rowAnnotation(df = data.frame(Cluster = rsplit), show_annotation_name = FALSE, show_legend = FALSE, col = list(Cluster = c19), annotation_legend_param = list(Cluster = list(grid_height = unit(.5, "cm"), grid_width = unit(.5, "cm"), title_gp = gpar(col = "black", fontsize = 14), labels_gp = gpar(col = "black", fontsize = 12)))) hm = Heatmap(matrix = hm_matrix, name = "Expr.", col = col_sc, cluster_columns = FALSE, cluster_rows = FALSE, rect_gp = gpar(col = "grey", lwd = 0), border = FALSE, show_row_names = FALSE, show_column_names = FALSE, column_title_gp = gpar(fontsize = 0), row_split = rsplit, #cluster_row_slices = TRUE, #column_split = obj@meta.data[[clustering]], column_split = rep(0:(nrow(hm_matrix) - 1), each = ntop), column_gap = unit(4, "mm"), row_title_gp = gpar(fontsize = 18), row_title_rot = 0, left_annotation = ha_col, bottom_annotation = ha_tre) return(hm) } ggtable <- function(d, p = NULL) { as.ggplot(tableGrob2(d, p)) } tableGrob2 <- function(d, p = NULL, labels = NULL) { d <- d[order(rownames(d)),] if (is.null(labels)) { tp <- tableGrob(d, theme=ttheme_default(base_size = 24)) } else { tp <- tableGrob(d, rows = labels, parse = TRUE, theme=ttheme_default(base_size = 24)) } if (is.null(p)) { return(tp) } pcol <- unique(ggplot_build(p)$data[[1]][["colour"]]) j <- which(tp$layout$name == "rowhead-fg") for (i in seq_along(pcol)) { tp$grobs[j][[i+1]][["gp"]] = gpar(col = pcol[i], fontsize = 24) } return(tp) } gseaScores <- getFromNamespace("gseaScores", "DOSE") gsInfo <- function(object, geneSetID) { geneList <- object@geneList if (is.numeric(geneSetID)) geneSetID <- object@result[geneSetID, "ID"] geneSet <- object@geneSets[[geneSetID]] exponent <- object@params[["exponent"]] df <- gseaScores(geneList, geneSet, exponent, fortify=TRUE) df$ymin=0 df$ymax=0 pos <- df$position == 1 h <- diff(range(df$runningScore))/20 df$ymin[pos] <- -h df$ymax[pos] <- h df$geneList <- geneList df$Description <- object@result[geneSetID, "Description"] return(df) } gseaplot2 <- function(x, geneSetID, title = "", color="green", base_size = 14, rel_heights=c(1.5, .5, 1), subplots = 1:3, pvalue_table = FALSE, ES_geom="line", nopadj = FALSE, labels = NULL) { ES_geom <- match.arg(ES_geom, c("line", "dot")) geneList <- position <- NULL ## to satisfy codetool if (length(geneSetID) == 1) { gsdata <- gsInfo(x, geneSetID) } else { gsdata <- do.call(rbind, lapply(geneSetID, gsInfo, object = x)) } p <- ggplot(gsdata, aes_(x = ~x)) + xlab(NULL) + #theme_classic(base_size) + theme_classic(base_size) + theme(panel.grid.major = element_line(colour = "grey92"), panel.grid.minor = element_line(colour = "grey92"), panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank()) + scale_x_continuous(expand=c(0,0)) if (ES_geom == "line") { es_layer <- geom_line(aes_(y = ~runningScore, color= ~Description), size=1) } else { es_layer <- geom_point(aes_(y = ~runningScore, color= ~Description), size=1, data = subset(gsdata, position == 1)) } p.res <- p + es_layer + theme(legend.position = c(.8, .8), legend.title = element_blank(), legend.background = element_rect(fill = "transparent")) p.res <- p.res + ylab("Running Enrichment Score") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.line.x=element_blank(), plot.margin=margin(t=.2, r = .2, b=0, l=.2, unit="cm")) i <- 0 for (term in unique(gsdata$Description)) { idx <- which(gsdata$ymin != 0 & gsdata$Description == term) gsdata[idx, "ymin"] <- i gsdata[idx, "ymax"] <- i + 1 i <- i + 1 } p2 <- ggplot(gsdata, aes_(x = ~x)) + geom_linerange(aes_(ymin=~ymin, ymax=~ymax, color=~Description)) + xlab(NULL) + ylab(NULL) + theme_classic(base_size) + theme(legend.position = "none", plot.margin = margin(t=-.1, b=0,unit="cm"), axis.ticks = element_blank(), axis.text = element_blank(), axis.line.x = element_blank()) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0)) if (length(geneSetID) == 1) { ## geneList <- gsdata$geneList ## j <- which.min(abs(geneList)) ## v1 <- quantile(geneList[1:j], seq(0,1, length.out=6))[1:5] ## v2 <- quantile(geneList[j:length(geneList)], seq(0,1, length.out=6))[1:5] ## v <- sort(c(v1, v2)) ## inv <- findInterval(geneList, v) v <- seq(1, sum(gsdata$position), length.out=9) inv <- findInterval(rev(cumsum(gsdata$position)), v) if (min(inv) == 0) inv <- inv + 1 col=c(rev(brewer.pal(5, "Blues")), brewer.pal(5, "Reds")) ymin <- min(p2$data$ymin) yy <- max(p2$data$ymax - p2$data$ymin) * .3 xmin <- which(!duplicated(inv)) xmax <- xmin + as.numeric(table(inv)[unique(inv)]) d <- data.frame(ymin = ymin, ymax = yy, xmin = xmin, xmax = xmax, col = col[unique(inv)]) p2 <- p2 + geom_rect( aes_(xmin=~xmin, xmax=~xmax, ymin=~ymin, ymax=~ymax, fill=~I(col)), data=d, alpha=.9, inherit.aes=FALSE) } ## p2 <- p2 + ## geom_rect(aes(xmin=x-.5, xmax=x+.5, fill=geneList), ## ymin=ymin, ymax = ymin + yy, alpha=.5) + ## theme(legend.position="none") + ## scale_fill_gradientn(colors=color_palette(c("blue", "red"))) df2 <- p$data #data.frame(x = which(p$data$position == 1)) df2$y <- p$data$geneList[df2$x] p.pos <- p + geom_segment(data=df2, aes_(x=~x, xend=~x, y=~y, yend=0), color="grey") p.pos <- p.pos + ylab("Ranked list metric") + xlab("Rank in Ordered Dataset") + theme(plot.margin=margin(t = -.1, r = .2, b=.2, l=.2, unit="cm")) #if (!is.null(title) && !is.na(title) && title != "") { p.res <- p.res + ggtitle(title) + theme(plot.title = element_text(size = 40, face = "bold", hjust = 0.5)) #} if (length(color) == length(geneSetID)) { p.res <- p.res + scale_color_manual(values=color) if (length(color) == 1) { p.res <- p.res + theme(legend.position = "none") p2 <- p2 + scale_color_manual(values = "black") } else { p2 <- p2 + scale_color_manual(values = color) } } if (pvalue_table) { if (nopadj) { pd <- x[geneSetID, c("Description", "pvalue", "NES")] } else { pd <- x[geneSetID, c("Description", "pvalue", "p.adjust", "NES")] } pd <- pd[order(pd[,1], decreasing=FALSE),] #pd$Description <- gsub("HALLMARK_", "", pd$Description) #rownames(pd) <- pd$Description pd <- pd[,-1] pd <- round(pd, 4) #tp <- tableGrob2(pd, p.res) if (is.null(labels)) { tp <- tableGrob(pd, rows = NULL, theme=ttheme_default(base_size = 26)) p.res <- p.res + theme(legend.position = "none") + annotation_custom(tp, xmin = quantile(p.res$data$x, .65), xmax = quantile(p.res$data$x, .95), ymin = quantile(p.res$data$runningScore, .5), ymax = quantile(p.res$data$runningScore, .9)) } else { tp <- tableGrob2(pd, p.res, labels) p.res <- p.res + theme(legend.position = "none") + annotation_custom(tp, xmin = quantile(p.res$data$x, .6), xmax = quantile(p.res$data$x, .9), ymin = quantile(p.res$data$runningScore, .8), ymax = quantile(p.res$data$runningScore, .95)) } } plotlist <- list(p.res, p2, p.pos)[subplots] n <- length(plotlist) plotlist[[n]] <- plotlist[[n]] + theme(axis.line.x = element_line(), axis.ticks.x=element_line(), axis.text.x = element_text()) if (length(subplots) == 1) { return(plotlist[[1]] + theme(plot.margin=margin(t=.2, r = .2, b=.2, l=.2, unit="cm"))) } if (length(rel_heights) > length(subplots)) { rel_heights <- rel_heights[subplots] } plot_grid(plotlist = plotlist, ncol=1, align="v", rel_heights=rel_heights) } ################# wdir <- "Birocchi_SciTraMed2022/results" out_dir <- "Birocchi_SciTraMed2022/Paper" ################## #### Figure 1 #### ################## # 1 - UMAP cluster res 0.6 fobj <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/")) full_umap_plot <- cbind(fobj@reductions$umap@cell.embeddings, fobj@meta.data[rownames(fobj@reductions$umap@cell.embeddings),]) label <- data.frame( UMAP_1 = c(-9.2, -7.5, -8, -6, -5.4, -4.8, -4, -2.5, -2.5, -0.9, -2, 1, 1.5, 4, 4, 7, 4, 5.2, 9), UMAP_2 = c(1, -1, 6, 1, -1.5, -3.3, 6.5, -1.7, -5.5, -3.3, -8, -6.5, -10, -3.5, -1, 2.5, 3, 6.5, 3), RNA_snn_res.0.6 = c(9, 3, 10, 4, 14, 5, 12, 18, 6, 17, 13, 8, 7, 11, 2, 1, 0, 15, 16), label = c("Naive~T~cells", "CD4", "NK", "CD8", "T~regs", "Prolif.~T~cells", "B~cells", "Mast~cells", "Unknown", "pDC", "DC1", "DC2", "mregDC", "Monocytes", "Macro~3", "Macro~2", "Macro~1", "Resid.like~Macro", "Prolif.~Macro") ) p <- ggplot(full_umap_plot, aes(UMAP_1, UMAP_2)) + geom_point(aes(colour = RNA_snn_res.0.6), size = 1, alpha = 1) + theme_nothing(font_size = 20) + #theme_bw(base_size = 20) + theme(legend.position = "none") + #geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) + geom_label(data = label, aes(label = label), parse = T, size = 5) + #guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) + xlab("") + ylab("") ggsave(filename = paste(out_dir, "Fig1_1-Full_fastMNN_Res06.png", sep = "/"), plot = p, width = 10, height = 10, dpi = 100) # 2 - Dotplot of enrichment pos HALLMARK ds_high <- readRDS(paste(wdir, "Full", "02-fastMNN_2-post-analysis", "CD_High-AB_Ctrl", "CD_High-AB_Ctrl-post-analysis.rds", sep = "/")) ds_high_hall_enrich <- ds_high$results$enrichment$enrichment_h.all.v6.1.symbols_mouse_pos@result[,c("ID", "GeneRatio", "p.adjust", "Count")] ds_high_hall_enrich$Contr <- "IFN-alpha" ds_low <- readRDS(paste(wdir, "Full", "02-fastMNN_2-post-analysis", "EF_Low-AB_Ctrl", "EF_Low-AB_Ctrl-post-analysis.rds", sep = "/")) ds_low_hall_enrich <- ds_low$results$enrichment$enrichment_h.all.v6.1.symbols_mouse_pos@result[,c("ID", "GeneRatio", "p.adjust", "Count")] ds_low_hall_enrich$Contr <- "IFN-alpha-DHFR" ds_hall_enrich <- rbind(ds_high_hall_enrich, ds_low_hall_enrich) gsize <- as.numeric(sub("/\\d+$", "", as.character(ds_hall_enrich$GeneRatio))) gcsize <- as.numeric(sub("^\\d+/", "", as.character(ds_hall_enrich$GeneRatio))) ds_hall_enrich$GeneRatio = gsize/gcsize ds_hall_enrich <- subset(ds_hall_enrich, p.adjust < 0.05) ds_hall_enrich$ID <- gsub("HALLMARK_", "", ds_hall_enrich$ID) #ds_hall_enrich$ID <- tolower(gsub("_", " ", ds_hall_enrich$ID)) print(head(ds_hall_enrich)) id_order <- ds_hall_enrich %>% select(ID, GeneRatio) %>% group_by(ID) %>% summarize_all(sum) %>% data.frame() id_order <- id_order[order(id_order$GeneRatio, decreasing = F),] ds_hall_enrich$ID <- factor(ds_hall_enrich$ID, levels = id_order$ID) p <- ggplot(ds_hall_enrich, aes(x = GeneRatio, y = ID, size = Count, color = p.adjust)) + theme_bw(base_size = 20) + scale_color_gradient(low = "red", high = "blue", guide = guide_colorbar(reverse = T)) + geom_point() + ylab("") + facet_grid(.~Contr, labeller = label_parsed) ggsave(plot = p, filename = paste(out_dir, "Fig1_2-Full_Hallmark_enrichment_pos.png", sep = "/"), width = 12, height = 8, dpi = 100) # 3/4 - GSEA hallmark ifna response High vs CTRL o Low vs CTRL ds_high <- readRDS(paste(wdir, "Full", "02-fastMNN_2-post-analysis", "CD_High-AB_Ctrl", "CD_High-AB_Ctrl-post-analysis.rds", sep = "/")) ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse@result[1,] p <- gseaplot2(ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse, geneSetID = 1, base_size = 18, title = expression(IFN-alpha), color = c("red"), pvalue_table = T) ggsave(filename = paste(out_dir, "Fig1_3-CD_High-AB_Ctrl_HALLMARK_INT_ALPHA_GSEA.png", sep = "/"), plot = p, width = 12, height = 8, dpi = 100) ds_low <- readRDS(paste(wdir, "Full", "02-fastMNN_2-post-analysis", "EF_Low-AB_Ctrl", "EF_Low-AB_Ctrl-post-analysis.rds", sep = "/")) p <- gseaplot2(ds_low$results$gsea$gsea_h.all.v6.1.symbols_mouse, geneSetID = 1, base_size = 18, title = expression(IFN-alpha-DHFR), color = c("red"), pvalue_table = T) ggsave(filename = paste(out_dir, "Fig1_4-EF_Low-AB_Ctrl_HALLMARK_INT_ALPHA_GSEA.png", sep = "/"), plot = p, width = 12, height = 8, dpi = 100) # 5 - Heatmap DEGs HIGH or LOW vs CTRL for each cluster fobj <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/")) clu_deg <- fobj@misc[["deg_clumarkers"]] full_res_genes <- c() for (clu in names(clu_deg)) { for (contr in c("EF_Low-AB_Ctrl", "CD_High-AB_Ctrl")) { clu_name <- paste0(clu, "-", contr) res <- subset(clu_deg[[clu]][[contr]], p_val_adj < 0.05 & abs(avg_logFC) >= 1) full_res_genes <- union(full_res_genes, rownames(res)) } } full_res <- data.frame() for (clu in as.character(seq(0, 18, 1))) { for (contr in c("EF_Low-AB_Ctrl", "CD_High-AB_Ctrl")) { clu_name <- paste0(clu, "-", contr) res <- clu_deg[[clu]][[contr]] res_filt <- res[rownames(res) %in% full_res_genes, "avg_logFC", drop = FALSE] colnames(res_filt) <- clu_name full_res <- merge(full_res, res_filt, by = 0, all = T) rownames(full_res) <- full_res$Row.names full_res$Row.names <- NULL } } nrow(full_res) full_res[is.na(full_res)] <- 0 sinfo <- data.frame(Dataset = colnames(full_res)) sinfo <- cbind(sinfo, str_split_fixed(sinfo$Dataset, "-", 3)) colnames(sinfo) <- c("Dataset", "Cluster", "Dose", "Ctrl") rownames(sinfo) <- sinfo$Dataset sinfo$Cluster <- factor(sinfo$Cluster, levels = as.character(seq(0, 18, 1))) sinfo$Ctrl <- NULL sinfo$Dataset <- NULL sinfo$Treatment <- ifelse(sinfo$Dose == "EF_Low", "IFN-a-DHFR", "IFN-a") sinfo$Dose <- NULL ctypes <- c("Macro1","Macro2", "Macro3","CD4","CD8","ProlifT","Unk", "mregDC","DC2","NaiveT","NK","Mono","Bcells","DC1","Tregs", "ResMacro","ProlifMacro","pDC","MastCells") levels(sinfo$Cluster) <- ctypes inf_alpha <- c("1700057G04Rik","AC125149.4","AC125149.5","AC132444.3","AC133103.4","AC133103.5","AC168977.1","Adar","B2m","Batf2","Bst2","C1s1","C1s2","Casp1", "Casp8","Ccrl2","Cd47","Cd74","Cmpk2","Cmtr1","Cnp","Csf1","Cxcl10","Cxcl11","Ddx60","Dhx58","Eif2ak2","Elf1","Epsti1","Gbp2","Gbp2b","Gmpr", "Helz2","Herc6","Ifi35","Ifi44","Ifi44l","Ifih1","Ifit2","Ifit3","Ifit3b","Ifitm1","Ifitm2","Ifitm3","Ifitm6","Ifitm7","Il15","Il4ra","Il7", "Irf1","Irf2","Irf7","Irf9","Isg15","Isg20","Lap3","Lgals3bp","Lpar6","Ly6e","Mov10","Mvb12a","Mx1","Mx2","Ncoa7","Nmi","Nub1","Oas1a","Oas1b", "Oas1c","Oas1d","Oas1e","Oas1f","Oas1g","Oas1h","Oasl1","Ogfr","Parp12","Parp14","Parp9","Plscr1","Plscr2","Pnpt1","Procr","Psma3","Psmb8","Psmb9", "Psme1","Psme2","Psme2b","Ripk2","Rnf31","Rsad2","Rtp4","Samd9l","Sell","Slc25a28","Sp110","Stat2","Tap1","Tdrd7","Tent5a","Tmem140","Trafd1", "Trim12a","Trim12c","Trim14","Trim21","Trim25","Trim26","Trim30a","Trim30c","Trim30d","Trim5","Txnip","Uba7","Ube2l6","Usp18","Wars") ifn_font <- ifelse(rownames(full_res) %in% inf_alpha, 2, 1) ifn_col <- ifelse(rownames(full_res) %in% inf_alpha, "red", "black") c19 <- hue_pal()(19) names(c19) <- levels(sinfo$Cluster) ha_col <- columnAnnotation(df = sinfo, show_annotation_name = FALSE, col = list(Treatment = c("IFN-a-DHFR" = "#6699FF", "IFN-a" = "#CC0066"), Cluster = c19), annotation_name_side = "right", annotation_legend_param = list(Treatment = list(at = c("IFN-a-DHFR", "IFN-a"), grid_height = unit(.5, "cm"), grid_width = unit(.5, "cm"), title_gp = gpar(col = "black", fontsize = 14), labels = c(expression(paste("IFN-", alpha, "-DHFR ")), expression(paste("IFN-", alpha, " "))), # labels = gt_render(c("IFN-&alpha-DHFR ", "IFN-alpha- ")), labels_gp = gpar(col = "black", fontsize = 12)), Cluster = list(grid_height = unit(.5, "cm"), grid_width = unit(.5, "cm"), title_gp = gpar(col = "black", fontsize = 14), labels_gp = gpar(col = "black", fontsize = 12)))) hm = Heatmap(matrix = as.matrix(full_res), col = colorRampPalette(rev(brewer.pal(n = 9, name ="RdYlBu")))(100), cell_fun = function(j, i, x, y, width, height, fill) { if(abs(as.matrix(full_res)[i, j]) > 0.5) grid.text(sprintf("%.1f", as.matrix(full_res)[i, j]), x, y, gp = gpar(fontsize = 9)) }, cluster_columns = FALSE, cluster_rows = TRUE, rect_gp = gpar(col = "grey", lwd = 0), border = FALSE, name = "LogFC", show_row_names = TRUE, row_names_gp = gpar(fontsize = 11, col = ifn_col), show_column_names = F, column_title_gp = gpar(col = "black", fontsize = 11), heatmap_legend_param = list(legend_height = unit(3, "cm"), title_gp = gpar(col = "black", fontsize = 14), labels_gp = gpar(col = "black", fontsize = 12)), column_gap = unit(1, "mm"), column_split = sinfo$Cluster, top_annotation = ha_col, #split = 2, row_dend_reorder = TRUE, width = 18, height = 10) png(paste(out_dir, "Fig1_5-Full_fastMNN_deg_clumarkers_abslogFC1_fdr05_05Lab_noColNames.png", sep = "/"), width = 1800, height = 1000, res = 100) draw(hm, merge_legends = T) dev.off() ######################### #### Figure 1 Suppl. #### ######################### # 1 - Heatmap of top marker mean expression hm <- doSCHeatmapMean(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"), "RNA_snn_res.0.6", 10, 3, c("Saa3", "Chil3", "Ly6i", "C1qc", "Selenop", "Ms4a7", "Il1b", "Dusp1", "Socs3", "Furin", "Penk", "Trbc2", "Cd8b1", "Cd8a", "Klrd1", "Ube2c", "Dut", "Stmn1", "Fscn1", "Il2ra", "Pkib", "Ifi27l2a", "Tmem123", "Socs2", "Ccl17", "Il1r2", "Plet1", "S1pr1", "Klf2", "Dusp10", "Prf1", "Car2", "Klrb1c", "Cxcl10", "Ifitm6", "Clec4e", "Igkc", "Igha", "Cd79a", "Cd24a", "Xcr1", "Naaa", "Izumo1r", "Cd2", "Ctla4", "Hexb", "Cd81", "Psap", "Cks1b", "Top2a", "Siglech", "Ccr9", "Rnase6", "Mcpt2", "Cma1", "Cpa3"), c("Macro1","Macro2", "Macro3","CD4","CD8","Prolif. T","Unknown", "mregDC","DC2","Naive T","NK","Monocytes", "B cells","DC1","T regs","Res. Macro","Prolif. Macro","pDC","Mast Cells")) png(paste(out_dir, "Fig1Suppl_1-Full_fastMNN_HM_top10_c3_clumean_topnames_2.png", sep = "/"), width = 1500, height = 800, res = 100) draw(hm) dev.off() # 2 - HALLMARK_IFNA_RESPONSE modulescore expression fobj <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/")) inf_alpha <- c("1700057G04Rik","AC125149.4","AC125149.5","AC132444.3","AC133103.4","AC133103.5","AC168977.1","Adar","B2m","Batf2","Bst2","C1s1","C1s2","Casp1", "Casp8","Ccrl2","Cd47","Cd74","Cmpk2","Cmtr1","Cnp","Csf1","Cxcl10","Cxcl11","Ddx60","Dhx58","Eif2ak2","Elf1","Epsti1","Gbp2","Gbp2b","Gmpr", "Helz2","Herc6","Ifi35","Ifi44","Ifi44l","Ifih1","Ifit2","Ifit3","Ifit3b","Ifitm1","Ifitm2","Ifitm3","Ifitm6","Ifitm7","Il15","Il4ra","Il7", "Irf1","Irf2","Irf7","Irf9","Isg15","Isg20","Lap3","Lgals3bp","Lpar6","Ly6e","Mov10","Mvb12a","Mx1","Mx2","Ncoa7","Nmi","Nub1","Oas1a","Oas1b", "Oas1c","Oas1d","Oas1e","Oas1f","Oas1g","Oas1h","Oasl1","Ogfr","Parp12","Parp14","Parp9","Plscr1","Plscr2","Pnpt1","Procr","Psma3","Psmb8","Psmb9", "Psme1","Psme2","Psme2b","Ripk2","Rnf31","Rsad2","Rtp4","Samd9l","Sell","Slc25a28","Sp110","Stat2","Tap1","Tdrd7","Tent5a","Tmem140","Trafd1", "Trim12a","Trim12c","Trim14","Trim21","Trim25","Trim26","Trim30a","Trim30c","Trim30d","Trim5","Txnip","Uba7","Ube2l6","Usp18","Wars") inf_alpha_filt <- inf_alpha[inf_alpha %in% rownames(fobj)] fobj <- AddModuleScore(object = fobj, features = list(IFN_ALPHA = inf_alpha_filt), ctrl = length(inf_alpha_filt), name = "IFN_ALPHA") umap_plot <- cbind(fobj@reductions$umap@cell.embeddings, fobj@meta.data[rownames(fobj@reductions$umap@cell.embeddings),]) umap_plot$RNA_Group2 <- factor(umap_plot$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low")) levels(umap_plot$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR") str(umap_plot) summary(umap_plot$IFN_ALPHA1) p <- ggplot(umap_plot, aes(x = UMAP_1, y = UMAP_2)) + geom_point(data = subset(umap_plot, IFN_ALPHA1 > 0.2 & IFN_ALPHA1 < 0.5), aes(colour = IFN_ALPHA1), size = 1, alpha = 0.8) + geom_point(data = subset(umap_plot, IFN_ALPHA1 <= 0.2 & IFN_ALPHA1 > 0.1), aes(colour = IFN_ALPHA1), size = 1, alpha = 0.8) + geom_point(data = subset(umap_plot, IFN_ALPHA1 >= 0.5 & IFN_ALPHA1 < 0.75), aes(colour = IFN_ALPHA1), size = 1, alpha = 0.8) + geom_point(data = subset(umap_plot, IFN_ALPHA1 <= 0.1), aes(colour = IFN_ALPHA1), size = 1, alpha = 1) + geom_point(data = subset(umap_plot, IFN_ALPHA1 >= 0.75), aes(colour = IFN_ALPHA1), size = 1, alpha = 1) + scale_color_gradientn(colours = c("#2166AC","#4393C3","#D1E5F0","#F4A582","#B2182B"), name = "Expr.", #values = rescale(c(-.1,mean(macro_umap_plot$IFNa1),.9)), values = rescale(c(-.1,.05,.2,.4,.7)), guide = "colorbar") + #theme_bw(base_size = 20) + #theme_nothing(font_size = 20) + theme_void(base_size = 20) + #ggtitle("INTERFERON ALPHA RESPONSE") + theme(legend.position = "right", strip.text.x = element_text(size = 40)) + #xlab("UMAP 1") + #ylab("UMAP 2") + facet_grid(.~RNA_Group2, labeller = label_parsed) ggsave(filename = paste(out_dir, "Fig1Suppl_2-Full_fastMNN_IFNa_split.png", sep = "/"), plot = p, width = 20, height = 8, dpi = 100) ################## #### Figure 2 #### ################## # 1 - UMAP plot T cells res 0.6 with labels tobj <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/")) t_umap_plot <- cbind(tobj@reductions$umap@cell.embeddings, tobj@meta.data[rownames(tobj@reductions$umap@cell.embeddings),]) label <- data.frame( UMAP_1 = c(0.5, 0, -5, 0, 3.5, 4, 2.2, -2.5), UMAP_2 = c(1.5, -3.5, -0.3, 4.3, -3, 2.5, 0.5, 0.5), RNA_snn_res.0.6 = c(0, 1, 2, 3, 4, 5, 6, 7), label = c("CD4 1", "CD8 1", "Proliferating", "CD4 2", "CD8 2", "Naive", "T regs", "CD4 3") ) p <- ggplot(t_umap_plot, aes(UMAP_1, UMAP_2)) + geom_point(aes(colour = RNA_snn_res.0.6), size = 1.5, alpha = 1) + theme_nothing(font_size = 20) + #theme_bw(base_size = 20) + theme(legend.position = "none") + #geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) + geom_label(data = label, aes(label = label), size = 8) + #guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) + xlab("") + ylab("") ggsave(filename = paste(out_dir, "Fig2_1-TcellFiltered_fastMNN_Res06.png", sep = "/"), plot = p, width = 10, height = 8, dpi = 100) # 2 - GSEA DEG High vs. CTRL e Low vs. CTRL Tcells, di signature Exhaustion_GSEA_new (Wherry et al., 2007) tobj <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/")) t_gmt <- read.gmt(paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "PaperSig.gmt", sep = "/")) for (contr in c("CD_High-AB_Ctrl", "EF_Low-AB_Ctrl")) { deg_res <- tobj@misc$dgemarkers[[contr]] deg_res.logfc <- as.vector(deg_res$avg_logFC) names(deg_res.logfc) <- row.names(deg_res) deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE) deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))] print(deg_res.logfc[1:5]) # GSEA contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = t_gmt, nPerm = 10000, pvalueCutoff = 0.2) num_id <- 0 for (p in seq(1, nrow(contr.enricher.gsea@result))) { if (as.character(contr.enricher.gsea@result[p, "ID"]) == "Exhaustion") { num_id <- p } } if (contr == "CD_High-AB_Ctrl") { p <- gseaplot2(contr.enricher.gsea, geneSetID = num_id, base_size = 18, title = expression(IFN-alpha), color = c("red"), pvalue_table = T, nopadj = T) } else { p <- gseaplot2(contr.enricher.gsea, geneSetID = num_id, base_size = 18, title = expression(IFN-alpha-DHFR), color = c("red"), pvalue_table = T, nopadj = T) } ggsave(filename = paste(out_dir, paste0("Fig2_2-TcellFiltered_", contr, "_Exhaustion_GSEA.png"), sep = "/"), plot = p, width = 12, height = 8, dpi = 100) } # 3 - UMAP plot DC with labels res 0.2 dcobj <- readRDS(paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/")) dc_umap_plot <- cbind(dcobj@reductions$umap@cell.embeddings, dcobj@meta.data[rownames(dcobj@reductions$umap@cell.embeddings),]) dc_umap_plot$RNA_snn_res.0.1 <- factor(dc_umap_plot$RNA_snn_res.0.1, levels = c(0, 1, 2, 3)) #levels(dc_umap_plot$RNA_snn_res.0.1) <- c("0", "1", "2", "3") dc_umap_plot$RNA_snn_res.0.1 label <- data.frame( UMAP_1 = c(-5, 5, -1, 1), UMAP_2 = c(2.5, 0, -7.5, 10), RNA_snn_res.0.1 = c(0, 1, 2, 3), label = c("mregDC", "DC2", "DC1", "pDC") ) p <- ggplot(dc_umap_plot, aes(UMAP_1, UMAP_2)) + geom_point(aes(colour = RNA_snn_res.0.1), size = 1.5, alpha = 1) + theme_nothing(font_size = 20) + #theme_bw(base_size = 20) + theme(legend.position = "none") + #geom_text(aes(label = RNA_snn_res.0.1), check_overlap = TRUE) + geom_label(data = label, aes(label = label), size = 8) + #guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) + xlab("") + ylab("") ggsave(filename = paste(out_dir, "Fig2_3-DendriticFilteredNEW_fastMNN_Res01.png", sep = "/"), plot = p, width = 10, height = 8, dpi = 100) # 4/5 - GSEA DEG High vs. CTRL e Low vs. CTRL DC di signatures “Unstimulated DC vs. 4h LPS DN” e “Hallmark Oxidative phosphorylation” dcobj <- readRDS(paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/")) c7_gmt <- read.gmt(paste(wdir, "DendriticFilteredNEW", "02-fastMNN_2-post-analysis", "c7.all.v6.1.symbols_mouse.gmt", sep = "/")) for (contr in c("CD_High-AB_Ctrl", "EF_Low-AB_Ctrl")) { deg_res <- dcobj@misc$dgemarkers[[contr]] deg_res.logfc <- as.vector(deg_res$avg_logFC) names(deg_res.logfc) <- row.names(deg_res) deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE) deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))] print(deg_res.logfc[1:5]) # GSEA contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = c7_gmt, nPerm = 10000, pvalueCutoff = 0.2) num_id <- 0 for (p in seq(1, nrow(contr.enricher.gsea@result))) { if (as.character(contr.enricher.gsea@result[p, "ID"]) == "GSE14000_UNSTIM_VS_4H_LPS_DC_DN") { num_id <- p } } if (contr == "CD_High-AB_Ctrl") { p <- gseaplot2(contr.enricher.gsea, geneSetID = num_id, base_size = 18, title = expression(IFN-alpha), color = c("red"), pvalue_table = T, nopadj = T) } else { p <- gseaplot2(contr.enricher.gsea, geneSetID = num_id, base_size = 18, title = expression(IFN-alpha-DHFR), color = c("red"), pvalue_table = T, nopadj = T) } ggsave(filename = paste(out_dir, paste0("Fig2_4-DendriticFilteredNEW_", contr, "_GSE14000_UNSTIM_VS_4H_LPS_DC_DN_GSEA.png"), sep = "/"), plot = p, width = 12, height = 8, dpi = 100) } # HALLMARK_OXIDATIVE_PHOSPHORYLATION ds_high <- readRDS(paste(wdir, "DendriticFilteredNEW", "02-fastMNN_2-post-analysis", "CD_High-AB_Ctrl", "CD_High-AB_Ctrl-post-analysis.rds", sep = "/")) num_id <- 0 for (p in seq(1, nrow(ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse@result))) { if (as.character(ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse@result[p, "ID"]) == "HALLMARK_OXIDATIVE_PHOSPHORYLATION") { num_id <- p } } ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse@result[num_id,3:8] p <- gseaplot2(ds_high$results$gsea$gsea_h.all.v6.1.symbols_mouse, geneSetID = num_id, base_size = 18, title = expression(IFN-alpha), color = c("red"), pvalue_table = T, nopadj = T) ggsave(filename = paste(out_dir, "Fig2_5-DendriticFilteredNEW_CD_High-AB_Ctrl_HALLMARK_OXIDATIVE_PHOSPHORYLATION_GSEA.png", sep = "/"), plot = p, width = 12, height = 8, dpi = 100) # HALLMARK_OXIDATIVE_PHOSPHORYLATION ds_low <- readRDS(paste(wdir, "DendriticFilteredNEW", "02-fastMNN_2-post-analysis", "EF_Low-AB_Ctrl", "EF_Low-AB_Ctrl-post-analysis.rds", sep = "/")) num_id <- 0 for (p in seq(1, nrow(ds_low$results$gsea$gsea_h.all.v6.1.symbols_mouse@result))) { if (as.character(ds_low$results$gsea$gsea_h.all.v6.1.symbols_mouse@result[p, "ID"]) == "HALLMARK_OXIDATIVE_PHOSPHORYLATION") { num_id <- p } } ds_low$results$gsea$gsea_h.all.v6.1.symbols_mouse@result[num_id,3:8] p <- gseaplot2(ds_low$results$gsea$gsea_h.all.v6.1.symbols_mouse, geneSetID = num_id, base_size = 18, title = expression(IFN-alpha-DHFR), color = c("red"), pvalue_table = T, nopadj = T) ggsave(filename = paste(out_dir, "Fig2_5-DendriticFilteredNEW_EF_Low-AB_Ctrl_HALLMARK_OXIDATIVE_PHOSPHORYLATION_GSEA.png", sep = "/"), plot = p, width = 12, height = 8, dpi = 100) ######################### #### Figure 2 Suppl. #### ######################### # 1 - Heatmap marker genes Tcells with labels hm <- doSCHeatmapMean(fname = paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/"), clustering = "RNA_snn_res.0.6", ntop = 15, num_color = 3, high_genes = c("Furin","Cd4","Cxcr6","Bhlhe40","Gzma","Cd8a","Irf8","Prf1","Stmn1","Pclaf","Hist1h2ap", "Tubb5","Ccl1","Penk","Xcl1","Ifng","Cd7","Gzmk","Irf7","Eomes","S1pr1","Dusp10","Ramp3", "Tcf7","Ikzf2","Izumo1r","Foxp3","Ctla4","Nme1","Lgals1","Pgk1","Lgals7"), clu_lab = c("CD4 1", "CD8 1", "Proliferating", "CD4 2", "CD8 2", "Naive", "T regs", "CD4 3")) png(paste(out_dir, "Fig2Suppl_1-TcellFiltered_fastMNN_RNA_snn_res.0.6_top15-color3-wider2.png", sep = "/"), width = 1500, height = 700, res = 150) draw(hm) dev.off() # 2 - Heatmap marker genes DC with labels # obj <- readRDS(paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/")) # Idents(obj) <- "RNA_snn_res.0.1" # marks <- FindAllMarkers(object = obj, min.pct = 0.1, logfc.threshold = 0.25, only.pos = FALSE, # min.cells.group = 10, test.use = "wilcox", return.thresh = 1e-06) # obj@misc$markers[["RNA_snn_res.0.1"]] <- marks # saveRDS(object = obj, file = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/")) hm <- doSCHeatmapMean(fname = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/"), clustering = "RNA_snn_res.0.1", ntop = 15, num_color = 3, high_genes = c("Fscn1", "Tmem123", "Ccr7", "Ms4a6c", "S100a6", "Fcer1g", "Naaa", "Xcr1", "Cd24a", "Ly6d", "Ccr9", "Siglech"), clu_lab = c("mregDC", "DC2", "DC1", "pDC")) png(paste(out_dir, "Fig2Suppl_2-DendriticFilteredNEW_fastMNN_RNA_snn_res.0.1_top15-color3-wider2.png", sep = "/"), width = 1500, height = 700, res = 150) draw(hm) dev.off() ################## #### Figure 3 #### ################## # 1 - UMAP plot Mono/macrophages split by treatment, res 0.6 mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) m_umap_plot <- cbind(mobj@reductions$umap@cell.embeddings, mobj@meta.data[rownames(mobj@reductions$umap@cell.embeddings),]) label <- data.frame( UMAP_1 = c(-1.5, 0, 4, 3, -2.5, -4, -4, 6, -5.5, -1.8, 7.5), UMAP_2 = c(-1.5, 3, -4, 0, 1.5, -3, 3, -1.5, 3.5, 6, 0), RNA_snn_res.0.6 = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10") ) p <- ggplot(m_umap_plot, aes(UMAP_1, UMAP_2)) + geom_point(aes(colour = RNA_snn_res.0.6), size = 1.5, alpha = 1) + theme_nothing(font_size = 20) + #theme_bw(base_size = 20) + theme(legend.position = "none") + #geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) + geom_label(data = label, aes(label = label), size = 8) + #guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) + xlab("") + ylab("") ggsave(filename = paste(out_dir, "Fig3_1-MacrophagesFiltered_fastMNN_Res06.png", sep = "/"), plot = p, width = 10, height = 8, dpi = 100) # 2 - Vulcano plot DEG Gene Therapy vs Control hgenes <- c("Ccl8", "Pf4", "Ccl7", "Selenop", "Clec4a1", "Cd72", "Cd81", "Eno1", "Mrc1", "Msr1", "Ly6c2", "Isg15", "Ier3", "Iigp1", "Dusp1", "Ier5", "Jun", "Gm4951") res_dge <- read.xlsx(xlsxFile = paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "GeneTherapy-Control-DGE.xlsx", sep = "/"), sheet = "GeneTherapy-Control", rowNames = T) vol <- plotVolcano(dge_res = res_dge, adj.pval.thr = 0.01, logfc.thr = 0, hgenes = hgenes) ggsave(filename = paste(out_dir, "Fig3_2-MacrophagesFiltered_GeneTherapy_Ctrl_Volcano.png", sep = "/"), plot = vol, width = 10, height = 9, dpi = 100) # 3 - UMAP plot Mono/macrophages split by treatment, res 0.6 mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) m_umap_plot <- cbind(mobj@reductions$umap@cell.embeddings, mobj@meta.data[rownames(mobj@reductions$umap@cell.embeddings),]) m_umap_plot$RNA_Group2 <- factor(m_umap_plot$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low")) levels(m_umap_plot$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR") label <- data.frame( UMAP_1 = c(-1.5, 0, 4, 3, -2.5, -4, -4, 6, -5.5, -1.8, 7.5), UMAP_2 = c(-1.5, 3, -4, 0, 1.5, -3, 3, -1.5, 3.5, 6, 0), RNA_snn_res.0.6 = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10") ) p <- ggplot(m_umap_plot, aes(UMAP_1, UMAP_2)) + geom_point(aes(colour = RNA_snn_res.0.6), size = 1.5, alpha = 1) + theme_nothing(font_size = 20) + #theme_bw(base_size = 20) + theme(legend.position = "none", strip.text.x = element_text(size = 40)) + #geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) + geom_label(data = label, aes(label = label), size = 8) + #guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) + xlab("") + ylab("") + facet_grid(.~RNA_Group2, labeller = label_parsed) ggsave(filename = paste(out_dir, "Fig3_3-MacrophagesFiltered_fastMNN_Res06_split.png", sep = "/"), plot = p, width = 28, height = 8, dpi = 100) # 4 - Featureplot genes # genes <- c("Pglyrp1", "Ace", "Plac8", "Ifitm6", "Stmn1", "Top2a", "Chil3", "Arg1", "Il1b", "Dusp1", "Tnf", "Ccl8", "Pf4", "Mrc1") # png(filename = paste(out_dir, "Fig3_4-MacrophagesFiltered_fastMNN_FeatureplotGenes.png", sep = "/"), width = 1200, height = 1200, res = 100) # FeaturePlot(object = mobj, features = genes, reduction = "umap", pt.size = 1, ncol = 4, cols = c("blue", "white", "red"), order = T) + NoLegend() # dev.off() ######################### #### Figure 3 Suppl. #### ######################### # 1 - GeneTherapy Dotplot of enrichment pos HALLMARK gt_enr_pos <- read.xlsx(xlsxFile = paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "GeneTherapy-Control_enrichment_pos.xlsx", sep = "/"), sheet = "enr_pos_HALLMARK") ds_hall_enrich <- gt_enr_pos[,c("ID", "GeneRatio", "p.adjust", "Count")] gsize <- as.numeric(sub("/\\d+$", "", as.character(ds_hall_enrich$GeneRatio))) gcsize <- as.numeric(sub("^\\d+/", "", as.character(ds_hall_enrich$GeneRatio))) ds_hall_enrich$GeneRatio = gsize/gcsize ds_hall_enrich <- subset(ds_hall_enrich, p.adjust < 0.05) ds_hall_enrich$ID <- gsub("HALLMARK_", "", ds_hall_enrich$ID) #ds_hall_enrich$ID <- tolower(gsub("_", " ", ds_hall_enrich$ID)) print(head(ds_hall_enrich)) id_order <- ds_hall_enrich %>% select(ID, GeneRatio) %>% group_by(ID) %>% summarize_all(sum) %>% data.frame() id_order <- id_order[order(id_order$GeneRatio, decreasing = F),] ds_hall_enrich$ID <- factor(ds_hall_enrich$ID, levels = id_order$ID) p <- ggplot(ds_hall_enrich, aes(x = GeneRatio, y = ID, size = Count, color = p.adjust)) + theme_bw(base_size = 20) + scale_color_gradient(low = "red", high = "blue", guide = guide_colorbar(reverse = T)) + geom_point() + ylab("") ggsave(plot = p, filename = paste(out_dir, "Fig3Suppl_1-MacrophagesFiltered_GT_Ctrl_Hallmark_enrichment_pos.png", sep = "/"), width = 10, height = 8, dpi = 150) # 2 - GSEA MacroPolarization macro_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) gmt.obj <- read.gmt("/TIGET/coltella_sc_new/reference/MacroPolarizationSig.gmt") # DGE Marker for (contr in c("EF_Low-AB_Ctrl", "CD_High-AB_Ctrl")) { deg_res <- macro_obj@misc$dgemarkers[[contr]] deg_res.logfc <- as.vector(deg_res$avg_logFC) names(deg_res.logfc) <- row.names(deg_res) deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE) deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))] print(deg_res.logfc[1:5]) # GSEA contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = gmt.obj, nPerm = 10000) contr.enricher.gsea@result[,1:10] ptitle <- expression(IFN-alpha) if (contr == "EF_Low-AB_Ctrl") { ptitle <- expression(IFN-alpha-DHFR) } p <- gseaplot2(contr.enricher.gsea, geneSetID = c(1,2), base_size = 18, title = ptitle, color = c("red", "forestgreen"), pvalue_table = T, nopadj = T, labels = c("M1 Macro. IL-4", bquote("M2 Macro. IFN"~gamma~"+ LPS"))) ggsave(filename = paste(out_dir, paste0("Fig3Suppl_2-MacrophagesFiltered_", contr, "_MacroPolarization_GSEA.png"), sep = "/"), plot = p, width = 16, height = 8, dpi = 100) } # 3 - Heatmap marker genes Macrophages with labels hm <- doSCHeatmapMean(fname = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"), clustering = "RNA_snn_res.0.6", ntop = 15, num_color = 3, high_genes = c("Saa3", "Ccl5", "Ly6i", "Tmem176b", "Cx3cr1", "Hexb", "Dusp1", "Socs3", "Il1b", "Plac8", "Il1r2", "Mgl2", "Ccl8", "Cd63", "Il18bp", "Chil3", "Arg1", "Ccl24", "Pf4", "Mrc1", "Ms4a4c", "Ly6c2", "Cd74", "C1qb", "Stmn1", "Pclaf", "Birc5", "Pglyrp1", "Gsr", "Clec4e"), clu_lab = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10"), duplicated_labs = T) png(paste(out_dir, "Fig3Suppl_3-MacrophagesFiltered_fastMNN_RNA_snn_res.0.6_top15-color3-wider2.png", sep = "/"), width = 1800, height = 700, res = 150) draw(hm) dev.off() # 4 - Feature plot hypoxia (modulescore) from https://www.gsea-msigdb.org/gsea/msigdb/cards/SEMENZA_HIF1_TARGETS mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) hypoxia_sig <- c("Adm","Adra1b","Ak3","Aldoart2","Aldoc","Bnip3","Car9","Cdkn1a","Cited2","Cp","Edn1","Eno1","Epo","Flt1","Gm20425", "Hk1","Hk2","Hmox1","Igf2","Igfbp1","Igfbp2","Igfbp3","Ldha","Nos2","P4ha2","Pfkl","Pgk1","Serpine1","Slc2a1","Slc2a3", "Tfrc","Tgfb3","Tpi1","Vegfa") hypoxia_sig_filt <- hypoxia_sig[hypoxia_sig %in% rownames(mobj)] mobj <- AddModuleScore(object = mobj, features = list(HYPOXIA = hypoxia_sig_filt), ctrl = length(hypoxia_sig_filt), name = "HYPOXIA") macro_umap_plot <- cbind(mobj@reductions$umap@cell.embeddings, mobj@meta.data[rownames(mobj@reductions$umap@cell.embeddings),]) macro_umap_plot$RNA_Group2 <- factor(macro_umap_plot$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low")) levels(macro_umap_plot$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR") str(umap_plot) p <- ggplot(macro_umap_plot, aes(x = UMAP_1, y = UMAP_2)) + #geom_point(aes(colour = HYPOXIA1), size = 1, alpha = 0.8) + geom_point(data = subset(macro_umap_plot, HYPOXIA1 >= -0.25 & HYPOXIA1 < 0.25), aes(colour = HYPOXIA1), size = 1.5, alpha = 0.8) + geom_point(data = subset(macro_umap_plot, HYPOXIA1 < -0.25), aes(colour = HYPOXIA1), size = 1.5, alpha = 1) + geom_point(data = subset(macro_umap_plot, HYPOXIA1 >= 0.25 & HYPOXIA1 < 0.5), aes(colour = HYPOXIA1), size = 1.5, alpha = 1) + geom_point(data = subset(macro_umap_plot, HYPOXIA1 >= 0.5), aes(colour = HYPOXIA1), size = 1.5, alpha = 1) + # scale_colour_gradient2(low="#2166ac", # mid="gray95", # high="#b2182b", # name = "Expr.", # midpoint = -0.05) + scale_color_gradientn(colours = c("#2166AC","#4393C3","#D1E5F0","#F4A582","#B2182B"), name = "Expr.", #values = rescale(c(-.1,mean(macro_umap_plot$IFNa1),.9)), values = rescale(c(-.5,-.4,-.1,.2,.7)), guide = "colorbar") + theme_void(base_size = 20) + #ggtitle("INTERFERON ALPHA RESPONSE") + theme(legend.position = "right", strip.text.x = element_text(size = 40)) #xlab("UMAP 1") + #ylab("UMAP 2") + #facet_grid(.~RNA_Group2, labeller = label_parsed) ggsave(filename = paste(out_dir, "Fig3Suppl_4-MacrophagesFiltered_Hypoxia.png", sep = "/"), plot = p, width = 10, height = 8, dpi = 100) ################## #### Figure 4 #### ################## # 1 - Monocle pseudotime cds_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "03-fastMNN_2-pseudotime", "Monocle_cds_objects.rds", sep = "/")) label <- data.frame( UMAP_1 = c(-1.5, 0, 4, 3, -2.5, -4, -4, 6, -5.5, -1.8, 7.5), UMAP_2 = c(-1.5, 3, -4, 0, 1.5, -3, 3, -1.5, 3.5, 6, 0), RNA_snn_res.0.6 = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10") ) # # Single plots # prepare_plot <- function(cds) { # S_matrix <- reducedDims(cds)$UMAP # data_df <- data.frame(S_matrix[,c(1,2)]) # data_df$sample_name <- row.names(data_df) # data_df <- as.data.frame(cbind(data_df, colData(cds))) # data_df$RNA_Group2 <- factor(data_df$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low")) # levels(data_df$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR") # #str(data_df) # ica_space_df <- t(cds@principal_graph_aux$UMAP$dp_mst) %>% # as.data.frame() %>% # dplyr::select_(prin_graph_dim_1 = 1, prin_graph_dim_2 = 2) %>% # dplyr::mutate(sample_name = rownames(.), # sample_state = rownames(.)) # dp_mst <- cds@principal_graph$UMAP # edge_df <- dp_mst %>% # igraph::as_data_frame() %>% # dplyr::select_(source = "from", target = "to") %>% # dplyr::left_join(ica_space_df %>% # dplyr::select_( # source="sample_name", # source_prin_graph_dim_1="prin_graph_dim_1", # source_prin_graph_dim_2="prin_graph_dim_2"), # by = "source") %>% # dplyr::left_join(ica_space_df %>% # dplyr::select_( # target="sample_name", # target_prin_graph_dim_1="prin_graph_dim_1", # target_prin_graph_dim_2="prin_graph_dim_2"), # by = "target") # p <- ggplot(data_df, aes(x = UMAP_1, y = UMAP_2)) + # geom_point(aes(colour = RNA_snn_res.0.6), size = 1.5, alpha = 1) + # theme_nothing(font_size = 20) + # #theme_bw(base_size = 20) + # theme(legend.position = "none", # strip.text.x = element_text(size = 40)) + # geom_segment(data = edge_df, # aes_string(x="source_prin_graph_dim_1", # y="source_prin_graph_dim_2", # xend="target_prin_graph_dim_1", # yend="target_prin_graph_dim_2"), # size = 1, # color = "black", # linetype = "solid", # na.rm = TRUE) + # #geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) + # geom_label(data = label, aes(label = label), size = 8) + # #guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) + # xlab("") + # ylab("") + # facet_grid(.~RNA_Group2, labeller = label_parsed) # return(p) # } # p <- prepare_plot(cds_obj$AB_Ctrl) # ggsave(filename = paste(out_dir, "Fig4_1-MacrophagesFiltered_Ctrl_Monocle_RNA_snn_res.0.6_umap_graph.png", sep = "/"), # plot = p, width = 10, height = 8, dpi = 100) # p <- prepare_plot(cds_obj$CD_High) # ggsave(filename = paste(out_dir, "Fig4_1-MacrophagesFiltered_High_Monocle_RNA_snn_res.0.6_umap_graph.png", sep = "/"), # plot = p, width = 10, height = 8, dpi = 100) # p <- prepare_plot(cds_obj$EF_Low) # ggsave(filename = paste(out_dir, "Fig4_1-MacrophagesFiltered_Low_Monocle_RNA_snn_res.0.6_umap_graph.png", sep = "/"), # plot = p, width = 10, height = 8, dpi = 100) # UMAP Data data_df_ctrl <- data.frame(reducedDims(cds_obj$AB_Ctrl)$UMAP[,c(1,2)]) data_df_ctrl$sample_name <- row.names(data_df_ctrl) data_df_ctrl <- as.data.frame(cbind(data_df_ctrl, colData(cds_obj$AB_Ctrl))) data_df_high <- data.frame(reducedDims(cds_obj$CD_High)$UMAP[,c(1,2)]) data_df_high$sample_name <- row.names(data_df_high) data_df_high <- as.data.frame(cbind(data_df_high, colData(cds_obj$CD_High))) data_df_low <- data.frame(reducedDims(cds_obj$EF_Low)$UMAP[,c(1,2)]) data_df_low$sample_name <- row.names(data_df_low) data_df_low <- as.data.frame(cbind(data_df_low, colData(cds_obj$EF_Low))) data_df <- rbind(data_df_ctrl, rbind(data_df_high, data_df_low)) data_df$RNA_Group2 <- factor(data_df$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low")) levels(data_df$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR") # Graph Data ica_space_df_ctrl <- t(cds_obj$AB_Ctrl@principal_graph_aux$UMAP$dp_mst) %>% as.data.frame() %>% dplyr::select_(prin_graph_dim_1 = 1, prin_graph_dim_2 = 2) %>% dplyr::mutate(sample_name = rownames(.), sample_state = rownames(.)) dp_mst_ctrl <- cds_obj$AB_Ctrl@principal_graph$UMAP edge_df_ctrl <- dp_mst_ctrl %>% igraph::as_data_frame() %>% dplyr::select_(source = "from", target = "to") %>% dplyr::left_join(ica_space_df_ctrl %>% dplyr::select_( source="sample_name", source_prin_graph_dim_1="prin_graph_dim_1", source_prin_graph_dim_2="prin_graph_dim_2"), by = "source") %>% dplyr::left_join(ica_space_df_ctrl %>% dplyr::select_( target="sample_name", target_prin_graph_dim_1="prin_graph_dim_1", target_prin_graph_dim_2="prin_graph_dim_2"), by = "target") ica_space_df_high <- t(cds_obj$CD_High@principal_graph_aux$UMAP$dp_mst) %>% as.data.frame() %>% dplyr::select_(prin_graph_dim_1 = 1, prin_graph_dim_2 = 2) %>% dplyr::mutate(sample_name = rownames(.), sample_state = rownames(.)) dp_mst_high <- cds_obj$CD_High@principal_graph$UMAP edge_df_high <- dp_mst_high %>% igraph::as_data_frame() %>% dplyr::select_(source = "from", target = "to") %>% dplyr::left_join(ica_space_df_high %>% dplyr::select_( source="sample_name", source_prin_graph_dim_1="prin_graph_dim_1", source_prin_graph_dim_2="prin_graph_dim_2"), by = "source") %>% dplyr::left_join(ica_space_df_high %>% dplyr::select_( target="sample_name", target_prin_graph_dim_1="prin_graph_dim_1", target_prin_graph_dim_2="prin_graph_dim_2"), by = "target") ica_space_df_low <- t(cds_obj$EF_Low@principal_graph_aux$UMAP$dp_mst) %>% as.data.frame() %>% dplyr::select_(prin_graph_dim_1 = 1, prin_graph_dim_2 = 2) %>% dplyr::mutate(sample_name = rownames(.), sample_state = rownames(.)) dp_mst_low <- cds_obj$EF_Low@principal_graph$UMAP edge_df_low <- dp_mst_low %>% igraph::as_data_frame() %>% dplyr::select_(source = "from", target = "to") %>% dplyr::left_join(ica_space_df_low %>% dplyr::select_( source="sample_name", source_prin_graph_dim_1="prin_graph_dim_1", source_prin_graph_dim_2="prin_graph_dim_2"), by = "source") %>% dplyr::left_join(ica_space_df_low %>% dplyr::select_( target="sample_name", target_prin_graph_dim_1="prin_graph_dim_1", target_prin_graph_dim_2="prin_graph_dim_2"), by = "target") edge_df_ctrl$RNA_Group <- "AB_Ctrl" edge_df_high$RNA_Group <- "CD_High" edge_df_low$RNA_Group <- "EF_Low" edge_df <- rbind(edge_df_ctrl, rbind(edge_df_high, edge_df_low)) edge_df$RNA_Group2 <- factor(edge_df$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low")) levels(edge_df$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR") p <- ggplot(data_df, aes(x = UMAP_1, y = UMAP_2)) + geom_point(aes(colour = RNA_snn_res.0.6), size = 1.5, alpha = 1) + theme_nothing(font_size = 20) + #theme_bw(base_size = 20) + theme(legend.position = "none", strip.text.x = element_text(size = 40)) + geom_segment(data = edge_df, aes_string(x="source_prin_graph_dim_1", y="source_prin_graph_dim_2", xend="target_prin_graph_dim_1", yend="target_prin_graph_dim_2"), size = 1, color = "black", linetype = "solid", na.rm = TRUE, inherit.aes=FALSE) + #geom_text(aes(label = RNA_snn_res.0.6), check_overlap = TRUE) + geom_label(data = label, aes(label = label), size = 8) + #guides(colour = guide_legend(ncol = 1, title = "Res. 0.6", override.aes = list(size = 3))) + xlab("") + ylab("") + facet_grid(.~RNA_Group2, labeller = label_parsed) ggsave(filename = paste(out_dir, "Fig4_1-MacrophagesFiltered_Monocle_RNA_snn_res.0.6_umap_graph.png", sep = "/"), plot = p, width = 28, height = 8, dpi = 100) # 3 - Feature plot IFN_ALPHA (modulescore). mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) inf_alpha <- c("1700057G04Rik","AC125149.4","AC125149.5","AC132444.3","AC133103.4","AC133103.5","AC168977.1","Adar","B2m","Batf2","Bst2","C1s1","C1s2","Casp1", "Casp8","Ccrl2","Cd47","Cd74","Cmpk2","Cmtr1","Cnp","Csf1","Cxcl10","Cxcl11","Ddx60","Dhx58","Eif2ak2","Elf1","Epsti1","Gbp2","Gbp2b","Gmpr", "Helz2","Herc6","Ifi35","Ifi44","Ifi44l","Ifih1","Ifit2","Ifit3","Ifit3b","Ifitm1","Ifitm2","Ifitm3","Ifitm6","Ifitm7","Il15","Il4ra","Il7", "Irf1","Irf2","Irf7","Irf9","Isg15","Isg20","Lap3","Lgals3bp","Lpar6","Ly6e","Mov10","Mvb12a","Mx1","Mx2","Ncoa7","Nmi","Nub1","Oas1a","Oas1b", "Oas1c","Oas1d","Oas1e","Oas1f","Oas1g","Oas1h","Oasl1","Ogfr","Parp12","Parp14","Parp9","Plscr1","Plscr2","Pnpt1","Procr","Psma3","Psmb8","Psmb9", "Psme1","Psme2","Psme2b","Ripk2","Rnf31","Rsad2","Rtp4","Samd9l","Sell","Slc25a28","Sp110","Stat2","Tap1","Tdrd7","Tent5a","Tmem140","Trafd1", "Trim12a","Trim12c","Trim14","Trim21","Trim25","Trim26","Trim30a","Trim30c","Trim30d","Trim5","Txnip","Uba7","Ube2l6","Usp18","Wars") inf_alpha_filt <- inf_alpha[inf_alpha %in% rownames(mobj)] mobj <- AddModuleScore(object = mobj, features = list(IFNa = inf_alpha_filt), ctrl = length(inf_alpha_filt), name = "IFNa") macro_umap_plot <- cbind(mobj@reductions$umap@cell.embeddings, mobj@meta.data[rownames(mobj@reductions$umap@cell.embeddings),]) macro_umap_plot$RNA_Group2 <- factor(macro_umap_plot$RNA_Group, levels = c("AB_Ctrl", "CD_High", "EF_Low")) levels(macro_umap_plot$RNA_Group2) <- c("Control", "IFN-alpha", "IFN-alpha-DHFR") summary(macro_umap_plot$IFNa1) p <- ggplot(macro_umap_plot, aes(x = UMAP_1, y = UMAP_2)) + #geom_point(aes(colour = IFNa1), size = 1.5, alpha = 1) + geom_point(data = subset(macro_umap_plot, IFNa1 >= 0.1 & IFNa1 < 0.3), aes(colour = IFNa1), size = 2, alpha = 0.8) + geom_point(data = subset(macro_umap_plot, IFNa1 < 0.1), aes(colour = IFNa1), size = 2, alpha = 1) + #geom_point(data = subset(macro_umap_plot, IFNa1 >= 0.25 & IFNa1 < 0.5), aes(colour = IFNa1), size = 1.5, alpha = 1) + geom_point(data = subset(macro_umap_plot, IFNa1 >= 0.3), aes(colour = IFNa1), size = 2, alpha = 1) + # scale_colour_gradient2(low="#2166ac", # mid="gray95", # high="#b2182b", # name = "Expr.", # midpoint = 0.25) + # scale_color_gradientn(colours = c("blue","blue","grey95","red","red"), # #values = rescale(c(-.1,mean(macro_umap_plot$IFNa1),.9)), # values = rescale(c(-.1,0,.2,.5,.65)), # guide = "colorbar") + scale_color_gradientn(colours = c("#2166AC","#4393C3","#D1E5F0","#F4A582","#B2182B"), name = "Expr.", #values = rescale(c(-.1,mean(macro_umap_plot$IFNa1),.9)), values = rescale(c(-.1,0,.2,.35,.65)), guide = "colorbar") + theme_void(base_size = 20) + #ggtitle("INTERFERON ALPHA RESPONSE") + theme(legend.position = "right", strip.text.x = element_text(size = 40)) + #xlab("UMAP 1") + #ylab("UMAP 2") + facet_grid(.~RNA_Group2, labeller = label_parsed) ggsave(filename = paste(out_dir, "Fig4_3-MacrophagesFiltered_IFNa.png", sep = "/"), plot = p, width = 28, height = 8, dpi = 100) # # 2 - GSEA signatures M1 e M2 DEG High vs. CTRL e Low vs. CTRL (macro) # #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/CD_High-AB_Ctrl/Alternatively_Activated_IL-4_Signature_GSEA) # #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/CD_High-AB_Ctrl/Classically_Activated_LPS-IFN_Signature_GSEA) # #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/EF_Low-AB_Ctrl/Alternatively_Activated_IL-4_Signature_GSEA) # #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/EF_Low-AB_Ctrl/Classically_Activated_LPS-IFN_Signature_GSEA) # mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) # pol_gmt <- read.gmt(paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "MacroPolarizationSig.gmt", sep = "/")) # for (contr in c("CD_High-AB_Ctrl", "EF_Low-AB_Ctrl")) { # deg_res <- mobj@misc$dgemarkers[[contr]] # deg_res.logfc <- as.vector(deg_res$avg_logFC) # names(deg_res.logfc) <- row.names(deg_res) # deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE) # deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))] # print(deg_res.logfc[1:5]) # # GSEA # contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = pol_gmt, nPerm = 10000, pvalueCutoff = 0.2) # for (p in seq(1, nrow(contr.enricher.gsea@result))) { # d <- as.character(contr.enricher.gsea@result[p, "ID"]) # if (contr == "CD_High-AB_Ctrl") { # p <- gseaplot2(contr.enricher.gsea, geneSetID = p, base_size = 18, # title = bquote(.(d)~": " ~ IFN-alpha ~ " vs. Control"), # color = c("red"), # pvalue_table = T) # } else { # p <- gseaplot2(contr.enricher.gsea, geneSetID = p, base_size = 18, # title = bquote(.(d)~": " ~ IFN-alpha-DHFR ~ " vs. Control"), # color = c("red"), # pvalue_table = T) # } # ggsave(filename = paste(out_dir, paste0("Fig4_2_MacrophagesFiltered_", contr, "_", d, "_GSEA.png"), sep = "/"), plot = p, width = 12, height = 8, dpi = 100) # } # } # # 4 - GSEA signatures M1 e M2 on pos marker genes of cluster 6 (res 0.6) # #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/ RNA_snn_res.0.6_Macro_polarization_GSEA/ clu6_Alternatively_Activated_IL-4_Signature_GSEA) # #(coltella_sc/Harmony/ MacrophagesFiltered/02-fastMNN_2-post-analysis/ RNA_snn_res.0.6_Macro_polarization_GSEA/clu6_Classically_Activated_LPS-IFN_Signature_GSEA) # mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) # pol_gmt <- read.gmt(paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "MacroPolarizationSig.gmt", sep = "/")) # deg_res <- subset(mobj@misc$markers_full$RNA_snn_res.0.6, cluster == "6") # deg_res.logfc <- as.vector(deg_res$avg_logFC) # names(deg_res.logfc) <- deg_res$gene # deg_res.logfc <- sort(deg_res.logfc, decreasing = TRUE) # deg_res.logfc <- deg_res.logfc[!duplicated(names(deg_res.logfc))] # print(deg_res.logfc[1:5]) # # GSEA # contr.enricher.gsea <- GSEA(deg_res.logfc, TERM2GENE = pol_gmt, nPerm = 10000, pvalueCutoff = 0.2) # for (cc in c(1, 2)) { # d <- as.character(contr.enricher.gsea@result[cc, "ID"]) # p <- gseaplot2(contr.enricher.gsea, geneSetID = cc, base_size = 18, # title = bquote(.(d)), # color = c("red"), # pvalue_table = T) # ggsave(filename = paste(out_dir, paste0("Fig4_4_MacrophagesFiltered_Clu6_", d, "_GSEA.png"), sep = "/"), plot = p, width = 12, height = 8, dpi = 100) # } # # 3 - Featureplot Ccl8, Mrc1 # mobj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) # ccl8_expr_plot <- cbind(as.data.frame(mobj@reductions$umap@cell.embeddings), # "Expr" = as.numeric(mobj@assays$RNA@data["Ccl8", rownames(mobj@reductions$umap@cell.embeddings)])) # p <- ggplot(ccl8_expr_plot, aes(x = UMAP_1, y = UMAP_2)) + # geom_point(aes(colour = Expr), size = 1, alpha = 0.8) + # scale_colour_gradient2(low="#2166ac", # mid="grey", # high="#b2182b", # name = "Expr.", # midpoint = (((max(ccl8_expr_plot$Expr)-min(ccl8_expr_plot$Expr))/2)-0.7)) + # theme_bw(base_size = 20) + # #ggtitle("Ccl8 Expression") + # theme(legend.position = "right") + # xlab("UMAP 1") + # ylab("UMAP 2") # ggsave(filename = paste(out_dir, "Fig5_3_MacrophagesFiltered_Ccl8_Expr.png", sep = "/"), # plot = p, width = 10, height = 8, dpi = 100) # mrc1_expr_plot <- cbind(as.data.frame(mobj@reductions$umap@cell.embeddings), # "Expr" = mobj@assays$RNA@data["Mrc1", rownames(mobj@reductions$umap@cell.embeddings)]) # p <- ggplot(mrc1_expr_plot, aes(x = UMAP_1, y = UMAP_2)) + # geom_point(data = subset(mrc1_expr_plot, Expr < 1), aes(colour = Expr), size = 1, alpha = 0.8) + # geom_point(data = subset(mrc1_expr_plot, Expr >= 1), aes(colour = Expr), size = 1.1, alpha = 1) + # scale_colour_gradient2(low="#2166ac", # mid="grey", # high="#b2182b", # name = "Expr.", # midpoint = (((max(mrc1_expr_plot$Expr)-min(mrc1_expr_plot$Expr))/2)-0.4)) + # theme_bw(base_size = 20) + # #ggtitle("Mrc1 Expression") + # theme(legend.position = "right") + # xlab("UMAP 1") + # ylab("UMAP 2") # ggsave(filename = paste(out_dir, "Fig5_3_MacrophagesFiltered_Mrc1_Expr.png", sep = "/"), # plot = p, width = 10, height = 8, dpi = 100)