Coverage calculation: long reads (RNA-seq)

Question

Say your aim is to calculate the coverage of an RNA-seq experiment generated with long-read sequencing (so, uneven read length).

Up to now, I relied on the Lander/Waterman equation:

$$C = L*N / G$$

where:

$C$ = final coverage
$G$ = haploid genome (or transcriptome) length
$L$ = read length
$N$ = number of reads

I have two major conceptual issues with this formula:

• $L$ is uneven in this case. Should we use the median/mean read length instead? If so, side question: if the RNA-seq experiment was executed following the size-fractionated protocol (PacBio), is it advisable to merge all the fractions prior to coverage calculation or should it be calculated separately for each fraction?
• Should $N$ consist of the mapped reads to the reference or simply the number of reads as found in the subreads whole dataset?

Answer 1

The read length is irrelevant when calculating the mean coverage statistic. It's simply the total number of bases sequenced divided by the target Xome length.

In the example provided in the question, $L*N$ is equivalently expressed as $\sum^{N} l_x$, or the total length of sequenced reads. For equal-length reads, this is easier to calculate as number of reads $*$ length of reads, but the full sum is necessary if read lengths are not equal.

For what it's worth, this value, $L*N$, is also equivalent to the average read length multiplied by the number of reads, i.e.

$$\frac{\sum^{N} l_x}{N} * N$$

... which makes it a bit more obvious why the total length of sequenced bases is very slightly easier to calculate.

In answer to the second part of the question, $N$ should ideally correspond to only the sequences that are from the reference. For an assembly, only the reads that were used to create the assembly should count when considering coverage. When there's no likelihood of contamination (and especially for RNASeq with variable coverage per transcript), a statistic using the total number of reads should be good enough for most purposes that I can think of. High-throughput sequencing is such a fuzzy process that a few tens of millions of bases here or there shouldn't make a difference.