--- 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 ```{r} # # install package dada2 & phyloseq # if (!require("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # # BiocManager::install(c("dada2", "phyloseq", "tidyverse")) # download packages library(dada2) library(phyloseq) library(tidyverse) library(ShortRead) # cutadapt for preprocessing data cutadapt <- "PATH_to_cutadapt" # check version system2(cutadapt, args = "--version") ``` # Determination of Amplicon Sequence Variant (ASV) The script is adapted from: https://github.com/ErnakovichLab/dada2_ernakovichlab#learn-the-error-rates https://benjjneb.github.io/dada2/tutorial.html ## Sample preparation ### Filtered out adapters using cutadapt This step is done with *cutadapt*. ```{r} # Set path to shared data folder and contents data.fp <- "Path_to_raw_sequences" project.fp <- "path_to_output_folder" # CHANGE ME to project directory; don't append with a "/" # Set up names of sub directories to stay organized preprocess.fp <- file.path(project.fp, "01_preprocess") filtN.fp <- file.path(preprocess.fp, "filtN") trimmed.fp <- file.path(preprocess.fp, "trimmed") filter.fp <- file.path(project.fp, "02_filter") table.fp <- file.path(project.fp, "03_tabletax") ``` ```{r} # Get full paths for all files and save them for downstream analyses # Forward and reverse fastq filenames have format: fnFs <- sort(list.files(data.fp, pattern="_R1.fastq", full.names = TRUE)) fnRs <- sort(list.files(data.fp, pattern="_R2.fastq", full.names = TRUE)) # Ambiguous bases will make it hard for cutadapt to find short primer sequences in the reads. To solve this problem, we will remove sequences with ambiguous bases (Ns) # Name the N-filtered files to put them in filtN/ subdirectory fnFs.filtN <- file.path(preprocess.fp, "filtN", basename(fnFs)) fnRs.filtN <- file.path(preprocess.fp, "filtN", basename(fnRs)) # Filter Ns from reads and put them into the filtN directory filterAndTrim(fnFs, fnFs.filtN, fnRs, fnRs.filtN, maxN = 0, multithread = TRUE) # Set up the primer sequences to pass along to cutadapt FWD <- "AGMGTTYGATYMTGGCTCAG" ## primer F for v12 #REV <- "ACTCCTACGGGAGGCAGC" ## primer R for v12 REV <- "GCTGCCTCCCGTAGGAGT" # Write a function that creates a list of all orientations of the primers allOrients <- function(primer) { # Create all orientations of the input sequence require(Biostrings) dna <- DNAString(primer) # The Biostrings works w/ DNAString objects rather than character vectors orients <- c(Forward = dna, Complement = complement(dna), Reverse = reverse(dna), RevComp = reverseComplement(dna)) return(sapply(orients, toString)) # Convert back to character vector } # Save the primer orientations to pass to cutadapt FWD.orients <- allOrients(FWD) REV.orients <- allOrients(REV) FWD.orients # Write a function that counts how many time primers appear in a sequence primerHits <- function(primer, fn) { # Counts number of reads in which the primer is found nhits <- vcountPattern(primer, sread(readFastq(fn)), fixed = FALSE) return(sum(nhits > 0)) } # Before running cutadapt, we will look at primer detection for the first sample, as a check rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn = fnFs.filtN[[1]]), FWD.ReverseReads = sapply(FWD.orients, primerHits, fn = fnRs.filtN[[1]]), REV.ForwardReads = sapply(REV.orients, primerHits, fn = fnFs.filtN[[1]]), REV.ReverseReads = sapply(REV.orients, primerHits, fn = fnRs.filtN[[1]])) # run cutadapt # Create directory to hold the output from cutadapt if (!dir.exists(trimmed.fp)) dir.create(trimmed.fp) fnFs.cut <- file.path(trimmed.fp, basename(fnFs)) fnRs.cut <- file.path(trimmed.fp, basename(fnRs)) # Save the reverse complements of the primers to variables FWD.RC <- dada2:::rc(FWD) REV.RC <- dada2:::rc(REV) ## Create the cutadapt flags ## # Trim FWD and the reverse-complement of REV off of R1 (forward reads) R1.flags <- paste("-g", FWD, "-a", REV.RC, "--minimum-length 50") # Trim REV and the reverse-complement of FWD off of R2 (reverse reads) R2.flags <- paste("-G", REV, "-A", FWD.RC, "--minimum-length 50") # Run Cutadapt for (i in seq_along(fnFs)) { system2(cutadapt, args = c("-j", 0, R1.flags, R2.flags, "-n", 2, # -n 2 required to remove FWD and REV from reads "-o", fnFs.cut[i], "-p", fnRs.cut[i], # output files fnFs.filtN[i], fnRs.filtN[i])) # input files } # As a sanity check, we will check for primers in the first cutadapt-ed sample: ## should all be zero! rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn = fnFs.cut[[1]]), FWD.ReverseReads = sapply(FWD.orients, primerHits, fn = fnRs.cut[[1]]), REV.ForwardReads = sapply(REV.orients, primerHits, fn = fnFs.cut[[1]]), REV.ReverseReads = sapply(REV.orients, primerHits, fn = fnRs.cut[[1]])) ``` ### Inspect reads quality All the following steps are done using DADA2. ```{r} # Put filtered reads into separate sub-directories for big data workflow dir.create(filter.fp) subF.fp <- file.path(filter.fp, "preprocessed_F") subR.fp <- file.path(filter.fp, "preprocessed_R") dir.create(subF.fp) dir.create(subR.fp) # Move R1 and R2 from trimmed to separate forward/reverse sub-directories fnFs.Q <- file.path(subF.fp, basename(fnFs)) fnRs.Q <- file.path(subR.fp, basename(fnRs)) file.copy(from = fnFs.cut, to = fnFs.Q) file.copy(from = fnRs.cut, to = fnRs.Q) # File parsing; create file names and make sure that forward and reverse files match filtpathF <- file.path(subF.fp, "filtered") # files go into preprocessed_F/filtered/ filtpathR <- file.path(subR.fp, "filtered") # ... fastqFs <- sort(list.files(subF.fp, pattern="fastq")) fastqRs <- sort(list.files(subR.fp, pattern="fastq")) if(length(fastqFs) != length(fastqRs)) stop("Forward and reverse files do not match.") ``` ### visualizing the quality profiles of the forward reads ```{r, fig.width= 8, fig.height=6} # If the number of samples is 20 or less, plot them all, otherwise, just plot 20 randomly selected samples if( length(fastqFs) <= 20) { fwd_qual_plots <- plotQualityProfile(paste0(subF.fp, "/", fastqFs)) rev_qual_plots <- plotQualityProfile(paste0(subR.fp, "/", fastqRs)) } else { rand_samples <- sample(size = 20, 1:length(fastqFs)) # grab 20 random samples to plot fwd_qual_plots <- plotQualityProfile(paste0(subF.fp, "/", fastqFs[rand_samples])) rev_qual_plots <- plotQualityProfile(paste0(subR.fp, "/", fastqRs[rand_samples])) } fwd_qual_plots ``` ** Comments** In gray-scale is a heat map of the frequency of each quality score at each base position. The mean quality score at each position is shown by the green line, and the quartiles of the quality score distribution by the orange lines. The red line shows the scaled proportion of reads that extend to at least that position (this is more useful for other sequencing technologies, as Illumina reads are typically all the same length, hence the flat red line) ### visualize the quality profile of the reverse reads ```{r, fig.width= 8, fig.height=6} rev_qual_plots ``` ### Filter and trim ```{r} filt_out <- filterAndTrim(fwd=file.path(subF.fp, fastqFs), filt=file.path(filtpathF, fastqFs), rev=file.path(subR.fp, fastqRs), filt.rev=file.path(filtpathR, fastqRs), truncLen=c(225,225), maxEE=c(2,2), truncQ=2, maxN=0, rm.phix=TRUE, compress=TRUE, verbose=TRUE, multithread=TRUE) # look at how many reads were kept head(filt_out) # summary of samples in filt_out by percentage filt_out %>% data.frame() %>% mutate(Samples = rownames(.), percent_kept = 100*(reads.out/reads.in)) %>% select(Samples, everything()) %>% summarise(min_remaining = paste0(round(min(percent_kept), 2), "%"), median_remaining = paste0(round(median(percent_kept), 2), "%"), mean_remaining = paste0(round(mean(percent_kept), 2), "%"), max_remaining = paste0(round(max(percent_kept), 2), "%")) ``` ```{r} if( length(fastqFs) <= 20) { remaining_samplesF <- fastqFs[ which(fastqFs %in% list.files(filtpathF))] # keep only samples that haven't been filtered out remaining_samplesR <- fastqRs[ which(fastqRs %in% list.files(filtpathR))] # keep only samples that haven't been filtered out fwd_qual_plots_filt <- plotQualityProfile(paste0(filtpathF, "/", remaining_samplesF)) rev_qual_plots_filt <- plotQualityProfile(paste0(filtpathR, "/", remaining_samplesR)) } else { remaining_samplesF <- fastqFs[rand_samples][ which(fastqFs[rand_samples] %in% list.files(filtpathF))] # keep only samples that haven't been filtered out remaining_samplesR <- fastqRs[rand_samples][ which(fastqRs[rand_samples] %in% list.files(filtpathR))] # keep only samples that haven't been filtered out fwd_qual_plots_filt <- plotQualityProfile(paste0(filtpathF, "/", remaining_samplesF)) rev_qual_plots_filt <- plotQualityProfile(paste0(filtpathR, "/", remaining_samplesR)) } fwd_qual_plots_filt ``` ```{r} rev_qual_plots_filt ``` ## Error rates ```{r} # File parsing filtFs <- list.files(filtpathF, pattern="fastq", full.names = TRUE) filtRs <- list.files(filtpathR, pattern="fastq", full.names = TRUE) # Sample names in order sample.names <- basename(filtFs) # doesn't drop fastq.gz sample.names <- gsub("_R1.fastq", "", sample.names) sample.namesR <- basename(filtRs) # doesn't drop fastq.gz sample.namesR <- gsub("_R2.fastq", "", sample.namesR) # Double check if(!identical(sample.names, sample.namesR)) stop("Forward and reverse files do not match.") names(filtFs) <- sample.names names(filtRs) <- sample.names ``` ```{r} set.seed(100) # set seed to ensure that randomized steps are replicatable ``` ### Learning error rates In the plots below, the red line represents what we should expect the learned error rates to look like for each of the 16 possible base transitions (A->A, A->C, A->G, etc.) and the black line and grey dots represent what the observed error rates are. If the black line and the red lines are very far off from each other, it may be a good idea to increase the nbases parameter. This allows the machine learning algorithm to train on a larger portion of your data and may help improve the fit. ```{r} # Learn forward error rates (Notes: randomize default is FALSE) errF <- learnErrors(filtFs, nbases = 1e14, multithread = TRUE, randomize = TRUE) # Learn reverse error rates errR <- learnErrors(filtRs, nbases = 1e14, multithread = TRUE, randomize = TRUE) saveRDS(errF, paste0(filtpathF, "/errF.rds")) saveRDS(errR, paste0(filtpathR, "/errR.rds")) ``` ```{r} # Original default recommended way (not optimal for NovaSeq data!) errF_plot <- plotErrors(errF, nominalQ = TRUE) errR_plot <- plotErrors(errR, nominalQ = TRUE) errF_plot errR_plot ``` ## Dereplication, sequence inference ```{r} #For complex communities when you want to preserve rare taxa alternative: swap pool = TRUE with pool = "pseudo" # Dereplicate forward reads derepF.p <- derepFastq(filtFs) names(derepF.p) <- sample.names # Infer sequences for forward reads dadaF.p <- dada(derepF.p, err = errF, multithread = TRUE, pool = "pseudo") names(dadaF.p) <- sample.names # Dereplicate reverse reads derepR.p <- derepFastq(filtRs) names(derepR.p) <- sample.names # Infer sequences for reverse reads dadaR.p <- dada(derepR.p, err = errR, multithread = TRUE, pool = "pseudo") names(dadaR.p) <- sample.names ``` ## Merge paired reads ```{r} # Merge reads together mergers <- mergePairs(dadaF.p, derepF.p, dadaR.p, derepR.p) ``` ## Construct sequence table ```{r} seqtab <- makeSequenceTable(mergers) dim(seqtab) # Inspect distribution of sequence lengths table(nchar(getSequences(seqtab))) seqtab2 <- seqtab[,nchar(colnames(seqtab)) %in% 150] table(nchar(getSequences(seqtab2))) # Save table as an r data object file dir.create(table.fp) saveRDS(seqtab, paste0(table.fp, "/seqtab.rds")) ``` ## Remove chimera ```{r} # Read in RDS st.all <- readRDS(paste0(table.fp, "/seqtab.rds")) # Remove chimeras seqtab.nochim <- removeBimeraDenovo(st.all, method="consensus", multithread=TRUE) save(seqtab.nochim, file = "./seqtab.nochim.RData") # Print percentage of our seqences that were not chimeric. 100*sum(seqtab.nochim)/sum(seqtab) ``` ## Track reads through the pipeline ```{r} getN <- function(x) sum(getUniques(x)) # function to grab sequence counts from output objects # tracking reads by counts filt_out_track <- filt_out %>% data.frame() %>% mutate(Sample = gsub("_R1.fastq","",rownames(.))) %>% dplyr::rename(input = reads.in, filtered = reads.out) rownames(filt_out_track) <- filt_out_track$Sample ddF_track <- data.frame(denoisedF = sapply(dadaF.p[sample.names], getN)) %>% mutate(Sample = row.names(.)) ddR_track <- data.frame(denoisedR = sapply(dadaR.p[sample.names], getN)) %>% mutate(Sample = row.names(.)) merge_track <- data.frame(merged = sapply(mergers, getN)) %>% mutate(Sample = row.names(.)) chim_track <- data.frame(nonchim = rowSums(seqtab.nochim)) %>% mutate(Sample = row.names(.)) track <- left_join(filt_out_track, ddF_track, by = "Sample") %>% left_join(ddR_track, by = "Sample") %>% left_join(merge_track, by = "Sample") %>% left_join(chim_track, by = "Sample") %>% replace(., is.na(.), 0) %>% select(Sample, everything()) row.names(track) <- track$Sample head(track) # tracking reads by percentage track_pct <- track %>% data.frame() %>% mutate(Sample = rownames(.), filtered_pct = ifelse(filtered == 0, 0, 100 * (filtered/input)), denoisedF_pct = ifelse(denoisedF == 0, 0, 100 * (denoisedF/filtered)), denoisedR_pct = ifelse(denoisedR == 0, 0, 100 * (denoisedR/filtered)), merged_pct = ifelse(merged == 0, 0, 100 * merged/((denoisedF + denoisedR)/2)), nonchim_pct = ifelse(nonchim == 0, 0, 100 * (nonchim/merged)), total_pct = ifelse(nonchim == 0, 0, 100 * nonchim/input)) %>% select(Sample, ends_with("_pct")) # summary stats of tracked reads averaged across samples track_pct_avg <- track_pct %>% summarize_at(vars(ends_with("_pct")), list(avg = mean)) head(track_pct_avg) track_pct_med <- track_pct %>% summarize_at(vars(ends_with("_pct")), list(avg = stats::median)) head(track_pct_med) # Plotting each sample's reads through the pipeline track_plot <- track %>% data.frame() %>% mutate(Sample = rownames(.)) %>% gather(key = "Step", value = "Reads", -Sample) %>% mutate(Step = factor(Step, levels = c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim"))) %>% ggplot(aes(x = Step, y = Reads)) + geom_line(aes(group = Sample), alpha = 0.2) + geom_point(alpha = 0.5, position = position_jitter(width = 0)) + stat_summary(fun.y = median, geom = "line", group = 1, color = "steelblue", size = 1, alpha = 0.5) + stat_summary(fun.y = median, geom = "point", group = 1, color = "steelblue", size = 2, alpha = 0.5) + stat_summary(fun.data = median_hilow, fun.args = list(conf.int = 0.5), geom = "ribbon", group = 1, fill = "steelblue", alpha = 0.2) + geom_label(data = t(track_pct_avg[1:5]) %>% data.frame() %>% dplyr::rename(Percent = 1) %>% mutate(Step = c("filtered", "denoisedF", "denoisedR", "merged", "nonchim"), Percent = paste(round(Percent, 2), "%")), aes(label = Percent), y = 1.1 * max(track[,2])) + geom_label(data = track_pct_avg[6] %>% data.frame() %>% dplyr::rename(total = 1), aes(label = paste("Total\nRemaining:\n", round(track_pct_avg[1,6], 2), "%")), y = mean(track[,6]), x = 6.5) + expand_limits(y = 1.1 * max(track[,2]), x = 7) + theme_classic() track_plot # Write results to disk saveRDS(track, paste0(project.fp, "/tracking_reads.rds")) saveRDS(track_pct, paste0(project.fp, "/tracking_reads_percentage.rds")) saveRDS(track_plot, paste0(project.fp, "/tracking_reads_summary_plot.rds")) ggsave(plot = track_plot, filename = paste0(project.fp, "/tracking_reads_summary_plot.png"), width = 10, height = 10, dpi = "retina") ``` ## Assign taxonomy Tha MIDAS database can be download in the following link: https://www.midasfieldguide.org/guide/downloads The database used for this analysis was MIDAS 5.3 ```{r} taxa <- assignTaxonomy(seqtab.nochim, "Path_to_MIDAS 5 DB", multithread=TRUE) taxa.print <- taxa # Removing sequence rownames for display only rownames(taxa.print) <- NULL head(taxa.print) ``` ```{r} # Flip table seqtab.t <- as.data.frame(t(seqtab.nochim)) # Pull out ASV repset rep_set_ASVs <- as.data.frame(rownames(seqtab.t)) rep_set_ASVs <- mutate(rep_set_ASVs, ASV_ID = 1:n()) rep_set_ASVs$ASV_ID <- sub("^", "ASV_", rep_set_ASVs$ASV_ID) rep_set_ASVs$ASV <- rep_set_ASVs$`rownames(seqtab.t)` rep_set_ASVs$`rownames(seqtab.t)` <- NULL # Add ASV numbers to table rownames(seqtab.t) <- rep_set_ASVs$ASV_ID # Add ASV numbers to taxonomy taxonomy <- as.data.frame(taxa) taxonomy$ASV <- as.factor(rownames(taxonomy)) taxonomy <- merge(rep_set_ASVs, taxonomy, by = "ASV") rownames(taxonomy) <- taxonomy$ASV_ID taxonomy_for_mctoolsr <- unite_(taxonomy, "taxonomy", c("Kingdom", "Phylum", "Class", "Order","Family", "Genus", "ASV_ID"), sep = ";") # Write repset to fasta file # create a function that writes fasta sequences writeRepSetFasta<-function(data, filename){ fastaLines = c() for (rowNum in 1:nrow(data)){ fastaLines = c(fastaLines, as.character(paste(">", data[rowNum,"name"], sep = ""))) fastaLines = c(fastaLines,as.character(data[rowNum,"seq"])) } fileConn<-file(filename) writeLines(fastaLines, fileConn) close(fileConn) } # Arrange the taxonomy dataframe for the writeRepSetFasta function taxonomy_for_fasta <- taxonomy %>% unite("TaxString", c("Kingdom", "Phylum", "Class", "Order","Family", "Genus", "ASV_ID"), sep = ";", remove = FALSE) %>% unite("name", c("ASV_ID", "TaxString"), sep = " ", remove = TRUE) %>% select(ASV, name) %>% dplyr::rename(seq = ASV) # write fasta file writeRepSetFasta(taxonomy_for_fasta, paste0(table.fp, "/repset.fasta")) # Merge taxonomy and table seqtab_wTax <- merge(seqtab.t, taxonomy_for_mctoolsr, by = 0) seqtab_wTax$ASV <- NULL # Set name of table in mctoolsr format and save out_fp <- paste0(table.fp, "/seqtab_wTax_mctoolsr.txt") names(seqtab_wTax)[1] = "#ASV_ID" write("#Exported for mctoolsr", out_fp) suppressWarnings(write.table(seqtab_wTax, out_fp, sep = "\t", row.names = FALSE, append = TRUE)) # Also export files as .txt write.table(seqtab.t, file = paste0(table.fp, "/seqtab_final.txt"), sep = "\t", row.names = TRUE, col.names = NA) write.table(taxonomy[, 3:ncol(taxonomy)], "./tax_final.csv") ```