---
title: "Relation of pest insect-killing and soilborne pathogen-inhibition abilities to species diversification in environmental Pseudomonas protegens. ISME J 17, 1369–1381 (2023)."
author: "Daniel Garrido-Sanz (daniel.garridosanz@unil.ch)"
affiliation: "Keel Lab | Department of fundamental microbiology | University of Lausanne"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: 
  rmdformats::html_clean:
    code_folding: hide
    toc: true
---

Code used in: Garrido-Sanz, D., Vesga, P., Heiman, C.M. et al. Relation of pest insect-killing and soilborne pathogen-inhibition abilities to species diversification in environmental Pseudomonas protegens. ISME J 17, 1369–1381 (2023). https://doi.org/10.1038/s41396-023-01451-8

#1. Genome fractions

##1.1 Load packages

```{r}
library("dplyr")
library("matrixStats")
library("ggplot2")
library(gridExtra)
library(reshape2)
library(tidyverse)
```

##1.2 Load Data
We only need one of the files produced by OrthoFinder analysis: "Orthogroups.GeneCount.tsv"
This file DOES NOT CONTAIN the singletons, which are present in the Orthogroups.txt. Manually add in the Orthogroups.Gene.tsv the Orthogroups that are only present in one genome (singletons). Save it as Orthologs_and_singletons.txt.

```{r}
orthologs <- read.table("Orthologs_and_singletons.txt", header=T, sep="\t")
rownames(orthologs)<-orthologs[, 1] 
orthologs<-orthologs[,-1]
#Convert all values greater than 1 into 1 (i.e., removing paralogs from the calculation)
orthologs[orthologs > 1] = 1
orthologs <- orthologs %>% mutate_all(as.integer)
head(orthologs)
```


##1.3 Index for random sampling

```{r}
strainsv = as.vector(colnames(orthologs))
strainsdf = NULL
strainsdf = as.data.frame(colnames(orthologs))

for (i in 1:100) { #100 Replicates
  strainsdf = cbind(strainsdf, sample(strainsv))
  tag = paste("model", i, sep="_")
  colnames(strainsdf)[i+1] <- tag
}
head(strainsdf)
#write.table(strainsdf, file="sampling_table.txt", sep="\t") #Save it into a file in case we need it latter.
```


##1.4 Genome fractions
The core genome is composed of OGs (orthogroups) present in a number X of strains. 
The pangenome is the total count of unique OGs.
The strain-specific genome is the novel OGs that appear over sampled genomes.
We are going to randomly sample all genomes and ask about how many OGs are present to obtain these fractions.


```{r}
outtable_core<-data.frame(a=1:length(rownames(strainsdf)))
outtable_core$a<- NULL

outtable_pan<-data.frame(a=1:length(rownames(strainsdf)))
outtable_pan$a<- NULL

outtable_strain<-data.frame(a=1:length(rownames(strainsdf)))
outtable_strain$a<- NULL

for (i in 1:length(strainsdf)) {
  print(paste("Loop", as.character(i), "of ", as.character(length(strainsdf))))
  coretable <- data.frame(a=1:length(rownames(orthologs)))
  rownames(coretable) <- rownames(orthologs)
  coretable$a<-NULL
  coreevol <- vector()
  panevol <- vector()
  strainspecevol <- vector()
  modelog <- as.vector(strainsdf[,i])
  for (strain in modelog) {
   coretable<- bind_cols(coretable, select(orthologs, strain))
   core<-nrow(filter(coretable, rowSums(coretable)==length(colnames(coretable))))
   coreevol <- c(coreevol, core)
   pan<-nrow(filter(coretable, rowSums(coretable)>0))
   panevol <- c(panevol, pan)
   strainspec<-nrow(coretable %>% filter((!!as.symbol(strain)) == 1 & rowSums(coretable)==1))
   strainspecevol <- c(strainspecevol, strainspec)
  }
  outtable_core <- cbind(outtable_core,coreevol)
  outtable_pan <- cbind(outtable_pan,panevol)
  outtable_strain <- cbind(outtable_strain,strainspecevol)
}
colnames(outtable_core)<-colnames(strainsdf)
colnames(outtable_pan)<-colnames(strainsdf)
colnames(outtable_strain)<-colnames(strainsdf)

#write.table(outtable_core, "CORE_EVOL_100.txt", quote = F, sep = "\t")
#write.table(outtable_pan, "PAN_EVOL_100.txt", quote = F, sep = "\t")
#write.table(outtable_strain, "STRAIN_EVOL_100.txt", quote = F, sep = "\t")
```

Load output again if restarting

```{r}
#outtable_core = read.table("All_Fractions/CORE_EVOL_100.txt", header = T, sep = "\t")
#outtable_pan = read.table("All_Fractions/PAN_EVOL_100.txt", header = T, sep = "\t")
#outtable_strain = read.table("All_Fractions/STRAIN_EVOL_100.txt", header = T, sep = "\t")
```

##1.5 Fit models

For fitting the models, use the output already written (CORE_EVOL_100.txt, PAN_EVOL_100.txt, SPECIFIC_EVOL_100.txt). Transform them into the output of get_homologues and use the plot_pancore_matrix.pl to obtain the fitted models. These values will be used to calculate the extrapolation to see if we expect to see an stabilization soon.

```{r}
corefit_will = function(genomes) {
  coregenes_will = -3669 + 10924*exp((-sqrt(sqrt(genomes))/ 6.28))
  return(coregenes_will)
}

corefit_tett = function(genomes) {
  coregenes_tett = 2737 + 2423 * exp((-genomes/41.31))
  return(coregenes_tett)
}

panfit = function(g) {
  pangenes = 6595 + 43.0 *(g - 1) + 302 *exp(-2/15.43)*(1 - exp(-(g - 1)/ 15.43))/(1 - exp(-1/ 15.43)) 
  return(pangenes)
}

specific_fit = function(g) {
  specificgenes = 6078 - 1.4(g-1) - 434411 * exp(.2/0.45)* (1-exp(-(g-1)/0.45))/1-exp(-1/0.45)
  return(specificgenes)
}
```

