library(Seurat) library(SeuratWrappers) library(ggplot2) library(dplyr) library(ggrepel) library(reshape) library(openxlsx) library(stringr) library(RColorBrewer) library(SingleR) library(ComplexHeatmap) library(circlize) plotVolcano <- function(dge_res, adj.pval.thr, logfc.thr) { 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),] highl <- head(subset(data, color != "NonSignificant"), 20) vol <- ggplot(data, aes(x = lfc, y = pval, colour = color)) + theme_bw(base_size = 12) + theme(legend.position = "right") + 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 = highl, aes(label = gene), show.legend = FALSE) return(vol) } wdir <- "Birocchi_SciTraMed2022/results" adj.pval.thr <- 1e-06 logfc.thr <- 0 # Full fmnn_full <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/")) Idents(fmnn_full) <- "RNA_Group" fmnn_full_mark_low <- FindMarkers(object = fmnn_full, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_full_mark_high <- FindMarkers(object = fmnn_full, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_full_mark_highlow <- FindMarkers(object = fmnn_full, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) fmnn_full_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_full_mark_low, "CD_High-AB_Ctrl" = fmnn_full_mark_high, "CD_High-EF_Low" = fmnn_full_mark_highlow) fmnn_full@misc[["dgemarkers"]] <- fmnn_full_dge_mark write.xlsx(x = fmnn_full_dge_mark, file = paste(wdir, "Full", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE) saveRDS(object = fmnn_full, file = paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/")) for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) { res_dge <- fmnn_full_dge_mark[[contr]] print(head(res_dge)) vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr) ggsave(filename = paste(wdir, "Full", "02-fastMNN_2-post-analysis", paste(contr, "volcano.png", sep = "-"), sep = "/"), plot = vol, width = 10, height = 9, dpi = 600) } # Macrophages fmnn_macro <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) Idents(fmnn_macro) <- "RNA_Group" fmnn_macro_mark_low <- FindMarkers(object = fmnn_macro, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_macro_mark_high <- FindMarkers(object = fmnn_macro, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_macro_mark_highlow <- FindMarkers(object = fmnn_macro, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) fmnn_macro_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_macro_mark_low, "CD_High-AB_Ctrl" = fmnn_macro_mark_high, "CD_High-EF_Low" = fmnn_macro_mark_highlow) fmnn_macro@misc[["dgemarkers"]] <- fmnn_macro_dge_mark write.xlsx(x = fmnn_macro_dge_mark, file = paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE) saveRDS(object = fmnn_macro, file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) { res_dge <- fmnn_macro_dge_mark[[contr]] print(head(res_dge)) vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr) ggsave(filename = paste(wdir, "MacrophagesFiltered", "02-fastMNN_2-post-analysis", paste(contr, "volcano.png", sep = "-"), sep = "/"), plot = vol, width = 10, height = 9, dpi = 600) } # T cell fmnn_tcell <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/")) Idents(fmnn_tcell) <- "RNA_Group" fmnn_tcell_mark_low <- FindMarkers(object = fmnn_tcell, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_tcell_mark_high <- FindMarkers(object = fmnn_tcell, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_tcell_mark_highlow <- FindMarkers(object = fmnn_tcell, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) fmnn_tcell_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_tcell_mark_low, "CD_High-AB_Ctrl" = fmnn_tcell_mark_high, "CD_High-EF_Low" = fmnn_tcell_mark_highlow) fmnn_tcell@misc[["dgemarkers"]] <- fmnn_tcell_dge_mark write.xlsx(x = fmnn_tcell_dge_mark, file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE) saveRDS(object = fmnn_tcell, file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/")) for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) { res_dge <- fmnn_tcell_dge_mark[[contr]] print(head(res_dge)) vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr) ggsave(filename = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", paste(contr, "volcano.png", sep = "-"), sep = "/"), plot = vol, width = 10, height = 9, dpi = 600) } # Cd4 - T cell fmnn_cd4 <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_final.rds", sep = "/")) Idents(fmnn_cd4) <- "RNA_Group" fmnn_cd4_mark_low <- FindMarkers(object = fmnn_cd4, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_cd4_mark_high <- FindMarkers(object = fmnn_cd4, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_cd4_mark_highlow <- FindMarkers(object = fmnn_cd4, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) fmnn_cd4_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_cd4_mark_low, "CD_High-AB_Ctrl" = fmnn_cd4_mark_high, "CD_High-EF_Low" = fmnn_cd4_mark_highlow) fmnn_cd4@misc[["dgemarkers"]] <- fmnn_cd4_dge_mark write.xlsx(x = fmnn_cd4_dge_mark, file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "Cd4sub", "DGE_markers.xlsx", sep = "/"), row.names = TRUE) saveRDS(object = fmnn_cd4, file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_final.rds", sep = "/")) for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) { res_dge <- fmnn_cd4_dge_mark[[contr]] print(head(res_dge)) vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr) ggsave(filename = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "Cd4sub", paste(contr, "volcano.png", sep = "-"), sep = "/"), plot = vol, width = 10, height = 9, dpi = 600) } # Cd8 - T cell fmnn_cd8 <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_final.rds", sep = "/")) Idents(fmnn_cd8) <- "RNA_Group" fmnn_cd8_mark_low <- FindMarkers(object = fmnn_cd8, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_cd8_mark_high <- FindMarkers(object = fmnn_cd8, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_cd8_mark_highlow <- FindMarkers(object = fmnn_cd8, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) fmnn_cd8_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_cd8_mark_low, "CD_High-AB_Ctrl" = fmnn_cd8_mark_high, "CD_High-EF_Low" = fmnn_cd8_mark_highlow) fmnn_cd8@misc[["dgemarkers"]] <- fmnn_cd8_dge_mark write.xlsx(x = fmnn_cd8_dge_mark, file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "Cd8sub", "DGE_markers.xlsx", sep = "/"), row.names = TRUE) saveRDS(object = fmnn_cd8, file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_final.rds", sep = "/")) for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) { res_dge <- fmnn_cd8_dge_mark[[contr]] print(head(res_dge)) vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr) ggsave(filename = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "Cd8sub", paste(contr, "volcano.png", sep = "-"), sep = "/"), plot = vol, width = 10, height = 9, dpi = 600) } # Dendritic NEW fmnn_dc <- readRDS(paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/")) Idents(fmnn_dc) <- "RNA_Group" fmnn_dc_mark_low <- FindMarkers(object = fmnn_dc, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_dc_mark_high <- FindMarkers(object = fmnn_dc, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_dc_mark_highlow <- FindMarkers(object = fmnn_dc, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) fmnn_dc_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_dc_mark_low, "CD_High-AB_Ctrl" = fmnn_dc_mark_high, "CD_High-EF_Low" = fmnn_dc_mark_highlow) fmnn_dc@misc[["dgemarkers"]] <- fmnn_dc_dge_mark write.xlsx(x = fmnn_dc_dge_mark, file = paste(wdir, "DendriticFilteredNEW", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE) saveRDS(object = fmnn_dc, file = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/")) for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) { res_dge <- fmnn_dc_dge_mark[[contr]] print(head(res_dge)) vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr) ggsave(filename = paste(wdir, "DendriticFilteredNEW", "02-fastMNN_2-post-analysis", paste(contr, "volcano.png", sep = "-"), sep = "/"), plot = vol, width = 10, height = 9, dpi = 600) } # T cell fmnn_tc <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/")) Idents(fmnn_tc) <- "RNA_Group" fmnn_tc_mark_low <- FindMarkers(object = fmnn_tc, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_tc_mark_high <- FindMarkers(object = fmnn_tc, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_tc_mark_highlow <- FindMarkers(object = fmnn_tc, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) fmnn_tc_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_tc_mark_low, "CD_High-AB_Ctrl" = fmnn_tc_mark_high, "CD_High-EF_Low" = fmnn_tc_mark_highlow) fmnn_tc@misc[["dgemarkers"]] <- fmnn_tc_dge_mark write.xlsx(x = fmnn_tc_dge_mark, file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE) saveRDS(object = fmnn_tc, file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/")) for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) { res_dge <- fmnn_tc_dge_mark[[contr]] print(head(res_dge)) vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr) ggsave(filename = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", paste(contr, "volcano.png", sep = "-"), sep = "/"), plot = vol, width = 10, height = 9, dpi = 600) } # T cell fmnn_tcell <- readRDS(paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/")) dim(fmnn_tcell) cd4_ssub <- subset(x = fmnn_tcell, Cd4 > Cd8a & Cd4 > 0) dim(cd4_ssub) FeaturePlot(object = cd4_ssub, features = c("Cd4", "Cd8a"), reduction = "umap", split.by = "RNA_Group") Idents(cd4_ssub) <- "RNA_Group" cd4_ssub_mark_low <- FindMarkers(object = cd4_ssub, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) cd4_ssub_mark_high <- FindMarkers(object = cd4_ssub, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) cd4_ssub_mark_highlow <- FindMarkers(object = cd4_ssub, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) cd4_ssub_dge_mark <- list("EF_Low-AB_Ctrl" = cd4_ssub_mark_low, "CD_High-AB_Ctrl" = cd4_ssub_mark_high, "CD_High-EF_Low" = cd4_ssub_mark_highlow) cd4_ssub@misc[["dgemarkers"]] <- cd4_ssub_dge_mark for (contr in names(cd4_ssub@misc$dgemarkers)) { deg_res <- cd4_ssub@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 = exah.gmt, nPerm = 10000) if (nrow(contr.enricher.gsea) > 0) { # Write res contr.enricher.gsea.res <- list(contr.enricher.gsea) names(contr.enricher.gsea.res) <- c(contr) write.xlsx(contr.enricher.gsea.res, paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", contr, "Cd4_Exhaustion_GSEA.xlsx", sep = "/")) # Plot for (t in seq(1, nrow(contr.enricher.gsea))) { d <- contr.enricher.gsea$Description[t] p <- gseaplot2(contr.enricher.gsea, geneSetID = t, #title = gsub("_", " " , d), color = c("forestgreen"), pvalue_table = T) fn <- paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", contr, paste0("Cd4_", d, "_GSEA.pdf"), sep = "/") ggsave(filename = fn, plot = p, width = 10, height = 7) } } } # for (contr in names(cd4_ssub@misc$dgemarkers)) { # deg_res <- cd4_ssub@misc$dgemarkers[[contr]] # print(nrow(deg_res)) # deg_res_top <- subset(deg_res, p_val_adj < 1e-06 & avg_logFC > 0) # cat(cdsub, contr, nrow(deg_res_top), "\n") # # Enrichment # contr.enrichement <- enricher(rownames(deg_res_top), TERM2GENE = exah.gmt, pvalueCutoff = 0.05, universe = rownames(cd4_ssub)) # #print(contr.enrichement) # if (nrow(contr.enrichement) > 0) { # print(contr.enrichement@result[,1:7]) # # # Write res # # contr.enricher.gsea.res <- list(contr.enricher.gsea) # # names(contr.enricher.gsea.res) <- c(contr) # # write.xlsx(contr.enricher.gsea.res, paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", cdsub, contr, "Exhaustion_GSEA.xlsx", sep = "/")) # # # Plot # # for (t in seq(1, nrow(contr.enricher.gsea))) { # # d <- contr.enricher.gsea$Description[t] # # p <- gseaplot2(contr.enricher.gsea, geneSetID = t, # # #title = gsub("_", " " , d), # # color = c("forestgreen"), # # pvalue_table = T) # # fn <- paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", cdsub, contr, paste0(d, "_GSEA.pdf"), sep = "/") # # ggsave(filename = fn, plot = p, width = 10, height = 7) # # } # } # } saveRDS(object = cd4_ssub, file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "TcellFiltered_fastMNN_Cd4subTest_final.rds", sep = "/")) cd8_ssub <- subset(x = fmnn_tcell, Cd8a > Cd4 & Cd8a > 0) dim(cd8_ssub) FeaturePlot(object = cd8_ssub, features = c("Cd4", "Cd8a"), reduction = "umap", split.by = "RNA_Group") Idents(cd8_ssub) <- "RNA_Group" cd8_ssub_mark_low <- FindMarkers(object = cd8_ssub, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) cd8_ssub_mark_high <- FindMarkers(object = cd8_ssub, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) cd8_ssub_mark_highlow <- FindMarkers(object = cd8_ssub, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) cd8_ssub_dge_mark <- list("EF_Low-AB_Ctrl" = cd8_ssub_mark_low, "CD_High-AB_Ctrl" = cd8_ssub_mark_high, "CD_High-EF_Low" = cd8_ssub_mark_highlow) cd8_ssub@misc[["dgemarkers"]] <- cd8_ssub_dge_mark for (contr in names(cd8_ssub@misc$dgemarkers)) { deg_res <- cd8_ssub@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 = exah.gmt, nPerm = 10000) if (nrow(contr.enricher.gsea) > 0) { # Write res contr.enricher.gsea.res <- list(contr.enricher.gsea) names(contr.enricher.gsea.res) <- c(contr) write.xlsx(contr.enricher.gsea.res, paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", contr, "Cd8_Exhaustion_GSEA.xlsx", sep = "/")) # Plot for (t in seq(1, nrow(contr.enricher.gsea))) { d <- contr.enricher.gsea$Description[t] p <- gseaplot2(contr.enricher.gsea, geneSetID = t, #title = gsub("_", " " , d), color = c("forestgreen"), pvalue_table = T) fn <- paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", contr, paste0("Cd8_", d, "_GSEA.pdf"), sep = "/") ggsave(filename = fn, plot = p, width = 10, height = 7) } } } # for (contr in names(cd8_ssub@misc$dgemarkers)) { # deg_res <- cd8_ssub@misc$dgemarkers[[contr]] # print(nrow(deg_res)) # deg_res_top <- subset(deg_res, p_val_adj < 1e-06 & avg_logFC > 0) # cat(cdsub, contr, nrow(deg_res_top), "\n") # # Enrichment # contr.enrichement <- enricher(rownames(deg_res_top), TERM2GENE = exah.gmt, pvalueCutoff = 0.05, universe = rownames(cd8_ssub)) # #print(contr.enrichement) # if (nrow(contr.enrichement) > 0) { # print(contr.enrichement@result[,1:7]) # # # Write res # # contr.enricher.gsea.res <- list(contr.enricher.gsea) # # names(contr.enricher.gsea.res) <- c(contr) # # write.xlsx(contr.enricher.gsea.res, paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", cdsub, contr, "Exhaustion_GSEA.xlsx", sep = "/")) # # # Plot # # for (t in seq(1, nrow(contr.enricher.gsea))) { # # d <- contr.enricher.gsea$Description[t] # # p <- gseaplot2(contr.enricher.gsea, geneSetID = t, # # #title = gsub("_", " " , d), # # color = c("forestgreen"), # # pvalue_table = T) # # fn <- paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", cdsub, contr, paste0(d, "_GSEA.pdf"), sep = "/") # # ggsave(filename = fn, plot = p, width = 10, height = 7) # # } # } # } saveRDS(object = cd8_ssub, file = paste(wdir, "TcellFiltered", "02-fastMNN_2-post-analysis", "TcellFiltered_fastMNN_Cd8subTest_final.rds", sep = "/")) # Monocytes fmnn_monoc <- readRDS(paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_fastMNN_final.rds", sep = "/")) Idents(fmnn_monoc) <- "RNA_Group" fmnn_monoc_mark_low <- FindMarkers(object = fmnn_monoc, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_monoc_mark_high <- FindMarkers(object = fmnn_monoc, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_monoc_mark_highlow <- FindMarkers(object = fmnn_monoc, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) fmnn_monoc_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_monoc_mark_low, "CD_High-AB_Ctrl" = fmnn_monoc_mark_high, "CD_High-EF_Low" = fmnn_monoc_mark_highlow) fmnn_monoc@misc[["dgemarkers"]] <- fmnn_monoc_dge_mark write.xlsx(x = fmnn_monoc_dge_mark, file = paste(wdir, "MonocytesFiltered", "02-fastMNN_2-post-analysis", "DGE_markers.xlsx", sep = "/"), row.names = TRUE) saveRDS(object = fmnn_monoc, file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_fastMNN_final.rds", sep = "/")) for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) { res_dge <- fmnn_monoc_dge_mark[[contr]] print(head(res_dge)) vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr) ggsave(filename = paste(wdir, "MonocytesFiltered", "02-fastMNN_2-post-analysis", paste(contr, "volcano.png", sep = "-"), sep = "/"), plot = vol, width = 10, height = 9, dpi = 600) } # Monocytes No3 fmnn_monoc <- readRDS(paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_No3_fastMNN_final.rds", sep = "/")) Idents(fmnn_monoc) <- "RNA_Group" fmnn_monoc_mark_low <- FindMarkers(object = fmnn_monoc, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_monoc_mark_high <- FindMarkers(object = fmnn_monoc, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_monoc_mark_highlow <- FindMarkers(object = fmnn_monoc, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) fmnn_monoc_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_monoc_mark_low, "CD_High-AB_Ctrl" = fmnn_monoc_mark_high, "CD_High-EF_Low" = fmnn_monoc_mark_highlow) fmnn_monoc@misc[["dgemarkers"]] <- fmnn_monoc_dge_mark write.xlsx(x = fmnn_monoc_dge_mark, file = paste(wdir, "MonocytesFiltered", "02-fastMNN_2-post-analysis", "DGE_markers_No3.xlsx", sep = "/"), row.names = TRUE) saveRDS(object = fmnn_monoc, file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_No3_fastMNN_final.rds", sep = "/")) for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) { res_dge <- fmnn_monoc_dge_mark[[contr]] print(head(res_dge)) vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr) ggsave(filename = paste(wdir, "MonocytesFiltered", "02-fastMNN_2-post-analysis", paste(contr, "_No3_volcano.png", sep = "-"), sep = "/"), plot = vol, width = 10, height = 9, dpi = 600) } # Full cluster markers fmnn_full <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/")) Idents(fmnn_full) <- "RNA_snn_res.0.6" clu_deg <- list() for (clu in unique(fmnn_full$RNA_snn_res.0.6)) { print(clu) clu_mark <- list() fmnn_full_clu <- subset(fmnn_full, ident = clu) Idents(fmnn_full_clu) <- "RNA_Group" clu_mark[["EF_Low-AB_Ctrl"]] <- FindMarkers(object = fmnn_full_clu, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) clu_mark[["CD_High-AB_Ctrl"]] <- FindMarkers(object = fmnn_full_clu, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) clu_mark[["CD_High-EF_Low"]] <- FindMarkers(object = fmnn_full_clu, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) clu_deg[[clu]] <- clu_mark } fmnn_full@misc[["deg_clumarkers"]] <- clu_deg saveRDS(object = fmnn_full, file = paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/")) wb <- createWorkbook() for (clu in names(clu_deg)) { addWorksheet(wb, clu) writeData(wb, clu, "EF_Low-AB_Ctrl", startCol = 1, startRow = 1) writeData(wb, clu, clu_deg[[clu]][["EF_Low-AB_Ctrl"]], startCol = 1, startRow = 2, rowNames = TRUE) writeData(wb, clu, "CD_High-AB_Ctrl", startCol = 10, startRow = 1) writeData(wb, clu, clu_deg[[clu]][["CD_High-AB_Ctrl"]], startCol = 10, startRow = 2, rowNames = TRUE) } ## Save workbook to working directory saveWorkbook(wb, paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_deg_clumarkers.xlsx", sep = "/"), overwrite = TRUE) # SingleR Annotation fmnn_full <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/")) organism.db <- list(ImmGenData = ImmGenData(), MouseRNAseqData = MouseRNAseqData()) for(el in names(organism.db)){ # Main lname <- paste("SingleR", el, sep = "_") fmnn_full@misc[[lname]] <- SingleR(test = fmnn_full@assays$RNA@data, ref = organism.db[[el]], labels = organism.db[[el]]$label.main) lname_lab <- paste(lname, "labels", sep = "_") fmnn_full@meta.data[[lname_lab]] <- fmnn_full@misc[[lname]]$labels # png(filename = paste(plot_dir, paste(sample, lname, "umap", "plot.png", sep = "_"), sep = "/"), # units = "in", width = 12, height = 9, res = 200) # print(DimPlot(object = fmnn_full, reduction = "umap", group.by = lname_lab, label = T)) # dev.off() # Refined lname_ref <- paste("SingleRrefined", el, sep = "_") fmnn_full@misc[[lname_ref]] <- SingleR(test = fmnn_full@assays$RNA@data, ref = organism.db[[el]], labels = organism.db[[el]]$label.fine) lname_ref_lab <- paste(lname_ref, "labels", sep = "_") fmnn_full@meta.data[[lname_ref_lab]] <- fmnn_full@misc[[lname_ref]]$labels # png(filename = paste(plot_dir, paste(sample, lname_ref, "umap", "plot.png", sep = "_"), sep = "/"), # units = "in", width = 12, height = 9, res = 200) # print(DimPlot(object = fmnn_full, reduction = "umap", group.by = lname_ref_lab, label = T)) # dev.off() } saveRDS(object = fmnn_full, file = paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/")) # B cells fmnn_bcells <- readRDS(file = paste(wdir, "BcellFiltered", "01-fastMNN_2", "BcellFiltered_fastMNN_final.rds", sep = "/")) Idents(fmnn_bcells) <- "RNA_Group" fmnn_bcells_mark_low <- FindMarkers(object = fmnn_bcells, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_bcells_mark_high <- FindMarkers(object = fmnn_bcells, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) fmnn_bcells_mark_highlow <- FindMarkers(object = fmnn_bcells, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) fmnn_bcells_dge_mark <- list("EF_Low-AB_Ctrl" = fmnn_bcells_mark_low, "CD_High-AB_Ctrl" = fmnn_bcells_mark_high, "CD_High-EF_Low" = fmnn_bcells_mark_highlow) fmnn_bcells@misc[["dgemarkers"]] <- fmnn_bcells_dge_mark write.xlsx(x = fmnn_bcells_dge_mark, file = paste(wdir, "BcellFiltered", "02-fastMNN_2-post-analysis", "DGE_markers_No3.xlsx", sep = "/"), row.names = TRUE) saveRDS(object = fmnn_bcells, file = paste(wdir, "BcellFiltered", "01-fastMNN_2", "BcellFiltered_fastMNN_final.rds", sep = "/")) for (contr in c("EF_Low-AB_Ctrl" , "CD_High-AB_Ctrl", "CD_High-EF_Low")) { res_dge <- fmnn_bcells_dge_mark[[contr]] print(head(res_dge)) vol <- plotVolcano(res_dge, adj.pval.thr, logfc.thr) ggsave(filename = paste(wdir, "BcellFiltered", "02-fastMNN_2-post-analysis", paste(contr, "_volcano.png", sep = "-"), sep = "/"), plot = vol, width = 10, height = 9, dpi = 600) } # # Macrophages cluster markers # fmnn_macro <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) # Idents(fmnn_macro) <- "RNA_snn_res.0.6" # clu_deg <- list() # for (clu in unique(fmnn_macro$RNA_snn_res.0.6)) { # print(clu) # clu_mark <- list() # fmnn_macro_clu <- subset(fmnn_macro, ident = clu) # Idents(fmnn_macro_clu) <- "RNA_Group" # clu_mark[["EF_Low-AB_Ctrl"]] <- FindMarkers(object = fmnn_macro_clu, ident.1 = "EF_Low", ident.2 = "AB_Ctrl", logfc.threshold = 0) # clu_mark[["CD_High-AB_Ctrl"]] <- FindMarkers(object = fmnn_macro_clu, ident.1 = "CD_High", ident.2 = "AB_Ctrl", logfc.threshold = 0) # clu_mark[["CD_High-EF_Low"]] <- FindMarkers(object = fmnn_macro_clu, ident.1 = "CD_High", ident.2 = "EF_Low", logfc.threshold = 0) # clu_deg[[clu]] <- clu_mark # } # fmnn_macro@misc[["deg_clumarkers"]] <- clu_deg # saveRDS(object = fmnn_macro, # file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/")) # # full_res <- data.frame() # 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) # res_filt <- res[, "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 # # wb <- createWorkbook() # for (clu in names(clu_deg)) { # addWorksheet(wb, clu) # writeData(wb, clu, "EF_Low-AB_Ctrl", startCol = 1, startRow = 1) # writeData(wb, clu, clu_deg[[clu]][["EF_Low-AB_Ctrl"]], startCol = 1, startRow = 2, rowNames = TRUE) # writeData(wb, clu, "CD_High-AB_Ctrl", startCol = 10, startRow = 1) # writeData(wb, clu, clu_deg[[clu]][["CD_High-AB_Ctrl"]], startCol = 10, startRow = 2, rowNames = TRUE) # } # ## Save workbook to working directory # saveWorkbook(wb, paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_deg_clumarkers.xlsx", sep = "/"), overwrite = TRUE) #