Computational protocol: Agricultural Freshwater Pond Supports Diverse and Dynamic Bacterial and Viral Populations

Similar protocols

Protocol publication

[…] 16S rRNA gene reads were initially screened for low quality bases and short read lengths (), paired reads were merged using PANDAseq (), de-multiplexed, trimmed of artificial barcodes and primers, and assessed for chimeras using UCHIME in de novo mode, as implemented in Quantitative Insights Into Microbial Ecology (QIIME; version 1.9.1-20150604) (). Quality-controlled reads were clustered at 97% de novo into operational taxonomic units (OTUs) with the SILVA 16S database () in QIIME (). All sequences taxonomically assigned to chloroplasts were removed from further downstream analysis. When appropriate, data was normalized using metagenomeSeq’s cumulative sum scaling to account for uneven sampling depth ().To visualize the relative abundance of bacterial phyla, stacked bar charts were generated using ggplot2 (). In addition, bacterial taxa were summarized at the genera level in QIIME (level = 6) and those with a maximum relative abundance greater than 1% in at least one sample were used to build a heatmap via R (ver. 3.3.2) and vegan heatplus ().Significance tests were conducted using a Tukey’s HSD Test between filter size fractions and among sampling dates. Additionally, Pearson correlation coefficients were calculated to identify associations between the water characteristics and the relative abundance of the bacterial phyla/genera.Alpha diversity was calculated using the R packages: Bioconductor (), metagenomeSeq (), vegan (), phyloseq (), fossil (), biomformat (), and ggplot2 () on data rarefied to an even sampling depth (55,307 sequences) and tested for significance using a Tukey’s HSD Test.Beta diversity was determined through principle coordinates analysis (PCoA) plots of Bray–Curtis, Weighted and Unweighted UniFrac distances calculated using the R package phyloseq and tested for significance with ANOSIM (9,999 permutations) (; ; ; ).Core bacterial OTUs for each filter fraction, and the sample as a whole, were defined as OTUs in 100% of samples determined with QIIME’s compute_core_microbiome.py script (). Core OTUs were then visualized using ggplot2, and Krona (; ). [...] The paired-end reads were quality trimmed using Trimmomatic (ver. 0.36) (), merged with FLASh (ver. 1.2.11) (), and assembled de novo with metaSPAdes (ver. 3.10.1) without read error correction (). Taxonomic classifications were assigned to contigs by searching predicted peptide open reading frames (ORFs) against the peptide SEED and Phage SEED databases (retrieved 11/17/2017) (). Briefly, peptide ORFs were predicted from virome contigs using MetaGene () and were searched against the SEED and Phage SEED databases using protein-protein BLAST (BLASTp ver. 2.6.0+) (E value ≤ 1e-3) (). Taxonomy was assigned to contigs using ORF BLASTp hits to SEED sequences with NCBI taxonomy IDs. A contig is assigned the NCBI taxonomy ID with the maximum sum bit score across all ORF BLASTp hits in the contig. Scripts performing these assignments are available at https://github.com/dnasko/baby_virome.Functional classifications were assigned to ORFs by searching predicted peptide ORFs against the same peptide SEED database. Peptide ORFs were searched against the SEED databases using BLASTp (E value ≤ 1e-3). ORFs were assigned to a SEED subsystem with the maximum sum bit score among all of the ORF’s hits. Only functions associated with viral hits were considered for this analysis.Counts for taxonomy and functions identified in each virome are based on normalized contig and peptide ORF abundances, respectively. The abundance for each contig was estimated by recruiting all quality controlled reads to the assembled contigs using Bowtie2 (ver. 2.3.3) in very sensitive local mode (). Coverage for each contig was calculated using Samtools depth () and a custom parser available at https://github.com/dnasko/baby_virome. ORF abundances are calculated by computing the coverage of the contig within the ORF’s start and stop coordinates. To compare abundance measurements between viromes they are normalized to account for sequencing effort and assembly/recruitment proficiency. Briefly, each contig/ORF’s coverage is divided by the giga base pairs (Gbp) of reads able to recruit back to contigs/ORFs. Taxonomic and functional data were visualized using ggplot2 (). Additionally, Pearson correlation coefficients were calculated to identify the associations between the water characteristics and the abundance of the predicted viral taxa.For the phylogenetic marker gene, predicted ORFs were queried against Pol I UniRef90 () clusters using protein-protein BLAST (BLASTp) () with an E value cutoff ≤ 1e-5. Significant hits were filtered based on length (≥ 100 amino acids) and then confirmed to be Pol I via NCBI’s Conserved Domain BLAST online tool (). The sequences were then aligned with MAFFT using the FFT-NS-i × 1000 algorithm (). To obtain biologically meaningful data on the Pol I-containing phage, a region of interest (I547 to N923 in E. coli IAI39) containing the Phe762 position relative to E. coli IAI39 was selected and used to create an unrooted maximum likelihood tree with 100 bootstrap replicates using Geneious 6.0.5 () with PhyML (). Those with a Phe762 or Tyr762 were defined as generally virulent phage, while those with a Leu762 were defined as generally temperate (). Abundances for each Pol I peptide were calculated as described above. […]

Pipeline specifications