##1.6 Plots

```{r}
core_melt = outtable_core
core_melt$index = rownames(core_melt)
core_melt = melt(core_melt, id.vars = "index")

#Filter outliers for the jitter
core_melt = core_melt %>% group_by(index) %>% mutate(value_filtered = case_when(as.numeric(value) - quantile(as.numeric(value))[4] > 1.5*IQR(as.numeric(value)) ~ NA_real_,
                                  quantile(as.numeric(value))[2] - as.numeric(value) > 1.5*IQR(as.numeric(value)) ~ NA_real_,
                                  TRUE ~ as.numeric(value)))


core_mean = core_melt %>% group_by(index) %>% drop_na() %>% summarise(mean = mean(value_filtered), sd = sd(value_filtered))

core = ggplot(core_melt, mapping = aes(x = reorder(index, as.numeric(index)), y = value)) +
  geom_jitter(aes(x = reorder(index, as.numeric(index)), y = value_filtered), size = 0.7, alpha= 0.4, color= "dodgerblue4", width = 0.3, shape = 16)+
  geom_boxplot(outlier.size = 0.5, outlier.alpha = 0.25, outlier.colour = "red", alpha = 0.2, fill="dodgerblue4", lwd=0.5, outlier.shape = 16) +
  stat_summary(data = core_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "line", color="dodgerblue2", group=1) +
  stat_summary(data = core_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "point",  color="dodgerblue2", size = 0.9, group=1, shape = 16) +
  geom_hline(yintercept = min(core_melt$value_filtered, na.rm =  T), color="red", linetype = "dashed") +
  theme_classic() + 
  scale_y_continuous(breaks = c(2756, 3000,4000,5000,6000,7000,8000,9000,10000,11000,12000,13000,14000,15000,16000), limits = c(2750,7000)) +
  scale_x_discrete(breaks = seq(from = 0, to = 115, by = 5)) +
  stat_function(fun = corefit_tett, color = "red") +
  stat_function(fun = corefit_will, color = "blue") +
  ylab("orthologous groups") + xlab("Sampled genomes")
```

```{r}
pan_melt = outtable_pan
pan_melt$index = rownames(pan_melt)
pan_melt = melt(pan_melt, id.vars = "index")

#Filter outliers for the jitter
pan_melt = pan_melt %>% group_by(index) %>% mutate(value_filtered = case_when(as.numeric(value) - quantile(as.numeric(value))[4] > 1.5*IQR(as.numeric(value)) ~ NA_real_,
                                  quantile(as.numeric(value))[2] - as.numeric(value) > 1.5*IQR(as.numeric(value)) ~ NA_real_,
                                  TRUE ~ as.numeric(value)))

pan_mean = pan_melt %>% group_by(index) %>% drop_na() %>% summarise(mean = mean(value_filtered), sd = sd(value_filtered))

pange = ggplot(pan_melt, mapping = aes(x = reorder(index, as.numeric(index)), y = value)) +
  geom_jitter(aes(x = reorder(index, as.numeric(index)), y = value_filtered), size = 0.7, alpha= 0.4, color= "springgreen4", width = 0.3, shape = 16)+
  geom_boxplot(outlier.size = 0.5, outlier.alpha = 0.25, outlier.colour = "red", alpha = 0.2, fill="springgreen4", lwd=0.5, shape = 16) +
  stat_summary(data = pan_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "line", color="springgreen2", group=1) +
  stat_summary(data = pan_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "point",  color="springgreen2", size = 0.9, group=1, shape = 16) +
  geom_hline(yintercept = max(pan_melt$value_filtered, na.rm =  T), color="red", linetype = "dashed") +
  theme_classic() +
  scale_y_continuous(breaks = c(3000,4000,5000,6000,7000,8000,9000,10000,11000,12000,13000,14000,15000,15620, 16000), limits = c(5000,16000)) +
  scale_x_discrete(breaks = seq(from = 0, to = 115, by = 5)) +
  stat_function(fun = panfit, color = "red") +
  ylab("orthologous groups") + xlab("genomes sampled")
```

