```{r}
library(tidyverse)
library(ggrepel)
library(circlize)
library(viridis)
library(patchwork)
library(ggpubr)
library(cowplot)
library(rstatix)
library(pheatmap)
library(Seurat)
library(fgsea)
library(msigdbr)
library(readxl)
library(lsr)
library(biomaRt)
library(janitor)
library(ggh4x)
library(scales)
library(pROC)
library(onewaytests)
library(nichenetr)
library(RColorBrewer)
library(Matrix)
library(DirichletReg)
library(openxlsx)

group1 <- c("pDC", "Mac2", "Mac8", "Mac5", "Mac1", "Mac3")
group1.color <- "#faaa5e"
group2 <- c("Mast", "Fibro", "Mac4", "T1P", "Endo", "Plasma")
group2.color <- "#d7003d"
group3 <- c("TProliferating", "Mac7", "TT1T17Pop3", "TMetallothionein", "TT1T17Pop4", "TCytotoxic4", "TT1T17Pop1", "TCytotoxic6", "TStemLike", "TCytotoxic2", "TSRRM2")
group3.color <- "#0160cc"
group4 <- c("Neutrophil", "TT1T17Pop2", "TTreg", "Mac9", "RBC", "T2P")
group4.color <- "#6ca659"
group5 <- c("TCytotoxic1", "TCytotoxic3", "TCytotoxic5", "B", "cDC", "Mac6", "TInterferon")
group5.color <- "#4e4ca0"
groupT <- c("TProliferating", "TT1T17Pop3", "TMetallothionein", "TT1T17Pop4", "TCytotoxic4", "TT1T17Pop1", "TCytotoxic6", "TStemLike", "TCytotoxic2", "TSRRM2", "TT1T17Pop2", "TTreg",
            "TCytotoxic1", "TCytotoxic3", "TCytotoxic5", "TInterferon")

font.size = 16
figdir <- "path-to-figures"
```

```{r Set up Seurat objects and harmonize metadata}
load("/Users/ctzouanas/Dropbox (MIT)/ConstantineData/SeuratObjects/MacrophageSeurat.Rdata")
load("/Users/ctzouanas/Dropbox (MIT)/ConstantineData/SeuratObjects/GranulomaTotalSeurat.Rdata")
load("/Users/ctzouanas/Documents/MIT/Shalek/TB_IntercellularInteraction/MetaMerge.Rdata")
load("/Users/ctzouanas/Dropbox (MIT)/ConstantineData/SeuratObjects/TcellSeuratUpdate.Rdata")
load("/Users/ctzouanas/Documents/MIT/Shalek/TB_IntercellularInteraction/CFUDataFinal_Update.Rdata")

MetaMerge <- full_join(x = MetaMerge,
                       y = CFUDataFinal_Update,
                       by = c("Granuloma"))

MacrophageSeurat <- UpdateSeuratObject(object = MacrophageSeurat)
TcellSeuratUpdate <- UpdateSeuratObject(object = TcellSeuratUpdate)
GranulomaTotalSeurat <- UpdateSeuratObject(object = GranulomaTotalSeurat)

lowburden.granulomas <- CFUDataFinal_Update %>% dplyr::filter(HighLowUpdate == 0)
lowburden.granulomas <- lowburden.granulomas$Granuloma %>% as.character()
highburden.granulomas <- CFUDataFinal_Update %>% dplyr::filter(HighLowUpdate == 1)
highburden.granulomas <- highburden.granulomas$Granuloma %>% as.character()

old.granulomas <- CFUDataFinal_Update %>% dplyr::filter(Age == "Old")
old.granulomas <- old.granulomas$Granuloma %>% as.character()
new.granulomas <- CFUDataFinal_Update %>% dplyr::filter(Age == "New")
new.granulomas <- new.granulomas$Granuloma %>% as.character()

seurat.cellnames <- colnames(GranulomaTotalSeurat) %>% data.frame()
colnames(seurat.cellnames) <- "CellID"
orig.metadata <- GranulomaTotalSeurat@meta.data %>% rownames_to_column("CellID")
joined.metadata <- left_join(x = seurat.cellnames, y = orig.metadata, by = "CellID")
joined.metadata <- left_join(x = joined.metadata, y = MetaMerge, by = "CellID")
joined.metadata <- joined.metadata %>% column_to_rownames("CellID")

GranulomaTotalSeurat <- AddMetaData(object = GranulomaTotalSeurat, metadata = joined.metadata)
GranulomaTotalSeurat$CellTypeGroup <- ""
GranulomaTotalSeurat$CellTypeGroup[GranulomaTotalSeurat$SpecificFinal %in% group1] <- "Group1"
GranulomaTotalSeurat$CellTypeGroup[GranulomaTotalSeurat$SpecificFinal %in% group2] <- "Group2"
GranulomaTotalSeurat$CellTypeGroup[GranulomaTotalSeurat$SpecificFinal %in% group3] <- "Group3"
GranulomaTotalSeurat$CellTypeGroup[GranulomaTotalSeurat$SpecificFinal %in% group4] <- "Group4"
GranulomaTotalSeurat$CellTypeGroup[GranulomaTotalSeurat$SpecificFinal %in% group5] <- "Group5"
GranulomaTotalSeurat <- GranulomaTotalSeurat[, GranulomaTotalSeurat$CellTypeGroup != ""]

GranulomaTotalSeurat$GranulomaBurden = ""
GranulomaTotalSeurat$GranulomaBurden[GranulomaTotalSeurat$Granuloma %in% lowburden.granulomas] <- "LowBurden"
GranulomaTotalSeurat$GranulomaBurden[GranulomaTotalSeurat$Granuloma %in% highburden.granulomas] <- "HighBurden"

GranulomaTotalSeurat$GranulomaAge = ""
GranulomaTotalSeurat$GranulomaAge[GranulomaTotalSeurat$Granuloma %in% old.granulomas] <- "Old"
GranulomaTotalSeurat$GranulomaAge[GranulomaTotalSeurat$Granuloma %in% new.granulomas] <- "New"

GranulomaTotalSeurat@reductions$umap <- CreateDimReducObject(embeddings = as.matrix(joined.metadata[,c("X_umap1", "X_umap2")]), key = "UMAP_", assay = "RNA")

DimPlot(GranulomaTotalSeurat, group.by = "GranulomaBurden") + xlim(c(-15, 15)) + ylim(c(-15, 15)) + theme(legend.position = "none")

macrophage.metadata <- GranulomaTotalSeurat@meta.data[grepl(pattern = "Mac[0-9]", x = GranulomaTotalSeurat$SpecificFinal),]
nrow(macrophage.metadata) == ncol(MacrophageSeurat)
T.metadata <- GranulomaTotalSeurat@meta.data[grepl(pattern = "^T[A-Z]", x = GranulomaTotalSeurat$SpecificFinal),]
nrow(T.metadata) == ncol(TcellSeuratUpdate)

macrophage.metadata <- macrophage.metadata[colnames(MacrophageSeurat),]
MacrophageSeurat <- AddMetaData(object = MacrophageSeurat, metadata = macrophage.metadata)
MacrophageSeurat@reductions$umap <- CreateDimReducObject(embeddings = as.matrix(MacrophageSeurat@meta.data[,c("X_umap1", "X_umap2")]), key = "UMAP_", assay = "RNA")
T.metadata <- T.metadata[colnames(TcellSeuratUpdate),]
TcellSeuratUpdate <- AddMetaData(object = TcellSeuratUpdate, metadata = T.metadata)
TcellSeuratUpdate@reductions$umap <- CreateDimReducObject(embeddings = as.matrix(TcellSeuratUpdate@meta.data[,c("X_umap1", "X_umap2")]), key = "UMAP_", assay = "RNA")

DimPlot(MacrophageSeurat, group.by = "SpecificFinal", reduction = "umap") + xlim(c(-15, 15)) + ylim(c(-15, 15)) + theme(legend.position = "none")
DimPlot(TcellSeuratUpdate, group.by = "SpecificFinal", reduction = "umap") + xlim(c(-15, 15)) + ylim(c(-15, 15)) + theme(legend.position = "none")
```

