How Many Reads In A Bam File Are Aligned To a Specific Region

Question

[h.rokny]

I have BAM files from Hi-C experiment. I'd like to know how many of my aligned reads fall in a particular region of a given chromosome (i have a bed file that contains list of regions).

bed file:

chr1 1000 2000  
chr1 3000 4000   
chrY 1000000 2000000

Thanks

Answer 1

[igor]

You can use bedtools coverage -a regions.bed -b reads.bam

[cmbkojic]

I was trying to solve the same problem with Bedtools, version 2.26.0.

I am using Bedtools coverage option to calculate coverage of my aligned reads in .bam format over 2kb windows of chromosome. So, I created a bed file for the chromosome separated in 2kb windows. The bed file has only 3 columns:chr1 0 2000; chr1 2000 4000... and so on. I am running the bedtools coverage option with -a genome.bed -b reads.bam. When I try to upload it to UCSC it tells me ''Error line 1 of custom track: Expecting + or - in strand'' How can I solve this?

The output of bedtools is like this:

chr1 0 2000 0 0 2000 0.0000000
chr1 2000 4000 0 0 2000 0.0000000

Answer 2

[swbarnes2]

samtools view -L can take a .bed file as input, and output reads that overlap from the .bam. samtools view -L -c will just return the number of reads.

Answer 3

[benformatics]

An R code version using GenomicRanges. Not super fast but will do the job and let you keep working with the information in R.

library(GenomicAlignments)
library(rtracklayer)

## import bam file
bamfile <- readGAlignments('sample.bam') ## assuming single end otherwise use readGAlignmentPairs

## import bed file with regions of interest
bedfile <- import.bed('regions.bed')

## count number of overlaps for each element of the bedfile object
overlap.counts <- countOverlaps(bedfile,bamfile)

Answer 4

[EagleEye]

featureCounts is very simple to use and it is fast,

http://bioinf.wehi.edu.au/featureCounts/

Make your bed file look like SAF file as mentioned in above link.

[glihm]

If you can make a Gff with your regions of interest (pretty simple and fast), I think like EagleEye => Faster solution is using FeatureCounts.

Answer 5

[vakul.mohanty]

pysam's (python module)fetch can be used to extract and count reads(or do just about anything) for a specific region in a BAM file. I find this a very flexible way to parsing BAM files.