```{r} library(Seurat) library(magrittr) library(tidyverse) library(EnhancedVolcano) library(msigdbr) library(fgsea) library(scales) library(cowplot) library(lsr) library(stringi) source("/Users/ctzouanas/Documents/MIT/Shalek/AncestralWisdom/fgsea_wrappers_CT.R") source("/Users/ctzouanas/Documents/MIT/Shalek/AncestralWisdom/AncestralWisdom_Terra/speciesconversion_biomart.R") library(biomaRt) library(ActivePathways) library(DirichletReg) library(monocle) library(SeuratObject) library(SeuratWrappers) library(SeuratDisk) library(qusage) library(pheatmap) ancestral.wisdom.path <- "/Users/ctzouanas/Documents/MIT/Shalek/AncestralWisdom/AncestralWisdom_Terra" code.files <- list.files(paste0(ancestral.wisdom.path))[grepl(pattern = ".R$", x = list.files(paste0(ancestral.wisdom.path)))] files.to.sources <- paste0(ancestral.wisdom.path, "/", code.files) for(ii in 1:length(files.to.sources)){ file.ii <- files.to.sources[[ii]] source(file.ii) } ``` ```{r Load necessary objects} BALrog.obj <- readRDS("/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/BAL_velo_srt_LOCKED.RDS") PBj.obj <- readRDS("/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/PBMC_srt_SEMI_LOCKED.rds") BAL_PB_myeloid_srt <- readRDS("/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/BAL_PBMC_myeloid_srt_LOCKED.rds") BAL_PB_full <- readRDS("/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Kwon_2020_Manuscript_rdata_objects_BAL_PBMC_integrated_LOCKED_srt.rds") ragon.id.mapping <- data.frame( ragon.id = c("851204", "936223", "217405", "608217", "911594", "445474", "899898", "343744", "733361"), pack.years = c(0, 0, 0, 0, 17, 28, 17, 6, 15), age = c(39, 32, 23, 22, 52, 55, 64, 57, 60) ) if(exists("BAL_PB_myeloid_srt")){ BAL_PB_myeloid_srt@meta.data$Smoking_Status[BAL_PB_myeloid_srt@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0118"))] = 'SMOKER' BAL_PB_myeloid_srt@meta.data$Smoking_Status[BAL_PB_myeloid_srt@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0208"))] = 'SMOKER' BAL_PB_myeloid_srt@meta.data$Smoking_Status[BAL_PB_myeloid_srt@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0425"))] = 'SMOKER' BAL_PB_myeloid_srt@meta.data$Smoking_Status[BAL_PB_myeloid_srt@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0513"))] = 'NON_SMOKER' BAL_PB_myeloid_srt@meta.data$Smoking_Status[BAL_PB_myeloid_srt@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0531"))] = 'SMOKER' BAL_PB_myeloid_srt@meta.data$Smoking_Status[BAL_PB_myeloid_srt@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "07XX"))] = 'SMOKER' BAL_PB_myeloid_srt@meta.data$Smoking_Status[BAL_PB_myeloid_srt@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0911"))] = 'NON_SMOKER' BAL_PB_myeloid_srt@meta.data$Smoking_Status[BAL_PB_myeloid_srt@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0924"))] = 'NON_SMOKER' BAL_PB_myeloid_srt@meta.data$Smoking_Status[BAL_PB_myeloid_srt@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "1130"))] = 'NON_SMOKER' BAL_PB_myeloid_srt@meta.data$Smoking_Status = as.factor(BAL_PB_myeloid_srt@meta.data$Smoking_Status) BAL_PB_myeloid_srt@meta.data$Smoking_Status = as.character(BAL_PB_myeloid_srt@meta.data$Smoking_Status) table(BAL_PB_myeloid_srt$Smoking_Status, BAL_PB_myeloid_srt$compartment) table(BAL_PB_myeloid_srt$Smoking_Status, BAL_PB_myeloid_srt$orig.ident) volunteer.id <- gsub(pattern = "BAL|PBMC", replacement = "volunteer_", x = BAL_PB_myeloid_srt$orig.ident) volunteer.id <- factor(x = volunteer.id, levels = c("volunteer_0118", "volunteer_0208", "volunteer_0425", "volunteer_0531", "volunteer_07XX", "volunteer_0513", "volunteer_0911", "volunteer_0924", "volunteer_1130"), ordered = TRUE) BAL_PB_myeloid_srt$volunteer.id <- volunteer.id } if(exists("BALrog.obj")){ BALrog.obj$pack.years <- factor( BALrog.obj$Ragon_ID, levels = ragon.id.mapping$ragon.id, labels = ragon.id.mapping$pack.years ) %>% as.character() %>% as.numeric() BALrog.obj@meta.data$Smoking_Status[BALrog.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0118"))] = 'SMOKER' BALrog.obj@meta.data$Smoking_Status[BALrog.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0208"))] = 'SMOKER' BALrog.obj@meta.data$Smoking_Status[BALrog.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0425"))] = 'SMOKER' BALrog.obj@meta.data$Smoking_Status[BALrog.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0513"))] = 'NON_SMOKER' BALrog.obj@meta.data$Smoking_Status[BALrog.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0531"))] = 'SMOKER' BALrog.obj@meta.data$Smoking_Status[BALrog.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "07XX"))] = 'SMOKER' BALrog.obj@meta.data$Smoking_Status[BALrog.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0911"))] = 'NON_SMOKER' BALrog.obj@meta.data$Smoking_Status[BALrog.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0924"))] = 'NON_SMOKER' BALrog.obj@meta.data$Smoking_Status[BALrog.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "1130"))] = 'NON_SMOKER' BALrog.obj@meta.data$Smoking_Status = as.factor(BALrog.obj@meta.data$Smoking_Status) BALrog.obj@meta.data$Smoking_Status = as.character(BALrog.obj@meta.data$Smoking_Status) table(BALrog.obj$Smoking_Status, BALrog.obj$orig.ident) volunteer.id <- gsub(pattern = "BAL|PBMC", replacement = "volunteer_", x = BALrog.obj$orig.ident) volunteer.id <- factor(x = volunteer.id, levels = c("volunteer_0118", "volunteer_0208", "volunteer_0425", "volunteer_0531", "volunteer_07XX", "volunteer_0513", "volunteer_0911", "volunteer_0924", "volunteer_1130"), ordered = TRUE) BALrog.obj$volunteer.id <- volunteer.id } if(exists("PBj.obj")){ PBj.obj@meta.data$Smoking_Status[PBj.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0118"))] = 'SMOKER' PBj.obj@meta.data$Smoking_Status[PBj.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0208"))] = 'SMOKER' PBj.obj@meta.data$Smoking_Status[PBj.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0425"))] = 'SMOKER' PBj.obj@meta.data$Smoking_Status[PBj.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0513"))] = 'NON_SMOKER' PBj.obj@meta.data$Smoking_Status[PBj.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0531"))] = 'SMOKER' PBj.obj@meta.data$Smoking_Status[PBj.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "07XX"))] = 'SMOKER' PBj.obj@meta.data$Smoking_Status[PBj.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0911"))] = 'NON_SMOKER' PBj.obj@meta.data$Smoking_Status[PBj.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "0924"))] = 'NON_SMOKER' PBj.obj@meta.data$Smoking_Status[PBj.obj@meta.data$orig.ident %in% (paste0(c("BAL", "PBMC"), "1130"))] = 'NON_SMOKER' PBj.obj@meta.data$Smoking_Status = as.factor(PBj.obj@meta.data$Smoking_Status) PBj.obj@meta.data$Smoking_Status = as.character(PBj.obj@meta.data$Smoking_Status) table(PBj.obj$Smoking_Status, PBj.obj$orig.ident) volunteer.id <- gsub(pattern = "BAL|PBMC", replacement = "volunteer_", x = PBj.obj$orig.ident) volunteer.id <- factor(x = volunteer.id, levels = c("volunteer_0118", "volunteer_0208", "volunteer_0425", "volunteer_0531", "volunteer_07XX", "volunteer_0513", "volunteer_0911", "volunteer_0924", "volunteer_1130"), ordered = TRUE) PBj.obj$volunteer.id <- volunteer.id } figdir.save <- "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Clean_Figdir" ``` ```{r Read in gene lists to test} presentation.substring = "ANTIGEN|PRESENTATION" phagocytosis.substring = "PHAGOCYTOSIS|PHAGOCYTIC|PHAGOLYSOSOME|PHAGOSOMAL" myeloid.substring = "MYELOID|MONOCYTE|MACROPHAGE" differentiation.substring = "MATURATION|DIFFERENTIATION" presentation.sets <- geneset.loader(gsea.category.in = NA, gsea.subcategory.in = NA, species.in = "Homo sapiens", oneoff.folder.path.in = NULL, substring.list.to.filter = list(presentation.substring)) phagocytosis.sets <- geneset.loader(gsea.category.in = NA, gsea.subcategory.in = NA, species.in = "Homo sapiens", oneoff.folder.path.in = NULL, substring.list.to.filter = list(phagocytosis.substring)) myeloid.diff.sets <- geneset.loader(gsea.category.in = NA, gsea.subcategory.in = NA, species.in = "Homo sapiens", oneoff.folder.path.in = NULL, substring.list.to.filter = list(myeloid.substring, differentiation.substring)) presentation.lists.to.exclude <- "T_CELL|BCELL|BREECH|HDAC|B_CELL|NMDA|TESTIS|TCELL" differentiation.lists.to.exclude <- "MICROGLIA" presentation.sets <- presentation.sets[!grepl(presentation.lists.to.exclude, names(presentation.sets))] myeloid.diff.sets <- myeloid.diff.sets[!grepl(differentiation.lists.to.exclude, names(myeloid.diff.sets))] chosen.patterns <- c("CYTOKINE", "CHEMOKINE", "IMMUNE", "MYELOID", "INFLAMMATION", "INFLAMMATORY", "LEUKOCYTE", "INTERLEUKIN", "IMMUNITY", "IMMUNE") immune.patterns <- paste(chosen.patterns, sep = "|", collapse = "|") immune.sets <- geneset.loader(species.in = "Homo sapiens", oneoff.folder.path.in = NULL, substring.list.to.filter = list(immune.patterns)) scale.fill.limits = c(-2.25, 2.25) ``` ```{r Compositional analyses on BAL cells} BALrog.obj.composition <- BALrog.obj@meta.data[, c("volunteer.id", "Smoking_Status", "celltypes_2")] %>% group_by(volunteer.id, Smoking_Status, celltypes_2) %>% dplyr::summarise(count = n()) %>% ungroup() %>% group_by(volunteer.id, Smoking_Status) %>% dplyr::mutate(Freq = count / sum(count)) %>% ungroup() %>% dplyr::mutate(Smoking_Status = if_else(Smoking_Status == "SMOKER", 1, 0)) %>% dplyr::select(!c(count)) %>% tidyr::pivot_wider(names_from = celltypes_2, values_from = Freq, values_fill = 0) %>% column_to_rownames(var = "volunteer.id") %>% ungroup() smoking.status.vec <- BALrog.obj.composition[, 1] %>% as.logical() wilcox.df <- data.frame() for(ii in 2:ncol(BALrog.obj.composition)){ celltype.ii <- colnames(BALrog.obj.composition)[[ii]] smoking.vec <- BALrog.obj.composition[smoking.status.vec, ii] nonsmoking.vec <- BALrog.obj.composition[!smoking.status.vec, ii] pval.ii <- wilcox.test(x = smoking.vec, y = nonsmoking.vec, exact = TRUE)$p.value wilcox.df.ii <- data.frame(celltype = celltype.ii, p.val = pval.ii) wilcox.df <- rbind(wilcox.df, wilcox.df.ii) } wilcox.df$p.val.adj <- p.adjust(p = wilcox.df$p.val, method = "BH") AL <- DR_data(BALrog.obj.composition[,colnames(BALrog.obj.composition) != "Smoking_Status"]) test2 <- DirichReg(AL ~ Smoking_Status, BALrog.obj.composition) summary.out <- summary(test2) df.out <- summary.out[["coef.mat"]] %>% as.data.frame() %>% rownames_to_column("variable.name") %>% dplyr::filter(grepl(pattern = "Smoking_Status", x = .$variable.name)) %>% dplyr::mutate(cluster = paste0(summary.out[["varnames"]])) %>% dplyr::select(cluster, Estimate, `Pr(>|z|)`) %>% dplyr::rename("dirich.pval" = 3, "dirich.Estimate" = 2) %>% arrange(dirich.pval) composition.plot.df <- BALrog.obj@meta.data[, c("volunteer.id", "Smoking_Status", "celltypes_2")] %>% group_by(volunteer.id, Smoking_Status, celltypes_2) %>% dplyr::summarise(count = n()) %>% ungroup() %>% group_by(volunteer.id, Smoking_Status) %>% dplyr::mutate(Freq = count / sum(count)) %>% ungroup() %>% tidyr::complete(volunteer.id, celltypes_2) %>% dplyr::mutate(count = if_else(is.na(count), 0, as.numeric(count))) %>% dplyr::mutate(Freq = if_else(is.na(Freq), 0, as.numeric(Freq))) %>% dplyr::mutate(Smoking_Status = if_else(volunteer.id %in% c(paste0("volunteer_", c("0118", "0208", "0425", "0531", "07XX"))), "SMOKER", "NON_SMOKER")) %>% inner_join(x = ., y = df.out, by = c("celltypes_2" = "cluster")) %>% dplyr::mutate(celltypes_2 = gsub(pattern = "_", replacement = " ", x = celltypes_2)) %>% dplyr::mutate(plot.title = paste0(celltypes_2)) %>% dplyr::mutate(Smoking_Status = factor(Smoking_Status, levels = c("NON_SMOKER", "SMOKER"), labels = c("Non\nSmoker", "Smoker"))) %>% dplyr::mutate(volunteer.id = as.character(volunteer.id)) font.size = 16 set.seed(42) p.dirich <- ggplot(composition.plot.df, aes(x = Smoking_Status, y = Freq)) + geom_boxplot(outlier.shape = NA) + geom_point(aes(color = volunteer.id), size = 3, position = position_jitter(w = 0.2, h = 0)) + ylab("Cell Type Proportion") + theme_classic() + facet_wrap(~ plot.title, scales = "free", nrow = 2) + theme(legend.position = "none", axis.title.x = element_blank(), strip.background = element_blank(), text = element_text(size = font.size), strip.text = element_text(size = font.size)) pdf(file = paste0(figdir.save, "/", "BALComposition_DirichletRegression_220813.pdf"), height = 12, width = 12) plot(p.dirich) dev.off() bothtests.pval.df <- inner_join( x = df.out, y = wilcox.df %>% dplyr::select(celltype, p.val.adj) %>% dplyr::rename(wilcox.p.val.adj = 2), by = c("cluster" = "celltype") ) CD93.expr <- BALrog.obj@assays$RNA@counts["CD93", ] > 0 BALrog.obj$CD93.Status <- if_else(CD93.expr, "CD93pos", "CD93neg") BALrog.obj$CD93.Status <- factor(BALrog.obj$CD93.Status, levels = c("CD93pos", "CD93neg"), labels = c("CD93: +", "CD93: -")) Idents(BALrog.obj) <- "celltypes_2" CD93.BALmono.marker <- FindMarkers(BALrog.obj, ident.1 = "CD93: +", group.by = "CD93.Status", subset.ident = "BAL_Monocytes", features = c("FCGR3A", "CD14"), logfc.threshold = 0) %>% rownames_to_column(var = "gene") %>% arrange(-avg_log2FC) p.balmono.cd14 <- VlnPlot(BALrog.obj[, BALrog.obj$celltypes_2 == "BAL_Monocytes"], features = c("CD14"), pt.size = 0, group.by = "CD93.Status") + ggtitle("Alveolar Monocytes:\nCD14 Expression") + ylab("CD14 Expression Level") + theme(legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 0, hjust = 0.5)) pdf(file = paste0(figdir.save, "/", "BALMonocyte_CD14Expression_by_CD93Status.pdf"), height = 3, width = 2) plot(p.balmono.cd14) dev.off() ``` ```{r Make supplementary QC plots} p1 <- VlnPlot(BALrog.obj, features = "nFeature_RNA", group.by = "celltypes_2", pt.size = 0) + labs(title = "Number of Genes Detected") + NoLegend() + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + scale_y_continuous(trans = "log10") p2 <- VlnPlot(BALrog.obj, features = "nCount_RNA", group.by = "celltypes_2", pt.size = 0) + labs(title = "Number of UMIs Detected") + NoLegend() + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + scale_y_continuous(trans = "log10") p3 <- VlnPlot(BALrog.obj, features = "percent.mt", group.by = "celltypes_2", pt.size = 0) + labs(title = "Percent of Mitochondrial Reads") + NoLegend() + theme(axis.title.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) plotlist.cluster <- list(p1, p2, p3) p4 <- cowplot::plot_grid(plotlist = plotlist.cluster, nrow = length(plotlist.cluster), align = "v", rel_heights = c(1, 1, 1.5)) # plot(p4) p5 <- VlnPlot(BALrog.obj, features = "nFeature_RNA", group.by = "Smoking_Status", pt.size = 0) + labs(title = "") + NoLegend() + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + scale_y_continuous(trans = "log10") p6 <- VlnPlot(BALrog.obj, features = "nCount_RNA", group.by = "Smoking_Status", pt.size = 0) + labs(title = "") + NoLegend() + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + scale_y_continuous(trans = "log10") p7 <- VlnPlot(BALrog.obj, features = "percent.mt", group.by = "Smoking_Status", pt.size = 0) + labs(title = "") + NoLegend() + theme(axis.title.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) plotlist.cluster <- list(p5, p6, p7) p8 <- cowplot::plot_grid(plotlist = plotlist.cluster, nrow = length(plotlist.cluster), align = "v", rel_heights = c(1, 1, 1.5)) # plot(p8) p9 <- VlnPlot(BAL_PB_myeloid_srt, features = "nFeature_RNA", group.by = "compartment", pt.size = 0) + labs(title = "") + NoLegend() + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + scale_y_continuous(trans = "log10") p10 <- VlnPlot(BAL_PB_myeloid_srt, features = "nCount_RNA", group.by = "compartment", pt.size = 0) + labs(title = "") + NoLegend() + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + scale_y_continuous(trans = "log10") p11 <- VlnPlot(BAL_PB_myeloid_srt, features = "percent.mt", group.by = "compartment", pt.size = 0) + labs(title = "") + NoLegend() + theme(axis.title.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) plotlist.cluster <- list(p9, p10, p11) p12 <- cowplot::plot_grid(plotlist = plotlist.cluster, nrow = length(plotlist.cluster), align = "v", rel_heights = c(1, 1, 1.5)) # plot(p12) p17 <- VlnPlot(BALrog.obj, features = "nFeature_RNA", group.by = "volunteer.id", pt.size = 0) + labs(title = "") + NoLegend() + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + scale_y_continuous(trans = "log10") p18 <- VlnPlot(BALrog.obj, features = "nCount_RNA", group.by = "volunteer.id", pt.size = 0) + labs(title = "") + NoLegend() + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + scale_y_continuous(trans = "log10") p19 <- VlnPlot(BALrog.obj, features = "percent.mt", group.by = "volunteer.id", pt.size = 0) + labs(title = "") + NoLegend() + theme(axis.title.x=element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) plotlist.cluster <- list(p17, p18, p19) p20 <- cowplot::plot_grid(plotlist = plotlist.cluster, nrow = length(plotlist.cluster), align = "v", rel_heights = c(1, 1, 1.5)) # plot(p20) widths <- c(2.5, 1, 1, 1.75) heights = c(1, 1, 1) big.plotlist <- list(p1, p5, p9, p17, p2, p6, p10, p18, p3, p7, p11, p19) p21 <- cowplot::plot_grid(plotlist = big.plotlist, ncol = 4, nrow = 3, align = "hv", rel_widths = widths, rel_heights = heights) plot(p21) pdf(file = paste0(figdir.save, "/QC by cluster volunteer compartment smoking.pdf"), width = 16, height = 18, useDingbats = FALSE) plot(p21) dev.off() p22 <- ggplot(BALrog.obj@meta.data, aes(x=BALrog.obj$volunteer.id, fill=eval(parse(text = "celltypes_2")))) + geom_bar(position = "fill", colour="black") + theme_classic() + NoLegend() + xlab("") + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + ylab("Cell Type Proportion") BALrog.obj$celltypes_4 <- factor(x = BALrog.obj$celltypes_2, levels = sort(unique(BALrog.obj$celltypes_2))) p23 <- ggplot(BALrog.obj@meta.data, aes(x=BALrog.obj$celltypes_4, fill=eval(parse(text = "Smoking_Status")))) + geom_bar(position = "fill", colour="black") + theme_classic() + NoLegend() + xlab("") + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + ylab("Cell Type Proportion") p24 <- cowplot::plot_grid(plotlist = list(p22, p.dirich + theme(text = element_text(size = 12), strip.text.x = element_text(size = 12))), ncol = 2, nrow = 1, align = "hv", axis = "tblr") pdf(file = paste0(figdir.save, "/QC cell type composition by donor and smoking 220817.pdf"), width = 16, height = 6, useDingbats = FALSE) plot(p24) dev.off() ``` ```{r Look at expression of CCR2/3/5} monocyte.obj <- BAL_PB_myeloid_srt[,BAL_PB_myeloid_srt$celltypes_3 %in% c("CD14_Monocytes", "CD16_Monocytes", "BAL_Monocytes")] receptors <- c("CCR2", "CCR3", "CCR5") monocyte.obj$monocyte.type <- monocyte.obj$celltypes_3 monocyte.obj$monocyte.type[monocyte.obj$monocyte.type %in% c("CD14_Monocytes", "CD16_Monocytes")] <- "Blood_Monocytes" monocyte.obj$monocyte.type <- factor(monocyte.obj$monocyte.type, levels = c("Blood_Monocytes", "BAL_Monocytes")) VlnPlot(monocyte.obj, features = receptors, group.by = "monocyte.type", split.by = "Smoking_Status", split.plot = T, ncol = 1) celltypes <- unique(monocyte.obj$monocyte.type) wilcox.df <- data.frame() for(ii in 1:length(celltypes)){ smoker.cell.names <- monocyte.obj$Smoking_Status == "SMOKER" & monocyte.obj$monocyte.type == celltypes[[ii]] nonsmoker.cell.names <- monocyte.obj$Smoking_Status == "NON_SMOKER" & monocyte.obj$monocyte.type == celltypes[[ii]] for(jj in 1:length(receptors)){ smoker.gene.vals <- monocyte.obj[receptors[jj], smoker.cell.names]@assays$RNA@data[1,] nonsmoker.gene.vals <- monocyte.obj[receptors[jj], nonsmoker.cell.names]@assays$RNA@data[1,] wilcox.result.ii <- wilcox.test(smoker.gene.vals, nonsmoker.gene.vals) wilcox.df.ii <- data.frame(gene = receptors[jj], celltype = celltypes[[ii]], pval = wilcox.result.ii$p.value, median.smoker = median(smoker.gene.vals), median.nonsmoker = median(nonsmoker.gene.vals), cohens.D = cohensD(x = smoker.gene.vals, y = nonsmoker.gene.vals)) wilcox.df <- rbind(wilcox.df, wilcox.df.ii) } } DefaultAssay(monocyte.obj) wilcox.df$pval.benhoch <- p.adjust(p = wilcox.df$pval, method = "BH") vln.collection <- VlnPlot(monocyte.obj, features = receptors, group.by = "monocyte.type", split.by = "Smoking_Status", split.plot = TRUE, pt.size = 0.5, combine = FALSE) for(ii in 1:length(vln.collection)){ p.ii <- vln.collection[[ii]] p.ii.updated <- p.ii + NoLegend() + theme(axis.title.x = element_blank(), axis.text.x = element_text(angle = 0, hjust = 0.5)) + ylab("log(TPK10 + 1)") plot(p.ii.updated) vln.collection[[ii]] <- p.ii.updated } p.merge <- cowplot::plot_grid(plotlist = vln.collection, nrow = 1, align = "hv") plot(p.merge) pdf(paste0(figdir.save, "/CCR_Plots.pdf"), height = 5, width = 16) plot(p.merge) dev.off() ``` ```{r Figure 1E - add module scores for terms related to antigen presentation, phagocytosis, and myeloid differentiation} module.wrapper <- function(seurat.in, modules.list.in, save.path.in){ pdf(file = save.path.in, useDingbats = FALSE) for (ii in 1:length(modules.list.in)){ name.access.ii = names(modules.list.in)[[ii]] print(name.access.ii) markers.in.ii = list(modules.list.in[[ii]]) name.ii = paste0(name.access.ii, "1") if (name.ii %in% colnames(seurat.in@meta.data)){(seurat.in[[name.ii]] <- NULL)} # print(name.access.ii) if(sum(markers.in.ii[[1]] %in% rownames(seurat.in)) == 0){ seurat.in[[name.ii]] <- 0 } else{ seurat.in <- AddModuleScore(seurat.in, markers.in.ii, name = name.access.ii) } p_ii <- VlnPlot(seurat.in, name.ii, group.by = "celltypes_3", pt.size = 0) + NoLegend() + theme(plot.title = element_text(size=8)) plot(p_ii) } dev.off() } bassler.mono.like.mphage.list <- list("SPP1", "CCL2", "CLEC5A", "EMP1", "CD84", "CHIT1", "SAMSN1") bassler.lipid.list <- list("PPARG", "AKR1C3", "ZDHHC2", "CYP51A1", "PPT1", "GNA13", "HNRNPK", "ALOX5", "HADHB", "THBS1", "CBR1", "GPX4", "ACADS", "ACOT2", "ACSL1", "CYP27A1", "LRP1", "KHSRP", "NPC1", "PLCB2", "LTA4H", "RXRA", "ARID1A", "NCOA1", "ACOT7", "S1PR4", "DECR1", "ACSM3", "CSK", "ZDHHC21", "MECP2", "JAK2", "MYO5A", "NCOA2", "POR", "ABCG1", "HILPDA", "FDFT1", "AP2B1", "PECR", "GPX1", "LDLR", "MSMO1", "STARD4", "DHCR24", "LPL", "ACAT2", "SQLE", "HPGDS", "ACOT4", "SLC25A20", "ACOT1", "GPX3", "PRKAR2B", "SC5D", "FABP3", "NPC2", "YKT6", "RAP2B", "ECHS1", "TREM2", "EBP", "CYP1B1", "LPCAT2", "SPTLC2", "AP2M1", "GOLGA7", "NR1H2", "PHB", "ZDHHC24", "ACADM", "CLTA", "ARF4", "ARL1", "FDX1", "LIPA", "PRKACB", "IL1B", "SERINC5", "ACOX1", "HSD17B4", "NR3C1", "C5AR1", "CALR", "LPGAT1", "CRLS1", "ACADVL", "CPT1A", "CPT2", "STUB1", "CYP4V2", "ATP1A1", "YWHAH", "ZDHHC7", "CLOCK", "ATP11B", "CLTC", "ACAA1", "GLUL", "PRKAR1A", "SORL1", "GNAQ", "PDPK1", "AP2A1", "AP2A2", "OSBPL11", "AHR", "SPTLC2", "ACAT1", "ZDHHC17", "ZDHHC5", "HADHA", "PNPLA8", "MIA3", "PLAA", "ZDHHC3", "PRKAR2A", "SERINC1", "ZDHHC20", "CAPN2", "PTAFR", "ZDHHC6", "ACSL3", "ACAA2", "SCARB2", "PLA2G12A", "HADH", "SCP2", "HPGD", "PLBD1", "MSR1", "SOAT1", "COLEC12", "REST", "GNPAT", "TBXAS1", "IDH1", "ACSL4", "DLAT", "SERINC3", "OSBPL8", "NCEH1", "TMEM30A", "CNPY2", "CD36", "SPTSSA", "GNA15", "PRKAA1", "LPCAT3", "UNC119", "FPR2", "STOML2", "ANXA2P2", "AP2S1", "ANXA2", "ZDHHC12") cibersort.lists <- read.GMT("/Users/ctzouanas/Documents/MIT/Shalek/AncestralWisdom/Gene Lists/GMT_Gene_Lists/Cibersort_LM22_ImmuneGeneMarkers.gmt") mono.list <- cibersort.lists[["MONOCYTES"]] lists.to.plot <- list( phagocytosis.sets[names(phagocytosis.sets) == "GOBP_PHAGOCYTOSIS"][[1]], presentation.sets[names(presentation.sets) == "GOBP_ANTIGEN_PROCESSING_AND_PRESENTATION"][[1]], bassler.mono.like.mphage.list %>% as.character(), bassler.lipid.list %>% as.character(), mono.list$gene %>% as.character() ) names(lists.to.plot) <- c("GO.Phagocytosis", "GO.Antigen.Processing.and.Presentation", "Bassler.Monocyte.Like.Macrophage", "Bassler.COPD.Linked.Lipid", "CIBERSORT.Monocytes") celltypes.to.keep <- c("BAL_Monocytes", "CXCL5_Macrophages", "IFI27_Macrophages", "Inflammatory_Macrophages", "Macrophage1", "Macrophage2", "Macrophage3", "Macrophage4", "Macrophage5", "Macrophage6") monotypes.to.keep <- c("CD14_Monocytes", "CD16_Monocytes") truth.keep.vec <- BAL_PB_myeloid_srt$celltypes_2 %in% celltypes.to.keep mono.truth.keep.vec <- BAL_PB_myeloid_srt$celltypes_3 %in% monotypes.to.keep BAL_PB_myeloid_subset <- BAL_PB_myeloid_srt[, truth.keep.vec] BAL_PB_myeloid_subset$celltypes_4 <- BAL_PB_myeloid_subset$celltypes_2 mono_subset <- BAL_PB_myeloid_srt[, mono.truth.keep.vec] mono_subset$celltypes_4 <- mono_subset$celltypes_3 BAL_PB_myeloid_subset <- merge(BAL_PB_myeloid_subset, y = mono_subset) myeloid.indicator <- rep("other", length(colnames(BAL_PB_myeloid_subset))) myeloid.indicator[BAL_PB_myeloid_subset$celltypes_4 == "BAL_Monocytes"] = "BAL_Monocytes" myeloid.indicator[BAL_PB_myeloid_subset$celltypes_4 %in% c("CD14_Monocytes", "CD16_Monocytes")] = "Blood_Monocytes" myeloid.indicator[grepl(pattern = "Macrophage", x = BAL_PB_myeloid_subset$celltypes_4)] = "Macrophage" BAL_PB_myeloid_subset$myeloid.indicator <- myeloid.indicator unique(BAL_PB_myeloid_subset$myeloid.indicator) BAL_PB_myeloid_subset$myeloid.indicator <- factor(BAL_PB_myeloid_subset$myeloid.indicator, levels = c("Blood_Monocytes", "BAL_Monocytes", "Macrophage")) BAL_PB_myeloid_subset$celltypes_4 <- factor(x = BAL_PB_myeloid_subset$celltypes_4, levels = sort(unique(BAL_PB_myeloid_subset$celltypes_4))) Idents(BAL_PB_myeloid_subset) <- "celltypes_4" for (ii in 1:length(lists.to.plot)){ BAL_PB_myeloid_subset <- AddModuleScore(BAL_PB_myeloid_subset, lists.to.plot[ii], name = names(lists.to.plot)[ii], assay = "RNA", slot = "scale.data") title.ii <- gsub("\\.", " ", names(lists.to.plot)[ii]) p1 <- VlnPlot(BAL_PB_myeloid_subset, features = paste0(names(lists.to.plot)[ii], "1"), pt.size = 0, group.by = "celltypes_4") + NoLegend() + xlab("") + ggtitle(label = title.ii) plot(p1) } p1 <- VlnPlot(BAL_PB_myeloid_subset, features = "GO.Phagocytosis1", pt.size = 0, group.by = "celltypes_4") + NoLegend() + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "GO Phagocytosis") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"celltypes_4"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), axis.text.y = element_text(size = 6), text = element_text(size = 6)) p2 <- VlnPlot(BAL_PB_myeloid_subset, features = "GO.Antigen.Processing.and.Presentation1", pt.size = 0, group.by = "celltypes_4") + NoLegend() + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "GO Antigen Processing and Presentation") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"celltypes_4"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), text = element_text(size = 6), axis.text.y = element_blank()) p3 <- VlnPlot(BAL_PB_myeloid_subset, features = "Bassler.Monocyte.Like.Macrophage1", pt.size = 0, group.by = "celltypes_4") + NoLegend() + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "Bassler et al. Mono-Like Macrophage Markers") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"celltypes_4"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), text = element_text(size = 6), axis.text.y = element_blank()) p4 <- VlnPlot(BAL_PB_myeloid_subset, features = "Bassler.COPD.Linked.Lipid1", pt.size = 0, group.by = "celltypes_4") + NoLegend() + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "Bassler et al. COPD-Linked Lipid-Associated Genes") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"celltypes_4"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), text = element_text(size = 6), axis.text.y = element_blank()) plot(p1) plot(p2) plot(p3) plot(p4) p5 <- cowplot::plot_grid(plotlist = list(p1, p2, p3, p4), align = "v", ncol = 4) plot(p5) BAL_PB_myeloid_subset$myeloid.type <- "Macrophage" BAL_PB_myeloid_subset$myeloid.type[BAL_PB_myeloid_subset$celltypes_3 %in% c("CD14_Monocytes", "CD16_Monocytes")] <- "Blood Monocyte" BAL_PB_myeloid_subset$myeloid.type[BAL_PB_myeloid_subset$celltypes_3 %in% c("BAL_Monocytes")] <- "BAL Monocyte" BAL_PB_myeloid_subset$myeloid.type <- factor(BAL_PB_myeloid_subset$myeloid.type, levels = c("Blood Monocyte", "BAL Monocyte", "Macrophage")) p1 <- VlnPlot(BAL_PB_myeloid_subset, features = "GO.Phagocytosis1", pt.size = 0, group.by = "myeloid.type") + NoLegend() + geom_boxplot(width = 0.3) + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "GO Phagocytosis") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"myeloid.type"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), axis.text.y = element_text(size = 6), text = element_text(size = 6)) p2 <- VlnPlot(BAL_PB_myeloid_subset, features = "GO.Antigen.Processing.and.Presentation1", pt.size = 0, group.by = "myeloid.type") + NoLegend() + geom_boxplot(width = 0.3) + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "GO Antigen Processing and Presentation") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"myeloid.type"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), text = element_text(size = 6), axis.text.y = element_blank()) p3 <- VlnPlot(BAL_PB_myeloid_subset, features = "Bassler.Monocyte.Like.Macrophage1", pt.size = 0, group.by = "myeloid.type") + NoLegend() + geom_boxplot(width = 0.3) + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "Bassler et al. Mono-Like Macrophage Markers") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"myeloid.type"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), text = element_text(size = 6), axis.text.y = element_blank()) p4 <- VlnPlot(BAL_PB_myeloid_subset, features = "Bassler.COPD.Linked.Lipid1", pt.size = 0, group.by = "myeloid.type") + NoLegend() + geom_boxplot(width = 0.3) + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "Bassler et al. COPD-Linked Lipid-Associated Genes") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"myeloid.type"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), text = element_text(size = 6), axis.text.y = element_blank()) plot(p1) plot(p2) plot(p3) plot(p4) p5 <- cowplot::plot_grid(plotlist = list(p1, p2, p3, p4), align = "v", ncol = 4) plot(p5) wilcox.test(x = BAL_PB_myeloid_subset$GO.Phagocytosis1[BAL_PB_myeloid_subset$myeloid.type == "Macrophage"], y = BAL_PB_myeloid_subset$GO.Phagocytosis1[BAL_PB_myeloid_subset$myeloid.type == "Blood Monocyte"])$p.value wilcox.test(x = BAL_PB_myeloid_subset$GO.Phagocytosis1[BAL_PB_myeloid_subset$myeloid.type == "BAL Monocyte"], y = BAL_PB_myeloid_subset$GO.Phagocytosis1[BAL_PB_myeloid_subset$myeloid.type == "Blood Monocyte"])$p.value wilcox.test(x = BAL_PB_myeloid_subset$GO.Antigen.Processing.and.Presentation1[BAL_PB_myeloid_subset$myeloid.type == "Macrophage"], y = BAL_PB_myeloid_subset$GO.Antigen.Processing.and.Presentation1[BAL_PB_myeloid_subset$myeloid.type == "Blood Monocyte"])$p.value wilcox.test(x = BAL_PB_myeloid_subset$GO.Antigen.Processing.and.Presentation1[BAL_PB_myeloid_subset$myeloid.type == "BAL Monocyte"], y = BAL_PB_myeloid_subset$GO.Antigen.Processing.and.Presentation1[BAL_PB_myeloid_subset$myeloid.type == "Blood Monocyte"])$p.value BAL_PB_myeloid_subset$myeloid.type.rev <- factor(BAL_PB_myeloid_subset$myeloid.type, levels = rev(c("Blood Monocyte", "BAL Monocyte", "Macrophage"))) pdf(file = paste0(figdir.save, "/Fig1_MyeloidMarkersDotPlot.pdf"), useDingbats = FALSE, width = 6, height = 4) DotPlot(BAL_PB_myeloid_subset, group.by = "myeloid.type.rev", features = c("VCAN", "FCN1", "FGL2", "SORL1", "S100A8", "S100A9", "APOC1", "CES1", "SCD", "SERPING1", "FABP4", "MME"), scale = FALSE) + scale_colour_gradient2(low = "#125175", mid = "white", high = "#ad2524") + theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.title = element_blank(), legend.position = "bottom", legend.justification = "center", legend.box = "vertical") dev.off() pdf(file = paste0(figdir.save, "/Fig1E_GO_and_Bassler_Vlns_split_by_myeloidtype.pdf"), useDingbats = FALSE, width = 10, height = 2) plot(p5) dev.off() p1 <- VlnPlot(BAL_PB_myeloid_subset, features = "GO.Phagocytosis1", pt.size = 0, group.by = "myeloid.indicator", split.by = "Smoking_Status", split.plot = TRUE) + NoLegend() + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "GO Phagocytosis") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"myeloid.indicator"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), axis.text.y = element_text(size = 6), text = element_text(size = 6)) p2 <- VlnPlot(BAL_PB_myeloid_subset, features = "GO.Antigen.Processing.and.Presentation1", pt.size = 0, group.by = "myeloid.indicator", split.by = "Smoking_Status", split.plot = TRUE) + NoLegend() + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "GO Antigen Processing and Presentation") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"myeloid.indicator"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), text = element_text(size = 6), axis.text.y = element_blank()) p3 <- VlnPlot(BAL_PB_myeloid_subset, features = "Bassler.Monocyte.Like.Macrophage1", pt.size = 0, group.by = "myeloid.indicator", split.by = "Smoking_Status", split.plot = TRUE) + NoLegend() + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "Bassler et al. Mono-Like Macrophage Markers") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"myeloid.indicator"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), text = element_text(size = 6), axis.text.y = element_blank()) p4 <- VlnPlot(BAL_PB_myeloid_subset, features = "Bassler.COPD.Linked.Lipid1", pt.size = 0, group.by = "myeloid.indicator", split.by = "Smoking_Status", split.plot = TRUE) + NoLegend() + xlab("") + ylab("Module Score (A.U.)") + ggtitle(label = "Bassler et al. COPD-Linked Lipid-Associated Genes") + scale_x_discrete(limits = rev(levels(BAL_PB_myeloid_subset@meta.data[,"myeloid.indicator"]))) + coord_flip() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 6), text = element_text(size = 6), axis.text.y = element_blank()) plot(p1) plot(p2) plot(p3) plot(p4) p5 <- cowplot::plot_grid(plotlist = list(p1, p2, p3, p4), align = "v", ncol = 4) plot(p5) ##### Do statistical testing for BAL monocytes, blood monocytes, or macrophages split by smoking status Idents(BAL_PB_myeloid_subset) <- "myeloid.indicator" unique.myeloid <- unique(as.character(BAL_PB_myeloid_subset$myeloid.indicator)) gene.lists <- c("CD93", "FCGR3A") storage.df <- data.frame() for(ii in 1:length(unique.myeloid)){ for(jj in 1:length(gene.lists)){ smoker.vec.ii <- BAL_PB_myeloid_subset@assays$RNA@data[gene.lists[[jj]], (BAL_PB_myeloid_subset$myeloid.indicator == unique.myeloid[[ii]] & BAL_PB_myeloid_subset$Smoking_Status == "SMOKER")] nonsmoker.vec.ii <- BAL_PB_myeloid_subset@assays$RNA@data[gene.lists[[jj]], (BAL_PB_myeloid_subset$myeloid.indicator == unique.myeloid[[ii]] & BAL_PB_myeloid_subset$Smoking_Status == "NON_SMOKER")] wilcox.ii <- wilcox.test(x = smoker.vec.ii, y = nonsmoker.vec.ii) cohensD.ii <- cohensD(x = smoker.vec.ii, y = nonsmoker.vec.ii) row.ii <- data.frame(geneset = gene.lists[[jj]], myeloid.type = unique.myeloid[[ii]], p.val = wilcox.ii$p.value, cohensD.val = cohensD.ii, median.smoke = median(smoker.vec.ii), median.nonsmoke = median(nonsmoker.vec.ii)) storage.df <- rbind(storage.df, row.ii) } } storage.df <- storage.df %>% mutate(p.val.adj = p.adjust(p.val, method = "bonferroni")) Idents(BAL_PB_myeloid_subset) <- "myeloid.indicator" unique.myeloid <- unique(as.character(BAL_PB_myeloid_subset$myeloid.indicator)) gene.lists <- c("CD93", "FCGR3A") storage.df <- data.frame() for(ii in 1:length(unique.myeloid)){ for(jj in 1:length(gene.lists)){ celltypeofint.vec.ii <- BAL_PB_myeloid_subset@assays$RNA@data[gene.lists[[jj]], (BAL_PB_myeloid_subset$myeloid.indicator == unique.myeloid[[ii]])] bloodmono.vec.ii <- BAL_PB_myeloid_subset@assays$RNA@data[gene.lists[[jj]], (BAL_PB_myeloid_subset$myeloid.indicator == "Blood_Monocytes")] wilcox.ii <- wilcox.test(x = celltypeofint.vec.ii, y = bloodmono.vec.ii) cohensD.ii <- cohensD(x = celltypeofint.vec.ii, y = bloodmono.vec.ii) row.ii <- data.frame(geneset = gene.lists[[jj]], myeloid.type = unique.myeloid[[ii]], p.val = wilcox.ii$p.value, cohensD.val = cohensD.ii, median.smoke = median(smoker.vec.ii), median.nonsmoke = median(nonsmoker.vec.ii)) storage.df <- rbind(storage.df, row.ii) } } storage.df <- storage.df %>% mutate(p.val.adj = p.adjust(p.val, method = "bonferroni")) anova.wrapper <- function(seurat.in, feature.in, metadata.in){ feature.isolated <- seurat.in@meta.data[, feature.in] metadata.isolated <- seurat.in@meta.data[, metadata.in] aov.df <- data.frame(feature = feature.isolated, metadata = metadata.isolated) aov.out <- aov(feature ~ metadata, data = aov.df) } gophago.myeloid.aov <- anova.wrapper(seurat.in = BAL_PB_myeloid_subset, feature.in = "GO.Phagocytosis1", metadata.in = "myeloid.indicator") %>% summary() goantigen.myeloid.aov <- anova.wrapper(seurat.in = BAL_PB_myeloid_subset, feature.in = "GO.Antigen.Processing.and.Presentation1", metadata.in = "myeloid.indicator") %>% summary() gophago.celltypes.aov <- anova.wrapper(seurat.in = BAL_PB_myeloid_subset, feature.in = "Bassler.Monocyte.Like.Macrophage1", metadata.in = "celltypes_4") %>% summary() goantigen.celltypes.aov <- anova.wrapper(seurat.in = BAL_PB_myeloid_subset, feature.in = "Bassler.COPD.Linked.Lipid1", metadata.in = "celltypes_4") %>% summary() myeloid.val <- BAL_PB_myeloid_subset@meta.data[BAL_PB_myeloid_subset$myeloid.indicator == "Macrophage","GO.Phagocytosis1"] BAL.val <- BAL_PB_myeloid_subset@meta.data[BAL_PB_myeloid_subset$myeloid.indicator == "BAL_Monocytes","GO.Phagocytosis1"] Blood.val <- BAL_PB_myeloid_subset@meta.data[BAL_PB_myeloid_subset$myeloid.indicator == "Blood_Monocytes","GO.Phagocytosis1"] cohensD(x = myeloid.val, y = BAL.val) cohensD(x = Blood.val, y = BAL.val) ``` ```{r Plots for showing scoring of cells for age and gender-related gene modules} celltypes.to.keep <- c("BAL_Monocytes", "CXCL5_Macrophages", "IFI27_Macrophages", "Inflammatory_Macrophages", "Macrophage1", "Macrophage2", "Macrophage3", "Macrophage4", "Macrophage5", "Macrophage6") monotypes.to.keep <- c("CD14_Monocytes", "CD16_Monocytes") truth.keep.vec <- BAL_PB_myeloid_srt$celltypes_2 %in% celltypes.to.keep mono.truth.keep.vec <- BAL_PB_myeloid_srt$celltypes_3 %in% monotypes.to.keep BAL_PB_myeloid_subset <- BAL_PB_myeloid_srt[, truth.keep.vec] BAL_PB_myeloid_subset$celltypes_4 <- BAL_PB_myeloid_subset$celltypes_2 mono_subset <- BAL_PB_myeloid_srt[, mono.truth.keep.vec] mono_subset$celltypes_4 <- mono_subset$celltypes_3 BAL_PB_myeloid_subset <- merge(BAL_PB_myeloid_subset, y = mono_subset) BAL_PB_myeloid_subset$celltypes_4 <- factor(x = BAL_PB_myeloid_subset$celltypes_4, levels = sort(unique(BAL_PB_myeloid_subset$celltypes_4))) Idents(BAL_PB_myeloid_subset) <- "celltypes_4" bassler.mono.like.mphage.list <- list("SPP1", "CCL2", "CLEC5A", "EMP1", "CD84", "CHIT1", "SAMSN1") bassler.lipid.list <- list("PPARG", "AKR1C3", "ZDHHC2", "CYP51A1", "PPT1", "GNA13", "HNRNPK", "ALOX5", "HADHB", "THBS1", "CBR1", "GPX4", "ACADS", "ACOT2", "ACSL1", "CYP27A1", "LRP1", "KHSRP", "NPC1", "PLCB2", "LTA4H", "RXRA", "ARID1A", "NCOA1", "ACOT7", "S1PR4", "DECR1", "ACSM3", "CSK", "ZDHHC21", "MECP2", "JAK2", "MYO5A", "NCOA2", "POR", "ABCG1", "HILPDA", "FDFT1", "AP2B1", "PECR", "GPX1", "LDLR", "MSMO1", "STARD4", "DHCR24", "LPL", "ACAT2", "SQLE", "HPGDS", "ACOT4", "SLC25A20", "ACOT1", "GPX3", "PRKAR2B", "SC5D", "FABP3", "NPC2", "YKT6", "RAP2B", "ECHS1", "TREM2", "EBP", "CYP1B1", "LPCAT2", "SPTLC2", "AP2M1", "GOLGA7", "NR1H2", "PHB", "ZDHHC24", "ACADM", "CLTA", "ARF4", "ARL1", "FDX1", "LIPA", "PRKACB", "IL1B", "SERINC5", "ACOX1", "HSD17B4", "NR3C1", "C5AR1", "CALR", "LPGAT1", "CRLS1", "ACADVL", "CPT1A", "CPT2", "STUB1", "CYP4V2", "ATP1A1", "YWHAH", "ZDHHC7", "CLOCK", "ATP11B", "CLTC", "ACAA1", "GLUL", "PRKAR1A", "SORL1", "GNAQ", "PDPK1", "AP2A1", "AP2A2", "OSBPL11", "AHR", "SPTLC2", "ACAT1", "ZDHHC17", "ZDHHC5", "HADHA", "PNPLA8", "MIA3", "PLAA", "ZDHHC3", "PRKAR2A", "SERINC1", "ZDHHC20", "CAPN2", "PTAFR", "ZDHHC6", "ACSL3", "ACAA2", "SCARB2", "PLA2G12A", "HADH", "SCP2", "HPGD", "PLBD1", "MSR1", "SOAT1", "COLEC12", "REST", "GNPAT", "TBXAS1", "IDH1", "ACSL4", "DLAT", "SERINC3", "OSBPL8", "NCEH1", "TMEM30A", "CNPY2", "CD36", "SPTSSA", "GNA15", "PRKAA1", "LPCAT3", "UNC119", "FPR2", "STOML2", "ANXA2P2", "AP2S1", "ANXA2", "ZDHHC12") cibersort.lists <- read.gmt("/Users/ctzouanas/Documents/MIT/Shalek/AncestralWisdom/Gene Lists/GMT_Gene_Lists/Cibersort_LM22_ImmuneGeneMarkers.gmt") mono.list <- cibersort.lists[grepl(pattern = "MONOCYTES", x = names(cibersort.lists))] lists.to.plot <- list( presentation.sets[names(presentation.sets) == "GOBP_ANTIGEN_PROCESSING_AND_PRESENTATION"][[1]], bassler.mono.like.mphage.list %>% as.character(), bassler.lipid.list %>% as.character() ) names(lists.to.plot) <- c("GO.Antigen.Processing.and.Presentation", "Bassler.Monocyte.Like.Macrophage", "Bassler.COPD.Linked.Lipid") BAL_PB_myeloid_subset$celltypes_4 <- factor(x = BAL_PB_myeloid_subset$celltypes_4, levels = sort(unique(BAL_PB_myeloid_subset$celltypes_4))) Idents(BAL_PB_myeloid_subset) <- "celltypes_4" for (ii in 1:length(lists.to.plot)){ BAL_PB_myeloid_subset <- AddModuleScore(BAL_PB_myeloid_subset, lists.to.plot[ii], name = names(lists.to.plot)[ii], assay = "RNA", slot = "scale.data") title.ii <- gsub("\\.", " ", names(lists.to.plot)[ii]) p1 <- VlnPlot(BAL_PB_myeloid_subset, features = paste0(names(lists.to.plot)[ii], "1"), pt.size = 0, group.by = "celltypes_4") + NoLegend() + xlab("") + ggtitle(label = title.ii) plot(p1) # cluster.median <- BAL_PB_myeloid_subset@meta.data %>% select(c(paste0(names(lists.to.plot)[ii], "1"), "celltypes_4")) %>% group_by(celltypes_4) %>% dplyr::summarise(median.val = median(!!as.name(paste0(names(lists.to.plot)[ii], "1")))) %>% arrange(-median.val) # print(cluster.median) } chow.age.genes <- read.csv(file = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Chow_LungAging_SuppTable7.csv") padj.thresh <- 0.001 chow.filtered.up <- chow.age.genes %>% dplyr::filter(category == "AgeUp") %>% dplyr::select(genes) chow.filtered.down <- chow.age.genes %>% dplyr::filter(category == "AgeDown") %>% dplyr::select(genes) BAL_PB_gene_mat <- BAL_PB_myeloid_subset@assays$RNA@data pct.expr.mat <- as.matrix(rowMeans(BAL_PB_gene_mat > 0))*100 pct.expr.thresh <- 25 BAL.genes.passing.thresh <- rownames(pct.expr.mat)[pct.expr.mat > pct.expr.thresh] morrow.AM.genes <- readxl::read_excel("/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/12931_2019_1032_MOESM3_ESM/Table_S5.xls", sheet = "smoking") morrow.filtered.up <- morrow.AM.genes %>% filter(padj < padj.thresh) %>% filter(log2FoldChange > 0) morrow.filtered.down <- morrow.AM.genes %>% filter(padj < padj.thresh) %>% filter(log2FoldChange < 0) yang.gender.disc <- read.csv("/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Yang_SciReports_2019_SuppTable_Discovery.csv") yang.gender.disc.df <- yang.gender.disc %>% filter(!(Chromosome %in% c("", "X", "Y"))) yang.gender.repl <- read.csv("/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Yang_SciReports_2019_SuppTable_Replication.csv") yang.gender.repl.df <- yang.gender.repl %>% filter(!(Chromosome %in% c("", "X", "Y"))) yang.gender.disc.up <- yang.gender.disc.df %>% filter(logFC > 0) yang.gender.disc.down <- yang.gender.disc.df %>% filter(logFC < 0) yang.gender.repl.up <- yang.gender.repl.df %>% filter(logFC > 0) yang.gender.repl.down <- yang.gender.repl.df %>% filter(logFC < 0) yang.gender.up.unique <- unique(c(yang.gender.disc.up$Gene, yang.gender.repl.up$Gene)) yang.gender.down.unique <- unique(c(yang.gender.disc.down$Gene, yang.gender.repl.down$Gene)) yang.gender.up.filtered <- yang.gender.up.unique[!(yang.gender.up.unique %in% yang.gender.down.unique)] yang.gender.down.filtered <- yang.gender.down.unique[!(yang.gender.down.unique %in% yang.gender.up.unique)] angelidis.age.genes <- read.csv("//Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Angelidis_SuppTable5.csv") angelidis.age.filtered <- angelidis.age.genes[, c("gene.name", colnames(angelidis.age.genes)[grepl("Alveolar_", colnames(angelidis.age.genes))])] %>% filter(Alveolar_macrophage.p_val_adj < padj.thresh) angelidis.filtered.up <- angelidis.age.filtered %>% dplyr::filter(Alveolar_macrophage.avg_logFC > 0) %>% dplyr::select(gene.name) angelidis.filtered.down <- angelidis.age.filtered %>% dplyr::filter(Alveolar_macrophage.avg_logFC < 0) %>% dplyr::select(gene.name) if(!exists("marts")){marts <- loadmarts(host.in = "https://dec2021.archive.ensembl.org")} angelidisup.human <- convertMouseGeneList(gene.list.in = angelidis.filtered.up$gene.name, mouse.mart.in = marts[[2]], human.mart.in = marts[[1]], uniqueRows.in = T) angelidisdown.human <- convertMouseGeneList(gene.list.in = angelidis.filtered.down$gene.name, mouse.mart.in = marts[[2]], human.mart.in = marts[[1]], uniqueRows.in = T) angelidisup.module <- angelidisup.human$HGNC.symbol[angelidisup.human$HGNC.symbol %in% BAL.genes.passing.thresh] angelidisdown.module <- angelidisdown.human$HGNC.symbol[angelidisdown.human$HGNC.symbol %in% BAL.genes.passing.thresh] chowup.module <- chow.filtered.up$genes[chow.filtered.up$genes %in% BAL.genes.passing.thresh] chowdown.module <- chow.filtered.down$genes[chow.filtered.down$genes %in% BAL.genes.passing.thresh] morrowup.module <- morrow.filtered.up$symbol[morrow.filtered.up$symbol %in% BAL.genes.passing.thresh] morrowdown.module <- morrow.filtered.down$symbol[morrow.filtered.down$symbol %in% BAL.genes.passing.thresh] yangup.module <- yang.gender.up.filtered[yang.gender.up.filtered %in% BAL.genes.passing.thresh] yangdown.module <- yang.gender.down.filtered[yang.gender.down.filtered %in% BAL.genes.passing.thresh] BAL_PB_myeloid_subset <- AddModuleScore(BAL_PB_myeloid_subset, list(chowup.module), name = "Chow.Aging.Up", assay = "RNA", slot = "scale.data") BAL_PB_myeloid_subset <- AddModuleScore(BAL_PB_myeloid_subset, list(morrowup.module), name = "Morrow.Smoking.Up", assay = "RNA", slot = "scale.data") BAL_PB_myeloid_subset <- AddModuleScore(BAL_PB_myeloid_subset, list(angelidisup.module), name = "Angelidis.Aging.Up", assay = "RNA", slot = "scale.data") BAL_PB_myeloid_subset <- AddModuleScore(BAL_PB_myeloid_subset, list(yangup.module), name = "Yang.Gender.Up", assay = "RNA", slot = "scale.data") BAL_PB_myeloid_subset <- AddModuleScore(BAL_PB_myeloid_subset, list(yangdown.module), name = "Yang.Gender.Down", assay = "RNA", slot = "scale.data") cor.chowup.monomphage <- cor.test(x = BAL_PB_myeloid_subset$Bassler.Monocyte.Like.Macrophage1, y = BAL_PB_myeloid_subset$Chow.Aging.Up1, alternative = "greater", method = "pearson") cor.angelidisup.monomphage <- cor.test(x = BAL_PB_myeloid_subset$Bassler.Monocyte.Like.Macrophage1, y = BAL_PB_myeloid_subset$Angelidis.Aging.Up1, alternative = "greater", method = "pearson") cor.yangup.monomphage <- cor.test(x = BAL_PB_myeloid_subset$Bassler.Monocyte.Like.Macrophage1, y = BAL_PB_myeloid_subset$Yang.Gender.Up1, alternative = "greater", method = "pearson") cor.yangdown.monomphage <- cor.test(x = BAL_PB_myeloid_subset$Bassler.Monocyte.Like.Macrophage1, y = BAL_PB_myeloid_subset$Yang.Gender.Down1, alternative = "greater", method = "pearson") cor.morrowup.monomphage <- cor.test(x = BAL_PB_myeloid_subset$Bassler.Monocyte.Like.Macrophage1, y = BAL_PB_myeloid_subset$Morrow.Smoking.Up1, alternative = "greater", method = "pearson") cor.mat <- c(cor.chowup.monomphage$estimate, cor.angelidisup.monomphage$estimate, cor.yangup.monomphage$estimate, cor.yangdown.monomphage$estimate, cor.morrowup.monomphage$estimate) %>% as.matrix() %>% t() colnames(cor.mat) <- c("Aging-Related Geneset\n(Chow et al.)", "Aging-Related Geneset\n(Angelidis et al.)", "Male-Upregulated Geneset\n(Yang et al.)", "Female-Upregulated Geneset\n(Yang et al.)", "Smoking-Related Geneset\n(Morrow et al.)") rownames(cor.mat) <- c("Mono-Like\nMPhage Geneset\n(Bassler et al.)") cor.mat.labels <- signif(cor.mat, digits = 2) paletteLength = 100 myColor = colorRampPalette(c("#256aae", "#f7f6f6","#b21f2c"))(paletteLength) progenyBreaks = c(seq(-0.1, 0, length.out=ceiling(paletteLength/2) + 1), seq(0.3/paletteLength, 0.3, length.out=floor(paletteLength/2))) p1 <- pheatmap(cor.mat, cluster_rows = FALSE, cluster_cols = FALSE, display_numbers = cor.mat.labels, cellwidth = 100, cellheight = 100, angle_col = 90, fontsize = 16, fontsize_number = 16, color=myColor, breaks = progenyBreaks) pdf(paste0(figdir.save, "/", "Age_Sex_Smoking_Corr_Heatmap.pdf"), width = 10, height = 5, useDingbats = FALSE) p1 dev.off() ``` ```{r USE ME - Code from Marc for Monocle2} HSMM <- readRDS('/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/HSMM_method2_2.rds') ordering_genes <- readRDS('/Users/ctzouanas/Downloads/HSMM_ordering_genes_method2_2.rds') plot_cell_trajectory(HSMM, color_by = "celltypes_2")+theme(legend.position = 'right')+facet_wrap('celltypes_2') diff_test_res <- differentialGeneTest(HSMM[ordering_genes,],fullModelFormulaStr = "~sm.ns(Pseudotime)") sig_gene_names <- row.names(diff_test_res)[order(diff_test_res$qval)][1:100] plot_pseudotime_heatmap(HSMM[sig_gene_names,], num_clusters = 4, cores = 1, show_rownames = T) # # my_genes <- row.names(subset(fData(HSMM),gene_short_name %in% c("CD93"))) HSMM.pheno.df <- HSMM@phenoData@data %>% rownames_to_column(var = "cellname") pseudotime.coords <- HSMM@reducedDimS %>% t() %>% data.frame() %>% rownames_to_column(var = "cellname") %>% dplyr::rename(Component1 = X1, Component2 = X2) HSMM.pheno.df <- left_join(x = HSMM.pheno.df, y = pseudotime.coords, by = "cellname") HSMM.pheno.df <- HSMM.pheno.df %>% mutate(BALmonovec = ifelse(celltypes_2 %in% c("BAL_Monocytes"), "BAL Monocytes", "Other")) HSMM.pheno.df <- HSMM.pheno.df %>% mutate(Bloodmonovec = ifelse(celltypes_2 %in% c("CD14_Monocytes", "CD16_Monocytes"), "Blood Monocytes", "Other")) HSMM.pheno.df <- HSMM.pheno.df %>% mutate(Mphagevec = ifelse(grepl(pattern = "Macrophage", x = HSMM.pheno.df$celltypes_2), "Macrophages", "Other")) CD93.expression.df <- data.frame(cellname = colnames(BAL_PB_full), CD93.expr = BAL_PB_full@assays$RNA@data["CD93", ]) %>% dplyr::mutate(CD93.scale.expr = (CD93.expr - mean(CD93.expr)) / sd(CD93.expr)) HSMM.pheno.df <- inner_join(x = HSMM.pheno.df, y = CD93.expression.df, by = "cellname") p.all <- ggplot(HSMM.pheno.df, aes(x = Component1, y = Component2, color = celltypes_2)) + geom_point() + theme_classic() + xlim(c(-10, 11.5)) + ylim(c(-7, 4)) + theme(axis.text = element_blank(), axis.title = element_blank(), axis.ticks = element_blank()) p.BALmono <- ggplot(HSMM.pheno.df %>% dplyr::filter(BALmonovec == "BAL Monocytes"), aes(x = Component1, y = Component2, color = celltypes_2)) + geom_point() + theme_classic() + xlim(c(-10, 11.5)) + ylim(c(-7, 4)) + scale_colour_manual(values = c("#f3766e")) + theme(axis.text = element_blank(), axis.title = element_blank(), axis.ticks = element_blank()) p.Bloodmono <- ggplot(HSMM.pheno.df %>% dplyr::filter(Bloodmonovec == "Blood Monocytes"), aes(x = Component1, y = Component2, color = celltypes_2)) + geom_point() + theme_classic() + xlim(c(-10, 11.5)) + ylim(c(-7, 4)) + scale_colour_manual(values = c("#e08b26", "#bc9d2f")) + theme(axis.text = element_blank(), axis.title = element_blank(), axis.ticks = element_blank()) # ggplot(HSMM.pheno.df %>% dplyr::filter(Mphagevec == "Macrophages"), aes(x = Component1, y = Component2, color = celltypes_2)) + geom_point() + theme_classic() + xlim(c(-10, 11.5)) + ylim(c(-7, 4)) p.pseudocoord <- ggplot(HSMM.pheno.df, aes(x = Component1, y = Component2, color = Pseudotime)) + geom_point() + theme_classic() + xlim(c(-10, 11.5)) + ylim(c(-7, 4)) + theme(axis.text = element_blank(), axis.title = element_blank(), axis.ticks = element_blank()) p.5 <- plot_grid(p.all + theme(legend.position = "none"), p.pseudocoord + theme(legend.position = "none"), p.BALmono + theme(legend.position = "none"), p.Bloodmono + theme(legend.position = "none"), ncol = 2, nrow = 2, align = "hv") p.5legend <- plot_grid(p.all, p.pseudocoord, p.BALmono, p.Bloodmono, ncol = 2, nrow = 2, align = "hv") height = 4 width = 6 pdf(file = paste0(figdir.save, "/", "Pseudotime_Recoloredandsplit_v1.pdf"), height = height, width = width, useDingbats = FALSE) plot(p.5) plot(p.5legend) dev.off() tiff(paste0(figdir.save, "/", "Pseudotime_Recoloredandsplit_v1.tif"), units="in", height = height, width = width, res=300) plot(p.5) dev.off() ``` ```{r Revamp 2F} BAL_PB_full$celltypes_4 <- factor(BAL_PB_full$celltypes_2, levels = c("CD14_Monocytes", "CD16_Monocytes", "BAL_Monocytes", "CXCL5_Macrophages", "IFI27_Macrophages", "Inflammatory_Macrophages", "Macrophage1", "Macrophage2", "Macrophage3", "Macrophage4", "Macrophage5", "Macrophage6", "Proliferating_Macrophages"), labels = c("CD14 Monocytes", "CD16 Monocytes", "BAL Monocytes", "CXCL5 Macrophages", "IFI27 Macrophages", "Inflammatory Macrophages", "Macrophage1", "Macrophage2", "Macrophage3", "Macrophage4", "Macrophage5", "Macrophage6", "Proliferating Macrophages")) gene.of.int <- c("CD14" = "CD14") p.cd93.vln <- VlnPlot(BAL_PB_full, features = names(gene.of.int), pt.size = 0, group.by = "celltypes_4", slot = "data", assay = "RNA") + NoLegend() + labs(title = gene.of.int) + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), text = element_text(size = 16), plot.title = element_text(hjust = 0.5)) p.cd93.dot <- DotPlot(BAL_PB_full, features = names(gene.of.int), group.by = "celltypes_4", assay = "RNA", scale = FALSE) + scale_colour_gradient2(low = "#125175", mid = "white", high = "#ad2524") + coord_flip() + theme(axis.title = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "bottom", legend.justification = "center", text = element_text(size = 16)) + scale_x_discrete(labels = gene.of.int) pdf(file = paste0(figdir.save, "/", "Plot2F_", gene.of.int, "_VlnandDotPlot.pdf"), width = 6, height = 10, useDingbats = FALSE) p.1 = cowplot::plot_grid(p.cd93.vln, p.cd93.dot, axis = "v", align = "lr", nrow = 2, rel_heights = c(1.75, 1)) plot(p.1) dev.off() ``` ```{r Revamp figures} if(!exists("BAL_PB_myeloid_subset")){ celltypes.to.keep <- c("BAL_Monocytes", "CXCL5_Macrophages", "IFI27_Macrophages", "Inflammatory_Macrophages", "Macrophage1", "Macrophage2", "Macrophage3", "Macrophage4", "Macrophage5", "Macrophage6") monotypes.to.keep <- c("CD14_Monocytes", "CD16_Monocytes") truth.keep.vec <- BAL_PB_myeloid_srt$celltypes_2 %in% celltypes.to.keep mono.truth.keep.vec <- BAL_PB_myeloid_srt$celltypes_3 %in% monotypes.to.keep BAL_PB_myeloid_subset <- BAL_PB_myeloid_srt[, truth.keep.vec] BAL_PB_myeloid_subset$celltypes_4 <- BAL_PB_myeloid_subset$celltypes_2 mono_subset <- BAL_PB_myeloid_srt[, mono.truth.keep.vec] mono_subset$celltypes_4 <- mono_subset$celltypes_3 BAL_PB_myeloid_subset <- merge(BAL_PB_myeloid_subset, y = mono_subset) } myeloid.indicator <- rep("other", length(colnames(BAL_PB_myeloid_subset))) myeloid.indicator[BAL_PB_myeloid_subset$celltypes_4 == "BAL_Monocytes"] = "BAL_Monocytes" myeloid.indicator[BAL_PB_myeloid_subset$celltypes_4 %in% c("CD14_Monocytes", "CD16_Monocytes")] = "Blood_Monocytes" myeloid.indicator[grepl(pattern = "Macrophage", x = BAL_PB_myeloid_subset$celltypes_4)] = "Macrophage" BAL_PB_myeloid_subset$myeloid.indicator <- myeloid.indicator unique(BAL_PB_myeloid_subset$myeloid.indicator) BAL_PB_monocyte_subset <- BAL_PB_myeloid_subset[, BAL_PB_myeloid_subset$myeloid.indicator != "Macrophage"] BAL_PB_monocyte_subset$myeloid.indicator <- factor(BAL_PB_monocyte_subset$myeloid.indicator, levels = c("Blood_Monocytes", "BAL_Monocytes")) BAL_PB_monocyte_subset$myeloid.indicator.rev <- factor(BAL_PB_monocyte_subset$myeloid.indicator, levels = rev(c("Blood_Monocytes", "BAL_Monocytes"))) BAL_PB_monocyte_subset$myeloid.smoking <- factor( paste0(BAL_PB_monocyte_subset$myeloid.indicator.rev, "_", BAL_PB_monocyte_subset$Smoking_Status), levels = rev(c("Blood_Monocytes_SMOKER", "Blood_Monocytes_NON_SMOKER", "BAL_Monocytes_SMOKER", "BAL_Monocytes_NON_SMOKER")), labels = rev(c("Blood Monocytes\nSmoker", "Blood Monocytes\nNon-Smoker", "BAL Monocytes\nSmoker", "BAL Monocytes\nNon-Smoker")) ) wilcox.test(x = BAL_PB_monocyte_subset@assays$RNA@data["CCR2", BAL_PB_monocyte_subset$myeloid.smoking == "BAL Monocytes\nSmoker"], y = BAL_PB_monocyte_subset@assays$RNA@data["CCR2", BAL_PB_monocyte_subset$myeloid.smoking == "BAL Monocytes\nNon-Smoker"]) pdf(file = paste0(figdir.save, "/", "PlotFigure4_CCRComparison.pdf"), width = 4, height = 4, useDingbats = FALSE) p.1 <- DotPlot(object = BAL_PB_monocyte_subset, features = c("CCR2", "CCR3", "CCR5"), group.by = "myeloid.smoking", scale = FALSE, assay = "RNA") + scale_colour_gradient2(low = "#125175", mid = "#888888", high = "#ad2524") + theme(legend.position = "bottom", legend.justification = "center", legend.box = "vertical", axis.title = element_blank()) + coord_flip() plot(p.1) dev.off() ``` ```{r Corleis et al.: Derive markers of airspace monocytes and their smoking response} mphage.names <- c("Proliferating_Macrophages", "Macrophage5", "Macrophage3", "Macrophage6", "Macrophage4", "IFI27_Macrophages", "BAL_Monocytes", "Macrophage1", "CXCL5_Macrophages", "Inflammatory_Macrophages", "Macrophage2") Idents(BALrog.obj) <- "celltypes_2" airspace.vs.all.df <- FindMarkers(BALrog.obj, ident.1 = "BAL_Monocytes", logfc.threshold = 0.25, only.pos = TRUE) %>% arrange(-avg_log2FC) airspace.vs.all.markers <- airspace.vs.all.df %>% slice_max(order_by = avg_log2FC, n = 10) %>% rownames_to_column(var = "gene") %>% pull(gene) airspace.smoking.markers.df <- FindMarkers( BALrog.obj, ident.1 = "SMOKER", ident.2 = "NON_SMOKER", group.by = "Smoking_Status", subset.ident = "BAL_Monocytes", logfc.threshold = 0, min.pct = 0.1 ) %>% dplyr::mutate(p_val_BH_adj = p.adjust(p = p_val, method = "BH")) # VlnPlot(BALrog.obj, features = "CCR5", group.by = "celltypes_2", split.by = "Smoking_Status", slot = "data", assay = "RNA") ``` ```{r ONLY RUN ONCE - HLCA: Subset to myeloid cells and add CD93 metadata} # if(!exists("hlca.expr.obj")){ # hlca.expr.obj <- readRDS("/Users/ctzouanas/Downloads/local.rds") # } # # hlca.myeloid.keep <- # hlca.expr.obj$ann_level_3 %in% c("Macrophages", "Monocytes") | # hlca.expr.obj$ann_level_4 %in% c("Alveolar macrophages", "Classical monocytes", "Interstitial macrophages", "Non-classical monocytes") # # myeloid.hlca.expr.obj <- hlca.expr.obj[, hlca.myeloid.keep] # rm(hlca.expr.obj) # gc() # # # myeloid.hlca.expr.obj$ann_level_4 <- factor( # myeloid.hlca.expr.obj$ann_level_4, # levels = rev(c("Classical monocytes", "Non-classical monocytes", "Interstitial macrophages", "Alveolar macrophages")) # ) # # # CD93 = ENSG00000125810 # myeloid.hlca.expr.obj$CD93_count <- myeloid.hlca.expr.obj@assays$RNA@counts["ENSG00000125810", ] # myeloid.hlca.expr.obj$CD93_data <- myeloid.hlca.expr.obj@assays$RNA@data["ENSG00000125810", ] # myeloid.hlca.expr.obj$CD93_status <- # if_else(condition = myeloid.hlca.expr.obj@assays$RNA@counts["ENSG00000125810", ] > 0, # true = "CD93.Pos", # false = "CD93.Neg") # # myeloid.hlca.expr.obj$original_ann_level_3_withCD93 <- # if_else(condition = myeloid.hlca.expr.obj$original_ann_level_3 == "Monocytes" & # myeloid.hlca.expr.obj@assays$RNA@counts["ENSG00000125810", ] > 0, # true = "Monocytes_CD93+", # false = myeloid.hlca.expr.obj$original_ann_level_3) # # # hlca.expr.obj@assays$RNA <- NULL # # cd93.df <- data.frame( # cell.barcode = colnames(myeloid.hlca.expr.obj), # CD93_count = myeloid.hlca.expr.obj$CD93_count, # CD93_data = myeloid.hlca.expr.obj$CD93_data, # CD93_status = myeloid.hlca.expr.obj$CD93_status, # original_ann_level_3_withCD93 = myeloid.hlca.expr.obj$original_ann_level_3_withCD93 # ) # # saveRDS(object = myeloid.hlca.expr.obj, # file = "/Users/ctzouanas/Downloads/HLCAMyeloid_withCD93status_withRNAexpr.rds") ``` ```{r HLCA: Read in all-cell metadata and myeloid expression data from Human Lung Cell Atlas} if(!exists("hlca.metadata.obj")){ hlca.metadata.path <- "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/HLCA_emb_and_metadata.h5ad" Convert(hlca.metadata.path, dest = "h5seurat", overwrite = TRUE) hlca.metadata.obj <- LoadH5Seurat("/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/HLCA_emb_and_metadata.h5seurat") } if(!exists("myeloid.hlca.expr.obj")){ myeloid.hlca.expr.obj <- readRDS("/Users/ctzouanas/Downloads/HLCAMyeloid_withCD93status_withRNAexpr.rds") } hlca.smoking.age.df <- hlca.metadata.obj@meta.data[, c("age", "smoking_status")] %>% distinct() %>% dplyr::mutate(smoking_status = as.character(smoking_status)) %>% dplyr::filter(smoking_status != "nan") %>% dplyr::group_by(smoking_status) %>% dplyr::summarise(percentile.25 = quantile(age, probs = c(0.25)), percentile.50 = quantile(age, probs = c(0.5)), percentile.75 = quantile(age, probs = c(0.75))) ``` ```{r HLCA: Set up function to find metadata-based cell type proportions} proportion.summary.fcn <- function(seurat.in, cluster.level, split.field, participant.field, plot.order, jitter.width, str.to.remove, color.palette){ metadata.df <- seurat.in@meta.data[,c(cluster.level, split.field, participant.field)] colnames(metadata.df) <- c("cluster.level", "split.field", "participant.field") metadata.df <- metadata.df %>% mutate(particip.split.combo = paste0(participant.field, "__", split.field)) n.clusters <- length(unique(metadata.df$cluster.level)) n.participants <- length(unique(metadata.df$particip.split.combo)) summary.df <- metadata.df %>% group_by(particip.split.combo, cluster.level, .drop=FALSE) %>% summarise(cell.count = n()) %>% ungroup() %>% complete(particip.split.combo, cluster.level) summary.df$cell.count[is.na(summary.df$cell.count)] = 0 resplit <- strsplit(summary.df$particip.split.combo, split = "__") %>% data.frame() %>% t() %>% data.frame() summary.df$participant.field <- resplit$X1 summary.df <- summary.df %>% group_by(particip.split.combo) %>% mutate(participant.split.prop = cell.count/sum(cell.count)) %>% separate(particip.split.combo, into = c("particip2", "split.field"), sep = "__", remove = FALSE) %>% dplyr::select(!particip2) count.test <- summary.df %>% ungroup() %>% dplyr::count(cluster.level, .drop = FALSE) Idents(seurat.in) <- cluster.level if(!is.null(plot.order)){ y.limits <- plot.order } else{ cluster.metadata <- seurat.in[[cluster.level]][[1]] if(is.factor(cluster.metadata)){ y.limits <- levels(cluster.metadata) } else if(is.character(cluster.metadata)){ y.limits <- unique(cluster.metadata) %>% factor() } else{ print("D'oh! Improperly formatted cluster metadata! Try again, champ.") break } } y.labels <- as.character(y.limits) %>% gsub(pattern = str.to.remove, replacement = "", x = .) %>% rev() if(length(unique(count.test$n)) != 1){ stop("Participants with 0 cells in some cluster have been omitted during the process of summarizing cell counts/proportions - results may be erroneous!") } p1 <- ggplot(summary.df, aes(x = cluster.level, y = participant.split.prop, fill = split.field)) + ggtitle("Cluster Proportions Across Participants") + geom_bar(stat = "summary", position = "dodge", aes(fill = split.field)) + geom_point(position = position_jitterdodge(jitter.width = jitter.width, dodge.width = 0.9)) + scale_x_discrete(labels = factor(y.labels), limits = rev(y.limits)) + ylab("Cell Type Frequency") + theme_classic() + scale_y_continuous(expand = c(0, 0)) + scale_fill_manual(values = rev(color.palette)) + theme(plot.title = element_text(hjust = 0.5), text = element_text(size=16), axis.title.y = element_blank()) + coord_flip() output.list <- list( "output.df" = summary.df, "output.plot" = p1) return(output.list) } ``` ```{r HLCA: Try out Dirichlet regression on metadata-based smoking status} if(!exists("seurat.in")){ seurat.in = hlca.metadata.obj[, hlca.metadata.obj$anatomical_region_level_1 != "nose" & hlca.metadata.obj$smoking_status != "nan"] } smokingstatus.colors <- c("nan" = "#666666", "never" = "#4b7369", "former" = "#faaa5e", "active" = "#c85e28") chosen.metadata.field <- "ann_level_4" smoking.prop <- proportion.summary.fcn( seurat.in = seurat.in, cluster.level = chosen.metadata.field, split.field = "smoking_status", participant.field = "subject_ID", plot.order = NULL, jitter.width = 0.05, str.to.remove = "", color.palette = smokingstatus.colors ) biosex.prop <- proportion.summary.fcn( seurat.in = seurat.in, cluster.level = chosen.metadata.field, split.field = "sex", participant.field = "subject_ID", plot.order = NULL, jitter.width = 0.05, str.to.remove = "", color.palette = smokingstatus.colors ) summary.df <- smoking.prop[["output.df"]] cluster.dirich.df <- summary.df %>% dplyr::select(cluster.level, participant.field, participant.split.prop, split.field) %>% dplyr::filter(split.field != "nan") %>% dplyr::mutate(split.field = factor(split.field, levels = c("active", "former", "never"), labels = c("2", "1", "0")) %>% as.character() %>% as.numeric()) %>% pivot_wider(names_from = cluster.level, values_from = participant.split.prop) %>% column_to_rownames("particip.split.combo") %>% dplyr::select(!participant.field) AL <- DR_data(cluster.dirich.df[,colnames(cluster.dirich.df) != "split.field"]) test2 <- DirichReg(AL ~ split.field, cluster.dirich.df) summary.out <- summary(test2) df.out <- summary.out[["coef.mat"]] %>% as.data.frame() %>% rownames_to_column("variable.name") %>% dplyr::filter(grepl(pattern = "split.field", x = .$variable.name)) %>% dplyr::mutate(cluster = paste0(summary.out[["varnames"]])) %>% dplyr::select(cluster, Estimate, `Pr(>|z|)`) %>% dplyr::rename("dirich.pval" = 3, "dirich.Estimate" = 2) %>% arrange(dirich.pval) %>% arrange(-dirich.Estimate) %>% dplyr::mutate(cluster = factor(cluster, levels = rev(cluster))) %>% dplyr::mutate(log10pval = -1*log10(dirich.pval)) summary.df <- inner_join(x = summary.df, y = df.out, by = c("cluster.level" = "cluster")) font.size.chosen <- 16 p.smoking.dirich <- ggplot(df.out, aes(x = dirich.Estimate, y = cluster, size = log10pval, color = dirich.Estimate)) + geom_point() + geom_vline(xintercept = 0) + scale_color_gradient2(low = "#125175", mid = "#cccccc", high = "#ad2524") + theme_classic() + ylab("Human Lung Cell Atlas\nCell Type Annotation Level 4") + xlab("Dirichlet Regression Coefficient:\nComposition Association with Smoking Status") + theme(legend.position = "bottom", text = element_text(size = font.size.chosen), axis.text.y = element_text(size = (font.size.chosen*0.75))) plot(p.smoking.dirich) figdir <- figdir.save pdf(file = paste0(figdir, "/", "HLCA_DirichletRegression_CellTypeProportions_AnnLevel4CD93_SmokingAssociation.pdf"), width = 10, height = 8, useDingbats = FALSE) plot(p.smoking.dirich) dev.off() chosen.mono <- c("Classical monocytes", "Non-classical monocytes", "Interstitial macrophages") chosen.mono.labels <- c("Classical\nmonocytes", "Non-classical\nmonocytes", "Interstitial\nmacrophages") summary.mono.df <- summary.df %>% dplyr::filter(cluster.level %in% chosen.mono) %>% dplyr::mutate(cluster.level = factor(as.character(cluster.level), levels = rev(chosen.mono), labels = rev(chosen.mono.labels))) p.smoking.monoprop <- ggplot(summary.mono.df, aes(x = cluster.level, y = participant.split.prop, fill = split.field)) + geom_bar(stat = "summary", position = "dodge", aes(fill = split.field)) + geom_point(position = position_jitterdodge(jitter.width = 0.05, dodge.width = 0.9)) + xlab("Human Lung Cell Atlas\nCell Type Annotation Level 4") + ylab("Cell Type Frequency") + theme_classic() + scale_y_continuous(expand = c(0, 0)) + scale_fill_manual(values = c(smokingstatus.colors)) + theme(plot.title = element_text(hjust = 0.5), text = element_text(size=font.size.chosen), legend.position = "bottom") + coord_flip() plot(p.smoking.monoprop) pdf(file = paste0(figdir, "/", "HLCA_MonocyteProportions_SmokingStatus.pdf"), width = 10, height = 8, useDingbats = FALSE) plot(p.smoking.monoprop) dev.off() active.vs.never.df <- data.frame() hlca.all.celltypes <- smoking.prop$output.df$cluster.level %>% as.character() %>% unique() for(ii in 1:length(hlca.all.celltypes)){ chosen.cell.ii <- hlca.all.celltypes[[ii]] active.vals <- smoking.prop$output.df %>% dplyr::filter(split.field == "active" & cluster.level == chosen.cell.ii) %>% pull(participant.split.prop) never.vals <- smoking.prop$output.df %>% dplyr::filter(split.field == "never" & cluster.level == chosen.cell.ii) %>% pull(participant.split.prop) p.val.ii <- t.test(x = active.vals, y = never.vals) active.vs.never.ii <- data.frame( celltype = chosen.cell.ii, statistic = p.val.ii$statistic, p.val = p.val.ii$p.value ) active.vs.never.df <- rbind(active.vs.never.df, active.vs.never.ii) } active.vs.never.df <- active.vs.never.df %>% dplyr::mutate(p.val.adj = p.adjust(p = p.val, method = "BH")) biosex.colors <- c( "male" = "#57cc99", "female" = "#38a3a5" ) biosex.malevsfemale.df <- data.frame() hlca.all.celltypes <- biosex.prop$output.df$cluster.level %>% as.character() %>% unique() for(ii in 1:length(hlca.all.celltypes)){ chosen.cell.ii <- hlca.all.celltypes[[ii]] male.vals <- biosex.prop$output.df %>% dplyr::filter(split.field == "male" & cluster.level == chosen.cell.ii) %>% pull(participant.split.prop) female.vals <- biosex.prop$output.df %>% dplyr::filter(split.field == "female" & cluster.level == chosen.cell.ii) %>% pull(participant.split.prop) p.val.ii <- t.test(x = male.vals, y = female.vals) biosex.malevsfemale.ii <- data.frame( celltype = chosen.cell.ii, statistic = p.val.ii$statistic, p.val = p.val.ii$p.value ) biosex.malevsfemale.df <- rbind(biosex.malevsfemale.df, biosex.malevsfemale.ii) } biosex.malevsfemale.df <- biosex.malevsfemale.df %>% dplyr::mutate(p.val.adj = p.adjust(p = p.val, method = "BH")) p.mono.biosex <- biosex.prop$output.df %>% dplyr::filter(cluster.level %in% c("Classical monocytes", "Non-classical monocytes")) %>% dplyr::mutate(cluster.level.sex = paste0(cluster.level, ".", split.field)) %>% dplyr::mutate(cluster.level.plot = factor(cluster.level, levels = rev(c("Classical monocytes", "Non-classical monocytes")), labels = rev(c("Classical\nMono.", "Non-class.\nMono.")))) %>% ggplot(aes(x = participant.split.prop, y = cluster.level.plot, color = split.field, fill = split.field)) + geom_point(position = position_jitterdodge(jitter.width = 0.1)) + scale_color_manual(values = biosex.colors) + scale_fill_manual(values = biosex.colors) + xlab("Cell Type Proportion") + ylab("Cell Type") + geom_bar(position = "dodge", stat = "summary", fun = "mean", alpha = 0.2) + theme_classic() + theme(legend.position = "none") t.test( x = biosex.prop$output.df %>% dplyr::filter(cluster.level %in% c("Classical monocytes") & split.field == "male") %>% pull(participant.split.prop), y = biosex.prop$output.df %>% dplyr::filter(cluster.level %in% c("Classical monocytes") & split.field == "female") %>% pull(participant.split.prop) ) t.test( x = biosex.prop$output.df %>% dplyr::filter(cluster.level %in% c("Non-classical monocytes") & split.field == "male") %>% pull(participant.split.prop), y = biosex.prop$output.df %>% dplyr::filter(cluster.level %in% c("Non-classical monocytes") & split.field == "female") %>% pull(participant.split.prop) ) ``` ```{r HLCA: Try out regression on age} age.prop <- proportion.summary.fcn( seurat.in = seurat.in, cluster.level = "ann_level_4", split.field = "age", participant.field = "subject_ID", plot.order = NULL, jitter.width = 0.05, str.to.remove = "", color.palette = smokingstatus.colors ) age.mono.prop <- age.prop[["output.df"]] %>% dplyr::filter(cluster.level %in% c("Classical monocytes", "Non-classical monocytes")) %>% dplyr::mutate(split.field = as.numeric(split.field)) %>% inner_join(x = ., y = smoking.prop[["output.df"]] %>% dplyr::select(participant.field, split.field) %>% unique(), by = "participant.field") %>% dplyr::mutate(log10prop = -1*log10(participant.split.prop)) age.classicalmono.prop <- age.mono.prop %>% dplyr::filter(cluster.level == "Classical monocytes") age.nonclassicalmono.prop <- age.mono.prop %>% dplyr::filter(cluster.level == "Non-classical monocytes") classicalmono.age.spearman <- cor.test(age.classicalmono.prop$participant.split.prop, age.classicalmono.prop$split.field.x, method = "spearman") p.age.classicalmono <- ggplot(age.classicalmono.prop, aes(x = split.field.x, y = participant.split.prop, color = split.field.y)) + geom_point() + scale_color_manual(values = smokingstatus.colors) + xlab("Age (years)") + ylab("Human Lung Cell Atlas\nProportion of Classical Monocytes") + labs(color = "Smoking Status") + annotate("text", x = 30, y = 0.3, label = paste0("Spearman's Rho = ", signif(classicalmono.age.spearman$estimate, digits = 2), "\n", "P-Value = ", signif(classicalmono.age.spearman$p.value, digits = 2))) + theme_classic() + theme(legend.position = "bottom") plot(p.age.classicalmono) nonclassicalmono.age.spearman <- cor.test(age.nonclassicalmono.prop$participant.split.prop, age.nonclassicalmono.prop$split.field.x, method = "spearman") p.age.nonclassicalmono <- ggplot(age.nonclassicalmono.prop, aes(x = split.field.x, y = participant.split.prop, color = split.field.y)) + geom_point() + scale_color_manual(values = smokingstatus.colors) + xlab("Age (years)") + ylab("Human Lung Cell Atlas\nProportion of Non-Classical Monocytes") + labs(color = "Smoking Status") + annotate("text", x = 35, y = 0.2, label = paste0("Spearman's Rho = ", signif(nonclassicalmono.age.spearman$estimate, digits = 2), "\n", "P-Value = ", signif(nonclassicalmono.age.spearman$p.value, digits = 2))) + theme_classic() + theme(legend.position = "bottom") plot(p.age.nonclassicalmono) ``` ```{r Gideon et al.: Set up 4 week Seurat object and harmonize metadata} if(!exists("Granuloma4WeekSeurat")){ load("/Users/ctzouanas/Dropbox (MIT)/Reanalysis2022/Granuloma4WeekSeurat.Rdata") Granuloma4WeekSeurat <- UpdateSeuratObject(object = Granuloma4WeekSeurat) Granuloma4WeekSeurat$Round1CellTypes.plotting <- factor(as.character(Granuloma4WeekSeurat$Round1CellTypes), levels = c("B", "T", "Proliferating", "Plasma", "Mphage", "Mast", "Neutrophil", "pDC", "Endo", "Fibro", "T2P", "Club"), labels = c("B", "T", "T", "Plasma", "Macrophage", "Mast", "Neutrophil", "pDC", "Endothelial", "Fibroblast", "T2P", "Club")) granuloma.metadata <- strsplit(as.character(Granuloma4WeekSeurat$CellID), split = "_") %>% data.frame() %>% t() %>% data.frame() Granuloma4WeekSeurat$Animal <- granuloma.metadata$X2 Granuloma4WeekSeurat$Granuloma <- paste0(granuloma.metadata$X1, "_", granuloma.metadata$X2) Granuloma4Week.Metadata <- read.table(file = "/Users/ctzouanas/Dropbox (MIT)/ConstantineData/4Week/4WeekCFU.csv", header = TRUE, sep = ",") %>% dplyr::mutate(formatted_name = paste0("Array", Granuloma, "_", Animal)) %>% dplyr::filter(formatted_name %in% unique(Granuloma4WeekSeurat$Granuloma)) %>% dplyr::select(formatted_name, logCFU) %>% dplyr::rename(Granuloma = 1) merged.4week.metadata <- inner_join( x = Granuloma4WeekSeurat@meta.data, y = Granuloma4Week.Metadata, by = c("Granuloma") ) rownames(merged.4week.metadata) <- rownames(Granuloma4WeekSeurat@meta.data) Granuloma4WeekSeurat@meta.data <- merged.4week.metadata } ``` ```{r Set up functions} cor.abundance.burden <- function(object.in, cluster.level, group.var, burden.var){ granuloma.burden.df <- object.in@meta.data[c(group.var, burden.var)] %>% unique() metadata.df <- object.in@meta.data[,c(cluster.level, group.var)] %>% group_by((!!sym(cluster.level)), (!!sym(group.var))) %>% dplyr::summarise(cluster.count = n()) %>% ungroup() %>% group_by((!!sym(group.var))) %>% dplyr::mutate(cluster.prop = 100*cluster.count / sum(cluster.count)) %>% ungroup() %>% complete((!!sym(cluster.level)), (!!sym(group.var)), fill = list(cluster.prop = 0, cluster.count = 0)) metadata.df <- inner_join(x = metadata.df, y = granuloma.burden.df, by = group.var) clusters <- unique(object.in[[cluster.level]]) %>% unlist() cluster.summary.df <- data.frame() cluster.raw.df <- data.frame() for(ii in 1:length(clusters)){ cluster.ii <- clusters[[ii]] subset.df <- metadata.df %>% dplyr::filter((!!sym(cluster.level)) == cluster.ii) cor.test.val <- cor.test(x = subset.df[,"cluster.prop"] %>% unlist(), subset.df[,burden.var] %>% unlist(), method = "spearman") cluster.df.ii <- data.frame(cluster = cluster.ii, cor = cor.test.val$estimate, pval = cor.test.val$p.value) cluster.summary.df <- rbind(cluster.summary.df, cluster.df.ii) cluster.raw.df <- rbind(cluster.raw.df, subset.df) } cluster.summary.df$p.adj <- p.adjust(p = cluster.summary.df$pval, method = "BH") cluster.summary.df <- inner_join(x = cluster.summary.df, y = cluster.raw.df, by = c("cluster" = cluster.level)) return(cluster.summary.df) } ``` ```{r Gideon et al.: Map 4 week granuloma to monocytes and examine composition} if(!exists("macrophage.4week.opt.obj")){ macrophage.4week.opt.obj <- readRDS("/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Gideon_4WeekGranuloma_Macrophages_Optimized.rds") } hallmark.genesets <- geneset.loader(gsea.category.in = "H", gsea.subcategory.in = NA, species.in = "Homo sapiens") macro.4week.opt.markers <- FindAllMarkers( object = macrophage.4week.opt.obj, assay = "RNA", logfc.threshold = 0.5, min.pct = 0.25, only.pos = TRUE ) %>% group_by(cluster) %>% dplyr::arrange(-avg_log2FC, .by_group = TRUE) macro.4week.specific.markers <- macro.4week.opt.markers %>% dplyr::filter(!grepl(pattern = "^LOC", x = gene)) %>% dplyr::mutate(pct.diff = pct.1 - pct.2) %>% dplyr::group_by(cluster) %>% dplyr::slice_max(order_by = avg_log2FC, n = 2) macrophage.4week.opt.obj <- AddModuleScore( object = macrophage.4week.opt.obj, features = list(airspace.vs.all.df %>% slice_max(order_by = avg_log2FC, n = 25) %>% rownames_to_column(var = "gene") %>% pull(gene)), name = "Corleis.AirspaceMonocytes" ) BAL_PB_myeloid_srt <- AddModuleScore( object = BAL_PB_myeloid_srt, features = list(airspace.vs.all.df %>% slice_max(order_by = avg_log2FC, n = 25) %>% rownames_to_column(var = "gene") %>% pull(gene)), name = "Corleis.AirspaceMonocytes" ) inflam.hallmark <- hallmark.genesets[["HALLMARK_INFLAMMATORY_RESPONSE"]][hallmark.genesets[["HALLMARK_INFLAMMATORY_RESPONSE"]] %in% rownames(BALrog.obj@assays$RNA)] %>% unique() macro.4wk.logCFU.limmaout.df <- limmatrend_CT(obj.in = macrophage.4week.opt.obj[, macrophage.4week.opt.obj$seurat_clusters %in% c("0")], method="single.cells", metadata.to.model = c("logCFU"), metadata.to.test = c("logCFU"), sample.metadata.field="Granuloma", assay.in="RNA", slot.in = "counts", frac.express.thresh.in = 0.1) Idents(BALrog.obj) <- "celltypes_2" airspacemono.smoking.hallmarkinflam.df <- FindMarkers( BALrog.obj, ident.1 = "SMOKER", ident.2 = "NON_SMOKER", group.by = "Smoking_Status", subset.ident = "BAL_Monocytes", logfc.threshold = 0, min.pct = 0.1, assay = "RNA", features = inflam.hallmark ) %>% rownames_to_column(var = "gene") %>% arrange(-avg_log2FC) %>% dplyr::mutate(p_val_adj_BH = p.adjust(p_val, method = "BH")) p.val.thresh <- 0.05 airspace.4wkgran.combo.df <- inner_join( x = airspacemono.smoking.hallmarkinflam.df %>% dplyr::filter(p_val_adj_BH < p.val.thresh) %>% dplyr::rename(airspace_avg_log2FC = 3) %>% dplyr::select(gene, airspace_avg_log2FC) , y = macro.4wk.logCFU.limmaout.df %>% dplyr::filter(logCFU_adj.P.Val < p.val.thresh) %>% dplyr::rename(logCFU_4wk_logFC = 2) %>% dplyr::select(gene.name, logCFU_4wk_logFC), by = c("gene" = "gene.name") ) airspace.4wkgran.filter.df <- airspace.4wkgran.combo.df %>% dplyr::filter(gene %in% hallmark.genesets[["HALLMARK_INFLAMMATORY_RESPONSE"]]) macrophage.4week.opt.obj$seurat_clusters_plot <- if_else( macrophage.4week.opt.obj$seurat_clusters == "0", true = "Monocyte-Like", false = paste0("Myeloid.", macrophage.4week.opt.obj$seurat_clusters) ) %>% factor(x = ., levels = c("Monocyte-Like", paste0("Myeloid.", 1:(length(unique(macrophage.4week.opt.obj$seurat_clusters))-1)))) p.gran4wk.airspacemono.vln <- VlnPlot( macrophage.4week.opt.obj, features = "Corleis.AirspaceMonocytes1", group.by = "seurat_clusters_plot", pt.size = 0) + NoLegend() + xlab("Tuberculosis Granuloma Myeloid Clusters\nGideon et al.") + ylab("Corleis et al.:\nAirspace Monocyte Markers\nModule Score (Arb. Units)") + ggtitle("Airspace Monocyte Marker Expression\nin Tuberculosis Granuloma Myeloid Cells") + theme(plot.title = element_text(hjust = 0.5, face = "plain")) p.gran4wk.airspacemono.dot <- DotPlot( macrophage.4week.opt.obj, features = "Corleis.AirspaceMonocytes1", group.by = "seurat_clusters_plot") + coord_flip() + scale_color_gradient2(low = "#125175", mid = "#ffffff", high = "#ad2524") + ylab("Tuberculosis Granuloma Myeloid Clusters\nGideon et al.") + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1), plot.title = element_blank(), legend.position = "bottom", legend.justification = "center") p.gran4wk.CD93.vln <- VlnPlot( macrophage.4week.opt.obj, features = "CD93", group.by = "seurat_clusters_plot", pt.size = 0) + NoLegend() + xlab("Tuberculosis Granuloma Myeloid Clusters\nGideon et al.") + ylab("CD93 Expression\nlog(TP10k + 1)") + ggtitle("CD93 Expression\nin Tuberculosis Granuloma Myeloid Cells") + theme(plot.title = element_text(hjust = 0.5, face = "plain")) p.gran4wk.cd93.dot <- DotPlot( macrophage.4week.opt.obj, features = "CD93", group.by = "seurat_clusters_plot") + coord_flip() + scale_color_gradient2(low = "#125175", mid = "#ffffff", high = "#ad2524") + ylab("Tuberculosis Granuloma Myeloid Clusters\nGideon et al.") + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1), plot.title = element_blank(), legend.position = "bottom", legend.justification = "center") mac.4wk.burden.cor <- cor.abundance.burden( object.in = macrophage.4week.opt.obj, cluster.level = "seurat_clusters_plot", group.var = "Granuloma", burden.var = "logCFU") %>% dplyr::mutate(Granuloma2 = Granuloma) %>% tidyr::separate(col = Granuloma2, into = c(NA, "Animal"), sep = "_") color.monolike.burdencor <- c( "#f27f55", "#e05a54" ) p.gran4wk.monolike.burdencor <- ggplot(mac.4wk.burden.cor %>% dplyr::filter(cluster == "Monocyte-Like"), aes(x = cluster.prop, y = logCFU, color = Animal)) + geom_point(aes(shape = Animal, size = 2)) + geom_smooth(aes(x = cluster.prop, y = logCFU), method = "lm", se = FALSE, color = "#999999") + annotate("text", x = 10, y = 5.75, size = 5, label = paste0("Spearman's Rho = ", signif(mac.4wk.burden.cor %>% dplyr::filter(cluster == "Monocyte-Like") %>% pull(cor) %>% unique(), digits = 2), "\n", "BH-Adj. P-Value = N.S.")) + xlab("Percent of Airspace Monocyte-Like Cells\nAmong Granuloma Myeloid Cells (%)") + scale_color_manual(values = color.monolike.burdencor) + ylab("Bacterial Burden\nlog10(Tuberculosis CFU)") + theme_classic() + theme(legend.position = "none") ylim.chosen <- airspace.4wkgran.filter.df$airspace_avg_log2FC %>% abs() %>% max() xlim.chosen <- airspace.4wkgran.filter.df$logCFU_4wk_logFC %>% abs() %>% max() p.hallmarkinflam.airspacevs4wkgran <- ggplot(airspace.4wkgran.filter.df, aes(x = logCFU_4wk_logFC, y = airspace_avg_log2FC, label = gene)) + geom_point() + geom_text_repel(max.overlaps = 50) + geom_hline(yintercept = 0, color = "#aaaaaa", linetype = "dashed") + geom_vline(xintercept = 0, color = "#aaaaaa", linetype = "dashed") + scale_x_continuous(expand = c(0, 0.1), limits = c(-1.1*xlim.chosen, 1.1*xlim.chosen)) + scale_y_continuous(expand = c(0, 0.3), limits = c(-1.1*ylim.chosen, 1.1*ylim.chosen)) + ggtitle("Hallmark Geneset: Inflammatory Response", subtitle = "All Genes with >10% Detection & BH-Adj. P-Value < 0.05") + xlab("Tuberculosis Burden Gene Expr. Change\nAirspace Monocyte-Like Cells (Gideon et al.)") + ylab("Smoking Gene Expr. Change\nAirspace Monocyte (Corleis et al.)") + theme_classic() + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5)) p.hallmarkglycolysis.airspacevs4wkgran <- ggplot(airspace.4wkgran.combo.df %>% dplyr::filter(gene %in% hallmark.genesets$HALLMARK_GLYCOLYSIS), aes(x = logCFU_4wk_logFC, y = airspace_avg_log2FC, label = gene)) + geom_point() + geom_text_repel(max.overlaps = 100) + geom_hline(yintercept = 0, color = "#aaaaaa", linetype = "dashed") + geom_vline(xintercept = 0, color = "#aaaaaa", linetype = "dashed") + scale_x_continuous(expand = c(0, 0.1), limits = c(-1.1*xlim.chosen, 1.1*xlim.chosen)) + scale_y_continuous(expand = c(0, 0.3), limits = c(-1.1*ylim.chosen, 1.1*ylim.chosen)) + ggtitle("Hallmark Geneset: Glycolysis", subtitle = "All Genes with >10% Detection & BH-Adj. P-Value < 0.05") + xlab("Tuberculosis Burden Gene Expr. Change\nAirspace Monocyte-Like Cells (Gideon et al.)") + ylab("Smoking Gene Expr. Change\nAirspace Monocyte (Corleis et al.)") + theme_classic() + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5)) ``` ```{r Corleis et al.: Pack-years or age vs. airspace monocyte abundance} corleis.smokingstatus.colors <- c("NON_SMOKER" = "#4b7369", "SMOKER" = "#c85e28") corleis.smoking.composition.list <- proportion.summary.fcn( seurat.in = BALrog.obj, cluster.level = "celltypes_2", split.field = "Smoking_Status", participant.field = "Ragon_ID", plot.order = NULL, jitter.width = 0, str.to.remove = "", color.palette = smokingstatus.colors[c(2,4)] %>% unname()) corleis.smoking.composition.df <- corleis.smoking.composition.list$output.df %>% dplyr::mutate(pack.years = factor( participant.field, levels = ragon.id.mapping$ragon.id, labels = ragon.id.mapping$pack.years ) %>% as.character() %>% as.numeric()) %>% dplyr::mutate(age = factor( participant.field, levels = ragon.id.mapping$ragon.id, labels = ragon.id.mapping$age ) %>% as.character() %>% as.numeric()) corleis.mono.smoking.comp.df <- corleis.smoking.composition.df %>% dplyr::filter(cluster.level == "BAL_Monocytes") mono.packyear.spearman <- cor.test(corleis.mono.smoking.comp.df$participant.split.prop, corleis.mono.smoking.comp.df$pack.years, method = "spearman") p.corleis.mono.composition <- ggplot( corleis.mono.smoking.comp.df, aes(x = pack.years, y = participant.split.prop)) + scale_color_manual(values = corleis.smokingstatus.colors) + geom_point(aes(color = split.field, size = 2)) + geom_smooth(method = "lm", se = FALSE, color = "black") + annotate("text", x = 10, y = 0.03, label = paste0("Spearman's Rho = ", signif(mono.packyear.spearman$estimate, digits = 2), "\n", "P-Value = ", signif(mono.packyear.spearman$p.value, digits = 2))) + xlab("Pack Years") + ylab("Corleis et al.\nProportion of Airspace Monocytes") + theme_classic() + theme(legend.position = "bottom", text = element_text(size = font.size.chosen)) pdf(file = paste0(figdir, "/", "Corleis_MonocyteProportions_Cor_PackYears.pdf"), width = 6, height = 4, useDingbats = FALSE) plot(p.corleis.mono.composition) dev.off() ``` ```{r Corleis et al.: Monocyte inflammation changes with smoking} hallmark.geneset <- geneset.loader(gsea.category.in = "H", gsea.subcategory.in = NA, species.in = "Homo sapiens") monocyte.obj <- BAL_PB_myeloid_srt[,BAL_PB_myeloid_srt$celltypes_3 %in% c("CD14_Monocytes", "CD16_Monocytes", "BAL_Monocytes")] cytokine.df <- read.table(file = "/Users/ctzouanas/Downloads/protein_class_Predicted (1).tsv", header = TRUE, sep = "\t") %>% dplyr::filter(grepl(pattern = "Cytokine", x = Molecular.function)) %>% dplyr::filter(Gene %in% rownames(BALrog.obj)) Idents(BALrog.obj) <- "celltypes_2" mono.smoking.markers <- FindMarkers( object = BALrog.obj, ident.1 = "SMOKER", ident.2 = "NON_SMOKER", group.by = "Smoking_Status", subset.ident = "BAL_Monocytes", logfc.threshold = 0, min.pct = 0.1 ) %>% arrange(-avg_log2FC) %>% rownames_to_column(var = "gene") %>% dplyr::mutate(p_val_adj = p.adjust(p = p_val, method = "BH")) mono.smoking.cytokines <- FindMarkers( object = BALrog.obj, ident.1 = "SMOKER", ident.2 = "NON_SMOKER", group.by = "Smoking_Status", subset.ident = "BAL_Monocytes", features = cytokine.df$Gene, logfc.threshold = 0, min.pct = 0.1 ) %>% arrange(-avg_log2FC) %>% rownames_to_column(var = "gene") %>% dplyr::mutate(p_val_adj = p.adjust(p = p_val, method = "BH")) monosmoking.hallmark.gsea <- GSEA.function( DEG.list.in = mono.smoking.markers, DEG.format.in = "V4", lists.to.test = hallmark.geneset, pval.adj.sig.thresh = 1)[[1]] %>% arrange(-NES) %>% dplyr::filter(padj < 0.05) %>% dplyr::filter(grepl(pattern = "TNF|IL6|INFLAMMATORY", x = pathway)) %>% dplyr::mutate(minuslog10.padj = -1*log10(padj)) %>% dplyr::mutate(plot.name = gsub(pattern = "HALLMARK_", replacement = "", x = pathway)) %>% dplyr::mutate(plot.name = gsub(pattern = "_", replacement = " ", x = plot.name)) %>% arrange(-NES) %>% dplyr::mutate(plot.name = factor(plot.name, levels = rev(plot.name))) p.smokingmono.gsea <- ggplot(monosmoking.hallmark.gsea, aes(x = NES, y = plot.name, color = NES, size = minuslog10.padj)) + geom_point() + scale_color_gradient(low = "#cccccc", high = "#2f6131", limits = c(0, 2)) + scale_size(limits = c(1, 3), range = c(1, 5)) + scale_x_continuous(limits = c(0, NA)) + ylab("Hallmark Gene Sets") + theme_classic() + theme(legend.position = "bottom", legend.justification = "center") plot(p.smokingmono.gsea) p.smokingmono.cytokine <- DotPlot(BALrog.obj[, BALrog.obj$celltypes_2 == "BAL_Monocytes"], features = c("IL1B", "TNF", "CCL2"), group.by = "Smoking_Status", scale = FALSE) + ylab("Corleis et al.\nSmoking Status") + scale_color_gradient(low = "#cccccc", high = "#2f6131") + scale_y_discrete(labels = c("SMOKER" = "Active\nSmoker", "NON_SMOKER" = "Never\nSmoker")) + theme(legend.position = "bottom", legend.justification = "center", axis.title.x = element_blank()) plot(p.smokingmono.cytokine) ``` ```{r Examine Corleis and HLCA for glycolysis and OXPHOS} mono.mac.obj <- BAL_PB_myeloid_srt[, !BAL_PB_myeloid_srt$celltypes_2 %in% c("DC_1", "Eosinophils")] mono.mac.obj$celltypes.ratelimitenz.plot <- factor( x = if_else(mono.mac.obj$celltypes_3 == "BAL_Monocytes", true = "Airspace\nMonocytes", false = if_else(mono.mac.obj$celltypes_3 == "CD14_Monocytes", true = "PBMC: CD14\nMonocytes", false = if_else(mono.mac.obj$celltypes_3 == "CD16_Monocytes", true = "PBMC: CD16\nMonocytes", false = "Lung\nMacrophages"))), levels = c("PBMC: CD14\nMonocytes", "PBMC: CD16\nMonocytes", "Airspace\nMonocytes", "Lung\nMacrophages")) mono.mac.obj$celltypes.ratelimitenz.limma <- as.character( if_else(condition = mono.mac.obj$celltypes.ratelimitenz.plot %in% c("PBMC: CD14\nMonocytes", "PBMC: CD16\nMonocytes"), true = "PBMC\nMonocytes", false = as.character(mono.mac.obj$celltypes.ratelimitenz.plot)) ) glycolysis.sets <- c( geneset.loader(gsea.category.in = "H", gsea.subcategory.in = NA, substring.list.to.filter = "GLYCOLYSIS")[[1]], geneset.loader(gsea.category.in = "C2", gsea.subcategory.in = "CP:KEGG", substring.list.to.filter = "GLYCOLYSIS")[[1]] ) %>% unique() oxphos.set <- geneset.loader(gsea.category.in = "H", gsea.subcategory.in = NA, substring.list.to.filter = "OXIDATIVE")[[1]] glycolysis.oxphos.set <- c(glycolysis.sets, oxphos.set) %>% unique() Idents(mono.mac.obj) <- "celltypes.ratelimitenz.limma" lungmac.vs.bloodmono.glycolysis.markers <- FindMarkers( object = mono.mac.obj, ident.1 = "Lung\nMacrophages", ident.2 = "PBMC\nMonocytes", assay = "RNA", features = glycolysis.oxphos.set[glycolysis.oxphos.set %in% rownames(mono.mac.obj)], min.pct = 0, logfc.threshold = 0, only.pos = FALSE ) %>% rownames_to_column("gene") %>% arrange(-avg_log2FC, .by_group = TRUE) %>% dplyr::mutate(comparison = "LungMacro_vs_BloodMono") airspacemono.vs.bloodmono.glycolysis.markers <- FindMarkers( object = mono.mac.obj, ident.1 = "Airspace\nMonocytes", ident.2 = "PBMC\nMonocytes", assay = "RNA", features = glycolysis.oxphos.set[glycolysis.oxphos.set %in% rownames(mono.mac.obj)], min.pct = 0, logfc.threshold = 0, only.pos = FALSE ) %>% rownames_to_column("gene") %>% arrange(-avg_log2FC, .by_group = TRUE) %>% dplyr::mutate(comparison = "AirspaceMono_vs_BloodMono") lungmac.airspacemono.chosen.markers <- c( lungmac.vs.bloodmono.glycolysis.markers %>% dplyr::filter(p_val_adj < 0.01 & abs(avg_log2FC) > 0.25) %>% pull(gene), airspacemono.vs.bloodmono.glycolysis.markers %>% dplyr::filter(p_val_adj < 0.01 & abs(avg_log2FC) > 0.25) %>% pull(gene) ) %>% unique() lungmac.airspacemono_vs_bloodmono.glycolysis <- inner_join( x = lungmac.vs.bloodmono.glycolysis.markers %>% dplyr::select(gene, avg_log2FC) %>% dplyr::rename(lungmac.vs.bloodmono_avg_log2FC = 2), y = airspacemono.vs.bloodmono.glycolysis.markers %>% dplyr::select(gene, avg_log2FC) %>% dplyr::rename(airspacemono.vs.bloodmono_avg_log2FC = 2), by = "gene" ) %>% dplyr::mutate(is.glycolysis.related = if_else( condition = gene %in% glycolysis.sets, true = "Yes", false = "No" )) %>% dplyr::mutate(is.oxphos.related = if_else( condition = gene %in% oxphos.set, true = "Yes", false = "No" )) %>% dplyr::filter(gene %in% lungmac.airspacemono.chosen.markers) p.corleis.glycolysis.oxphos <- DotPlot(object = mono.mac.obj, features = c("FBP1", "AKR1A1", "ALDH3A2", "UGP2", "PCK2", "HK2", "PFKL", "PGK1", "GNPDA1", "ENO1", "PFKP"), scale = TRUE, assay = "RNA", group.by = "celltypes.ratelimitenz.plot") + ylab("Corleis et al.\nAnnotation") + theme(axis.title.x = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "bottom", legend.box = "horizontal", legend.justification = "center") + scale_color_gradient2(low = "#125175", mid = "#cccccc", high = "#ad2524") ``` ```{r Corleis et al. Milliplex data} milliplex.df <- readxl::read_xlsx( path = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/2018-12-03 25plex AMs BSL3_for revisison_CT_Formatted.xlsx") milliplex.df[milliplex.df == "OOR <"] <- "-Inf" milliplex.df[milliplex.df == "OOR >"] <- "Inf" milliplex.df.long <- milliplex.df %>% dplyr::mutate(across(where(is.double), as.character)) %>% pivot_longer(cols = !c(Patient, Acronym, Condition), names_to = "Ligand", values_to = "Conc") %>% dplyr::mutate(OOR.Tracker = if_else(Conc %in% c("-Inf", "Inf"), true = "OOR", false = "Value")) %>% dplyr::mutate(Conc = as.numeric(Conc)) milliplex.df.oor.long <- milliplex.df.long milliplex.df.oor.long$Conc[milliplex.df.long$Conc == -Inf] <- min(milliplex.df.long$Conc[milliplex.df.long$Conc != -Inf]) milliplex.df.oor.long$Conc[milliplex.df.long$Conc == Inf] <- max(milliplex.df.long$Conc[milliplex.df.long$Conc != Inf]) milliplex.df.oor.wide <- milliplex.df.oor.long %>% tidyr::pivot_wider(names_from = Condition, values_from = Conc) %>% dplyr::mutate(Mtb.vs.Ctrl = Mtb / Ctrl) %>% dplyr::arrange(-Mtb.vs.Ctrl) %>% dplyr::group_by(Acronym, Ligand) %>% dplyr::summarise(Mean.Mtb.vs.Ctrl = mean(Mtb.vs.Ctrl)) %>% ungroup() %>% dplyr::mutate(log10.Mean.Mtb.vs.Ctrl = log10(Mean.Mtb.vs.Ctrl)) milliplex.df.oor.wide.wide <- milliplex.df.oor.wide %>% dplyr::select(!Mean.Mtb.vs.Ctrl) %>% tidyr::pivot_wider(names_from = Acronym, values_from = log10.Mean.Mtb.vs.Ctrl) milliplex.df.oor.wide.wide.filter <- milliplex.df.oor.wide.wide %>% dplyr::filter(AirspaceMono > AlveolarMphage & sign(AirspaceMono) == sign(BloodMono)) %>% arrange(-AirspaceMono) %>% relocate(BloodMono, .after = Ligand) chosen.milliplex.ligands <- c("IL1a", "IL1b", "GMCSF", "TNFa", "CCL4", "CCL3") milliplex.df.oor.plot <- milliplex.df.oor.wide.wide.filter %>% dplyr::filter(Ligand %in% chosen.milliplex.ligands) %>% tidyr::pivot_longer(!Ligand, names_to = "CellType", values_to = "log10.Mean.Mtb.vs.Ctrl") milliplex.df.oor.long.plot <- milliplex.df.oor.long %>% dplyr::filter(Ligand %in% chosen.milliplex.ligands) %>% dplyr::mutate(Acronym.Condition = paste0(Acronym, ".", Condition)) %>% dplyr::mutate(Acronym.Condition = factor(Acronym.Condition, levels = c("BloodMono.Ctrl", "BloodMono.Mtb", "AirspaceMono.Ctrl", "AirspaceMono.Mtb", "AlveolarMphage.Ctrl", "AlveolarMphage.Mtb"))) %>% dplyr::mutate(Acronym = factor(Acronym, levels = c("BloodMono", "AirspaceMono", "AlveolarMphage"))) %>% dplyr::mutate(OOR.Tracker = factor(OOR.Tracker, levels = c("Value", "OOR"))) p.milliplex <- ggplot(milliplex.df.oor.long.plot, aes(x = Acronym, y = Conc, fill = Condition)) + geom_bar(position = "dodge", stat = "summary", fun = "mean", alpha = 0.2) + ylab("Concentration (pg/mL)") + geom_point(aes(color = Condition, shape = OOR.Tracker), position = position_jitterdodge(jitter.width = 0)) + facet_wrap(facets = vars(Ligand), nrow = 2, scales = "free_y") + scale_fill_manual(values = c("Ctrl" = "#888888", "Mtb" = "#ad2524")) + scale_color_manual(values = c("Ctrl" = "#888888", "Mtb" = "#ad2524")) + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.title.x = element_blank(), strip.background = element_blank(), legend.position = "bottom") pdf(file = paste0(figdir, "/", "MilliplexData_MoreMolecules.pdf"), width = 8, height = 5) plot(p.milliplex) dev.off() ``` ```{r Set up UMAP with cells colored by distance from airspace monocytes} DimPlot(BALrog.obj, group.by = "celltypes_2") umap.coords <- BALrog.obj@reductions$umap@cell.embeddings airspacemono.umap1 <- umap.coords[BALrog.obj$celltypes_2 == "BAL_Monocytes", 1] %>% median() airspacemono.umap2 <- umap.coords[BALrog.obj$celltypes_2 == "BAL_Monocytes", 2] %>% median() BALrog.obj$Airspacemono.Distance = sqrt( (umap.coords[, 1] - airspacemono.umap1)^2 + (umap.coords[, 2] - airspacemono.umap2)^2 ) # BALrog.subset <- BALrog.obj[, umap.coords[, 1] > -6 & # umap.coords[, 1] < 10 & # umap.coords[, 2] > -9 & # umap.coords[, 2] < 9.5] FeaturePlot(BALrog.obj, features = "Airspacemono.Distance") spliced.matrix <- GetAssay(BALrog.obj, assay = "spliced") unspliced.matrix <- GetAssay(BALrog.obj, assay = "unspliced") spanning.matrix <- GetAssay(BALrog.obj, assay = "spanning") RNA.matrix <- GetAssay(BALrog.obj, assay = "RNA") SCT.matrix <- GetAssay(BALrog.obj, assay = "SCT") genes.spliced <- rownames(spliced.matrix) genes.unspliced <- rownames(unspliced.matrix) genes.spanning <- rownames(spanning.matrix) genes.RNA <- rownames(RNA.matrix) genes.SCT <- rownames(SCT.matrix) genes.overlap <- Reduce(intersect, list(genes.spliced, genes.unspliced, genes.spanning, genes.RNA, genes.SCT)) BALrog.obj <- BALrog.obj[genes.overlap,] # DefaultAssay(BALrog.obj) <- "RNA" # SaveH5Seurat(BALrog.obj, filename = "/Users/ctzouanas/Dropbox (MIT)/Nails & Hammers/Smoking.Bronch.2019/Manuscript/DE_Analysis/MarcBAL_withAirspacemonoDistance_231017.h5Seurat") # Convert("/Users/ctzouanas/Dropbox (MIT)/Nails & Hammers/Smoking.Bronch.2019/Manuscript/DE_Analysis/MarcBAL_withAirspacemonoDistance_231017.h5Seurat", dest = "h5ad") ``` ```{r Set up BAL object for Alexandria SCP upload} # counts.BAL <- BALrog.obj@assays$RNA@counts # writeMM(obj = counts.BAL, file="/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/BAL_countsmatrix.mtx") # data.BAL <- BALrog.obj@assays$RNA@data # writeMM(obj = data.BAL, file="/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/BAL_datamatrix.mtx") # # # save genes and cells names # write(x = rownames(counts.BAL), file = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/BAL_countsfeatures.tsv") # write(x = colnames(counts.BAL), file = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/BAL_countsbarcodes.tsv") # # write(x = rownames(counts.BAL), file = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/BAL_datafeatures.tsv") # write(x = colnames(counts.BAL), file = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/BAL_databarcodes.tsv") BALrog.obj$CellID <- colnames(BALrog.obj) metadata.BAL.scp <- BALrog.obj@meta.data[,c("CellID", #NAME "orig.ident", #biosample_id "volunteer.id", #donor_id ##species ##species__ontology_label ##disease "Smoking_Status", #disease__ontology_label ## organ ## organ__ontology_label ##library_preparation_protocol ## library_preparation_protocol__ontology_label "Biological_Sex", "nCount_RNA", "nFeature_RNA", "celltypes_2", "pack.years")] %>% dplyr::rename(NAME = CellID, biosample_id = orig.ident, donor_id = volunteer.id, CellTypeAnnotations = celltypes_2, sex = Biological_Sex, pack_years = pack.years) %>% dplyr::mutate(sex = factor(sex, levels = c("Female", "Male", "XXXX"), labels = c("female", "male", "unknown"))) %>% mutate(across(everything(), as.character)) metadata.BAL.scp$species <- "NCBITaxon_9606" metadata.BAL.scp$species__ontology_label <- "Homo sapiens" metadata.BAL.scp$disease <- "PATO_0000461" metadata.BAL.scp$disease__ontology_label <- "normal" metadata.BAL.scp$paired_ends <- "True" metadata.BAL.scp$end_bias <- "3 prime end bias" metadata.BAL.scp$biosample_type <- "PrimaryBioSample_Tissue" metadata.BAL.scp$library_preparation_protocol <- "EFO_0008919" metadata.BAL.scp$library_preparation_protocol__ontology_label <- "Seq-Well" metadata.BAL.scp$organ <- "UBERON_0002048" metadata.BAL.scp$organ__ontology_label <- "lung" scp.BAL.metadata.header <- data.frame(matrix(ncol = ncol(metadata.BAL.scp), nrow = 0)) colnames(scp.BAL.metadata.header) <- colnames(metadata.BAL.scp) scp.BAL.metadata.header[1,] <- colnames(metadata.BAL.scp) scp.BAL.metadata.header[2, colnames(scp.BAL.metadata.header) %in% c("nCount_RNA", "nFeature_RNA", "pack_years")] <- "numeric" scp.BAL.metadata.header[2, is.na(scp.BAL.metadata.header[2,])] <- "group" scp.BAL.metadata.export <- rbind(scp.BAL.metadata.header, metadata.BAL.scp) %>% relocate(NAME) scp.BAL.metadata.export[2,1] <- "TYPE" write.table(x = scp.BAL.metadata.export, file = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/BAL_alexandria_structured_metadata3.txt", row.names = FALSE, col.names = FALSE, sep = "\t", quote = FALSE) ``` ```{r Set up PBMC object for Alexandria SCP upload} counts.PBMC <- PBj.obj@assays$RNA@counts writeMM(obj = counts.PBMC, file="/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/PBMC_countsmatrix.mtx") data.PBMC <- PBj.obj@assays$RNA@data writeMM(obj = data.PBMC, file="/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/PBMC_datamatrix.mtx") # save genes and cells names write(x = rownames(counts.PBMC), file = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/PBMC_countsfeatures.tsv") write(x = colnames(counts.PBMC), file = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/PBMC_countsbarcodes.tsv") write(x = rownames(counts.PBMC), file = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/PBMC_datafeatures.tsv") write(x = colnames(counts.PBMC), file = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/PBMC_databarcodes.tsv") ported.metadata <- BALrog.obj@meta.data[, c("volunteer.id", "Biological_Sex", "pack.years")] %>% distinct() PBj.obj$CellID <- colnames(PBj.obj) metadata.PBMC.scp <- PBj.obj@meta.data[,c("CellID", #NAME "orig.ident", #biosample_id "volunteer.id", #donor_id ##species ##species__ontology_label ##disease "Smoking_Status", #disease__ontology_label ## organ ## organ__ontology_label ##library_preparation_protocol ## library_preparation_protocol__ontology_label "nCount_RNA", "nFeature_RNA")] %>% inner_join(x = ., y = ported.metadata, by = "volunteer.id") %>% dplyr::rename(NAME = CellID, biosample_id = orig.ident, donor_id = volunteer.id, sex = Biological_Sex, pack_years = pack.years) %>% dplyr::mutate(sex = factor(sex, levels = c("Female", "Male", "XXXX"), labels = c("female", "male", "unknown"))) %>% mutate(across(everything(), as.character)) metadata.PBMC.scp$species <- "NCBITaxon_9606" metadata.PBMC.scp$species__ontology_label <- "Homo sapiens" metadata.PBMC.scp$disease <- "PATO_0000461" metadata.PBMC.scp$disease__ontology_label <- "normal" metadata.PBMC.scp$paired_ends <- "True" metadata.PBMC.scp$end_bias <- "3 prime end bias" metadata.PBMC.scp$biosample_type <- "PrimaryBioSample_Tissue" metadata.PBMC.scp$library_preparation_protocol <- "EFO_0008919" metadata.PBMC.scp$library_preparation_protocol__ontology_label <- "Seq-Well" metadata.PBMC.scp$organ <- "UBERON_0000178" metadata.PBMC.scp$organ__ontology_label <- "blood" scp.PBMC.metadata.header <- data.frame(matrix(ncol = ncol(metadata.PBMC.scp), nrow = 0)) colnames(scp.PBMC.metadata.header) <- colnames(metadata.PBMC.scp) scp.PBMC.metadata.header[1,] <- colnames(metadata.PBMC.scp) scp.PBMC.metadata.header[2, colnames(scp.PBMC.metadata.header) %in% c("nCount_RNA", "nFeature_RNA", "pack_years")] <- "numeric" scp.PBMC.metadata.header[2, is.na(scp.PBMC.metadata.header[2,])] <- "group" scp.PBMC.metadata.export <- rbind(scp.PBMC.metadata.header, metadata.PBMC.scp) %>% relocate(NAME) scp.PBMC.metadata.export[2,1] <- "TYPE" write.table(x = scp.PBMC.metadata.export, file = "/Users/ctzouanas/Documents/MIT/Shalek/BAL_Bjorn_Marc/Code_Data_Upload/Alexandria_SCP_Upload/PBMC_alexandria_structured_metadata.txt", row.names = FALSE, col.names = FALSE, sep = "\t", quote = FALSE) ```