Extract Mapped Reads grouped by positions

[JdeBruin]

Hi,

I have mapped reads to a reference genome using Bowtie2 and BWA, using samtools I have created bam files etc.

What I would like to achieve is to group reads based on the mapping location, a new group is formed when a large gap in the mapping results are observed.
For example, all reads that map to position 1-10000 will be in group 1 then a gap in the mapping may be present and the next mapped read will start at 15000 - 18000 followed by a gap this will then form group2.

Any help would be appreciated.

[sphil]

Do you know where your "gaps" are present?

[JdeBruin]

Hi thanks for the reply.

The gaps can be seen in the mapping results using IGV/picard thus the gaps can be defined.

[sphil]

samtools view is a method of choice, if you have defined gaps you can go like this:

samtools view indexed_sorted.bam chr2:1,000,000-2,000,000

for example to view just the alignments on chr2 between region 1000000 and 2000000. Note the coordinate is 1-based.
If you need more infos this page might help: samtools FAQ
Every read which is returned for your specific position is in that group. Region to look at == group.

Hope that helps.

[swbarnes2]

You can also use BEDTools intersectBed to make a new .bam which is the intersect of your .bam file and a .bed file of coordinates.

[JdeBruin]

Thanks for all the solutions!!

[JdeBruin]

Hi,

when I use the following:
samtools view indexed_sorted.bam chr2:1,000,000-2,000,000
I get the following error:
[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
[main_samview] fail to read the header from "indexed_sorted.bam".

Would the view not require a flag as to specify a bam file?

[JdeBruin]

I have been able to solve the problem