Bio::Align DNAStatistics
SummaryIncluded librariesPackage variablesSynopsisDescriptionGeneral documentationMethods
Toolbar
WebCvs
Summary
Bio::Align::DNAStatistics - Calculate some statistics for a DNA alignment
Package variables
Privates (from "my" definitions)
@t = split '', "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG"
Included modules
Bio::Align::PairwiseStatistics
Bio::Matrix::PhylipDist
Bio::Tools::IUPAC
Inherit
Bio::Align::StatisticsI Bio::Root::Root
Synopsis
  use Bio::AlignIO;
use Bio::Align::DNAStatistics;
my $stats = Bio::Align::DNAStatistics->new(); my $alignin = Bio::AlignIO->new(-format => 'emboss', -file => 't/data/insulin.water'); my $aln = $alignin->next_aln; my $jcmatrix = $stats->distance(-align => $aln, -method => 'Jukes-Cantor'); print $jcmatrix->print_matrix; ## and for measurements of synonymous /nonsynonymous substitutions ## my $in = Bio::AlignIO->new(-format => 'fasta', -file => 't/data/nei_gojobori_test.aln'); my $alnobj = $in->next_aln; my ($seq1id,$seq2id) = map { $_->display_id } $alnobj->each_seq; my $results = $stats->calc_KaKs_pair($alnobj, $seq1id, $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.
Several different distance method calculations are supported. Listed
in brackets are the pattern which will match
   *(1)
   JukesCantor [jc|jukes|jukescantor|jukes-cantor]
   *(2)
   Uncorrected [jcuncor|uncorrected]
   *(3)
   F81 [f81|felsenstein]
   *(4)
   Kimura [k2|k2p|k80|kimura]
   *(5)
   Tamura [t92|tamura|tamura92]
   *(6)
   F84 [f84|felsenstein84]
   *(7)
   TajimaNei [tajimanei|tajima\-nei]
   *(8)
   JinNei [jinnei|jin\-nei] (not implemented)
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
Bio::Align::Utilities/aa_to_dna_aln 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:
   *(9)
   S_d - Number of synonymous mutations between the 2 sequences.
   *(10)
   N_d - Number of non-synonymous mutations between the 2 sequences.
   *(11)
   S - Mean number of synonymous sites in both sequences.
   *(12)
   N - mean number of synonymous sites in both sequences.
   *(13)
   P_s - proportion of synonymous differences in both sequences given by
P_s = S_d/S.
   *(14)
   P_n - proportion of non-synonymous differences in both sequences given
by P_n = S_n/S.
   *(15)
   D_s - estimation of synonymous mutations per synonymous site (by
Jukes-Cantor).
   *(16)
   D_n - estimation of non-synonymous mutations per non-synonymous site (by
Jukes-Cantor).
   *(17)
   D_n_var - estimation of variance of D_n .
   *(18)
   D_s_var - estimation of variance of S_n.
   *(19)
   z_value - calculation of z value.Positive value indicates D_n > D_s,
negative value indicates D_s > D_n.
The statistics returned by calc_average_KaKs are:
   *(20)
   D_s - Average number of synonymous mutations/synonymous site.
   *(21)
   D_n - Average number of non-synonymous mutations/non-synonymous site.
   *(22)
   D_s_var - Estimated variance of Ds from bootstrapped alignments.
   *(23)
   D_n_var - Estimated variance of Dn from bootstrapped alignments.
   *(24)
   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.
Much of the DNA distance code is based on implementations in EMBOSS
(Rice et al, www.emboss.org) [distmat.c] and PHYLIP (J. Felsenstein et
al) [dnadist.c]. Insight also gained from Eddy, Durbin, Krogh, &
Mitchison.
Methods
BEGIN Code
newDescriptionCode
distanceDescriptionCode
available_distance_methodsDescriptionCode
D_JukesCantorDescriptionCode
D_F81DescriptionCode
D_UncorrectedDescriptionCode
D_KimuraDescriptionCode
D_Kimura_varianceDescriptionCode
D_TamuraDescriptionCode
D_F84DescriptionCode
D_TajimaNeiDescriptionCode
D_JinNeiDescriptionCode
transversionsDescriptionCode
transitionsDescriptionCode
_trans_count_helper
No description
Code
_build_nt_matrix
No description
Code
_check_ambiguity_nucleotide
No description
Code
_check_arg
No description
Code
pairwise_statsDescriptionCode
calc_KaKs_pairDescriptionCode
calc_all_KaKs_pairsDescriptionCode
calc_average_KaKsDescriptionCode
_run_bootstrap
No description
Code
_resample
No description
Code
_get_av_ds_dn
No description
Code
jk
No description
Code
jk_var
No description
Code
analyse_mutations
No description
Code
count_diffs
No description
Code
get_syn_changesDescriptionCode
dnds_pattern_numberDescriptionCode
count_syn_sites
No description
Code
get_syn_sites
No description
Code
_make_codons
No description
Code
get_codons
No description
Code
mean
No description
Code
variance
No description
Code
sampling_variance
No description
Code
Methods description
newcode    nextTop
 Title   : new
Usage : my $obj = Bio::Align::DNAStatistics->new();
Function: Builds a new Bio::Align::DNAStatistics object
Returns : Bio::Align::DNAStatistics
Args : none
distancecodeprevnextTop
 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 : Bio::Matrix::PhylipDist object
Args : -align => Bio::Align::AlignI object
-method => String specifying specific distance method
(implementing class may assume a default)
See also: Bio::Matrix::PhylipDist
available_distance_methodscodeprevnextTop
 Title   : available_distance_methods
Usage : my @methods = $stats->available_distance_methods();
Function: Enumerates the possible distance methods
Returns : Array of strings
Args : none
D_JukesCantorcodeprevnextTop
 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 : Bio::Matrix::PhylipDist
Args : Bio::Align::AlignI of DNA sequences
double - gap penalty
D_F81codeprevnextTop
 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.
Relaxes the assumption of equal base frequencies that is
in JC.
Returns : Bio::Matrix::PhylipDist
Args : Bio::Align::AlignI of DNA sequences
D_UncorrectedcodeprevnextTop
 Title   : D_Uncorrected
Usage : my $d = $stats->D_Uncorrected($aln)
Function: Calculate a distance D, no correction for multiple substitutions
is used. In rare cases where sequences may not overlap, 'NA' is
substituted for the distance.
Returns : Bio::Matrix::PhylipDist
Args : Bio::Align::AlignI (DNA Alignment)
[optional] gap penalty
D_KimuracodeprevnextTop
 Title   : D_Kimura
Usage : my $d = $stat->D_Kimura($aln)
Function: Calculates D (pairwise distance) between all pairs of sequences
in an alignment using the Kimura 2 parameter model.
Returns : Bio::Matrix::PhylipDist
Args : Bio::Align::AlignI of DNA sequences
D_Kimura_variancecodeprevnextTop
 Title   : D_Kimura
