---
title: "Multiple extracellular polymeric substances pathways expressed by Accumulibacter and the flanking community during aerobic granule formation and after influent modification"
author: "LBE - Laëtitia Cardona"
date: "12.12.2024"
output:
  html_document:
    number_sections: yes
    toc: yes
    toc_depth: 4
  pdf_document:
    number_sections: yes
    toc: yes
    toc_depth: 4
  word_document:
    number_sections: yes
    toc: yes
    toc_depth: 4
header-includes:
  \usepackage{dcolumn}
---

# General information

## Authors & affiliation
Laëtitia Cardona, Jaspreet Singh Saini, Pilar Natalia Rodilla Ramírez, Aline Adler, Christof Holliger

*Laboratory for Environmental Biotechnology, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland

## Corresponding author
laetitia.cardona@epfl.ch

## Complementary information

*preprint* https://doi.org/10.1101/2024.09.04.611157 
*zenodo* https://zenodo.org/records/13693480
*16S rRNA gene raw data* PRJNA1125294
*raw and de novo assembly metatranscriptomic* PRJNA1144857

# Abastract

Aerobic granular sludge is a biological wastewater treatment process in which a microbial community forms a granular biofilm. The role of Candidatus Accumulibacter in the production of a biofilm matrix composed of extracellular polymeric substances was studied in a sequencing batch reactor enriched with polyphosphate-accumulating organisms. The metabolisms of the microbial populations were investigated using de novo metatranscriptomics analysis. Finally, the effect of decreasing the influent phosphate concentration was investigated. 
A few weeks after the reactor start-up, the microbial community was dominated by Accumulibacter. Up to nine species were active in parallel. However, the most active species differed according to sampling time. Reducing the phosphate concentration led to a dominance of the glycogen-accumulating organism Propionivibrio, with some Accumulibacter species still abundant. De novo metatranscriptomics analysis indicated a high diversity of potential extracellular substances produced mainly by Azonexus, Accumulibacter, Candidatus Contendobacter, and Propionivibrio. Moreover, the results suggest that Azonexus, Contendobacter and Propionivibrio recycle the neuraminic acid produced by Accumulibacter. Changes in the microbial community did not cause the granules to disintegrate, indicating that a Propionivibrio-dominated community can maintain stable granules.

# Keywords

metatranscriptomics, sequencing batch reactor, Propionivibrio, 16S rRNA gene amplicon sequencing, Enhanced biological phosphate removal

# Figures

```{r global_options, include=FALSE}
# Clean the Global Environment
rm(list=ls(all=TRUE))

# Data
load("./Data.RData")

# description of the files
## ASV.count = count table from the DADA2 pipeline
## ASV.taxa = taxonomic information for each ASV
## ASV.data = metadata of the samples from 16S rRNA gene amplicon sequencing

## FC.count = gene count from de novo metatranscriptomic analysis
## FC.annot = annotation of the genes from DRAM and eggNOG combined
## FC.annot.no0.GoI = 
## FC.data = metadata of the samples from the metatransriptomic analysis

## FC.genes.list = list of genes of interest

## mOTUs.count = count table from mOTUs analysis
## mOTUs.data = metadata of the samples from the metatranscriptomic analysis
## mOTUs.taxa = taxonomic information from mOTUs analysis

## PERF.melt = results of the chemical analysis in long table format

## download the packages
### from bioconductor
library(phyloseq)
library(edgeR)
library(DESeq2)

### from CRAN
library(RColorBrewer)
library(reshape2)
library(ggpubr) 
library(plyr)
library(tidyverse)

# remove legend from graph
rm_legend <- function(p){p + theme(legend.position = "none")}

```

## Figure 1. Nutrient removal and efficiency

For this part, the file used are the AppendixB_16S-rRNA-gene-amplicon, where:
*16S_metadata* contains the performance results
*16S_count taxa* group the count and taxonomy results