```{r}
specific_melt = outtable_strain
specific_melt$index = rownames(specific_melt)
specific_melt = melt(specific_melt, id.vars = "index")

specific_melt = specific_melt %>% group_by(index) %>% mutate(value_filtered = case_when(as.numeric(value) - quantile(as.numeric(value))[4] > 1.5*IQR(as.numeric(value)) ~ NA_real_,
                                  quantile(as.numeric(value))[2] - as.numeric(value) > 1.5*IQR(as.numeric(value)) ~ NA_real_,
                                  TRUE ~ as.numeric(value)))

specific_mean = specific_melt %>% group_by(index) %>% drop_na() %>% summarise(mean = mean(value_filtered), sd = sd(value_filtered))

ggplot(specific_melt, mapping = aes(x = reorder(index, as.numeric(index)), y = value)) +
  geom_jitter(aes(x = reorder(index, as.numeric(index)), y = value_filtered), size = 0.7, alpha= 0.4, color= "darkorchid4", width = 0.3,shape = 16)+
  geom_boxplot(outlier.size = 0.5, outlier.alpha = 0.25, outlier.colour = "red", alpha = 0.2, fill="darkorchid4", lwd=0.5,shape = 16) +
  stat_summary(data = specific_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "line", color="darkorchid2", group=1) +
  stat_summary(data = specific_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "point",  color="darkorchid2", size = 0.9, group=1) +
  scale_x_discrete(breaks = seq(from = 0, to = 115, by = 5)) +
  theme_classic() + 
  ylab("orthologous groups") + xlab("genomes sampled")

s1 = ggplot(specific_melt, mapping = aes(x = reorder(index, as.numeric(index)), y = value)) +
  geom_jitter(aes(x = reorder(index, as.numeric(index)), y = value_filtered), size = 0.7, alpha= 0.4, color= "darkorchid4", width = 0.3,shape = 16)+
  geom_boxplot(outlier.size = 0.5, outlier.alpha = 0.25, outlier.colour = "red", alpha = 0.2, fill="darkorchid4", lwd=0.5,shape = 16) +
  stat_summary(data = specific_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "line", color="darkorchid2", group=1) +
  stat_summary(data = specific_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "point",  color="darkorchid2", size = 0.9, group=1) +
  theme_classic() +
  coord_cartesian(ylim = c(5700, 6800)) +
  scale_y_continuous(breaks=c(6000,6500)) +
  scale_x_discrete(breaks = seq(from = 0, to = 115, by = 5)) +
  xlab("") + ylab("") + 
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
  theme(plot.margin=unit(c(5.5,5.5,-10, 5.5), "pt"), axis.line.y=element_line())

s2 = ggplot(specific_melt, mapping = aes(x = reorder(index, as.numeric(index)), y = value)) +
  geom_jitter(aes(x = reorder(index, as.numeric(index)), y = value_filtered), size = 0.7, alpha= 0.4, color= "darkorchid4", width = 0.3,shape = 16)+
  geom_boxplot(outlier.size = 0.5, outlier.alpha = 0.25, outlier.colour = "red", alpha = 0.2, fill="darkorchid4", lwd=0.5,shape = 16) +
  stat_summary(data = specific_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "line", color="darkorchid2", group=1) +
  stat_summary(data = specific_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "point",  color="darkorchid2", size = 0.9, group=1) +
  theme_classic() + 
  coord_cartesian(ylim = c(0, 1700))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100,120,140,160,180,200,400,600,800,1000,1200,1400,1600)) +
  scale_x_discrete(breaks = seq(from = 0, to = 115, by = 5)) +
  stat_function(fun= specific_fit, color="red")+
  xlab("Sampled genomes") + ylab("Orthologous groups") +
  theme(axis.line.x=element_line(), axis.line.y=element_line())

grid_s = grid.arrange(s1, s2, heights = c(1.5,3))

ggplot(specific_melt, mapping = aes(x = reorder(index, as.numeric(index)), y = value)) +
  geom_jitter(aes(x = reorder(index, as.numeric(index)), y = value_filtered), size = 0.7, alpha= 0.4, color= "darkorchid4", width = 0.3,shape = 16)+
  geom_boxplot(outlier.size = 0.5, outlier.alpha = 0.25, outlier.colour = "red", alpha = 0.2, fill="darkorchid4", lwd=0.5,shape = 16) +
  stat_summary(data = specific_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "line", color="darkorchid2", group=1) +
  stat_summary(data = specific_mean, aes(x = reorder(index, as.numeric(index)), y = mean), geom = "point",  color="darkorchid2", size = 0.9, group=1) +
  theme_classic() + 
  coord_cartesian(ylim = c(0, 200), xlim=c(100,115))+
  scale_y_continuous(breaks=c(0,20,40,60,80,100,120,140,160,180,200)) +
  scale_x_discrete(breaks = seq(from = 0, to = 115, by = 5)) +
  xlab("Sampled genomes") + ylab("Orthologous groups") +
  theme(axis.line.x=element_line(), axis.line.y=element_line())
```

```{r}
grid_final = grid.arrange(core, pange, s2, ncol = 3)
```


##1.7 Genome fractions per group
Load table with the subset of strains belonging to a specific species cluster.

```{r}
clusters = read.table("clusters_no.txt", header = T, sep = "\t")
head(clusters)
```

Calculate fractions

```{r}
SC.results_core = data.frame()
SC.results_pan = data.frame()
SC.results_strain = data.frame()

for (i in sort(unique(clusters$SpeciesCluster))) { #iterate through the different species clusters
  
  #1. create OG subtable with columns of the species cluster
  SC.names = clusters %>% select(-OrderTree, -GenomeName) %>% filter(SpeciesCluster == i)
  SC.og = orthologs[as.vector(SC.names$NumericName)]
  
  #2. create random index of 100 replicates
  SC.strainsv = as.vector(colnames(SC.og))
  SC.strainsdf = NULL
  SC.strainsdf = as.data.frame(colnames(SC.og))
  for (j in 1:100) { #loop for the index
    SC.strainsdf = cbind(SC.strainsdf, sample(SC.strainsv))
    tag = paste("model", j, sep="_")
    colnames(SC.strainsdf)[j+1] <- tag
    
    #3. Create empty dfs to store results
    SC_core<-data.frame(a=1:length(rownames(SC.strainsdf)))
    SC_core$a<- NULL
    
    SC_pan<-data.frame(a=1:length(rownames(SC.strainsdf)))
    SC_pan$a<- NULL
    
    SC_strain<-data.frame(a=1:length(rownames(SC.strainsdf)))
    SC_strain$a<- NULL
  }
    
  #4. Calculate genome fractions
  for (k in 1:length(SC.strainsdf)) { 
    #print(paste("Analyzing", i, "Loop", as.character(k), "of", as.character(length(SC.strainsdf))))
    SC.coretable = data.frame(a=1:length(rownames(SC.og)))
    rownames(SC.coretable) = rownames(SC.og)
    SC.coretable$a=NULL
    SC.coreevol=vector()
    SC.panevol=vector()
    SC.strainevol = vector()
    modelog = as.vector(SC.strainsdf[,k])
    for (strain in modelog) { 
      SC.coretable = bind_cols(SC.coretable, select(SC.og, strain))
      SC.core = nrow(filter(SC.coretable, rowSums(SC.coretable)==length(colnames(SC.coretable))))
      SC.coreevol = c(SC.coreevol, SC.core)
      SC.pan = nrow(filter(SC.coretable, rowSums(SC.coretable)>0))
      SC.panevol = c(SC.panevol, SC.pan)
      SC.strain = nrow(SC.coretable %>% filter((!!as.symbol(strain)) == 1 & rowSums(SC.coretable)==1))
      SC.strainevol = c(SC.strainevol, SC.strain)
    }
    #add species to first row
    SC_core = cbind(SC_core, SC.coreevol)
    SC_pan = cbind(SC_pan, SC.panevol)
    SC_strain = cbind(SC_strain, SC.strainevol)
  }
  colnames(SC_core) = colnames(SC.strainsdf)
  colnames(SC_pan) = colnames(SC.strainsdf)
  colnames(SC_strain) = colnames(SC.strainsdf)
    
  #Final merge of all results per species
  SC.results_core = rbind(SC.results_core, SC_core)
  SC.results_pan = rbind(SC.results_pan, SC_pan)
  SC.results_strain = rbind(SC.results_strain, SC_strain)

}

SC.results_core$Species = sort(clusters$SpeciesCluster)
SC.results_pan$Species = sort(clusters$SpeciesCluster)
SC.results_strain$Species = sort(clusters$SpeciesCluster)
```

