This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How To Find Number Of Mapped Reads For A Series Of Windows

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

reads plink windows
• 6,061 views • link
updated by Alex Reynolds Alex Reynolds • written by roll roll

2 answers

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 views • link
written by Alex Reynolds Alex Reynolds

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 views • link
written by roll roll

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.

• 0 views • link
written by Sean Davis 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, ....)

• 108 views • link
written by Sean Davis Sean Davis

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 views • link
written by roll roll

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 views • link
written by Sean Davis Sean Davis

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 views • link
written by roll roll

Log in to answer this question.

More posts like this