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

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

alignment next-gen

• 11,078 views • link

updated 5.3 years ago by Ram

• written 9.9 years ago by h.rokny •

5 answers

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

• 0 views • link

written 9.9 years ago by igor

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

• 0 views • link

updated 5.3 years ago by Ram

• written 9.8 years ago by cmbkojic •

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.

• 0 views • link

written 9.9 years ago by swbarnes2

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)

• 0 views • link

written 9.9 years ago by benformatics

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.

• 0 views • link

written 9.9 years ago by EagleEye

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

• 0 views • link

written 9.9 years ago by glihm

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.

• 0 views • link

written 9.9 years ago by vakul.mohanty