---
title: "CodeForImmunity"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

### Figure 1B - Boxplot Log10 CFU/Granuloma by Animal
```{r,}
CFUAnimalBoxplot <- ggplot(GranulomaMetadata,aes(x=as.factor(Animal),y=logCFU,fill=as.factor(Animal))) + geom_boxplot(outlier.shape=NA) +  geom_jitter(shape=16, position=position_jitter(0.2)) +theme_classic() +theme_classic() + theme(legend.position = "none")

Plot <- GranulomaMetadata %>% ggplot(aes(x=as.factor(Animal), y=logCFU, fill=as.factor(Animal))) + stat_boxplot(geom ='errorbar',width=0.3,coef=10) + geom_boxplot(outlier.shape = NA) + geom_point(position=position_jitterdodge(jitter.width = 0.75),alpha=0.8,shape=21,size=5) + scale_fill_manual(values=AnimalColors)  + theme_classic() + theme(
axis.line.x = element_line(colour = 'black', size=0.0, linetype='solid'),
axis.line.y = element_line(colour = 'black', size=0.0, linetype='solid')) +  theme(legend.position="none",axis.text.x=element_text(angle=45, hjust=1)) + ylab("log10(CFU)") + xlab("Animal") + scale_x_discrete(labels = c("1","2","3","4")) +theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"), panel.border = element_rect(linetype = "solid", colour = "black", fill = "NA", size = 0.5))

```

### Figure 1C - Boxplot Log10 CFU/Granuloma by Burden (High vs. Low)
```{r,}
CEQPlot <- ggplot(GranulomaMetadata,aes(x=HighLow,y=Log10.CFU,group=as.factor(HighLow))) + geom_boxplot(outlier.shape=NA) +  geom_jitter(shape=16, position=position_jitter(0.2)) 
```


### Figure 1D - Boxplot Log10 CEQ/Granuloma by Burden (High vs. Low)
```{r,}
CEQPlot <- ggplot(GranulomaMetadata,aes(x=HighLow,y=Log10.CEQ,group=as.factor(HighLow))) + geom_boxplot(outlier.shape=NA) +  geom_jitter(shape=16, position=position_jitter(0.2)) 
```


### Figure 1E - Boxplot Log10 CFU/CEQ per Granuloma by Burden (High vs. Low)
```{r,}

CFUCEQBoxplot <- ggplot(GranulomaMetadata,aes(x=as.factor(HighLow),y=Log10.CFU.CEQ,group=as.factor(HighLow))) + geom_boxplot(outlier.shape=NA) +  geom_jitter(shape=16, position=position_jitter(0.2)) 
```

### Figure 1G - Boxplot Log10 CFU/Granuloma by Time of Detection (Early vs. Late)
```{r,}
AgeCFUBoxplot <- ggplot(GranulomaMetadata,aes(x=as.factor(Age),y=logCFU,group=as.factor(Age))) + geom_boxplot(outlier.shape=NA) +  geom_jitter(shape=16, position=position_jitter(0.2)) 
```

### Figure 1H - Boxplot Log10 CEQ/Granuloma by Time of Detection (Early vs. Late)
```{r,}
load("Metadata/CFUDataFinal_Update.Rdata")
AgeCEQBoxplot <- ggplot(GranulomaMetadata,aes(x=Age,y=Log10.CEQ,group=as.factor(Age))) + geom_boxplot(outlier.shape=NA) +  geom_jitter(shape=16, position=position_jitter(0.2)) +theme_classic()
```

### Figure 1I - Boxplot Log10 CFU/CEQ per Granuloma by Time of Detection (Early vs. Late)
```{r,}
AgeCFUCEQBoxplot <- ggplot(GranulomaMetadata,aes(x=as.factor(Age),y=Log10.CFU.CEQ,group=as.factor(Age))) + geom_boxplot(outlier.shape=NA) +  geom_jitter(shape=16, position=position_jitter(0.2)) +theme_classic()
```
 
