On the optimal trimming of high-throughput mRNA sequence data

The widespread and rapid adoption of high-throughput sequencing technologies has afforded researchers the opportunity to gain a deep understanding of genome level processes that underlie evolutionary change, and perhaps more importantly, the links between genotype and phenotype. In particular, researchers interested in functional biology and adaptation have used these technologies to sequence mRNA transcriptomes of specific tissues, which in turn are often compared to other tissues, or other individuals with different phenotypes. While these techniques are extremely powerful, careful attention to data quality is required. In particular, because high-throughput sequencing is more error-prone than traditional Sanger sequencing, quality trimming of sequence reads should be an important step in all data processing pipelines. While several software packages for quality trimming exist, no general guidelines for the specifics of trimming have been developed. Here, using empirically derived sequence data, I provide general recommendations regarding the optimal strength of trimming, specifically in mRNA-Seq studies. Although very aggressive quality trimming is common, this study suggests that a more gentle trimming, specifically of those nucleotides whose Phred score <2 or <5, is optimal for most studies across a wide variety of metrics.


Introduction 1
The popularity of genome-enabled biology has increased dramatically over the last few years. While 2 researchers involved in the study of model organisms have had the ability to leverage the power of 3 genomics for nearly a decade, this power is only now available for the study of non-model organisms. 4 For many, the primary goal of these newer works is to better understand the genomic underpinnings of complexity of the graph, and thus increase both RAM requirements and runtime. 25 In addition to sequence error correction, which has been shown to improve accuracy of the de novo 26 assembly (MacManes and Eisen, 2013), low quality (=high probability of error) nucleotides are 27 commonly removed from the sequencing reads prior to assembly, using one of several available tools (http://www.biopieces.org/)). These tools typically use either a sliding window approach, 31 discarding nucleotides falling below a given (user selected) average quality threshold, or trimming of 32 low-quality nucleotides at one or both ends of the sequencing read. Though the absolute number will 33 surely be decreased in the trimmed dataset, aggressive quality trimming may remove a substantial 34 portion of the total read dataset, which in transcriptome studies may disproportionately effect lower 35 expression transcripts. 36 Although the process of nucleotide quality trimming is commonplace, particularly in the implementation has not been well defined. Though the rigor with which trimming is performed may be 41 guided by the design of the experiment, a deeper understanding of the effects of trimming is desirable.

42
As transcriptome-based studies of functional genomics continue to become more popular,  Here, I provide recommendations regarding the efficient trimming of high-throughput sequence reads, 51 specifically for mRNASeq reads from the Illumina platform. To do this, I used publicly available 52 datasets containing Illumina reads derived from Mus musculus. Subsets of these data (10 million, 20 53 million, 50 million, 75 million, 100 million reads) were randomly chosen, trimmed to various levels of 54 stringency, assembled then analyzed for assembly error and content. In addition to this, I develop a set 55 of metrics that may be generally useful in evaluating the quality of transcriptome assemblies. These 56 results aim to guide researchers through this critical aspect of the analysis of high-throughput 57 sequence data. While the results of this paper may not be applicable to all studies, that so many 58 researchers are interested in the genomics of adaptation and phenotypic diversity, particularly in 59 non-model organisms suggests its widespread utility.

61
Because I was interested in understanding the effects of sequence read quality trimming on the quality 62 of vertebrate transcriptome assembly, I elected to analyze a publicly available (SRR797058) paired-end

107
Quality trimming of sequence reads had a relatively large effect on the total number of errors 108 contained in the final assembly (Figure 1), which was reduced by between 9 and 26% when comparing 109 the assemblies of untrimmed versus Phred=20 trimmed sequence reads. Most of the improvement in 110 accuracy is gained when trimming at the level of Phred=5 or greater, with modest improvements 111 potentially garnered with more aggressive trimming at certain coverage levels (Table 1).

116
Trimming at Phred = 5 may be optimal, given the potential untoward effects of more stringent 117 quality trimming. 10M, 20M, 50M, 75M, 100M refer to the subsamples size. 10M replicate is the 118 technical replicate, 10M alt. dataset is the secondary dataset. Note that to enhance clarity, the 119 Y-axis does not start at zero.

120
In de Bruijn graph-based assemblers, the kmer is the fundamental unit of assembly. Even in 121 transcriptome datasets, unique kmers are likely to be formed as a results of sequencing error, and 122 therefore may be removed during the trimming process. Figure 2A shows the pattern of unique kmer 123 loss across the various trimming levels and read datasets. What is apparent, is that trimming at 124 Phred=5 removes a large fraction of unique kmers, with either less-or more-aggressive trimming 125 resulting in smaller effects. In contrast to the removal of unique kmers, those kmers whose frequency 126 is >1 are more likely to be real, and therefore should be retained. Figure 2B shows that while 127 Phred=5 removes unique kmers, it may also reduce the number of non-unique kmers, which may 128 hamper the assembly process. trimming and non-unique kmers, whose pattern is similar to that of unique kmers.

134
In addition to looking at nucleotide error and kmer distributions, assembly quality may be measured by 135 the the proportion of sequencing reads that map concordantly to a given transcriptome assembly 136 (Hunt et al., 2013). As such, the analysis of assembly quality includes study of the mapping rates.     Of note, all assembly files will be deposited in Dryad upon acceptance for publication. Until then, they 169 can be accessed via https://www.dropbox.com/sh/oiem0v5jgr5c5ir/TYQdGcpYwP.

171
Although the process of nucleotide quality trimming is commonplace in HTS analysis pipelines, 172 particularly those involving assembly, its optimal implementation has not been well defined. Though 173 the rigor with which trimming is performed seems to vary, there is a bias towards stringent trimming  based on this dataset contained more error (e.g. Figure 1), aggressive trimming did not improve quality 202 for any of the assessed metrics, though like other datasets, the absolute number of errors were reduced.

203
In addition to the specific dataset, the subsampling procedure may have resulted in undetected biases.

204
To address these concerns, a technical replicate of the original dataset was produced at the 10M 205 subsampling level. This level was selected as a smaller sample of the total dataset is more likely to with ultra-deep sequencing. Indeed, when looking at the differences in the magnitude of negative 223 effects in the datasets presented here, it is apparent that trimming at Phred=20 is 'less bad' in the 224 100M read dataset than it is in the 10M read datasets. For instance, Figure 2B demonstrates that one 225 of the untoward effects of trimming, the reduction of non-unique kmers, is reduced as the depth of 226 sequencing is increased. Figures 3 and 4  studies that revolve around assembly, quality control of sequence data has been thought to be a 243 crucial step. Though the removal of erroneous nucleotides is the goal, how best to accomplish this is 244 less clear. As described above, quality trimming has been a common method, but in its commonplace 245 usage, may be detrimental to assembly. What if, instead of relying on quality scores, we instead rely 246 on the distribution of kmers to guide our quality control endeavors? In transcriptomes of typical 247 complexity, sequenced to even moderate coverage, it is reasonable to expect that all but the most 248 exceptionally rare mRNA molecules are sequenced at a depth >1. Following this, all kmer whose 249 frequency is <2 are putative errors, and should be removed before assembly, though this process may 250 result in the loss of kmers from extremely low abundance transcripts or isoforms. This idea and its 251 implementation are fodder for future research.

252
In summary, the process of nucleotide quality trimming is commonplace in many HTS analysis 253 pipelines, but its optimal implementation has not been well defined. A very aggressive strategy, where 254 sequence reads are trimmed when Phred scores fall below 20 is common. My analyses suggest that 255 for studies whose primary goal is transcript discovery, that a more gentle trimming strategy (e.g.

256
Phred=2 or Phred=5) that removes only the lowest quality bases is optimal.