```{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)
```