### Figure 2A - Global UMAP
```{r,}
GlobalUMAPCellType_labeled <- ggplot(GranulomaCellLevelMetadata,aes(x=X_umap1,y=X_umap2,color=as.factor(GenericFinal))) + theme(axis.line = element_blank())+ xlab("UMAP 1") + ylab("UMAP 2")+ geom_point_rast(size=0.25)  + theme(panel.grid = element_blank(),axis.text = element_blank(),axis.ticks = element_blank(),panel.background = element_blank()) +  scale_color_manual(values=  GenericColors3)
GlobalUMAPCellType_labeled
 
```

### Figure 2B - Global Heatmap
```{r,}


### Gene Selection
GenericMarkersTop4 <- GenericMarkers %>% group_by(cluster) %>% top_n(n =4, wt = avg_logFC)
GenericMarkersTop4 <- GenericMarkersTop4[order(as.character(GenericMarkersTop4$cluster)),]
GeneAnnotation <- as.character(GenericMarkersTop4$gene)x

### Code to Downsample Cells 
CellTypeList <- as.character(unique(GranulomaMinus23Metadata$GenericFinal))
CellTypeList <- CellTypeList[order(CellTypeList)]

DownSampleMatrix <- function(Metadata){
nCells <- min(table(Metadata$GenericFinal))
IDsSubset <-c()
for (i in 1:length(unique(Metadata$GenericFinal))){
  CellTypeUse <- CellTypeList[i]
  CellsUse <- as.character(Metadata$CellID[Metadata$GenericFinal == CellTypeUse])
  IDsUse <- sample(CellsUse,nCells)
  IDsSubset <- c(IDsSubset,IDsUse)
}
assign("IDsSubset",IDsSubset,envir=.GlobalEnv)
}

DownSampleMatrix(GranulomaCellLevelMetadata)

### Format Column Data - Pheatmap
ColumnData <- as.data.frame(MetadataSubset$GenericFinal)
colnames(ColumnData) <- "CellType"
row.names(ColumnData) <- MetadataSubset$CellID

Breaks <- c()
for(i in 1:13){
Break <- sum(summary(as.factor(ColumnData$CellType))[1:i])
 Breaks <- c(Breaks,Break)
}

PrettyHeatmap <- function(NormData,GeneList,Clusters,Name){ 
  #NormData <- log10(NormData + 1)
  NormDataUse <- NormData[GeneList,]
 Clusters$CellID <- row.names(Clusters)
  Clusters <- as.data.frame(Clusters[order(Clusters$CellType),])
  row.names(Clusters) <- Clusters$CellID
  Clusters$CellID <- NULL
  NormDataUse <- NormDataUse[,row.names(Clusters)]
  Clusters$CellType <- as.factor(Clusters$CellType)
  #NormDataUse <- as.data.frame(t(NormDataUse))
  #NormDataUse <- as.data.frame(t(scale(NormDataUse)))
  PlotName <- paste(Name,"Heatmap.pdf",sep="_")
   pdf(file=PlotName,height=8,width=12)
  my_heatmap <- pheatmap(NormDataUse,annotation_col=Clusters,border_col=NA,show_colnames = F, show_rownames = T,legend = TRUE,
             cluster_rows=F, cluster_cols=F, gaps_col=Breaks
)
  dev.off()
}

PrettyHeatmap(GranulomaNormal,GeneAnnotation,ColumnData,"OutputDirectory")



```


