read coverage and transcripts

Question

[qudrat]

Can somebody suggest me how to obtain read coverage of a particular transcript in a transcript assembly done by TopHat/Cufflinks?

[WouterDeCoster]

a transcript assembly done by TopHat/Cufflinks?

Correct me if I'm wrong, but are you sure this is an assembly? TopHat is an aligner, not an assembler.

[qudrat]

@WouterDeCoster, TopHat is an aligner and Cufflinks is used for assembly.

[WouterDeCoster]

Right, but it doesn't modify the alignment, does it? You want to assess the coverage in the alignment. An assembly doesn't have a coverage.

[qudrat]

Yes I want to assess the coverage in alignment. Actually in my question, I wanted to make it clear that what software I used

[WouterDeCoster]

Alright very well. Excuse me for my pedantic nitpicking here, but correct terminology is quite important for quickly getting the right answers.

Answer 1

[Renesh]

Here, you need to find the number of reads that are mapped to the given transcripts. For this purpose you need mapped BAM file, the transcript chromosome number and it's start-end co-ordinates. If you have this information, you can find read coverage using samtools as follows

samtools view BAM_file Chromosome:start-end

It will give you all reads that mapped to given transcript region. To count the read coverage, you can pipe wc -l command

samtools view BAM_file Chromosome:start-end | wc -l

[qudrat]

Renesh, what if a gene has more than one transcript and each transcript has same co-ordinate as in case of exon skipping or intron retention. How one would find the coverage for each different transcript?

Answer 2

[Hussain Ather]

Alongside samtools, bedtools offers ways of reading Tophat/cufflinks output too http://bedtools.readthedocs.io/en/latest/index.html

Answer 3

[WouterDeCoster]

In addition to the tools mentioned earlier, commonly used tools for counting reads in a genomic interval (gene, exon,...) are featureCounts and htseq-counts. I would recommend featureCounts, because it's very fast, has convenient options, but htseq-counts also works fine. Both are nicely documented.