How would you calculate the coverage for RNAseq?

[LeonDK]

I am thinking that the most important quantification parameter must be how many reads are mapped within annotated regions?

Hence it is not mapped reads per whole genome, but rather mapped reads per annotated region, like so:
coverage = n_mapped_reads * read_length / (n_annotated_regions * average_length_of_annotated_regions)

This way we get average number of times each annotated nucleotide was sequenced?

L

[GenoMax]

So if a region is not annotated you are going to consider reads mapping there not important/real?

[LeonDK]

So if a region is not annotated you are going to consider reads mapping there not important/real?

Hmm... No and yes... Of course if something maps it (in theory) means that a transcript originated from there. But if I am to use the annotations for GRCh38 for analysis of DEG, than the reads mapping in non-annotated regions are not considered?

[GenoMax]

If you are explicitly interested in the "annotated" regions then that is a limitation you are imposing. Reads mapping elsewhere will still be real.

Take a look at these: http://bedtools.readthedocs.org/en/l.../coverage.html http://bedtools.readthedocs.org/en/l...genomecov.html

[jwfoley]

Why do you want to calculate coverage in the first place? If you're using it for gene-expression profiling, RNA-seq is a read-counting application, not a base-counting application, so for most purposes the read count is a more appropriate measure of sequencing depth.

[LeonDK]

Why do you want to calculate coverage in the first place? If you're using it for gene-expression profiling, RNA-seq is a read-counting application, not a base-counting application, so for most purposes the read count is a more appropriate measure of sequencing depth.

Good point!

So in your opinion, it is sufficient to report the number of mapped reads? E.g. "we included samples with at least 10 mio. mapped reads" or something similar?

[jwfoley]

Good point!

So in your opinion, it is sufficient to report the number of mapped reads? E.g. "we included samples with at least 10 mio. mapped reads" or something similar?

Yes, that is the single number that wraps it up. For reference, ENCODE's old standard for gene-expression profiling was 30 million paired-end reads (for human samples), but that's not the greatest recommendation anyway.

[GenoMax]

So in your opinion, it is sufficient to report the number of mapped reads? E.g. "we included samples with at least 10 mio. mapped reads" or something similar?

What kind of "report" are we referring to here? How do you decide on a number for cut-off? Is the distribution of mapped reads assured to be equally representative for all samples, if you are only going to report a total number for mapped reads? And so on ...

[LeonDK]

Ok, so of course I know very well, that there is no magic number. In fact pretty much every cut-off we use is arbitrary, especially the infamous 0.05 - Anyhoo...

For publication-purposes it is necessary to account for how 'bad' samples were excluded.

Hart et al. (10.1089/cmb.2012.0283) states that: "In other words, a sequencing depth of 10 million reads will ensure that approximately 90% of all genes will be covered by at least 10 reads."

That's why a suggested 10 mio. as a 'nice' arbitrary number