number of .bam reads mapping to genes in gff format

Question

[Floydian_slip]

Hi,

I have a gff file containing a gene set and I also have a .bam file with read mappings to the genome. How can I get the number of reads that map to each gene (which as Start/Stop codon, CDS, intron specified in the gff file).

Any help will be greatly appreciated.

Thanks,
-n

Answer 1

[Devon Ryan]

You can use featureCounts, htseq-count, or Rsamtools with GenomicFeatures in R (you'd need to use the findOverlaps() function). featureCounts and htseq-count are probably the easier options.

[poisonAlien]

+1 for featureCounts speed.

Answer 2

[Bert Overduin]

You can use the intersect sub-command of bedtools with the -c option.