Calculate the mean and SD using the 100 replicates.

```{r statistics}
#Create iteration index
cn=clusters %>% group_by(SpeciesCluster) %>% summarize(count=n())

iter=data.frame()
n=0

for (i in cn$SpeciesCluster) {
  n= n+1
  Iteration=as.numeric(seq(from=1, to=cn$count[n]))
  Species=rep(i, cn$count[n])
  res= cbind(Iteration,Species)
  iter=rbind(iter, res)
}

#Core genome
SC.core_mean = data.frame(Mean=rowMeans(SC.results_core[,1:101]))
SC.core_sd = data.frame(SD=rowSds(as.matrix(SC.results_core[,1:101])))
SC.core_stats = cbind(Strain=SC.results_core[102], Interation=iter$Iteration, SC.core_mean, SC.core_sd)
SC.core_stats$Species = as.numeric(gsub("species", "", SC.core_stats$Species))
#write.table(SC.core_stats, "SC.core_stats.txt", quote = F, sep = "\t")

#Pan genome
SC.pan_mean = data.frame(Mean=rowMeans(SC.results_pan[,1:101]))
SC.pan_sd = data.frame(SD=rowSds(as.matrix(SC.results_pan[,1:101])))
SC.pan_stats = cbind(Strain=SC.results_pan[102], Interation=iter$Iteration, SC.pan_mean, SC.pan_sd)
SC.pan_stats$Species = as.numeric(gsub("species", "", SC.pan_stats$Species))
#write.table(SC.pan_stats, "SC.pan_stats.txt", quote = F, sep = "\t")

#Strain-specific genome
SC.strain_mean = data.frame(Mean=rowMeans(SC.results_strain[,1:101]))
SC.strain_sd = data.frame(SD=rowSds(as.matrix(SC.results_strain[,1:101])))
SC.strain_stats = cbind(Strain=SC.results_strain[102], Interation=iter$Iteration, SC.strain_mean, SC.strain_sd)
SC.strain_stats$Species = as.numeric(gsub("species", "", SC.strain_stats$Species))
#write.table(SC.strain_stats, "SC.strain_stats.txt", quote = F, sep = "\t")
```


Plot

```{r}
SCcor = ggplot(SC.core_stats, aes(x=as.numeric(Interation), y = Mean, group=as.factor(Species), ymin=Mean-SD, ymax=Mean+SD)) +
  geom_line(aes(color=as.factor(Species)), size=0.5) +
  geom_ribbon(aes(fill=as.factor(Species)), alpha=0.07) +
  ggtitle("Core genome") + 
  scale_color_manual(values=c("#2f0f3e", "#000000", "#ea785a", "#000000", "#000000", "#000000", "#feedb0", "#d3d398", "#c5dbc1", "#6cc5a3", "#4ea0a1", "#427c99")) +
  scale_fill_manual(values=c("#2f0f3e", "#000000", "#ea785a", "#000000", "#000000", "#000000", "#feedb0", "#d3d398", "#c5dbc1", "#6cc5a3", "#4ea0a1", "#427c99")) +
  theme_classic() +
  theme(axis.line.x=element_line(), axis.line.y=element_line())+
  xlab("Genomes sampled") + ylab("Core genome size (OGs)") + scale_x_continuous(breaks=c(1,5,10,15,20,25,30,35,40,45,50,55,60,65,70), limits = c(1,70))
  
SCpan = ggplot(SC.pan_stats, aes(x=as.numeric(Interation), y = Mean, group=as.factor(Species), ymin=Mean-SD, ymax=Mean+SD)) +
  geom_line(aes(color=as.factor(Species)), size=0.5) +
  geom_ribbon(aes(fill=as.factor(Species)), alpha=0.07) +
  ggtitle("Pangenome") + 
  scale_color_manual(values=c("#2f0f3e", "#000000", "#ea785a", "#000000", "#000000", "#000000", "#feedb0", "#d3d398", "#c5dbc1", "#6cc5a3", "#4ea0a1", "#427c99")) +
  scale_fill_manual(values=c("#2f0f3e", "#000000", "#ea785a", "#000000", "#000000", "#000000", "#feedb0", "#d3d398", "#c5dbc1", "#6cc5a3", "#4ea0a1", "#427c99")) +
  theme_classic() +
  theme(axis.line.x=element_line(), axis.line.y=element_line())+
  xlab("Genomes sampled") + ylab("Pangenome size (OGs)") + scale_x_continuous(breaks=c(1,5,10,15,20,25,30,35,40,45,50,55,60,65,70), limits = c(1,70)) +
  scale_y_continuous(limits = c(5500, 10000))

SCstr = ggplot(SC.strain_stats, aes(x=as.numeric(Interation), y = Mean, group=as.factor(Species), ymin=Mean-SD, ymax=Mean+SD)) +
  geom_line(aes(color=as.factor(Species)), size=0.5) +
  geom_ribbon(aes(fill=as.factor(Species)), alpha=0.07) +
  ggtitle("Genome-specific genome") + 
  scale_color_manual(values=c("#2f0f3e", "#000000", "#ea785a", "#000000", "#000000", "#000000", "#feedb0", "#d3d398", "#c5dbc1", "#6cc5a3", "#4ea0a1", "#427c99")) +
  scale_fill_manual(values=c("#2f0f3e", "#000000", "#ea785a", "#000000", "#000000", "#000000", "#feedb0", "#d3d398", "#c5dbc1", "#6cc5a3", "#4ea0a1", "#427c99")) +
  theme_classic() +
  theme(axis.line.x=element_line(), axis.line.y=element_line())+
  xlab("Genomes sampled") + ylab("Genome-specific genome size (OGs)") + scale_x_continuous(breaks=c(1,5,10,15,20,25,30,35,40,45,50,55,60,65,70), limits = c(1,70))


grid.arrange(SCcor, SCpan, SCstr, ncol = 3)
```


