SupplementalMethods_COGs.Rmd

---
title: \"Genomic Stability and Genetic Defense Systems in *Dolosigranulum pigrum* a Candidate Beneficial Bacterium from the Human Microbiome\" Supplemental Methods
output: rmarkdown::github_document
---

```{r message=FALSE, include=FALSE}
library(kableExtra)
library(tidyverse)
library(readr)
library(forcats)
library(rstatix)
library(ggpubr)
library(cowplot)
library(grid)
library(RColorBrewer)
knitr::opts_chunk$set(message = FALSE)
```

# COG ANALYSIS

## Data Import

We import the output of `anvi-summarize` and select the most relevant variables for the functional analysis:

```{r message=FALSE}
DpigPangenome <- read_delim("analysis_Anvio7/Pangenomic_Results_Dpig/Dpig-PAN-SUMMARY/PAN_DPIG_prokka_gene_clusters_summary.txt.gz", "\t")
DpigPangenome <- DpigPangenome %>%
  select(-functional_homogeneity_index, -geometric_homogeneity_index, -combined_homogeneity_index, -aa_sequence)
```

In the new variable "accessory_vs_core" we define "Soft/Core" as "MC_core"+"SC_core"+"soft_core" and "accessory" as "shell"+"cloud":

```{r}
DpigPangenome <- DpigPangenome %>%
  mutate(accessory_vs_core=ifelse(grepl("Core", bin_name), "Core", "Accessory"))
```

"Core" is used in the code to avoid problems with the "/" symbol and later replaced with "Soft/Core" for plotting.

```{r echo=FALSE, message=FALSE}
vAccesory <- nrow(DpigPangenome %>% group_by(gene_cluster_id) %>% filter(accessory_vs_core =="Accessory") %>% summarise)
vCore <- nrow(DpigPangenome %>% group_by(gene_cluster_id) %>% filter(accessory_vs_core =="Core") %>% summarise)

vAccesory.p <- round(100*vAccesory/(vAccesory+vCore), 1)
vCore.p <- round(100*vCore/(vAccesory+vCore), 1)
```

There are `r vAccesory` gene clusters (GC) (`r vAccesory.p`%) in the "Accessory" vs. `r vCore` (`r vCore.p`%) in the "Soft/Core" at the pangenome level

## COG Analysis at the Gene Level

We define a new variable `COGs` to use in the plots. This variable is based on `COG20_CATEGORY` but with a cleaner definition of unclassified, uninformative, or mixed assignments:

-   COG categories "Function Unknown" and "General function predictions only" were considered as "Uninformative".
-   If the COG category is mix (e.g., G\|S\|M) it gets labeled as "Ambiguous".
-   Missing values (NA) are labeled as "Unclassified".

```{r, message=FALSE}
DpigPangenome$COGs <- DpigPangenome$COG20_CATEGORY_ACC
DpigPangenome$COGs[DpigPangenome$COGs =="S"]<- "Uninformative"
DpigPangenome$COGs[DpigPangenome$COGs =="R"]<- "Uninformative"
DpigPangenome$COGs[grepl('|', DpigPangenome$COGs,fixed=TRUE)]<-"Ambiguous"
DpigPangenome$COGs[is.na(DpigPangenome$COGs)]<-"Unclassified"
```

Summary of COG annotated genes:

```{r echo=FALSE}
TableGene <- data.frame(
  "Genes" = c("Total in Dpig Pangenome", 
              "COG Category Uninformative = Function Unknown", 
              "COG Category Uninformative = General function prediction only",
              "COG Category Ambiguous (Mixed COG Category)",
              "COG Category Unclassified (Non-assigned)", 
              "Informative COGs (Total - Uninformative, Ambiguous & Unclassified)"),
  "Count" = c(nrow(DpigPangenome), 
              nrow(DpigPangenome %>% filter(COG20_CATEGORY_ACC =="S")),
              nrow(DpigPangenome %>% filter(COG20_CATEGORY_ACC =="R")),
              nrow(DpigPangenome %>% filter(COGs =="Ambiguous")),
              nrow(DpigPangenome %>% filter(COGs =="Unclassified")), 
              nrow(DpigPangenome %>% filter(COGs !="Uninformative" & COGs !="Ambiguous" & COGs !="Unclassified"))
              )
)
TableGene$Percentage <- round(100*(TableGene$Count/49412),1)

kable(TableGene)
```

