Please also reference the journal article:¶

Trimming and decontamination of metagenomic data can significantly impact assembly and binning metrics, phylogenomic and functional analysis¶

Jason M. Whitham and Amy M. Grunden, 2021¶

[email protected]✉ and [email protected]¶

North Carolina State University, 4550A Thomas Hall, Box 7615, Raleigh NC, 27695, United States of America¶

ABSTRACT¶

Background¶

Investigators using metagenomic sequencing to study their microbiomes are often provided data that has been trimmed and decontaminated or do it themselves without knowing the effect these procedures can have on their downstream analyses. Here we evaluated the impact that JGI trimming and decontamination procedures had on assembly and binning metrics, placement of metagenome assembled genomes into species trees, and functional profiles of metagenome-assembled genomes (MAGs) extracted from twenty three complex rhizosphere metagenomes. We also investigated how more aggressive trimming impacts these binning metrics.¶

Results¶

We found that JGI trimmed and decontamination of input reads had some significant impacts in assembly and binning metrics compared to raw reads, and that differences in placement of MAGs in species trees increased with decreasing completeness and contamination thresholds. More aggressive trimming beyond those used by JGI were found to reduce MAG counts.¶

Conclusions¶

Mild trimming and decontamination of metagenomics reads prior to assembly can change an investigator’s answer to the questions, “Who is there and what are they doing? However, mild trimming and decontamination of metagenomic reads with high quality scores is recommended for those who elect to do so.¶

INTRODUCTION¶

Metagenomics is the study of communities of microorganisms. Investigators in this field often ask the questions, "Who is there, and what are they doing?" They answer these questions by comparing microorganisms' genomes from one environment to ones found in other environments, in space and time. Differences in sequences of homologous universal genetic markers are used to construct phylogenetic (one marker) or phylogenomic (multiple marker) trees, enabling inference of identities, phenotypes, and functions of microorganisms based on those of neighboring organisms. Functional potential of community members is further evaluated by comparing non-universal genes to homologous genes of characterized enzymes, their domains, and the metabolic pathways they belong to. Leading up to these analyses, is a pipeline of trimming and decontamination, assembly, and binning steps.¶

Assemblies are long fragments of DNA, contigs, constructed from smaller fragments, known as reads, which are commonly acquired by shotgun sequencing of the DNA found in environmental samples. Since assemblies are made of contigs, some key metrics for evaluating assembly quality are total length of all assembly contigs, total number of assembly contigs, the number of bases of the largest contig in an assembly, and GC content [1]. When comparing metagenome assemblies generated with the same datasets, algorithm and parameters, assemblies with longer total lengths are generally better than ones with short total lengths. Conversely, total length can be inflated with misassemblies, and misassembly frequency can vary with different assemblers and datasets [2]. Having the largest contig relative to other assemblies is typically a favorable quality since longer contigs are good starting points for generating gap-free MAGs. Unlike the first two metrics mentioned, interpreting a high number for total contigs is a little less straightforward. An assembly with more total contigs relative to other assemblies could mean that it is larger, but it could also mean that it is more fragmented. So, the total contigs metric should be interpreted in the context metrics describing the size of contigs, such as the quantity of 1k, 10k, 100k bp contigs. GC content can be useful for simple assembly sequence comparisons.¶

Binning contigs into MAGs is a common step in characterizing metagenomes. Depending on the study, investigators might be interested in generating as many MAGs as possible or MAGs of a certain quality threshold in accordance with accepted recommendations [3]. The quality of MAGs is often defined by their completeness and contamination; each of which can be calculated by presence and absence of lineage-specific single-copy marker genes [4]. An important consideration when evaluating MAGs with this method is estimation bias. The real completeness is increasingly overestimated, and real contamination is increasingly underestimated with lower completeness and higher contamination estimates. The estimation bias is due to contaminant markers being counted as completeness markers. However, this bias is less than 3% when MAGs are estimated to be greater than 70% complete and less than or equal to 10% contaminated [4]. In other words, at these estimate limits, the MAGs might actually be 67% complete and 13% contaminated. This level of completeness is often good enough for investigators answering the questions "Who's there, and what are they doing?" Investigators considering quality in their MAG selection often even use recommended medium quality thresholds; completeness of greater than or equal to 50% and less than 10% contamination [3]. In comparison to the low level of bias found with 70% completeness and 10% contamination, one can expect as much as 8% bias with medium quality thresholds, resulting in some MAGs with 42% completeness and 18% contamination.¶

Prior to assembly and binning, investigators often perform trimming and decontamination of reads with the purpose of removing low quality bases, presumed errors generated in the sequencing process, contaminants like adaptors used in sequencing library preparation, and foreign DNA introduced during DNA extraction and library construction. The conventional wisdom is that trimming and decontamination of metagenomes are important to do before assembly and binning. The United States Department of Energy (DOE) Joint Genome Institute (JGI) follows this workflow and many bioinformatic forums like BioStars and SeqAnswers advocate for some level of trimming and decontamination, offering parameter suggestions for commonly used applications. However, we have not seen any published research that addresses the consequences of trimming and decontamination on assembly and binning metrics.¶