```{r Set up function to evaluate p-values and correlations across cell types}
cor.abundance.burden <- function(object.in, cluster.level, group.var, burden.var, split.var){
  granuloma.burden.df <- object.in@meta.data[c(group.var, burden.var, split.var)] %>% unique()
  
  metadata.df <- object.in@meta.data[,c(cluster.level, group.var)] %>%
    group_by((!!sym(cluster.level)), (!!sym(group.var))) %>%
    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)
}

dirichwrapper.highlowburden <- function(cluster.storage.df.in, color.palette.in = c("#056c9a", "#ae2425"), nrow.in = 2, font.size = 16){
  cluster.dirich.df <- cluster.storage.df.in %>%
    dplyr::select(cluster, Granuloma, cluster.prop, GranulomaBurden) %>%
    dplyr::mutate(GranulomaBurden = factor(GranulomaBurden, levels = c("LowBurden", "HighBurden"), labels = c("0", "1")) %>% as.character() %>% as.numeric()) %>%
    dplyr::mutate(cluster.prop = cluster.prop / 100) %>%
    pivot_wider(names_from = cluster, values_from = cluster.prop) %>%
    column_to_rownames("Granuloma")
  AL <- DR_data(cluster.dirich.df[,colnames(cluster.dirich.df) != "GranulomaBurden"])
  # test1 <- DirichReg(AL ~ 1, cluster.dirich.df)
  test2 <- DirichReg(AL ~ GranulomaBurden, cluster.dirich.df)
  # anova(test1, test2)
  summary.out <- summary(test2)
  df.out <- summary.out[["coef.mat"]] %>% as.data.frame() %>% rownames_to_column("variable.name") %>% 
    dplyr::filter(grepl(pattern = "GranulomaBurden", 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)
  cluster.storage.df.in <- inner_join(x = cluster.storage.df.in, y = df.out, by = "cluster")
  cluster.storage.df.in$plot.title <- paste0(cluster.storage.df.in$cluster, 
                                             "\nDirich. Est. = ", signif(cluster.storage.df.in$dirich.Estimate, digits = 2),
                                             "\nDirich. P-val = ", signif(cluster.storage.df.in$dirich.pval, digits = 2))
  cluster.storage.df.in$GranulomaBurden <- factor(cluster.storage.df.in$GranulomaBurden, levels = c("LowBurden", "HighBurden"), 
                                                  labels = c("Low", "High"))
  p.dirich.out <- cluster.storage.df.in %>% 
    ggplot(aes(x = GranulomaBurden, y = cluster.prop, color = GranulomaBurden)) + 
    geom_boxplot(aes(color = GranulomaBurden), outlier.shape = NA) + 
    geom_point(size = 3, position = "jitter") + 
    scale_color_manual(values = color.palette.in) +
    xlab("Granuloma Burden") + 
    facet_wrap(~ plot.title, nrow = nrow.in, scales = "free") + 
    theme_classic() + 
    ylab("Percent of Cells in Granuloma (%)") +
    theme(text = element_text(size = font.size), legend.position = "none",
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          strip.background = element_blank(), strip.text = element_text(size = font.size))
  dirich.output <- vector(mode = "list", length = 2)
  dirich.output[[1]] <- p.dirich.out
  dirich.output[[2]] <- cluster.storage.df.in
  return(dirich.output)
}

GranulomaTotalSeurat$SpecificFinal_GroupedTT1T17 <- GranulomaTotalSeurat$SpecificFinal
GranulomaTotalSeurat$SpecificFinal_GroupedTT1T17[grepl(pattern = "^TT1T17", x = GranulomaTotalSeurat$SpecificFinal)] <- "TT1T17"

cluster.storage.SpecificFinal.CFU <- cor.abundance.burden(object.in = GranulomaTotalSeurat,
                                                          cluster.level = "SpecificFinal_GroupedTT1T17",
                                                          group.var = "Granuloma",
                                                          burden.var = "logCFU",
                                                          split.var = "GranulomaBurden")

cluster.storage.SpecificFinal.CFUCEQ <- cor.abundance.burden(object.in = GranulomaTotalSeurat,
                                                          cluster.level = "SpecificFinal_GroupedTT1T17",
                                                          group.var = "Granuloma",
                                                          burden.var = "Log10.CFU.CEQ",
                                                          split.var = "GranulomaBurden")
cluster.storage.GenericFinal.CFU <- cor.abundance.burden(object.in = GranulomaTotalSeurat,
                                                          cluster.level = "GenericFinal",
                                                          group.var = "Granuloma",
                                                          burden.var = "logCFU",
                                                          split.var = "GranulomaBurden")

cluster.storage.GenericFinal.CFUCEQ <- cor.abundance.burden(object.in = GranulomaTotalSeurat,
                                                          cluster.level = "GenericFinal",
                                                          group.var = "Granuloma",
                                                          burden.var = "Log10.CFU.CEQ",
                                                          split.var = "GranulomaBurden")

cluster.storage.SpecificFinal.CFU.highresTT1T17 <- cor.abundance.burden(object.in = GranulomaTotalSeurat,
                                                                        cluster.level = "SpecificFinal",
                                                                        group.var = "Granuloma",
                                                                        burden.var = "logCFU",
                                                                        split.var = "GranulomaBurden")
cluster.storage.SpecificFinal.CFUCEQ.highresTT1T17 <- cor.abundance.burden(object.in = GranulomaTotalSeurat,
                                                                           cluster.level = "SpecificFinal",
                                                                           group.var = "Granuloma",
                                                                           burden.var = "Log10.CFU.CEQ",
                                                                           split.var = "GranulomaBurden")
##### Generic Cell Types: Composition vs. CFU or CFU/CEQ #####
cluster.storage.GenericFinal.CFU$cluster <- factor(cluster.storage.GenericFinal.CFU$cluster,
                                       levels = c("B", "cDC", "Endo", "Fibro", "Mast", "Mphage", "Neutrophil", "pDC", "Plasma", "RBC", "T", "T1P", "T2P"),
                                       labels = c("B", "cDC", "Endothelial", "Fibroblast", "Mast", "Macrophage", "Neutrophil", "pDC", "Plasma", "RBC", "T/NK", "T1P", "T2P"))
cluster.storage.GenericFinal.CFUCEQ$cluster <- factor(cluster.storage.GenericFinal.CFUCEQ$cluster,
                                       levels = c("B", "cDC", "Endo", "Fibro", "Mast", "Mphage", "Neutrophil", "pDC", "Plasma", "RBC", "T", "T1P", "T2P"),
                                       labels = c("B", "cDC", "Endothelial", "Fibroblast", "Mast", "Macrophage", "Neutrophil", "pDC", "Plasma", "RBC", "T/NK", "T1P", "T2P"))

p.generic.CFU <- cluster.storage.GenericFinal.CFU %>%
  dplyr::mutate(title = paste0(cluster, "\nRho = ", round(cor, digits = 2), " q = ", signif(p.adj, digits = 2))) %>%
  ggplot(aes(x = logCFU, y = cluster.prop, color = GranulomaBurden)) + 
  scale_color_manual(values = c("#f7931e", "#029151")) + 
  geom_point() + 
  ggtitle("Cell Type Composition vs. CFU") + 
  facet_wrap(~ title, nrow = 2, scales = "free") + theme_classic() + 
  xlab(expression(~log[10]~CFU~per~Granuloma)) + 
  ylab("Percent of All Cells in Granuloma (%)") + 
  theme(strip.background = element_blank(), text = element_text(size = 12), 
        strip.text = element_text(size = 12), legend.position = "none", plot.title = element_text(hjust = 0.5))
p.generic.CFU.CEQ <- cluster.storage.GenericFinal.CFUCEQ %>%
  dplyr::mutate(title = paste0(cluster, "\nRho = ", round(cor, digits = 2), " q = ", signif(p.adj, digits = 2))) %>%
  ggplot(aes(x = Log10.CFU.CEQ, y = cluster.prop, color = GranulomaBurden)) + 
  scale_color_manual(values = c("#f7931e", "#029151")) + 
  geom_point() + 
  ggtitle("Cell Type Composition vs. CFU/CEQ") + 
  facet_wrap(~ title, nrow = 2, scales = "free") + theme_classic() + 
  xlab(expression(~log[10]~CFU/CEQ~per~Granuloma)) + 
  ylab("Percent of All Cells in Granuloma (%)") + 
  theme(strip.background = element_blank(), text = element_text(size = 12), 
        strip.text = element_text(size = 12), legend.position = "none", plot.title = element_text(hjust = 0.5))

##### Macrophages: Composition vs. CFU or CFU/CEQ #####
mac.corr.CFU.plot <- rbind(cluster.storage.SpecificFinal.CFU, cluster.storage.GenericFinal.CFU) %>%
  dplyr::filter(cluster != "RBC") %>%
  dplyr::filter(grepl(pattern = "^Mac", x = .$cluster)) %>%
  dplyr::mutate(cluster = factor(cluster,
                                 levels = c("Macrophage", paste0("Mac", 1:9)),
                                 labels = c("All Macrophages", paste0("Mac", 1:9)))) %>%
  arrange(cluster)
mac.corr.CFU.plot$p.adj[grepl(pattern = "^Mac", x = mac.corr.CFU.plot$cluster)] <-
  p.adjust(p = mac.corr.CFU.plot$pval[grepl(pattern = "^Mac", x = mac.corr.CFU.plot$cluster)], method = "BH")
  
p.macrophage.specific.CFU <- mac.corr.CFU.plot %>%
  dplyr::mutate(title = paste0(cluster, "\nRho = ", round(cor, digits = 2), " q = ", signif(p.adj, digits = 2))) %>%
  ggplot(aes(x = logCFU, y = cluster.prop, color = GranulomaBurden)) + 
  scale_color_manual(values = c("#f7931e", "#029151")) + 
  geom_point() + 
  ggtitle("Macrophage Subclusters vs. CFU") + 
  facet_wrap(~ title, ncol = 5, scales = "free") + theme_classic() + 
  xlab(expression(~log[10]~CFU~per~Granuloma)) + 
  ylab("Percent of All Cells in Granuloma (%)") + 
  theme(strip.background = element_blank(), text = element_text(size = 12), 
        strip.text = element_text(size = 12), legend.position = "none", plot.title = element_text(hjust = 0.5))

mac.corr.CFUCEQ.plot <- rbind(cluster.storage.SpecificFinal.CFUCEQ, cluster.storage.GenericFinal.CFUCEQ) %>%
  dplyr::filter(cluster != "RBC") %>%
  dplyr::filter(grepl(pattern = "^Mac", x = .$cluster)) %>%
  dplyr::mutate(cluster = factor(cluster,
                                 levels = c("Macrophage", paste0("Mac", 1:9)),
                                 labels = c("All Macrophages", paste0("Mac", 1:9)))) %>%
  arrange(cluster)
mac.corr.CFUCEQ.plot$p.adj[grepl(pattern = "^Mac", x = mac.corr.CFUCEQ.plot$cluster)] <-
  p.adjust(p = mac.corr.CFUCEQ.plot$pval[grepl(pattern = "^Mac", x = mac.corr.CFUCEQ.plot$cluster)], method = "BH")
p.macrophage.specific.CFUCEQ <- mac.corr.CFUCEQ.plot %>%
  dplyr::mutate(title = paste0(cluster, "\nRho = ", round(cor, digits = 2), " q = ", signif(p.adj, digits = 2))) %>%
  ggplot(aes(x = Log10.CFU.CEQ, y = cluster.prop, color = GranulomaBurden)) + 
  scale_color_manual(values = c("#f7931e", "#029151")) + 
  geom_point() + 
  ggtitle("Macrophage Subclusters vs. CFU/CEQ") + 
  facet_wrap(~ title, ncol = 5, scales = "free") + theme_classic() + 
  xlab(expression(~log[10]~CFU/CEQ~per~Granuloma)) + 
  ylab("Percent of All Cells in Granuloma (%)") + 
  theme(strip.background = element_blank(), text = element_text(size = 12), 
        strip.text = element_text(size = 12), legend.position = "none", plot.title = element_text(hjust = 0.5))

##### T: Composition vs. CFU or CFU/CEQ #####
T.corr.CFU.plot <- rbind(cluster.storage.SpecificFinal.CFU, cluster.storage.GenericFinal.CFU) %>%
  dplyr::filter(cluster != "RBC") %>%
  dplyr::filter(grepl(pattern = "^T", x = .$cluster)) %>%
  dplyr::filter(!cluster %in% c("T1P", "T2P")) %>%
  dplyr::mutate(cluster = factor(cluster,
                                 levels = c("T/NK", paste0("TCytotoxic", 1:6), "TStemLike", "TTreg", "TInterferon", "TMetallothionein", "TSRRM2", "TProliferating", "TT1T17"),
                                 labels = c("All T/NK", paste0("Cytotoxic", 1:6), "Stem-Like", "Treg", "IFN Responsive", "Metallothionein", "SRRM2+ T Cells", "Proliferating", "T1-T17 Cells"))) %>%
  arrange(cluster)
T.corr.CFU.plot$p.adj[T.corr.CFU.plot$cluster != "All T/NK"] <-
  p.adjust(p = T.corr.CFU.plot$pval[T.corr.CFU.plot$cluster != "All T/NK"], method = "BH")

p.T.specific.CFU <- T.corr.CFU.plot %>%
  dplyr::mutate(title = paste0(cluster, "\nRho = ", round(cor, digits = 2), " q = ", signif(p.adj, digits = 2))) %>%
  ggplot(aes(x = logCFU, y = cluster.prop, color = GranulomaBurden)) + 
  scale_color_manual(values = c("#f7931e", "#029151")) + 
  geom_point() + 
  ggtitle("T Cell Subclusters vs. CFU") + 
  facet_wrap(~ title, ncol = 7, scales = "free") + theme_classic() + 
  xlab(expression(~log[10]~CFU~per~Granuloma)) + 
  ylab("Percent of All Cells in Granuloma (%)") + 
  theme(strip.background = element_blank(), text = element_text(size = 12), 
        strip.text = element_text(size = 12), legend.position = "none", plot.title = element_text(hjust = 0.5))

T.corr.CFUCEQ.plot <- rbind(cluster.storage.SpecificFinal.CFUCEQ, cluster.storage.GenericFinal.CFUCEQ) %>%
  dplyr::filter(cluster != "RBC") %>%
  dplyr::filter(grepl(pattern = "^T", x = .$cluster)) %>%
  dplyr::filter(!cluster %in% c("T1P", "T2P")) %>%
  dplyr::mutate(cluster = factor(cluster,
                                 levels = c("T/NK", paste0("TCytotoxic", 1:6), "TStemLike", "TTreg", "TInterferon", "TMetallothionein", "TSRRM2", "TProliferating", "TT1T17"),
                                 labels = c("All T/NK", paste0("Cytotoxic", 1:6), "Stem-Like", "Treg", "IFN Responsive", "Metallothionein", "SRRM2+ T Cells", "Proliferating", "T1-T17 Cells"))) %>%
  arrange(cluster)
T.corr.CFUCEQ.plot$p.adj[T.corr.CFUCEQ.plot$cluster != "All T/NK"] <-
  p.adjust(p = T.corr.CFUCEQ.plot$pval[T.corr.CFUCEQ.plot$cluster != "All T/NK"], method = "BH")

p.T.specific.CFUCEQ <- T.corr.CFUCEQ.plot %>%
  dplyr::mutate(title = paste0(cluster, "\nRho = ", round(cor, digits = 2), " q = ", signif(p.adj, digits = 2))) %>%
  ggplot(aes(x = Log10.CFU.CEQ, y = cluster.prop, color = GranulomaBurden)) + 
  scale_color_manual(values = c("#f7931e", "#029151")) + 
  geom_point() + 
  ggtitle("T Cell Subclusters vs. CFU/CEQ") + 
  facet_wrap(~ title, ncol = 7, scales = "free") + theme_classic() + 
  xlab(expression(~log[10]~CFU/CEQ~per~Granuloma)) + 
  ylab("Percent of All Cells in Granuloma (%)") + 
  theme(strip.background = element_blank(), text = element_text(size = 12), 
        strip.text = element_text(size = 12), legend.position = "none", plot.title = element_text(hjust = 0.5))

# comparison.list <- list(p.generic.CFU, p.generic.CFU.CEQ, p.macrophage.specific.CFU, p.macrophage.specific.CFUCEQ, p.T.specific.CFU, p.T.specific.CFUCEQ)
# pdf(file = paste0(figdir, "/", "CellTypeCompositionvsBurdenMetrics.pdf"), height = 24, width = 14)
# cowplot::plot_grid(plotlist = comparison.list, nrow = 6, align = "hv", axis = "tblr", rel_heights = c(1,1,1,1,1,1), labels = "AUTO")
# dev.off()
# 
# comparison.list <- list(p.generic.CFU, p.generic.CFU.CEQ, p.macrophage.specific.CFU, p.macrophage.specific.CFUCEQ, p.T.specific.CFU, p.T.specific.CFUCEQ)
# pdf(file = paste0(figdir, "/", "CellTypeCompositionvsBurdenMetrics_layout.pdf"), height = 12, width = 28)
# cowplot::plot_grid(plotlist = comparison.list, nrow = 3, align = "hv", axis = "tblr", rel_heights = c(1,1,1), rel_widths = c(1,1), labels = "AUTO")
# dev.off()
##### T1-T17: Composition vs. CFU or CFU/CEQ #####
T1T17.corr.CFU.plot <- rbind(cluster.storage.SpecificFinal.CFU, cluster.storage.SpecificFinal.CFU.highresTT1T17) %>%
  dplyr::filter(cluster != "RBC") %>%
  dplyr::filter(grepl(pattern = "^TT1T17", x = .$cluster)) %>%
  dplyr::mutate(cluster = factor(cluster,
                                 levels = c("TT1T17", paste0("TT1T17Pop", 1:4)),
                                 labels = c("All T1-T17 Cells", paste0("T1-T17 Pop. ", 1:4)))) %>%
  arrange(cluster)
T1T17.corr.CFU.plot$p.adj[T1T17.corr.CFU.plot$cluster != "All T1-T17 Cells"] <-
  p.adjust(p = T1T17.corr.CFU.plot$pval[T1T17.corr.CFU.plot$cluster != "All T1-T17 Cells"], method = "BH")

p.T1T17.specific.CFU <- T1T17.corr.CFU.plot %>%
  dplyr::mutate(title = paste0(cluster, "\nRho = ", round(cor, digits = 2), " q = ", signif(p.adj, digits = 2))) %>%
  ggplot(aes(x = logCFU, y = cluster.prop, color = GranulomaBurden)) + 
  scale_color_manual(values = c("#f7931e", "#029151")) + 
  geom_point() + 
  ggtitle("T1-T17 Subclusters vs. CFU") + 
  facet_wrap(~ title, ncol = 7, scales = "free") + theme_classic() + 
  xlab(expression(~log[10]~CFU~per~Granuloma)) + 
  ylab("Percent of All Cells in Granuloma (%)") + 
  theme(strip.background = element_blank(), text = element_text(size = 12), 
        strip.text = element_text(size = 12), legend.position = "none", plot.title = element_text(hjust = 0.5))

T1T17.corr.CFUCEQ.plot <- rbind(cluster.storage.SpecificFinal.CFUCEQ, cluster.storage.SpecificFinal.CFUCEQ.highresTT1T17) %>%
  dplyr::filter(cluster != "RBC") %>%
  dplyr::filter(grepl(pattern = "^TT1T17", x = .$cluster)) %>%
  dplyr::mutate(cluster = factor(cluster,
                                 levels = c("TT1T17", paste0("TT1T17Pop", 1:4)),
                                 labels = c("All T1-T17 Cells", paste0("T1-T17 Pop. ", 1:4)))) %>%
  arrange(cluster)
T1T17.corr.CFU.plot$p.adj[T1T17.corr.CFU.plot$cluster != "All T1-T17 Cells"] <-
  p.adjust(p = T1T17.corr.CFU.plot$pval[T1T17.corr.CFU.plot$cluster != "All T1-T17 Cells"], method = "BH")

p.T1T17.specific.CFUCEQ <- T1T17.corr.CFUCEQ.plot %>%
  dplyr::mutate(title = paste0(cluster, "\nRho = ", round(cor, digits = 2), " q = ", signif(p.adj, digits = 2))) %>%
  ggplot(aes(x = Log10.CFU.CEQ, y = cluster.prop, color = GranulomaBurden)) + 
  scale_color_manual(values = c("#f7931e", "#029151")) + 
  geom_point() + 
  ggtitle("T1-T17 Subclusters vs. CFU/CEQ") + 
  facet_wrap(~ title, ncol = 7, scales = "free") + theme_classic() + 
  xlab(expression(~log[10]~CFU/CEQ~per~Granuloma)) + 
  ylab("Percent of All Cells in Granuloma (%)") + 
  theme(strip.background = element_blank(), text = element_text(size = 12), 
        strip.text = element_text(size = 12), legend.position = "none", plot.title = element_text(hjust = 0.5))

comparison.list <- list(p.generic.CFU, p.generic.CFU.CEQ, p.macrophage.specific.CFU, p.macrophage.specific.CFUCEQ, p.T.specific.CFU, p.T.specific.CFUCEQ, p.T1T17.specific.CFU, p.T1T17.specific.CFUCEQ)
pdf(file = paste0(figdir, "/", "CellTypeCompositionvsBurdenMetrics.pdf"), height = 32, width = 14)
cowplot::plot_grid(plotlist = comparison.list, nrow = 8, align = "hv", axis = "tblr", rel_heights = c(1,1,1,1,1,1,0.75,0.75), labels = "AUTO")
dev.off()

comparison.list <- list(p.generic.CFU, p.generic.CFU.CEQ, p.macrophage.specific.CFU, p.macrophage.specific.CFUCEQ, p.T.specific.CFU, p.T.specific.CFUCEQ, p.T1T17.specific.CFU, p.T1T17.specific.CFUCEQ)
pdf(file = paste0(figdir, "/", "CellTypeCompositionvsBurdenMetrics_4by2layout.pdf"), height = 16, width = 28)
cowplot::plot_grid(plotlist = comparison.list, nrow = 4, align = "hv", axis = "tblr", rel_heights = c(1,1,1,0.75), rel_widths = c(1,1), labels = "AUTO")
dev.off()

dirich.generic <- dirichwrapper.highlowburden(cluster.storage.df.in = cluster.storage.GenericFinal.CFU, font.size = 12, nrow = 2)
dirich.specific <- dirichwrapper.highlowburden(cluster.storage.df.in = cluster.storage.SpecificFinal.CFU, font.size = 12, nrow = 4)
dirich.specific.splitt1t17 <- dirichwrapper.highlowburden(cluster.storage.df.in = cluster.storage.SpecificFinal.CFU.highresTT1T17, 
                                                            font.size = 12, nrow = 4)
p.dirich.generic <- dirich.generic[[1]]
p.dirich.specific <- dirich.specific[[1]]
p.dirich.specific.splitt1t17 <- dirich.specific.splitt1t17[[1]]

df.dirich.generic <- dirich.generic[[2]]
df.dirich.specific <- dirich.specific[[2]]
df.dirich.specific.splitt1t17 <- dirich.specific.splitt1t17[[2]]
# plot(p.dirich.specific)
# macro.storage.SpecificFinal.CFU <- cor.abundance.burden(object.in = MacrophageSeurat,
#                                                         cluster.level = "SpecificFinal",
#                                                         group.var = "Granuloma",
#                                                         burden.var = "logCFU",
#                                                         split.var = "GranulomaBurden")
# p.dirich.macrophage <- dirichwrapper.highlowburden(cluster.storage.df.in = macro.storage.SpecificFinal.CFU, 
#                                                    font.size = 12, nrow = 2)
# plot(p.dirich.macrophage)
# TcellSeuratUpdate$SpecificFinal_GroupedTT1T17 <- GranulomaTotalSeurat$SpecificFinal
# TcellSeuratUpdate$SpecificFinal_GroupedTT1T17[grepl(pattern = "^TT1T17", x = TcellSeuratUpdate$SpecificFinal)] <- "TT1T17"
# T.storage.SpecificFinal.CFU <- cor.abundance.burden(object.in = TcellSeuratUpdate,
#                                                     cluster.level = "SpecificFinal_GroupedTT1T17",
#                                                     group.var = "Granuloma",
#                                                     burden.var = "logCFU",
#                                                     split.var = "GranulomaBurden")
# p.dirich.T <- dirichwrapper.highlowburden(cluster.storage.df.in = T.storage.SpecificFinal.CFU, 
#                                           font.size = 12, nrow = 2)
# plot(p.dirich.T)
# TT1T17.storage.SpecificFinal.CFU <- cor.abundance.burden(object.in = TcellSeuratUpdate[,grepl(pattern = "^TT1T17", x = TcellSeuratUpdate$SpecificFinal)],
#                                                     cluster.level = "SpecificFinal",
#                                                     group.var = "Granuloma",
#                                                     burden.var = "logCFU",
#                                                     split.var = "GranulomaBurden")
# p.dirich.T1T17 <- dirichwrapper.highlowburden(cluster.storage.df.in = TT1T17.storage.SpecificFinal.CFU, 
#                                               font.size = 12, nrow = 2)
dirich.list <- list(p.dirich.generic, p.dirich.specific, p.dirich.specific.splitt1t17)
pdf(file = paste0(figdir, "/", "CellTypeCompositionDirichletRegression.pdf"), height = 26, width = 20)
cowplot::plot_grid(plotlist = dirich.list, nrow = 3, align = "hv", axis = "tblr", rel_heights = c(1,2,2), labels = "AUTO", label_size = 24)
dev.off()

celltype.cor.dirich.list <- list(df.dirich.generic, df.dirich.specific, df.dirich.specific.splitt1t17)
for(ii in 1:length(celltype.cor.dirich.list)){
  df.ii <- celltype.cor.dirich.list[[ii]]
  df.ii <- df.ii %>% dplyr::select(!plot.title) %>%
    dplyr::rename("Cluster Name" = "cluster", 
                  "Spearman Correlation" = "cor",
                  "Spearman p-value" = "pval",
                  "Spearman p-value with Benj-Hoch mult. testing correction" = "p.adj",
                  "Cell Count" = "cluster.count",
                  "Cell Proportion" = "cluster.prop",
                  "Log10 CFU" = "logCFU",
                  "Granuloma Burden Bin" = "GranulomaBurden",
                  "Dirichlet coefficient" = "dirich.Estimate",
                  "Dirichlet p-value" = "dirich.pval")
  celltype.cor.dirich.list[[ii]] <- df.ii
}
names(celltype.cor.dirich.list) <- c("Canonical Cell Types", "Subclustered Cell Types", "Subclust w Sep T1-T17 Subpop")

write.xlsx(x = celltype.cor.dirich.list,
           file = paste0(figdir, "/", "Table S5-Association with bacterial burden.xlsx"),
           overwrite = TRUE)
```

