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

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 <- "insert-figure-directory-here"
```

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

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

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

```{r Generate composition heatmaps}
celltype.sample.df <- GranulomaTotalSeurat@meta.data[,c("SpecificFinal", "Granuloma")]
granuloma.celltype.abundance.df <- celltype.sample.df %>% group_by(SpecificFinal, Granuloma) %>% summarise(celltypecount = n()) %>% ungroup() %>%
  group_by(Granuloma) %>% mutate(celltypeprop = celltypecount / sum(celltypecount)) %>%
  dplyr::select(!celltypecount) %>% 
  pivot_wider(names_from = Granuloma, values_from = celltypeprop) %>% column_to_rownames("SpecificFinal")
granuloma.celltype.abundance.df[is.na(granuloma.celltype.abundance.df)] <- 0
SpecificCellTypeCor_Pearson <- cor(t(granuloma.celltype.abundance.df), method = "pearson")

correlation.heatmap <- pheatmap(SpecificCellTypeCor_Pearson)
Groups <- sort(cutree(correlation.heatmap$tree_row, k=5))

group.anno.col <- data.frame(Group = Groups[colnames(SpecificCellTypeCor_Pearson)]) %>% 
  dplyr::mutate(Group = factor(as.character(Group), levels = c("1", "2", "3", "4", "5")))
group.anno.col.colors <- c("1" = "#faaa5e", "2" = "#ca5e4a", "3" = "#4e88b9", "4" = "#6ca659", "5" = "#4e4ca0")
# pheatmap's internal numbering doesn't go 1 to 5 from left to right on the heatmap, so have to do some weird shuffling
group.anno.col.colors <- c("1" = "#4e4ca0", "2" = "#ca5e4a", "3" = "#faaa5e", "4" = "#4e88b9", "5" = "#6ca659")

paletteLength <- 50
myColor <- colorRampPalette(c("#6d8fb9", "white", "#cb604c"))(paletteLength)
# use floor and ceiling to deal with even/odd length pallettelengths
myBreaks <- c(seq(min(SpecificCellTypeCor_Pearson), 0, length.out=ceiling(paletteLength/2) + 1), 
              seq(max(SpecificCellTypeCor_Pearson)/paletteLength, max(SpecificCellTypeCor_Pearson), length.out=floor(paletteLength/2)))

correlation.heatmap <- pheatmap(SpecificCellTypeCor_Pearson, color = myColor, breaks = myBreaks,
                 filename = paste0(figdir, "/", "CorrelationHeatmap_WithLegend.pdf"),
                 cutree_rows = 5,cutree_cols = 5, border_color = NA,
                 annotation_col = group.anno.col, treeheight_col = 0,
                 annotation_colors = list(Group = group.anno.col.colors), show_colnames = FALSE)
correlation.heatmap <- pheatmap(SpecificCellTypeCor_Pearson, color = myColor, breaks = myBreaks,
                                filename = paste0(figdir, "/", "CorrelationHeatmap_NoLegend.pdf"),
                                cutree_rows = 5,cutree_cols = 5, border_color = NA, legend = FALSE, annotation_legend = FALSE,
                                annotation_col = group.anno.col, treeheight_col = 0,
                                annotation_colors = list(Group = group.anno.col.colors), show_colnames = FALSE)

granuloma.by.celltypegroups <- GranulomaTotalSeurat@meta.data[,c("Granuloma", "CellTypeGroup", "GranulomaBurden")] %>% 
  group_by(Granuloma, CellTypeGroup, GranulomaBurden) %>% summarise(cellcount = n()) %>% ungroup() %>%
  group_by(Granuloma) %>% mutate(cellprop = cellcount / sum(cellcount)) %>% ungroup()
granuloma.by.celltypegroups <- granuloma.by.celltypegroups %>% 
  dplyr::mutate(GranulomaLabel = as.numeric(gsub(pattern = "Granuloma", replacement = "", x = Granuloma)))
granuloma.by.celltypegroups <- granuloma.by.celltypegroups %>% arrange(GranulomaLabel) %>% arrange(GranulomaBurden) %>% 
  dplyr::mutate(GranulomaLabel = as.character(GranulomaLabel)) %>% 
  dplyr::mutate(GranulomaLabel = factor(GranulomaLabel, levels = unique(GranulomaLabel)))

allgranuloma.barplot <- ggplot(granuloma.by.celltypegroups, aes(fill = CellTypeGroup, y=cellcount, x = GranulomaLabel)) + 
  geom_bar(position = "fill", stat = "identity") + theme_classic() + 
  scale_y_continuous(expand = c(0,0), breaks = c(0, 1)) + 
  scale_fill_manual(values = c("#faaa5e", "#ca5e4a", "#4e88b9", "#6ca659", "#4e4ca0")) + 
  ylab("Group Proportion") + xlab("Granuloma") + 
  theme(text = element_text(size = 16), axis.text.x = element_blank(), 
        axis.line.x = element_blank(), axis.ticks.x = element_blank(), legend.position = "none", axis.line.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank(), axis.text.y = element_blank())

granuloma.metadata.df <- CFUDataFinal_Update %>% 
  dplyr::mutate(Granuloma = factor(Granuloma, levels = paste0("Granuloma", levels(granuloma.by.celltypegroups$GranulomaLabel)))) %>%
  arrange(Granuloma) %>% dplyr::select(Granuloma, HighLow, Animal)
colnames(granuloma.metadata.df) <- c("Granuloma", "Burden", "Animal")
granuloma.metadata.df <- granuloma.metadata.df %>% pivot_longer(!Granuloma, names_to = "variable", values_to = "value") %>% 
  dplyr::mutate(value = as.character(value))  %>%
  dplyr::mutate(GranulomaLabel = gsub(pattern = "Granuloma", replacement = "", x = Granuloma)) %>%
  dplyr::mutate(GranulomaLabel = factor(as.character(GranulomaLabel), levels = unique(as.character(GranulomaLabel))))

granuloma.metadata.plot <- ggplot(granuloma.metadata.df, aes(x = GranulomaLabel, y = variable)) + 
  geom_tile(aes(fill = value), color = "white", size = 0.5) + theme_classic() + 
  theme(text = element_text(size = 16)) + 
  scale_fill_manual(values = c("#046c9a", "#ad2524", "#ffaa7f", "#abcfff", "#dd775a", "#c4bedf")) + 
  scale_y_discrete(position = "right") + 
  theme(axis.title.y = element_blank(), legend.position = "none", axis.line.y = element_blank(), axis.ticks.y = element_blank(),
        axis.title.x = element_blank(), axis.text.x = element_blank(), axis.line.x = element_blank(), axis.ticks.x = element_blank())