Usage : my $d = $stat->D_Kimura_variance($aln)
Function: Calculates D (pairwise distance) between all pairs of sequences
in an alignment using the Kimura 2 parameter model.
Returns : array of 2 Bio::Matrix::PhylipDist,
the first is the Kimura distance and the second is
a matrix of variance V(K)
Args : Bio::Align::AlignI of DNA sequences
D_TamuracodeprevnextTop
 Title   : D_Tamura
Usage : Calculates D (pairwise distance) between 2 sequences in an
alignment using Tamura 1992 distance model.
Returns : Bio::Matrix::PhylipDist
Args : Bio::Align::AlignI of DNA sequences
D_F84codeprevnextTop
 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 : Bio::Matrix::PhylipDist
Args : Bio::Align::AlignI of DNA sequences
[optional] double - gap penalty
D_TajimaNeicodeprevnextTop
 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 : Bio::Matrix::PhylipDist
Args : Bio::Align::AlignI of DNA sequences
D_JinNeicodeprevnextTop
 Title   : D_JinNei
Usage : my $d = $stat->D_JinNei($aln)
Function: Calculates D (pairwise distance) between 2 sequences in an
alignment using the Jin-Nei 1990 distance model.
Returns : Bio::Matrix::PhylipDist
Args : Bio::Align::AlignI of DNA sequences
transversionscodeprevnextTop
 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
transitionscodeprevnextTop
 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
pairwise_statscodeprevnextTop
 Title   : pairwise_stats
Usage : $obj->pairwise_stats($newval)
Function:
Returns : value of pairwise_stats
Args : newvalue (optional)
calc_KaKs_paircodeprevnextTop
 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_pairscodeprevnextTop
 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_KaKscodeprevnextTop
 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_changescodeprevnextTop
 Title   : get_syn_changes
Usage : Bio::Align::DNAStatitics->get_syn_changes
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 synonymous
0 non-syn
-1 either codon is a stop codon
Args : none
dnds_pattern_numbercodeprevnextTop
 Title   : dnds_pattern_number