```{r Produce Fig. S2H}
GranulomaTotalSeurat.4017 <- GranulomaTotalSeurat[,GranulomaTotalSeurat$Animal.x == "4017"]

metadata.4017.df <- GranulomaTotalSeurat.4017@meta.data[,c("GenericFinal", "Granuloma")]
summary.4017.df <- metadata.4017.df %>%
  group_by(GenericFinal, Granuloma) %>%
  summarise(celltypecount = n()) %>%
  ungroup() %>%
  complete(GenericFinal, Granuloma, fill = list(celltypecount = 0)) %>%
  group_by(Granuloma) %>% 
  dplyr::mutate(celltypeprop = 100 * celltypecount / sum(celltypecount)) %>%
  ungroup() %>% 
  dplyr::mutate(GranulomaAge = factor(.$Granuloma, 
                                      levels = CFUDataFinal_Update$Granuloma,
                                      labels = CFUDataFinal_Update$Age))

celltypes.4017 <- unique(summary.4017.df$GenericFinal)
pval.4017.df <- data.frame()
for(ii in 1:length(celltypes.4017)){
  celltype.ii <- celltypes.4017[[ii]]
  summary.4017.df.ii <- summary.4017.df %>% dplyr::filter(GenericFinal == celltype.ii)
  oldburden.vals.ii <- summary.4017.df.ii %>% dplyr::filter(GranulomaAge == "Old") %>% pull(celltypeprop)
  newburden.vals.ii <- summary.4017.df.ii %>% dplyr::filter(GranulomaAge == "New") %>% pull(celltypeprop)
  wilcox.ii <- wilcox.test(x = oldburden.vals.ii,
                           y = newburden.vals.ii)$p.value
  pval.df.ii <- data.frame(celltype = celltype.ii, p.val = wilcox.ii)
  pval.4017.df <- rbind(pval.4017.df, pval.df.ii)
}

pval.4017.df$p.adj <- p.adjust(p = pval.4017.df$p.val,
                               method = "BH")
summary.4017.df$GranulomaAge <- factor(summary.4017.df$GranulomaAge,
                                       levels = c("New", "Old"),
                                       labels = c("Late", "Early"))
summary.4017.df <- inner_join(x = summary.4017.df,
                              y = pval.4017.df,
                              by = c("GenericFinal" = "celltype"))
summary.4017.df$GenericFinal <- factor(summary.4017.df$GenericFinal,
                                       levels = c("B", "cDC", "Endo", "Fibro", "Mast", "Mphage", "Neutrophil", "pDC", "Plasma", "RBC", "T", "T1P", "T2P"),
                                       labels = c("B", "cDC", "Endothelial", "Fibroblast", "Mast", "Macrophage", "Neutrophil", "pDC", "Plasma", "RBC", "T/NK", "T1P", "T2P"))

summary.4017.df$plot.title <- paste0(summary.4017.df$GenericFinal, "\nq = ", signif(summary.4017.df$p.adj, digits = 2))
granuloma.colors <- c("#9b2e29",
                      "#fdd377",
                      "#f7911e",
                      "#3d57a7",
                      "#015395",
                      "#901b4f",
                      "#7956a4",
                      "#ee3324",
                      "#7ad2f7",
                      "#259859")

dirich.4017.df <- summary.4017.df %>%
  dplyr::select(GenericFinal, Granuloma, celltypeprop, GranulomaAge) %>%
  dplyr::mutate(GranulomaNum = as.numeric(as.character(factor(GranulomaAge, levels = c("Late", "Early"), labels = c("0", "1"))))) %>%
  dplyr::mutate(celltypeprop = celltypeprop / 100) %>%
  pivot_wider(names_from = GenericFinal, values_from = celltypeprop) %>%
  column_to_rownames("Granuloma") %>%
  dplyr::select(-GranulomaAge)
AL.4017 <- DR_data(dirich.4017.df[,colnames(dirich.4017.df) != "GranulomaNum"])
test1.4017 <- DirichReg(AL.4017 ~ 1, dirich.4017.df)
test2.4017 <- DirichReg(AL.4017 ~ GranulomaNum, dirich.4017.df)
anova(test1.4017, test2.4017)
summary.4017 <- summary(test2.4017)
df.4017.out <- summary.4017[["coef.mat"]] %>% as.data.frame() %>% rownames_to_column("variable.name") %>% 
    dplyr::filter(grepl(pattern = "GranulomaNum", x = .$variable.name)) %>%
    dplyr::mutate(cluster = summary.4017[["varnames"]]) %>%
    dplyr::select(cluster, Estimate, `Pr(>|z|)`) %>%
    dplyr::rename("dirich.pval" = 3, "dirich.Estimate" = 2) %>%
    arrange(dirich.pval)
summary.4017.df.in <- inner_join(x = summary.4017.df, y = df.4017.out, by = c("GenericFinal" = "cluster")) %>% rename("cluster" = 1)
summary.4017.df.in$plot.title <- paste0(summary.4017.df.in$cluster, 
                                        "\nEst. = ", signif(summary.4017.df.in$dirich.Estimate, digits = 2),
                                        "\nP-val = ", signif(summary.4017.df.in$dirich.pval, digits = 2))
summary.4017.df.in$GranulomaBurden <- factor(summary.4017.df.in$GranulomaAge, levels = c("Late", "Early"), 
                                             labels = c("Late", "Early"))
set.seed(11)
p.dirich.4017.earlylate <- summary.4017.df.in %>% 
  ggplot(aes(x = GranulomaAge, y = celltypeprop, color = Granuloma)) + 
  geom_boxplot(aes(color = GranulomaAge), outlier.shape = NA) + 
  geom_point(size = 3, position = "jitter") + 
  scale_color_manual(values = granuloma.colors) +
  xlab("Granuloma Age") + 
  facet_wrap(~ plot.title, nrow = 1, scales = "free") + 
  theme_classic() + 
  ylab("Percent of Cells in Granuloma (%)") +
  theme(text = element_text(size = 12), legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1),
        strip.background = element_blank(), strip.text = element_text(size = 12))

set.seed(11)
font.size = 16
p.celltypes.props.4017 <- summary.4017.df %>% 
  ggplot(aes(x = GranulomaAge, y = celltypeprop, color = Granuloma)) + 
  geom_boxplot(aes(color = GranulomaAge), outlier.shape = NA) + 
  geom_point(size = 3, position = "jitter") + 
  scale_color_manual(values = granuloma.colors) +
  facet_wrap(~ plot.title, ncol = length(celltypes.4017), scales = "free") + 
  theme_classic() + 
  ylab("Percent of Cells in Granuloma (%)") +
  theme(text = element_text(size = font.size), legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1),
        strip.background = element_blank(), strip.text = element_text(size = 16))

p.granuloma.cfu.4017 <- CFUDataFinal_Update %>% dplyr::filter(Animal == "4017") %>%
  dplyr::mutate(GranulomaAge = factor(.$Age, levels = c("New", "Old"), labels = c("Late", "Early"))) %>%
  ggplot(aes(x = GranulomaAge, y = logCFU, color = Granuloma)) + 
  geom_boxplot(aes(color = GranulomaAge), outlier.shape = NA) + 
  geom_point(size = 3, position = "jitter") + 
  scale_color_manual(values = granuloma.colors) +
  theme_classic() + 
  ylab(expression(~log[10]~CFU~per~Granuloma)) +
  theme(text = element_text(size = font.size), legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))

plotlist.4017 <- list(p.granuloma.cfu.4017, p.celltypes.props.4017)
pdf(file = paste0(figdir, "/", "Animal4017_Analyses.pdf"), width = 18, height = 6)
cowplot::plot_grid(plotlist = plotlist.4017, ncol = 2, rel_widths = c(0.1, 1), axis = "tblr", align = "hv")
dev.off()

pdf(file = paste0(figdir, "/", "Animal4017_DirichletAnalyses.pdf"), width = 18, height = 6)
plot(p.dirich.4017.earlylate)
dev.off()


```

