[Bioperl-l] Getting pileup consensus from BAM files using Bio::DB::Sam

[Lincoln Stein]

Samtools is running MAQ on the pileup. You could either implement MAQ in
perl, or come up with your own consensus caller.

Lincoln

On Tue, Aug 3, 2010 at 1:29 AM, Gowthaman Ramasamy <
gowthaman.ramasamy at seattlebiomed.org> wrote:

> Hi List,
> I am trying to find out the consensus using pileup via Bio::DB::Sam. Using
> the following script I could parse out the ref_base and different bases from
> reads at that position. Though, I am not able to find a method to derive
> consensus. Similar to the values produced by "samtools pileup -c -f
> xxxxxx.fasta yyyyyyy.bam".
>
> The script I use now retrives ref base, query bases for each position. How
> do I improve it to get the consensus?
>
> Thanks very much in advance,
> Gowthaman
>
>
> use Bio::DB::Sam;
>
> my $bam = Bio::DB::Sam->new(-bam => 'something.bam',
>                            -fasta => 'something.fasta'
>                           );
>
> my $cb = sub {
>                        my ($seqid, $pos, $pileups) = @_;
>                        my $refBase = $bam->segment($seqid, $pos,
> $pos)->dna;
>                        print "\n$pos\t$refBase=>";
>                        for my $pileup (@$pileups){
>                                my $al = $pileup->alignment;
>                                my $qBase = substr($al->qseq, $pileup->qpos,
> 1);
>                                print "$qBase,";
>                                }
>                        };
>
> $bam->pileup('Lin.chr10i', $cb);
>
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at lists.open-bio.org
> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>



-- 
Lincoln D. Stein
Director, Informatics and Biocomputing Platform
Ontario Institute for Cancer Research
101 College St., Suite 800
Toronto, ON, Canada M5G0A3
416 673-8514
Assistant: Renata Musa <Renata.Musa at oicr.on.ca>