Usage : my $patterns = $stats->dnds_pattern_number($alnobj);
Function: Counts the number of codons with no gaps in the MSA
Returns : Number of codons with no gaps ('patterns' in PAML notation)
Args : A Bio::Align::AlignI compliant object such as a
Bio::SimpleAlign object.
Methods code
BEGINTop
BEGIN {
    $GapChars = '[\.\-]';
    $GCChhars = '[GCS]';
    @Nucleotides = qw(A G T C);
    $SeqCount = 2;
    $Precision = 5;
    
    # these values come from EMBOSS distmat implementation
%NucleotideIndexes = ( 'A' => 0, 'T' => 1, 'C' => 2, 'G' => 3, 'AT' => 0, 'AC' => 1, 'AG' => 2, 'CT' => 3, 'GT' => 4, 'CG' => 5, # these are wrong now
# 'S' => [ 1, 3],
# 'W' => [ 0, 4],
# 'Y' => [ 2, 3],
# 'R' => [ 0, 1],
# 'M' => [ 0, 3],
# 'K' => [ 1, 2],
# 'B' => [ 1, 2, 3],
# 'H' => [ 0, 2, 3],
# 'V' => [ 0, 1, 3],
# 'D' => [ 0, 1, 2],
); $DefaultGapPenalty = 0; # could put ambiguities here?
%DNAChanges = ( 'Transversions' => { 'A' => [ 'T', 'C'], 'T' => [ 'A', 'G'], 'C' => [ 'A', 'G'], 'G' => [ 'C', 'T'], }, 'Transitions' => { 'A' => [ 'G' ], 'G' => [ 'A' ], 'C' => [ 'T' ], 'T' => [ 'C' ], }, ); %DistanceMethods = ( 'jc|jukes|jukescantor|jukes\-cantor' => 'JukesCantor', 'jcuncor|uncorrected' => 'Uncorrected', 'f81|felsenstein81' => 'F81', 'k2|k2p|k80|kimura' => 'Kimura', 't92|tamura|tamura92' => 'Tamura', 'f84|felsenstein84' => 'F84', 'tajimanei|tajima\-nei' => 'TajimaNei', 'jinnei|jin\-nei' => 'JinNei');
}
newdescriptionprevnextTop
sub new {
     my ($class,@args) = @_;
    my $self = $class->SUPER::new(@args);
    
    $self->pairwise_stats( Bio::Align::PairwiseStatistics->new());

    return $self;
}
distancedescriptionprevnextTop
sub distance {
   my ($self,@args) = @_;
   my ($aln,$method) = $self->_rearrange([qw(ALIGN METHOD)],@args);
   if( ! defined $aln || ! ref ($aln) || ! $aln->isa('Bio::Align::AlignI') ) { 
       $self->throw("Must supply a valid Bio::Align::AlignI for the -align parameter in distance");
   }
   $method ||= 'JukesCantor';
   foreach my $m ( keys %DistanceMethods ) {
       if(defined $m &&  $method =~ /$m/i ) {
	   my $mtd = "D_$DistanceMethods{$m}";
	   return $self->$mtd($aln);
       }
   }
   $self->warn("Unrecognized distance method $method must be one of [".
	       join(',',$self->available_distance_methods())."]");
   return;
}
available_distance_methodsdescriptionprevnextTop
sub available_distance_methods {
   my ($self,@args) = @_;
   return values %DistanceMethods;
}
D_JukesCantordescriptionprevnextTop
sub D_JukesCantor {
   my ($self,$aln,$gappenalty) = @_;
   return 0 unless $self->_check_arg($aln);
   $gappenalty = $DefaultGapPenalty unless defined $gappenalty;
   # ambiguities ignored at this point
my (@seqs,@names,@values,%dist); my $seqct = 0; foreach my $seq ( $aln->each_seq) { push @names, $seq->display_id; push @seqs, uc $seq->seq(); $seqct++; } my $precisionstr = "%.$Precision"."f"; for(my $i = 0; $i < $seqct-1; $i++ ) { # (diagonals) distance is 0 for same sequence
$dist{$names[$i]}->{$names[$i]} = [$i,$i]; $values[$i][$i] = sprintf($precisionstr,0); for( my $j = $i+1; $j < $seqct; $j++ ) { my ($matrix,$pfreq,$gaps) = $self->_build_nt_matrix($seqs[$i], $seqs[$j]); # just want diagonals
my $m = ( $matrix->[0]->[0] + $matrix->[1]->[1] + $matrix->[2]->[2] + $matrix->[3]->[3] ); my $D = 1 - ( $m / ($aln->length - $gaps + ( $gaps * $gappenalty)));
my $d = (- 3 / 4) * log ( 1 - (4 * $D/ 3)); # fwd and rev lookup
$dist{$names[$i]}->{$names[$j]} = [$i,$j]; $dist{$names[$j]}->{$names[$i]} = [$i,$j]; $values[$j][$i] = $values[$i][$j] = sprintf($precisionstr,$d); # (diagonals) distance is 0 for same sequence
$dist{$names[$j]}->{$names[$j]} = [$j,$j]; $values[$j][$j] = sprintf($precisionstr,0); } } return Bio::Matrix::PhylipDist->new(-program => 'bioperl_DNAstats', -matrix =>\% dist, -names =>\@ names, -values =>\@ values);
}
D_F81descriptionprevnextTop
sub D_F81 {
   my ($self,$aln,$gappenalty) = @_;
   return 0 unless $self->_check_arg($aln);
   $gappenalty = $DefaultGapPenalty unless defined $gappenalty;
   # ambiguities ignored at this point
my (@seqs,@names,@values,%dist); my $seqct = 0; foreach my $seq ( $aln->each_seq) { push @names, $seq->display_id;; push @seqs, uc $seq->seq(); $seqct++; } my $precisionstr = "%.$Precision"."f"; for(my $i = 0; $i < $seqct-1; $i++ ) { # (diagonals) distance is 0 for same sequence
$dist{$names[$i]}->{$names[$i]} = [$i,$i]; $values[$i][$i] = sprintf($precisionstr,0); for( my $j = $i+1; $j < $seqct; $j++ ) { my ($matrix,$pfreq,$gaps) = $self->_build_nt_matrix($seqs[$i], $seqs[$j]); # just want diagonals
my $m = ( $matrix->[0]->[0] + $matrix->[1]->[1] + $matrix->[2]->[2] + $matrix->[3]->[3] ); my $D = 1 - ( $m / ($aln->length - $gaps + ( $gaps * $gappenalty)));
my $d = (- 3 / 4) * log ( 1 - (4 * $D/ 3)); # fwd and rev lookup
$dist{$names[$i]}->{$names[$j]} = [$i,$j]; $dist{$names[$j]}->{$names[$i]} = [$i,$j]; $values[$j][$i] = $values[$i][$j] = sprintf($precisionstr,$d); # (diagonals) distance is 0 for same sequence
$dist{$names[$j]}->{$names[$j]} = [$j,$j]; $values[$j][$j] = sprintf($precisionstr,0); } } return Bio::Matrix::PhylipDist->new(-program => 'bioperl_DNAstats', -matrix =>\% dist, -names =>\@ names, -values =>\@ values);
}
D_UncorrecteddescriptionprevnextTop
sub D_Uncorrected {
   my ($self,$aln,$gappenalty) = @_;
   $gappenalty = $DefaultGapPenalty unless defined $gappenalty;
   return 0 unless $self->_check_arg($aln);
   # ambiguities ignored at this point
my (@seqs,@names,@values,%dist); my $seqct = 0; foreach my $seq ( $aln->each_seq) { push @names, $seq->display_id; push @seqs, uc $seq->seq(); $seqct++; } my $precisionstr = "%.$Precision"."f"; my $len = $aln->length; for( my $i = 0; $i < $seqct-1; $i++ ) { # (diagonals) distance is 0 for same sequence
$dist{$names[$i]}->{$names[$i]} = [$i,$i]; $values[$i][$i] = sprintf($precisionstr,0); for( my $j = $i+1; $j < $seqct; $j++ ) { my ($matrix,$pfreq,$gaps) = $self->_build_nt_matrix($seqs[$i], $seqs[$j]); my $m = ( $matrix->[0]->[0] + $matrix->[1]->[1] + $matrix->[2]->[2] + $matrix->[3]->[3] ); my $denom = ( $len - $gaps + ( $gaps * $gappenalty)); $self->warn("No distance calculated between $names[$i] and $names[$j], inserting -1") unless $denom; my $D = $denom ? 1 - ( $m / $denom) : -1;
# fwd and rev lookup
$dist{$names[$i]}->{$names[$j]} = [$i,$j]; $dist{$names[$j]}->{$names[$i]} = [$i,$j]; $values[$j][$i] = $values[$i][$j] = $denom ? sprintf($precisionstr,$D) : sprintf("%-*s", $Precision + 2, $D); # (diagonals) distance is 0 for same sequence
$dist{$names[$j]}->{$names[$j]} = [$j,$j]; $values[$j][$j] = sprintf($precisionstr,0); } } return Bio::Matrix::PhylipDist->new(-program => 'bioperl_DNAstats', -matrix =>\% dist, -names =>\@ names, -values =>\@ values); } # M Kimura, J. Mol. Evol., 1980, 16, 111.
}
D_KimuradescriptionprevnextTop
sub D_Kimura {
   my ($self,$aln) = @_;
   return 0 unless $self->_check_arg($aln);
   # ambiguities ignored at this point
my (@names,@values,%dist); my $seqct = 0; foreach my $seq ( $aln->each_seq) { push @names, $seq->display_id; $seqct++; } my $precisionstr = "%.$Precision"."f"; for( my $i = 0; $i < $seqct-1; $i++ ) { # (diagonals) distance is 0 for same sequence
$dist{$names[$i]}->{$names[$i]} = [$i,$i]; $values[$i][$i] = sprintf($precisionstr,0); for( my $j = $i+1; $j < $seqct; $j++ ) { my $pairwise = $aln->select_noncont($i+1,$j+1); my $L = $self->pairwise_stats->number_of_comparable_bases($pairwise); unless( $L ) { $L = 1; } my $P = $self->transitions($pairwise) / $L;
my $Q = $self->transversions($pairwise) / $L;
my $K = 0; my $denom = ( 1 - (2 * $P) - $Q); if( $denom == 0 ) { $self->throw("cannot find distance for ",$i+1, ",",$j+1," $P, $Q\n"); } my $a = 1 / ( 1 - (2 * $P) - $Q);
my $b = 1 / ( 1 - 2 * $Q );
if( $a < 0 || $b < 0 ) { $K = -1; } else{ $K = (1/2) * log ( $a ) + (1/4) * log($b); } # fwd and rev lookup
$dist{$names[$i]}->{$names[$j]} = [$i,$j]; $dist{$names[$j]}->{$names[$i]} = [$i,$j]; $values[$j][$i] = $values[$i][$j] = sprintf($precisionstr,$K); # (diagonals) distance is 0 for same sequence
$dist{$names[$j]}->{$names[$j]} = [$j,$j]; $values[$j][$j] = sprintf($precisionstr,0); } } return Bio::Matrix::PhylipDist->new(-program => 'bioperl_DNAstats', -matrix =>\% dist, -names =>\@ names, -values =>\@ values);
}
D_Kimura_variancedescriptionprevnextTop
sub D_Kimura_variance {
   my ($self,$aln) = @_;
   return 0 unless $self->_check_arg($aln);
   # ambiguities ignored at this point
my (@names,@values,%dist,@var); my $seqct = 0; foreach my $seq ( $aln->each_seq) { push @names, $seq->display_id; $seqct++; } my $precisionstr = "%.$Precision"."f"; for( my $i = 0; $i < $seqct-1; $i++ ) { # (diagonals) distance is 0 for same sequence
$dist{$names[$i]}->{$names[$i]} = [$i,$i]; $values[$i][$i] = sprintf($precisionstr,0); for( my $j = $i+1; $j < $seqct; $j++ ) { my $pairwise = $aln->select_noncont($i+1,$j+1); my $L = $self->pairwise_stats->number_of_comparable_bases($pairwise); unless( $L ) { $L = 1; } my $P = $self->transitions($pairwise) / $L;
my $Q = $self->transversions($pairwise) / $L;
my ($a,$b,$K,$var_k); my $a_denom = ( 1 - (2 * $P) - $Q); my $b_denom = 1 - 2 * $Q; unless( $a_denom > 0 && $b_denom > 0 ) { $a = 1; $b = 1; $K = -1; $var_k = -1; } else { $a = 1 / $a_denom;
$b = 1 / $b_denom;
$K = (1/2) * log ( $a ) + (1/4) * log($b); # from Wu and Li 1985 which in turn is from Kimura 1980
my $c = ( $a - $b ) / 2;
my $d = ( $a + $b ) / 2;
$var_k = ( $a**2 * $P + $d**2 * $Q - ( $a * $P + $d * $Q)**2 ) / $L;
} # fwd and rev lookup
$dist{$names[$i]}->{$names[$j]} = [$i,$j]; $dist{$names[$j]}->{$names[$i]} = [$i,$j]; $values[$j][$i] = $values[$i][$j] = sprintf($precisionstr,$K); # (diagonals) distance is 0 for same sequence
$dist{$names[$j]}->{$names[$j]} = [$j,$j]; $values[$j]->[$j] = sprintf($precisionstr,0); $var[$j]->[$i] = $var[$i]->[$j] = sprintf($precisionstr,$var_k); $var[$j]->[$j] = $values[$j]->[$j]; } } return ( Bio::Matrix::PhylipDist->new(-program => 'bioperl_DNAstats', -matrix =>\% dist, -names =>\@ names, -values =>\@ values), Bio::Matrix::PhylipDist->new(-program => 'bioperl_DNAstats', -matrix =>\% dist, -names =>\@ names, -values =>\@ var) ); } # K Tamura, Mol. Biol. Evol. 1992, 9, 678.
}
D_TamuradescriptionprevnextTop
sub D_Tamura {
   my ($self,$aln) = @_;
   return 0 unless $self->_check_arg($aln);
   # ambiguities ignored at this point
my (@seqs,@names,@values,%dist,$i,$j); my $seqct = 0; my $length = $aln->length; foreach my $seq ( $aln->each_seq) { push @names, $seq->display_id;; push @seqs, uc $seq->seq(); $seqct++; } my $precisionstr = "%.$Precision"."f"; my (@gap,@gc,@trans,@tranv,@score); $i = 0; for my $t1 ( @seqs ) { $j = 0; for my $t2 ( @seqs ) { $gap[$i][$j] = 0; for( my $k = 0; $k < $length; $k++ ) { my ($c1,$c2) = ( substr($seqs[$i],$k,1), substr($seqs[$j],$k,1) ); if( $c1 =~ /^$GapChars$/ || $c2 =~ /^$GapChars$/ ) { $gap[$i][$j]++; } elsif( $c2 =~ /^$GCChhars$/i ) { $gc[$i][$j]++; } } $gc[$i][$j] = ( $gc[$i][$j] /
(
$length - $gap[$i][$j]) );
$j++; } $i++; } for( $i = 0; $i < $seqct-1; $i++ ) { # (diagonals) distance is 0 for same sequence
$dist{$names[$i]}->{$names[$i]} = [$i,$i]; $values[$i][$i] = sprintf($precisionstr,0); for( $j = $i+1; $j < $seqct; $j++ ) { my $pairwise = $aln->select_noncont($i+1,$j+1); my $L = $self->pairwise_stats->number_of_comparable_bases($pairwise); my $P = $self->transitions($pairwise) / $L;
my $Q = $self->transversions($pairwise) / $L;
my $C = $gc[$i][$j] + $gc[$j][$i]- ( 2 * $gc[$i][$j] * $gc[$j][$i] ); if( $P ) { $P = $P / $C;
} my $d = -($C * log(1- $P - $Q)) -(0.5* ( 1 - $C) * log(1 - 2 * $Q)); # fwd and rev lookup
$dist{$names[$i]}->{$names[$j]} = [$i,$j]; $dist{$names[$j]}->{$names[$i]} = [$i,$j]; $values[$j][$i] = $values[$i][$j] = sprintf($precisionstr,$d); # (diagonals) distance is 0 for same sequence
$dist{$names[$j]}->{$names[$j]} = [$j,$j]; $values[$j][$j] = sprintf($precisionstr,0); } } return Bio::Matrix::PhylipDist->new(-program => 'bioperl_DNAstats', -matrix =>\% dist, -names =>\@ names, -values =>\@ values);
}
D_F84descriptionprevnextTop
sub D_F84 {
   my ($self,$aln,$gappenalty) = @_;
   return 0 unless $self->_check_arg($aln);
   $self->throw_not_implemented();
   # ambiguities ignored at this point
my (@seqs,@names,@values,%dist); my $seqct = 0; foreach my $seq ( $aln->each_seq) { # if there is no name,
my $id = $seq->display_id; if( ! length($id) || # deal with empty names
$id =~ /^\s+$/ ) { $id = $seqct+1; } push @names, $id; push @seqs, uc $seq->seq(); $seqct++; } my $precisionstr = "%.$Precision"."f"; for( my $i = 0; $i < $seqct-1; $i++ ) { # (diagonals) distance is 0 for same sequence
$dist{$names[$i]}->{$names[$i]} = [$i,$i]; $values[$i][$i] = sprintf($precisionstr,0); for( my $j = $i+1; $j < $seqct; $j++ ) { } } } # Tajima and Nei, Mol. Biol. Evol. 1984, 1, 269.
# Tajima-Nei correction used for multiple substitutions in the calc
# of the distance matrix. Nucleic acids only.
#
# D = p-distance = 1 - (matches/(posns_scored + gaps)
#
# distance = -b * ln(1-D/b)
#
}
D_TajimaNeidescriptionprevnextTop
sub D_TajimaNei {
   my ($self,$aln) = @_;
   return 0 unless $self->_check_arg($aln);
   # ambiguities ignored at this point
my (@seqs,@names,@values,%dist); my $seqct = 0; foreach my $seq ( $aln->each_seq) { # if there is no name,
push @names, $seq->display_id; push @seqs, uc $seq->seq(); $seqct++; } my $precisionstr = "%.$Precision"."f"; my ($i,$j,$bs); # pairwise
for( $i =0; $i < $seqct -1; $i++ ) { $dist{$names[$i]}->{$names[$i]} = [$i,$i]; $values[$i][$i] = sprintf($precisionstr,0); for ( $j = $i+1; $j <$seqct;$j++ ) { my ($matrix,$pfreq,$gaps) = $self->_build_nt_matrix($seqs[$i], $seqs[$j]); my $pairwise = $aln->select_noncont($i+1,$j+1); my $slen = $self->pairwise_stats->number_of_comparable_bases($pairwise); my $fij2 = 0; for( $bs = 0; $bs < 4; $bs++ ) { my $fi = 0; map {$fi += $matrix->[$bs]->[$_] } 0..3; my $fj = 0; # summation
map { $fj += $matrix->[$_]->[$bs] } 0..3; my $fij = ( $fi && $fj ) ? ($fi + $fj) /( 2 * $slen) : 0;
$fij2 += $fij**2; } my ($pair,$h) = (0,0); for( $bs = 0; $bs < 3; $bs++ ) { for(my $bs1 = $bs+1; $bs1 <= 3; $bs1++ ) { my $fij = $pfreq->[$pair++] / $slen;
if( $fij ) { my ($ci1,$ci2,$cj1,$cj2) = (0,0,0,0); map { $ci1 += $matrix->[$_]->[$bs] } 0..3; map { $cj1 += $matrix->[$bs]->[$_] } 0..3; map { $ci2 += $matrix->[$_]->[$bs1] } 0..3; map { $cj2 += $matrix->[$bs1]->[$_] } 0..3; if( $fij ) { $h += ( ($fij**2) / 2 ) / ( ( ( $ci1 + $cj1 ) / (2 * $slen) ) *
( (
$ci2 + $cj2 ) / (2 * $slen) ) ); } $self->debug( "slen is $slen h is $h fij = $fij ci1 =$ci1 cj1=$cj1 ci2=$ci2 cj2=$cj2\n"); } } } # just want diagonals which are matches (A matched A, C -> C)
my $m = ( $matrix->[0]->[0] + $matrix->[1]->[1] + $matrix->[2]->[2] + $matrix->[3]->[3] ); my $D = 1 - ( $m / $slen);
my $d; if( $h == 0 ) { $d = -1; } else { my $b = (1 - $fij2 + (($D**2)/$h)) / 2; my $c = 1- $D/ $b;
if( $c < 0 ) { $d = -1; } else { $d = (-1 * $b) * log ( $c); } } # fwd and rev lookup
$dist{$names[$i]}->{$names[$j]} = [$i,$j]; $dist{$names[$j]}->{$names[$i]} = [$i,$j]; $values[$j][$i] = $values[$i][$j] = sprintf($precisionstr,$d); # (diagonals) distance is 0 for same sequence
$dist{$names[$j]}->{$names[$j]} = [$j,$j]; $values[$j][$j] = sprintf($precisionstr,0); } } return Bio::Matrix::PhylipDist->new(-program => 'bioperl_DNAstats', -matrix =>\% dist, -names =>\@ names, -values =>\@ values); } # Jin and Nei, Mol. Biol. Evol. 82, 7, 1990.
}
D_JinNeidescriptionprevnextTop
sub D_JinNei {
   my ($self,@args) = @_;
   $self->warn("JinNei implementation not completed");
   return;
}
transversionsdescriptionprevnextTop
sub transversions {
   my ($self,$aln) = @_;
   return $self->_trans_count_helper($aln, $DNAChanges{'Transversions'});
}
transitionsdescriptionprevnextTop
sub transitions {
   my ($self,$aln) = @_;
   return $self->_trans_count_helper($aln, $DNAChanges{'Transitions'});
}
_trans_count_helperdescriptionprevnextTop
sub _trans_count_helper {
    my ($self,$aln,$type) = @_;
    return 0 unless( $self->_check_arg($aln) );
    if( ! $aln->is_flush ) { $self->throw("must be flush") }
    my (@tcount);
    my ($first,$second) = ( uc $aln->get_seq_by_pos(1)->seq(),
			    uc $aln->get_seq_by_pos(2)->seq() );
    my $alen = $aln->length; 
    for (my $i = 0;$i<$alen; $i++ ) { 
	my ($c1,$c2) = ( substr($first,$i,1),
			 substr($second,$i,1) );
	if( $c1 ne $c2 ) { 
	    foreach my $nt ( @{$type->{$c1}} ) {
		if( $nt eq $c2) {
		   $tcount[$i]++;
	       }
	    }
	}
    }
    my $sum = 0;
    map { if( $_) { $sum += $_} } @tcount;
    return $sum;
}