```{r Produce Fig. S6}
mac.color <- c("Mac1" = "#216fa4",
               "Mac2" = "#1a9497",
               "Mac3" = "#289245",
               "Mac4" = "#669840",
               "Mac5" = "#a5b238",
               "Mac6" = "#fcb442",
               "Mac7" = "#f69a1e",
               "Mac8" = "#ed6b23",
               "Mac9" = "#c8262e")
if(!exists("fraction.express")){
  MacrophageSeurat$SpecificFinal <- factor(MacrophageSeurat$SpecificFinal,
                                           levels = names(mac.color))
  font.size.in = 12
  
  # Quick detour to calculate percent expression of genes in macrophage object
  binary.matrix <- as.matrix(MacrophageSeurat@assays[["RNA"]]@counts)
  binary.matrix <- as.matrix(rowSums(binary.matrix > 0))
  fraction.express <- binary.matrix / ncol(MacrophageSeurat)
  fraction.express <- rownames_to_column(as.data.frame(fraction.express))
  colnames(fraction.express) <- c("gene.name", "frac.express")
}

set.seed(100)
chosen.genes <- fraction.express[fraction.express$frac.express > 0.25, "gene.name"]

cronan.supp <- read.csv(file = "/Users/ctzouanas/Documents/MIT/Shalek/TB_IntercellularInteraction/UsefulFiles/CronanSupplemental.csv") %>% clean_names()
colnames(cronan.supp) <- gsub(pattern = "x[0-9]+_", replacement = "", x = colnames(cronan.supp))
colnames(cronan.supp)[[1]] <- "gene"
cronan.supp <- data.frame(cronan.supp)
avg_epithelialized_macrophage <- rowMeans(cronan.supp[,grepl(pattern = "epithelialized_macrophages", x = colnames(cronan.supp))])
avg_epithelialized_macrophage <- data.frame(gene = cronan.supp$gene, avg_epithelialized_macrophage = avg_epithelialized_macrophage)
avg_all_others <- rowMeans(cronan.supp[,!grepl(pattern = "epithelialized_macrophages|gene", x = colnames(cronan.supp))])
avg_all_others <- data.frame(gene = cronan.supp$gene, avg_all_others = avg_all_others)
avg_join_df <- inner_join(x = avg_epithelialized_macrophage,
                          y = avg_all_others,
                          by = "gene") %>% 
  dplyr::mutate(avg_log2FC = log2(avg_epithelialized_macrophage / avg_all_others)) %>%
  arrange(-avg_log2FC) %>%
  dplyr::filter(avg_log2FC > 0.75)
cronan.epithelioid <- avg_join_df$gene

if(!exists("marts.rhesus")){
  marts.rhesus <- useMart('ensembl', dataset = 'mmulatta_gene_ensembl')
  marts.zebrafish <- useMart('ensembl', dataset = 'drerio_gene_ensembl')
}

cronan.m1 <- c("nos2b", "il1b", "il12ba", "il6", "cxcl11.1")

all.zebrafish <- c(cronan.epithelioid, cronan.m1, cronan.m2)
all.rhesus = getLDS(attributes = c("external_gene_name"), filters = "external_gene_name", values = all.zebrafish, mart = marts.zebrafish, 
                 attributesL = c("external_gene_name"), martL = marts.rhesus, 
                 uniqueRows=FALSE)

rhesus.epithelioid <- all.rhesus %>% dplyr::filter(Gene.name %in% cronan.epithelioid) %>% 
  dplyr::filter(Gene.name.1 %in% chosen.genes) %>%
  pull(Gene.name.1) %>% unique()
rhesus.m1 <- all.rhesus %>% dplyr::filter(Gene.name %in% cronan.m1) %>% pull(Gene.name.1) %>% unique()

MacrophageSeurat <- AddModuleScore(object = MacrophageSeurat,
                                   features = list(rhesus.epithelioid),
                                   name = "Cronan.Epithelioid")
MacrophageSeurat <- AddModuleScore(object = MacrophageSeurat,
                                   features = list(rhesus.m1),
                                   name = "Cronan.M1")

genes.of.int <- c("CDH1", "TJP1", "STAT6", "ARG1")

p.cronanepi.vln <- VlnPlot(MacrophageSeurat, features = "Cronan.Epithelioid1", group.by = "SpecificFinal",
                           pt.size = 0, cols = rev(mac.color)) + NoLegend() +
  ylab("Module Score (Arb. Units)") +
  xlab("This Study: Macrophage Cluster") +
  scale_x_discrete(limits = rev(names(mac.color))) + 
  ggtitle("Cronan et al.\nEpithelioid Macrophage Module") +
  coord_flip() + 
  theme(plot.title = element_text(face = "plain", hjust = 0.5),
        text = element_text(size = font.size.in),
        axis.text.y = element_text(color = rev(mac.color)),
        axis.text.x = element_text(angle = 0, hjust = 0.5))

p.cronanepi.dot <- DotPlot(object = MacrophageSeurat, 
                          features = c("CDH1", "TJP1", "STAT6", 
                                       "CD36", "ALDH9A1", "MGST3", "FABP5", "APOC1", 
                                       "MRC1", "MS4A7",
                                       "ANXA1", "ANXA2", "MMP14", "FN1", "ITGA6", "ECM1", 
                                       "CD93"),
                          group.by = "SpecificFinal",
                          cols = c("lightgrey", "#ad2524")) + 
  scale_y_discrete(limits = rev(names(mac.color))) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.text.y = element_text(color = rev(mac.color)),
        legend.position = "bottom", legend.justification = "center", legend.box = "vertical",
        axis.title = element_blank(), text = element_text(size = font.size.in)) +
  scale_colour_gradient2(low = "#125175", mid = "white", high = "#ad2524")

mac4.features <- c("CXCL8", "AREG", "NFKBIA", "IFI27", "IFI30",
                   "IL1B", "LYZ", "HIF1A", 
                   "GPR183", "INSIG1",  "EREG", 
                   "C1QA", "C1QB", "C1QC", 
                   "B2M", "CD74", 
                   "OLR1",
                   "CTSH", "CTSB", "CTSC", "CTSS", "CTSL", "CSTB")
p.mac4.dotplot <- DotPlot(object = MacrophageSeurat, 
                          features = mac4.features,
                          group.by = "SpecificFinal",
                          cols = c("lightgrey", "#ad2524")) + 
  scale_y_discrete(limits = rev(names(mac.color))) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.text.y = element_text(color = rev(mac.color)), 
        legend.position = "bottom", legend.justification = "center", legend.box = "vertical",
        axis.title = element_blank(), text = element_text(size = font.size.in)) +
  scale_colour_gradient2(low = "#125175", mid = "white", high = "#ad2524")

mphage.metadata.df <- MacrophageSeurat@meta.data[,c("SpecificFinal", "GranulomaBurden", "Cronan.Epithelioid1", "Cronan.M11")] 
mphage.module.plot <- mphage.metadata.df %>%
  group_by(SpecificFinal) %>%
  summarise(mean.Cronan.Epithelioid1 = mean(Cronan.Epithelioid1),
            sd.Cronan.Epithelioid1 = sd(Cronan.Epithelioid1),
            mean.Cronan.M11 = mean(Cronan.M11),
            sd.Cronan.M11 = sd(Cronan.M11))
p.cronan.epivsm1 <- ggplot(mphage.module.plot,
                           aes(x = mean.Cronan.M11,
                               y = mean.Cronan.Epithelioid1, label = SpecificFinal, color = SpecificFinal)) +
  xlab("Cronan et al.\nM1 Module Score (Arb. Units)") + 
  ylab("Cronan et al.\nEpithelioid Module Score (Arb. Units)") + 
  geom_point(size = 3) + theme_classic() + geom_text_repel() + 
  scale_color_manual(values = mac.color) + 
  theme(text = element_text(size = font.size.in), legend.position = "none", axis.text = element_blank())
cor.test(mphage.module.plot$mean.Cronan.M11,
         mphage.module.plot$mean.Cronan.Epithelioid1,
         method = "spearman")

mphage.metadata.df <- MacrophageSeurat@meta.data[,c("SpecificFinal", "Granuloma", "GranulomaBurden", "logCFU")] %>%
  group_by(SpecificFinal, Granuloma, GranulomaBurden, logCFU) %>%
  summarise(cluster.count = n()) %>%
  ungroup() %>%
  group_by(Granuloma, GranulomaBurden, logCFU) %>%
  dplyr::mutate(cluster.prop = 100*cluster.count / sum(cluster.count)) %>%
  complete(SpecificFinal, Granuloma, fill = list(cluster.prop = 0, cluster.count = 0))

axis.break = 10
mac4.cluster.storage.SpecificFinal.CFU <- cor.abundance.burden(object.in = GranulomaTotalSeurat,
                                                          cluster.level = "SpecificFinal",
                                                          group.var = "Granuloma",
                                                          burden.var = "logCFU",
                                                          split.var = "GranulomaBurden")
mac4.mphage.metadata.df <- mac4.cluster.storage.SpecificFinal.CFU %>% dplyr::filter(cluster == "Mac4") %>%
  dplyr::mutate(axis.break.field = case_when(cluster.prop < axis.break ~ 1, cluster.prop > axis.break ~ 0))
p.mac4.cfu <- ggplot(data = mac4.mphage.metadata.df, aes(x = logCFU, y = cluster.prop, color = GranulomaBurden)) + 
  geom_point(size = 3) + 
  facet_grid(rows = vars(axis.break.field), scales = "free") + 
  force_panelsizes(rows = c(0.75, 1)) + 
  theme_classic() + 
  xlab(expression(~log[10]~CFU~per~Granuloma)) + 
  ylab("Percent of All Macrophages (%)") + 
  scale_color_manual(values = c("#f7931e", "#029151")) + 
  scale_x_continuous(expand = c(0,0)) + scale_y_continuous(expand = c(0,0)) + 
  coord_cartesian(clip = "off") + ggtitle("Mac4: Correlation with\nGranuloma CFU") + 
  theme(strip.text.y = element_blank(), legend.position = "none", 
        text = element_text(size = font.size.in), plot.title = element_text(hjust = 0.5))
cor.test(mac4.mphage.metadata.df$logCFU,
         mac4.mphage.metadata.df$cluster.prop,
         method = "spearman")

esaulova.mac.plotlist <- list(p.mac4.cfu, p.mac4.dotplot, NULL)
cronan.mac.plotlist <- list(p.cronanepi.vln, p.cronanepi.dot, p.cronan.epivsm1)
p.esaulova <- cowplot::plot_grid(plotlist = esaulova.mac.plotlist, ncol = 3, rel_widths = c(1,1,1), 
                   axis = "tblr", align = "hv", labels = c("A", "B", NA), label_size = 16)
p.cronan <- cowplot::plot_grid(plotlist = cronan.mac.plotlist, ncol = 3, rel_widths = c(1,1,1), 
                   axis = "tblr", align = "hv", labels = c("C", "D", "E"), label_size = 16)
all.plotlist <- list(p.mac4.cfu, p.mac4.dotplot, NULL, p.cronanepi.vln, p.cronanepi.dot, p.cronan.epivsm1)
pdf(paste0(figdir, "/", "EsaulovaCronan_MacrophageComparison.pdf"), height = 12, width = 16)
cowplot::plot_grid(plotlist = all.plotlist, nrow = 2, ncol = 3, labels = c("A", "B", NA, "C", "D", "E"), 
                   rel_heights = c(1,1), rel_widths = c(1,1.2,1), axis = "tblr", align = "hv")
dev.off()
```

