bioperl for samtools

[m_elena_bioinfo]

Hi users,
I'm using Bio:B::Sam to analyse my data from Illumina.
I have bam, sam and fasta files,but I don't know how use bio-samtools and its modules to interface to the pileup routines.
Anyone could help me?!
thanks a lot,
ME

[krobison]

Which of the examples from the online documentation have you tried? The SNP caller example is a good illustration of writing a "walker" which will visit regions of the alignment and compute something.

[m_elena_bioinfo]

Hi Krobison, thanx for oyur quick reply! I'm using the indications from this site: http://search.cpan.org/~lds/Bio-SamT.../Bio/DB/Sam.pm . But the examples are not so easy to understand...How can I extract pileup from bam? Or from my pileup file, how can I extract the information? I really can't understand the example or the module!

[bbl]

Hello there,
I recently have encoutered the same problem. Have you got it sorted? Perhaps, you can divulge a bit information... thanx

[krobison]

Okay, here is one possible exam answer :-) The $sam argument needs to be the Perl object created from the .bam file. The callback scheme is probably not how I would have designed things, but that's a matter of personal design taste & one must work with the way the Perl interface is.

This code could be a little rusty; I mostly use another tool now for SNP calling. Wouldn't shock me if there is a more compact & clever way to code this either.

package SnpCaller;
## please give Keith Robison, Infinity Pharmaceuticals & Omics! Omics! credit
## if you find this useful. 
## If you don't please post your improvement or criticism!

sub new
{
    my ($class)=@_;
    my $this={'minQ'=>25, 'minMinor'=>2, 'minDepth'=>4};
    bless $this,$class;
    return $this;
}

### variation on "extremely stupid" SNP caller from Sam.pm docs

sub callSnps
{
    my ($this,$sam,$region)=@_;
    # region in form chr2:1-10
    my @SNPs=();
    my $snpCaller=sub 
    {
        my ($seqid,$pos,$p)=@_;
        # any speed gain from fetching long segment at start?
        my $refbase=$sam->segment($seqid,$pos,$pos)->dna;

        $refbase=~tr/a-z/A-Z/;
        my %bases=();
        my %maxQual=();
        return if (scalar(@{$p})<$this->{'minDepth'});
        for my $pileup(@{$p})
        {
            my $b=$pileup->b;
            next if ($pileup->indel); # doesn't handle indels; need to in future
            my $qbase=substr($b->qseq,$pileup->qpos,1);
            next unless ($qbase=~/[ATCG]/i);
            my $qscore=ord(substr($b->qscore,$pileup->qpos,1))-33;
            next unless ($qscore>$this->{'minQ'});
            $bases{$qbase}++;
            $maxQual{$qbase}=$qscore if ($qscore>$maxQual{$qbase});
        }
        my $store=1;
        if (scalar(keys %bases)==1)
        {
            $store=0 if (defined $bases{$refbase});
        }
        else
        {
            my @counts=sort {$b<=>$a} values %bases;
            $store=0 if ($counts[1]<$this->{'minMinor'});
        }
        push(@SNPs,{'pos'=>$pos,'ref'=>$refbase, 'bases'=>\%bases,'maxQual'=>\%maxQual}) if ($store); 
    };
    $sam->pileup($region,$snpCaller);
    return @SNPs;
}