Getting pileup consensus from BAM files using Bio::DB::Sam

[ragowthaman]

Hi List,
I am trying to find out the consensus using pileup via Bio:B::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:B::Sam;

my $bam = Bio:B::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);

[lh3]

There is no perl binding for samtools-like consensus calling.

[ragowthaman]

Oh! Thanks. I'll be good reading the pileup outfile then.