```{r Gene shuffling to assess Mac5 significance}
seurat.in <- MacrophageSeurat
geneset.oi.in <- rhesus.epithelioid
n.rep <- 100
ii = 1
avgexpr.geneset.ctrl <- data.frame()
ctrl.gene.df <- matched.expression.ctrl.gene.module.CT(seurat.in = MacrophageSeurat, gene.list.in = geneset.oi.in, nbin.in = 24, ctrl.in = n.rep)
set.seed(54)
for(ii in 1:n.rep){
  if(ii %% 5 == 0){print(paste0("Done with ", ii, " iterations!"))}
  geneset.ctrl.ii <- ctrl.gene.df[ii,] %>% as.character()
  seurat.in <- AddModuleScore(object = seurat.in, features = list(geneset.ctrl.ii), name = "Control.Module.ii")
  avgexpr.geneset.ctrl.ii <- seurat.in@meta.data %>% dplyr::select(SpecificFinal, Control.Module.ii1) %>% 
    group_by(SpecificFinal) %>% summarise(meanModule = mean(Control.Module.ii1)) %>% dplyr::mutate(iter.count = ii)
  avgexpr.geneset.ctrl <- rbind(avgexpr.geneset.ctrl, avgexpr.geneset.ctrl.ii)
  seurat.in$Control.Module.ii1 <- NULL
}
avgexpr.geneset.real <- seurat.in@meta.data %>% dplyr::select(SpecificFinal, eval("Cronan.Epithelioid1")) %>% 
  dplyr::rename("ModuleScore" = 2) %>%
  group_by(SpecificFinal) %>% 
  summarise(RealDataMeanModule = mean(ModuleScore))

avgexpr.geneset.ctrl.summary <- avgexpr.geneset.ctrl %>%
  ungroup() %>%
  group_by(SpecificFinal) %>%
  summarise(meanmeanModule = mean(meanModule), sdmeanModule = sd(meanModule))

avgexpr.geneset <- inner_join(x = avgexpr.geneset.ctrl.summary, 
                              y = avgexpr.geneset.real,
                              by = "SpecificFinal")
avgexpr.geneset <- avgexpr.geneset %>%
  dplyr::mutate(z.score = (RealDataMeanModule - meanmeanModule) / sdmeanModule) %>%
  dplyr::mutate(pval = pnorm(z.score, lower.tail = FALSE)) %>%
  dplyr::mutate(p.adj = p.adjust(p = pval, method = "BH"))
```

