samtools view region vs. bamtools filter region
1
1
Entering edit mode
7.5 years ago
Lídia ▴ 80

When extracting reads aligned to a certain region, what is the difference between using samtools view and bamtools filter region?

Here is the code and the counts of the number of reads recovered. I don't understand why the number of reads recovered are so different between the two approaches. Any help??

samtools index IV75_realigned_reads.bam
samtools view -hb IV75_realigned_reads.bam 3R:15,000,000-15,100,000 > Control_region.bam
bamtools count -in Control_region.bam
# 28034
bamtools filter -region 3R:15,000,000-15,100,000 -in IV75_realigned_reads.bam -out test.bam
bamtools count -in test.bam
# 7523775
samtools NGS bam bamtools • 5.8k views
ADD COMMENT
0
Entering edit mode

It is not clear what is more maddening

  1. that bamtools uses a different character to indicate the range or
  2. that it still produces a count when you provide the incorrect one
ADD REPLY
0
Entering edit mode

The addition of both terms + the interaction... ;-)

ADD REPLY
2
Entering edit mode
7.5 years ago

as far as I can read in the manual of bamtools, region are defined with two dots and not using a hyphen. You should try (and don't play with fire, remove the commas):

3R:15000000..15100000

REGION string format: A proper REGION string can be formatted like any of the following examples (where 'chr1' is the name of a reference (not its ID)and '' is any valid integer position within that reference.):
-region chr1
only alignments on (entire) reference 'chr1'
-region chr1:500
only alignments overlapping the region starting at chr1:500 and continuing to the end of chr1
-region chr1:500..1000
only alignments overlapping the region starting at chr1:500 and continuing to chr1:1000
-region chr1:500..chr3:750
only alignments overlapping the region starting at chr1:500 and continuing to chr3:750. T

ADD COMMENT
0
Entering edit mode

Sorry for such a silly mistake... You were completely right!

filter -region 3R:15000000..15100000 -in IV75_realigned_reads.bam -out test.bam
bamtools count -in test.bam
28034

Thank you very much for your answer

ADD REPLY

Login before adding your answer.

Traffic: 1803 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6