### Figure 2C - Individual Spearman's Rho Scatterplots
```{r,}
load("GenericProportions.Rdata")
load("GranulomaMetadata.Rdata")

GenericProportionsMerge <- merge(GenericProportions,GranulomaMetadata,by="Granuloma")

BcellPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=B)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("B cells")


cDCPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=cDC)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("cDCs")

EndoPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=Endo)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Endothelial cells")

FibroPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=Fibro)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Fibroblasts")

MastPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=Mast)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Mast cells")

MphagePlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=Mphage)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Macrophages")

NeutrophilPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=Neutrophil)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Neutrophils")


pDCPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=pDC)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("pDCs")

PlasmaPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=Plasma)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Plasma cells")

RBCPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=RBC)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Erythrocytes")

TcellPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=T)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("T cells")

T1PPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=T1P)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Type 1 Pneumocytes")

T2PPlot <- ggplot(GenericProportionsMerge,aes(x=logCFU,y=T2P)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Type 2 Pneumocytes") + theme(plot.title = element_text(hjust = 0.5))

pdf(file="FigurePanels/GenericSupplement/GenericProportionsScatter.pdf",height=11,width=38)
print(ggarrange(BcellPlot,cDCPlot,EndoPlot,FibroPlot,MastPlot,MphagePlot,NeutrophilPlot,pDCPlot,PlasmaPlot,RBCPlot,TcellPlot,T1PPlot,T2PPlot,ncol=7,nrow=2))
dev.off()

```

### Figure 3A - T cell UMAP
```{r,}
TcellTotalUMAP_Labeled <- ggplot(TcellMetaFinal,aes(x=UMAP1,y=UMAP2,color=as.factor(TcellPhenoFinal))) + theme(axis.line = element_blank()) + theme(axis.line = element_blank())+ xlab("UMAP 1") + ylab("UMAP 2")+ geom_point_rast(size=0.50)  + theme(panel.grid = element_blank(),axis.text = element_blank(),axis.ticks = element_blank(),panel.background = element_blank()) + scale_color_manual(values=TcellColors)

```

### Figure 3B - Barplots by Cluster
```{r,}
load("CD4CD8Proportions.Rdata") 
barplot(t(CD4CD8Proportions))
```

### Figure 3C - T cell Dotplot
```{r,}

load("TcellSeurat.Rdata")
TcellMarkersuse <- c("CCR6","FURIN","COTL1","TOP2A","MKI67","PNISR","SRRM2","LOC102140106","LOC102137613","IFI6","ISG15","MX1","GATA3","FOXP3","TIGIT","JUNB","FOS","TCF7","GZMK","CX3CR1","TGFBR3","S100A10","XCL1","KLRC1","CD69","EGR1","KLRD1","GZMA","GZMH","KLRC2","TYROBP","GNLY")

pdf(file="TcellDotplot.pdf",height=11,width=8,useDingbats = FALSE)
DotPlot(TcellSeuratUpdate,features=rev(TcellMarkersuse),cols=c("lightgrey","darkred")) + RotatedAxis() + coord_flip()
dev.off()
```