```{r Re-examine stem-like object}
Idents(TcellSeuratUpdate) <- "SpecificFinal"
project.name = "Granuloma_TStemLike"
T_stemlike_figdir <- "/Users/ctzouanas/Documents/MIT/Shalek/TB_IntercellularInteraction/Revisions/UsefulFiles/Subclustering/T_StemLike"

T.stemlike.obj <- TcellSeuratUpdate[,TcellSeuratUpdate$SpecificFinal == "TStemLike"]
T.stemlike.obj = SCTransform(object = T.stemlike.obj)
T.stemlike.obj = NormalizeData(T.stemlike.obj, assay = "RNA")
T.stemlike.obj = RunPCA(object = T.stemlike.obj)
ElbowPlot(T.stemlike.obj)
T.stemlike.obj@misc$nPCs <- 1:13
T.stemlike.obj = FindNeighbors(object = T.stemlike.obj, dims = T.stemlike.obj@misc$nPCs, k.param=ceiling(0.5*sqrt(ncol(T.stemlike.obj))))
T.stemlike.obj = RunUMAP(srobj, T.stemlike.obj = T.stemlike.obj@misc$nPCs)

DimPlot(T.stemlike.obj)
gene.of.int <- c("CCR7", "TCF7", "IL7R", "TXNIP", 
                 "CD69", "ITGB1")

plot.list <- list()
for(ii in 1:length(gene.of.int)){
  gene.ii <- gene.of.int[[ii]]
  p.ii <- FeaturePlot(object = T.stemlike.obj, features = gene.ii, order = FALSE, cols = c("lightgrey", "#ad2524"), pt.size = 0.5) + 
    xlab("UMAP1") + ylab("UMAP2") + 
    labs(color = "Expression Level\nlog2(TP10k+1)") + 
    theme(legend.position = "bottom", axis.text = element_blank(), legend.justification = "center", legend.title = element_text(hjust = 0.5))
  plot.list[[ii]] <- p.ii
}

pdf(file = paste0(T_stemlike_figdir, "/", "TStemLike_Markers.pdf"), width = 12, height = 8)
cowplot::plot_grid(plotlist = plot.list, nrow = 1, align = "hv", axis = "tblr")
dev.off()
```