```{r 1_PERF, include=TRUE, fig.height=5, fig.width=6}
# phosphate plot
plot.P <- 
  PERF.melt %>% filter(Analysis == "P-PO[4] (mg/L)") %>%
  ggplot(aes(x = as.numeric(Day), y = value, color = Phase, shape = Phase)) +
  geom_point(size = 3, na.rm = TRUE) +
  geom_line(na.rm = FALSE) +
  scale_shape_manual("Phase", values = c(15, 17, 6), breaks = c("Influent", "Anaerobic", "Aerobic")) +
  scale_colour_manual("Phase", values = c("#888888", "#999933", "#DDCC77"), 
                      breaks = c("Influent", "Anaerobic", "Aerobic")) +
  scale_x_continuous(name= "Time (days)", breaks = seq(0, 182, 20)) +
  geom_vline(xintercept = 182, color = "white") +
  geom_vline(xintercept = 103, color = "black", lty = 2) +
  scale_y_continuous(breaks = seq(0, 300, 50)) +
  labs(y = bquote(P-PO[4] (mg/L))) +
  theme_bw() +
  theme(axis.text.y = element_text(size = 10, face = "bold"),
        axis.text.x = element_text(size = 10, face = "bold"),
        axis.title.y = element_text(size = 12),
        axis.line = element_line(color = 'black'),
        plot.margin = margin(0.1, 0.1, 0.1, 0.1, "cm"))

# COD plot
plot.CODrem <- 
  PERF.melt %>% filter(Analysis == "C[removal]") %>%
  ggplot(aes(x = as.numeric(Day), y = value)) +
  geom_point(size = 3, na.rm = TRUE, color = "#999933") +
  geom_line(na.rm = FALSE, color = "#999933") +
  scale_x_continuous(name= "", breaks = seq(0, 182, 20))+
  geom_vline(xintercept = 182, color = "white") +
  geom_vline(xintercept = 103, color = "black", lty = 2) +
  labs(y = "C removal (%)") +
  theme_bw() +
  theme(axis.text.y = element_text(size = 10, face = "bold"),
        axis.text.x = element_text(size = 10, face = "bold"),
        axis.title.y = element_text(size = 12),
        axis.line = element_line(color = 'black'),
        plot.margin = margin(0.1, 0.1, 0.1, 0.1, "cm"))

ggarrange(rm_legend(plot.CODrem), rm_legend(plot.P), 
          ncol = 1, nrow = 2,
          align = "v",
          heights = c(0.2,0.2))
```

## Figure 2. Relative abundance of the most abundant microorganisms

```{r 2_16S, include=TRUE, fig.height=4, fig.width=10}
# Create phyloseq object
ASV.physeq <- phyloseq(otu_table(ASV.count, taxa_are_rows = FALSE),
                        sample_data(ASV.data),
                        tax_table(as.matrix(ASV.taxa)))

# Some samples seemed to be outliers, for the purpose of the figure we removed them
ASV.physeq.sub <- subset_samples(ASV.physeq, !Names %in% c("RAen2_26", "RAen2_111", "RAen2_131"))

# agglomerate at genus level
ASV.physeq.sub.glom <- tax_glom(ASV.physeq.sub, taxrank = "Genus")

# To make the figure readable by removing the low abundant ASVs we filtered out those lower than 3%
minper = 3
filter.param = genefilter_sample(ASV.physeq.sub.glom, filterfun_sample(function(x) x/sum(x)*100 >= minper))
ASV.sub.glom.psmelt <- ASV.physeq.sub.glom %>%
  transform_sample_counts(function(x) x/sum(x)) %>%
  prune_taxa(filter.param, .) %>%
  psmelt()

# here we ordered manually the genera to group them by class
ASV.sub.glom.psmelt$Genus <- factor(ASV.sub.glom.psmelt$Genus, 
                                       levels = c("Microtrichaceae",
                                                  "Ca_Phosphoribacter","Propionicimonas", 
                                                  "UKL13-1", "Micavibrionales", "Rhodobacter", "Defluviimonas", "Novosphingobium", 
                                                  "env.OPS_17", "Flavobacterium", "Flavobacteriaceae", "Chryseobacterium", "Saprospiraceae", 
                                                  "Herpetosiphon",
                                                  "Run-SP154", "Ca_Contendobacter", "Ca_Competibacter", "Ca_Accumulibacter", "Propionivibrio", "Azonexus", "Thiothrix", "Pseudoxanthomonas",
                                                  "Nitrospira", 
                                                  "Salinispira"))

# we modify the naming
mic.labels <- c("Microtrichaceae", 
                "Ca. Phosphoribacter","Propionicimonas", 
                "UKL13-1", "Micavibrionales", "Rhodobacter", "Defluviimonas", "Novosphingobium", 
                "env.OPS_17", "Flavobacterium", "Flavobacteriaceae", "Chryseobacterium", "Saprospiraceae",
                "Herpetosiphon",
                "Run-SP154", "Ca. Contendobacter", "Ca. Competibacter", "Ca. Accumulibacter", "Propionivibrio", "Azonexus", "Thiothrix", "Pseudoxanthomonas",
                "Nitrospira", "Salinispira")

# define specific colors
col.BAC = c("mediumpurple4",
            "khaki", "khaki3", 
            "#084594", "#4292C6", "#9ECAE1", "#2171B5", "#DEEBF7", 
            "#005A32", "#2c743f", "#41AB5D", "#A1D99B", "#E5F5E0", 
            "#737373", 
            "red4", "tomato4", "tan4", "#C994C7","#FFD92F", "tan2", "wheat", "snow2", 
            "#E7298A", 
            "azure4")

# plot
ASV.sub.glom.psmelt %>%
  filter(!Day %in% c("128", "127")) %>%
  ggplot(., aes(x = Day, y = Abundance, fill = Genus)) +
  geom_area(alpha = 0.85) +
  labs(y="Relative abundance (%)",x="Time (days)", fill = "") +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(labels = mic.labels, values=col.BAC, limits = force) +
  scale_x_continuous(breaks = seq(0, 182, 20))+
  geom_vline(xintercept = 182, color = "white") +
  geom_vline(xintercept = 103, color = "black", lty = 2) +
  theme_bw()+
  theme(legend.position = "right",
        axis.text.y = element_text(size = 10, face = "bold"),
        axis.text.x = element_text(size = 10, face = "bold"),
        axis.title.y = element_text(size = 12),
        axis.line = element_line(color = 'black'),
        plot.margin = margin(0.1, 0.1, 0.1, 0.1, "cm")) 
```