highvslow.cellcounts.df <- GranulomaTotalSeurat@meta.data[,c("CellTypeGroup", "GranulomaBurden")] %>% 
  group_by(CellTypeGroup, GranulomaBurden) %>% summarise(cellcount = n()) %>% ungroup() %>%
  group_by(GranulomaBurden) %>% mutate(cellprop = cellcount / sum(cellcount)) %>% ungroup() %>% 
  dplyr::mutate(GranulomaBurdenLabel = gsub(pattern = "Burden", replacement = "", x = GranulomaBurden))
granuloma.summary.plot <- ggplot(highvslow.cellcounts.df, aes(fill = CellTypeGroup, y=cellcount, x = GranulomaBurdenLabel)) + 
  geom_bar(position = "fill", stat = "identity") + theme_classic() + 
  scale_y_continuous(expand = c(0,0), breaks = c(0, 1)) + 
  scale_fill_manual(values = c("#faaa5e", "#ca5e4a", "#4e88b9", "#6ca659", "#4e4ca0")) + 
  ylab("Group Proportion") + 
  theme(text = element_text(size = 16), legend.position = "none", axis.title.x = element_blank(), axis.line.x = element_blank(), axis.ticks.x = element_blank(),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

plotlist <- list(correlation.heatmap[[4]], granuloma.summary.plot + theme(axis.text.x = element_blank()), allgranuloma.barplot + theme(axis.title.x = element_blank()),
                 NULL, NULL, granuloma.metadata.plot)
compiled.plot <- plot_grid(plotlist = plotlist, align = "hv", axis = "tblr", ncol = 3, 
                           rel_widths = c(3, 1, 3), rel_heights = c(4,1))
pdf(paste0(figdir, "/", "CompositionBarplots_and_CorrHeatmap_NoLegend.pdf"), width = 12, height = 6, useDingbats = FALSE)
plot(compiled.plot)
dev.off()

plotlist <- list(correlation.heatmap[[4]], granuloma.summary.plot, allgranuloma.barplot,
                 NULL, NULL, granuloma.metadata.plot)
compiled.plot <- plot_grid(plotlist = plotlist, align = "hv", axis = "tblr", ncol = 3, 
                           rel_widths = c(3, 1, 3), rel_heights = c(4,1))
pdf(paste0(figdir, "/", "CompositionBarplots_and_CorrHeatmap_WithLegend.pdf"), width = 12, height = 6, useDingbats = FALSE)
plot(compiled.plot)
dev.off()
 
group2.highburden.props <- granuloma.by.celltypegroups %>% 
  dplyr::filter(GranulomaBurden == "HighBurden" & CellTypeGroup == "Group2") %>% 
  dplyr::select(cellprop) %>% pull()
group2.lowburden.props <- granuloma.by.celltypegroups %>% 
  dplyr::filter(GranulomaBurden == "LowBurden" & CellTypeGroup == "Group2") %>% 
  dplyr::select(cellprop) %>% pull()
pval.group2.highvslow <- wilcox.test(x = group2.highburden.props,
                                     y = group2.lowburden.props)$p.value

group3.highburden.props <- granuloma.by.celltypegroups %>% 
  dplyr::filter(GranulomaBurden == "HighBurden" & CellTypeGroup == "Group3") %>% 
  dplyr::select(cellprop) %>% pull()
group3.lowburden.props <- granuloma.by.celltypegroups %>% 
  dplyr::filter(GranulomaBurden == "LowBurden" & CellTypeGroup == "Group3") %>% 
  dplyr::select(cellprop) %>% pull()
pval.group3.highvslow <- wilcox.test(x = group3.highburden.props,
                                     y = group3.lowburden.props)$p.value
```

```{r Analyses to look at more global signaling patterns}
travis.rdata.path <- "/Users/ctzouanas/Dropbox (MIT)/ConstantineData/Interactions/InteractionsMergeTotal.Rdata"

if(!exists("InteractionsMergeTotal")){
  load(travis.rdata.path)
  InteractionsMergeTotal <- InteractionsMergeTotal %>% 
    dplyr::select(c(RLTypes, HighLowPVal, LowPermPvals, ComboReceptorMeans_High, ComboReceptorMeans_Low)) %>%
    dplyr::rename(ComboReceptorHighLowPVal = HighLowPVal) %>%
    dplyr::mutate(ComboReceptorHighLowFC = ComboReceptorMeans_High / ComboReceptorMeans_Low) %>% 
    tidyr::replace_na(list(ComboReceptorHighLowPVal = 1, ComboReceptorHighLowFC = 0))
  interact.df <- InteractionsMergeTotal %>% dplyr::filter(ComboReceptorHighLowPVal < 0.05 & LowPermPvals < 0.05)
  interact.df$ComboReceptorHighLowFC[interact.df$ComboReceptorHighLowFC == Inf] = 
    max(interact.df$ComboReceptorHighLowFC[interact.df$ComboReceptorHighLowFC != Inf])
  interact.df$ComboReceptorHighLowFC[interact.df$ComboReceptorHighLowFC == 0] = 
    min(interact.df$ComboReceptorHighLowFC[interact.df$ComboReceptorHighLowFC != 0])
  interact.df <- interact.df %>% tidyr::separate(col = RLTypes, into = c("Sender", "Receiver", "Ligand", "Receptor"), sep = "_") %>%
    dplyr::mutate(Sender_Ligand = paste0(Sender, "_", Ligand), 
                  Receiver_Receptor = paste0(Receiver, "_", Receptor))
}

poslogFC.cellcell.count <- interact.df %>% ungroup() %>% dplyr::filter(ComboReceptorHighLowFC > 1) %>% group_by(Sender, Receiver) %>% summarise(count.interacts = n())
neglogFC.cellcell.count <- interact.df %>% ungroup() %>% dplyr::filter(ComboReceptorHighLowFC < 1) %>% group_by(Sender, Receiver) %>% summarise(count.interacts = n())

thresh.method <- "proportion"
if(thresh.method == "proportion"){
  pos.min.thresh <- 0.00*sum(poslogFC.cellcell.count$count.interacts)
  neg.min.thresh <- 0.00*sum(neglogFC.cellcell.count$count.interacts)
} else if(thresh.method == "absolute"){
  pos.min.thresh <- 4
  neg.min.thresh <- 4
}

poslogFC.cellcell.count <- poslogFC.cellcell.count %>% ungroup() %>% dplyr::filter(count.interacts > pos.min.thresh) %>% dplyr::mutate(GranulomaBurden = "HighBurden")
neglogFC.cellcell.count <- neglogFC.cellcell.count %>% ungroup() %>% dplyr::filter(count.interacts > neg.min.thresh) %>% dplyr::mutate(GranulomaBurden = "LowBurden")

dataframe.choice <- rbind(poslogFC.cellcell.count, neglogFC.cellcell.count)
dataframe.choice <- dataframe.choice[,c("Receiver", "Sender", "count.interacts", "GranulomaBurden")]

dataframe.choice$Sender <- paste0("S.", dataframe.choice$Sender)
dataframe.choice$Receiver <- paste0("R.", dataframe.choice$Receiver)

dataframe.choice$ReceiverGroup <- ""
dataframe.choice$ReceiverGroup[dataframe.choice$Receiver %in% paste0("R.", group1)] <- "Group1"
dataframe.choice$ReceiverGroup[dataframe.choice$Receiver %in% paste0("R.", group2)] <- "Group2"
dataframe.choice$ReceiverGroup[dataframe.choice$Receiver %in% paste0("R.", group3)] <- "Group3"
dataframe.choice$ReceiverGroup[dataframe.choice$Receiver %in% paste0("R.", group4)] <- "Group4"
dataframe.choice$ReceiverGroup[dataframe.choice$Receiver %in% paste0("R.", group5)] <- "Group5"
dataframe.choice$SenderGroup <- ""
dataframe.choice$SenderGroup[dataframe.choice$Sender %in% paste0("S.", group1)] <- "Group1"
dataframe.choice$SenderGroup[dataframe.choice$Sender %in% paste0("S.", group2)] <- "Group2"
dataframe.choice$SenderGroup[dataframe.choice$Sender %in% paste0("S.", group3)] <- "Group3"
dataframe.choice$SenderGroup[dataframe.choice$Sender %in% paste0("S.", group4)] <- "Group4"
dataframe.choice$SenderGroup[dataframe.choice$Sender %in% paste0("S.", group5)] <- "Group5"

highburden.groupsummarize <- dataframe.choice %>% dplyr::filter(GranulomaBurden == "HighBurden") %>% ungroup() %>% 
  group_by(ReceiverGroup, SenderGroup) %>% summarise(group.count.interacts = sum(count.interacts))
highburden.groupsummarize$ReceiverGroup <- paste0("R.", highburden.groupsummarize$ReceiverGroup)
highburden.groupsummarize$SenderGroup <- paste0("S.", highburden.groupsummarize$SenderGroup)
highburden.groupsummarize <- highburden.groupsummarize %>% pivot_wider(names_from = ReceiverGroup, values_from = group.count.interacts) %>% 
  pivot_longer(!SenderGroup, names_to = "ReceiverGroup", values_to = "group.count.interacts")
colnames(highburden.groupsummarize) <- c("SenderGroup", "ReceiverGroup", "interact.counts.highburden")

lowburden.groupsummarize <- dataframe.choice %>% dplyr::filter(GranulomaBurden == "LowBurden") %>% ungroup() %>% 
  group_by(ReceiverGroup, SenderGroup) %>% summarise(group.count.interacts = sum(count.interacts))
lowburden.groupsummarize$ReceiverGroup <- paste0("R.", lowburden.groupsummarize$ReceiverGroup)
lowburden.groupsummarize$SenderGroup <- paste0("S.", lowburden.groupsummarize$SenderGroup)
lowburden.groupsummarize <- lowburden.groupsummarize %>% pivot_wider(names_from = ReceiverGroup, values_from = group.count.interacts) %>% 
  pivot_longer(!SenderGroup, names_to = "ReceiverGroup", values_to = "group.count.interacts")
colnames(lowburden.groupsummarize) <- c("SenderGroup", "ReceiverGroup", "interact.counts.lowburden")

allburden.groupsummarize <- full_join(x = highburden.groupsummarize, y = lowburden.groupsummarize, by = c("SenderGroup", "ReceiverGroup"))
allburden.groupsummarize[is.na(allburden.groupsummarize)] <- 0

allburden.table <- as.table(rbind(
  allburden.groupsummarize$interact.counts.highburden,
  allburden.groupsummarize$interact.counts.lowburden
))
dimnames(allburden.table) <- list(
  Burden = c("HighBurden", "LowBurden"),
  Interaction = paste0(allburden.groupsummarize$SenderGroup, "__", allburden.groupsummarize$ReceiverGroup)
)

## Summarise intercellular interaction topologies
fisher_test(allburden.table, detailed = TRUE, simulate.p.value = TRUE)
chisq.test(allburden.table)

highburden.summarizesenders <- allburden.groupsummarize %>% group_by(SenderGroup) %>% summarise(sum.sender.interacts = sum(interact.counts.highburden)) %>% 
  dplyr::mutate(proportion.sender.interacts = sum.sender.interacts/sum(sum.sender.interacts))
lowburden.summarizesenders <- allburden.groupsummarize %>% group_by(SenderGroup) %>% summarise(sum.sender.interacts = sum(interact.counts.lowburden)) %>% 
  dplyr::mutate(proportion.sender.interacts = sum.sender.interacts/sum(sum.sender.interacts))
allburden.groupsummarize <- allburden.groupsummarize %>% ungroup() %>% group_by(ReceiverGroup) %>% 
  dplyr::mutate(highburden.group.sum.interact.counts = sum(interact.counts.highburden), lowburden.group.sum.interact.counts = sum(interact.counts.lowburden)) %>% 
  dplyr::mutate(highburden.proportion.signals.received = interact.counts.highburden / highburden.group.sum.interact.counts, 
                lowburden.proportion.signals.received = interact.counts.lowburden / lowburden.group.sum.interact.counts)
```

```{r Fuller circos plots but overlaying all interactions onto same plot}
travis.rdata.path <- "/Users/ctzouanas/Dropbox (MIT)/ConstantineData/Interactions/InteractionsMergeTotal.Rdata"
if(!exists("InteractionsMergeTotal")){
  load(travis.rdata.path)
  InteractionsMergeTotal <- InteractionsMergeTotal %>% 
    dplyr::select(c(RLTypes, HighLowPVal, LowPermPvals, ComboReceptorMeans_High, ComboReceptorMeans_Low)) %>%
    dplyr::rename(ComboReceptorHighLowPVal = HighLowPVal) %>%
    dplyr::mutate(ComboReceptorHighLowFC = ComboReceptorMeans_High / ComboReceptorMeans_Low) %>% 
    tidyr::replace_na(list(ComboReceptorHighLowPVal = 1, ComboReceptorHighLowFC = 0)) %>%
    dplyr::mutate(ComboReceptorHighLowFC = replace(ComboReceptorHighLowFC, is.infinite(ComboReceptorHighLowFC), 1.5))
  interact.df <- InteractionsMergeTotal %>% dplyr::filter(ComboReceptorHighLowPVal < 0.05 & LowPermPvals < 0.05) %>%
    tidyr::separate(col = RLTypes, into = c("Sender", "Receiver", "Ligand", "Receptor"), sep = "_") %>%
    dplyr::mutate(Sender_Ligand = paste0(Sender, "_", Ligand), 
                  Receiver_Receptor = paste0(Receiver, "_", Receptor))
}

poslogFC.cellcell.count <- interact.df %>% ungroup() %>% dplyr::filter(ComboReceptorHighLowFC > 1) %>% group_by(Sender, Receiver) %>% summarise(count.interacts = n())
neglogFC.cellcell.count <- interact.df %>% ungroup() %>% dplyr::filter(ComboReceptorHighLowFC < 1) %>% group_by(Sender, Receiver) %>% summarise(count.interacts = n())

thresh.method <- "proportion"
thresh.val <- 0.00
if(thresh.method == "proportion"){
  pos.min.thresh <- thresh.val*sum(poslogFC.cellcell.count$count.interacts)
  neg.min.thresh <- thresh.val*sum(neglogFC.cellcell.count$count.interacts)
} else if(thresh.method == "absolute"){
  pos.min.thresh <- 4
  neg.min.thresh <- 4
}

poslogFC.cellcell.count <- poslogFC.cellcell.count %>% ungroup() %>% dplyr::filter(count.interacts > pos.min.thresh) %>% dplyr::mutate(GranulomaBurden = "HighBurden")
neglogFC.cellcell.count <- neglogFC.cellcell.count %>% ungroup() %>% dplyr::filter(count.interacts > neg.min.thresh) %>% dplyr::mutate(GranulomaBurden = "LowBurden")

dataframe.choice <- rbind(poslogFC.cellcell.count, neglogFC.cellcell.count)
dataframe.choice <- dataframe.choice[,c("Receiver", "Sender", "count.interacts", "GranulomaBurden")]

dataframe.choice$Sender <- paste0("S.", dataframe.choice$Sender)
dataframe.choice$Receiver <- paste0("R.", dataframe.choice$Receiver)

dataframe.choice$ReceiverGroup <- ""
dataframe.choice$ReceiverGroup[dataframe.choice$Receiver %in% paste0("R.", group1)] <- "Group1"
dataframe.choice$ReceiverGroup[dataframe.choice$Receiver %in% paste0("R.", group2)] <- "Group2"
dataframe.choice$ReceiverGroup[dataframe.choice$Receiver %in% paste0("R.", group3)] <- "Group3"
dataframe.choice$ReceiverGroup[dataframe.choice$Receiver %in% paste0("R.", group4)] <- "Group4"
dataframe.choice$ReceiverGroup[dataframe.choice$Receiver %in% paste0("R.", group5)] <- "Group5"
dataframe.choice$SenderGroup <- ""
dataframe.choice$SenderGroup[dataframe.choice$Sender %in% paste0("S.", group1)] <- "Group1"
dataframe.choice$SenderGroup[dataframe.choice$Sender %in% paste0("S.", group2)] <- "Group2"
dataframe.choice$SenderGroup[dataframe.choice$Sender %in% paste0("S.", group3)] <- "Group3"
dataframe.choice$SenderGroup[dataframe.choice$Sender %in% paste0("S.", group4)] <- "Group4"
dataframe.choice$SenderGroup[dataframe.choice$Sender %in% paste0("S.", group5)] <- "Group5"

highburden.groupsummarize <- dataframe.choice %>% dplyr::filter(GranulomaBurden == "HighBurden") %>% ungroup() %>% 
  group_by(ReceiverGroup, SenderGroup) %>% summarise(group.count.interacts = sum(count.interacts)) %>% dplyr::mutate(GranulomaBurden = "HighBurden")
highburden.groupsummarize$ReceiverGroup <- paste0("R.", highburden.groupsummarize$ReceiverGroup)
highburden.groupsummarize$SenderGroup <- paste0("S.", highburden.groupsummarize$SenderGroup)

lowburden.groupsummarize <- dataframe.choice %>% dplyr::filter(GranulomaBurden == "LowBurden") %>% ungroup() %>% 
  group_by(ReceiverGroup, SenderGroup) %>% summarise(group.count.interacts = sum(count.interacts)) %>% dplyr::mutate(GranulomaBurden = "LowBurden")
lowburden.groupsummarize$ReceiverGroup <- paste0("R.", lowburden.groupsummarize$ReceiverGroup)
lowburden.groupsummarize$SenderGroup <- paste0("S.", lowburden.groupsummarize$SenderGroup)

dataframe.choice <- rbind(highburden.groupsummarize, lowburden.groupsummarize)

order <- c("Group1", "Group2", "Group3", "Group4", "Group5")

mapping.df <- data.frame(sender.group = paste0("S.", order),
                         receiver.group = paste0("R.", order),
                         colors = c(group1.color, group2.color, group3.color, group4.color, group5.color))

sender.order <- paste0("S.", order)
sender.order <- sender.order[sender.order %in% dataframe.choice$SenderGroup]
receiver.order <- paste0("R.", rev(order))
receiver.order <- receiver.order[receiver.order %in% dataframe.choice$ReceiverGroup]
order.vec <- c(rev(sender.order), rev(receiver.order))

##### Do work to get colors into the right format - vector where names are cell types, and entries are colors #####
colors.groups <- mapping.df$colors
names(colors.groups) <- order
sender.colors.groups <- colors.groups
names(sender.colors.groups) <- paste0("S.", order)
sender.colors.groups <- sender.colors.groups[names(sender.colors.groups) %in% sender.order]

receiver.colors.groups <- rev(colors.groups)
names(receiver.colors.groups) <- paste0("R.", rev(order))
receiver.colors.groups <- receiver.colors.groups[names(receiver.colors.groups) %in% receiver.order]

combined.colors <- c(sender.colors.groups, receiver.colors.groups)

pdf(paste0(figdir, "/", "CircosPlot_AllInteracts_SplitbyBurden_HighBurden.pdf"), width = 10, height = 10, useDingbats = FALSE)
circos.clear()
circos.par(start.degree = 270)
chordDiagram(dataframe.choice, directional = -1,
             order = order.vec, big.gap = 15,
             link.visible = ((dataframe.choice$GranulomaBurden == "HighBurden")),
             grid.col = combined.colors, link.zindex = rank(dataframe.choice[[3]]),
             direction.type = c("arrows", "diffHeight"), diffHeight  = 0, link.arr.type = "big.arrow",
             annotationTrack = c("grid"), annotationTrackHeight = c(0.04, 0.01),
             preAllocateTracks = list(track.height = mm_h(6), track.margin = c(mm_h(4), 0)))
for(si in get.all.sector.index()) {
  label.ii <- substr(si, 3, nchar(si))
  xlim = get.cell.meta.data("xlim", sector.index = si, track.index = 2)
  ylim = get.cell.meta.data("ylim", sector.index = si, track.index = 2)
  circos.text(mean(xlim), mean(ylim), label.ii, sector.index = si, track.index = 2, 
              facing = "bending.inside", niceFacing = TRUE, col = "white")
}
highlight.sector(dataframe.choice$SenderGroup, track.index = 1, col = "#f46d5a", 
    text = "Sender", cex = 1.2, text.col = "white", facing = "bending.inside", niceFacing = FALSE)
highlight.sector(dataframe.choice$ReceiverGroup, track.index = 1, col = "#ffb1a6", 
    text = "Receiver", cex = 1.2, facing = "bending.inside", text.col = "white", niceFacing = FALSE)
title("Network of Interactions\nStrengthened in High Burden Granulomas")
dev.off()

pdf(paste0(figdir, "/", "CircosPlot_AllInteracts_SplitbyBurden_LowBurden.pdf"), width = 10, height = 10, useDingbats = FALSE)
circos.clear()
circos.par(start.degree = 270)
chordDiagram(dataframe.choice, directional = -1,
             order = order.vec, big.gap = 15,
             link.visible = ((dataframe.choice$GranulomaBurden == "LowBurden")),
             grid.col = combined.colors, link.zindex = rank(dataframe.choice[[3]]),
             direction.type = c("arrows", "diffHeight"), diffHeight  = 0, link.arr.type = "big.arrow",
             annotationTrack = c("grid"), annotationTrackHeight = c(0.04, 0.01),
             preAllocateTracks = list(track.height = mm_h(6), track.margin = c(mm_h(4), 0)))
for(si in get.all.sector.index()) {
  label.ii <- substr(si, 3, nchar(si))
  xlim = get.cell.meta.data("xlim", sector.index = si, track.index = 2)
  ylim = get.cell.meta.data("ylim", sector.index = si, track.index = 2)
  circos.text(mean(xlim), mean(ylim), label.ii, sector.index = si, track.index = 2, 
              facing = "bending.inside", niceFacing = TRUE, col = "white")
}
highlight.sector(dataframe.choice$SenderGroup, track.index = 1, col = "#f46d5a", 
    text = "Sender", cex = 1.2, text.col = "white", facing = "bending.inside", niceFacing = FALSE)
highlight.sector(dataframe.choice$ReceiverGroup, track.index = 1, col = "#ffb1a6", 
    text = "Receiver", cex = 1.2, facing = "bending.inside", text.col = "white", niceFacing = FALSE)
title("Network of Interactions\nStrengthened in Low Burden Granulomas")
dev.off()

circos.low.fcn <- function(){
  circos.clear()
  circos.par(start.degree = 270)
  chordDiagram(dataframe.choice, directional = -1,
               order = order.vec, big.gap = 15,
               link.visible = ((dataframe.choice$GranulomaBurden == "LowBurden")),
               grid.col = combined.colors, link.zindex = rank(dataframe.choice[[3]]),
               direction.type = c("arrows", "diffHeight"), diffHeight  = 0, link.arr.type = "big.arrow",
               annotationTrack = c("grid"), annotationTrackHeight = mm_h(4),
               preAllocateTracks = list(track.height = mm_h(6), track.margin = c(mm_h(2), 0)))
  for(si in get.all.sector.index()) {
    label.ii <- substr(si, 3, nchar(si))
    xlim = get.cell.meta.data("xlim", sector.index = si, track.index = 2)
    ylim = get.cell.meta.data("ylim", sector.index = si, track.index = 2)
    circos.text(mean(xlim), mean(ylim), label.ii, sector.index = si, track.index = 2, 
                facing = "bending.inside", niceFacing = TRUE, col = "white", cex = 0.8)
  }
  highlight.sector(dataframe.choice$SenderGroup, track.index = 1, col = "#f46d5a", 
                   text = "Sender", cex = 1.2, text.col = "white", facing = "bending.inside", niceFacing = FALSE)
  highlight.sector(dataframe.choice$ReceiverGroup, track.index = 1, col = "#ffb1a6", 
                   text = "Receiver", cex = 1.2, facing = "bending.inside", text.col = "white", niceFacing = FALSE)
  title("Network of Interactions\nStrengthened in Low Burden Granulomas")
}

circos.high.fcn <- function(){
  circos.clear()
  circos.par(start.degree = 270)
  chordDiagram(dataframe.choice, directional = -1,
               order = order.vec, big.gap = 15,
               link.visible = ((dataframe.choice$GranulomaBurden == "HighBurden")),
               grid.col = combined.colors, link.zindex = rank(dataframe.choice[[3]]),
               direction.type = c("arrows", "diffHeight"), diffHeight  = 0, link.arr.type = "big.arrow",
               annotationTrack = c("grid"), annotationTrackHeight = mm_h(4),
               preAllocateTracks = list(track.height = mm_h(6), track.margin = c(mm_h(2), 0)))
  for(si in get.all.sector.index()) {
    label.ii <- substr(si, 3, nchar(si))
    xlim = get.cell.meta.data("xlim", sector.index = si, track.index = 2)
    ylim = get.cell.meta.data("ylim", sector.index = si, track.index = 2)
    circos.text(mean(xlim), mean(ylim), label.ii, sector.index = si, track.index = 2, 
                facing = "bending.inside", niceFacing = TRUE, col = "white", cex = 0.8)
  }
  highlight.sector(dataframe.choice$SenderGroup, track.index = 1, col = "#f46d5a", 
                   text = "Sender", cex = 1.2, text.col = "white", facing = "bending.inside", niceFacing = FALSE)
  highlight.sector(dataframe.choice$ReceiverGroup, track.index = 1, col = "#ffb1a6", 
                   text = "Receiver", cex = 1.2, facing = "bending.inside", text.col = "white", niceFacing = FALSE)
  title("Network of Interactions\nStrengthened in High Burden Granulomas")
}

circos.agrob.low <- cowplot::as_grob(circos.low.fcn)
circos.agrob.high <- cowplot::as_grob(circos.high.fcn)

# Code for figuring out angle needed for circos plot legend
number.interacts.received <- dataframe.choice[,"group.count.interacts"] %>% sum()
number.interacts.low <- dataframe.choice[dataframe.choice$GranulomaBurden == "LowBurden","group.count.interacts"] %>% sum()
number.interacts.high <- dataframe.choice[dataframe.choice$GranulomaBurden == "HighBurden","group.count.interacts"] %>% sum()
group2.interacts <- sum(as.matrix(dataframe.choice[(dataframe.choice$ReceiverGroup == "R.Group2"),"group.count.interacts"]))
group2.angle <- 330*group2.interacts / number.interacts.received
interacts.for.legend <- 200
angle.for.legend <- interacts.for.legend / group2.interacts * group2.angle

```

```{r Interaction count bar plots for ecosystem analyses}
travis.rdata.path <- "/Users/ctzouanas/Dropbox (MIT)/ConstantineData/Interactions/InteractionsMergeTotal.Rdata"
if(!exists("InteractionsMergeTotal")){
  load(travis.rdata.path)
  InteractionsMergeTotal <- InteractionsMergeTotal %>% 
    dplyr::select(c(RLTypes, HighLowPVal, LowPermPvals, ComboReceptorMeans_High, ComboReceptorMeans_Low)) %>%
    dplyr::rename(ComboReceptorHighLowPVal = HighLowPVal) %>%
    dplyr::mutate(ComboReceptorHighLowFC = ComboReceptorMeans_High / ComboReceptorMeans_Low) %>% 
    tidyr::replace_na(list(ComboReceptorHighLowPVal = 1, ComboReceptorHighLowFC = 0)) %>%
    dplyr::mutate(ComboReceptorHighLowFC = replace(ComboReceptorHighLowFC, is.infinite(ComboReceptorHighLowFC), 1.5))
  interact.df <- InteractionsMergeTotal %>% dplyr::filter(ComboReceptorHighLowPVal < 0.05 & LowPermPvals < 0.05) %>%
    tidyr::separate(col = RLTypes, into = c("Sender", "Receiver", "Ligand", "Receptor"), sep = "_") %>%
    dplyr::mutate(Sender_Ligand = paste0(Sender, "_", Ligand), 
                  Receiver_Receptor = paste0(Receiver, "_", Receptor))
}

interact.sender.low.df <- interact.df %>% dplyr::filter(ComboReceptorHighLowFC < 1) %>% group_by(Sender) %>% summarise(neg_count = n())
interact.sender.high.df <- interact.df %>% dplyr::filter(ComboReceptorHighLowFC > 1) %>% group_by(Sender) %>% summarise(pos_count = n())
interact.sender.count.df <- full_join(x = interact.sender.low.df, y = interact.sender.high.df, by = "Sender", fill = 0)
interact.sender.count.df[is.na(interact.sender.count.df)] <- 0

interact.sender.count.df <- interact.sender.count.df %>% mutate(net_interact = pos_count - neg_count)
interact.sender.count.df$SenderGroup <- ""
interact.sender.count.df$SenderGroup[interact.sender.count.df$Sender %in% group1] <- "Group1"
interact.sender.count.df$SenderGroup[interact.sender.count.df$Sender %in% group2] <- "Group2"
interact.sender.count.df$SenderGroup[interact.sender.count.df$Sender %in% group3] <- "Group3"
interact.sender.count.df$SenderGroup[interact.sender.count.df$Sender %in% group4] <- "Group4"
interact.sender.count.df$SenderGroup[interact.sender.count.df$Sender %in% group5] <- "Group5"
interact.sender.count.df$SenderGroup <- factor(interact.sender.count.df$SenderGroup, levels = c("Group1", "Group2", "Group3", "Group4", "Group5"))

plot.df <- interact.sender.count.df
plot.df$Sender <- fct_reorder(plot.df$Sender, plot.df$pos_count, .desc = TRUE)
pdf(paste0(figdir, "/", "HighCFUInteractCounts_BarChart.pdf"), width = 6, height = 8, useDingbats = FALSE)
p.highcfu.interact.counts <- ggplot(plot.df, aes(x = Sender, y = pos_count, fill = SenderGroup)) + geom_bar(stat = "identity") + 
  theme_classic() + scale_fill_manual(values = c(group1.color, "#ca5e4a", "#4e88b9", group4.color, group5.color)) + 
  theme(text = element_text(size = font.size), legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1, size = 0.5*font.size)) + 
  scale_y_continuous(expand = c(0,0)) + ylab("Number of\nHigh CFU-Linked\nInteractions")