### Figure 3D - T cell Scatterplots
```{r,}
load("Metadata/GranulomaMinus23Metadata.Rdata")

Counts <- table(GranulomaCellLevelMetadata$Specific,as.character(GranulomaCellLevelMetadata$Granuloma))
Proportions <- scale(Counts,scale=colSums(Counts),center=FALSE)*100 
Proportions <- as.data.frame(Proportions)


Proportions <- Proportions[Proportions$Var1 %in% TcellPopulations,]
colnames(Proportions) <- c("CellType","Granuloma","Frequency")
Proportions <- merge(Proportions,CFUDataSubset,by="Granuloma")
Proportions$Group <- paste(Proportions$CellType,Proportions$Tertile,sep="_")

# T cell Cluster 1
Proportions_1 <- Proportions[Proportions$CellType == "Cytotoxic1",]
T1_Plot <- ggplot(Proportions_1,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Cytotoxic1")

# T cell Cluster 2
Proportions_2 <- Proportions[Proportions$CellType == "Cytotoxic2",]
T2_Plot <- ggplot(Proportions_2,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Cytotoxic2")

# T cell Cluster 3
Proportions_3 <- Proportions[Proportions$CellType == "Cytotoxic3",]
T3_Plot <- ggplot(Proportions_3,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Cytotoxic3")


# T cell Cluster 4
Proportions_4 <- Proportions[Proportions$CellType == "Cytotoxic4",]
T4_Plot <- ggplot(Proportions_4,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Cytotoxic4")

Proportions_5 <- Proportions[Proportions$CellType == "Cytotoxic5",]
T5_Plot <- ggplot(Proportions_5,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Cytotoxic5")

Proportions_6 <- Proportions[Proportions$CellType == "Cytotoxic6",]
T6_Plot <- ggplot(Proportions_6,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Cytotoxic6")

Proportions_7 <- Proportions[Proportions$CellType == "Stem-like",]
T7_Plot <- ggplot(Proportions_7,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Stem-like")

Proportions_8 <- Proportions[Proportions$CellType == "T_Treg",]
T8_Plot <- ggplot(Proportions_8,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("T_Treg")

Proportions_9 <- Proportions[Proportions$CellType == "T_Interferon",]
T9_Plot <- ggplot(Proportions_9,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("T_Interferon")

Proportions_10 <- Proportions[Proportions$CellType == "T_Metallothionein",]
T10_Plot <- ggplot(Proportions_10,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("T_Metallothionein")

Proportions_11 <- Proportions[Proportions$CellType == "T_SRRM2",]
T11_Plot <- ggplot(Proportions_11,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("T_SRRM2")

Proportions_12 <- Proportions[Proportions$CellType == "T_Proliferating",]
T12_Plot <- ggplot(Proportions_12,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("T_Proliferating")

Proportions_13 <- Proportions[Proportions$CellType == "T_T1T17",]
T13_Plot <- ggplot(Proportions_13,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("T_T1T17")


```



### Figure 4 - T1T17 Sub-clustering
### Figure 4A - T cell Expression Overlays
```{r,}
load("Metadata/Tcell/TcellMetaFinal.Rdata")
T1T17HighlightUMAP <- ggplot(TcellMetaFinal,aes(x=UMAP1,y=UMAP2,color=as.character(T1T17))) + theme(axis.line = element_blank())+ xlab("UMAP 1") + ylab("UMAP 2")+ geom_point_rast(size=0.25) + scale_color_manual(values = c("lightgrey","darkred"))  + theme(legend.position = "none",panel.grid = element_blank(),axis.text = element_blank(),axis.ticks = element_blank(),panel.background = element_blank())

### Overlay of Individual Genes

GenePlot <- function(Metadata,NormData,Gene){
  GeneUse <- as.data.frame(t(NormData[Gene,]))
  colnames(GeneUse) <- Gene
  GeneUse$CellID <- row.names(GeneUse)
  Metadata <- merge(Metadata,GeneUse,by='CellID')
AssignName <- paste(Gene,"Plot",sep="_")
Metadata <- Metadata[order(Metadata[,Gene]),]
Plot <- ggplot(Metadata,aes(x=UMAP1,y=UMAP2,color=eval(parse(text=Gene)))) + theme(axis.line = element_blank())+ xlab("FA 1") + ylab("FA 2")+ geom_point_rast(size=0.5) + ggtitle(Gene)  + theme(legend.position = "none",panel.grid = element_blank(),axis.text = element_blank(),axis.ticks = element_blank(),panel.background = element_blank()) + scale_color_gradient(low='lightgrey',high='darkred',name=NULL)
assign(AssignName,Plot,envir=.GlobalEnv)

  PlotName <- paste("Directory/",Gene,"_Plot.pdf",sep="")

pdf(file=PlotName,height=5,width=5)
print(Plot)
dev.off()
}

  GenePlot(TcellMetaFinal,TcellNorm,"IL17A")
  GenePlot(TcellMetaFinal,TcellNorm,"IL17F")
  GenePlot(TcellMetaFinal,TcellNorm,"RORC")
  GenePlot(TcellMetaFinal,TcellNorm,"RORA")
  GenePlot(TcellMetaFinal,TcellNorm,"CCR6")
  GenePlot(TcellMetaFinal,TcellNorm,"IL23R")
  GenePlot(TcellMetaFinal,TcellNorm,"TBX21")
  GenePlot(TcellMetaFinal,TcellNorm,"FOXP3")
  GenePlot(TcellMetaFinal,TcellNorm,"GATA3")
  GenePlot(TcellMetaFinal,TcellNorm,"CD4")
  GenePlot(TcellMetaFinal,TcellNorm,"CD8A")
  GenePlot(TcellMetaFinal,TcellNorm,"CD8B")
  GenePlot(TcellMetaFinal,TcellNorm,"CXCR3")
  GenePlot(TcellMetaFinal,TcellNorm,"CCL20")
  GenePlot(TcellMetaFinal,TcellNorm,"FURIN")
  GenePlot(TcellMetaFinal,TcellNorm,"COTL1")
  GenePlot(TcellMetaFinal,TcellNorm,"BHLHE40")
  GenePlot(TcellMetaFinal,TcellNorm,"HIF1A")
  GenePlot(TcellMetaFinal,TcellNorm,"RBPJ")
  GenePlot(TcellMetaFinal,TcellNorm,"CCR6")
  GenePlot(TcellMetaFinal,TcellNorm,"IFNG")
  GenePlot(TcellMetaFinal,TcellNorm,"TNF")
  GenePlot(TcellMetaFinal,TcellNorm,"IL2RA")
```

