NAME

SYNOPSIS

  # 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
  # (i.e. 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();  # e.g. for DNA sequence
  foreach $base (sort keys %$hash_ref) {
      print "Number of bases of type ", $base, "= ", 
         %$hash_ref->{$base},"\n";
  }

  # 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";

  # Calculate mean Kyte-Doolittle hydropathicity (aka "gravy" score)
  my $prot = Bio::PrimarySeq->new(-seq=>'MSFVLVAPDMLATAAADVVQIGSAVSAGS',
                                  -alphabet=>'protein');
  my $gravy = Bio::Tools::SeqStats->hydropathicity($seqobj);
  print "might be hydropathic" if $gravy > 1;

DESCRIPTION

Nota that nucleotide sequences in bioperl do not strictly separate RNA and DNA sequences. By convention, 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();
  $gravy = $seq_stats->hydropathicity();

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()
• The mean hydropathicity ("gravy" score) of a protein: hydropathicity()

For DNA and RNA sequences single-stranded weights are returned. The molecular weights are calculated for neutral, or 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 (e.g. Informax VectorNTI).

Note that since sequences may contain ambiguous monomers (e.g. "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.

DEVELOPERS NOTES

Heikki made tiny adjustments (+/- 0.01 daltons) to amino acid molecular weights to have the output match values in SWISS-PROT.

Torsten added hydropathicity calculation.

FEEDBACK

  bioperl-l@bioperl.org                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

  https://github.com/bioperl/bioperl-live/issues

AUTHOR - Peter Schattner

CONTRIBUTOR - Torsten Seemann

APPENDIX

 Title   : count_monomers
 Usage   : $rcount = $seq_stats->count_monomers();
           or $rcount = $seq_stats->Bio::Tools::SeqStats->($seqobj);
 Function: Counts the number of each type of monomer (amino acid or
                base) in the sequence.
           Ts are counted as Us in RNA sequences.
 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.

 Title   : get_mol_wt
 Usage   : $wt = $seqobj->get_mol_wt() or
           $wt = Bio::Tools::SeqStats ->get_mol_wt($seqobj);
 Function: Calculate molecular weight of sequence
           Ts are counted as Us in RNA sequences.
 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.
 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.

 Title   : count_codons
 Usage   : $rcount = $seqstats->count_codons() or
           $rcount = Bio::Tools::SeqStats->count_codons($seqobj)
 Function: Counts the number of each type of codons for a dna or rna 
           sequence, starting at the 1st triple of the input 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 sequence object
 Throws  : an exception if type of sequence is unknown or protein.

 Title   : hydropathicity
 Usage   : $gravy = $seqstats->hydropathicity(); or
           $gravy = Bio::Tools::SeqStats->hydropathicity($seqobj);

 Function: Calculates the mean Kyte-Doolittle hydropathicity for a
           protein sequence. Also known as the "gravy" score. Refer to 
           Kyte J., Doolittle R.F., J. Mol. Biol. 157:105-132(1982). 
 Example :
 Returns : float 
 Args    : None or reference to sequence object

 Throws  : an exception if type of sequence is not protein.

 Title   :  _is_alphabet_strict
 Usage   :
 Function: internal function to determine whether there are
           any ambiguous elements in the current sequence
 Example :
 Returns : 1 if strict alphabet is being used,
           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.

 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