#2 Functional Characterization

NMDS analysis of the Functional BRITE KEEG categories with the orthologous groups indentified by OrthoFinder.

##2.1 Load packages

```{r}
library(tidyverse)
library(rentrez)
library(stringr)
library(devtools)
#install_github("vqv/ggbiplot")
library(ggbiplot)
library(AMR)
library(dplyr)
#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("KEGGREST")
library(KEGGREST)
library(vegan)
library(factoextra)
library(pals)
```

##2.2 Load the Orthogroups.txt file

```{r}
orthogroups = read.csv("Orthogroups.txt", sep=":", header = F, col.names = c("OG", "Prot_list")) 
head(orthogroups)
```


##2.3 Get OG protein IDs

The OG groups contains orthologous sequences. These contain the protein acc. no. of the different genomes. We are going to retrieve the sequence of either P. protegens CHA0, because the genome is highly curated. If not available, the sequence of the first genome it can find will be retrieved.

```{r}
prots =  as.data.frame(matrix(ncol = 2, nrow = 0))

for (i in 1:as.numeric(length(orthogroups$Prot_list))) {
  spl = strsplit(orthogroups$Prot_list[i], "[[:space:]]")
  for (j in 1:as.numeric(length(spl[[1]]))) {
    if (grepl("0000\\|", spl[[1]][j])) { #this is CHA0
      prot = spl[[1]][j]
      f = c(orthogroups$OG[i], prot)
      prots = rbind(prots, f)
      break
    } else if (grepl("00..\\|", spl[[1]][j])) { #This will select the first protein 
      prot = spl[[1]][j]
      f = c(orthogroups$OG[i], prot)
      prots = rbind(prots, f)
      break
    } else if (grepl("\\|", spl[[1]][j])) { #This will select the first protein
      prot = spl[[1]][j]
      f = c(orthogroups$OG[i], prot)
      prots = rbind(prots, f)
      break
    }
  }
}

colnames(prots) = c("OG", "protein")
```

Separate the species code from the actual protein ID

```{r}
prots2 = as_tibble(prots) %>% separate(protein, sep = "\\|", c("species_id", "Protein"))
head(prots2)
```

##2.4 Retrieve protein sequences

```{r, DO NOT EXECUTE IF ALREADY DONE!!!}
#SET API KEY
set_entrez_key("a0ef24798beb1779c91b955913b2c605ad08")
Sys.getenv("ENTREZ_KEY")

prot_seq = as.data.frame(matrix(ncol=2, nrow=0))
 for (i in 1:as.numeric(length(prots2$Protein))) {
   while(TRUE){
     ids2 = try(entrez_fetch(db="protein", id = prots2$Protein[i], rettype="fasta"))
     if(!is(ids2, 'try-error')) break
   }
   sequence = str_remove(ids2, "(.*)\n")
   sequence = str_remove_all(sequence, "\n")
   
   perc = i * 100 / as.numeric(length(prots2$Protein)) #calculate the percentage done.
   print(paste("Iteration no.", i, " ", prots2$Protein[i], "->", perc, "%")) #print the percentage done
   seq = c(prots2$Protein[i], sequence)
   prot_seq = rbind(prot_seq, seq)
}
 
colnames(prot_seq) = c("ID", "Sequence")
#write.table(prot_seq, "prot_seqFinal.txt", sep = "\t", row.names = F, quote = F)
#prot_seq = read.csv("prot_seqFinal.txt", sep = "\t", header = T, quote = "", na.strings = "")
```

As a sanity check, check that the number of observations in this prot_seq is the same as in orthogroups!

##2.5 Generate multifasta for functional annotation.

```{r}
multifasta = 0
count = 0
for (i in 1:length(prot_seq$Sequence)) {
  id = prot_seq$ID[i]
  seq = prot_seq$Sequence[i]
  count = count + 1
  fasta = paste(">", id, "\n", seq, sep = "") #this will transform into multifasta.
  multifasta = c(multifasta, fasta)
}
 
print(paste("Multifasta done. ", i, " sequences written. Now generating the file...", sep = ""))


multifasta = as.data.frame(multifasta[-1])
#write.table(multifasta, "multifasta.fa", row.names = F, col.names = F, quote = F)
```


##2.6 Annotate with eggNOG

