Bio::T...qStats(3pm)

💖 Bio::Tools::SeqStats(3pm) 💖

ok so bestie, this is a User Contributed Perl Documentation page for Bio::Tools::SeqStats — a lil Bioperl module that does math on your sequences. it's giving "calculator but for DNA". let's get into it 🧬✨

---

🌸 NAME 🌸

Bio::Tools::SeqStats — an object that holds statistics for one specific sequence. just the one. it's monogamous like that 💅

---

✨ SYNOPSIS ✨

ok so first you build a primary nucleic acid or protein sequence object somehow, and then you build a statistics object from it. like so:

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

want a hash of counts of each type of monomer (that's amino or nucleic acid, iykyk)? bet:

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

orrr you can skip the whole statistics-object thing and just get the count directly, no cap:

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

now codon counts in a nucleic acid sequence, same energy:

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

orrr again:

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

ok now the molecular weight. here's the tea ☕ — since a sequence can have ambiguous monomers in it, the molecular weight comes back as a (reference to a) two element array, holding the greatest lower bound (GLB) and the least upper bound (LUB) of the molecular weight:

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

orrrr:

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

and finally, calculate mean Kyte-Doolittle hydropathicity (aka the "gravy" score, iconic name tbh 🍓):

# 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 🦋

Bio::Tools::SeqStats is a lightweight object for calculating simple statistical and numerical properties of a sequence. and when I say "lightweight", I mean it only handles "primary" sequences, periodt. if you want stats on a sequence feature, your calling script has to build the appropriate primary sequence and hand it to SeqStats yourself. same deal if you want a codon count for a frame-shifted sequence and/or a negative strand sequence — you make that sequence, you pass it to the SeqStats object. it's not gonna do it for you, bestie 💅

not me noticing the typo — ok so, Nota that nucleotide sequences in bioperl do not strictly separate RNA and DNA sequences. by convention, sequences from RNA molecules are shown as if they were DNA. objects are supposed to make the distinction when needed. this class is one of the few where that distinction needs to be made. internally, it changes all Ts into Us before weight and monomer count. slay 🔥

SeqStats can be called in two distinct manners. if you only need a single computation on a given sequence object, just call the method directly on the SeqStats object:

$weight = Bio::Tools::SeqStats->get_mol_wt($seqobj);

but if you're gonna do several computations on the same sequence object, build yourself an "instance" statistics object and reuse it for all 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 account for the additional OH on the phosphate of the 5' residue and the additional H on the sugar ring of the 3' residue. heads up though — this leads to a difference of 18 in calculated molecular weights compared to some other available programs (e.g. Informax VectorNTI). so if your numbers don't match, that's why, don't panic 😭

and 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. so like, it'll tell you. pay attention 🫶

---

🧚 DEVELOPERS NOTES 🧚

Ewan moved it from Bio::SeqStats to Bio::Tools::SeqStats.

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 💌

📬 Mailing Lists

user feedback is an integral part of the evolution of this and other Bioperl modules. send your comments and suggestions preferably to one of the Bioperl mailing lists. your participation is much appreciated, fr fr 💖

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. the more detail the better, bestie ✨

🐛 Reporting Bugs

report bugs to the Bioperl bug tracking system to help us keep track the bugs and their resolution. bug reports can be submitted the web:

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

---

👑 AUTHOR - Peter Schattner 👑

email schattner AT alum.mit.edu

---

🌈 CONTRIBUTOR - Torsten Seemann 🌈

email torsten.seemann AT infotech.monash.edu.au

---

(part 1 of 2 — more coming, don't go anywhere 💜)

🧾 APPENDIX 🧾

okay bestie the rest of this doc is a lil deep dive into every single object method, one by one. internal methods (the sneaky ones) usually start with an underscore _, so if you see that you know it's backstage crew only 💅

---

🍓 count_monomers

• Title — count_monomers
• Usage — $rcount = $seq_stats->count_monomers(); or $rcount = $seq_stats->Bio::Tools::SeqStats->($seqobj);
• Function — counts how many of each type of monomer (that's an amino acid or a base, depending on what you're working with) are chillin' in the sequence. also, heads up: in RNA sequences the Ts get counted as Us. periodt.
• Example — (none given, it's giving mysterious ✨)
• Returns — a reference to a hash, where the keys are the letters of the genetic alphabet being used and the values are how many times that letter shows up in the sequence 💖
• Args — none, or a reference to a sequence object
• Throws — it throws an exception if the type of sequence is unknown (i.e. not amino and not nucleic), or if there's an unknown letter in the alphabet. ambiguous elements? those are allowed, no stress 🫶

---

🌸 get_mol_wt

• Title — get_mol_wt
• Usage — $wt = $seqobj->get_mol_wt() or $wt = Bio::Tools::SeqStats ->get_mol_wt($seqobj);
• Function — calculates the molecular weight of a sequence. and yes, same energy as before: Ts are counted as Us in RNA sequences 🧬
• Example — (blank, we love a suspense moment)
• Returns — a reference to a two element array holding the lower and upper bounds of the molecule's molecular weight. for DNA and RNA sequences, single-stranded weights are what you get. and if the sequence has no ambiguous elements in it, both entries in that array are equal to the molecular weight of the molecule ⭐
• Args — none, or a reference to a sequence object
• Throws — exception if the type of sequence is unknown (i.e. not amino and not nucleic) or if there's an unknown letter in the alphabet. ambiguous elements are allowed 💜

---

🦋 count_codons

• Title — count_codons
• Usage — $rcount = $seqstats->count_codons() or $rcount = Bio::Tools::SeqStats->count_codons($seqobj)
• Function — counts how many of each type of codon there are for a dna or rna sequence, and it starts counting right at the 1st triple of the input sequence. no skipping, no warm-up 👑
• Example — (empty, iconic silence)
• Returns — a reference to a hash where the keys are the codons of the genetic alphabet used and the values are how many times those codons appear in the sequence. and btw, all codons with "ambiguous" bases get counted together, so don't come for me about it 🙈
• Args — none, or a sequence object
• Throws — an exception if the type of sequence is unknown, or if it's protein 💀

---

🎀 hydropathicity

• Title — hydropathicity
• Usage — $gravy = $seqstats->hydropathicity(); or $gravy = Bio::Tools::SeqStats->hydropathicity($seqobj);
• Function — calculates the mean Kyte-Doolittle hydropathicity for a protein sequence. it also goes by the "gravy" score, which honestly? iconic name. if you wanna read the receipts, refer to Kyte J., Doolittle R.F., J. Mol. Biol. 157:105-132(1982) 📚
• Example — (blank, as usual)
• Returns — float
• Args — none, or a reference to a sequence object
• Throws — an exception if the type of sequence is not protein 🌈

---

🧚 _is_alphabet_strict

• Title — _is_alphabet_strict
• Usage — (nothing here, it's an internal girly)
• Function — an internal function that figures out whether there are any ambiguous elements hanging out in the current sequence 🔍
• Example — (nada)
• Returns — 1 if a strict alphabet is being used, 0 if ambiguous elements are present
• Args — (none)
• Throws — an exception if the type of sequence is unknown (i.e. amino or nucleic) or if there's an unknown letter in the alphabet. ambiguous monomers are allowed, bestie 💖

---

💅 _print_data

• Title — _print_data
• Usage — $seqobj->_print_data() or Bio::Tools::SeqStats->_print_data();
• Function — displays dna / rna parameters (this one's for debugging, so it's giving developer main character energy) 🛠️
• Returns — 1
• Args — None

used for debugging. that's the whole vibe.

---

aaaand that's a wrap on the whole page, bestie!! 🎉 this doc came to you from perl v5.32.1 and it was last dated 2021-08-15 — so yeah, she's been around a minute and she still ate 👑💖 now go count some monomers fr fr ✨🦋