/ 
bioperl-1.4
/ Bio::Align::DNAStatistics

NAME

Bio::Align::DNAStatistics - Calculate some statistics for a DNA alignment

SYNOPSIS

  use Bio::AlignIO;
  use Bio::Align::DNAStatistics;

  my $stats = new Bio::Align::DNAStatistics;
  my $alignin = new Bio::AlignIO(-format => 'emboss',
                                 -file   => 't/data/insulin.water');
  my $aln = $alignin->next_aln;
  my $jc = $stats->distance(-align => $aln, 
                            -method => 'Jukes-Cantor');
  foreach my $d ( @$jc )  {
      print "\t";
      foreach my $r ( @$d ) {
	  print "$r\t";
      } 
      print "\n";
  }
  ## and for measurements of synonymous /nonsynonymous substitutions ##

  my $in = new Bio::AlignIO(-format => 'fasta',
                            -file   => 't/data/nei_gojobori_test.aln');
  my $alnobj = $in->next_aln;
  my ($seqid,$seq2id) = map { $_->display_id } $alnobj->each_seq;
  my $results = $stats->calc_KaKs_pair($alnobj, $seqid, $seq2id);
  print "comparing ".$results->[0]{'Seq1'}." and ".$results->[0]{'Seq2'}."\n";
  for (sort keys %{$results->[0]} ){
      next if /Seq/;
      printf("%-9s %.4f \n",$_ , $results->[0]{$_});
  }

  my $results2 = $stats->calc_all_KaKs_pairs($alnobj);
  for my $an (@$results2){
      print "comparing ". $an->{'Seq1'}." and ". $an->{'Seq2'}. " \n";
      for (sort keys %$an ){
	  next if /Seq/;
	  printf("%-9s %.4f \n",$_ , $an->{$_});
      }
      print "\n\n";
  }

  my $result3 = $stats->calc_average_KaKs($alnobj, 1000);
  for (sort keys %$result3 ){
      next if /Seq/;
      printf("%-9s %.4f \n",$_ , $result3->{$_});
  }

DESCRIPTION

This object contains routines for calculating various statistics and distances for DNA alignments. The routines are not well tested and do contain errors at this point. Work is underway to correct them, but do not expect this code to give you the right answer currently! Use dnadist/distmat in the PHLYIP or EMBOSS packages to calculate the distances.

There are also three methods to calculate the ratio of synonymous to non-synonymous mutations. All are implementations of the Nei-Gojobori evolutionary pathway method and use the Jukes-Cantor method of nucleotide substitution. This method works well so long as the nucleotide frequencies are roughly equal and there is no significant transition/transversion bias. In order to use these methods there are several pre-requisites for the alignment.

  1. DNA alignment must be based on protein alignment. Use the subroutine aa_to_dna_aln in Bio::Align::Utilities to achieve this.

  2. Therefore alignment gaps must be in multiples of 3 (representing an aa deletion/insertion) and at present must be indicated by a '-' symbol.

  3. Alignment must be solely of coding region and be in reading frame 0 to achieve meaningful results

  4. ALignment must therefore be a multiple of 3 nucleotides long.

  5. All sequences must be the same length (including gaps). This should be the case anyway if the sequences have been automatically aligned using a program like Clustal.

  6. Only the standard codon alphabet is supported at present.

calc_KaKs_pair() calculates a number of statistics for a named pair of sequences in the alignment.

calc_all_KaKs_pairs() calculates these statistics for all pairwise comparisons in an MSA. The statistics returned are:

S_d

Number of synonymous mutations between the 2 sequences.

N_d

Number of non-synonymous mutations between the 2 sequences.

S

Mean number of synonymous sites in both sequences.

N

mean number of synonymous sites in both sequences.

P_s

proportion of synonymous differences in both sequences given by P_s = S_d/S.

P_n

proportion of non-synonymous differences in both sequences given by P_n = S_n/S.

D_s