We are going to use eggNOG mapper (https://github.com/eggnogdb/eggnog-mapper/wiki/eggNOG-mapper-v2.1.5#Installation) to annotate all these proteins.


```{bash}
#download main annotation databases:
download_eggnog_data.py --data_dir /mnt/f/Databases/eggnog/
#run it
emapper.py -i multifasta.fa -o eggnog_annotation -m diamond --cpu 44 --data_dir /mnt/f/Databases/eggnog/ -d bact
```

##2.7 Retrieve BRITE classification

```{r}
eggnog = read.csv("eggnog_annotation.emapper.annotations", 
                  header = T, 
                  sep = "\t", 
                  quote = "", 
                  skip = 4, 
                  na.strings = "")
eggnog_filt = eggnog %>% select(X.query, BRITE)
colnames(eggnog_filt) = c("query_name", "BRITE")
head(eggnog_filt)
```

And process

```{r}
brite = as.data.frame(matrix(ncol = 4, nrow = 0))
for (i in 1:length(head(eggnog_filt$BRITE,-3))) {
  id = eggnog_filt$query_name[i]
  print(paste("Seq", i, id, "-->", i * 100/length(eggnog_filt$BRITE) ))
  if (eggnog_filt$BRITE[i] == "-") {
    lastbri = "-"
    bri_annot = "-"
    bri = c(id, lastbri, bri_annot)
    brite = rbind(brite, bri)
    next
  }
  lastbri = sapply(strsplit(eggnog_filt$BRITE[i], ","), tail, 1)
  while(T){
     bri_annot = try(keggFind("brite", lastbri))
     if(!is(bri_annot, "try-error")) break
  }
  bri = c(id, lastbri, bri_annot)
  brite = rbind(brite, bri)
}

colnames(brite) = c("ID", "lastBRITE", "BRITE_ANNOTATION")
#write.table(brite, "brite.txt", quote = F, sep = "\t")
```

Gather info

Add the OG, speicies_ID, etc to the brite table

```{r}
brite_table = left_join(prots2, brite, by = c("Protein" = "ID"))
brite_table = brite_table %>% select(OG, species_id, Protein, lastBRITE, BRITE_ANNOTATION)
head(brite_table)
```

Get the table with the paralos and the total number of proteins per genome

```{r}
og_sing_parag = read.csv("Orthologs_and_singletons.txt",
                         header = T, sep = "\t", quote = "")
total_prots = as.data.frame(colSums(og_sing_parag[,2:length(og_sing_parag)]))
colnames(total_prots) = "TotalProts"
```

What we need is a table of the form:

BRITE genome1 genome2 genome3 ... 
T        1      5       1 ... 
V        3      3       3 ... 
H        1      2       2 ... 
...

We can match the BRITE with the OG thanks to the og_sing_parag table... and then group by BRITE.

```{r}
og_brite = left_join(brite_table, og_sing_parag, by = "OG")
og_brite = og_brite %>% select(-OG, -lastBRITE, -species_id, -Protein) #remove these columns

#Group by brite
table_brite_final = og_brite %>% group_by(BRITE_ANNOTATION) %>% summarise_all(list(sum))

#Remove NAs and proteins not annotated
brite_filt = table_brite_final %>% subset(BRITE_ANNOTATION != "NA" & !grepl("^WP", BRITE_ANNOTATION) & !grepl("^KEGG", BRITE_ANNOTATION) & !grepl("-", BRITE_ANNOTATION))

#Add BRITE as row names instead of columns
bri = as.data.frame(brite_filt)
rownames(bri) = bri$BRITE_ANNOTATION
bri = bri[,-1]
```

We need the grouping variable for the genomes, in species groups. Load this table containing the SpeciesCluster and the GenomeName

```{r}
names = read.table("Names.txt", header = T, sep = "\t")
head(names)
```

##2.8 Calculate relative abundances per genome

Aout of the total number of assigned BRITEs per genome, to avoid genome-size inference.

```{r}
#this takes into account the genome size by divided the proteins assiged to brite cagetory by the total proteins in the geome
brite2perc.t = apply(bri, 1, FUN = "/", t(total_prots/100))
rownames(brite2perc.t) = names$GenomeName
```

##2.9 Filter genomes by contigs

Calculate number of contigs per genome and filter

```{r}
#exceute in bash
#grep -c '>' *.fna > conting_count.txt

contigs = read.table("contig_count.txt", sep = "\t", col.names = c("GenomeName", "Contigs"))

#Filter genomes having more or 85 contigs
contigs_filt = contigs %>% subset(Contigs <= 85)
head(contigs_filt)
```

From 115 genomes we now have 80

Filter the "brite2perc.t" table

```{r}
briteperc.t_filt = as.data.frame(brite2perc.t) %>% rownames_to_column("GenomeName") %>% inner_join(contigs_filt, by = "GenomeName")
briteperc.t_filt = column_to_rownames(briteperc.t_filt, "GenomeName") %>% select(-Contigs)
head(briteperc.t_filt)
```

##2.10 NMDS

```{r}
bridist_filt = vegdist(briteperc.t_filt, method = "bray")

set.seed(1) 
bri_filtnmds = metaMDS(bridist_filt, distance = "bray", k = 3, trymax = 100)
bri_filtnmds$stress

goodness(bri_filtnmds)
stressplot(bri_filtnmds)
plot(bri_filtnmds, "sites")
orditorp(bri_filtnmds, "sites")
```

Obtain the scores of the NMDS to further plot this points with ggplot.

```{r}
sites_filt.scrs = as.data.frame(scores(bri_filtnmds, display="sites"))

#add names
sites_filt.scrs = rownames_to_column(sites_filt.scrs, "GenomeName")
sites_filt = left_join(sites_filt.scrs, names)
head(sites_filt)
```

Plot

```{r}
nmds_plot = ggplot(sites_filt, aes(x=NMDS1, y=NMDS2)) +
  geom_point(aes(colour = as.factor(SpeciesCluster)), size=4, shape=21, alpha=0.75) +
  geom_point(aes(colour = as.factor(SpeciesCluster)), size=4, alpha=0.5) +
  scale_y_continuous(limits = c(-0.03,0.025), breaks = c(-0.03,-0.020,-0.01,0.00, 0.01, 0.02))+
  annotate("text",x=0.02, y=-0.03, label = paste0("Stress: ", format(round(bri_filtnmds$stress, 5), nsmall = 5))) +
  geom_hline(yintercept = 0.0, color= "grey", linetype="dashed") +
  geom_vline(xintercept = 0.0, color= "grey", linetype="dashed") +
  #ggtitle("With CHA0") +
  #geom_text(aes(label=sites_filt$GenomeName, hjust=0.5, vjust=-1.2), size = 2, alpha = 0.8) +
  scale_color_manual(values=c(ocean.matter(4), rev(ocean.deep(8))))+
  theme_classic() 

nmds_plot
```

Retrieve variables

```{r}
nmds.test.fit_filt = envfit(bri_filtnmds, briteperc.t_filt, permutations = 999)

plot(bri_filtnmds)
plot(nmds.test.fit_filt)

test.scrs_filt = as.data.frame(scores(nmds.test.fit_filt, display = "vectors")) * ordiArrowMul(nmds.test.fit_filt)
test.scrs_filt = cbind(test.scrs_filt, Test= rownames(test.scrs_filt)) #add test names
test.scrs_filt = cbind(test.scrs_filt, pval = nmds.test.fit_filt$vectors$pvals) #add p-values to dataframe
nmds.test.fit_filt
test.scrs_filt

```

Remove not significant variables

```{r}
test.scrs_filt2 = test.scrs_filt %>% filter(pval <= 0.05)
test.scrs_filt2
#save to retrieve information later
#write.table(test.scrs_filt2, "scr_filt2.txt", quote = F, sep = "\t", row.names = F)
```

Plot

```{r}
library(ggpubr)

#Assign each variable to a number
test.scrs_filt2$No= seq(1,dim(test.scrs_filt2)[1])

ggplot(sites_filt, aes(x=NMDS1, y=NMDS2)) +
  geom_segment(data=test.scrs_filt2, 
               aes(x=0, xend = NMDS1, y = 0, yend = NMDS2), 
               lwd = 0.25, alpha= 0.75, color="cornsilk4", 
               arrow = arrow(length = unit(0.005, "npc"))) +
  scale_y_continuous(limits = c(-0.025,0.025), breaks = c(-0.03,-0.020,-0.01,0.00, 0.01, 0.02))+
  scale_x_continuous(limits = c(-0.025, 0.06)) +
  geom_hline(yintercept = 0.0, color= "grey", linetype="dashed") +
  geom_vline(xintercept = 0.0, color= "grey", linetype="dashed") +
  scale_color_manual(values=c(ocean.matter(4), rev(ocean.deep(8))))+
  ggrepel::geom_text_repel(data=test.scrs_filt2, aes(x=NMDS1, y=NMDS2, label = No), max.overlaps = Inf,
                           cex = 3, direction = "both", segment.size = 0.25, box.padding = 0.01,
                           segment.alpha=0.5, segment.color="grey") +
  coord_flip()+
  #theme_transparent()+
  labs(color="P. protegens\nspecies cluster no.")


ggplot(sites_filt, aes(x=NMDS1, y=NMDS2)) +
  geom_point(aes(colour = as.factor(SpeciesCluster)), size=4, shape=21, alpha=0.75) +
  geom_segment(data=test.scrs_filt2, 
               aes(x=0, xend = NMDS1, y = 0, yend = NMDS2), 
               lwd = 0.25, alpha= 0.75, color="cornsilk4", 
               arrow = arrow(length = unit(0.005, "npc"))) +
  #geom_text(data=test.scrs_filt2, aes(label=Test, hjust=0.5, vjust=-1.2), size = 2, alpha = 1, color="cornsilk4") +
  geom_point(aes(colour = as.factor(SpeciesCluster)), size=4, alpha=0.5) +
  scale_y_continuous(limits = c(-0.025,0.025), breaks = c(-0.03,-0.020,-0.01,0.00, 0.01, 0.02))+
  scale_x_continuous(limits = c(-0.025, 0.06)) +
  annotate("text", size =3, y=-0.02, x=-0.024, label = paste0("Stress: ", format(round(bri_filtnmds$stress, 5), nsmall = 5))) +
  geom_hline(yintercept = 0.0, color= "grey", linetype="dashed") +
  geom_vline(xintercept = 0.0, color= "grey", linetype="dashed") +
  #geom_text(data=test.scrs_filt2, aes(label=No, hjust=0.5, vjust=-1.2), size = 2, alpha = 0.8) +
  scale_color_manual(values=c(ocean.matter(4), rev(ocean.deep(8))))+
  ggrepel::geom_text_repel(data=test.scrs_filt2, aes(x=NMDS1, y=NMDS2, label = No), max.overlaps = Inf,
                           cex = 3, direction = "both", segment.size = 0.25, box.padding = 0.01,
                           segment.alpha=0.5, segment.color="grey") +
  coord_flip()+
  theme_classic()+
  labs(color="P. protegens\nspecies cluster no.")

```


##2.11 Significance (MANOVA)


```{r}
com = as.matrix(bridist_filt)
```

```{r}
ano = anosim(com, sites_filt$SpeciesCluster, distance = "bray", permutations = 9999)
ano
```

The result tell us that there are indeed differences between the groups. Now we need to know which groups, using a Post-hoc pairwise comparisons. We are going to use the pairwise.adonis

```{r}
#if (!require(pairwiseAdonis)) install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")
library(pairwiseAdonis)

set.seed(123)
pairwise_adonis_res = pairwise.adonis(x = com,
                    factors= as.numeric(sites_filt$SpeciesCluster), sim.method = 'bray',
                    p.adjust.m='fdr', perm = 9999)

pairwise_adonis_res

#write.table(pairwise_adonis_res, "pairwise_adonis_test.txt", quote = F, sep = "\t", row.names = F)
```


#3. Hierarchical clustering of intergenomic distances

##3.1 Load packages

```{r}
library(mclust)
library(fpc)
library(reshape2)
library(ggpubr)
library(ggplot2)
```

## 3.2 Load inputs

```{r}
dist_prot = read.table("Species_protegens/Distance_matrix_Protegens.txt", header = F, sep = "\t")
taxa_names = scan("Species_protegens/taxa_names_protegens.txt", character())
partition = scan("Species_protegens/protegens_partition.txt")
head(dist_prot)
```

## 3.3. Clustering at Species level:

Calculation of single-, average-, and complete-linkeage (F) across a range of distance threshold (T) from 0 to the maximum value of the distance matrix. Calculate different statistic indexes based on the reference partition, including the Adapted Rand Index (ARI), at T and F.

```{r warning=FALSE}
matr = as.dist(dist_prot) #this is the matrix as a dist object
step_size = 0.0005
min_thres = 0
max_thresh = max(matr)

#Create the hierarchical clustering based on single, average and complete linkage
clust_single = hclust(matr, method = 'single')
clust_averag = hclust(matr, method = 'average')
clust_complet = hclust(matr, method = 'complete')

#Create empty vectors to store the results
results = data.frame(matrix(ncol = 10, nrow = 0))
results_names = c("Threshold_T", "Clust_F0.0", "ARI_F0.0", "Dunn2_F0.0",
                  "Clust_F0.5", "ARI_F0.5", "Dunn2_F0.5",
                  "Clust_F1", "ARI_F1", "Dunn2_F1")

for (i in seq(min_thres, max_thresh, step_size)) {
  cut_single = cutree(clust_single, h=i)
  clusters_no_single = max(cut_single)
  ARI_single = adjustedRandIndex(partition, cut_single)
  if (clusters_no_single > 1) {
  single_all_tests = cluster.stats(matr, cut_single, partition)
  dunn2_single = single_all_tests$dunn2
  } else {dunn2_single = "-"}
  
  cut_average = cutree(clust_averag, h=i)
  clusters_no_average = max(cut_average)
  ARI_average = adjustedRandIndex(partition, cut_average)
  if (clusters_no_average > 1) {
  average_all_tests = cluster.stats(matr, cut_average, partition)
  dunn2_average = average_all_tests$dunn2
  } else {dunn2_average = "-"}
  
  cut_complete = cutree(clust_complet, h=i)
  clusters_no_complete = max(cut_complete)
  ARI_complete = adjustedRandIndex(partition, cut_complete)
  if (clusters_no_complete > 1) {
  complete_all_tests = cluster.stats(matr, cut_complete, partition)
  dunn2_complete = complete_all_tests$dunn2
  } else {dunn2_complete = "-"}
  
  threshold = i
  results_all = cbind(threshold, as.numeric(clusters_no_single), as.numeric(ARI_single), as.numeric(dunn2_single), 
                      as.numeric(clusters_no_average), as.numeric(ARI_average), as.numeric(dunn2_average), 
                      as.numeric(clusters_no_complete), as.numeric(ARI_complete), as.numeric(dunn2_complete))
  results = rbind(results, results_all)
}

colnames(results) = results_names

#write.table(results, "Clustering_results_species.txt", quote = F, sep = "\t", row.names = F)
```

Plot 

```{r}
plot_sp_clust = ggplot(results, aes(x=Threshold_T)) +
  geom_line(aes(y=Clust_F0.0, color = "F = 0.0 (single)")) +
  geom_line(aes(y=Clust_F0.5, color = "F = 0.5 (average)")) +
  geom_line(aes(y=Clust_F1, color = "F = 1 (complete)")) +
  scale_color_manual(values = c("khaki3", "green", "steelblue", "gray40", "darkgreen", "black")) +
  
  geom_line(aes(y=ARI_F0.0*100, colour = "ARI w/F = 0.0 (single")) +
  geom_line(aes(y=ARI_F0.5*100, colour = "ARI w/F = 0.5 (average)")) +
  geom_line(aes(y=ARI_F1*100, colour = "ARI w/F = 1 (complete)"))+
  

  scale_y_continuous(sec.axis = sec_axis( trans=~./100, name="Adapted Rand Index"), name = "Clusters no.")+
  geom_vline(xintercept = 0.036, color = "black", alpha = 0.3, linetype = "longdash") +
  geom_hline(yintercept = 12, color = "red", alpha = 0.3)+
  labs(color = "Linkage F") +
  xlab(" Distance Threshold T") + 
  ggtitle("Species clustering", subtitle = "Distance = 0.036, dDDH% = 70%, clusters = 12")+
  theme_gray()
  
plot_sp_clust
```


#4. Biolog analysis 

Heatmap based on the Biolog results.

##4.1 Load packages

```{r}
library(ggplot2)
library(tidyverse)
library(ComplexHeatmap)
library(scico)
```

##4.2Load Data

As input data we are gonna use the Biolog data in which the first two columns are the category and the compounds, and the rest of the columns are each strain tested. Each row is a individual compound and the OD obtained.

```{r}
data = read.table("biolog_data.txt", header = T, sep = "\t")
info = data %>% select(Category, Compounds)
head(info)

data = column_to_rownames(data, "Compounds") %>% select(-Category)
head(data)
```

Create the df for insecticidal or not strain

```{r}
insec = data.frame(Strain = colnames(data), "Insecticidal" = c("yes", "yes", "yes", "yes", "no", "no", "yes", "yes", "no", "no", "yes", "yes"))
insec = column_to_rownames(insec, "Strain")
insec
```


##4.3 Filtering

Remove  rows in which the sum of all compounds is 0 (meaning no differences between the strains) and also applying a minimum number to consider growth, in this case 0.1

```{r}
data.filt = data
threshold = 0.1
data.filt[data.filt < threshold] = 0

data.filt = data.filt %>% filter(rowSums(.) > 0)
head(data.filt)
```

##4.4 Heatmap

```{r}
insect = HeatmapAnnotation(df = insec, col = list(Insecticidal = c("yes" = "firebrick3", "no" = "grey80")), border = TRUE)

Heatmap(data.filt, clustering_distance_columns = "euclidean", clustering_distance_rows = "euclidean",
        column_dend_reorder = F, row_dend_reorder = T,
        column_names_side = "top", column_dend_side = "bottom", column_names_rot = 40,
        row_names_gp = gpar(fontsize=8),
        #row_title = "Compounds",
        #row_km = 2, row_km_repeats = 10)
        row_split = info$Category, row_gap = unit(2, "mm"), border = TRUE,
        top_annotation = insect, 
        col = scico(palette = "lajolla", n=50))

```