In this study, we evaluated how trimming and decontamination impacted total assembly length, total number of assembly contigs, largest contig length, the number of 1k, 10k, 100k bp contigs, GC content, total counts of bins, their completeness and contamination, those of "good” and “medium quality" bins. We did this with select data from a large sequencing effort completed by JGI in 2018 to enable investigators to identify beneficial ecofunctional genes selected by rhizospheres of three important biofuel crops - corn, switchgrass, and Miscanthus xgiganteus [5]. Rhizosphere samples were taken at seven plot areas at the Kellogg Biological Station (KBS) in Michigan during five seasons of crop establishment. Some of the data generated early on was used to create a pipeline for faster identification of short-subunit rRNA gene fragments and enabled unsupervised operational taxonomic unit community analysis with shotgun data [6]. Shotgun sequencing was compared to amplicon sequencing, demonstrating bias in the later for specific community members. More recently, the data was repurposed and combined with data from other sequencing efforts to evaluate the prevalence and activity of trace gas oxidizers of soil microbial communities [7]. In line with our research interests, we focused on metagenomes from Miscanthus xgiganteus, selecting 23 metagenomes that spanned all sampling seasons and plots and clustered well in principal component analysis (PCA) by taxon.¶

Researchers and students alike are increasingly interested in metagenomics. So, we made our analysis highly accessible to active and new metagenomics investigators by using KBase [8], an open science data platform funded by the DOE. This platform has a variety of command-line bioinformatics applications that have been simplified into modules that can be arranged into Jupitor Notebook [9] narratives. This format removes many of the learning curves associated with working in the terminal and provides visual displays for data and analyses, making research and findings more interoperable, reproducible, and repurposable. While the platform is not perfect, it is under heavy development, with increasing modular applications, and a funded help desk, with staff that are quick to respond to bugs and inquires. In addition to research, it is also being used to teach students in the classroom and asynchronously online, especially since the beginning of the coronavirus pandemic.¶

METHODS¶

Trimming and Decontamination¶

Select raw fastq files ("raw") Miscanthus xgiganteus rhizosphere metagenomes from the JGI Project "Metagenomic analysis of the rhizosphere of three biofuel crops at the KBS intensive site" (Proposal ID: 1296) (Tiedje 2013) were imported from JGI's Integrated Microbial Genomes and Microbiomes (IMG/M) website [10-11]. Selection of metagenomes was based on clustering in PCA plots by taxon, accomplished with the genome clustering tool on the JGI IMG/M website. Metadata for these metagenomes is provided in Table 1.¶

JGI trimmed and decontaminated reads (referred to as "qc" in this study) were imported directly into KBase for subsequent assembly and binning. JGI reports indicated that raw fastq file used in this study contained contaminant reads including PhiX (< 0.3%), DNA Spike-ins (< 3%), E. coli (< 0.7%), Mitocondria (< 4%), Chloroplast (< 4%), and rRNA (< 3%). The following parameters for trimming and decontamination were used to generate qc libraries from raw reads. BBDuk [12] adapter trimming (mink=11, k=23, hdist=1, hdist2=1, ktrim=r, tpe, tbo, minlen=40, minlenfraction=0.6, ftm=5) was used to remove known Illumina adapters. The reads were then processed using BBDuk filtering and trimming (maq=8, maxns=1, minlen=40, minlenfraction=0.6, k=27, hdist=1, trimq=12, qtrim=rl). At this stage reads ends were trimmed where quality values were less than 12. Read pairs containing more than three 'N', or with quality scores (before trimming) averaging less than 3 over the read, or length under 51bp after trimming, as well as reads matching illumina artifact, spike-ins and phiX were discarded. Remaining reads were mapped to a masked version of human HG19 with BBMap (fast local minratio=0.84 maxindel=6 tipsearch=4 bw=18 bwr=0.18 usemodulo printunmappedcount idtag minhits=1), discarding all hits over 93% identity.¶

Additional trimming and decontamination were performed on six randomly chosen raw fastq files with BBDuk (version 37.62) on Henry2, the high performance computing cluster at North Carolina State University, because BBDuk was not available on KBase at the time of this study. Six files were processed fourteen different ways with varying aggressiveness. Parameters used for each step were derived from author recommendations provided on the JGI website [13] and bioinformatics forums, SeqAnswers [14] and BioStars [15-16]. BBDuk trimmed adaptors of paired 151 bp raw reads; targeted kmers sized between 11 (mink=11) and 23 (k=23) and ones with a single substitution (hdist =1), all bases to the right of 3' adapters (ktrim=r), and paired reads to the same length in events where the adapter kmer was only detected in one of the paired reads (tpe). Paired overlapping reads were also merged (tpo). Reads that were quality trimmed were either trimmed within the recommend Q8-Q12 range (trimq=10) or to Q20 (trimq=20) from the left and right (qtrim=rl); reads shorter than 100 base pair were discarded (minlen=100). Illumina artifacts, spike-ins and phiX were also discarded from libraries that were decontaminated. Force trimming was also performed on some libraries since the last two (ftr=149) and first six (ftl = 7) were found to be low quality. Six files were chosen for additional trimming and decontamination in accordance with the minimum five to six groups recommended for mixed effects models [17].¶