### Figure 4B - UMAP T1T17 Subclusters
```{r,}


GenePlot <- function(Metadata,NormData,Gene){
  GeneUse <- as.data.frame(t(NormData[Gene,]))
  colnames(GeneUse) <- Gene
  GeneUse$CellID <- row.names(GeneUse)
  Metadata <- merge(Metadata,GeneUse,by='CellID')
AssignName <- paste(Gene,"Plot",sep="_")
Metadata <- Metadata[order(Metadata[,Gene]),]
Plot <- ggplot(Metadata,aes(x=UMAP1,y=UMAP2,color=eval(parse(text=Gene)))) + theme(axis.line = element_blank())+ xlab("FA 1") + ylab("FA 2")+ geom_point_rast(size=0.5) + ggtitle(Gene)  + theme(legend.position = "none",panel.grid = element_blank(),axis.text = element_blank(),axis.ticks = element_blank(),panel.background = element_blank()) + scale_color_gradient(low='lightgrey',high='darkred',name=NULL)
assign(AssignName,Plot,envir=.GlobalEnv)

  PlotName <- paste("Directory/",Gene,"_Plot.pdf",sep="")

pdf(file=PlotName,height=5,width=5)
print(Plot)
dev.off()
}

  GenePlot(T1T17Metadata,T1T17cellNorm,"BHLHE40")
  GenePlot(T1T17Metadata,T1T17cellNorm,"RORC")
  GenePlot(T1T17Metadata,T1T17cellNorm,"IL23R")
  GenePlot(T1T17Metadata,T1T17cellNorm,"RORA")
  GenePlot(T1T17Metadata,T1T17cellNorm,"CD4")
  GenePlot(T1T17Metadata,T1T17cellNorm,"IL7R")
  GenePlot(T1T17Metadata,T1T17cellNorm,"PDE4D")
  GenePlot(T1T17Metadata,T1T17cellNorm,"TXNIP")
  GenePlot(T1T17Metadata,T1T17cellNorm,"CD8B")
  GenePlot(T1T17Metadata,T1T17cellNorm,"GZMA")
  GenePlot(T1T17Metadata,T1T17cellNorm,"PRF1")
  GenePlot(T1T17Metadata,T1T17cellNorm,"GZMK")
  GenePlot(T1T17Metadata,T1T17cellNorm,"CD8A")
  GenePlot(T1T17Metadata,T1T17cellNorm,"TNF")
  GenePlot(T1T17Metadata,T1T17cellNorm,"IFNG")
  GenePlot(T1T17Metadata,T1T17cellNorm,"IL2RA")
```

