How To Find Number Of Mapped Reads For A Series Of Windows
2
2
Entering edit mode
9.0 years ago
roll ▴ 330

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

3
Entering edit mode
9.0 years ago

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.)

0
Entering edit mode

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

0
Entering edit mode

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.

1
Entering edit mode
9.0 years ago

Try bedtools makewindows:

Followed by coverageBed:

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

0
Entering edit mode

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?

0
Entering edit mode

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.

0
Entering edit mode

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?

0
Entering edit mode

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