estimation of synonymous mutations per synonymous site (by Jukes-CAntor).

D_n

estimation of non-synonymous mutations per non-synonymous site (by Jukes-CAntor).

D_n_var

estimation of variance of D_n .

D_s_var

estimation of variance of S_n.

z_value

calculation of z value.Positive value indicates D_n > D_s, negative value indicates D_s > D_n.

calc_average_KaKs

The statistics returned by calc_average_KaKs are:

D_s

Average number of synonymous mutations/synonymous site.

D_n

Average number of non-synonymous mutations/non-synonymous site.

D_s_var

Estimated variance of Ds from bootstrapped alignments.

D_n_var

Estimated variance of Dn from bootstrapped alignments.

z_score

calculation of z value. POsitive value indicates D_n >D_s, negative values vice versa.

The design of the code is based around the explanation of the Nei-Gojobori algorithm in the excellent book "Molecular Evolution and Phylogenetics" by Nei and Kumar, published by Oxford University Press. The methods have been tested using the worked example 4.1 in the book, and reproduce those results. If people like having this sort of analysis in BioPerl other methods for estimating Ds and Dn can be provided later.

FEEDBACK

Mailing Lists

User feedback is an integral part of the evolution of this and other Bioperl modules. Send your comments and suggestions preferably to the Bioperl mailing list. Your participation is much appreciated.

bioperl-l@bioperl.org              - General discussion
http://bioperl.org/MailList.shtml  - About the mailing lists

Reporting Bugs

Report bugs to the Bioperl bug tracking system to help us keep track of the bugs and their resolution. Bug reports can be submitted via email or the web:

bioperl-bugs@bioperl.org
http://bugzilla.bioperl.org/

AUTHOR - Jason Stajich

Email jason@bioperl.org

Describe contact details here

CONTRIBUTORS

Additional contributors names and emails here Richard Adams, richard.adams@ed.ac.uk

APPENDIX

The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _

new

Title   : new
Usage   : my $obj = new Bio::Align::DNAStatistics();
Function: Builds a new Bio::Align::DNAStatistics object 
Returns : Bio::Align::DNAStatistics
Args    : none

distance

 Title   : distance
 Usage   : my $distance_mat = $stats->distance(-align  => $aln, 
		 			       -method => $method);
 Function: Calculates a distance matrix for all pairwise distances of
           sequences in an alignment.
 Returns : Array ref
 Args    : -align  => Bio::Align::AlignI object
           -method => String specifying specific distance method 
                      (implementing class may assume a default)

available_distance_methods

Title   : available_distance_methods
Usage   : my @methods = $stats->available_distance_methods();
Function: Enumerates the possible distance methods
Returns : Array of strings
Args    : none

D - distance methods

D_JukesCantor

Title   : D_JukesCantor
Usage   : my $d = $stat->D_JukesCantor($aln)
Function: Calculates D (pairwise distance) between 2 sequences in an 
          alignment using the Jukes-Cantor 1 parameter model. 
Returns : ArrayRef of all pairwise distances of all sequence pairs in the alignment
Args    : Bio::Align::AlignI of DNA sequences
          double - gap penalty

D_F81

Title   : D_F81
Usage   : my $d = $stat->D_F81($aln)
Function: Calculates D (pairwise distance) between 2 sequences in an 
          alignment using the Felsenstein 1981 distance model. 
Returns : ArrayRef of a 2d array of all pairwise distances in the alignment
Args    : Bio::Align::AlignI of DNA sequences

D_Kimura

Title   : D_Kimura
Usage   : my $d = $stat->D_Kimura($aln)
Function: Calculates D (pairwise distance) between 2 sequences in an 
          alignment using the Kimura 2 parameter model.
Returns : ArrayRef of pairwise distances between all sequences in alignment
Args    : Bio::Align::AlignI of DNA sequences

D_Tamura

Title   : D_Tamura
Usage   :
Function:
Returns : 
Args    :

D_F84