### Figure 4C - Dotplot
```{r,}
load("T1T17Seurat.Rdata")
T1T17Genes <- c("CD4","CD8A","CD8B","IL7R","CD6","TXNIP","PDE4D","ZFP36L2","ITGB1","CCR6","CXCR3","GZMA","GZMH","GZMK","GNLY","PRF1","KLRC1","LTA","LTB","IFNG","TNF","CTLA4","GADD45B","SLA","TPI1","PKM","HSPA5","ENO1","DNAJB1","HSPH1")
  
pdf(file="T1T17Dotplot.pdf",height=11,width=8,useDingbats = FALSE)
DotPlot(T1T17Seurat,features=rev(T1T17Genes),cols=c("lightgrey","darkred")) + RotatedAxis() + coord_flip()
dev.off()
```

### Figure 4D - Expxression Overlay T1T17 
```{r,}
GenePlot <- function(Metadata,NormData,Gene){
  GeneUse <- as.data.frame(t(NormData[Gene,]))
  colnames(GeneUse) <- Gene
  GeneUse$CellID <- row.names(GeneUse)
  Metadata <- merge(Metadata,GeneUse,by='CellID')
  Metadata <- Metadata[order(Metadata[,Gene]),]
AssignName <- paste(Gene,"Plot",sep="_")
Plot <- ggplot(Metadata,aes(x=UMAP1,y=UMAP2,color=eval(parse(text=Gene)))) + theme(axis.line = element_blank())+ xlab("FA 1") + ylab("FA 2")+ geom_point_rast(size=0.5) + ggtitle(Gene)  + theme(legend.position = "none",panel.grid = element_blank(),axis.text = element_blank(),axis.ticks = element_blank(),panel.background = element_blank()) + scale_color_gradient(low='lightgrey',high='darkred',name=NULL)
assign(AssignName,Plot,envir=.GlobalEnv)

PlotName <- paste("Directory/",Gene,"_Plot.pdf",sep="")

pdf(file=PlotName,height=5,width=5)
print(Plot)
dev.off()
}

  GenePlot(T1T17Metadata,TcellNormUse,"CD4")
  GenePlot(T1T17Metadata,TcellNormUse,"RORA")
  GenePlot(T1T17Metadata,TcellNormUse,"RORC")
  GenePlot(T1T17Metadata,TcellNormUse,"IL17A")
  GenePlot(T1T17Metadata,TcellNormUse,"IL17F")
  GenePlot(T1T17Metadata,TcellNormUse,"RBPJ")
  GenePlot(T1T17Metadata,TcellNormUse,"IL23R")
  GenePlot(T1T17Metadata,TcellNormUse,"BHLHE40")
  GenePlot(T1T17Metadata,TcellNormUse,"TBX21")
  GenePlot(T1T17Metadata,TcellNormUse,"TBLYM")
  GenePlot(T1T17Metadata,TcellNormUse,"CCR6")
  GenePlot(T1T17Metadata,TcellNormUse,"CXCR3")
```

