Getting mapping statistics from a .bed file

[ruchi1621]

Hello,
I have a .bed file that has the aligned sequence reads for some regions of interest. I want to find the number of reads aligned within the target region. Can anybody please help me on how that is to be done?
Thank you!

[Apexy]

Hi,
I have never stored short read alingments in a bed file. I have used the SAM format to store my reads, then use a bed file which have coordinates of those regions I want to extract reads from. So if you can get your alignment from a .bed to SAM, then this approach could be one way forward. The details of how to execute this protocol is found here

Extracting SAM entries mapping to a specific loci

Quite commonly, we want to extract all reads that map within a specific genomic loci.

First make a BED file, which in its simplest form is a tab-delimited file with 3 columns. Say you want to extract reads mapping within the genomic coordinates 250,000 to 260,000 on chromosome 1. Then your BED file (which I will call test.bed) will be:

chr1 250000 260000


Then run samtools as such:

samtools view -L test.bed test.bam | awk '$2 != 4 {print}'
OR
samtools view -S -L test.bed test.sam | awk '$2 != 4 {print}'
OR
samtools view -F 4 -S -L test.bed test.sam


Hope this helps

[crazyhottommy]

have a look at tabix http://samtools.sourceforge.net/tabix.shtml#4