ReQON: a Bioconductor package for recalibrating quality scores from next-generation sequencing data

Input

The input to ReQON is an indexed and sorted BAM file [9] with quality scores reported in the QUAL field. The reads can be aligned to any genome using any alignment algorithm. The recalibration results will be dependent on the accuracy of these alignments.

Algorithm

ReQON uses logistic regression to recalibrate the quality scores. The following parameters must be specified:

• Region: Genomic coordinates of the region that the regression model is trained on. This training region must be large enough to obtain accurate coefficients for the regression model. On the other hand, specifying a larger than necessary region will increase run time and may overfit the model to the training set. We recommend training on one of the smaller chromosomes, or specifying MaxTrain, described next.

• MaxTrain (optional): This allows the user to train on a fixed number of bases from a large training region. For example, training on the first 10 million bases of the genomic region given in Region is achieved by setting MaxTrain = 10,000,000. Typically, results are consistent when training on at least 10 million bases and do not improve when the training size is larger than 25 million bases.

• RefSeq: File containing the reference sequence corresponding to the training set.

• SNP (optional): File of known variant locations to remove from the training set before recalibration.

• nerr: The maximum number of errors tolerated at a genomic position (default = 2). Positions with more than nerr errors may likely be true variants, so bases from these positions are removed from the training region.

• nraf: The maximum non-reference allele frequency at a genomic position that is allowed (default = 0.05). Positions with non-reference allele frequency greater than nraf are removed from the training set for the same reason as nerr.

The first step is to read the training region, specified by Region and MaxTrain, from the input BAM file. Positions that are listed in SNP are removed from the training set. Next, sequencing errors are identified as bases that do not match RefSeq. In reality, some of these identified bases will not be errors but instead correct calls, such as novel variants or mapping errors. In an attempt to remedy this issue, two different thresholds are set to remove positions most likely to contain false error calls. At each genomic position, thresh = max{nerr, nraf × coverage} is calculated. This threshold sets a maximum number of allowable called errors for positions with low coverage, and a maximum frequency of non-reference bases for positions with high coverage. The default settings allow up to 2 errors per position if the coverage is less than 60x, and more errors (0.05 × coverage) for positions with at least 60x coverage. Positions with more called errors than thresh are identified and removed from the training set. Note that when we say that “a position is removed from the training set,” this is different than keeping the bases at this position in the training set but switching the call from error to non-error. Instead, all bases at this position are removed from the training set due to low confidence of the error call.

Every base in the training set has now been classified as either error or non-error. There is a strong relationship between errors and their position in the read. Most errors occur at the end of the read, due to technical aspects of the sequencing process [10]. This same pattern is present in Figure 1A, which shows the distribution of sequencing errors in the training set by their read position. The rest of Figure 1 is described in more detail in the next section. To account for these high-error read positions in the ReQON regression model, the algorithm flags read positions that have more errors than a specified threshold, set at 1.5 times the average number of errors per read position. This threshold is visualized by the dashed cyan line in Figure 1A. For notational purposes, assume that there are k flagged positions, denoted as fp ₁ fp ₂ , … fp _k .

