Bio Sequence Properties

Version Compatibility

Reference examples tested with: BioPython 1.83+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Sequence Properties

Calculate physical and chemical properties of nucleotide and protein sequences using Biopython.

"Calculate GC content" -> Compute the fraction of G+C bases in a nucleotide sequence.

• Python: gc_fraction(seq) (Bio.SeqUtils) - returns a FRACTION 0-1, multiply by 100 for percent.

"Compute a primer melting temperature" -> Estimate Tm for hybridization or PCR.

• Python: MeltingTemp.Tm_NN(seq) (Bio.SeqUtils) - nearest-neighbor, the accurate method for primers.

"Analyze protein properties" -> Compute MW, pI, stability, hydrophobicity from an amino-acid sequence.

• Python: ProteinAnalysis(str_seq) (Bio.SeqUtils.ProtParam).

The governing principle

Most of these functions return a plausible number for any input, so the danger is silent wrongness, not crashes. Three defaults bite hardest: gc_fraction returns a FRACTION (not the percent the legacy GC() returned), molecular_weight defaults to a SINGLE strand (~half a duplex), and Tm_Wallace/Tm_GC are composition-only methods that are wrong for real primers. Pick the function to match the question and verify its units, not just that it ran.

Required Imports

from Bio.Seq import Seq
from Bio.SeqUtils import gc_fraction, molecular_weight, GC123, GC_skew, MeltingTemp, nt_search, seq1, seq3
from Bio.SeqUtils.ProtParam import ProteinAnalysis

DNA/RNA Properties

GC Content

gc_fraction() returns a fraction in [0, 1]. The legacy GC() returned a percent in [0, 100] and was REMOVED in BioPython 1.82 (from Bio.SeqUtils import GC now raises ImportError).

from Bio.SeqUtils import gc_fraction

seq = Seq('ATGCGATCGATCGATCGATCG')
gc = gc_fraction(seq)        # 0.476... (FRACTION, not percent)
gc_percent = gc * 100        # 47.6 - multiply for percent

Factor-of-100 trap: porting GC(seq) to gc_fraction(seq) without * 100 silently underreports 100x, so downstream filters like "GC > 40" reject everything.

Ambiguity-default trap: the new default ambiguous='remove' strips ambiguity codes before computing, but legacy GC() counted them in the length only, which equals the new ambiguous='ignore'. The faithful drop-in replacement is gc_fraction(seq, ambiguous='ignore') * 100. The modes only diverge on sequences that actually contain ambiguity codes (so clean test fixtures hide the difference, real data exposes it).

gc_fraction(seq, ambiguous='remove')    # default: ambiguity codes stripped (neither numerator nor denominator)
gc_fraction(seq, ambiguous='ignore')    # counts ambiguous in denominator only - matches legacy GC()
gc_fraction(seq, ambiguous='weighted')  # each code contributes its mean GC probability (N/X = 0.5)

GC at Codon Positions (GC123)

GC123() returns FOUR PERCENTAGES (0-100) - total GC plus GC at codon positions 1, 2, 3 (position 3 is the wobble base, most free to vary under codon bias). Note the unit inconsistency with gc_fraction: these are percentages, not fractions. GC123 does not handle ambiguity codes.

from Bio.SeqUtils import GC123

gc_total, gc_pos1, gc_pos2, gc_pos3 = GC123(Seq('ATGCGATCGATCGATCGATCG'))  # all 0-100

GC Skew

GC_skew(seq, window=100) returns (G-C)/(G+C) for each non-overlapping window. A window with no G or C returns 0 (the zero-division is guarded), which can be misread as "no skew" rather than "no data".

from Bio.SeqUtils import GC_skew

skew_values = GC_skew(seq, window=1000)  # list of per-window skew values

Biology: cumulative GC skew has a global MINIMUM at the replication origin (oriC) and a MAXIMUM at the terminus on circular bacterial chromosomes - the leading strand is G-enriched from strand-asymmetric mutation/repair. This is the basis of in-silico origin prediction. Compute the cumulative skew by taking the running sum of GC_skew().

xGC_skew() is a GRAPHICS routine (its docstring literally says "GRAPHICS !!!") that draws on a Tkinter canvas. It raises a loud TclError/ImportError in headless environments. For headless cumulative-skew analysis, sum GC_skew() yourself instead.

Molecular Weight

molecular_weight(seq, seq_type='DNA', double_stranded=False, circular=False, monoisotopic=False).

from Bio.SeqUtils import molecular_weight