## Figure 3. Active microbial population at different stages of granulation

For this part of the code, we use files form AppendixC_BioInfoSummary, where:
*geneCounts_Raw* scorrespond in the code to FC.count
*geneCount_filtered-cpm* corresponds was obtained after 0 and low cont filtering followed by normalisation (in the script = edgeR.list.lfc.norm.cpm)

The annotation file is the AppendixD_geneAnnotation.
*annotation of all genes* corresponds to FC.annot

### Filtration and normalisaiton of the data

```{r}
#################### Filtering structural 0 ####################
# remove all 0
lowcount <- rowSums(FC.count) > 0
table(lowcount)
FC.count.no0 <- FC.count[lowcount,]

# get annotation
FC.annot.no0 <- FC.annot[FC.annot$X %in% rownames(FC.count.no0), ]

# put samples in order
FC.count.no0 <- FC.count.no0[ , FC.data$Sample]

#################### Filtering low count ####################

# determine on which data we work
counts = FC.count.no0
group = FC.data$Name

# put in DGE list for edgeR
edgeR.list <- DGEList(counts = counts, group = group)

# gene filtering: remove low count as low evidence of a real DGE
keep <- filterByExpr(y = edgeR.list, group = group, min.count = 15)
edgeR.list.lfc <- edgeR.list[keep, , keep.lib.sizes=FALSE]
# dim(edgeR.list.lfc)
# colSums(edgeR.list.lfc$counts)

melt(data.frame(names = colnames(FC.count.no0), 
           before_lfc = colSums(FC.count.no0),
           after_lfc = colSums(edgeR.list.lfc$counts))) %>%
  ggplot(aes(x = names, y = value, fill = variable)) +
  geom_bar(stat="identity", color="black", position=position_dodge()) + 
  coord_flip() +
  theme_minimal()

# check filtering
## average libray size in millions => min log-cpm value for each sample is then log2(2/L)  = -3.5072
L <- mean(edgeR.list$samples$lib.size) * 1e-6
M <- median(edgeR.list$samples$lib.size) * 1e-6
# c(L, M)
lcpm.cutoff <- log2(15/M + 2/L)

nsamples <- ncol(edgeR.list.lfc)
col <- brewer.pal(nsamples, "Paired")
lcpm <- cpm(edgeR.list, log=TRUE)
plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.4), las=2, main="", xlab="")
title(main="A. Raw data", xlab="Log-cpm")
abline(v=lcpm.cutoff, lty=3)
for (i in 2:nsamples){
  den <- density(lcpm[,i])
  lines(den$x, den$y, col=col[i], lwd=2)
}
lcpm <- cpm(edgeR.list.lfc, log=TRUE)
plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.4), las=2, main="", xlab="")
title(main="B. Filtered data", xlab="Log-cpm")
abline(v=lcpm.cutoff, lty=3)
for (i in 2:nsamples){
  den <- density(lcpm[,i])
  lines(den$x, den$y, col=col[i], lwd=2)
}

# filter the annotation file
FC.annot.no0.lfc <- FC.annot.no0[FC.annot.no0$X %in% rownames(edgeR.list.lfc$counts), ]

#################### Normalisation ####################

# normalisation to account for library size
edgeR.list.lfc.norm <- normLibSizes(edgeR.list.lfc, method = "TMM")
edgeR.list.lfc.norm.cpm <- cpm(edgeR.list.lfc.norm)

```

### fig 3.A) Relative proportion of gene expression per genus from the de novo metatranscriptomics analysis

