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

[Scott Cain]

Hi Gowthaman,

I don't see a method to extract the consensus.  You are welcome to
submit a patch :-)

Scott


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
>



-- 
------------------------------------------------------------------------
Scott Cain, Ph. D.                                   scott at scottcain dot net
GMOD Coordinator (http://gmod.org/)                     216-392-3087
Ontario Institute for Cancer Research