Genome size in X coverage calculation for RNASeq experiments

[LauraGP]

Hello everyone,

Traditionally in sequencing, X coverage is determined as the product of the number of sequenced reads and their length, divided by the genome size.

Regarding RNASeq experiments and using purified poly(A)+ RNA as starting material for library preparation, wouldn't be more accurate to use the product of the number of known genes (as in Ensembl) and tha average transcript size for the given organism, instead of the whole genome size, as the denominator in X coverage calculation?

At the end, because we are working with an mRNA-enriched library for sequencing it is obvious that there will be regions in the genome which we will not cover at all. So maybe we could just exclude them from our X coverage calculation?

Thanks in advanced for any help!

[kmcarr]

Hello everyone,

Traditionally in sequencing, X coverage is determined as the product of the number of sequenced reads and their length, divided by the genome size.

Regarding RNASeq experiments and using purified poly(A)+ RNA as starting material for library preparation, wouldn't be more accurate to use the product of the number of known genes (as in Ensembl) and tha average transcript size for the given organism, instead of the whole genome size, as the denominator in X coverage calculation?

At the end, because we are working with an mRNA-enriched library for sequencing it is obvious that there will be regions in the genome which we will not cover at all. So maybe we could just exclude them from our X coverage calculation?

Thanks in advanced for any help!

In one sense that is more accurate, but calculating coverage like this for RNA-Seq experiments is meaningless since it doesn't take into account the wide variation in transcript abundance. Within a single data set for an mRNA library you may find coverage of > 1000X for highly expressed transcripts and < 1X for rare transcripts.

If you are performing a differential expression study then most researchers think in terms of the number of reads (counts) per mRNA library, regardless of read length or total transcriptome length. On the other hand if your goal is de novo transcriptome assembly/discovery then you do need to think in terms of total Gbp of raw sequence data, but you acknowledge that you will sequence the abundant transcripts to a ridiculous depth in order to get enough data so that rarer transcripts gain a reasonable level of coverage for assembly.

[Brian Bushnell]

You can calculate coverage from total number of bases sequenced divided by the total transcriptome size. But that's not very useful in organisms (such as eukaryotes) with differential splicing. I don't think it is useful to use calculate average coverage for RNA-seq in a eukaryote. It makes more sense to decide whether you have an adequate number of reads to statistically determine expression.

[LauraGP]

Hello,

Thanks for your replies!

I believe we are dealing at the end with the dilemma: to increase the number of sequenced reads, that is the sequencing coverage/depth, or to increase the number of biological replicates. Ideally, we would like to increase both, but the budget is not always on our side.

For differential gene expression analysis, our goal is to statistically determine if an observed difference in read counts for each gene, is due to a particular condition or to random variation. In DESeq2 package, the use of shrunken log fold change estimates allows to deal with genes with low read count, basically filtering those which are unlikely to be found after a user-defined threshold for FDR. It seems to me that the identified DE genes after comparison, are just the top of the iceberg as soon as the sequencing depth or the number of sequenced reads are lower.

My main question is then, how to accurately estimate the number of reads we need to sequence per sample, in order to identify at least the 95% of the genes which are altered by our experimental condition?
So far I have just found in the literature previous experimental evidences and recomendations, but not a formal description of how to perform this kind of estimation.
Also, how to determine the variance in sequencing depth (X coverage) across the genome/transcriptome, if not by ultra-deep sequencing approaches? Is there a way to estimate a priori the variance in coverage, for instance analyzing the GC content across the known reference genome?

[nucacidhunter]

As pointed earlier for transcriptome coverage is not applicable as the size of transcriptome is variable and not known. If you have many samples (treatments, replicates) then you can sequence some of libraries deeper and down sample to find the optimum read number. For instance, for human DE 25M read is recommended but you can sequence some to 40M reads and if increasing read number gives better results for your experiment then do a top up sequencing for the rest.

[LauraGP]

As pointed earlier for transcriptome coverage is not applicable as the size of transcriptome is variable and not known. If you have many samples (treatments, replicates) then you can sequence some of libraries deeper and down sample to find the optimum read number. For instance, for human DE 25M read is recommended but you can sequence some to 40M reads and if increasing read number gives better results for your experiment then do a top up sequencing for the rest.

Thanks for your suggestions! got it.