63.8% of the gene calls are Informative.

## COG Analysis at the Gene Cluster Level

This analysis was done at the pangenomic gene cluster level (GC). Since many gene clusters had mixed COG category assignments a solution is to assign each individual gene call to their corresponding Genome/Accessory_vs_Core/COG grouping weighting their contribution by dividing their count by the number of genes in their GC.

### GCs by COG Category and Genome

The table "GCsbyCOG_Genome" groups the genes by genome; and inside genomes by "Accessory" vs. "Soft/Core" status, and nested inside as the COG category. But, in this case, instead of counting the elements in each group we calculated the sum of 1/`num_genes_in_gene_cluster`.

```{r, message=FALSE}
GCsbyCOG_Genome <- DpigPangenome %>%
  group_by(genome_name, accessory_vs_core, COGs) %>%
  summarise(num_corrected_genes=sum(1/num_genes_in_gene_cluster))
```

The total sum of all values in the `num_corrected_genes` variable should add up to the number of CGs:

```{r, message=FALSE}
sum(GCsbyCOG_Genome$num_corrected_genes)
nrow(DpigPangenome %>% group_by(gene_cluster_id) %>% summarise)
```

Adding extra column to label the gray scale portion of the plots:

```{r}
GCsbyCOG_Genome <- GCsbyCOG_Genome %>%
  mutate(Assignment=ifelse(COGs!="Uninformative" & COGs!="Ambiguous" & COGs!="Unclassified", "Informative", COGs))
```

#### Summary of GOC annotated GCs in the Accessory vs. Soft/Core :

```{r message=FALSE, warning=TRUE}
TableGC <- GCsbyCOG_Genome %>% 
  group_by(accessory_vs_core, Assignment) %>%
  summarize(corrected_genes=sum(num_corrected_genes))

TableGC$Percentages <- round(100*TableGC$corrected_genes/sum(TableGC$corrected_genes), 1)

kable(TableGC)
```

#### Summary of GOC annotated GCs in the Accessory:

```{r message=FALSE, warning=TRUE}
TableGCAccessory <- GCsbyCOG_Genome %>% 
  filter(accessory_vs_core =="Accessory") %>%
  group_by(accessory_vs_core, Assignment) %>%
  summarize(corrected_genes=sum(num_corrected_genes))

TableGCAccessory$Percentages <- round(100*TableGCAccessory$corrected_genes/sum(TableGCAccessory$corrected_genes), 1)

kable(TableGCAccessory)
```

#### Summary of GOC annotated GCs in the Soft/Core:

```{r message=FALSE, warning=TRUE}
TableGCCore <- GCsbyCOG_Genome %>% 
  filter(accessory_vs_core =="Core") %>%
  group_by(accessory_vs_core, Assignment) %>%
  summarize(corrected_genes=sum(num_corrected_genes))

TableGCCore$Percentages <- round(100*TableGCCore$corrected_genes/sum(TableGCCore$corrected_genes), 1)

kable(TableGCCore)
```

#### Summary of GOC annotated GCs by Genome in the Accessory vs. Soft/Core :

```{r message=FALSE, warning=TRUE}
TableGenomes <- GCsbyCOG_Genome %>% 
  group_by(genome_name, accessory_vs_core) %>% 
  summarize(corrected_genes=sum(num_corrected_genes))

kable(TableGenomes)
```

#### Renaming and ordering variables factor levels for plotting:

```{r}
GCsbyCOG_Genome$accessory_vs_core <- factor(GCsbyCOG_Genome$accessory_vs_core, levels =c("Core", "Accessory"))

GCsbyCOG_Genome$COGs <- recode_factor(GCsbyCOG_Genome$COGs, "Q"="Secondary metabolites biosynthesis, transport, and catabolism","P"="Inorganic ion transport and metabolism","I"="Lipid transport and metabolism","H"="Coenzyme transport and metabolism","G"="Carbohydrate transport and metabolism","F"="Nucleotide transport and metabolism","E"="Amino acid transport and metabolism","C"="Energy production and conversion","X"="Mobilome: prophages, transposons","L"="Replication, recombination and repair","K"="Transcription","J"="Translation, ribosomal structure and biogenesis","V"="Defense mechanisms","U"="Intracellular trafficking, secretion, and vesicular transport","T"="Signal transduction mechanisms","O"="Post-translational modification, protein turnover, and chaperones","N"="Cell Motility","M"="Cell wall/membrane/envelope biogenesis","D"="Cell cycle control, cell division, chromosome partitioning","Uninformative"="Uninformative","Ambiguous"="Ambiguous","Unclassified"="Unclassified", .ordered = TRUE)

GCsbyCOG_Genome$Assignment <- recode_factor(GCsbyCOG_Genome$Assignment,  "Informative"=" ", "Uninformative"="Uninformative", "Ambiguous"="Ambiguous", "Unclassified"="Unclassified", .ordered = TRUE)

GCsbyCOG_Genome$genome_name <- recode_factor(GCsbyCOG_Genome$genome_name, "ATCC_51524"="ATCC 51524", "KPL3250"="KPL3250", "KPL1939_CDC4792_99"="CDC 4792-99","KPL1934_CDC4709_98"="CDC 4709-98", "KPL1922_CDC39_95"="CDC 39-95", "KPL3264"="KPL3264", "KPL3256"="KPL3256", "KPL3033"="KPL3033", "KPL1933_CDC4545_98"="CDC 4545-98", "KPL1930_CDC2949_98"="CDC 2949-98", "KPL3069"="KPL3069", "KPL3052"="KPL3052", "KPL3090"="KPL3090", "KPL3086"="KPL3086", "KPL3065"="KPL3065", "KPL3043"="KPL3043", "KPL3911"="KPL3911", "KPL3084"="KPL3084", "KPL3070"="KPL3070",
"KPL3246"="KPL3246", "KPL1937_CDC4199_99"="CDC 4199-99","KPL3274"="KPL3274","KPL3050"="KPL3050","KPL1938_CDC4791_99"="CDC 4791-99", "KPL1932_CDC4420_98"="CDC 4420-98", "KPL3077"="KPL3077", "KPL1931_CDC4294_98"="CDC 4294-98", "KPL1914"="KPL1914", .ordered = TRUE)
```

### GCs by COG Category

The table "GCsbyCOG" groups the genes by "Accessory" vs. "Soft/Core" status, and nested inside as the COG category.

```{r, message=FALSE}
GCsbyCOG <- DpigPangenome %>%
  group_by(accessory_vs_core, COGs) %>%
  summarise(num_corrected_genes=sum(1/num_genes_in_gene_cluster))
```

#### Renaming and ordering variables factor levels for plotting:

```{r}
GCsbyCOG$COGs <- recode_factor(GCsbyCOG$COGs, "Q"="Secondary metabolites biosynthesis, transport, and catabolism",
                               "P"="Inorganic ion transport and metabolism",
                               "I"="Lipid transport and metabolism",
                               "H"="Coenzyme transport and metabolism",
                               "G"="Carbohydrate transport and metabolism",
                               "F"="Nucleotide transport and metabolism",
                               "E"="Amino acid transport and metabolism",
                               "C"="Energy production and conversion",
                               "X"="Mobilome: prophages, transposons",
                               "L"="Replication, recombination and repair",
                               "K"="Transcription",
                               "J"="Translation, ribosomal structure and biogenesis",
                               "V"="Defense mechanisms",
                               "U"="Intracellular trafficking, secretion, and vesicular transport",
                               "T"="Signal transduction mechanisms",
                               "O"="Post-translational modification, protein turnover, and chaperones",
                               "N"="Cell Motility",
                               "M"="Cell wall/membrane/envelope biogenesis",
                               "D"="Cell cycle control, cell division, chromosome partitioning",
                               "Uninformative"="Uninformative",
                               "Ambiguous"="Ambiguous",
                               "Unclassified"="Unclassified", .ordered = TRUE)
```