plot(p.highcfu.interact.counts)
dev.off()

plot.df$Sender <- fct_reorder(plot.df$Sender, plot.df$neg_count, .desc = TRUE)
pdf(paste0(figdir, "/", "LowCFUInteractCounts_BarChart.pdf"), width = 6, height = 8, useDingbats = FALSE)
p.lowcfu.interact.counts <- ggplot(plot.df, aes(x = Sender, y = abs(neg_count), fill = SenderGroup)) + geom_bar(stat = "identity") + 
  theme_classic() + scale_fill_manual(values = c(group1.color, "#ca5e4a", "#4e88b9", group4.color, group5.color)) + 
  theme(text = element_text(size = font.size), legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1, size = 0.5*font.size)) + 
  ylab("Number of\nLow CFU-Linked\nInteractions") + scale_y_continuous(expand = c(0,0))
plot(p.lowcfu.interact.counts)
dev.off()

group.interact.sender.count.df <- interact.sender.count.df %>% group_by(SenderGroup) %>% summarise(neg_count = sum(neg_count), pos_count = sum(pos_count))

```

```{r Analyze what cell type + ligand combos are most associated with high or low granuloma burden}
pospleioneg.colorvec <- c("#ad2524", "#572b88", "#046c9a")
font.size = 16

##### Find where top ligands come from, in each burden condition #####
average.sender.ligand.strength <- interact.df %>% dplyr::group_by(Sender_Ligand) %>% 
  summarise(mean_FC = mean(ComboReceptorHighLowFC), sd_FC = sd(ComboReceptorHighLowFC), count_FC = n())
