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

[Gowthaman Ramasamy]

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);