### Figure 4E - Scatterplot of T1T17 Proportions
```{r,}
T1T17_1_Plot <- ggplot(Proportions_T1T17_0,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Sub-Population 1")

T1T17_2_Plot <- ggplot(Proportions_T1T17_1,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+xlim(0,5) + theme_classic2()+ theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 2)") + ggtitle("Sub-Population 2")

T1T17_3_Plot <- ggplot(Proportions_T1T17_3,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ xlim(0,5) + theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Sub-Population 3")

T1T17_4_Plot <- ggplot(Proportions_T1T17_4,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ xlim(0,5) + theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("Sub-Population 4")

T1T17_Total_Plot <- ggplot(Proportions,aes(x=as.numeric(as.character(logCFU)),y=Frequency)) +geom_point_rast(size=5.0) + stat_smooth(method="lm", se=FALSE,col="red",size=0.5)+ theme_classic2()+ xlim(0,5) + theme(axis.title.y = element_text(size=10))+ theme(legend.position='none',panel.border = element_rect(fill = NA, colour = "black"),panel.background = element_rect(fill = NA, color = "black"))+ ylab("Proportion of Total cells") +xlab("log(CFU + 1)") + ggtitle("T1T17 Total")
```

### Figure 5
```{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
```

### Figure 6
### Figure 6A - Correlation heatmap of cell type proportions
```{r,}
load("Proportions.Rdata")
CellTypeCor_Spearman <- cor(t(Proportions),method='spearman')
CellTypeCor_Pearson <- cor(t(Proportions),method='pearson')

library(pheatmap)

#Cut 5
pdf(file="Spearman_Heatmap.pdf",height=7,width=7)
pheatmap(CellTypeCor_Spearman,cutree_rows = 5,cutree_cols = 5)
dev.off()
pdf(file="Pearson_Heatmap.pdf",height=7,width=7)
pheatmap(CellTypeCor_Pearson,cutree_rows = 5,cutree_cols = 5)
dev.off()

```

### Figure 6B - Barplot of cell-type clusters
```{r,}
library(pheatmap)

Test <- pheatmap(CellTypeCor_Pearson,cutree_rows = 5,cutree_cols = 5)
Groups <- sort(cutree(Test$tree_row, k=5))
Group1 <- names(Groups[Groups == 1])
Group2 <- names(Groups[Groups == 2])
Group3 <- names(Groups[Groups == 3])
Group4 <- names(Groups[Groups == 4])
Group5 <- names(Groups[Groups == 5])


ProportionsGroups1 <- as.data.frame(Proportions[Group1,])
ProportionsGroups2 <- as.data.frame(Proportions[Group2,])
ProportionsGroups3 <- as.data.frame(Proportions[Group3,])
ProportionsGroups4 <- as.data.frame(Proportions[Group4,])
ProportionsGroups5 <- as.data.frame(Proportions[Group5,])

ProportionsGroups1 <- colSums(ProportionsGroups1)
ProportionsGroups2 <- colSums(ProportionsGroups2)
ProportionsGroups3 <- colSums(ProportionsGroups3)
ProportionsGroups4 <- colSums(ProportionsGroups4)
ProportionsGroups5<- colSums(ProportionsGroups5)


CollapsedProrpotions <- as.data.frame(rbind(ProportionsGroups1,ProportionsGroups2,ProportionsGroups3,ProportionsGroups4,ProportionsGroups5))

GranulomaOrder <- as.character(GranulomaMetadata$Granuloma[order(GranulomaMetadata$CFU)])
CollapsedProrpotions <- CollapsedProrpotions[,GranulomaOrder]

HighGranulomas <- as.character(GranulomaMetadata$Granuloma[GranulomaMetadata$HighLow == 1])
LowGranulomas <- as.character(GranulomaMetadata$Granuloma[GranulomaMetadata$HighLow == 0])

CollapsedProrpotionsHigh <- CollapsedProrpotions[,HighGranulomas]
CollapsedProrpotionsLow <- CollapsedProrpotions[,LowGranulomas]

CollapsedProrpotionsHighLow <- cbind(rowMeans(CollapsedProrpotionsLow),rowMeans(CollapsedProrpotionsHigh))

pdf(file="FinalFigures/Figure6 - ReceptorLigand/HierarchicalCompositionBarplotHighLow_Cut5.pdf",height=5,width=10)
print(barplot(as.matrix(CollapsedProrpotionsHighLow)))
dev.off()

```
