Test for read coverage uniformity of genes?

[travelk]

Does anyone know of a program/package/method to test coverage uniformity of reads over individual genes?

My differential expression analysis shows 1500 genes... but when I look closer at IGV, I see that some of the genes that come up are artefacts in that all the reads of a particular sample map only to the UTR or within an intron. (I'm analysing single cell data so it's a bit messy that way.) Is there a way to test whether the reads cover all (or most) of the exons of a gene, rather than one by one in IGV?

Thanks for your help!

[blancha]

I see that some of the genes that come up are artefacts in that all the reads of a particular sample map only to the UTR or within an intron

Reads aligning to the introns will not be counted, so they will not affect the differential expression analysis.

Reads aligning to the UTR will be counted, obviously. However, you should be aware that some library preps will results in a very strong 3' bias, when using poly-A enrichment. We have uniform coverage at my institute even when doing poly-A enrichment. I do occasionally get data from other institutes, with a different library prep protocol, who do show this bias. If your genes show a systematic bias, you wouldn't want to exclude genes with a majority of reads aligning to the 3' UTR. If the majority of your genes have uniform coverage, you don't need to worry about this bias.

Is there a way to test whether the reads cover all (or most) of the exons of a gene, rather than one by one in IGV?

You could count the reads for individual exons with featureCounts. It will be difficult to have an algorithm to determine uniform coverage, though, in my opinion. Some exons are spliced out, and will therefore not have any reads align to them. There is also always some variation in the number of reads that align to exons.

Obviously, you should eliminate genes with low count reads, which is generally the most effective manner of removing artifacts.

If you still have artifacts after removing genes with low counts, I would rather investigate more the source of the artifacts than try and determine if the coverage is uniform.
For example, do the peaks correspond to repeat regions, or miRNAs for example?

[WhatsOEver]

I use RSeqQC to check for potential amplification biases in single cell data. The geneBody_coverage script is the one you are probably looking for. The nice thing is that the script is linked to R scripts in order to produce some easy-to-interpret graphics. The following code will calculate the coverage over a set of housekepping genes supplied with RSeqQC (But you are of course free to supply your own bed file):

python /usr/local/RSeQC-2.6.1/scripts/geneBody_coverage.py --refgene /usr/local/RSeQC-2.6.1/resources/hg19.HouseKeepingGenes.bed --input=/txtFileWithBamFilesToAnalyse --out-prefix ./genebodyCov_hg19housekeeping

[travelk]

Ok, I've tested it, and RSeQC is exactly the type of thing I was looking for. I need to modify it slightly to suit my purposes but it lets me automate the verification of my genes of interest.

Thanks for the suggestion!

[maxsalm]

Hi there! You may also find mRIN useful: http://www.ncbi.nlm.nih.gov/pubmed/26234653