Exon and junctions count per gene

Question

[rhonddaskl]

Hi,

I have aligned bam file from STAR. I want to count how many exon per gene, by using --quantMode GeneCounts, I can only get the count number per geneid, I want to get count number per geneid per exon_id. And I also want to get junction counts per gene. Can star do this kind of count for you or you need other count method?

Thanks.

Answer 1

[lieven.sterck]

Not sure if STAR can do this (might be), but other counting tools certainly can do this (HTseq-count, FeatureCounts, ...)

Key here is the GFF file you input, if it has the exon features present in it you can count on it (you might want to be a bit more lenient than usual as many reads will only match an exon partially (the splicing you mention indeed)

[rhonddaskl]

Thanks! I tried htseq-count and featurecounts, both can do exon-level count, but I haven't figure out how to do junction count, any suggestions? I'm using the "Comprehensive gene annotation" from genecode which contains exon information.

[lieven.sterck]

Great and yes, I kinda deliberately left that from my previous reply as I don't have an answer ready on that issue :-)

I'm pretty sure that is a count feature that will not be present as such in most tools, so I think it will require some post-processing of initially obtained counts.

I know that back in the days TopHat used to output a junctions file that contained all spliced mapped reads but as said I can't recall seeing any recent tool that outputs a junction count.

[lieven.sterck]

after some thought :

do you need those junctions counts on a per junction basis or on gene basis?

if the latter I was thinking you could try to do the count on exon basis for each gene, then the counts on gene basis for each gene (mRNA or CDS) and then subtract both numbers. That should give you a count of the junctions (as the "junction reads" won't be counted in the exon approach, but are in the gene approach, so the subtraction should give the junction counts)

[rhonddaskl]

thanks! it's really helpful. I need junction reads per gene basis, I'm gonna try to use TopHat.