Bio::A...tistics(3pm)

📦
version 1.7.8 · 3 versions
, in debian bookworm, debian trixie, ubuntu jammy, ubuntu noble

🎀 Bio::Align::DNAStatistics 🎀

Okay so this one is from the BioPerl user-contributed Perl documentation, the 3pm section — like, Bio::A...tistics(3pm) if your terminal is being dramatic and truncating. It's the whole page, and bestie, it's a lot. Sit down, hydrate, we're doing this. 💖

✨ NAME ✨

Bio::Align::DNAStatistics — it calculates a bunch of statistics for a DNA alignment. That's the whole job description. 🌸

✨ SYNOPSIS ✨

Here's the vibe in code form, exactly as it is:

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 ✨

So this object is full of routines that crunch numbers for DNA alignments — statistics, distances, the works. BUT — and this is the part where I hold your hand — the routines are not well tested and they do contain errors at this point. Work is underway to fix them, but like… do not expect this code to hand you the right answer right now, fr fr. 😭 If you actually want correct distances, go use dnadist/distmat in the PHLYIP or EMBOSS packages instead. Periodt.

Several different distance method calculations are supported. The pattern that each one matches is listed in brackets:

  • JukesCantor — [jc|jukes|jukescantor|jukes-cantor]
  • Uncorrected — [jcuncor|uncorrected]
  • F81 — [f81|felsenstein]
  • Kimura — [k2|k2p|k80|kimura]
  • Tamura — [t92|tamura|tamura92]
  • F84 — [f84|felsenstein84]
  • TajimaNei — [tajimanei|tajima\-nei]
  • JinNei — [jinnei|jin\-nei] (not implemented) 🙈

There are also three methods for calculating the ratio of synonymous to non-synonymous mutations. All three are implementations of the Nei-Gojobori evolutionary pathway method, and they use the Jukes-Cantor method of nucleotide substitution. It works well so long as the nucleotide frequencies are roughly equal and there's no significant transition/transversion bias. ✨

But to use these methods, your alignment has to meet some pre-requisites, and bestie, there are six of them:

  1. The DNA alignment must be based on a 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. The alignment must be solely of coding region and be in reading frame 0 to achieve meaningful results.
  4. The 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 same 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.

And 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 they reproduce those results. 👑 If people like having this sort of analysis in BioPerl, other methods for estimating Ds and Dn can be provided later. No pressure though.

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. 🌈

✨ REFERENCES ✨

  • 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.

  • D_Tamura

K Tamura, Mol. Biol. Evol. 1992, 9, 678.

  • D_Kimura

M Kimura, J. Mol. Evol., 1980, 16, 111.

  • JinNei

Jin and Nei, Mol. Biol. Evol. 82, 7, 1990.

  • D_TajimaNei

Tajima and Nei, Mol. Biol. Evol. 1984, 1, 269.

✨ 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/wiki/Mailing_lists  - About the mailing lists

🩷 Support

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 responsive 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 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 the web:

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

✨ AUTHOR - Jason Stajich ✨

Email jason-AT-bioperl.org 💌

✨ CONTRIBUTORS ✨

Richard Adams, richard.adams@ed.ac.uk 🌸

✨ APPENDIX ✨

okay bestie, this is the part where we go method by method, one by one, no skipping 💅 the rest of the docs break down every single object method for you. and heads up: the internal ones usually start with an underscore _, so if you see that leading _ you already know it's giving private, stay-in-your-lane energy 🙈

new

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

new object, who dis 💖 no args, no drama, just vibes.

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: 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

this is the one that does the whole matrix, every pairwise distance in your alignment, all at once. iconic behavior, honestly 👑

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

basically the menu 🍓 — ask it what distance methods exist and it hands you back an array of strings. no args needed.

🔥 D - distance methods 🔥

okay so now the actual distance models, all the D's, they're all giving main character energy. let's go 👑

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: Bio::Matrix::PhylipDist
  • Args:
    • Bio::Align::AlignI of DNA sequences
    • double - gap penalty

one parameter, simple, clean, JC model. the classic ✨

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. Relaxes the assumption of equal base frequencies that is in JC.
  • Returns: Bio::Matrix::PhylipDist
  • Args: Bio::Align::AlignI of DNA sequences

Felsenstein 1981 said "not me assuming every base is equally likely 🙄" and relaxed that JC constraint. we love growth.

D_Uncorrected

  • 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

no correction at all, it's just raw like that 💀 and btw in rare cases where the sequences may not overlap, 'NA' gets subbed in for the distance, so watch out for that.

D_Kimura

  • 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

Kimura 2 parameter model, and yes it does all the pairs, every last one 🫶

D_Kimura_variance

  • 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

same Kimura 2 parameter model but it also hands you the variance. two matrices, one array — first slot is the Kimura distance, second slot is a matrix of variance V(K). double the slay, double the returns 💖

D_Tamura

  • 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

Tamura 1992, that's the model, that's the vibe. (the original usage line literally says "Calculates D..." — we keep it fr fr, no cap.)

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: Bio::Matrix::PhylipDist
  • Args:
    • Bio::Align::AlignI of DNA sequences
    • [optional] double - gap penalty

Felsenstein 1984, and there's an optional gap penalty you can throw in if you're feeling fancy.

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: Bio::Matrix::PhylipDist
  • Args: Bio::Align::AlignI of DNA sequences

TajimaNei 1984, carrying the whole pairwise distance situation on its back 🎀

D_JinNei

  • 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

Jin-Nei 1990, the finale of the D squad 🌈

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

counts up the transversions between two sequences. returns you a plain old integer, no matrix, no fuss. (yes the usage line says transversion not transversions, that's just how it is bestie, we don't fix the source 🙈)

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

and this one counts the transitions in your alignment. integer out, done, periodt 💅

📊 Data Methods 📊

moving on to the data-flavored methods, they're giving very "numbers go brrr" 🧚

pairwise_stats

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

(yeah the Function line is literally empty in the source, we leave it exactly as-is, no inventing facts up in here 🫶) you can pass a new value if you want, it's optional, otherwise it just returns the pairwise_stats value.

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.

Nei-Gojobori stats for a single pairwise comparison. give it an alignment object plus two sequence names, get back a hashref of stats. (the spelling Useage is straight from the source, we keep it verbatim ✨)

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.

every single pairwise combination, all of them, no sequence left behind 😭 you get an arrayref of hashrefs of stats.

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.

average stats across the whole alignment, and you can set the bootstrap iterations — default is 1000 if you don't say anything. bootstrapping is the moment 💅

get_syn_changes

  • 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. 1 — synonymous
  2. 0 — non-syn
  3. -1 — either codon is a stop codon
  • Args: none

this one builds you a symmetric matrix of every codon pair that differs by exactly 1, and tells you what kind of change it is. 1 is synonymous, 0 is non-syn, and -1 means either codon is a stop codon, so tread carefully around those 🚨

dnds_pattern_number

  • 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.

counts the codons in your MSA that have zero gaps. in PAML notation those are the "patterns" 🌸

and that's the whole page bestie, we made it to the end 🎀 this was pulled from perl v5.32.1, dated 2021-08-15, and the page itself is Bio::A...tistics(3pm). go forth and compute your distances, you absolute icon 👑✨

📦 where this page comes from