#### Summary of GOC annotated GCs GCs by COG Category:

New table "GCsbyCOG_CorevsAcc" in wide format. % of each category relative to the "Accessory" or "Soft/Core" was calculated (pTotal. variables). Total GCs for each COG category calculated, and % of GCs in the "Accessory" and "Soft/Core" relative to each category (p. values) were calculated as well. The ratio between the number of GC in the "Accessory" vs. the "Soft/Core" is calculated for each COG:

```{r}
GCsbyCOG_CorevsAcc <- spread(GCsbyCOG, accessory_vs_core, num_corrected_genes)
GCsbyCOG_CorevsAcc$pTotal.Accessory <- round(100*GCsbyCOG_CorevsAcc$Accessory/sum(GCsbyCOG_CorevsAcc$Accessory), 1)
GCsbyCOG_CorevsAcc$pTotal.Core <- round(100*GCsbyCOG_CorevsAcc$Core/sum(GCsbyCOG_CorevsAcc$Core), 1)
GCsbyCOG_CorevsAcc$total <- GCsbyCOG_CorevsAcc$Accessory + GCsbyCOG_CorevsAcc$Core
GCsbyCOG_CorevsAcc$pTotal.total <- round(100*GCsbyCOG_CorevsAcc$total/sum(GCsbyCOG_CorevsAcc$total), 1)
GCsbyCOG_CorevsAcc$p.accessory <- round(100*(GCsbyCOG_CorevsAcc$Accessory/GCsbyCOG_CorevsAcc$total), 1)
GCsbyCOG_CorevsAcc$p.core <- round(100*(GCsbyCOG_CorevsAcc$Core/GCsbyCOG_CorevsAcc$total), 1)
GCsbyCOG_CorevsAcc$ratio <- round(GCsbyCOG_CorevsAcc$Accessory/GCsbyCOG_CorevsAcc$Core, 2)

kable(GCsbyCOG_CorevsAcc)
```

## Plots

Color Palettes

```{r}
getPalette <- colorRampPalette(brewer.pal(8, "Set1"))
CountTotalCOGs <- length(unique(GCsbyCOG_Genome$COGs))

palette1 <- c("grey60", "grey40", "grey20", getPalette(CountTotalCOGs-3)) # 22 elements: Colors + Grays
palette2 <- getPalette(CountTotalCOGs-3) # 19 elements: Colors
palette3 <- c("grey60", "grey40", "grey20", "white") # 4 elements: White + Grays
```

### Plots Accessory vs. Soft/Core

Panel A in main figure:

```{r, TotalGCs.accessory_vs_core}
pA <- ggplot(GCsbyCOG_Genome, aes(x = accessory_vs_core, y = num_corrected_genes, fill = fct_rev(COGs))) +
  stat_summary(fun=sum ,geom="bar", position = "stack") +
  scale_x_discrete(labels = c("Soft/Core", "Accessory")) +
  scale_fill_manual(values = palette1) +
  scale_y_continuous(expand = c(0,0), breaks=seq(0, 1500, by = 250)) +
  labs(fill="COG Categories", x=" ", y= "Number of Gene Clusters") +
  theme_classic() +
  theme(axis.title = element_text(size = 9), axis.text = element_text(size=7), plot.margin=unit(c(10,0,10,20),"pt"), legend.position = "none") 
```

```{r, echo=FALSE,fig.height=4, fig.width=4}
pA
```

This plot is used for the grayscale legend:

```{r, InformativeGCs.accessory_vs_core}
pB <- ggplot(GCsbyCOG_Genome, aes(x = accessory_vs_core, y = num_corrected_genes, fill = fct_rev(Assignment))) +
  stat_summary(fun=sum ,geom="bar", position = "stack") +
  scale_x_discrete(labels = c("Soft/Core", "Accessory")) +
  scale_fill_manual(values = palette3) +
  scale_y_continuous(expand = c(0,0), breaks=seq(0, 1500, by = 250)) +
  labs(fill=" ", x=" ", y= "Number of Gene Clusters") +
  theme_classic() +
  theme(axis.title = element_text(size = 9), axis.text = element_text(size=7), plot.margin=unit(c(10,0,10,20),"pt"), legend.key.size = unit(0.7, "line"), legend.text = element_text(size = 7), legend.box.margin = margin(10,20,10,10)) +
  guides(fill=guide_legend(ncol=1, title.position = "top", title.hjust = 0.5))
```

