Hi
How can I get a read count per exon per transcript in a particular gene from BAM file? I tried using htseq-count with id_attr=transcript_id but this approach does not give me the read count per exon per transcript. I tried using DEXSeq to get the read count, but it does not get what I want. Is there any way to get the read count? I only got one gene of interest to analyse which TPM1. This gene has more than 40 isoforms, so it is possible to get each of the isoforms their read count for each exon. Thanks!
•
written
by
wan.fahmi •
There are many choices and tools for that. How about using coveragebed from bedtools?
bedtools coveragebed -s -abam bamFile -b exons.bed > yourResultFile
-s: only if you have stranded RNA-seq data.
To compare the expression level of exons, divide by tag-per-million and also normalize by length (to take varying exon length into account)
•
written
by
Chirag Nepal
I tried but it need 2 bam files which is not the case. I don't want to compare, I just want a read count of each exon for each transcript. I think it quite arbitrary as it quite hard to distinguished which read mapped to which isoform. How do you think?
Paste your exact command, so we can probably see where you might have wrong.
bedtools multicov
Your first step is to get a bed file in place which defines coordinates for each exon of each transcript.
If you do not have that file, consider using DEXSeq. It has python scripts to modify standard gtf file to make exon level counts.
what is the difference between multicov and coveragebd ? Multicov can take couple of bam files, right ?
bedtools multicov has options for selecting, % of overlapping bases to make a count, minimum mapping quality, considering duplicate reads etc and it gives you a matrix which can be fed in to DE analysis packages etc.
•
written
by
GouthamAtla
Exactly. That is the difference I knew before.
•
written
by
Chirag Nepal
I tried it too but please see my above comment.
•
written
by
wan.fahmi •
Log in to answer this question.