```{r 3_activeCom, fig.height=3, fig.width=5}
#################### looking at genus level based on Kaiju results ####################
FC.Expr <- melt(t(edgeR.list.lfc.norm.cpm)) %>% left_join(FC.annot.no0.lfc, join_by(Var2 == X)) %>%
  left_join(FC.data, join_by(Var1 == Sample)) 

FC.Expr.taxo.sum <- ddply(FC.Expr, c("Genus", "Phase", "Day", "Cycle"), summarise,
                          sum = sum(value), # sum data to taxo level
                          ngenes = length(value)) %>% # determine number of genes per taxo 
  filter(!is.na(.$Genus))

FC.Expr.taxo.sum$Phase <- factor(FC.Expr.taxo.sum$Phase, levels = c("Feeding", "Aerobic"))
FC.Expr.taxo.sum$Day <- factor(FC.Expr.taxo.sum$Day, levels = c("13", "26", "103", "182"))

# microbial abundance
FC.Expr.taxo.med <- FC.Expr.taxo.sum %>%
  group_by(Day, Genus) %>%
  summarise(median = median(sum)) %>%
  group_by(Day) %>%
  mutate(totalSum = sum(median)) %>%
  ungroup() %>%
  filter(Genus %in% unique(subset(., .$median >= 3000)$Genus)) %>%
  mutate(perc = median/totalSum,
         Genus = gsub("NA", "Unknown Genus", Genus))

# get color palette
col.genus = c("brown3","indianred1","darkseagreen", "#BDBDBD", "#C994C7","#1B9E77", "#4292C6", 
                "#BF812D", "#525252", "burlywood3", "darkolivegreen", "#FFD92F", "lightskyblue3", "beige")

FC.Expr.taxo.med %>% 
  ggplot(aes(x = Day, y = perc, fill = Genus)) +
  geom_bar(stat = "identity") +
    scale_fill_manual(values = col.genus,
                    limits = force, labels = c("Alkalilimnicola",
                              "Aurantibacter",
                              "Azoarcus",
                              "Azonexus",
                              "Accumulibacter",
                              "Contendobacter",
                              "Cloacibacterium",
                              "Flavobacterium",
                              "Kaistella",
                              "Nitrosomonas",
                              "Pedobacter",
                              "Propionivibrio",
                              "Thiothrix",
                              "Unknown Genus")) +
  labs(y="Relative expression",x ="Days") +
  theme_minimal() +
  theme(legend.key.size = unit(0.5, "cm"), 
        legend.title = element_blank(),
        legend.text =  element_text(size = 10, face = "bold"),
        axis.text.y = element_text(size = 11, face = "bold"),
        axis.text.x = element_text(size = 11, face = "bold"),
        axis.title.y = element_text(size = 13),
        axis.line = element_line(color = 'black'))
```

### fig 3.B) Shannon diversity index calculated using the level of expression of each genus from the de novo metatranscriptomics analysis

```{r 4_Shannon, fig.height=3, fig.width=5}
#Function to calculate shannon diversity
shannon <- function(x){
  
  rabund <- x[x>0]/sum(x)
  -sum(rabund * log (rabund))
}

FC.exp.div <- FC.Expr.taxo.sum %>%
  group_by(Day, Cycle, Phase) %>%
  summarise(Shannon = shannon(sum)) 

FC.exp.div %>% 
  ggplot(aes(x = Day, y = Shannon, fill = Day)) +
  geom_boxplot() + geom_point() +
  scale_fill_manual(values = c("#DDCC77", "#999933", "#117733", "#661100")) +
  labs(y = "Shannon index", x = "Days") +
  scale_y_continuous(limits = c(0, NA)) +
  #scale_y_break(c(0.7, 0.8), space = 0.5, ) +
  theme_minimal() +
  theme(legend.key.size = unit(1, "cm"), 
        legend.title = element_blank(),
        legend.text =  element_text(size = 11, face = "bold"),
        axis.text.y = element_text(size = 11, face = "bold"),
        axis.text.x = element_text(size = 11, face = "bold"),
        axis.title.y = element_text(size = 13),
        axis.line = element_line(color = 'black'))

```

### fig 3.C) Expression profile at species level of the most abundant microorganisms using mOTUs 