split.strings <- str_split(string = average.sender.ligand.strength$Sender_Ligand, pattern = "_") %>% data.frame() %>% t()
colnames(split.strings) <- c("Sender", "Ligand")
average.sender.ligand.strength <- cbind(average.sender.ligand.strength, split.strings)
average.sender.ligand.strength$SenderGroup <- ""
average.sender.ligand.strength$SenderGroup[average.sender.ligand.strength$Sender %in% group1] <- "Group1"
average.sender.ligand.strength$SenderGroup[average.sender.ligand.strength$Sender %in% group2] <- "Group2"
average.sender.ligand.strength$SenderGroup[average.sender.ligand.strength$Sender %in% group3] <- "Group3"
average.sender.ligand.strength$SenderGroup[average.sender.ligand.strength$Sender %in% group4] <- "Group4"
average.sender.ligand.strength$SenderGroup[average.sender.ligand.strength$Sender %in% group5] <- "Group5"
average.sender.ligand.strength <- average.sender.ligand.strength %>% arrange(-mean_FC)
average.sender.ligand.strength$rank <- 1:nrow(average.sender.ligand.strength)
zero.truth.vec <- average.sender.ligand.strength$mean_FC == 0
average.sender.ligand.strength$mean_FC[zero.truth.vec] = min(average.sender.ligand.strength$mean_FC[!zero.truth.vec])

