There is more (length) to Ion Torrent reads than meets the eye (and is Ion Torrent hiding it?)

A quick summary for the impatient:

The sff files from the E coli Ion Torrent runs released by EdgeBio show much longer raw reads than the trimmed reads in the corresponing fastq/fasta files. The quality of those extra bases, however, is very low. This shows the potential for longer reads from the Ion Torrent platform.
The sff file released by Ion Torrent through their Dev Community site has these extra bases masked, which makes one wondering what if are trying to hide something…

Part 1: EdgeBio’s data
When EdgeBio released six runs with E coli DH10B Ion Torrent data (see http://www.edgebio.com/blog/?p=191), I decided to have a look inside the sff files they provided. I downloaded the data from data.edgebio.com, and used Roche’s sffinfo command to ‘peek inside’. The sffinfo command, accompanying the 454 Life Science software suite, will list the content of the binary sff file in text format (see the post on my other blog). Other, open source/access tools, such the ones my mention on my blog, might do this as well.
Here are the ‘header’ (manifest) part and the first read of DH10B-run01:

Common Header:
Magic Number: 0x2E736666
Version: 0001
Index Offset: 0
Index Length: 0
# of Reads: 320872
Header Length: 256
Key Length: 4
# of Flows: 220
Flowgram Code: 1
Flow Chars: TACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACG
TACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACG
TACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACG
TACG
Key Sequence: TCAG
>8UBVS:4:16
Read Header Len: 32
Name Length: 10
# of Bases: 141
Clip Qual Left: 5
Clip Qual Right: 83
Clip Adap Left: 0
Clip Adap Right: 116
Flowgram: 1.00 0.03 0.97 0.00 0.09 0.94 0.00 0.93 0.98 0.04 0.00 1.05 0.00 0.93 0.00 1.84 0.00 0.00 1.01
...
Flow Indexes: 1 3 6 8 9 12 14 16 16 19 20 20 22 22 22 23 25 26 29 29 32 33 35 36 39 41 43 46 48 48 50 50 50 51 51 54
...
Bases: tcagTGAGGCGGAAACTATTGTCGCTCAGGAAACCAGCGAATTTATGGCGTGGCTGCGAGCACAAAGC
CGCCAGCGAAACCAActtcgcggagtatacgccgagccagggaccagactggagactgaccaggacacagac
Quality Scores: 32 32 32 32 32 32 32 32 22 26 26 22 26 26 14 26 26 12 20 10 20 26 26 26 26 26 26 26 26 22 32 29 17 29 19 26
...

For sake of space, I cut short a few very long lines (indicated with ‘…’). A full explanation of the output is given in my blogpost.
Of note is the line that says

# of Flows: 220

Each flow is one base, in the order TACG, so there are 55 cycles of 4 flows. Note that the first instrument from 454, GS20, ran 42 cycles (168 flows), the GS FLX ran 100 (400 flows) and the current GS FLS Titanium runs 200 cycles (800 flows). The coming upgrade (when, Roche?) is supposed to give 400 cycles, 1600 flows…

I once looked into the number of bases one gets on average per cycle (of four flows), and this is remarkably close to 2.5. I basically chopped up E coli and human chromosomes into ‘cycles’ using flows in the order ‘TACG’ and counted how many T’s, A,’s C’s and G’s were actually read each cycle (I might write about this in another blog post). For the 454 instrument, the read length distributions, slightly depending on GC content, actually peaks around the then expected number of bases:
– the GS20, at 42 cycles, peaked around 100 bases
– GS FLX, 100 cycles, had a ~250 bases read length
– GS FLX Titanium, 200 cycles, peaks nicely at 500 nt.
Note that the expected output of the upcoming upgrade, supposedly at 400 cycles, actually seems to peak at 850 nt. Why this is I might have to look into once data become available…

Now, for the IonTorrent, at 220 flows and 55 cycles, one would expect a peak around 135.7. But the edgeBio data peaks at 107.5 bases!
If we look at the sff file data for the first read, it says

# of Bases: 141
Clip Qual Left: 5
Clip Qual Right: 83

So, the clipped read is 79 bases long (from position 5 to 83, including), but the raw read is indeed much longer. For this particular read, there is an adapter detected at the end, but for the majority of the reads in this file this is not the case.
Note that the reported sequence has the clipped part in lower case, and the remainder (trimmed part) in upper case.

So, I plotted this ‘raw read length’ versus the Trimmed one for DH10B-RUN01, see the graph.

The raw length peaks at 132 nt, close to, but below the expected 137.5 nt based on 55 cycles. The trimmed reads peak at 108. Also note the steep drop in read length beyond 108 bases. This drop is quite unusual, you never see it for 454 read, for example.

So, Ion Torrent reads are actually much longer than what is being reported. The key question, though, is whether there is any useful information in these extra bases. For this, I extracted the fasta sequences and quality scores, and uploaded them to the excellent tool PrinSeq, http://edwards.sdsu.edu/prinseq_beta. I first uploaded the data as they were with the original trimpoints:

The figure shows the quality value distribution along the read length in bins of 2 bases, steadily declining towards the end of the reads.

The same analysis for the retrimmed data:

Here the bins are 9 bases. It shows that the quality reaches its minimum at the 117 bp bin. This explains why the trimming is leading to the much shorter reads: beyond say 115 nt the quality of the bases is too low to be meaningful.

I would say that the analysis of EdgeBio’s file shows that the platform indeed has a potential to grow beyond the 100 base read length. ‘All it takes’ is increase the quality of the bases beyond cycle 27.

Part 2: Ion Torrent’s own data
Ion Torrent released their own run of E coli DH10B through their Torrent Dev community site http://lifetech-it.hosted.jivesoftware.com/index.jspa (check out the excellent analysis by Nick Loman on his blog http://pathogenomics.bham.ac.uk/blog/2011/05/first-look-at-ion-torrent-data-de-novo-assembly/
So, I naturally had a look at the information in this sff file. Here are the header and the first read, this time represented completely.

Common Header:
Magic Number: 0x2E736666
Version: 0001
Index Offset: 0
Index Length: 0
# of Reads: 522099
Header Length: 296
Key Length: 4
# of Flows: 260
Flowgram Code: 1
Flow Chars: TACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACG
TACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACG
TACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACG
TACG
Key Sequence: TCAG
>IB8EW:4:73
Read Header Len: 32
Name Length: 10
# of Bases: 190
Clip Qual Left: 5
Clip Qual Right: 94
Clip Adap Left: 0
Clip Adap Right: 0
Flowgram: 0.91 0.00 0.89 0.00 0.00 0.93 0.02 0.98 0.00 0.00 0.00 0.00 0.99 0.00 0.01 0.00 0.80 0.00 0.69 0.01 0.11 1.90 0.11 0.00 0.91 0.94 2.67 0.00 0.10 0.09 2.01 0.02 0.97 0.99 0.08 0.11 3.12 1.39 0.21 0.47 1.40 1.41 0.10 0.00 4.40 1.05 0.01 0.98 0.21 0.96 0.04 0.19 0.94 1.86 1.99 0.20 0.04 0.10 0.05 0.18 1.16 1.07 0.00 1.89 0.08 0.84 0.15 0.15 3.61 1.14 0.18 0.08 1.77 1.96 0.07 0.14 1.08 2.72 0.25 0.24 1.16 0.18 0.03 2.77 0.36 0.29 0.95 1.14 0.09 0.99 0.22 1.13 0.21 0.20 0.03 0.21 1.16 0.19 1.03 0.95 0.14 1.13 0.07 0.15 2.85 2.12 0.24 0.17 1.31 0.22 0.19 0.16 1.83 1.10 0.21 0.33 1.06 0.29 1.09 0.35 0.24 0.15 0.23 0.15 1.01 3.47 1.00 0.06 0.28 1.87 0.25 1.03 0.19 0.13 2.71 0.10 1.04 0.43 0.30 1.92 3.67 0.18 0.27 0.28 1.26 0.36 1.80 0.34 1.19 0.46 0.18 1.21 1.03 0.30 0.31 1.27 0.43 0.49 0.33 0.30 1.02 0.30 1.12 0.37 0.53 1.01 2.03 1.11 0.30 0.36 0.24 0.34 1.62 0.15 1.07 0.28 1.42 0.21 1.10 0.43 0.39 2.83 0.33 1.34 1.18 0.30 0.26 1.13 0.47 0.48 0.95 0.38 0.58 1.15 0.45 1.03 0.46 1.37 0.72 1.92 0.49 1.31 1.09 1.05 1.54 1.27 0.42 1.28 0.89 1.18 0.36 0.78 1.22 0.77 1.12 1.03 0.37 0.65 0.42 0.45 2.12 2.13 1.83 1.35 0.88 1.14 1.33 0.58 0.64 1.62 1.21 0.87 0.74 0.61 0.54 0.49 1.77 3.64 0.66 1.47 1.39 0.84 0.53 1.43 0.86 0.96 1.94 1.12 0.78 1.49 0.66 0.61 0.63 0.54 0.06 0.76 0.77 0.85 0.66 0.91
Flow Indexes: 1 3 6 8 13 17 19 22 22 25 26 27 27 27 31 31 33 34 37 37 37 38 41 42 45 45 45 45 46 48 50 53 54 54 55 55 61 62 64 64 66 69 69 69 69 70 73 73 74 74 77 78 78 78 81 84 84 84 87 88 90 92 97 99 100 102 105 105 105 106 106 109 113 113 114 117 119 125 126 126 126 127 130 130 132 135 135 135 137 140 140 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141 141
Bases: tcagATGCCTGTTTGGTATTTATCAAAAGACTCCGGCACCATTTTATTCCAGGGTAAAGAGATCGATT
TCCATTCTGCAAAGAAGCCCTGGAAAnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
Quality Scores: 28 23 23 23 22 16 16 16 15 23 10 14 14 6 9 12 23 23 23 23 6 10 9 10 10 10 10 5 10 20 20 20 20 14 20 14 20 20 20 14 10 10 10 10 5 10 10 10 12 14 12 12 12 7 12 17 17 17 17 17 20 20 20 20 20 17 17 17 9 17 12 8 12 14 15 15 15 8 8 8 3 8 15 12 14 14 14 7 14 14 8 10 10 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Note that the number of flows for this run is 260, equaling 65 cycles. But note the many lower case n’s at the end of the reported sequence! And the zeros in the reported quality values!
How to explain the difference between the sff files? I can only imagine the ‘masking’ in the file from the Ion Torrent website to be done deliberately on the part of Life Technologies. If so, were they trying to hide the low quality bases beyond the trimpoints?

Code used:
Read length distribution: I use a home-made perl script, called fasta_length.pl, which just plots the length of the fasta files it gets. The short version goes like this:

$/=">"; # set the record separator to the '>' symbol
; # remove the empty first 'sequence'
while (){
chomp; # remove the trailing '>' symbol
my @lines = split(/\n/,$_); # split the entry into individual lines based on the newline character
my $header = shift @lines; # the header is the first line (now without the '>' symbol)
my $seq = join "", @lines;
print length $seq; print "\n";
}

To get the distributions, I use the script in combination with sffinfo, awk and sort:

sffinfo -s DH10B-RUN01.sff |fasta_length | awk '{cnt[$0]++}END{for (x in cnt){print x"\t"cnt[x]}}' |sort -n -k 1,1

Plotting was done in excel (which I still use for these simple plots).
Getting the raw reads was done by making use of the trimpoint resetting function of the sfffile command. The first trimpoint is always 5, the first base beyond the four-base key sequence, for the second trimpoint I took the number reported as ‘# of bases’ (I later saw that I could just have set it to ‘0’, that would also have forced the max read length and would have made it much easier).

sffinfo DH10B-RUN01.sff|awk '{
if ($0 ~ />/){id=$0}
if ($0 ~/# of Bases/){l=$4}
if ($0 ~/Clip Qual Left/){print id" 5 "l}
}' |sed 's/>//' >DH10B-RUN01_trimpoints.txt

‘id’ becomes the read name, ‘l’ becomes the raw length.
To remake the sff file:

sfffile -o DH10B-RUN01_trimpoints_reset.sff -tr DH10B-RUN01_trimpoints.txt DH10B-RUN01.sff

To generate the fasta and qual files:

sffinfo -s DH10B-RUN01_trimpoints_reset.sff >DH10B-RUN01_trimpoints_reset.fna

sffinfo -sq DH10B-RUN01_trimpoints_reset.sff >DH10B-RUN01_trimpoints_reset.qual

I then manually (using the ‘nano’ program) removed read 8UBVS:1145:967, which was exceptionally long (and consisted of long stretches of TTTTTTTAAAAAACCCCCCGGGGGG), to tidy up the prinseq graph. Both fna and qual files were uploaded to prinseq after gzipping them.

Advertisements

2 thoughts on “There is more (length) to Ion Torrent reads than meets the eye (and is Ion Torrent hiding it?)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s