How To Find Number Of Mapped Reads For A Series Of Windows

Question

[roll]

Hi,

I am trying to find the number of mapped reads in a given window. And to be more precise i want to find number of mapped reads for a series of windows. I divided the genome into 100bp and would like to know what are the number of mapped reads in these window. The final aim would be to put these into plink format (as phenotypes).

Is there any tool to do this?

thanks

Answer 1

[Alex Reynolds]

There's a usage example on the BEDOPS site that demonstrates how to use awk and bedmap to make bins (windows) and count reads across them.

The basic idea is:

$ bedmap --faster --echo --count windows.bed reads.bed > answer.bed

The file answer.bed is a BED file that contains coordinates of windows and the number (or count) of reads across those windows.

The usage script creates the windows, but you can easily replace that with your own sorted windows file. It also uses bam2bed to convert reads in BAM format to BED format, ready for use with bedmap.

(If your reads are paired-end, omit the --faster option.)

[roll]

Hi Thanks. I was trying the bam2bed feature yesterday and it took a lot of time to convert. Is this expected? I couldn't try out the counting the reads bit yet. The --count command, does it count the mapped reads or all reads? thanks

[Sean Davis]

Yes, the conversion can take a lot of time (enough time to read the file and write another one, at least). The count command counts reads overlapping with the windows.bed windows; unmapped reads do not overlap any windows, so they are not counted.

Answer 2

[Sean Davis]

Try bedtools makewindows:

http://bedtools.readthedocs.org/en/latest/content/tools/makewindows.html

Followed by coverageBed:

http://bedtools.readthedocs.org/en/latest/content/tools/coverage.html

There are many other alternative approaches using different software packages (GenomicRanges, bedops, ....)

[roll]

Hi Sean Davis, thanks. I currently only have bam files. So if i get it right i need to convert them to bed? And then use the coverageBed to extract the number of mapped reads?

[Sean Davis]

Makewindows makes a bed file with windows across the genome of the size you want. The, coverageBed takes this bed file and your bam file to do the counting. No, you do not need to convert your bam files to bed files.

[roll]

Can you please give an example command line how will it use bam and bed file to get the number of mapped reads? Can i do this for 50 bam files? or does it have to be one by one?

[Sean Davis]

If you have multiple bam files, you might also look at:

http://bedtools.readthedocs.org/en/latest/content/tools/multicov.html