dna = Seq('ATGCGATCG')
mw_ss = molecular_weight(dna)                          # single-stranded (DEFAULT)
mw_ds = molecular_weight(dna, double_stranded=True)    # full duplex mass
mw_circ = molecular_weight(dna, circular=True)         # no terminal phosphate adjustment
mw_rna = molecular_weight(Seq('AUGCGAUCG'), seq_type='RNA')
mw_prot = molecular_weight(Seq('MRCRS'), seq_type='protein')
mw_mono = molecular_weight(dna, monoisotopic=True)     # most-abundant isotope, for high-res MS

Double-stranded trap (~2x, silent): the default double_stranded=False returns a single-strand mass - roughly half a duplex. For genomic dsDNA pass double_stranded=True. It is NOT exactly half, because the complementary strand has a different base composition, so a non-self-complementary single-strand number cannot be "fixed later" by doubling. ng-to-molecule, copy-number, and molarity conversions all come out ~2x wrong.

Monoisotopic vs average: the default is average mass (bulk, spectrophotometric). Pass monoisotopic=True to match high-resolution mass spec (ESI/MALDI); the wrong choice is a silent systematic offset that grows with mass. Ambiguous letters raise ValueError (loud - good).

Melting Temperature

Goal: Estimate the Tm of an oligo, choosing a method that matches its length and use.

Approach: Use Tm_NN for any PCR primer (nearest-neighbor, sequence-order aware). Reserve Tm_Wallace for very short probes and Tm_GC only when a composition-only estimate is acceptable.

from Bio.SeqUtils import MeltingTemp as mt

primer = Seq('ACGGTCAGGTCAGGTACGGT')

tm = mt.Tm_NN(primer, strict=True)                       # accurate primer Tm
tm_salt = mt.Tm_NN(primer, Na=50, dnac1=250, dnac2=250)  # 50 mM Na+, 250 nM each strand
tm_mg = mt.Tm_NN(primer, Mg=1.5, dNTPs=0.2, saltcorr=7)  # Mg/dNTPs ONLY honored at saltcorr 6 or 7

Tm_NN defaults: nn_table=None which selects DNA_NN3 (Allawi & SantaLucia 1997), saltcorr=5, strand concentrations dnac1=dnac2=25 nM. saltcorr ranges 1-7, but only 6 (Owczarzy 2004) and 7 (Owczarzy 2008) actually use Mg2+/dNTPs - setting Mg/dNTPs with saltcorr<=5 silently ignores them. Keep strict=True for primer work: it raises on ambiguous or unsupported nearest-neighbor pairs, whereas strict=False silently skips them and underestimates Tm.

IUPAC-Aware Search (nt_search)

from Bio.SeqUtils import nt_search

result = nt_search('ATGCGATCGATCGATNGATC', 'GATNGATC')  # ['GAT[GATC]GATC', 4] - result[0] is the EXPANDED regex, result[1:] are 0-based starts

Protein Properties

Goal: Compute biophysical properties of a protein from its amino-acid sequence.

Approach: Create one ProteinAnalysis object and call its methods. Non-standard residues (B, Z, X, U, *, -) are absent from the parameter tables and raise KeyError, so sanitize first.

from Bio.SeqUtils.ProtParam import ProteinAnalysis

clean = 'MAEGEITTFTALTEKFNLPPGNYKKPKLLYCSNG'.replace('*', '').replace('X', '')
protein = ProteinAnalysis(clean)

mw = protein.molecular_weight()              # protein average MW (Daltons)
pi = protein.isoelectric_point()             # pI from linear-sequence pKa tables
charge = protein.charge_at_pH(7.0)           # net charge at a given pH
ii = protein.instability_index()             # Guruprasad: > 40 => predicted unstable
gravy = protein.gravy()                       # mean Kyte-Doolittle hydropathy (neg = hydrophilic)
arom = protein.aromaticity()                 # relative frequency of F + W + Y
helix, turn, sheet = protein.secondary_structure_fraction()
eps_reduced, eps_oxidized = protein.molar_extinction_coefficient()  # 280 nm, (reduced, cystine)
flex = protein.flexibility()                 # per-residue, fixed window of 9

Interpretation caveats:

• isoelectric_point() and charge_at_pH() use fixed pKa tables on the LINEAR sequence - they ignore 3D environment and post-translational modifications, so a measured pI can differ by a full pH unit or more.
• gravy(scale='KyteDoolitle') - the default scale literal is MISSPELLED 'KyteDoolitle' (one 't'). Passing the correctly-spelled 'KyteDoolittle' raises KeyError. A single whole-protein average also collapses local topology (a TM helix plus hydrophilic loops can average near 0); use a windowed hydropathy profile for membrane topology.
• secondary_structure_fraction() is a composition propensity estimate, not a structure prediction.
• instability_index() uses Guruprasad's dipeptide method; > 40 predicts an unstable protein.