Coverage calculation w/genome

[ccard28]

Hello,
I am trying to calculate depth of coverage for some RNA-Seq data we have in lab. I have tried reading through posts but can't seem to get a grasp on what should be a fairly simple formula I think. I have paired-end reads (18,750,000 reads per file) at 100 bp and am using the bovine genome(UMD_3.1/btau6) which appears to be 2,670,422,299 for a total sequence length. Would the total coverage be (37,500,000 x 100)/2,670,422,299? This would give only like a 1.42x coverage which seems really low. Is this accurate or am I doing this calculation wrong? Also, what is an acceptable depth of coverage? This is also the amount of reads pre-alignment, for an accurate depth of coverage would I need to take (# reads aligned x 100)/2,670,422,299? Thank you in advance for help clarifying this issue.

[BAMseek]

Hi ccard28,

Since you are dealing with RNA-Seq data, you might want to look at the average coverage within the transcriptome. It's good to align to the genome, but you are probably interested in the average number of reads within the gene regions. Also, multiplying the number of aligned reads by the read length is a good back-of-the-envelope calculation, but it might not be quite right due to adapter trimming, quality trimming, soft clipping, reads that overhang your target regions . . . - all of which could cause the aligned bases of a read to be less than the read length.

Here would be my suggestion for calculating the total number of bases aligned to a target region:

samtools depth -b target.bed in.bam | awk '{s=s+$3};END{print s}'

So this would give you the sum of the aligned bases within your target (target.bed, which in your case would be the transcripts). You would then divide this by the total number of bases in the transcriptome (careful not to double count a base because it belongs to more than one transcript). I think samtools depth has a max depth of 8000, so that is one caveat when using it.

Justin