count reads or count base-pairs

[yuelics]

Hi all,

I am new to the field. Here is a silly question:

For quantification of defined genomic regions (let it be genes, windows, etc) based on sequencing data, why do people rarely count base-pairs (a convenient function for doing that is "coverage" in R) but rather count reads?

Thanks in advance,
Yue

[westerman]

By "quantification" I presume you mean something like a transcriptome project where you want to find out the expression level of genes. For such a project every biological transcript should equal to exactly one read. E.g., if a cell is producing 50 mRNAs for protein 'A' and 10 mRNAs for protein 'B' then there will be a total of 60 reads generated, 50 of which will map to gene 'A' and 10 to gene 'B'. Depending on the sequencing technology, the reads may vary in length. Therefore the number of reads is a better count of the transcriptomes than the base pair coverage.

For other projects you may want (or at least get by) with a base-pair count. For example when discovering SNPs via whole-genome-mapping or exon-mapping then, since they are (duh!) single-base the number of reads will equal the number of base-pairs. So either number could suffice. However for SNP quality control, it is useful to know the number of *unique* reads covering the SNP thus one might as well use the read count instead of the base-pair coverage count. (BTW: the number of unique reads should roughly equal the number of reads otherwise you potentially have problems.)

For other projects, such as general mapping or de-novo assembly then the coverage via base-pairs is an important statistic.

Hope that the above helps.

[yuelics]

Thanks westerman. Based on your analysis, if we know the transcripts along the genome, then it is better to use read coverage.

In de novo approach, where a sliding window of fixed size is used, should the base coverage be equivalent to window size of 1?

I have never seen anyone set this window size to 1 (25 bp at least in some literature, but usually 1kb). But using window size of 1 eliminates the concern for reads spanning across two windows. But the complication may be that we introduce highly dependent event for bases within a read, which might make it hard to model with simple Markov chain of limited order (e.g., HMM of order 1).

If using read coverage, what would be an effective window width, and should the adjacent windows of large size (e.g., 1kb) be overlapping each other to some extent (e.g., 500 bp)?

Thanks for the helps!

[westerman]

If you are talking about de-novo transcriptomes (which I do a lot of, mainly from the 454 platform) here too I would use read counts. Where appropriate I bundle all of the reads together to create the de-novo transcripts. I then figure how how many of each sample's reads map to the transcripts. This is the count that we use. In many ways this work is no different than when you actually know the transcripts ahead of time ... except for de-novo you have no idea (without extra work) of what the transcripts actually represent or if they actually exist in nature.

I can not see using a sliding window on transcripts. Except, perhaps, to determine if the transcripts are likely to exist. In other words transcripts that have a large number of bases supporting them are much more likely to be correct than those without a number of bases. The sliding window count of bases would even out disparities in read coverage and could allow one to differentiate between false transcripts versus those just expressed at a very low level.


But then you may be doing a more in-depth analysis (or an entirely different analysis) than what I normally do. Basically what we give our customers is a list of their generated transcripts (plus blast/GO terms) with the associated read counts for each of their samples. It is up to them to figure out if those transcripts and counts signify anything biologically important or how to model them further.