[Bioperl-l] Bio::Tools::SeqStats count_codons

Govind Chandra govind.chandra at bbsrc.ac.uk
Tue Sep 27 11:15:51 EDT 2005


The count_codons method in Bio::Tools::SeqStats seems to be case
sensitive. See script below. Is this intended? Or am I doing something
silly.

### script begin ###
use Bio::SeqIO;
use Bio::Seq;
use Bio::Tools::SeqStats;
$seqio=Bio::SeqIO->new('-file' => 'Ssc.fas',
                       '-format' => 'fasta');

$seqobj=$seqio->next_seq();

open(ORFS, "<ordered_all_orfs.out");

while(<ORFS>) {
chomp();
($strand, $start, $end, $nt_len, $aa_len)=split(/\s+/, $_);
#print("$strand   $start   $end   $nt_len   $aa_len\n");

if($strand==1) {
$cdsobj=$seqobj->trunc($start, $end);
}
elsif($strand==-1) {
$cdsobj=$seqobj->trunc($start, $end)->revcom();
}

$tseq=uc($cdsobj->seq());  # works
#$tseq=$cdsobj->seq(); # does not work. given and ambiguous count only

$codobj=Bio::Seq->new('-seq' => $tseq,
		      '-alphabet' => 'dna'
		      );

print(">$strand $aa_len\n",$cdsobj->seq(),"\n");

$ss=Bio::Tools::SeqStats->new(-seq => $codobj);

$refcount = $ss->count_codons();

foreach $codon (keys(%$refcount)) {
print("$codon   $refcount->{$codon}\n");
}
}

### script end ###

Cheers

Govind
John Innes Centre
Norwich UK.




More information about the Bioperl-l mailing list