Table 1: Data sources for this study. Starred(*) fastq files were randomly selected for BBDuk trimming and decontamination.¶

Assembly and Binning¶

Assembly and binning were performed in KBase. Currently, KBase has the leading metagenome assemblers, MEGAHIT [18], HipMER [19], SPAdes [20], and IDBA-UD [21] for assembling metagenomes. We used MEGAHIT in this study because it was the fastest, required the least amount of RAM, and had the longest assemblies [22]. KBase also currently has two of the leading binning tools, MetaBAT2 [23] and MaxBin2 [24]. Of these, we chose MetaBAT2 to bin all of our assemblies because it recently outperformed MaxBin2, generating more bins from the high-complexity Critical Assessment of Metagenome Interpretation dataset [25]. We binned some assemblies though with MaxBin2 to demonstrate similar results regardless of binning algorithm or min contig length.¶

Reads were assembled using MEGAHIT v1.2.9 with the meta-large parameter preset and a minimum contig length of 1000 bp. Assembly quality metrics including number of contigs, length of contigs, longest contig, and the number of contigs with sizes greater than 1k, 10k, and 100k bp were assessed with QUAST v4.4. MetaBAT2 v1.7 was used to bin minimum length contigs of 2500 bp. This minimum length was chosen because the MetBAT2 application occasionally failed when specifying shorter minimum contig lengths. Quality metrics including number of bins, bin completeness, bin contamination, single-copy and multi-copy markers were evaluated with CheckM v1.0.18.¶

Phylogenomic and Functional Analyses¶

Good quality MAGs that were greater than 70% complete and less than or equal to 10% contaminated were extracted from select raw and qc libraries as assembly sets with the KBase "Extract Bins as Assemblies from BinnedContigs" utility v.1.0.2. Assembly sets were annotated with RASTtk [26] v1.073. Annotated genome sets were inserted into species trees based on multiple universal phylogenetic markers with the KBase Species Tree application v2.2.0 [27]. Refseq (GCF) IDs were removed to condense tree figures, and the final trees visualized in ETE toolkit [28].¶

Functional assessment including COG [29], Pfam [30], and TIGRFAM [31] annotations of MAG genes was performed with the "View Function Profile for Genomes" application v1.4.0 MAGs and "Annotate Domains in a GenomeSet" utility v1.4.0. Defaults were used except that the lower limit threshold for TIGRFAM annotations was decreased from 10% to 9% of genes for the 9053.2 (M3-2) medium quality raw and qc genomesets to accommodate one pair of poorly annotated MAGs. Categories common to genomesets were statistically compared. Less than 0.2% of genes were not included in the analysis because they were absent in some genomesets.¶

Statistical Analyses, Data and Code Availability¶

Trimming, decontamination, assembly, binning, and functional analysis data used in this study was saved to Excel files. Python was used to combine Excel worksheets into Pandas dataframes, converted to lists, and then copy and pasted into KBase for statistical analysis in narrative code cells. Python coding was used for statistical analyses with the exception of calculating effect sizes, which were performed with the R package effsize [32]. The Excel files, data wrangling code, effect size and power calculation code are hosted at GitHub [33], and the code for all other statistical analyses are in the narratives [34] and [35] for complete transparency and reproducibility.¶

Paired Student's t-test and the nonparametric Wilcoxon signed-rank test were used to determine if there were significant differences in the means and mean signed-ranks of average and median assembly lengths, total contigs, the largest contig length, GC%, and the number contigs greater than or equal to 1k, 10k, 100k bp of raw and qc assemblies (α = 0.05). These tests were also used for binning metrics- MAG completeness, contamination, single-copy and multi-copy marker counts, and total MAG counts. Ordinary least squares and linear mixed effects models were used to determine significant correlations in read and base counts of raw and trimmed and decontaminated reads with MAG counts, completeness and contamination.¶

Assembly of raw and qc reads¶

Modules used in 9117.5*fastq processing¶

Modules used in 10158.8*fastq processing¶

Modules used in 11263.1*fastq processing¶

Modules used in 11306.3*fastq processing¶

Modules used in 11306.1*fastq processing¶

Modules used in 11260.6*fastq processing¶

Modules used in 11260.5*fastq processing¶

Modules used in 9053.2*fastq processing¶

Modules used in 9672.8*fastq processing¶

Modules used in 9053.4*fastq processing¶

Modules used in 9053.3*fastq processing¶

Modules used in 9053.5*fastq processing¶

Modules used in 9108.2*fastq processing¶

Modules used in 9108.2*fastq processing (raw assembly duplicate)¶

Modules used in 10186.3*fastq processing¶

Modules used in 10186.3*fastq processing (raw assembly duplicate)¶

Modules used in 7333.1*fastq processing¶

Modules used in 9117.8*fastq processing¶