##### Set up bar plots of where the top ligands (whether most enriched in high or low burden) come from #####
n.ligand <- round(0.1*nrow(average.sender.ligand.strength))
high.ligands.group.count <- average.sender.ligand.strength %>% dplyr:: filter(rank <= n.ligand) %>% group_by(SenderGroup) %>% summarise(ligand.count = n()) %>% dplyr::mutate(ligand.prop = 100*ligand.count / sum(ligand.count), GranulomaBurden = "HighBurden")
low.ligands.group.count <- average.sender.ligand.strength %>% dplyr:: filter(rank > (nrow(average.sender.ligand.strength) - n.ligand)) %>% group_by(SenderGroup) %>% summarise(ligand.count = n()) %>% dplyr::mutate(ligand.prop = 100*ligand.count / sum(ligand.count), GranulomaBurden = "LowBurden")
all.ligands.group.count <- rbind(high.ligands.group.count, low.ligands.group.count)
pdf(paste0(figdir, "/", "Top10PercentLigands_HighBurdenGranulomas.pdf"), useDingbats = FALSE, width = 2, height = 10)
p.topligands.highburden <- ggplot(all.ligands.group.count %>% dplyr::filter(GranulomaBurden == "HighBurden"), aes(fill = SenderGroup, x = GranulomaBurden, y = ligand.prop)) +
  geom_bar(stat = "identity", color = "black") + theme_classic() + theme(text = element_text(size = font.size), axis.text.x = element_text(angle = 45, hjust = 1, size = 0.5*font.size), legend.position = "none") +
  scale_fill_manual(values = c(group1.color, "#ca5e4a", "#4e88b9", group4.color, group5.color)) +
  scale_y_continuous(expand = c(0,0)) + ylab('Percent of Top Decile Ligands\nSecreted from Group') + xlab('') + 
  scale_x_discrete(labels = c("High Burden\nGranuloma"))