# this will generate a matrix which records across the row, the number
# of DNA subst
#
}
_build_nt_matrixdescriptionprevnextTop
sub _build_nt_matrix {
    my ($self,$seqa,$seqb) = @_;
    

    my $basect_matrix = [ [ qw(0 0 0 0) ],  # number of bases that match
[ qw(0 0 0 0) ], [ qw(0 0 0 0) ], [ qw(0 0 0 0) ] ]; my $gaps = 0; # number of gaps
my $pfreq = [ qw( 0 0 0 0 0 0)]; # matrix for pair frequency
my $len_a = length($seqa); for( my $i = 0; $i < $len_a; $i++) { my ($ti,$tj) = (substr($seqa,$i,1),substr($seqb,$i,1)); $ti =~ tr/U/T/; $tj =~ tr/U/T/; if( $ti =~ /^$GapChars$/) { $gaps++; next; } if( $tj =~ /^$GapChars$/) { $gaps++; next } my $ti_index = $NucleotideIndexes{$ti}; my $tj_index = $NucleotideIndexes{$tj}; if( ! defined $ti_index ) { $self->warn("ti_index not defined for $ti\n"); next; } $basect_matrix->[$ti_index]->[$tj_index]++; if( $ti ne $tj ) { $pfreq->[$NucleotideIndexes{join('',sort ($ti,$tj))}]++; } } return ($basect_matrix,$pfreq,$gaps);
}
_check_ambiguity_nucleotidedescriptionprevnextTop
sub _check_ambiguity_nucleotide {
    my ($base1,$base2) = @_;
    my %iub = Bio::Tools::IUPAC->iupac_iub();
    my @amb1 = @{ $iub{uc($base1)} };
    my @amb2 = @{ $iub{uc($base2)} };    
    my ($pmatch) = (0);
    for my $amb ( @amb1 ) {
	if( grep { $amb eq $_ } @amb2 ) {
	    $pmatch = 1;
	    last;
	}
    }
    if( $pmatch ) { 
	return (1 / scalar @amb1) * (1 / scalar @amb2);
    } else { 
	return 0;
    }
}
_check_argdescriptionprevnextTop
sub _check_arg {
    my($self,$aln ) = @_;
    if( ! defined $aln || ! $aln->isa('Bio::Align::AlignI') ) {
	$self->warn("Must provide a Bio::Align::AlignI compliant object to Bio::Align::DNAStatistics");
	return 0;
    } elsif( $aln->get_seq_by_pos(1)->alphabet ne 'dna' ) { 
	$self->warn("Must provide a DNA alignment to Bio::Align::DNAStatistics, you provided a " . $aln->get_seq_by_pos(1)->alphabet);
	return 0;
    }
    return 1;
}
pairwise_statsdescriptionprevnextTop
sub pairwise_stats {
   my ($self,$value) = @_;
   if( defined $value) {
      $self->{'_pairwise_stats'} = $value;
    }
    return $self->{'_pairwise_stats'};
}
calc_KaKs_pairdescriptionprevnextTop
sub calc_KaKs_pair {
    my ( $self, $aln, $seq1_id, $seq2_id) = @_;
    $self->throw("Needs 3 arguments - an alignment object, and 2 sequence ids") 
	if @_!= 4;
    $self->throw ("This calculation needs a Bio::Align::AlignI compatible object, not a [ " . ref($aln) . " ]object") unless $aln->isa('Bio::Align::AlignI');
    my @seqs = (
		#{id => $seq1_id, seq =>($aln->each_seq_with_id($seq1_id))[0]->seq},
#{id => $seq2_id, seq =>($aln->each_seq_with_id($seq2_id))[0]->seq}
{id => $seq1_id, seq => uc(($aln->each_seq_with_id($seq1_id))[0]->seq)}, {id => $seq2_id, seq => uc(($aln->each_seq_with_id($seq2_id))[0]->seq)} ) ; if (length($seqs[0]{'seq'}) != length($seqs[1]{'seq'})) { $self->throw(" aligned sequences must be of equal length!"); } my $results = []; $self->_get_av_ds_dn(\@seqs, $results); return $results;
}
calc_all_KaKs_pairsdescriptionprevnextTop
sub calc_all_KaKs_pairs {
#returns a multi_element_array with all pairwise comparisons
my ($self,$aln) = @_; $self->throw ("This calculation needs a Bio::Align::AlignI compatible object, not a [ " . ref($aln) . " ]object") unless $aln->isa('Bio::Align::AlignI'); my @seqs; for my $seq ($aln->each_seq) { push @seqs, {id => $seq->display_id, seq=>$seq->seq}; } my $results ; $results = $self->_get_av_ds_dn(\@seqs, $results); return $results;
}
calc_average_KaKsdescriptionprevnextTop
sub calc_average_KaKs {
#calculates global value for sequences in alignment using bootstrapping
#this is quite slow (~10 seconds per 3 X 200nt seqs);
my ($self, $aln, $bootstrap_rpt) = @_; $bootstrap_rpt ||= 1000; $self->throw ("This calculation needs a Bio::Align::AlignI compatible object, not a [ " . ref($aln) . " ]object") unless $aln->isa('Bio::Align::AlignI'); my @seqs; for my $seq ($aln->each_seq) { push @seqs, {id => $seq->display_id, seq=>$seq->seq}; } my $results ; my ($ds_orig, $dn_orig) = $self->_get_av_ds_dn(\@seqs); #print "ds = $ds_orig, dn = $dn_orig\n";
$results = {D_s => $ds_orig, D_n => $dn_orig}; $self->_run_bootstrap(\@seqs, $results, $bootstrap_rpt); return $results; } ############## primary internal subs for alignment comparisons ########################
}
_run_bootstrapdescriptionprevnextTop
sub _run_bootstrap {
    ### generates sampled sequences, calculates Ds and Dn values,
### then calculates variance of sampled sequences and add results to results hash
###
my ($self,$seq_ref, $results, $bootstrap_rpt) = @_; my @seqs = @$seq_ref; my @btstrp_aoa; # to hold array of array of nucleotides for resampling
my %bootstrap_values = (ds => [], dn =>[]); # to hold list of av values
#1st make alternative array of codons;
my $c = 0; while ($c < length $seqs[0]{'seq'}) { for (0..$#seqs) { push @{$btstrp_aoa[$_]}, substr ($seqs[$_]{'seq'}, $c, 3); } $c+=3; } for (1..$bootstrap_rpt) { my $sampled = _resample (\@btstrp_aoa); my ($ds, $dn) = $self->_get_av_ds_dn ($sampled) ; # is array ref
push @{$bootstrap_values{'ds'}}, $ds; push @{$bootstrap_values{'dn'}}, $dn; } $results->{'D_s_var'} = sampling_variance($bootstrap_values{'ds'}); $results->{'D_n_var'} = sampling_variance($bootstrap_values{'dn'}); $results->{'z_score'} = ($results->{'D_n'} - $results->{'D_s'}) /
sqrt(
$results->{'D_s_var'} + $results->{'D_n_var'} );
#print "bootstrapped var_syn = $results->{'D_s_var'} \n" ;
#print "bootstrapped var_nc = $results->{'D_n_var'} \n";
#print "z is $results->{'z_score'}\n"; ### end of global set up of/perm look up data
}
_resampledescriptionprevnextTop
sub _resample {
    my $ref = shift;
    my $codon_num = scalar (@{$ref->[0]});
    my @altered;
    for (0..$codon_num -1) {	#for each codon
my $rand = int (rand ($codon_num)); for (0..$#$ref) { push @{$altered[$_]}, $ref->[$_][$rand]; } } my @stringed = map {join '', @$_}@altered; my @return; #now out in random name to keep other subs happy
for (@stringed) { push @return, {id=>'1', seq=> $_}; } return\@ return;
}
_get_av_ds_dndescriptionprevnextTop
sub _get_av_ds_dn {
    # takes array of hashes of sequence strings and ids   #
my $self = shift; my $seq_ref = shift; my $result = shift if @_; my @caller = caller(1); my @seqarray = @$seq_ref; my $bootstrap_score_list; #for a multiple alignment considers all pairwise combinations#
my %dsfor_average = (ds => [], dn => []); for (my $i = 0; $i < scalar @seqarray; $i++) { for (my $j = $i +1; $j<scalar @seqarray; $j++ ){ # print "comparing $i and $j\n";
if (length($seqarray[$i]{'seq'}) != length($seqarray[$j]{'seq'})) { $self->warn(" aligned sequences must be of equal length!"); next; } my $syn_site_count = count_syn_sites($seqarray[$i]{'seq'}, $synsites); my $syn_site_count2 = count_syn_sites($seqarray[$j]{'seq'}, $synsites); # print "syn 1 is $syn_site_count , syn2 is $syn_site_count2\n";
my ($syn_count, $non_syn_count, $gap_cnt) = analyse_mutations($seqarray[$i]{'seq'}, $seqarray[$j]{'seq'}); #get averages
my $av_s_site = ($syn_site_count + $syn_site_count2)/2;
my $av_ns_syn_site = length($seqarray[$i]{'seq'}) - $gap_cnt- $av_s_site ; #calculate ps and pn (p54)
my $syn_prop = $syn_count / $av_s_site;
my $nc_prop = $non_syn_count / $av_ns_syn_site ;
#now use jukes/cantor to calculate D_s and D_n, would alter here if needed a different method
my $d_syn = $self->jk($syn_prop); my $d_nc = $self->jk($nc_prop); #JK calculation must succeed for continuation of calculation
#ret_value = -1 if error
next unless $d_nc >=0 && $d_syn >=0; push @{$dsfor_average{'ds'}}, $d_syn; push @{$dsfor_average{'dn'}}, $d_nc; #if not doing bootstrap, calculate the pairwise comparisin stats
if ($caller[3] =~ /calc_KaKs_pair/ || $caller[3] =~ /calc_all_KaKs_pairs/) { #now calculate variances assuming large sample
my $d_syn_var = jk_var($syn_prop, length($seqarray[$i]{'seq'}) - $gap_cnt ); my $d_nc_var = jk_var($nc_prop, length ($seqarray[$i]{'seq'}) - $gap_cnt); #now calculate z_value
#print "d_syn_var is $d_syn_var,and d_nc_var is $d_nc_var\n";
#my $z = ($d_nc - $d_syn) / sqrt($d_syn_var + $d_nc_var);
my $z = ($d_syn_var + $d_nc_var) ? ($d_nc - $d_syn) / sqrt($d_syn_var + $d_nc_var) : 0;
# print "z is $z\n";
push @$result , {S => $av_s_site, N=>$av_ns_syn_site, S_d => $syn_count, N_d =>$non_syn_count, P_s => $syn_prop, P_n=>$nc_prop, D_s => @{$dsfor_average{'ds'}}[-1], D_n => @{$dsfor_average{'dn'}}[-1], D_n_var =>$d_nc_var, D_s_var => $d_syn_var, Seq1 => $seqarray[$i]{'id'}, Seq2 => $seqarray[$j]{'id'}, z_score => $z, }; $self->warn (" number of mutations too small to justify normal test for $seqarray[$i]{'id'} and $seqarray[$j]{'id'}\n- use Fisher's exact, or bootstrap a MSA") if ($syn_count < 10 || $non_syn_count < 10 ) && $self->verbose > -1 ; }#endif
} } #warn of failure if no results hashes are present
#will fail if Jukes Cantor has failed for all pairwise combinations
#$self->warn("calculation failed!") if scalar @$result ==0;
#return results unless bootstrapping
return $result if $caller[3]=~ /calc_all_KaKs/ || $caller[3] =~ /calc_KaKs_pair/; #else if getting average for bootstrap
return( mean ($dsfor_average{'ds'}),mean ($dsfor_average{'dn'})) ;
}
jkdescriptionprevnextTop
sub jk {
    my ($self, $p) = @_;
    if ($p > 0.75) {
	$self->warn( " Jukes Cantor won't  work -too divergent!");
	return -1;
    }
    return -1 * (3/4) * (log(1 - (4/3) * $p));
}

#works for large value of n (50?100?)
}
jk_vardescriptionprevnextTop
sub jk_var {
    my ($p, $n) = @_;
    return (9 * $p * (1 -$p))/(((3 - 4 *$p) **2) * $n);
} # compares 2 sequences to find the number of synonymous/non
# synonymous mutations between them
}
analyse_mutationsdescriptionprevnextTop
sub analyse_mutations {
    my ($seq1, $seq2) = @_;
    my %mutator = ( 2=> {0=>[[1,2],  # codon positions to be altered 
[2,1]], # depend on which is the same
1=>[[0,2], [2,0]], 2=>[[0,1], [1,0]], }, 3=> [ [0,1,2], # all need to be altered
[1,0,2], [0,2,1], [1,2,0], [2,0,1], [2,1,0] ], ); my $TOTAL = 0; # total synonymous changes
my $TOTAL_n = 0; # total non-synonymous changes
my $gap_cnt = 0; my %input; my $seqlen = length($seq1); for (my $j=0; $j< $seqlen; $j+=3) { $input{'cod1'} = substr($seq1, $j,3); $input{'cod2'} = substr($seq2, $j,3); #ignore codon if beeing compared with gaps!
if ($input{'cod1'} =~ /\-/ || $input{'cod2'} =~ /\-/){ $gap_cnt += 3; #just increments once if there is a pair of gaps
next; } my ($diff_cnt, $same) = count_diffs(\%input); #ignore if codons are identical
next if $diff_cnt == 0 ; if ($diff_cnt == 1) { $TOTAL += $synchanges{$input{'cod1'}}{$input{'cod2'}}; $TOTAL_n += 1 - $synchanges{$input{'cod1'}}{$input{'cod2'}}; #print " \nfordiff is 1 , total now $TOTAL, total n now $TOTAL_n\n\n"
} elsif ($diff_cnt ==2) { my $s_cnt = 0; my $n_cnt = 0; my $tot_muts = 4; #will stay 4 unless there are stop codons at intervening point
OUTER:for my $perm (@{$mutator{'2'}{$same}}) { my $altered = $input{'cod1'}; my $prev= $altered; # print "$prev -> (", $t[$CODONS->{$altered}], ")";
for my $mut_i (@$perm) { #index of codon mutated
substr($altered, $mut_i,1) = substr($input{'cod2'}, $mut_i, 1); if ($t[$CODONS->{$altered}] eq '*') { $tot_muts -=2; #print "changes to stop codon!!\n";
next OUTER; } else { $s_cnt += $synchanges{$prev}{$altered}; # print "$altered ->(", $t[$CODONS->{$altered}], ") ";
} $prev = $altered; } # print "\n";
} if ($tot_muts != 0) { $TOTAL += ($s_cnt/($tot_muts/2)); $TOTAL_n += ($tot_muts - $s_cnt)/ ($tot_muts / 2); } } elsif ($diff_cnt ==3 ) { my $s_cnt = 0; my $n_cnt = 0; my $tot_muts = 18; #potential number of mutations
OUTER: for my $perm (@{$mutator{'3'}}) { my $altered = $input{'cod1'}; my $prev= $altered; # print "$prev -> (", $t[$CODONS->{$altered}], ")";
for my $mut_i (@$perm) { #index of codon mutated
substr($altered, $mut_i,1) = substr($input{'cod2'}, $mut_i, 1); if ($t[$CODONS->{$altered}] eq '*') { $tot_muts -=3; # print "changes to stop codon!!\n";
next OUTER; } else { $s_cnt += $synchanges{$prev}{$altered}; # print "$altered ->(", $t[$CODONS->{$altered}], ") ";
} $prev = $altered; } # print "\n";
}#end OUTER loop
#calculate number of synonymous/non synonymous mutations for that codon
# and add to total
if ($tot_muts != 0) { $TOTAL += ($s_cnt / ($tot_muts /3)); $TOTAL_n += 3 - ($s_cnt / ($tot_muts /3)); } } #endif $diffcnt = 3
} #end of sequencetraversal
return ($TOTAL, $TOTAL_n, $gap_cnt);
}
count_diffsdescriptionprevnextTop
sub count_diffs {
    #counts the number of nucleotide differences between 2 codons
# returns this value plus the codon index of which nucleotide is the same when 2
#nucleotides are different. This is so analyse_mutations() knows which nucleotides
# to change.
my $ref = shift; my $cnt = 0; my $same= undef; #just for 2 differences
for (0..2) { if (substr($ref->{'cod1'}, $_,1) ne substr($ref->{'cod2'}, $_, 1)){ $cnt++; } else { $same = $_; } } return ($cnt, $same);
}
get_syn_changesdescriptionprevnextTop
sub get_syn_changes {
#hash of all pairwise combinations of codons differing by 1
# 1 = syn, 0 = non-syn, -1 = stop
my %results; my @codons = _make_codons (); my $arr_len = scalar @codons; for (my $i = 0; $i < $arr_len -1; $i++) { my $cod1 = $codons[$i]; for (my $j = $i +1; $j < $arr_len; $j++) { my $diff_cnt = 0; for my $pos(0..2) { $diff_cnt++ if substr($cod1, $pos, 1) ne substr($codons[$j], $pos, 1); } next if $diff_cnt !=1; #synon change
if($t[$CODONS->{$cod1}] eq $t[$CODONS->{$codons[$j]}]) { $results{$cod1}{$codons[$j]} =1; $results{$codons[$j]}{$cod1} = 1; } #stop codon
elsif ($t[$CODONS->{$cod1}] eq '*' or $t[$CODONS->{$codons[$j]}] eq '*') { $results{$cod1}{$codons[$j]} = -1; $results{$codons[$j]}{$cod1} = -1; } # nc change
else { $results{$cod1}{$codons[$j]} = 0; $results{$codons[$j]}{$cod1} = 0; } } } return %results;
}
dnds_pattern_numberdescriptionprevnextTop
sub dnds_pattern_number {
    my ($self, $aln) = @_;
    return ($aln->remove_gaps->length)/3;
}
count_syn_sitesdescriptionprevnextTop
sub count_syn_sites {
    #counts the number of possible synonymous changes for sequence
my ($seq, $synsite) = @_; __PACKAGE__->throw("not integral number of codons") if length($seq) % 3 != 0; my $S = 0; for (my $i = 0; $i< length($seq); $i+=3) { my $cod = substr($seq, $i, 3); next if $cod =~ /\-/; #deal with alignment gaps
$S += $synsite->{$cod}{'s'}; } #print "S is $S\n";
return $S;
}
get_syn_sitesdescriptionprevnextTop
sub get_syn_sites {
    #sub to generate lookup hash for the number of synonymous changes per codon
my @nucs = qw(T C A G); my %raw_results; for my $i (@nucs) { for my $j (@nucs) { for my $k (@nucs) { # for each possible codon
my $cod = "$i$j$k"; my $aa = $t[$CODONS->{$cod}]; #calculate number of synonymous mutations vs non syn mutations
for my $i (qw(0 1 2)){ my $s = 0; my $n = 3; for my $nuc (qw(A T C G)) { next if substr ($cod, $i,1) eq $nuc; my $test = $cod; substr($test, $i, 1) = $nuc ; if ($t[$CODONS->{$test}] eq $aa) { $s++; } if ($t[$CODONS->{$test}] eq '*') { $n--; } } $raw_results{$cod}[$i] = {'s' => $s , 'n' => $n }; } } #end analysis of single codon
} } #end analysis of all codons
my %final_results; for my $cod (sort keys %raw_results) { my $t = 0; map{$t += ($_->{'s'} /$_->{'n'})} @{$raw_results{$cod}};
$final_results{$cod} = { 's'=>$t, 'n' => 3 -$t}; } return\% final_results;
}
_make_codonsdescriptionprevnextTop
sub _make_codons {
#makes all codon combinations, returns array of them
my @nucs = qw(T C A G); my @codons; for my $i (@nucs) { for my $j (@nucs) { for my $k (@nucs) { push @codons, "$i$j$k"; } } } return @codons;
}
get_codonsdescriptionprevnextTop
sub get_codons {
 #generates codon translation look up table#
my $x = 0; my $CODONS = {}; for my $codon (_make_codons) { $CODONS->{$codon} = $x; $x++; } return $CODONS; } #########stats subs, can go in another module? Here for speed. ###
}
meandescriptionprevnextTop
sub mean {
    my $ref = shift;
    my $el_num = scalar @$ref;
    my $tot = 0;
    map{$tot += $_}@$ref;
    return ($tot/$el_num);
}
variancedescriptionprevnextTop
sub variance {
    my $ref = shift;
    my $mean = mean($ref);
    my $sum_of_squares = 0;
    map{$sum_of_squares += ($_ - $mean) **2}@$ref;
    return $sum_of_squares;
}
sampling_variancedescriptionprevnextTop
sub sampling_variance {
    my $ref = shift;
    return variance($ref) / (scalar @$ref -1);
} 1;
}
General documentation
REFERENCESTop
   *(1)
   D_JukesCantor
   "Phylogenetic Inference", Swoffrod, Olsen, Waddell and Hillis, in
