How to get the number of reads for a set of coordinates from a bam file.

Question

[hins.2026]

I have Whole exome data and aligned them to the genome, so I also have bam files. I also have the corresponding chromosome coordinates for which I need to find the number of reads. So, the question is that how can I count the number of reads that map to each coordinate separately in the genome. in other word, I want to count the number of reads that map to the list of coordinates which I have. do you guys know how to do that? Thanks in advance

[Juke34]

Using Python and the pysam lib you can do

samfile = pysam.AlignmentFile(bam_file, "rb")  
samfile.fetch(chromosome_name, start, stop)

To get all reads mapping at coordinate of the chosen chromosome.

or using samtools you convert it in sam and then extract the corresponding line with awk. Otherwise I'm sure tools like samtools, BEDOPS or Subread package have command dedicated for that purpose. Have a look here: Extract Reads From A Bam File That Fall Within A Given Region

Answer 1

[ATpoint]

Please use the search function and google, this has been asked dozens of times, e.g:

To Calculate The Exact Total Number Of Mapped Reads In Exome Regions