plot(p.topligands.highburden)
dev.off()

highburden.ligands.group.count <- all.ligands.group.count %>% dplyr::filter(GranulomaBurden == "HighBurden")
binom.test(x = highburden.ligands.group.count$ligand.count[highburden.ligands.group.count$SenderGroup == "Group2"],
           n = sum(highburden.ligands.group.count$ligand.count),
           p = 1/nrow(highburden.ligands.group.count))

pdf(paste0(figdir, "/", "Top10PercentLigands_LowBurdenGranulomas.pdf"), useDingbats = FALSE, width = 2, height = 10)
p.topligands.lowburden <- ggplot(all.ligands.group.count %>% dplyr::filter(GranulomaBurden == "LowBurden"), aes(fill = SenderGroup, x = GranulomaBurden, y = ligand.prop)) +
  geom_bar(stat = "identity", color = "black") + theme_classic() + theme(text = element_text(size = font.size), axis.text.x = element_text(angle = 45, hjust = 1, size = 0.5*font.size), legend.position = "none") +
  scale_fill_manual(values = c(group1.color, "#ca5e4a", "#4e88b9", group4.color, group5.color)) +
  scale_y_continuous(expand = c(0,0)) + ylab('Percent of Top Decile Ligands\nSecreted from Group') + xlab('') + 
  scale_x_discrete(labels = c("Low Burden\nGranuloma"))
