Hello
Im trying to use Bio:: DB::Sam to get list of reads covering certain position .Sam file using pileup method using this code
So this code only get names of reads which cover position 1240 in reference from $ARGV[0].sam file. The problem is that the output list of read neames is definitely not complete. I.e. pileup ignores some reads.
Trying to find the root of problem i analysed reads which alignment starts from prosition 1209 (just for example) and found that only the first N of such reads assessed by puleup and the last M of reads are ignored. But after those not ignored N reads it doesnt ignore all reads starting from other positions (i.e. 1150,1048,1231 and any others). So as i understand the problem must be in some kind of sorting, but i havent found any information about sorting should be done in manual and i actually do sort .bam file.
Also i found that using $sam->pileup("GBA:1239-1241",$callback) instead of $sam->pileup("GBA",$callback) makes another list of reads, which differs from the previous, but i could found any logic in these differrences
So what's wrong?
Here are link to some test data: Sam file - https://cloud.mail.ru/public/d9e9de482926/In.sam
Reference - https://cloud.mail.ru/public/4394275..._gla_ids.fasta
$: perl test.pl In gaa_gba_gla_ids.fasta > out
Im trying to use Bio:: DB::Sam to get list of reads covering certain position .Sam file using pileup method using this code
Code:
use strict; use warnings; use Bio::DB::Sam; `samtools view -Sb $ARGV[0].sam > $ARGV[0].bam`; `samtools sort $ARGV[0].bam $ARGV[0].sorted`; `mv $ARGV[0].sorted.bam $ARGV[0].bam`; `samtools index $ARGV[0].bam`; my $sam = Bio::DB::Sam->new( -bam => "$ARGV[0].bam", -fasta => "$ARGV[1]", ); my $callback = sub { my ($seqid,$pos,$p) = @_; if ($pos eq "1240") { for my $pileup (@$p) { my $b = $pileup->alignment; print $b->qname,"\n"; } } }; $sam->pileup("GBA",$callback); exit;
Trying to find the root of problem i analysed reads which alignment starts from prosition 1209 (just for example) and found that only the first N of such reads assessed by puleup and the last M of reads are ignored. But after those not ignored N reads it doesnt ignore all reads starting from other positions (i.e. 1150,1048,1231 and any others). So as i understand the problem must be in some kind of sorting, but i havent found any information about sorting should be done in manual and i actually do sort .bam file.
Also i found that using $sam->pileup("GBA:1239-1241",$callback) instead of $sam->pileup("GBA",$callback) makes another list of reads, which differs from the previous, but i could found any logic in these differrences
So what's wrong?
Here are link to some test data: Sam file - https://cloud.mail.ru/public/d9e9de482926/In.sam
Reference - https://cloud.mail.ru/public/4394275..._gla_ids.fasta
$: perl test.pl In gaa_gba_gla_ids.fasta > out