### Plots by Genome

Panel A in supplemental figure:

```{r, TotalGCs.accesory.byGenome}
pC <- ggplot(filter(GCsbyCOG_Genome, accessory_vs_core == "Accessory"), aes(x=genome_name, y=num_corrected_genes, fill = fct_rev(COGs))) +
  stat_summary(fun=sum ,geom="bar", position = "stack") +
  scale_fill_manual(values = palette1) +
  scale_y_continuous(expand = c(0,0)) + 
  labs(fill="COG Assignment", x="", y= "Number of Gene Clusters") +
  theme_classic() + 
  theme(axis.text.y = element_text(size=7), axis.text.x = element_text(size=8, angle=75, hjust=1)) +
  theme(legend.position = "none", plot.margin=unit(c(15,15,-10,20),"pt")) 
```

```{r, echo=FALSE,fig.height=4, fig.width=8}
pC
```

Panel B in supplemental figure:

```{r, InformativeGCs.accesory.byGenome}
pD <- ggplot(filter(GCsbyCOG_Genome %>% filter(COGs != "Uninformative", COGs !="Ambiguous", COGs != "Unclassified"), accessory_vs_core == "Accessory"), aes(x=genome_name, y=num_corrected_genes, fill = fct_rev(COGs))) +
  stat_summary(fun=sum ,geom="bar", position = "stack") +
  scale_y_continuous(expand = c(0,0)) + 
  scale_fill_manual(values = palette2) + 
  labs(fill="COG Categories", x="", y= "Number of Informative Gene Clusters") +
  theme_classic() + 
  theme(axis.text.y = element_text(size=7), axis.text.x = element_text(size=8, angle=75, hjust=1)) +
  theme(legend.position="bottom", legend.key.size = unit(0.7, "line"), legend.text = element_text(size = 8), plot.margin=unit(c(0,15,0,20),"pt")) +
  guides(fill=guide_legend(ncol=2, title.position = "top", title.hjust = 0.5)) 
```

```{r, echo=FALSE,fig.height=6.5, fig.width=8}
pD
```

This plot is used for the grayscale legend:

```{r, TotalGCs.accesory.byGenome.legend}
pE <- ggplot(filter(GCsbyCOG_Genome, accessory_vs_core == "Accessory"), aes(x=genome_name, y=num_corrected_genes, fill = fct_rev(Assignment))) +
  stat_summary(fun=sum ,geom="bar", position = "stack") +
  scale_fill_manual(values = palette3) +
  scale_y_continuous(expand = c(0,0)) + 
  labs(fill="Accessory Genome COG Assignment", x="", y= "Number of Gene Clusters") +
  theme_classic() + 
  theme(axis.text.y = element_text(size=7), axis.text.x = element_text(size=8, angle=75, hjust=1)) +
  theme(legend.position="bottom", legend.key.size = unit(0.7, "line"), legend.text = element_text(size = 8), legend.title = element_text(face="bold", size = 12), plot.margin=unit(c(15,15,0,20),"pt")) +
  guides(fill=guide_legend(nrow=1, title.position = "top", title.hjust = -5))
```

This is used for the clade labels in supplemental figure:

