[Bioperl-l] bio::db::sam coverage problem

[Peter Huang]

Hi there,

I am trying to use Bio::DB::SAM to process my sam file generated from samtools. However, I encountered one problem.

Let's say I want to extract a base from a specific position that contains a SNP. I'd like to obtain count information regarding the individual bases at this position for later analysis.
If I use samtools view
/share/apps/software/samtools/samtools   view  LibG_sort.bam   ARP6:5134-5134   -o filteredresults.sam
less filteredresults.sam | wc
  32333  614327 17194532

However, when I tried to use the Bio::DB::SAM, the coverage is only 8088.

Here is what I have in my perl script:
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
my $sam = Bio::DB::Sam->new( -bam => $bamFile, -fasta => $reference );
my $total;

my $cb = sub {
    my ($seqid, $site, $pileups) = @_;
if( $site == $position ) {
    $total = scalar @$pileups;
        my $refBase = $sam->segment($seqid, $site, $site)->dna;

        for my $pileup (@$pileups){
            my $al = $pileup->alignment;
            my $qSeq = $al->qseq;

            my $qStr = substr($qSeq, $pileup->qpos, 1);
        }
    }
}
$sam->pileup("ARP6:5134-5134", $cb);
print "total is $total\n";
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
And the output for the total is 8088

I am curious if there is anything I did wrong or I missed anything. Your help is greatly appreciated! Thanks!

Best,

Peter