```{r 3_mOTUs, fig.height= 7.5, fig.width=6}
# build phyloseq object
mOTUs.physeq <- phyloseq(otu_table(mOTUs.count, taxa_are_rows = TRUE),
                            sample_data(mOTUs.data),
                            tax_table(mOTUs.taxa))

# filter low abundant strains
minper = 3
Fil0 = genefilter_sample(mOTUs.physeq, filterfun_sample(function(x) x*100 >= minper))
mOTUs.physeq.F = prune_taxa(Fil0, mOTUs.physeq)

# plot
col.species = c("#67001F","#980043","#CE1256", "#E7298A", "#F7F4F9", "#E7E1EF", "#D4B9DA", "#C994C7", "#DF65B0",   
                "#BDBDBD", "#FFD92F", "#1B9E77", "#BF812D", "#969696", "#525252", "#B2DF8A", "#4292C6", "#FFF5F0")

mOTUs.physeq.F %>% psmelt() %>%
  group_by(Species, Day) %>%
  summarise(median = median(Abundance), .groups = 'drop') %>%
  ggplot(aes(x = Day, y = median, fill = Species)) +
  geom_bar(stat = "identity") +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = col.species,
                    limits = force, labels = c("Accumulibacter phosphatis",
                              "Accumulibacter regalis (BA-93)",
                              "Accumulibacter vicinus (SK-01)",
                              "Accumulibacter cognatus (SK-02)",
                              "Accumulibacter",
                              "Accumulibacter affinis",
                              "Accumulibacter delftensis",
                              "Accumulibacter propinquus",
                              "Accumulibacter proximus",
                              "Azonexus sp016714975",
                              "Propionivibrio aalborgensis",
                              "Burkholderiales",
                              "Flavobacterium",
                              "Kaistella",
                              "Kaistella beijingensis",
                              "Geothrix odensis",
                              "Micavibrionales",
                              "Pseudoxanthomonas")) +
  labs(y = "Expression profile (%)", x = "Days") +
  theme_minimal() +
  theme(legend.key.size = unit(0.5, "cm"), 
        legend.title = element_blank(),
        legend.text =  element_text(size = 10, face = "bold"),
        axis.text.y = element_text(size = 11, face = "bold"),
        axis.text.x = element_text(size = 11, face = "bold"),
        axis.title.y = element_text(size = 13),
        axis.line = element_line(color = 'black')) 

```

## Differential gene expression analysis using edgeR

The counts were summed at KO and genus level for each sample before DGE analysis.

```{r}
# here we combine KO info from DRAM and eggNOG, if the KO info is missing from DRAM then add the one from eggNOG
FC.annot.no0.lfc$KO = FC.annot.no0.lfc$ko_id
for(i in 1:nrow(FC.annot.no0.lfc)){
  if(is.na(FC.annot.no0.lfc$KO)[i]){FC.annot.no0.lfc$KO[i] = FC.annot.no0.lfc$KEGG_ko[i]}
}

# sum the count at genus and KO level
DGE.time <- melt(edgeR.list.lfc$counts) %>%
  left_join(FC.annot.no0.lfc, join_by(Var1 == X)) %>%
  group_by(Genus, KO, Var2) %>%
  summarise(Expression = sum(value), .groups = 'drop') %>% 
  filter(!is.na(KO), Genus != " NA")

# get the count table 
DGE.time.count <- dcast(DGE.time, formula = Genus+KO~Var2, value.var = "Expression") %>%
  add_column(Gene = paste0(.$Genus, "_", .$KO)) %>%
  select(! c("Genus", "KO")) %>%
  column_to_rownames(var = "Gene") 

# put in DGE list for edgeR
edgeR.list.DGEtime <- DGEList(counts = DGE.time.count, group = group)

# TMM normalisation 
edgeR.list.DGEtime.norm <- normLibSizes(edgeR.list.DGEtime, method = "TMM")

# define conditions that can be compared and the replicates
edgeR.design <- model.matrix(~0+group)
colnames(edgeR.design) <- levels(as.factor(group))

# define the comparisons we want to perform
edgeR.contrast <- makeContrasts(Fe26Fe13 = Feeding_26 - Feeding_13,
                                Fe103Fe26 = Feeding_103 - Feeding_26,
                                Fe182Fe103 = Feeding_182 - Feeding_103,
                                Ae26Ae13 = Aerobic_26 - Aerobic_13,
                                Ae103Ae26 = Aerobic_103 - Aerobic_26, 
                                Ae182Ae103 = Aerobic_182 - Aerobic_103, 
                                levels = edgeR.design)

# dispersion estimation
edgeR.disp <- estimateDisp(edgeR.list.DGEtime, edgeR.design)
plotBCV(edgeR.disp)

edgeR.model <- glmQLFit(edgeR.disp,edgeR.design)
plotQLDisp(edgeR.model)

```

```{r}
# vector with the name of the different contrasts
vec.contrast <- colnames(edgeR.contrast)

## using the previous vector, create for each contrast a table with all the results of the modelling
edgeR.test <- list() 
edger.test.adjust <- list()

for(i in seq_along(vec.contrast)) {
  edgeR.test[[i]] <- glmQLFTest(edgeR.model, contrast = edgeR.contrast[, vec.contrast[i]])
  names(edgeR.test)[i] <- vec.contrast[i]
  edger.test.adjust[[i]] <- topTags(edgeR.test[[i]], n = nrow(DGE.time.count))
  names(edger.test.adjust)[i] <- vec.contrast[i]
}

# define the value that will be used to determine if a gene is a DGE
logFCval = 0
pvalue = 0.01

# plot the volcan plot for each comparison
edgeR.DGE <- list()

for(i in 1:length(edgeR.test)) {
  df <- as.data.frame(edger.test.adjust[[i]]$table)
  title.graph <- names(edger.test.adjust[i])
  # define genes if they are UP are DOWN regulated based on Pvalue and logFC
  df$diffexpressed <- "NO"
  df$diffexpressed[df$logFC > logFCval & df$FDR < pvalue] <- "UP"
  df$diffexpressed[df$logFC < -logFCval & df$FDR < pvalue]<- "DOWN"
  edgeR.DGE[[i]] <- df
  names(edgeR.DGE)[i] <- vec.contrast[i]
  plotdf <- ggplot(data=df, aes(x=logFC, y=-log10(FDR), col=diffexpressed)) +
    geom_point() +
    scale_color_manual(values = c("blue", "black", "red")) +
    geom_vline(xintercept=c(-logFCval, logFCval), col="red") +
    geom_hline(yintercept=-log10(pvalue), col="red") +
    theme_minimal() +
    ggtitle(title.graph)
  print(plotdf)
}

```

