# build a primary nucleic acid or protein sequence object somehow
# then build a statistics object from the sequence object
$seqobj = Bio::PrimarySeq->new(-seq=>'ACTGTGGCGTCAACTG',
-alphabet=>'dna',
-id=>'test');
$seq_stats = Bio::Tools::SeqStats->new(-seq=>$seqobj);
# obtain a hash of counts of each type of monomer
# (ie amino or nucleic acid)
print "\nMonomer counts using statistics object\n";
$seq_stats = Bio::Tools::SeqStats->new(-seq=>$seqobj);
$hash_ref = $seq_stats->count_monomers(); # eg for DNA sequence
foreach $base (sort keys %$hash_ref) {
print "Number of bases of type ", $base, "= ", %$hash_ref->{$base},"\n";
}
# or obtain the count directly without creating a new statistics object
print "\nMonomer counts without statistics object\n";
$hash_ref = Bio::Tools::SeqStats->count_monomers($seqobj);
foreach $base (sort keys %$hash_ref) {
print "Number of bases of type ", $base, "= ", %$hash_ref->{$base},"\n";
}
# obtain hash of counts of each type of codon in a nucleic acid sequence
print "\nCodon counts using statistics object\n";
$hash_ref = $seq_stats-> count_codons(); # for nucleic acid sequence
foreach $base (sort keys %$hash_ref) {
print "Number of codons of type ", $base, "= ", %$hash_ref->{$base},"\n";
}
# or
print "\nCodon counts without statistics object\n";
$hash_ref = Bio::Tools::SeqStats->count_codons($seqobj);
foreach $base (sort keys %$hash_ref) {
print "Number of codons of type ", $base, "= ", %$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
$weight = $seq_stats->get_mol_wt();
print "\nMolecular weight (using statistics object) of sequence ", $seqobj->id(),
" is between ", $$weight[0], " and " ,
$$weight[1], "\n";
# or
$weight = Bio::Tools::SeqStats->get_mol_wt($seqobj);
print "\nMolecular weight (without statistics object) of sequence ", $seqobj->id(),
" is between ", $$weight[0], " and " ,
$$weight[1], "\n";
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.
Nota that nucleotide sequences in bioperl do not strictly separate RNA
and DNA sequences. By convension, sequences from RNA molecules are
shown as is they were DNA. Objects are supposed to make the
distinction when needed. This class is one of the few where this
distinctions needs to be made. Internally, it changes all Ts into Us
before weight and monomer count.
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 = 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 = Bio::Tools::SeqStats->new($seqobj);
$monomers = $seq_stats->count_monomers();
$codons = $seq_stats->count_codons();
$weight = $seq_stats->get_mol_wt();
As currently implemented the object can return the following values
from a sequence:
*(1)
The molecular weight of the sequence: get_mol_wt()
*(2)
The number of each type of monomer present: count_monomers()
*(3)
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.
BEGIN { %Alphabets = (
'dna' => [ qw(A C G T R Y M K S W H B V D X N) ],
'rna' => [ qw(A C G U R Y M K S W H B V D X N) ],
'protein' => [ qw(A R N D C Q E G H I L K M F
P S T W X Y V B Z *) ], );
%Alphabets_strict = (
'dna' => [ qw( A C G T ) ],
'rna' => [ qw( A C G U ) ],
'protein' => [ qw(A R N D C Q E G H I L K M F
P S T W Y V) ],
);
my $amino_A_wt = 89.09;
my $amino_C_wt = 121.15;
my $amino_D_wt = 133.1;
my $amino_E_wt = 147.13;
my $amino_F_wt = 165.19;
my $amino_G_wt = 75.07;
my $amino_H_wt = 155.16;
my $amino_I_wt = 131.18;
my $amino_K_wt = 146.19;
my $amino_L_wt = 131.18;
my $amino_M_wt = 149.22;
my $amino_N_wt = 132.12;
my $amino_P_wt = 115.13;
my $amino_Q_wt = 146.15;
my $amino_R_wt = 174.21;
my $amino_S_wt = 105.09;
my $amino_T_wt = 119.12;
my $amino_V_wt = 117.15;
my $amino_W_wt = 204.22;
my $amino_Y_wt = 181.19;
$amino_weights = {
'A' => [$amino_A_wt, $amino_A_wt], 'B' => [$amino_N_wt, $amino_D_wt], 'C' => [$amino_C_wt, $amino_C_wt], 'D' => [$amino_D_wt, $amino_D_wt], 'E' => [$amino_E_wt, $amino_E_wt], 'F' => [$amino_F_wt, $amino_F_wt], 'G' => [$amino_G_wt, $amino_G_wt], 'H' => [$amino_H_wt, $amino_H_wt], 'I' => [$amino_I_wt, $amino_I_wt], 'K' => [$amino_K_wt, $amino_K_wt], 'L' => [$amino_L_wt, $amino_L_wt], 'M' => [$amino_M_wt, $amino_M_wt], 'N' => [$amino_N_wt, $amino_N_wt], 'P' => [$amino_P_wt, $amino_P_wt], 'Q' => [$amino_Q_wt, $amino_Q_wt], 'R' => [$amino_R_wt, $amino_R_wt], 'S' => [$amino_S_wt, $amino_S_wt], 'T' => [$amino_T_wt, $amino_T_wt], 'V' => [$amino_V_wt, $amino_V_wt], 'W' => [$amino_W_wt, $amino_W_wt], 'X' => [$amino_G_wt, $amino_W_wt], 'Y' => [$amino_Y_wt, $amino_Y_wt], 'Z' => [$amino_Q_wt, $amino_E_wt], };
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
$rna_A_wt $rna_C_wt $rna_G_wt $rna_U_wt));
use vars ( qw($dna_weights $rna_weights %Weights));
$C = 12.01;
$O = 16.00;
$N = 14.01;
$H = 1.01;
$P = 30.97;
$water = 18.015;
$adenine = 5 * $C + 5 * $N + 5 * $H;
$guanine = 5 * $C + 5 * $N + 1 * $O + 5 * $H;
$cytosine = 4 * $C + 3 * $N + 1 * $O + 5 * $H;
$thymine = 5 * $C + 2 * $N + 2 * $O + 6 * $H;
$uracil = 4 * $C + 2 * $N + 2 * $O + 4 * $H;
$ribose_phosphate = 5 * $C + 7 * $O + 9 * $H + 1 * $P; $deoxyribose_phosphate = 5 * $C + 6 * $O + 9 * $H + 1 * $P;
$dna_A_wt = $adenine + $deoxyribose_phosphate - $water;
$dna_C_wt = $cytosine + $deoxyribose_phosphate - $water;
$dna_G_wt = $guanine + $deoxyribose_phosphate - $water;
$dna_T_wt = $thymine + $deoxyribose_phosphate - $water;
$rna_A_wt = $adenine + $ribose_phosphate - $water;
$rna_C_wt = $cytosine + $ribose_phosphate - $water;
$rna_G_wt = $guanine + $ribose_phosphate - $water;
$rna_U_wt = $uracil + $ribose_phosphate - $water;
$dna_weights = {
'A' => [$dna_A_wt,$dna_A_wt], 'C' => [$dna_C_wt,$dna_C_wt], 'G' => [$dna_G_wt,$dna_G_wt], 'T' => [$dna_T_wt,$dna_T_wt], 'M' => [$dna_C_wt,$dna_A_wt], 'R' => [$dna_A_wt,$dna_G_wt], 'W' => [$dna_T_wt,$dna_A_wt], 'S' => [$dna_C_wt,$dna_G_wt], 'Y' => [$dna_C_wt,$dna_T_wt], 'K' => [$dna_T_wt,$dna_G_wt], 'V' => [$dna_C_wt,$dna_G_wt], 'H' => [$dna_C_wt,$dna_A_wt], 'D' => [$dna_T_wt,$dna_G_wt], 'B' => [$dna_C_wt,$dna_G_wt], 'X' => [$dna_C_wt,$dna_G_wt], 'N' => [$dna_C_wt,$dna_G_wt], };
$rna_weights = {
'A' => [$rna_A_wt,$rna_A_wt], 'C' => [$rna_C_wt,$rna_C_wt], 'G' => [$rna_G_wt,$rna_G_wt], 'U' => [$rna_U_wt,$rna_U_wt], 'M' => [$rna_C_wt,$rna_A_wt], 'R' => [$rna_A_wt,$rna_G_wt], 'W' => [$rna_U_wt,$rna_A_wt], 'S' => [$rna_C_wt,$rna_G_wt], 'Y' => [$rna_C_wt,$rna_U_wt], 'K' => [$rna_U_wt,$rna_G_wt], 'V' => [$rna_C_wt,$rna_G_wt], 'H' => [$rna_C_wt,$rna_A_wt], 'D' => [$rna_U_wt,$rna_G_wt], 'B' => [$rna_C_wt,$rna_G_wt], 'X' => [$rna_C_wt,$rna_G_wt], 'N' => [$rna_C_wt,$rna_G_wt], };
%Weights = (
'dna' => $dna_weights,
'rna' => $rna_weights,
'protein' => $amino_weights,
); } |