Regression coefficients (βs) are obtained from the following logistic regression model,$logit\left(\mathrm{Probability}\mathrm{of}Y=1\right)={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+{\beta }_4{X}_4+{\beta }_5{X}_5+{\beta }_6{X}_6+{\beta }_7{X}_7+{\beta }_8{X}_8+···+{\beta }_{7+k}{X}_{7+k}$,

where Y is the (n x 1) response vector and the X s are (n x 1) vectors of explanatory variables, as defined below.

• Y(i) = I{base i is a sequencing error (i.e., does not match RefSeq)}.

· X₁(i) = original quality score of base i.

• X₂(i) = I{X₁(i) = 0}. Many base calling algorithms assign a quality of zero to indicate bad or randomly called bases, so these bases are treated separately.

• X₃(i) = average quality score across all bases in the read containing base i. This helps identify reads for which all bases are assigned low quality scores.

• X₄(i) = Read position of base i.

• X₅(i) = I{base i = A}.

• X₆(i) = I{base i = C}.

• X₇(i) = I{base i = G}.

• X_7+m(i) = I{X₄(i) = fp _m } for m = 1, …, k. These indicator variables identify bases originating from one of the flagged read positions.

Additional covariates were originally considered, but those that did not improve model performance were dropped from the final model.

Due to memory constraints, it is not practical to fit the model when the training set contains tens of millions of bases. Instead, the training set is split into smaller subsets of no more than 10 million bases and the regression model is run on all subsets. Then, for each coefficient, the median is calculated over all regression models. Once these median regression coefficients are obtained, the model is applied to each base in the input BAM file to obtain predicted error probability values. These probabilities are transformed to the Phred scale [6] and rounded down to the nearest integer.

Output and visualizations

ReQON outputs a BAM file with original quality scores replaced by the recalibrated scores. Additionally, ReQON produces diagnostic plots (Figure 1) which show the effectiveness of the quality score recalibration on the training set. Plot A shows the distribution of errors in the training set by their read position. Any read position above the threshold (dashed cyan line) is given an additional indicator variable in the regression model, discussed in the previous section. Plot B shows the relative frequency distribution of quality scores both before (solid blue line) and after (dashed red line) recalibration.

The bottom two plots show the reported quality score versus the empirical quality score before (plot C) and after (plot D) recalibration. The empirical quality score is calculated by computing the observed error rate for all bases in the training set that are assigned a specific quality score. This error rate is then converted to the log-transformed Phred scale. If the quality scores are accurate, which occurs when the observed and reported sequencing error rates match, the points will fall on the 45-degree line.

The bottom plots also report the Frequency-Weighted Squared Error (FWSE), a measure of how close the points lie to the 45-degree line. FWSE is calculated by squaring the error (vertical distance between the point and 45-degree line), weighting this squared error by its relative frequency (shown in plot B and represented by shading of the point) and summing across all quality scores. FWSE will be close to zero if the quality scores accurately reflect the probability of a sequencing error.

Accuracy

As an example, two replicates of RNA from the U87 glioblastoma cell line [11] were sequenced using Illumina’s Genome Analyzer II, representing identical sequence runs but of slightly differing quality. For both cell line replicates, ReQON was run using the default settings with the model trained on chromosome 10.

Figure 1 shows the diagnostic output of ReQON after recalibration. Plot A shows that a majority of the sequencing errors occur at the ends of the read. Plot B shows that the majority of quality scores before recalibration were larger than 35, with almost 50% of the bases receiving a quality score near 40. After recalibration, the quality scores are assigned smaller values, with most quality scores falling between 25 and 35. Plot C confirms that the original quality scores are not very accurate because the quality scores with the largest frequencies (shaded dark blue) are far from the 45-degree line and, thus, FWSE is large. For example, plot B shows that approximately 50% of the bases are assigned a quality score of 40. In plot C, the empirical quality of this reported score is around 30. Thus, this single quality score contributes approximately 0.5(40 – 30)² = 50 to the total FWSE of 69.36. Plot D shows that, after recalibration, the quality scores do a much better job of representing true sequencing error rates, reflected by the 95% decrease in FWSE. Some of the lower quality scores lie away from the 45-degree line, but because few bases are assigned these scores (represented by unshaded circles), they contribute little to FWSE. The high-frequency quality scores (shaded dark red) are very close to the 45-degree line, resulting in a total FWSE that is very small.

Discrimination

Another desirable property of quality scores is the ability to separate sequencing errors from true variants. To perform this analysis, all bases that do not match RefSeq (hg 19) are analyzed. These non-reference bases are further classified as true variants (reported in dbSNP version 132 [12]) or sequencing errors (not in dbSNP132). Similar to the training set, bases identified as sequencing errors are removed if there are more than two errors at a position (or allele frequency greater than 0.05 for high coverage positions) as these may represent novel variants or systematic alignment errors. To increase confidence in the true variant calls, positions identified as true variants with less than 3 non-reference bases (or allele frequency less than 0.05 for high coverage positions) are also removed, as these may actually be sequencing errors.

Classification performance is measured by the area under the corresponding receiver operating characteristic (ROC) curve, or AUC. Figure 2 plots the relative frequency distributions of quality scores for non-reference bases in chromosome 20, an independent testing set. The red curve plots the distribution of bases belonging to positions in dbSNP (13,491 bases) and the blue curve plots the distribution of bases identified as sequencing errors (83,180 bases). Plot A shows the distribution of original quality scores before recalibration. The AUC for the original quality scores is 0.764, representing reasonable separation between the two classes. Plot B shows the distribution of quality scores after recalibration with ReQON. The two curves now have very different distributions. The quality scores of most bases at positions in dbSNP (red) have high quality scores (above 25). It is possible that many of the lower quality bases at positions in dbSNP may be sequencing errors, as bases belonging to any third alleles were not filtered out. In contrast, the quality scores of sequencing errors (blue) are mostly below 25. This better separation between the two classes is supported by the increased AUC (0.881 vs. 0.764). This analysis shows that, in addition to providing quality scores that more accurately reflect the probability of a sequencing error, the recalibrated quality scores also do a better job of distinguishing sequencing errors from true variants. The AUC statistics for both cell line replicates are shown in Table 2.

Comparison to GATK and BAQ

There are other available computational tools that recalibrate base quality scores. For example, some variant calling tools, such as Atlas-SNP2 [13] and SOAPsnp [4], recalibrate the quality scores as part of their algorithm. These tools incorporate the recalibrated qualities into their method, but do not output a file with the recalibrated scores. Therefore, we are not able to compare their performance to ReQON.

Figure 3 shows an example of a position in the training set where the bases do not match the reference sequence and the position is not listed in dbSNP version 132. Plot A visualizes the position and surrounding region using the Integrative Genomics Viewer (IGV) [14]. At the position, the reference sequence is a T, but all of the 103 bases mapped to this position are C (colored blue). This position most likely represents a novel variant or systematic mapping error. These bases are removed from the training set by ReQON because this position has a higher non-reference allele frequency (100%) than nraf (set at 5%). However, GATK calls all of these non-reference bases as sequencing errors when training its model. Plot B shows boxplots of the quality scores of the bases at this position after recalibrating with ReQON and GATK. The quality scores assigned by ReQON are significantly higher than the quality scores assigned by GATK (two-sided paired t-test, p = 8.36 × 10^-9).

In order to assess the accuracy of all three methods, Table 1 compares the FWSE statistic between ReQON, GATK and BAQ for two separate chromosomes, 10 and 20, on the two cell line replicates. Note that while GATK and BAQ train their models on the entire genome, ReQON only trained its model on chromosome 10. In this regard, chromosome 20 can be viewed as an independent testing set when assessing model performance for ReQON but not GATK or BAQ. Error calls were made in the same manner as the ReQON algorithm, which we feel more accurately identifies true sequencing errors, as described previously. In every case, FWSE is much lower for quality scores recalibrated with ReQON than GATK or BAQ. From this, we can conclude that the recalibrated quality scores from ReQON more accurately represent the probability that a base is a sequencing error.

In order to determine how well GATK and BAQ quality scores distinguish between non-reference bases classified as sequencing errors and non-reference bases belonging to positions of known variants, the discrimination analysis was repeated using the GATK and BAQ recalibrated quality scores. Table 2 shows the AUC for the original quality scores and the recalibrated quality scores. All of the recalibration methods outperform the original quality scores for both cell line replicates. The discrimination performance of GATK and ReQON is approximately the same, with GATK performing better in one replicate and ReQON performing better in the other replicate, and both outperforming BAQ.