Mol. Systematics, 2nd ed, 1996, Ch 11. Derived from "Evolution of
Protein Molecules", Jukes & Cantor, in Mammalian Prot. Metab., III,
1969, pp. 21-132.
   *(2)
   D_Tamura
   K Tamura, Mol. Biol. Evol. 1992, 9, 678.
   *(3)
   D_Kimura
   M Kimura, J. Mol. Evol., 1980, 16, 111.
   *(4)
   JinNei
   Jin and Nei, Mol. Biol. Evol. 82, 7, 1990.
   *(5)
   D_TajimaNei
   Tajima and Nei, Mol. Biol. Evol. 1984, 1, 269.
FEEDBACKTop
Mailing ListsTop
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/wiki/Mailing_lists - About the mailing lists
Support Top
Please direct usage questions or support issues to the mailing list:
bioperl-l@bioperl.org
rather than to the module maintainer directly. Many experienced and
reponsive experts will be able look at the problem and quickly
address it. Please include a thorough description of the problem
with code and data examples if at all possible.
Reporting BugsTop
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 the
web:
  https://redmine.open-bio.org/projects/bioperl/
AUTHOR - Jason StajichTop
Email jason-AT-bioperl.org
CONTRIBUTORSTop
Richard Adams, richard.adams@ed.ac.uk
APPENDIXTop
The rest of the documentation details each of the object methods.
Internal methods are usually preceded with a _
D - distance methodsTop
Data MethodsTop