Read count per exon per transcript
2
0
Entering edit mode
6.7 years ago
wan.fahmi • 0

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=transcipt_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!

read count RNA-Seq bam file dexseq • 4.3k views
0
Entering edit mode
6.7 years ago
Chirag Nepal ★ 2.3k

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)

0
Entering edit mode

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?

0
Entering edit mode

Paste your exact command, so we can probably see where you might have wrong.

0
Entering edit mode
6.7 years ago

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.

0
Entering edit mode

what is the difference between multicov and coveragebd ? Multicov can take couple of bam files, right ?

0
Entering edit mode

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.

0
Entering edit mode

Exactly. That is the difference i knew before.

0
Entering edit mode

I tried it too but please see my above comment.