### EBPR genes
The annotation file was merged with the list of genes of interest, to associate a metabolism to the genes, then filtered to keep only the genes with an associated metabolism before further analyses.

The annotation file here is the AppendixD_geneAnnotation.
*annotation of genes of interest* corresponds to FC.annot.no0.GoI.

For the figure, we used the count table after normalisation, so we have to recalculate to KO and genus level fro each sample.

```{r, fig.height=10, fig.width=7}
# get the normalised-transformed genes that correspond to a gene of interest for plotting
edgeR.list.norm.lfc.cpm.GoI <- edgeR.list.lfc.norm.cpm[unique(FC.annot.no0.GoI$X), ] 

# group information from count, annotation and metadata 
FC.ExprGoI <- melt(t(edgeR.list.norm.lfc.cpm.GoI)) %>% 
  left_join(FC.annot.no0.GoI, join_by(Var2 == X)) %>%
  left_join(FC.data, join_by(Var1 == Sample), relationship = "many-to-many") %>%
  filter(!GeneID == "no")

# order per day and phase
FC.ExprGoI$Day = factor(FC.ExprGoI$Day, levels = c("13", "26", "103", "182"))
FC.ExprGoI$Phase = factor(FC.ExprGoI$Phase, levels = c("Feeding", "Aerobic"))

```

#### EBPR genes check 

Here we plot the transcription of the EBPR genes for all genera. As it can observed, a lot of microbes transcribed really few genes in each metabolism. 

```{r, fig.height=10, fig.width=30}
# list of the metabolisms we want to look at
vec.metabo = unique((FC.genes.list %>% filter(Metabolism %in% c("Polyphosphate", "PHA", "Glycogen", "Volatile fatty acid")))$Metabolism)

# get the information for each genus and KO if DGE for feeding
test <- melt(t(data.frame(row.names = rownames(edgeR.DGE[[1]]),
                   Feeding_26 = edgeR.DGE[[1]]$diffexpressed,
                   Feeding_103 = edgeR.DGE[[2]]$diffexpressed,
                   Feeding_182 = edgeR.DGE[[3]]$diffexpressed)))

# # get the information for each genus and KO if DGE for Aerobic
# test <- melt(t(data.frame(row.names = rownames(edgeR.DGE[[1]]),
#                    Aerobic_26 = edgeR.DGE[[4]]$diffexpressed,
#                    Aerobic_103 = edgeR.DGE[[5]]$diffexpressed,
#                    Aerobic_182 = edgeR.DGE[[6]]$diffexpressed)))

# add the information of DGE to the global dataframe (count, annotation, metadata)
FC.ExprGoI.DGEtime <- FC.ExprGoI %>%
  add_column(MicKO = paste0(.$Genus, "_", .$KO)) %>%
  left_join(test, join_by(MicKO == Var2, Name == Var1))

# Simplify the dataframe to the metabolisms of interest 
## sum the value to geneID level and 
### calculate the median of the replicates (cycles)
FC.ExprGoI.DGEtime.sum <- FC.ExprGoI.DGEtime %>% 
  filter(Metabolism %in% vec.metabo) %>% 
  group_by(Genus, GeneID, Day, Phase, Metabolism, Cycle, value.y) %>% 
  summarise(Expression = sum(value.x), .groups = 'drop') %>%
  group_by(Genus, Day, Phase, Metabolism, GeneID, value.y) %>% 
  summarise(Expression = median(Expression), .groups = 'drop') 

# as we plot the result of the comparison from day x to day x-1, the first point is not represented (for plotting purpose we add day 13 instead of NA value)
FC.ExprGoI.DGEtime.sum$value.y <- FC.ExprGoI.DGEtime.sum$value.y %>% replace_na("13")

# order the geneIDs
FC.ExprGoI.DGEtime.sum$GeneID <- factor(FC.ExprGoI.DGEtime.sum$GeneID, levels = c("acs", "ackA", "pta", "actP", "prpE", 
                                                                  "glgA", "glgB", "glgC", "glgP", "glgX", "glgY", 
                                                                  "malQ", "malZ",
                                                                  "atoB", "phaA", "phaB", "phaC", "phaE", "phaZ", 
                                                                  "pap", "ppk2", "ppx", "ppk", 
                                                                  "phoAB", "phoD", "phoR", "phoU", 
                                                                  "pitA", "pstA", "pstB", "pstC", "pstS"))

# plot for either feeding or aerobic phases
FC.ExprGoI.DGEtime.sum %>% 
  filter(Phase == "Feeding") %>% 
  ggplot(aes(x = Genus, y = GeneID, fill = Day, group = Day)) +
  geom_point(aes(size = log(Expression+1), shape = value.y), position = position_dodge(0.8)) +
  scale_shape_manual("Differentially expressed\nbetween days", values = c(21, 25, 21, 24)) +
  scale_fill_manual(values = c("#DDCC77", "#999933", "#117733", "#661100")) +
  labs(y = "", x = "") +
  facet_grid(Metabolism~., scales = "free", space = "free") +
  guides(shape = guide_legend(override.aes = list(size = 4))) +
  guides(fill = guide_legend(override.aes = list(size = 4))) +
  scale_size(name = "Expression\nlog(cpm)") +
  theme_bw() +
  theme(legend.position = "",  
        legend.title = element_text(size = 12),
        legend.text =  element_text(size = 11),
        axis.text.y = element_text(size = 12, hjust = 0, face = "italic"),
        axis.text.x = element_text(size = 13,  angle = 45, hjust = 1),
        axis.title.y = element_text(size = 13),
        axis.line = element_line(color = 'black'),
        strip.text.y = element_text(angle = 0),
        strip.text = element_text(size = 12),
        strip.background = element_rect(fill = "white"))

```

