[Bioperl-l] Bio::DB::Sam strange behaviour for read pairs
Frank Schwach
fs5 at sanger.ac.uk
Wed Jan 6 22:16:13 UTC 2010
I'm trying to extract paired reads from a BAM file that span a given
region. I would then like to get the two read ends of the sequenced
clone that spans the region.
I use Bio::DB::Sam->get_features_by_location for this and it does give
me the correct read pairs as a region match but it doesn't give me both
read pairs in all cases.
Here is the script:
#!/usr/bin/perl
use Bio::DB::Sam;
my $usage = "usage: $0 BAMFILE CHROMOSOME STARTPOS ENDPOS\n" ;
my ($bam_file,$chrom,$start,$end) = @ARGV ;
die $usage unless $bam_file && $chrom && $start && $end;
my $bam = Bio::DB::Sam->new(-bam => $bam_file);
my @pairs = $bam->get_features_by_location(
-type => 'read_pair',
-seq_id => $chrom,
-start => $start,
-end => $end);
print "region: $chrom:$start..$end\n" ;
foreach my $pair (@pairs) {
print " pair: id: ".$pair->id.", start".$pair->start.',
end:'.$pair->end."\n";
my ($first_mate,$second_mate) = $pair->get_SeqFeatures;
print " first_mate: start:".$first_mate->start.',
end:'.$first_mate->end."\n";
if ($second_mate){
print " second_mate: start:".$second_mate->start.',
end:'.$second_mate->end."\n";
} else {
print " no second mate\n";
}
}
And here are the matching pairs that it produces with one of my files
for the region tal12:22479..29232:
region:
tal12:22479..29232
pair: id: tal-2446c08, start17496,
end:29423
first_mate: start:28540,
end:29423
no second
mate
pair: id: tal-2463d10, start23534,
end:31363
first_mate: start:23534,
end:24448
no second
mate
pair: id: tal-2371c09, start20860,
end:28230
first_mate: start:27604,
end:28230
no second
mate
pair: id: tal-2440b06, start19232,
end:27099
first_mate: start:26025,
end:27099
no second
mate
pair: id: tal-2327g09, start18909,
end:26129
first_mate: start:25354,
end:26129
no second
mate
pair: id: tal-2381b05, start25658,
end:35054
first_mate: start:25658,
end:26295
no second
mate
pair: id: tal-2377c11, start20898,
end:28230
first_mate: start:27473,
end:28230
no second
mate
pair: id: tal-2426e12, start21975,
end:27562
first_mate: start:21975,
end:23008
second_mate: start:26396,
end:27562
pair: id: tal-2365h10, start22843,
end:31944
first_mate: start:22843,
end:23184
no second
mate
pair: id: tal-2388h09, start19016,
end:28238
first_mate: start:27475,
end:28238
no second mate
So it finds a lot of pairs that span the region and the start/end from
the pair is also correct but it only gives me both individual mates in
one case:
pair: id: tal-2426e12, start21975,
end:27562
first_mate: start:21975,
end:23008
second_mate: start:26396, end:27562
In this case, both pairs are actually inside the query region (at least
partially) whereas in the other cases, one of the mates is not inside,
e.g. this one:
pair: id: tal-2388h09, start19016,
end:28238
first_mate: start:27475,
end:28238
no second mate
> get this read pair from the BAM file:
$ samtools view clones.bam | grep tal-2388h09
tal-2388h09 99 tal12 19016 205
36H9M1D14M1D664M1D16M1D21M1D28M1D15M1D10M1D12M1D7M1D8M1D5M =
27475 9223
CTTTGGATGAAATAGTTTTTAAATAATACTTATTAAATATTAAATATATAACACATAAATAAGTATTGATGCAAATTTTAAAGTATTATAGAAAACTAGGTTTGATTATATTGTTATACTGTACTTTAAGAGGAGAGAGATAAGATATCTTTGCTCTTTTAATATATAAATTTAGATAAATATTCGTTAAATTTTCTACATAGTTATTTTTTATCTTATATATTATACTGCTATAGTTATCAATGTATATACATTCAAATAATTTATTAAAAATTCTATATTATATTAATTCTATGATAAAATAATCCTGTTTGTGATTTAAAAAATGATGATTCAATAAAAACTAATAATATAATACGAGTTAATATGGAATAATAAAATGGCATTTAACATGAATTTAGTCTTTAACCTTTTCTTTGTTTGTCAAGTTTTTTAAAACATAAAACCACACATTTCAAAATGGATTTTTAGCAAATATATAAAAATTATACATTTATAATGTATTGTTATGCGTCTTTTCGATAGAATCAATATTTAATTATATGAAGTTTCCACAATAAAATAATATTTAATATTATTTATTAGTAGAGTATTTGATTATATATATAGGCATATAATAATAACTCTAGTTCTATCTACCATATTATTTATAATTATTATAACAAAATGTGATATGAAATTTTATTATATACTTATATTATTTTTTTAACTATTTTAAAATATATTTATTTATACCTCAAAACTATAAAATTGAAATTATTAATAATAATCTAATATATACCTTTATAAAAATAAACGTATAAACTAAT
><:4/+1+*)+4>BEH=9-,,66IIIIIIIIEDA>>>>A at DDFFIHHHHHITIIIIIHIIHHHHHHIYYYYYTTTYDDDHDDDDDDDIIINNTNHHHHHIYYYYIIIIIINNNNTTIIIIIIIIIIITTNTTTTTYYYTTTTTTYYNNNNNNLLLLLLLLLLLNNNNNNTTTTTTTTTTTTTNNTNNNTTTYYTLLLLLLTTTTTTTTYTTTNNNNNNTTTTTTTNNTTNNTTTTTTTTTTYYTTTTTTTNNNNNNTTTTTTTYYTTTTTTYYYTTTTTTTTTTTYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYTTTTTTTTTYYYYYYYYYYYYYYYYYYYYYYYTTTTTYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYTTTTYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYTTNNINTTYYYYYYYYYYYYYYYTTTTTTYYYYYYYYYYYYTTTTTTYYYYYYYYYYYYYTTTTTOOOIFFFIFIIOICC>>II@>>>>>>C>>>>>>CIBECCCHIIOOOOOOOOTTTIIFDDEIQQA:55839AA>99>@IIIIII>>::;;I;>>CC>>>>>@III<::=>AAA<>>>>I>:>>99:>842225006824855;5>68844//.//00:>::338:99<:/-+*-./0)((((+00+..,++(((+-()(*((((()*)***))3)''')*..+*++((*1++--+*''''((+/)*42.((***)+,+('*'''*((''''((,'%%''''''''(
AS:i:614 MS:i:50
tal-2388h09 147 tal12 27475 205 1H764M40H =
19016 -9223
ATTAAATCGGTATCGCCAACACAATGAGTATAATCATTGTCAAATATGCGTTTGTAAGTATATTCATTGTCACATTCACGTTTGTAAGTATATTCATTGTCACATTCACGTTTGTAAGTATAATCATTGGAACGTTCATTTTTGTAAGTATAATCATTGGTACGTTCATTTTTGTGAGTATAATCATTGGTACGTTCATTTTTATGAGTATAATCATTGGTACGTTCATTTTTGTGAGTATAATCATTGGTACGTTCATTTTTATGAGTATAATCATTGGAACGTTCATTTTTGTAAGTATAATCATTGGTACGTTCATTTTTGTGAGTATAATCATTGGTACGTTCATTTTTATGAGTATAATCATTGGAACGTTCATTTTTATGAGTATAATCATTGGAACGTTCATTTTTGTAAGTATAATCATTGGTACGTTCATTTTTATGAGTATAATCATTGGTACGTTCATTTTTATGAGTATAATCATTGGAACGTTCATTTTTATGAGTATAATCATTGGAACGTTCATTTTTGTAAGTATAATCATTGGAACGTTCATTTTTATGAGTATAATCATTGGTACGTTCATTTTTGTGAGTATAATCATTGGTACGTTCATTTTTGTGAGTATAATCATTGGTACGTTCATTTTTATGAGTATAATCATTGGAACGTTCATTTTTATGAGTATAATCATTGGTACGTTCATTTTTGTGAGTATAATCATTGGTACGTTAATTTTTGTGAGTATAATCATTGGAACG
(((0))*,-1-../2((())03---03266300271+*.-0-*''''+*,+/+))*-05330+)..4>7=77273911**((+20+03688633:93036<8;::5:<99379>>::>>>:57:<:7--)))1435::333228>::>II>::>A>>3/.958677AA=AA:>:==IIII8338<>A>>>>IIIIIIIIYYYYYKKYYYMIFFFFEIIIMI::4..8AIIC>9>=EIQQQMCAAAAAACIIIIAICIIIOOYTIIIMOQQMIIIIC>>AAABCCCCCEAI>C>>IQQIIIIIIIIIIKKYYYYYYYYYYYYYYYYYYYYYTIIIIIIYYYYTNINNNTYYYYYYYYYYYYYTTTTTYYYYYYYYYYYYYYYYYYYYYYYYYTTTTTTYYYYYYYYYYYYYYYSSYYYYYYYTTTTTTYYYYYYYYYYYYYYYYYYYYYYYYTTTTTTTTTTTTTTTTTTYYYYYTTTTTTYYYTTNNNNTTYYYYYYYYYTTTTLLTTNNTTTTYTTTTTTYYYYYYYYYYYTTOOKKKLKOOTYYYYYYYYYYYYYYYYTNNNNNNNNNTTTNYNNNNTNNNNTTYYYYYYYYTTNNNNTTYNNNNNITTTTTYYYYYYYYYYTTNNIIIIIDIIIIHTNNNNTTYYYYTNNNIIIIIITTTINIIIINNNNTTTYYYYIHHHDDHHDDIHDDGDFFFTIIINTTYYYYTTTTHHHHCCIIIHIHHHHCAI9:++**1168>ACCIIDDDDDDI>>>>>?NNN
AS:i:688 MS:i:50
So the read in the first line starts before the start of the query
region and is not accessible via $pair->get_SeqFeatures although this is
a valid pair.
Am I doing something wrong, is this the desired behaviour or is it a bug?
Thanks for your help!
--
The Wellcome Trust Sanger Institute is operated by Genome Research
Limited, a charity registered in England with number 1021457 and a
company registered in England with number 2742969, whose registered
office is 215 Euston Road, London, NW1 2BE.
More information about the Bioperl-l
mailing list