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