[pmiguel]

Here is the problem:



Green is "Moore's Law", or cost of a unit of computer power, as projected by a halving in cost every 18 months. Blue is cost per megabase of raw sequence from:

Wetterstrand KA. DNA Sequencing Costs: Data from the NHGRI Large-Scale Genome Sequencing Program Available at: www.genome.gov/sequencingcosts. Accessed 9/7/2011

Whether deliberately or not, one will choose a percentage of a budget for analysis of data collected by an experiment. If you are not carefully watching the dramatic fall in the price of sequencing, you might be stuck in a 2006 mind-set. That is, back in 2006 a mythical unit of computation would be 8X its current price. Whereas a million bases of sequence would cost 3000X more.

--
Phillip

[rskr]

Runtime complexity? I love it when people trot out these raw numbers of gigabase vs. compute time charts. They never seem to mention that as the quantity of the data grows certain algorithms don't scale linearly. Take for example short read mapping. The longer the reads are the more mismatches are required to get accurate alignments, but the more mismatches per read and the more bases per read the more complex the search is. So having a gigbase of 50 mers is not the same as a gigabase of 100mers, which is why I am guessing so many labs still use 36mer data even though it is horrendously inaccurate for mapping or anything else, but they don't have the compute resources to deal with a petabyte of 100mer data, even though it is possible to routinely produce the 100mer data, at a cheaper price per base than 36mer data.

[HESmith]

Hi rskr,

I agree with your general point, but there are a number of applications (e.g., RNA-Seq or ChIP-Seq in model organisms) for which 36bp reads are perfectly acceptable rather than "horrendously inaccurate" for mapping. Not every project requires a petabyte of data (but I want lots of cores and RAM for the ones that do!).

[rskr]

Hi rskr,

I agree with your general point, but there are a number of applications (e.g., RNA-Seq or ChIP-Seq in model organisms) for which 36bp reads are perfectly acceptable rather than "horrendously inaccurate" for mapping. Not every project requires a petabyte of data (but I want lots of cores and RAM for the ones that do!).

Well, if you can settle for mapping 70% of what you are studying uniquely, 36mers will only guarantee something in that range. Sure the genomes are more than 70% unique, but when you factor in error and mapability(IE things that are nearly unique), forget studying indels forget studying splice variants, don't trim adapters etc.

Might as well throw out all the insignificant Illumina studies done before 2010, because that is what you call cheap science.

[pmiguel]

Runtime complexity?

Is this your catch-phrase or something? You are invoking it here without any context I can discern.

I love it when people trot out these raw numbers of gigabase vs. compute time charts. They never seem to mention that as the quantity of the data grows certain algorithms don't scale linearly. Take for example short read mapping. The longer the reads are the more mismatches are required to get accurate alignments, but the more mismatches per read and the more bases per read the more complex the search is. So having a gigbase of 50 mers is not the same as a gigabase of 100mers, which is why I am guessing so many labs still use 36mer data even though it is horrendously inaccurate for mapping or anything else, but they don't have the compute resources to deal with a petabyte of 100mer data, even though it is possible to routinely produce the 100mer data, at a cheaper price per base than 36mer data.

Not sure what you mean. Are you saying the amount of computer power you need is even higher than one would naively presume from the graph above?

If so, then by your argument it might be more cost-effective to use the 36mers even though the 100mers map more accurately. As the amount of data per dollar mounts, the budget going towards computer cycles/storage approaches 100%. Controlling costs then would nearly entirely be a matter of managing computational resources.

--
Phillip

[rskr]

Is this your catch-phrase or something? You are invoking it here without any context I can discern.

Hook line and sinker.

Not sure what you mean. Are you saying the amount of computer power you need is even higher than one would naively presume from the graph above?

Yes, depending on what one is doing the algorithms will probably scale super linearly in the amount of data. For example just sorting the reads is super-linear N*log(N) in the number of reads, and usually the time complexity in read length is greater than the time complexity in the read number O(l)>O(N), for the non computational that is big Oh length is greater than big OH read count.

If so, then by your argument it might be more cost-effective to use the 36mers even though the 100mers map more accurately. As the amount of data per dollar mounts, the budget going towards computer cycles/storage approaches 100%. Controlling costs then would nearly entirely be a matter of managing computational resources.

My point is that cheaper in the short read sense is also lower quality, in fact so low you might not even bother, though their are many a PhD with a career vested in 36 mers. maybe because their perl doesn't scale, other than that I can't figure out why they still do sequence 36mers, because they are more expensive per base than 100mers, and don't cover a fraction of the genomes or transcriptomes, find splicing, indels or more than 3SNPs in a row.