#### Figure 4 and supp Figure 2. Expression of Enhanced Biological Phosphate Removal related genes at the feeding (aerobic) phase

The microorganisms with only one or two genes were not selected. To be considered the microorganisms should transcribed more than one or two genes in the different metabolisms.

```{r, fig.height=10, fig.width=9}
# selection of the microorganisms 
FC.ExprGoI.morg.list <- c(" Candidatus Accumulibacter", " Candidatus Competibacter", " Candidatus Contendobacter", " Propionivibrio", " Flavobacterium", " Azonexus")

FC.ExprGoI.DGEtime.sum %>% 
  filter(Phase == "Feeding") %>% #here replace by Aerobic to obtain the results from sup fig 3
  filter(Genus %in% FC.ExprGoI.morg.list) %>%
  ggplot(aes(x = Genus, y = GeneID, fill = Day, group = Day)) +
  geom_point(aes(size = log(Expression+1), shape = value.y), position = position_dodge(0.8)) +
  scale_shape_manual("Differentially Transcribed\nbetween days", values = c(21, 25, 21, 24)) +
  scale_fill_manual(values = c("#DDCC77", "#999933", "#117733", "#661100")) +
  labs(y = "", x = "") +
  scale_x_discrete(labels = c(" Candidatus Accumulibacter" = " Accumulibacter", 
                              " Candidatus Contendobacter" = " Contendobacter",
                              " Candidatus Competibacter" = " Competibacter")) +
  facet_grid(Metabolism~., scales = "free", space = "free") +
  guides(shape = guide_legend(override.aes = list(size = 4))) +
  guides(fill = guide_legend(override.aes = list(size = 4))) +
  scale_size(name = "Transcription\nlog(cpm)") +
  theme_bw() +
  theme(legend.position = "right",  
        legend.title = element_text(size = 12),
        legend.text =  element_text(size = 11),
        axis.text.y = element_text(size = 12, hjust = 0, face = "italic"),
        axis.text.x = element_text(size = 13,  angle = 45, hjust = 1),
        axis.title.y = element_text(size = 13),
        axis.line = element_line(color = 'black'),
        strip.text.y = element_text(angle = 0),
        strip.text = element_text(size = 12),
        strip.background = element_rect(fill = "white"))

```

### EPS genes

#### EPS check

