Ion Torrent's Data Quality Is Pretty Good (and Better Than Ion Claims)

One of the key questions around Ion Torrent, as with all new platforms, is what is the sequence data quality like.  Now, that can be a loaded question but I'll ask a slight variant on it: how truthful (or accurate) is Ion at estimating base quality?

Quality scores are a useful adjunct to sequencing data and are commonly expressed as phred scores, which are the integer part of -10*log10 of the error probability.  Any base caller needs to estimate these and many downstream programs, from aligners to assemblers to variant callers, rely on these quality values for their operations.  In many cases, the individual quality scores are combined to generate some joint estimate of the error (I built one such model at Codon).  These error probabilities come not from an infallible source, but are rather estimated from aspects of the raw data.  

Now, if we have a sequencing dataset for which we know truth, we can compare the actual error rates with those predicted by the quality scores.  Ion Torrent uses sequencing of E.coli DH10B as such a dataset.   DH10B is useful for this because (assuming perfection in the aligner) it should have no polymorphisms; deviations from the reference should all be errors.  I've had access to a few such datasets for over a month, but from a third party who made me swear not to share any details about them.  Luckily, Ion now has on their site what is apparently a recent run on the 314 chip with well over 500K reads.  In addition to the raw data, they provide alignments generated with their TMAP program.

I'm going to spew out a bunch of plots, but first I'll describe them.  Each plot is comparing the difference between the observed quality score and the called quality score.  If this delta is positive, it means the real quality is better than observed.  I call this case "pessimistic prediction"; the base caller was undercalling the quality.  On the other hand, if the value is negative then we have an "optimistic prediction"; the caller thinks it is doing better than it is.  A perfect caller would always score 0.   One final note: to compute the error frequencies a pseudocount of 1 was added to the error counts; otherwise bins with no observed errors blow up in the log calculation.  

In the plots, the X-axis is sequence read position and the Y-axis is the quality delta value.  Each plot is further trellised by the quality value.  Now, not only does this show the delta plots, we can also see a lot about the quality profile across a run.  Also, only bins with at least 200 events are shown (and point sizes are proportional to the log2 of the number of events).

So first, let's look at this plot for TMAP.  The first thing I noted is that most of the values are above 0; the caller is being too pessimistic.  The caller is often being extremely pessimistic; delta values are often 20 or more.  So, for example, this means that some of the lowest quality calls (qs=3) near the beginning of the read are actually more like 15-25.  The high end is less affected, but still so; the maximum quality of 33 really should be more like 36 in places.
A curious bit is the structure of some of these plots.  I would have expected smoother curves, though perhaps that is naive.  Some have discontinuities or multiple peaks.  Presumably much of this structure are artifacts of the base caller's quality prediction algorithm.
Also note that very high quality values are seen only at certain parts of the run: the highest quality value is 33 and isn't seen after position 32 in a run.

Is Ion's caller ever overestimating quality?  This is easier to see with the plot filtered in Spotfire so that delta values greater than -1 are removed (above; -1 is pretty much accounted for by the rounding error in generating integer quality scores).  According to these plots it is rare and mostly (and is most intense) at the beginning of the reads.  That's a suspicious position, which we'll look at more below.  The one exception to this pattern is quality score 30 at position 24, which is is really more like 28.  Not off by much and perhaps this is just the sort of noise this fine tooth comb will find.

Now, when I saw the same beginning of read effect in the earlier data, it was suggested that maybe this was an alignment issue related to homopolymer runs; a gap in the beginning of the sequence might not be detected, causing the beginning of the read to be misaligned.  Indeed, I was using a non-gapping aligner (Bowtie) in that first analysis & a little customer end-adjuster confirmed that inserting a few gaps could drastically improve things.  What if we try a better aligner?  TMAP is meant to be optimized for Ion Torrent data, but is it?

My aligner I use in this situation is BWA-SW.  On this dataset, it aligns slightly fewer reads than TMAP (485707 vs.508989), or about 3% (there's also 5 reads aligned by BWA-SW but not TMAP)..  The plots are now remarkably better (BWA-SW in gold).  The end effect goes away entirely, and the plots overall are in BWA-SW's favor; the gold dots are nearly always above the blue dots.
It looks like it could be a major project to really scrutinize the differences between the TMAP and BWASW alignments.  Spot checking a few, I am seeing some TMAP alignments that are curiously aggressive. For example, in one the initial base is matched and then the next 4 bases of the read are deleted, followed by 95 matching bases; BWA-SW calls for the first base to be skipped (unaligned) and then 95 bases of match.  Another small observation is that in the BWA-SW alignments which leave the first base unaligned and it is part of a homopolymer run, an extra T seems to be called more often than an extra G and each more than an extra A or C.  I was guessing extra Gs would be most common, given that the last base of the key sequence is a G and there might be a problem counting after G.  The pattern holds for cases where BWA-SW thinks the first base should be skipped but it is not followed by an identical base; but over 50% of the time it is a T which should be skipped over and the pattern holds even if you remove homopolymers.  Very strange; one guess is that this is an artifact of the library preparation protocol (which I am not familiar with in depth) and that the T really is there, but didn't come from DH10B.

One last observation.  This dataset had 260 flows (including the 4-base key sequence) and a handful of reads were as long as 112, but there were few of them. With BWA-SW, the observed quality scores at positions of 100 or greater, for those bins with at least 500 counts, were well above 20.  This would argue that the run was really stopped well short of its maximum potential; additional flows would have yielded quality data (though  I won't attempt to estimate how many would make sense)
This are some preliminary observations; it wouldn't shock me if I find some condition I haven't accounted for which contributes to this and which isn't interesting.  Also, the impact of these observations on an experiment will depend on the experiment.  For example, in a non-barcoded amplicon experiment the precise quality of first 20 or so bases is nearly irrelevant, as these are part of the targeting sequences and therefore any variants seen are far more likely to be oligonucleotide synthesis errors than real polymorphisms.  For a barcoded design, on the other hand, the barcode is likely in that beginning stretch and must be designed to account for these issues.  Fusion primer designs would avoid the "skip T" issue described above, if that is an artifact of the ligation process.  On the other hand, if that really is an artifact of the sequencing platform it will show up again.  Applications such as de novo assembly, which have the least information to diagnose problems, should be the most cautious.   In any case, a good rule is to make sure any variant you are calling is supported by a diversity of read positions; variants always called by the beginning or end of a read are suspect (I've seen this in another experiment).