plot(p.topligands.lowburden)
dev.off()

lowburden.ligands.group.count <- all.ligands.group.count %>% dplyr::filter(GranulomaBurden == "LowBurden")
binom.test(x = lowburden.ligands.group.count$ligand.count[lowburden.ligands.group.count$SenderGroup == "Group3"],
           n = sum(lowburden.ligands.group.count$ligand.count),
           p = 1/nrow(lowburden.ligands.group.count))

##### Set up bar plot of top ligands, averaged across all senders #####
average.ligand.strength <- interact.df %>% dplyr::group_by(Ligand) %>%
  dplyr::mutate(HighBurdenvsLowBurden_LigandStrength_FC = mean(ComboReceptorHighLowFC), 
                sd_FC = sd(ComboReceptorHighLowFC), 
                total_interactions = n(),
                pos_interactions = sum(ComboReceptorHighLowFC > 1),
                neg_interactions = sum(ComboReceptorHighLowFC < 1)) %>% 
  dplyr::select(Ligand, HighBurdenvsLowBurden_LigandStrength_FC, sd_FC, total_interactions, pos_interactions, neg_interactions) %>% 
    unique() %>%
  dplyr::mutate(log10_ligand_mean_FC = log10(HighBurdenvsLowBurden_LigandStrength_FC)) %>% arrange(-HighBurdenvsLowBurden_LigandStrength_FC)
write.table(average.ligand.strength,
            file = paste0(figdir, "/", "LigandStrength_HighBurdenvsLowBurden_SummarizeAcrossAllSenders.csv"),
            sep = ",", row.names = FALSE)

chosen.ligands <- c("TNFSF11", "TNFSF15", "FGF7", "IL33", "TGFB2", "TGFB3", "WNT4", "RSPO3", "JAG2", "CX3CL1", "IL4", "IL5", "IL13", "IL24", "PDGFB", "CD5L", "IL12A", "VEGFC", "VEGFB", "CCL11", "ANGPTL4", "DCN", "FN1", "VCAM1", "NRG1", "RLN3", "NTS", "IL12B",  "IL15", "B2M", "TNFSF13", "CCL19", "IL1B", "IL1A", "IFNG", "CXCL10", "LYZ", "TGFB1", "GZMB", "SELPLG", "ICAM2", "ICAM4", "FASLG", "CD40LG", "CD58", "CXCL16", "CD80", "CD48", "IL18", "CXCL9", "CXCL11", "CCL20", "CCL3", "CD86", "LTA", "CCL5", "LTB", "TNF", "CXCL10", "CCL17", "CCL22", "IL2", "IL27", "UCN3", "RTN4", "NRG2")

average.ligand.strength <- average.ligand.strength %>% dplyr::filter(Ligand %in% chosen.ligands)

plot.heatmap.df <- interact.df %>% group_by(Sender, Ligand) %>% summarise(HighvsLowFC_AverageAcrossReceivers = mean(ComboReceptorHighLowFC)) %>%
  dplyr::mutate(log10_ligand_mean_FC = log10(HighvsLowFC_AverageAcrossReceivers))
write.table(plot.heatmap.df, 
            file = paste0(figdir, "/", "LigandStrength_HighBurdenvsLowBurden_SummarizeByEachSender.csv"), 
            sep = ",", row.names = FALSE)
plot.heatmap.df <- plot.heatmap.df %>% dplyr::filter(Ligand %in% chosen.ligands) %>% 
  group_by(Ligand) %>%
  dplyr::mutate(total_count = n(), 
                pos_count = sum(log10_ligand_mean_FC > 0),
                neg_count = sum(log10_ligand_mean_FC < 0))
ligand.labels <- plot.heatmap.df %>% 
  dplyr::select(Ligand, total_count, pos_count, neg_count) %>%
  unique() %>%
  dplyr::mutate(Sign = case_when(total_count == pos_count ~ "Positive",
                                 total_count == neg_count ~ "Negative",
                                 TRUE ~ "Pleiotropic"))