```{r, fig.height=15, fig.width=30}
vec.metabo = unique((FC.genes.list %>% filter(Sub.metabolism %in% c("alginate", "Curli biogenesis", "nonulosonic acid pathways (NeuAc)", "nonulosonic acid pathways (leg)","nonulosonic acid pathways (pse)", "Pel operon", "PNAG", "Cellulose biosynthesis", "polysaccharide biosynthesis/export protein")))$Sub.metabolism)

FC.ExprGoI.DGEtime.sum <- FC.ExprGoI.DGEtime %>% 
  filter(Sub.metabolism %in% vec.metabo) %>% 
  group_by(Genus, GeneID, Day, Phase, Sub.metabolism, Cycle, value.y) %>% 
  summarise(Expression = sum(value.x), .groups = 'drop') %>%
  group_by(Genus, Day, Phase, Sub.metabolism, GeneID, value.y) %>% 
  summarise(Expression = median(Expression), .groups = 'drop') %>%
  filter(GeneID != "neuA")

FC.ExprGoI.DGEtime.sum$value.y <- FC.ExprGoI.DGEtime.sum$value.y %>% replace_na("13")

EPS.labels <- c("Alginate", "Curli biogenesis", "Legionaminic acid\npathway (leg)", 
                "N-acetylneuraminic aic\npathay (NeuAc)", "Pseudaminic acid\npathway (pse)", 
                "Pel operon", "poly-N-acetylglucosamine\n(PNAG)", 
                "Cellulose biosynthesis", "polysaccharide biosynthesis\nexport protein")
names(EPS.labels) <- c("alginate", "Curli biogenesis", "nonulosonic acid pathways (leg)", 
                "nonulosonic acid pathways (NeuAc)", "nonulosonic acid pathways (pse)", 
                "Pel operon", "PNAG", 
                "Cellulose biosynthesis", "polysaccharide biosynthesis/export protein")


FC.ExprGoI.DGEtime.sum %>% 
  filter(Phase == "Feeding") %>% 
  ggplot(aes(x = Genus, y = GeneID, fill = Day, group = Day)) +
  geom_point(aes(size = log(Expression+1), shape = value.y), position = position_dodge(0.8)) +
  scale_shape_manual("Differentially expressed\nfrom the day before", values = c(21, 25, 21, 24)) +
  scale_fill_manual(values = c("#DDCC77", "#999933", "#117733", "#661100")) +
  labs(y = "", x = "") +
  facet_grid(Sub.metabolism~., scales = "free", labeller = labeller(Sub.metabolism = EPS.labels), space = "free") +
  scale_size(name = "Expression\nlog(cpm)") +
  scale_x_discrete(labels = c(" Candidatus Accumulibacter" = " Accumulibacter", 
                              " Candidatus Contendobacter" = " Contendobacter")) +
  guides(shape = guide_legend(override.aes = list(size = 4))) +
  guides(fill = guide_legend(override.aes = list(size = 4))) +
  theme_bw() +
  theme(legend.position = "right",  
        legend.title = element_text(size = 12),
        legend.text =  element_text(size = 11),
        axis.text.y = element_text(size = 12, hjust = 0, face = "italic"),
        axis.text.x = element_text(size = 13,  angle = 45, hjust = 1),
        axis.title.y = element_text(size = 13),
        axis.line = element_line(color = 'black'),
        strip.text.y = element_text(angle = 0),
        strip.text = element_text(size = 12),
        strip.background = element_rect(fill = "white"))
```

#### Figure 5 and supp Figure 3. Expression of biofilm related genes at the feeding (aerobic) phase

```{r, fig.height=15, fig.width=12}
FC.ExprGoI.morg.list <- c(" Candidatus Accumulibacter", " Candidatus Contendobacter", " Propionivibrio", " Flavobacterium", " Azonexus", " Rhodoferax", " Zoogloea", " Kaistella")

FC.ExprGoI.DGEtime.sum %>% 
  filter(Phase == "Feeding") %>% #here replace ba Aerobic to obtain the results from sup fig 3
  filter(Genus %in% FC.ExprGoI.morg.list) %>%
  ggplot(aes(x = Genus, y = GeneID, fill = Day, group = Day)) +
  geom_point(aes(size = log(Expression+1), shape = value.y), position = position_dodge(0.8)) +
  scale_shape_manual("Differentially Transcription\nfrom the day before", values = c(21, 25, 21, 24)) +
  scale_fill_manual(values = c("#DDCC77", "#999933", "#117733", "#661100")) +
  labs(y = "", x = "") +
  facet_grid(Sub.metabolism~., scales = "free", labeller = labeller(Sub.metabolism = EPS.labels), space = "free") +
  scale_size(name = "Transcription\nlog(cpm)") +
  scale_x_discrete(labels = c(" Candidatus Accumulibacter" = " Accumulibacter", 
                              " Candidatus Contendobacter" = " Contendobacter")) +
  guides(shape = guide_legend(override.aes = list(size = 4))) +
  guides(fill = guide_legend(override.aes = list(size = 4))) +
  theme_bw() +
  theme(legend.position = "right",  
        legend.title = element_text(size = 12),
        legend.text =  element_text(size = 11),
        axis.text.y = element_text(size = 12, hjust = 0, face = "italic"),
        axis.text.x = element_text(size = 13,  angle = 45, hjust = 1),
        axis.title.y = element_text(size = 13),
        axis.line = element_line(color = 'black'),
        strip.text.y = element_text(angle = 0),
        strip.text = element_text(size = 12),
        strip.background = element_rect(fill = "white"))
```

```{r}

sessionInfo()

```