[Bioperl-l] Jitterbug Issue 943
Roger Hall
Roger Hall" <roger@iosea.com
Thu, 7 Jun 2001 11:41:04 -0700
This is a multi-part message in MIME format.
------=_NextPart_000_006C_01C0EF46.BBBE6460
Content-Type: multipart/alternative;
boundary="----=_NextPart_001_006D_01C0EF46.BBBE6460"
------=_NextPart_001_006D_01C0EF46.BBBE6460
Content-Type: text/plain;
charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable
I enhanced SeqStats to use either a custom alphabet or a custom =
aplhabet/weight hash
The incorrect testing script in SeqStats has been replaced with the code =
from seqstats.pl, which includes the corrected original tests and =
enhancement tests.
Please find attached SeqStats.pm and seqstats.pl.
Did we want this enhancement and is this effective? :}
Thanks!
Roger Hall
G1440
------=_NextPart_001_006D_01C0EF46.BBBE6460
Content-Type: text/html;
charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=3DContent-Type content=3D"text/html; =
charset=3Diso-8859-1">
<META content=3D"MSHTML 5.50.4134.600" name=3DGENERATOR>
<STYLE></STYLE>
</HEAD>
<BODY bgColor=3D#ffffff>
<DIV><FONT face=3DArial size=3D2>
<DIV><FONT face=3DArial size=3D2><FONT face=3DArial size=3D2>I enhanced =
SeqStats to use=20
either a custom alphabet or a custom aplhabet/weight hash</FONT>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>The incorrect =
testing script in=20
SeqStats has been replaced with the code from seqstats.pl, which =
includes=20
the corrected original tests and enhancement tests.</FONT></DIV>
<DIV> </DIV>
<DIV><FONT face=3DArial size=3D2>Please find attached SeqStats.pm and=20
seqstats.pl.</FONT></DIV>
<DIV> </DIV>
<DIV>Did we want this enhancement and is this effective? :}</DIV>
<DIV> </DIV>
<DIV>Thanks!</DIV>
<DIV> </DIV>
<DIV>Roger Hall</DIV>
<DIV>G1440</DIV>
<DIV> </DIV></FONT></DIV></BODY></HTML>
------=_NextPart_001_006D_01C0EF46.BBBE6460--
------=_NextPart_000_006C_01C0EF46.BBBE6460
Content-Type: application/octet-stream;
name="SeqStats.pm"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="SeqStats.pm"
# $Id: SeqStats.pm,v 1.10.2.1 2001/03/03 08:28:58 heikki Exp $
#
# BioPerl module for Bio::Tools::SeqStats
#
# Cared for by
#
# Copyright Peter Schattner
#
# You may distribute this module under the same terms as perl itself
# POD documentation - main docs before the code
=3Dhead1 NAME
Bio::Tools::SeqStats - Object holding statistics for one particular =
sequence
=3Dhead1 SYNOPSIS
# build a primary nucleic acid or protein sequence object somehow
# then build a statistics object from the sequence object=20
$seqobj =3D Bio::PrimarySeq->new(-seq=3D>'ACTGTGGCGTCAACTG',=20
-moltype=3D>'dna',=20
-id=3D>'test');
$seq_stats =3D Bio::Tools::SeqStats->new(-seq=3D>$seqobj);
# obtain a hash of counts of each type of monomer=20
# (ie amino or nucleic acid)
print "\nMonomer Counts using statistics object\n";
$seq_stats =3D Bio::Tools::SeqStats->new(-seq=3D>$seqobj);
$hash_ref =3D $seq_stats->count_monomers(); # eg for DNA sequence
foreach $base (sort keys %$hash_ref) {
print "Number of bases of type ", $base, "=3D ", =
%$hash_ref->{$base},"\n";
}
# or obtain the count directly without creating a new statistics object
print "\nMonomer Counts without statistics object\n";
$hash_ref =3D Bio::Tools::SeqStats->count_monomers($seqobj);
foreach $base (sort keys %$hash_ref) {
print "Number of bases of type ", $base, "=3D ", =
%$hash_ref->{$base},"\n";
}
# obtain hash of counts of each type of codon in a nucleic acid =
sequence
print "\nCodon Counts usin statistics object\n";
$hash_ref =3D $seq_stats-> count_codons(); # for nucleic acid sequence
foreach $base (sort keys %$hash_ref) {
print "Number of codons of type ", $base, "=3D ", =
%$hash_ref->{$base},"\n";
}
# or
print "\nCodon Counts without statistics object\n";
$hash_ref =3D Bio::Tools::SeqStats->count_codons($seqobj);
foreach $base (sort keys %$hash_ref) {
print "Number of codons of type ", $base, "=3D ", =
%$hash_ref->{$base},"\n";
}
# Obtain the molecular weight of a sequence. Since the sequence may =
contain # ambiguous monomers, the molecular weight is returned as a =
(reference to) a # two element array containing greatest lower bound =
(GLB) and least upper bound # (LUB) of the molecular weight=20
$weight =3D $seq_stats->get_mol_wt();
print "\nMolecular weight (using statistics object) of sequence ", =
$seqobj->id(),=20
" is between ", $$weight[0], " and " ,=20
$$weight[1], "\n";
# or
$weight =3D Bio::Tools::SeqStats->get_mol_wt($seqobj);
print "\nMolecular weight (without statistics object) of sequence ", =
$seqobj->id(),=20
" is between ", $$weight[0], " and " ,=20
$$weight[1], "\n";
# use a custom alphabet/weight hash
$seqobj =3D Bio::PrimarySeq->new ( -seq =3D> =
'BACDDCDDCACADDAABBDADABBD',
-id =3D> 'TestFragment-12',
-accession_number =3D> 'Test',
-moltype =3D> 'rna'
);
print "\nMonomer Counts using custom alphabet\n";
my @alphabet =3D qw(A B C D);
$seq_stats =3D =
Bio::Tools::SeqStats->new(-seq=3D>$seqobj,-alphabet=3D>\@alphabet);
$hash_ref =3D $seq_stats->count_monomers(); # eg for DNA sequence
foreach $base (sort keys %$hash_ref) {
print "Number of bases of type ", $base, "=3D ", =
%$hash_ref->{$base},"\n";
}
# Obtain the molecular weight of a sequence with a custom alphabet
my %weights =3D =
('A'=3D>[9,10],'B'=3D>[8,10],'C'=3D>[11,11],'D'=3D>[12,13]);
$seq_stats =3D =
Bio::Tools::SeqStats->new(-seq=3D>$seqobj,-weights=3D>\%weights);
$weight =3D $seq_stats->get_mol_wt();
print "\nMolecular weight (using custom weights) of sequence ", =
$seqobj->id(),=20
" is between ", $$weight[0], " and " ,=20
$$weight[1], "\n";
=3Dhead1 DESCRIPTION
Bio::Tools::SeqStats is a lightweight object for the calculation of
simple statistical and numerical properties of a sequence. By
"lightweight" I mean that only "primary" sequences are handled by the
object. The calling script needs to create the appropriate primary
sequence to be passed to SeqStats if statistics on a sequence feature
are required. Similarly if a codon count is desired for a
frame-shifted sequence and/or a negative strand sequence, the calling
script needs to create that sequence and pass it to the SeqStats
object.
SeqStats can be called in two distinct manners. If only a single
computation is required on a given sequence object, the method can be
called easily using the SeqStats object directly:
$weight =3D Bio::Tools::SeqStats->get_mol_wt($seqobj);
Alternately, if several computations will be required on a given
sequence object, an "instance" statistics object can be constructed
and used for the method calls:
$seq_stats =3D Bio::Tools::SeqStats->new($seqobj);
$monomers =3D $seq_stats->count_monomers();
$codons =3D $seq_stats->count_codons();
$weight =3D $seq_stats->get_mol_wt();
As currently implemented the object can return the following values from =
a sequence:
* The molecular weight of the sequence: get_mol_wt()
* The number of each type of monomer present: count_monomers()
* The number of each codon present in a nucleic acid sequence: =
count_codons()
For dna (and rna) sequences, single-stranded weights are returned. The
molecular weights are calculated for neutral - ie not ionized -
nucleic acids. The returned weight is the sum of the
base-sugar-phosphate residues of the chain plus one weight of water to
to account for the additional OH on the phosphate of the 5' residue
and the additional H on the sugar ring of the 3' residue. Note that
this leads to a difference of 18 in calculated molecular weights
compared to some other available programs (eg Informax VectorNTI).
Note that since sequences may contain ambiguous monomers (eg "M"
meaning "A" or "C" in a nucleic acid sequence), the method get_mol_wt
returns a two-element array containing the greatest lower bound and
least upper bound of the molecule. (For a sequence with no ambiguous
monomers, the two elements of the returned array will be equal.) The
method count_codons() handles ambiguous bases by simply counting all
ambiguous codons together and issuing a warning to that effect.
=3Dhead1 DEVELOPERS NOTES
Ewan moved it from Bio::SeqStats to Bio::Tools::SeqStats
=3Dhead1 FEEDBACK
=3Dhead2 Mailing Lists
User feedback is an integral part of the evolution of this and other
Bioperl modules. Send your comments and suggestions preferably to one
of the Bioperl mailing lists. Your participation is much appreciated.
bioperl-l@bioperl.org - General discussion
http://bio.perl.org/MailList.html - About the mailing lists
=3Dhead2 Reporting Bugs
Report bugs to the Bioperl bug tracking system to help us keep track
the bugs and their resolution.
Bug reports can be submitted via email or the web:
bioperl-bugs@bio.perl.org
http://bio.perl.org/bioperl-bugs/
=3Dhead1 AUTHOR - Peter Schattner
Email schattner@alum.mit.edu
=3Dhead1 APPENDIX
The rest of the documentation details each of the object methods. =
Internal methods are
usually preceded with a _
=3Dcut
package Bio::Tools::SeqStats;
use strict;
use vars qw(@ISA %Alphabets %Alphabets_strict $amino_weights=20
$rna_weights $dna_weights %Weights );
use Bio::Seq;
use Bio::Root::RootI;
@ISA =3D qw(Bio::Root::RootI);
BEGIN {=20
%Alphabets =3D (
'dna' =3D> [ qw(A C G T R Y M K S W H B V D X N) ],
'rna' =3D> [ qw(A C G U R Y M K S W H B V D X N) ],
'protein' =3D> [ qw(A R N D C Q E G H I L K M F
P S T W X Y V B Z *) ], # sac: added B, Z
);
# SAC: new strict alphabet: doesn't allow any ambiguity characters.
%Alphabets_strict =3D (
'dna' =3D> [ qw( A C G T ) ],
'rna' =3D> [ qw( A C G U ) ],
'protein' =3D> [ qw(A R N D C Q E G H I L K M F
P S T W Y V) ],
);
# IUPAC-IUB SYMBOLS FOR NUCLEOTIDE NOMENCLATURE:
# Cornish-Bowden (1985) Nucl. Acids Res. 13: 3021-3030.
# Amino Acid alphabet
# ------------------------------------------
# Symbol Meaning
# ------------------------------------------
my $amino_A_wt =3D 89.09;
my $amino_C_wt =3D 121.15;
my $amino_D_wt =3D 133.1;
my $amino_E_wt =3D 147.13;
my $amino_F_wt =3D 165.19;
my $amino_G_wt =3D 75.07;
my $amino_H_wt =3D 155.16;
my $amino_I_wt =3D 131.18;
my $amino_K_wt =3D 146.19;
my $amino_L_wt =3D 131.18;
my $amino_M_wt =3D 149.22;
my $amino_N_wt =3D 132.12;
my $amino_P_wt =3D 115.13;
my $amino_Q_wt =3D 146.15;
my $amino_R_wt =3D 174.21;
my $amino_S_wt =3D 105.09;
my $amino_T_wt =3D 119.12;
my $amino_V_wt =3D 117.15;
my $amino_W_wt =3D 204.22;
my $amino_Y_wt =3D 181.19;
$amino_weights =3D {
'A' =3D> [$amino_A_wt, $amino_A_wt], # Alanine
'B' =3D> [$amino_N_wt, $amino_D_wt], # Aspartic Acid, Asparagine
'C' =3D> [$amino_C_wt, $amino_C_wt], # Cystine
'D' =3D> [$amino_D_wt, $amino_D_wt], # Aspartic Acid
'E' =3D> [$amino_E_wt, $amino_E_wt], # Glutamic Acid
'F' =3D> [$amino_F_wt, $amino_F_wt], # Phenylalanine
'G' =3D> [$amino_G_wt, $amino_G_wt], # Glycine
'H' =3D> [$amino_H_wt, $amino_H_wt], # Histidine
'I' =3D> [$amino_I_wt, $amino_I_wt], # Isoleucine
'K' =3D> [$amino_K_wt, $amino_K_wt], # Lysine
'L' =3D> [$amino_L_wt, $amino_L_wt], # Leucine
'M' =3D> [$amino_M_wt, $amino_M_wt], # Methionine
'N' =3D> [$amino_N_wt, $amino_N_wt], # Asparagine
'P' =3D> [$amino_P_wt, $amino_P_wt], # Proline
'Q' =3D> [$amino_Q_wt, $amino_Q_wt], # Glutamine
'R' =3D> [$amino_R_wt, $amino_R_wt], # Arginine
'S' =3D> [$amino_S_wt, $amino_S_wt], # Serine
'T' =3D> [$amino_T_wt, $amino_T_wt], # Threonine
'V' =3D> [$amino_V_wt, $amino_V_wt], # Valine
'W' =3D> [$amino_W_wt, $amino_W_wt], # Tryptophan
'X' =3D> [$amino_G_wt, $amino_W_wt], # Unknown
'Y' =3D> [$amino_Y_wt, $amino_Y_wt], # Tyrosine
'Z' =3D> [$amino_Q_wt, $amino_E_wt], # Glutamic Acid, Glutamine
};
# Extended Dna / Rna alphabet
use vars ( qw($C $O $N $H $P $water) );
use vars ( qw($adenine $guanine $cytosine $thymine =
$uracil));
use vars ( qw($ribose_phosphate $deoxyribose_phosphate $ppi));
use vars ( qw($dna_A_wt $dna_C_wt $dna_G_wt $dna_T_wt =20
$rna_A_wt $rna_C_wt $rna_G_wt $rna_U_wt));
use vars ( qw($dna_weights $rna_weights %Weights));
$C =3D 12.01;
$O =3D 16.00;
$N =3D 14.01;
$H =3D 1.01;
$P =3D 30.97;
$water =3D 18.015;
=20
$adenine =3D 5 * $C + 5 * $N + 5 * $H;
$guanine =3D 5 * $C + 5 * $N + 1 * $O + 5 * $H;
$cytosine =3D 4 * $C + 3 * $N + 1 * $O + 5 * $H;
$thymine =3D 5 * $C + 2 * $N + 2 * $O + 6 * $H;
$uracil =3D 4 * $C + 2 * $N + 2 * $O + 4 * $H;
=20
$ribose_phosphate =3D 5 * $C + 7 * $O + 9 * $H + 1 * $P; =
#neutral (unionized) form
$deoxyribose_phosphate =3D 5 * $C + 6 * $O + 9 * $H + 1 * $P;
=20
# the following are single strand molecular weights / base
$dna_A_wt =3D $adenine + $deoxyribose_phosphate - $water;
$dna_C_wt =3D $cytosine + $deoxyribose_phosphate - $water;
$dna_G_wt =3D $guanine + $deoxyribose_phosphate - $water;
$dna_T_wt =3D $thymine + $deoxyribose_phosphate - $water;
=20
$rna_A_wt =3D $adenine + $ribose_phosphate - $water;
$rna_C_wt =3D $cytosine + $ribose_phosphate - $water;
$rna_G_wt =3D $guanine + $ribose_phosphate - $water;
$rna_U_wt =3D $uracil + $ribose_phosphate - $water;
=20
$dna_weights =3D {
'A' =3D> [$dna_A_wt,$dna_A_wt], # Adenine
'C' =3D> [$dna_C_wt,$dna_C_wt], # Cytosine
'G' =3D> [$dna_G_wt,$dna_G_wt], # Guanine
'T' =3D> [$dna_T_wt,$dna_T_wt], # Thymine
'M' =3D> [$dna_C_wt,$dna_A_wt], # A or C
'R' =3D> [$dna_A_wt,$dna_G_wt], # A or G
'W' =3D> [$dna_T_wt,$dna_A_wt], # A or T
'S' =3D> [$dna_C_wt,$dna_G_wt], # C or G
'Y' =3D> [$dna_C_wt,$dna_T_wt], # C or T
'K' =3D> [$dna_T_wt,$dna_G_wt], # G or T
'V' =3D> [$dna_C_wt,$dna_G_wt], # A or C or G
'H' =3D> [$dna_C_wt,$dna_A_wt], # A or C or T
'D' =3D> [$dna_T_wt,$dna_G_wt], # A or G or T
'B' =3D> [$dna_C_wt,$dna_G_wt], # C or G or T
'X' =3D> [$dna_C_wt,$dna_G_wt], # G or A or T or =
C
'N' =3D> [$dna_C_wt,$dna_G_wt], # G or A or T or =
C
};
=20
$rna_weights =3D {
'A' =3D> [$rna_A_wt,$rna_A_wt], # Adenine
'C' =3D> [$rna_C_wt,$rna_C_wt], # Cytosine
'G' =3D> [$rna_G_wt,$rna_G_wt], # Guanine
'U' =3D> [$rna_U_wt,$rna_U_wt], # Uracil
'M' =3D> [$rna_C_wt,$rna_A_wt], # A or C
'R' =3D> [$rna_A_wt,$rna_G_wt], # A or G
'W' =3D> [$rna_U_wt,$rna_A_wt], # A or U
'S' =3D> [$rna_C_wt,$rna_G_wt], # C or G
'Y' =3D> [$rna_C_wt,$rna_U_wt], # C or U
'K' =3D> [$rna_U_wt,$rna_G_wt], # G or U
'V' =3D> [$rna_C_wt,$rna_G_wt], # A or C or G
'H' =3D> [$rna_C_wt,$rna_A_wt], # A or C or U
'D' =3D> [$rna_U_wt,$rna_G_wt], # A or G or U
'B' =3D> [$rna_C_wt,$rna_G_wt], # C or G or U
'X' =3D> [$rna_C_wt,$rna_G_wt], # G or A or U or =
C
'N' =3D> [$rna_C_wt,$rna_G_wt], # G or A or U or =
C
};
=20
%Weights =3D (
'dna' =3D> $dna_weights,
'rna' =3D> $rna_weights,
'protein' =3D> $amino_weights,
);
}
sub new {
my($class,@args) =3D @_;
my $self =3D $class->SUPER::new(@args);
# Add Weights param (J943)
# my ($seqobj) =3D $self->_rearrange([qw(SEQ)],@args); #(J943)
my ($seqobj, $alphabet, $alphaweights) =3D $self->_rearrange([qw(SEQ =
ALPHABET WEIGHTS)],@args); #(J943)
unless ($seqobj->isa("Bio::PrimarySeqI")) {
$self->throw(" SeqStats works only on PrimarySeqI objects \n");
}
# If Weights param, set %Alphabet and %Weights (J943)
if ($alphaweights) { #(J943)
$Weights{$seqobj->moltype} =3D $alphaweights; #(J943)
my @alphabet =3D sort keys %{ $alphaweights }; #(J943)
$Alphabets{$seqobj->moltype} =3D [ @alphabet ]; #(J943)
$Alphabets_strict{$seqobj->moltype} =3D [ @alphabet ]; #(J943)
} #(J943)
# If just Alphabet param, set %Alphabet (J943)
if ($alphabet && !$alphaweights) { #(J943)
$Alphabets{$seqobj->moltype} =3D [ @$alphabet ]; #(J943)
$Alphabets_strict{$seqobj->moltype} =3D [ @$alphabet ]; #(J943)
} #(J943)
if ( !defined $seqobj->moltype || ! defined =
$Alphabets{$seqobj->moltype}) {
$self->throw("Must have a valid moltype defined for seq (".
join(",",keys %Alphabets));
}
$self->{'_seqref'} =3D $seqobj;
$self->{'_is_strict'} =3D _is_alphabet_strict($seqobj); # check the =
letters in the sequence
return $self;=20
}
=3Dhead2 count_monomers
Title : count_monomers
Usage : $rcount =3D $seq_stats->count_monomers();=20
or $rcount =3D $seq_stats->Bio::Tools::SeqStats->($seqobj);
Function: Counts the number of each type of monomer (amino acid or
base) in the sequence.
Example :
Returns : Reference to a hash in which keys are letters of the
genetic alphabet used and values are number of occurrences
of the letter in the sequence.
Args : None or reference to sequence object
Throws : Throws an exception if type of sequence is unknown (ie amino
or nucleic)or if unknown letter in alphabet. Ambiguous
elements are allowed.
=3Dcut
sub count_monomers{
my $rcount;
my %count =3D ();
my $seqobj;
my $_is_strict;
my $element =3D '';
my $_is_instance =3D 1 ;
my $self =3D shift @_;
my $object_argument =3D shift @_;
# First we need to determine if the present object is an instance
# object or if the sequence object has been passed as an argument
if (defined $object_argument) {
$_is_instance =3D 0;
}
# If we are using an instance object...
if ($_is_instance) {
if ($rcount =3D $self->{'_monomer_count'}) {
return $rcount; # return count if previously calculated
}
$_is_strict =3D $self->{'_is_strict'}; # retrieve "strictness"
$seqobj =3D $self->{'_seqref'};
} else {
# otherwise...
$seqobj =3D $object_argument;
=09
# Following two lines lead to error in "throw" routine
$seqobj->isa("Bio::PrimarySeqI") ||
$self->throw(" SeqStats works only on PrimarySeqI objects \n");
# is alphabet OK? Is it strict?
$_is_strict =3D _is_alphabet_strict($seqobj);=20
}
=09
my $alphabet =3D $_is_strict ? $Alphabets_strict{$seqobj->moltype} =
:
$Alphabets{$seqobj->moltype} ; # get array of allowed letters
=09
# convert everything to upper case to be safe
my $seqstring =3D uc $seqobj->seq(); =20
# For each letter, count the number of times it appears in
# the sequence
LETTER:
foreach $element (@$alphabet) {
# skip terminator symbol which may confuse regex
next LETTER if ($element eq '*');=20
$count{$element} =3D ( $seqstring =3D~ s/$element/$element/g);
}
=20
$rcount =3D \%count;
=20
if ($_is_instance) {
$self->{'_monomer_count'} =3D $rcount; # Save in case called again =
later
}
=20
return $rcount;
}
=3Dhead2 get_mol_wt
Title : get_mol_wt
Usage : $wt =3D $seqobj->get_mol_wt() or=20
$wt =3D Bio::Tools::SeqStats ->get_mol_wt($seqobj);
Function: Calculate molecular weight of sequence
Example :
Returns : Reference to two element array containing lower and upper
bounds of molecule molecular weight. (For dna (and rna)
sequences, single-stranded weights are returned.) If
sequence contains no ambiguous elements, both entries in
array are equal to molecular weight of molecule. =20
Args : None or reference to sequence object
Throws : Exception if type of sequence is unknown (ie not amino or
nucleic) or if unknown letter in alphabet. Ambiguous
elements are allowed.
=3Dcut
sub get_mol_wt {
my $seqobj;
my $_is_strict;
my $element =3D '';
my $_is_instance =3D 1 ;
my $self =3D shift @_;
my $object_argument =3D shift @_;
my ($weight_array, $rcount);
=20
if (defined $object_argument) {
$_is_instance =3D 0;
}
=20
if ($_is_instance) {=09
if ($weight_array =3D $self->{'_mol_wt'}) {
# return mol. weight if previously calculated
return $weight_array; =20
}
$seqobj =3D $self->{'_seqref'};
$rcount =3D $self->count_monomers();
} else {
$seqobj =3D $object_argument;
$seqobj->isa("Bio::PrimarySeqI") ||
die("Error: SeqStats works only on PrimarySeqI objects \n");
$_is_strict =3D _is_alphabet_strict($seqobj); # is alphabet OK?
$rcount =3D $self->count_monomers($seqobj);
}
# We will also need to know what type of monomer we are dealing with
=20
my $moltype =3D $seqobj->moltype();
# In general,the molecular weight is bounded below by the sum of the
# weights of lower bounds of each alphabet symbol times the number of
# occurrences of the symbol in the sequence. A similar upper bound on
# the weight is also calculated.
# Note that for "strict" (ie unambiguous) sequences there is an
# inefficiency since the upper bound =3D the lower bound (and is
# calculated twice). However, this decrease in performance will be
# minor and leads to (IMO) significantly more readable code.
my $weight_lower_bound =3D 0;
my $weight_upper_bound =3D 0;
my $weight_table =3D $Weights{$moltype};
# my $water =3D 18.015;
# compute weight of all the residues
foreach $element (keys %$rcount) {
$weight_lower_bound +=3D $$rcount{$element} * =
$$weight_table{$element}->[0];
$weight_upper_bound +=3D $$rcount{$element} * =
$$weight_table{$element}->[1];
}
if ($moltype =3D~ /protein/) {
# remove of H2O during peptide bond formation.
$weight_lower_bound -=3D $water * ($seqobj->length - 1);
$weight_upper_bound -=3D $water * ($seqobj->length - 1);
} else {
# Correction because phosphate of 5' residue has additional OH and
# sugar ring of 3' residue has additional H
$weight_lower_bound +=3D $water;
$weight_upper_bound +=3D $water;
}
$weight_lower_bound =3D sprintf("%.0f", $weight_lower_bound);
$weight_upper_bound =3D sprintf("%.0f", $weight_upper_bound);
$weight_array =3D [$weight_lower_bound, $weight_upper_bound];
if ($_is_instance) {
$self->{'_mol_wt'} =3D $weight_array; # Save in case called again =
later
}
return $weight_array;
}
=3Dhead2 count_codons
Title : count_codons
Usage : $rcount =3D $seqstats->count_codons (); or=20
$rcount =3D Bio::Tools::SeqStats->count_codons($seqobj);
Function: Counts the number of each type of codons in a given frame=20
for a dna or rna sequence.
Example :
Returns : Reference to a hash in which keys are codons of the genetic
alphabet used and values are number of occurrences of the
codons in the sequence. All codons with "ambiguous" bases
are counted together.
Args : None or reference to sequence object
Throws : an exception if type of sequence is unknown or protein.
=3Dcut
sub count_codons {
my $rcount =3D {};
my $codon ;
my $seqobj;
my $_is_strict;
my $element =3D '';
my $_is_instance =3D 1 ;
my $self =3D shift @_;
my $object_argument =3D shift @_;
=20
if (defined $object_argument) {
$_is_instance =3D 0;
}
=20
if ($_is_instance) {
if ($rcount =3D $self->{'_codon_count'}) {
return $rcount; # return count if previously calculated
}
$_is_strict =3D $self->{'_is_strict'}; # retrieve "strictness"
$seqobj =3D $self->{'_seqref'};
} else {
$seqobj =3D $object_argument;
$seqobj->isa("Bio::PrimarySeqI") ||
die(" Error: SeqStats works only on PrimarySeqI objects \n");
$_is_strict =3D _is_alphabet_strict($seqobj);
}
=20
# Codon counts only make sense for nucleic acid sequences
my $moltype =3D $seqobj->moltype();
=20
unless ($moltype =3D~ /[dr]na/) {
$seqobj->throw(" Codon counts only meaningful for dna or rna, not for =
$moltype sequences. \n");
}
=20
# If sequence contains ambiguous bases, warn that codons containing them =
will all be
# lumped together in the count.
=20
if (!$_is_strict ) {
$seqobj->warn(" Sequence $seqobj contains ambiguous bases. \n All =
codons with ambiguous bases will be added together in count. \n");
}
=20
my $seq =3D $seqobj->seq();
=20
# Now step through the string by threes and count the codons
=20
CODON:
while (length($seq) > 2) {
$codon =3D substr($seq,0,3);
$seq =3D substr($seq,3);
if ($codon =3D~ /[^ACTGU]/) {
$$rcount{'ambiguous'}++; #lump together ambiguous codons
next CODON;
}
if (!defined $$rcount{$codon}) {
$$rcount{$codon}=3D 1 ;
next CODON;
}
$$rcount{$codon}++; # default
}
=20
=20
if ($_is_instance) {
$self->{'_codon_count'} =3D $rcount; # Save in case called again later
}
=20
=20
=20
=20
return $rcount;
}
=3Dhead2 _is_alphabet_strict
Title : _is_alphabet_strict
Usage : =20
Function: internal function to determine whether there are=20
any ambiguous elements in the current sequence
Example :
Returns : 1 if strict alphabet is being used,=20
0 if ambiguous elements are present
Args :
Throws : an exception if type of sequence is unknown (ie amino or
nucleic) or if unknown letter in alphabet. Ambiguous
monomers are allowed.
=3Dcut
sub _is_alphabet_strict {
my ($seqobj) =3D @_;
my $moltype =3D $seqobj->moltype();
# convert everything to upper case to be safe
my $seqstring =3D uc $seqobj->seq(); =20
# First we check if only the 'strict' letters are present in the
# sequence string If not, we check whether the remaining letters are
# ambiguous monomers or whether there are illegal letters in the
# string
# $alpha_array is a ref to an array of the 'strictly' allowed letters
my $alpha_array =3D $Alphabets_strict{$moltype} ;
=20
# $alphabet contains the allowed letters in string form
my $alphabet =3D join ('', @$alpha_array) ;
=20
unless ($seqstring =3D~ /[^$alphabet]/) {
return 1 ;
}
# Next try to match with the alphabet's ambiguous letters
=20
$alpha_array =3D $Alphabets{$moltype} ;
$alphabet =3D join ('', @$alpha_array) ;
=20
unless ($seqstring =3D~ /[^$alphabet]/) {
return 0 ;
}
=20
# If we got here there is an illegal letter in the sequence
=20
$seqobj->throw(" Alphabet not OK for $seqobj \n");
=20
}
=3Dhead2 _print_data
Title : _print_data
Usage : $seqobj->_print_data() or =
Bio::Tools::SeqStats->_print_data();
Function: Displays dna / rna parameters (used for debugging)
Returns : 1
Args : None
Used for debugging.
=3Dcut
sub _print_data {
print "\n adenine =3D : $adenine \n";
print "\n guanine =3D : $guanine \n";
print "\n cytosine =3D : $cytosine \n";
print "\n thymine =3D : $thymine \n";
print "\n uracil =3D : $uracil \n";
print "\n dna_A_wt =3D : $dna_A_wt \n";
print "\n dna_C_wt =3D : $dna_C_wt \n";
print "\n dna_G_wt =3D : $dna_G_wt \n";
print "\n dna_T_wt =3D : $dna_T_wt \n";
print "\n rna_A_wt =3D : $rna_A_wt \n";
print "\n rna_C_wt =3D : $rna_C_wt \n";
print "\n rna_G_wt =3D : $rna_G_wt \n";
print "\n rna_U_wt =3D : $rna_U_wt \n";
return 1;
}
------=_NextPart_000_006C_01C0EF46.BBBE6460
Content-Type: application/x-perl;
name="seqstats.pl"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="seqstats.pl"
#!C:\perl\bin\perl
use Bio::PrimarySeq;
use Bio::Tools::SeqStats;
$seqobj =3D Bio::PrimarySeq->new(-seq=3D>'ACTGTGGCGTCAACTG',=20
-moltype=3D>'dna',=20
-id=3D>'test');
$seq_stats =3D Bio::Tools::SeqStats->new(-seq=3D>$seqobj);
# obtain a hash of counts of each type of monomer=20
# (ie amino or nucleic acid)
print "\nMonomer Counts using statistics object\n";
$seq_stats =3D Bio::Tools::SeqStats->new(-seq=3D>$seqobj);
$hash_ref =3D $seq_stats->count_monomers(); # eg for DNA sequence
foreach $base (sort keys %$hash_ref) {
print "Number of bases of type ", $base, "=3D ", =
%$hash_ref->{$base},"\n";
}
# or obtain the count directly without creating a new statistics object
print "\nMonomer Counts without statistics object\n";
$hash_ref =3D Bio::Tools::SeqStats->count_monomers($seqobj);
foreach $base (sort keys %$hash_ref) {
print "Number of bases of type ", $base, "=3D ", =
%$hash_ref->{$base},"\n";
}
# obtain hash of counts of each type of codon in a nucleic acid sequence
print "\nCodon Counts usin statistics object\n";
$hash_ref =3D $seq_stats-> count_codons(); # for nucleic acid sequence
foreach $base (sort keys %$hash_ref) {
print "Number of codons of type ", $base, "=3D ", =
%$hash_ref->{$base},"\n";
}
# or
print "\nCodon Counts without statistics object\n";
$hash_ref =3D Bio::Tools::SeqStats->count_codons($seqobj);
foreach $base (sort keys %$hash_ref) {
print "Number of codons of type ", $base, "=3D ", =
%$hash_ref->{$base},"\n";
}
# Obtain the molecular weight of a sequence. Since the sequence may =
contain # ambiguous monomers, the molecular weight is returned as a =
(reference to) a # two element array containing greatest lower bound =
(GLB) and least upper bound # (LUB) of the molecular weight=20
$weight =3D $seq_stats->get_mol_wt();
print "\nMolecular weight (using statistics object) of sequence ", =
$seqobj->id(),=20
" is between ", $$weight[0], " and " ,=20
$$weight[1], "\n";
# or
$weight =3D Bio::Tools::SeqStats->get_mol_wt($seqobj);
print "\nMolecular weight (without statistics object) of sequence ", =
$seqobj->id(),=20
" is between ", $$weight[0], " and " ,=20
$$weight[1], "\n";
# use a custom alphabet/weight hash
$seqobj =3D Bio::PrimarySeq->new ( -seq =3D> =
'BACDDCDDCACADDAABBDADABBD',
-id =3D> 'TestFragment-12',
-accession_number =3D> 'Test',
-moltype =3D> 'rna'
);
print "\nMonomer Counts using custom alphabet\n";
my @alphabet =3D qw(A B C D);
$seq_stats =3D =
Bio::Tools::SeqStats->new(-seq=3D>$seqobj,-alphabet=3D>\@alphabet);
$hash_ref =3D $seq_stats->count_monomers(); # eg for DNA sequence
foreach $base (sort keys %$hash_ref) {
print "Number of bases of type ", $base, "=3D ", =
%$hash_ref->{$base},"\n";
}
# Obtain the molecular weight of a sequence with a custom alphabet
my %weights =3D =
('A'=3D>[9,10],'B'=3D>[8,10],'C'=3D>[11,11],'D'=3D>[12,13]);
$seq_stats =3D =
Bio::Tools::SeqStats->new(-seq=3D>$seqobj,-weights=3D>\%weights);
$weight =3D $seq_stats->get_mol_wt();
print "\nMolecular weight (using custom weights) of sequence ", =
$seqobj->id(),=20
" is between ", $$weight[0], " and " ,=20
$$weight[1], "\n";
exit(0);
------=_NextPart_000_006C_01C0EF46.BBBE6460--