min.finite <- min(plot.heatmap.df$log10_ligand_mean_FC[is.finite(plot.heatmap.df$log10_ligand_mean_FC)])
max.finite <- max(plot.heatmap.df$log10_ligand_mean_FC[is.finite(plot.heatmap.df$log10_ligand_mean_FC)])
plot.heatmap.df$log10_ligand_mean_FC[plot.heatmap.df$log10_ligand_mean_FC == -Inf] = 1.05*min.finite
plot.heatmap.df$log10_ligand_mean_FC[plot.heatmap.df$log10_ligand_mean_FC == Inf] = 1.05*max.finite

plot.df <- average.ligand.strength %>% arrange(-log10_ligand_mean_FC)
plot.df <- inner_join(x = plot.df, y = ligand.labels, by = "Ligand") %>% 
  arrange(-log10_ligand_mean_FC)
plot.df$Ligand <- factor(plot.df$Ligand, levels = plot.df$Ligand)
plot.df$Sign = factor(plot.df$Sign, levels = c("Negative", "Pleiotropic", "Positive"))
p.bar <- ggplot(plot.df, aes(x = Ligand, y = log10_ligand_mean_FC, fill = Sign)) + geom_bar(stat = "identity") + 
  theme_classic() + 
  theme(text = element_text(size = font.size), legend.position = "none", axis.title.x = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) + scale_fill_manual(values = rev(pospleioneg.colorvec)) +
  ylab("logFC(Interaction Strength)\nAcross High vs. Low\nBurden Granulomas") + scale_y_continuous(expand = c(0,0))
plot.heatmap.df$Ligand <- factor(plot.heatmap.df$Ligand, levels = plot.df$Ligand)
order <- c(group1, group2, group3, group4, group5)
plot.heatmap.df$Sender <- factor(plot.heatmap.df$Sender, levels = rev(order))
plot.heatmap.df$SenderGroup <- ""
plot.heatmap.df$SenderGroup[plot.heatmap.df$Sender %in% group1] <- "Group1"
plot.heatmap.df$SenderGroup[plot.heatmap.df$Sender %in% group2] <- "Group2"
plot.heatmap.df$SenderGroup[plot.heatmap.df$Sender %in% group3] <- "Group3"
plot.heatmap.df$SenderGroup[plot.heatmap.df$Sender %in% group4] <- "Group4"
plot.heatmap.df$SenderGroup[plot.heatmap.df$Sender %in% group5] <- "Group5"
p.heatmap <- ggplot(plot.heatmap.df, aes(x = Ligand, y = Sender, fill = log10_ligand_mean_FC)) + geom_tile() + theme_classic() +
  scale_x_discrete(position = "top", expand = c(0,0)) + scale_y_discrete(expand = c(0,0)) + 
  scale_fill_gradient2(low = "#016099", mid = "#ffffff", high = "#ca5e4a") + facet_grid(SenderGroup~., scales = "free", space = "free") + #facet by group
  theme(text = element_text(size = font.size), legend.position = "bottom", axis.title.x = element_blank(), 
        axis.text.x = element_blank(), axis.title.y = element_blank(),
        strip.background = element_blank(), #remove background for facet labels
        panel.border = element_rect(colour = "black", fill = NA), #add black border
        panel.spacing = unit(0, "lines"))

p.ligands.combined <- cowplot::plot_grid(plotlist = list(p.bar, p.heatmap), align = "v", 
                                         nrow = 2, ncol = 1, axis = "lr", rel_heights = c(1, 1.25))
plot(p.ligands.combined)
pdf(paste0(figdir, "/", "StrongestLigandsbyCellType_BarChartHeatmap.pdf"), useDingbats = FALSE, width = 14, height = 10)
plot(p.ligands.combined)
dev.off()

```

```{r Pull ecosystem plots together}
height.in = 30
width.in = 15
heights.in = c(1.1, 0.75, 1.25, 0.75, 1)
widths.interact.in = c(1, 0.25, 1, 0.25)
widths.gsea.in = c(1,1)
nrow.in = 5

plotlist <- list(correlation.heatmap[[4]], granuloma.summary.plot + theme(axis.text.x = element_blank()), allgranuloma.barplot + theme(axis.title.x = element_blank()),
                 NULL, NULL, granuloma.metadata.plot + theme(axis.text.y = element_blank()))
compiled.plot <- plot_grid(plotlist = plotlist, align = "hv", axis = "tblr", ncol = 3, 
                           rel_widths = c(3, 1, 3), rel_heights = c(4,1))
pdf(paste0(figdir, "/", "CompositionBarplots_and_CorrHeatmap_NoLegend.pdf"), width = 12, height = 6, useDingbats = FALSE)
plot(compiled.plot)
dev.off()

pdf(paste0(figdir, "/", "EcosystemCompiled_NoLabels.pdf"), width = width.in, height = height.in, useDingbats = FALSE)
ggarrange(
  compiled.plot,
  ggarrange(p.highcfu.interact.counts, p.topligands.highburden, p.lowcfu.interact.counts, p.topligands.lowburden, ncol = 4, widths = widths.interact.in, align = "hv"),
  ggarrange(circos.agrob.high, circos.agrob.low, ncol = 2, align = "hv"),
  p.bar + theme(axis.text.y = element_blank(), axis.title.y = element_blank()), 
  p.heatmap + theme(legend.position = "none", axis.text.y = element_blank(), axis.title.y = element_blank(), strip.text.y = element_blank()),
  nrow = nrow.in, heights = heights.in
)
dev.off()

plotlist.labeled <- list(correlation.heatmap[[4]], granuloma.summary.plot, allgranuloma.barplot,
                         NULL, NULL, granuloma.metadata.plot)
compiled.plot.labeled <- plot_grid(plotlist = plotlist.labeled, align = "hv", axis = "tblr", ncol = 3, 
                                   rel_widths = c(3, 1, 3), rel_heights = c(4,1))
pdf(paste0(figdir, "/", "CompositionBarplots_and_CorrHeatmap_WithLegend.pdf"), width = 12, height = 6, useDingbats = FALSE)
plot(compiled.plot.labeled)
dev.off()

pdf(paste0(figdir, "/", "EcosystemCompiled_WithLabels.pdf"), width = width.in, height = height.in, useDingbats = FALSE)
ggarrange(
  compiled.plot.labeled,
  ggarrange(p.highcfu.interact.counts, p.topligands.highburden, p.lowcfu.interact.counts, p.topligands.lowburden, ncol = 4, widths = widths.interact.in, align = "hv"),
  ggarrange(circos.agrob.high, circos.agrob.low, ncol = 2, align = "hv"),
  p.bar, 
  p.heatmap,
  nrow = nrow.in, heights = heights.in
)
dev.off()
```
