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