Title   : D_F84
Usage   : my $d = $stat->D_F84($aln)
Function: Calculates D (pairwise distance) between 2 sequences in an 
          alignment using the Felsenstein 1984 distance model. 
Returns : Distance value
Args    : Bio::Align::AlignI of DNA sequences
          double - gap penalty

D_TajimaNei

Title   : D_TajimaNei
Usage   : my $d = $stat->D_TajimaNei($aln)
Function: Calculates D (pairwise distance) between 2 sequences in an 
          alignment using the TajimaNei 1984 distance model. 
Returns : Distance value
Args    : Bio::Align::AlignI of DNA sequences

K - sequence substitution methods

K_JukesCantor

Title   : K_JukesCantor
Usage   : my $k = $stats->K_JukesCantor($aln)
Function: Calculates K - the number of nucleotide substitutions between 
          2 seqs - according to the Jukes-Cantor 1 parameter model
          This only involves the number of changes between two sequences.
Returns : double
Args    : Bio::Align::AlignI

K_TajimaNei

Title   : K_TajimaNei
Usage   : my $k = $stats->K_TajimaNei($aln)
Function: Calculates K - the number of nucleotide substitutions between 
          2 seqs - according to the Kimura 2 parameter model.
          This does not assume equal frequencies among all the nucleotides.
Returns : ArrayRef of 2d matrix which contains pairwise K values for 
          all sequences in the alignment
Args    : Bio::Align::AlignI

transversions

Title   : transversions
Usage   : my $transversions = $stats->transversion($aln);
Function: Calculates the number of transversions between two sequences in 
          an alignment
Returns : integer
Args    : Bio::Align::AlignI

transitions

Title   : transitions
Usage   : my $transitions = Bio::Align::DNAStatistics->transitions($aln);
Function: Calculates the number of transitions in a given DNA alignment
Returns : integer representing the number of transitions
Args    : Bio::Align::AlignI object

Data Methods

pairwise_stats

Title   : pairwise_stats
Usage   : $obj->pairwise_stats($newval)
Function: 
Returns : value of pairwise_stats
Args    : newvalue (optional)

calc_KaKs_pair

Title    : calc_KaKs_pair
Useage   : my $results = $stats->calc_KaKs_pair($alnobj,
           $name1, $name2).
Function : calculates Nei-Gojobori statistics for pairwise 
           comparison.
Args     : A Bio::Align::AlignI compliant object such as a 
           Bio::SimpleAlign object, and 2 sequence name strings.
Returns  : a reference to a hash of statistics with keys as 
           listed in Description.

calc_all_KaKs_pairs

Title    : calc_all_KaKs_pairs
Useage   : my $results2 = $stats->calc_KaKs_pair($alnobj).
Function : Calculates Nei_gojobori statistics for all pairwise
           combinations in sequence.
Arguments: A Bio::Align::ALignI compliant object such as
           a Bio::SimpleAlign object.
Returns  : A reference to an array of hashes of statistics of
           all pairwise comparisons in the alignment.

calc_average_KaKs

Title    : calc_average_KaKs.  
Useage   : my $res= $stats->calc_average_KaKs($alnobj, 1000).
Function : calculates Nei_Gojobori stats for average of all 
           sequences in the alignment.
Args     : A Bio::Align::AlignI compliant object such as a
           Bio::SimpleAlign object, number of bootstrap iterations
           (default 1000).
Returns  : A reference to a hash of statistics as listed in Description.

get_syn_changes

Title   : get_syn_changes
Usage   : Bio::Align::DNAStatitics->get_syn_chnages
Function: Generate a hashref of all pairwise combinations of codns
          differing by 1
Returns : Symetic matrix using hashes
          First key is codon
          and each codon points to a hashref of codons
          the values of which describe type of change.
          my $type = $hash{$codon1}->{$codon2};
          values are :
            1   synonomous
            0   non-syn
           -1   either codon is a stop codon
Args    : none