```{r Generate main text Figure 5}
if(!exists("Granuloma4WeekSeurat")){
  load("/Users/ctzouanas/Dropbox (MIT)/Reanalysis2022/Granuloma4WeekSeurat.Rdata")
  Granuloma4WeekSeurat <- UpdateSeuratObject(object = Granuloma4WeekSeurat)
}

figdir.preprocess <- "/Users/ctzouanas/Documents/MIT/Shalek/TB_IntercellularInteraction/Revisions/UsefulFiles/4WeekAnalyses"

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)

font.size.4wk = 16

BurdenScoring <- function(obj.4wk = Granuloma4WeekSeurat,
                          obj.10wk = GranulomaTotalSeurat,
                          ident.4wk = "Round1CellTypes.plotting",
                          ident.10wk = "GenericFinal",
                          celltype.4wk,
                          celltype.10wk,
                          assay.in = "RNA",
                          palette.10wkvs4wk = c("#056c9a", "#ae2425", "#cc6328"),
                          min.cells = 10,
                          min.pct.in = 0.4,
                          logfc.thresh = 0.25){
  DefaultAssay(obj.4wk) <- assay.in
  DefaultAssay(obj.10wk) <- assay.in
  
  Idents(obj.10wk) <- ident.10wk
  celltype.10wk.obj <- obj.10wk[,obj.10wk[[ident.10wk]] == celltype.10wk]
  celltype.4wk.obj <- obj.4wk[,obj.4wk[[ident.4wk]] == celltype.4wk]
  celltype.10wk.highvslow <- FindMarkers(object = celltype.10wk.obj,
                                         ident.1 = "HighBurden", ident.2 = "LowBurden", group.by = "GranulomaBurden",
                                         min.pct = min.pct.in, logfc.threshold = logfc.thresh) %>% arrange(-avg_log2FC)
  
  celltype.10wk.upmodule <- celltype.10wk.highvslow %>% dplyr::filter(avg_log2FC > logfc.thresh) %>% rownames_to_column("gene") %>% dplyr::pull(gene)
  celltype.10wk.downmodule <- celltype.10wk.highvslow %>% dplyr::filter(avg_log2FC < -1*logfc.thresh) %>% rownames_to_column("gene") %>% dplyr::pull(gene)
  celltype.10wk.module <- unique(c(celltype.10wk.upmodule, celltype.10wk.downmodule))
  celltype.10wk.module <- celltype.10wk.module[celltype.10wk.module %in% rownames(celltype.4wk.obj)]
  
  celltype.4wk.obj <- AddModuleScore(object = celltype.4wk.obj, features = list(celltype.10wk.upmodule), name = "CellType.10Wk.UpModule")
  celltype.4wk.obj <- AddModuleScore(object = celltype.4wk.obj, features = list(celltype.10wk.downmodule), name = "CellType.10Wk.DownModule")
  celltype.4wk.obj$CellType.10Wk.DiffModule1 <- celltype.4wk.obj$CellType.10Wk.UpModule1 - celltype.4wk.obj$CellType.10Wk.DownModule1
  
  celltype.10wk.obj <- AddModuleScore(object = celltype.10wk.obj, features = list(celltype.10wk.upmodule), name = "CellType.10Wk.UpModule")
  celltype.10wk.obj <- AddModuleScore(object = celltype.10wk.obj, features = list(celltype.10wk.downmodule), name = "CellType.10Wk.DownModule")
  celltype.10wk.obj$CellType.10Wk.DiffModule1 <- celltype.10wk.obj$CellType.10Wk.UpModule1 - celltype.10wk.obj$CellType.10Wk.DownModule1
  
  summary.10wk.df <- celltype.10wk.obj@meta.data[,c("CellType.10Wk.DiffModule1", "GranulomaBurden", "Granuloma")] %>%
    group_by(Granuloma, GranulomaBurden) %>% 
    dplyr::summarise(medianModule = median(CellType.10Wk.DiffModule1),
                     cellcount = n())
  
  summary.4wk.df <- celltype.4wk.obj@meta.data[,c("CellType.10Wk.DiffModule1", "Granuloma")] %>%
    group_by(Granuloma) %>% 
    dplyr::summarise(medianModule = median(CellType.10Wk.DiffModule1),
                     cellcount = n()) %>%
    dplyr::mutate(GranulomaBurden = "Granuloma4Week")
  
  summary.df <- rbind(summary.10wk.df, summary.4wk.df) %>%
    dplyr::mutate(GranulomaBurden = factor(GranulomaBurden,
                                           levels = c("LowBurden", "HighBurden", "Granuloma4Week"),
                                           labels = c("10 Week:\nLow Burden", "10 Week:\nHigh Burden", "4 Week"))) %>%
    dplyr::filter(cellcount >= min.cells)
  
  scores.10wkhigh <- summary.df %>% dplyr::filter(GranulomaBurden == "10 Week:\nHigh Burden") %>% pull(medianModule)
  scores.10wklow <- summary.df %>% dplyr::filter(GranulomaBurden == "10 Week:\nLow Burden") %>% pull(medianModule)
  scores.4wk <- summary.df %>% dplyr::filter(GranulomaBurden == "4 Week") %>% pull(medianModule)
  
  print("Cohen's D for 10-week high vs. 10-week low")
  out <- cohensD(x = scores.10wkhigh, y = scores.10wklow)
  print(out)
  print("Cohen's D for 10-week high vs. 4 week")
  out <- cohensD(x = scores.10wkhigh, y = scores.4wk)
  print(out)
  print("Cohen's D for 10-week low vs. 4 week")
  out <- cohensD(x = scores.10wklow, y = scores.4wk)
  print(out)
  print("P Values")
  pval.df <- data.frame(test = c("10wkHighvs10wkLow", "10wkHighvs4wk", "10wkLowvs4wk"),
                        pval = c(wilcox.test(x = scores.10wkhigh, y = scores.10wklow)$p.value,
                                 wilcox.test(x = scores.10wkhigh, y = scores.4wk)$p.value,
                                 wilcox.test(x = scores.10wklow, y = scores.4wk)$p.value))
  pval.df$padj <- p.adjust(p = pval.df$pval, method = "BH")
  print(pval.df$test)
  print(pval.df$padj)
  p.celltype.modulescore <- summary.df %>% ggplot(aes(x = GranulomaBurden, y = medianModule, color = GranulomaBurden)) +
    geom_boxplot(lwd = 1) + geom_jitter(size = 2, width = 0.1) + theme_classic() + 
    scale_color_manual(values = palette.10wkvs4wk) + 
    xlab("GranulomaType") + ylab("Burden-Associated Module Score (A.U.)") + 
    ggtitle(paste0(celltype.4wk, "\nBurden-Associated Gene Set")) + 
    theme(legend.position = "none", plot.title = element_text(hjust = 0.5), text = element_text(size = font.size.4wk))
  plot(p.celltype.modulescore)
  return(p.celltype.modulescore)
}

cluster.mapping <- c("T" = "T",
                     "Macrophage" = "Mphage")

p.T.modulescore <- BurdenScoring(obj.4wk = Granuloma4WeekSeurat,
                                 obj.10wk = GranulomaTotalSeurat,
                                 ident.4wk = "Round1CellTypes.plotting",
                                 ident.10wk = "GenericFinal",
                                 celltype.4wk = names(cluster.mapping)[[1]],
                                 celltype.10wk = cluster.mapping[[1]],
                                 assay.in = "RNA",
                                 palette.10wkvs4wk = c("#056c9a", "#ae2425", "#cc6328"),
                                 min.cells = 10,
                                 min.pct.in = 0.4,
                                 logfc.thresh = 0.25)
p.macro.modulescore <- BurdenScoring(obj.4wk = Granuloma4WeekSeurat,
                                     obj.10wk = GranulomaTotalSeurat,
                                     ident.4wk = "Round1CellTypes.plotting",
                                     ident.10wk = "GenericFinal",
                                     celltype.4wk = names(cluster.mapping)[[2]],
                                     celltype.10wk = cluster.mapping[[2]],
                                     assay.in = "RNA",
                                     palette.10wkvs4wk = c("#056c9a", "#ae2425", "#cc6328"),
                                     min.cells = 10,
                                     min.pct.in = 0.4,
                                     logfc.thresh = 0.25)

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) %>%
  rename(Granuloma = 1) %>% 
  dplyr::mutate(GranulomaBurden = "Granuloma4Week")

Granuloma10Week.Metadata <- CFUDataFinal_Update %>% dplyr::select(Granuloma, logCFU, HighLowUpdate) %>%
  dplyr::mutate(GranulomaBurden = case_when(HighLowUpdate == 1 ~ "HighBurden", HighLowUpdate == 0 ~ "LowBurden")) %>%
  dplyr::select(!HighLowUpdate)

Granuloma4and10Week.Metadata <- rbind(Granuloma4Week.Metadata, Granuloma10Week.Metadata) %>%
  dplyr::mutate(GranulomaBurden = factor(GranulomaBurden,
                                         levels = c("LowBurden", "HighBurden", "Granuloma4Week"),
                                         labels = c("10 Week:\nLow Burden", "10 Week:\nHigh Burden", "4 Week")))
palette.10wkvs4wk = c("#056c9a", "#ae2425", "#cc6328")
p.granuloma.cfu <- Granuloma4and10Week.Metadata %>% ggplot(aes(x = GranulomaBurden, y = logCFU, color = GranulomaBurden)) +
    geom_boxplot(lwd = 1) + geom_jitter(size = 2, width = 0.1) + theme_classic() + 
    scale_color_manual(values = palette.10wkvs4wk) + 
    xlab("GranulomaType") + ylab(expression(~log[10]~CFU~per~Granuloma)) + 
    ggtitle(paste0("Bacterial Burden\nAcross Granulomas")) + 
    theme(legend.position = "none", plot.title = element_text(hjust = 0.5), text = element_text(size = font.size.4wk))
cfu.10wkhigh <- Granuloma4and10Week.Metadata %>% dplyr::filter(GranulomaBurden == "10 Week:\nHigh Burden") %>% pull(logCFU)
cfu.10wklow <- Granuloma4and10Week.Metadata %>% dplyr::filter(GranulomaBurden == "10 Week:\nLow Burden") %>% pull(logCFU)
cfu.4wk <- Granuloma4and10Week.Metadata %>% dplyr::filter(GranulomaBurden == "4 Week") %>% pull(logCFU)
pval.df <- data.frame(test = c("10wkHighvs10wkLow", "10wkHighvs4wk", "10wkLowvs4wk"),
                      pval = c(wilcox.test(x = cfu.10wkhigh, y = cfu.10wklow)$p.value,
                               wilcox.test(x = cfu.10wkhigh, y = cfu.4wk)$p.value,
                               wilcox.test(x = cfu.10wklow, y = cfu.4wk)$p.value))
pval.df$padj <- p.adjust(p = pval.df$pval, method = "BH")

Idents(Granuloma4WeekSeurat) <- "Round1CellTypes.plotting"
markers.chosen <- c("MS4A1", "CD79B", "CD3D", "ETS1", "CD160", "JCHAIN", "MPEG1", "CSF1R", "CPA3", "GATA2", "CXCL8", "CSF3R", "LILRA4", "TCF4", "VWF", "FCN3", "COL3A1", "DCN", "SFTPC", "SFTPB", "SCGB3A1", "PIGR") %>% rev()

p.4wk.dotplot <- DotPlot(Granuloma4WeekSeurat, features = markers.chosen, group.by = "Round1CellTypes.plotting", 
                          assay = "RNA", scale = TRUE, cols = c("lightgrey", "#ad2524")) + 
  ggtitle("4-Week Granuloma scRNA-seq Markers") + 
  theme(legend.position = "bottom", legend.justification = "center", legend.margin=margin(), legend.box = "vertical", 
        plot.title = element_text(face = "plain", hjust = 0.5),
        axis.title = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1), text = element_text(size = font.size.4wk)) + 
  coord_flip() + 
  scale_colour_gradient2(low = "#125175", mid = "white", high = "#ad2524")
palette.4wk <- c("B" = "#584e9d",
                 "T" = "#377bb4",
                 "Plasma" = "#4da9a6",
                 "Macrophage" = "#84c69f",
                 "Mast" = "#bbd89e",
                 "Neutrophil" = "#e2e6a0",
                 "pDC" = "#f6e492",
                 "Endothelial" = "#f08d55",
                 "Fibroblast" = "#e15d49", 
                 "T2P" = "#c6374f",
                 "Club" = "#951f44")
p.4wk.umap <- DimPlot(Granuloma4WeekSeurat, group.by = "Round1CellTypes.plotting", label = TRUE, cols = palette.4wk) +
  xlab("UMAP1") + ylab("UMAP2") + ggtitle("4-Week Granuloma scRNA-seq") + 
  theme(axis.text = element_blank(), plot.title = element_text(face = "plain"), text = element_text(size = font.size.4wk),
        legend.position = "bottom", legend.justification = "center", legend.box = "vertical")

row1.4wk <- cowplot::plot_grid(NULL, p.granuloma.cfu, p.4wk.umap, ncol = 3, align = "hv", axis = "tblr", labels = c("A", "B", "C"), rel_widths = c(1,1,1.75), label_size = 24)
row2.4wk <- cowplot::plot_grid(p.4wk.dotplot, p.macro.modulescore, p.T.modulescore, ncol = 3, rel_widths = c(1.5,1,1), align = "hv", axis = "tblr", labels = c("D", "E", "F"), label_size = 24)
combined.4wk <- cowplot::plot_grid(row1.4wk, row2.4wk,
                                   nrow = 2, rel_heights = c(1, 1), align = "hv", axis = "tblr")
pdf(paste0(figdir, "/", "GranulomaPhenotype_Necrotic_4Week.pdf"), height = 16, width = 16, useDingbats = FALSE)
plot(combined.4wk)
dev.off()

p.macro.modulescore
```