[HESmith]

My previous post must have been too terse to understand; here's some data to clarify (C. elegans ChIP-Seq, identical library, two lanes each indicated read length, percent uniquely aligned):

Single-end 36bp: 85.90%, 85.92%
Single-end 76bp: 86.22%, 86.58%
Paired-end 36bp: 86.41%, 86.69%

In this example, doubling the read length increased the number of mappable reads by (wait for it) - less than 1%. Since all I care about in this experiment is the position of the read, the per base price is not the relevant metric - it's the marginal increase in the number of mappable reads. At current reagent costs for the HiSeq, I just spent $170 to map an additional 700K reads. Not much bang for the buck. Plus, that doesn't factor in the lost opportunity cost of running the sequencer for an extra two days that could be better spent on sequencing other samples.

I stand by my previous claim; for some applications, 36mers are fine.

[pmiguel]

Yes, depending on what one is doing the algorithms will probably scale super linearly in the amount of data. For example just sorting the reads is super-linear N*log(N) in the number of reads, and usually the time complexity in read length is greater than the time complexity in the read number O(l)>O(N), for the non computational that is big Oh length is greater than big OH read count.

I guess when I posted the graph up-thread the reaction I expected was "the vastly increasing amount of sequence data is going to bankrupt me computationally."

But your response, if I am now parsing your correctly, is "It is even worse than the graph implies. Because longer read lengths are desirable, the computational load engendered by those vastly increasing numbers of bases is greater than it would have been had only the number of reads been increasing."

Sure. Could be. I mean, I don't insist on being the most pessimistic person in the room...

I'm just saying that the OP's boss, who set aside less than generous amount of the budget aside for informatics, is not necessarily the monster he or she is being painted in the thread. I think in most domains of biology it is not common that collecting data will be a tiny percentage of the cost of analyzing it.

--
Phillip

[pmiguel]

[...]
I stand by my previous claim; for some applications, 36mers are fine.

Yes, and it stands to reason. 36 nt offers 2^72 possible combinations. You lose some of those to sequencing errors. But usually the reason a 36 nt can't be uniquely mapped (after quality filtering) will be duplication. And, in most cases duplicated segments will be larger than 100 bases, so longer reads will not necessarily save the day.

--
Phillip

[rskr]

Yes, and it stands to reason. 36 nt offers 2^72 possible combinations. You lose some of those to sequencing errors. But usually the reason a 36 nt can't be uniquely mapped (after quality filtering) will be duplication. And, in most cases duplicated segments will be larger than 100 bases, so longer reads will not necessarily save the day.

--
Phillip

However it you look at the quality of the mappings, indels, multiple, SNPs, splicing etc. Will fail so if you are doing RNA-seq on a non-homozygous species or on anything other than a homozygous species that was used to sequence the genome, you are going to have problems dealing with the variation, for anything but the simplest of SNPs, and even then, you are going to have problems with sensitivity, because mapping degrades near the end of the reads, so be prepared to throw away 12 bases on both ends of the reads. Or better yet, get longer sequences.

[HESmith]

However it you look at the quality of the mappings, indels, multiple, SNPs, splicing etc. Will fail so if you are doing RNA-seq on a non-homozygous species or on anything other than a homozygous species that was used to sequence the genome, you are going to have problems dealing with the variation, for anything but the simplest of SNPs, and even then, you are going to have problems with sensitivity, because mapping degrades near the end of the reads, so be prepared to throw away 12 bases on both ends of the reads. Or better yet, get longer sequences.

That's funny; I don't recall recommending 36mers for any of the applications/situations that you mention...

[rskr]

That's funny; I don't recall recommending 36mers for any of the applications/situations that you mention...

You did report mapping statistics, yet failed to report bench marking statistics, or address my concerns about the quality or your alignments, so say whatever you want.

[frozenlyse]

rskr - your argument boils down to "well 36bp is too short for MY application(s) so **** them right off"

cool story bro

[rskr]

rskr - your argument boils down to "well 36bp is too short for MY application(s) so **** them right off"

cool story bro

No I believe if you read into my argument a little more it boils down to that 36bp are most like too short for YOUR application therefore if you use 36bp I am going to write you off Sure there are applications for 36mers out there maybe siRNA or something.

[pmiguel]

[...] and even then, you are going to have problems with sensitivity, because mapping degrades near the end of the reads, so be prepared to throw away 12 bases on both ends of the reads. Or better yet, get longer sequences.

Does it? On 100 nt reads it is possible that the sequence quality dips low enough that the last 12 bases are poor. But I would not expect this to be an issue for 36 nt reads.

--
Phillip