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

ragowthaman · 08-02-2010, 11:13 PM

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

lh3 · 08-03-2010, 07:19 AM

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

ragowthaman · 08-03-2010, 08:21 AM

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

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Consensus FASTA from BAM files	mixter	Bioinformatics	8	11-14-2012 11:20 AM
Pileup from BAM or SAM	nilmot13	General	7	02-03-2012 02:54 AM
Consensus part from sequence read(fastq) and align(BAM) files	culmen	Bioinformatics	5	12-21-2010 03:57 AM
Max Mapping / SNP / Consensus Quality in pileup files	agc	Bioinformatics	1	11-22-2010 06:42 AM

08-02-2010, 11:13 PM	#1
ragowthaman Member Location: Seattle, USA Join Date: Nov 2009 Posts: 12	Getting pileup consensus from BAM files using Bio::DB::Sam 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);

08-03-2010, 07:19 AM	#2
lh3 Senior Member Location: Boston Join Date: Feb 2008 Posts: 611	There is no perl binding for samtools-like consensus calling.

08-03-2010, 08:21 AM	#3
ragowthaman Member Location: Seattle, USA Join Date: Nov 2009 Posts: 12	Oh! Thanks. I'll be good reading the pileup outfile then.