```{r}
pclades <- ggplot() +
  scale_y_continuous(limits = c(-1.5, 0.5), breaks = c(-1, 0)) +
  geom_segment(aes(x=0,xend=2.9,y=0,yend=0), color="#2c9b51ff") +
  geom_segment(aes(x=3,xend=4.9,y=0,yend=0), color="#a851a8ff") +
  geom_segment(aes(x=5,xend=9.9,y=0,yend=0), color="#e17139ff") +
  geom_segment(aes(x=10,xend=28,y=0,yend=0), color="#1b1d86ff") +
  annotate("text", x = 1.5, y = -1, label = "C1", fontface="bold", color="#2c9b51ff")+
  annotate("text", x = 4, y = -1, label = "C2", fontface="bold", color="#a851a8ff")+
  annotate("text", x = 7.3, y = -1, label = "C3", fontface="bold", color="#e17139ff")+
  annotate("text", x = 19, y = -1, label = "C4", fontface="bold", color="#1b1d86ff")+
  theme_classic() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_blank(), plot.margin=unit(c(0,0,5,20),"pt")) 
```

### Plots by COG Category

In order to represent the Soft/Core on the left of the plot with absolute values per COG category we create `core.neg`; a negative version of the `core` variable in GCsbyCOG_CorevsAcc. Table converted to the long format for plotting.

```{r}
GCsbyCOG_CorevsAcc$core.neg <- -GCsbyCOG_CorevsAcc$Core
GCsbyCOG_CorevsAccLong <- gather(GCsbyCOG_CorevsAcc, accessory_vs_core, plotting, core.neg, Accessory)
```

Panel B in main figure:

```{r, InformativeGCs.byCOG}
pF <- ggplot(filter(GCsbyCOG_CorevsAccLong, COGs != "Uninformative", COGs != "Ambiguous", COGs != "Unclassified"), aes(x = COGs, y = plotting, fill = COGs)) +
  geom_bar(stat="identity") + 
  scale_fill_manual(values = rev(palette2)) + 
  scale_x_discrete(position = "top") +
  labs(title= "COG Categories", x="", y= "Number of Gene Clusters") +
  coord_flip() +
  scale_y_continuous(limits = c(-180, 180), breaks = c(-150, -100, -50, 0, 50, 100, 150), label = c(150, 100, 50, 0, 50, 100, 150)) +
  geom_segment(aes(x=0,xend=19.5,y=0,yend=0), linetype=3, size=0.1) +
  geom_label(aes(x = 20.5, y = -95, label = "      Soft/Core       "), fontface="bold", size=3, fill = "grey90", label.size=NA, label.padding = unit(0.3, "lines")) +
  geom_label(aes(x = 20.5, y = 95, label = "     Accessory      "), fontface="bold", size=3, fill = "grey90", label.size=NA, label.padding = unit(0.3, "lines")) +
  theme_classic() +
  theme(axis.title = element_text(size = 9), axis.text.x = element_text(size=7), axis.ticks.y = element_blank(), axis.line.y = element_blank(), legend.position = "none", plot.margin=unit(c(5,10,10,25),"pt"), plot.title=element_text(face="bold", hjust=3, vjust=-3.9)) 

gpF <- ggplotGrob(pF)
gpF$layout$clip[gpF$layout$name=="panel"] <- "off"
```

```{r, echo=FALSE,fig.height=4, fig.width=8}
pF
```

## Final Figures

### Main figure

```{r eval=TRUE}
pMain <- ggarrange(ggarrange(get_legend(pB), pA, ncol = 1, heights = c(0.2, 1)),
                   gpF, ncol = 2, labels = c("A", "B"), hjust=-0.5, vjust=2, widths = c(0.7, 2))

ggsave("analysis_COGs/Fig4_COG_summary.tiff", pMain, width = 9, height = 4, dpi = 300, units="in")
```

```{r, echo=FALSE,fig.height=4, fig.width=9}
pMain
```

### Supplemental Figure

```{r eval=TRUE}
pSupple <- ggarrange(get_legend(pE),
                      ggarrange(pC+theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()), pD+theme(legend.position="none"), ncol = 1,  align = "v", labels = c("i", "ii"), hjust=-0.5, vjust=1, heights = c(1, 1)),
                      pclades, 
                      get_legend(pD), ncol = 1, heights = c(0.2, 2, 0.2, 0.6))

ggsave("analysis_COGs/FigS1D_COG_byGenome.tiff", pSupple, width = 8, height = 10, dpi = 150)
```

```{r, echo=FALSE,fig.height=10, fig.width=8}
pSupple
```