Bio.SeqUtils

Core Principle: Bio.SeqUtils provides analytical utilities for calculating sequence GC proportions, molecular weights, oligonucleotide melting temperatures ($T_m$), and codon usage metrics.

GC content

from Bio.SeqUtils import gc_fraction, GC123, gc_skew

gc_fraction(seq)                        # 0.0–1.0  ← A FRACTION
gc_fraction(seq, ambiguous="ignore")    # exclude N from the denominator
gc_fraction(seq, ambiguous="weighted")  # weight ambiguous codes by their possibilities

GC123(seq)                # (overall, GC1, GC2, GC3): by codon position
gc_skew(seq, window=100)  # (G-C)/(G+C) in sliding windows

SeqUtils.GC() was removed in Biopython 1.82 and it returned a percentage (0–100), while gc_fraction returns a fraction (0–1).

from Bio.SeqUtils import GC          # ImportError
GC(seq)                               # was 41.5
gc_fraction(seq)                      # is 0.415

A straight substitution introduces a silent 100× error. Any old code you port needs gc_fraction(seq) * 100 to preserve behaviour, or a review of every downstream threshold. See Biopython Deprecations.

GC123 is the codon-position breakdown: GC3 (third codon position) is the standard codon-bias measure, because third positions are largely synonymous and therefore reflect mutational and translational pressure rather than protein constraint.

gc_skew is used to locate bacterial replication origins: the skew switches sign at the origin and terminus.

Molecular weight

from Bio.SeqUtils import molecular_weight

molecular_weight(seq, seq_type="DNA")                      # double-stranded default
molecular_weight(seq, seq_type="DNA", double_stranded=False)
molecular_weight(seq, seq_type="RNA")
molecular_weight(seq, seq_type="protein")
molecular_weight(seq, seq_type="protein", monoisotopic=True)   # for mass spec

monoisotopic=True matters in proteomics: average and monoisotopic masses differ by roughly 0.05% and that is far more than mass-spec resolution.

Melting temperature

from Bio.SeqUtils import MeltingTemp as mt

mt.Tm_Wallace(primer)        # crude 2AT+4GC rule; short oligos only
mt.Tm_GC(primer)             # GC-content-based
mt.Tm_NN(primer)             # nearest-neighbour thermodynamics: use this
mt.Tm_NN(primer, Na=50, Mg=1.5, dNTPs=0.6, dnac1=250, nn_table=mt.DNA_NN4)
mt.chem_correction(tm, DMSO=5)

Tm_NN is the accurate one and the only one worth using for primer design. It accounts for stacking energies, and the salt/Mg²⁺/dNTP corrections matter: a Tm computed with default salt can be several degrees off from your actual PCR buffer.

Tm_Wallace (the "2+4 rule") is a rule of thumb for oligos under ~14 bp and is not appropriate for real design.

Protein properties

from Bio.SeqUtils.ProtParam import ProteinAnalysis

pa = ProteinAnalysis(str(protein_seq))

pa.molecular_weight()
pa.aromaticity()
pa.instability_index()               # >40 suggests an unstable protein in vitro
pa.isoelectric_point()
pa.gravy()                            # grand average of hydropathy
pa.secondary_structure_fraction()     # (helix, turn, sheet)
pa.count_amino_acids()
pa.get_amino_acids_percent()
pa.protein_scale(ProtParamData.kd, window=9)    # Kyte-Doolittle hydropathy profile
pa.molar_extinction_coefficient()     # for A280 quantification

ProteinAnalysis requires a plain str, not a Seq, and raises on non-standard residues (X, *, U). Clean the sequence first.

The extinction coefficient is the practically useful one: it is what you need to convert an A280 reading into a protein concentration.

Codon adaptation

from Bio.SeqUtils.CodonUsage import CodonAdaptationIndex

cai = CodonAdaptationIndex()
cai.generate_index("highly_expressed_genes.fasta")
cai.cai_for_gene(str(gene_seq))

CAI measures how closely a gene's codon usage matches that of highly expressed genes in the same organism: a proxy for translational efficiency, and a standard step in codon optimisation for heterologous expression.

The old CodonUsage and CodonUsageIndices module-level functions were removed in 1.82; the CodonAdaptationIndex class is the replacement.

Sequence utilities

from Bio.SeqUtils import seq1, seq3, nt_search, six_frame_translations

seq1("MetAlaIleVal")          # 'MAIV' : three-letter → one-letter
seq3("MAIV")                  # 'MetAlaIleVal'
nt_search(str(seq), "GAATTC") # find a pattern, expanding ambiguity codes
print(six_frame_translations(seq))    # human-readable six-frame display

seq1/seq3 are genuinely useful when parsing HGVS protein notation, which uses three-letter codes:

# p.Arg175His → R175H
import re
m = re.match(r"p\.([A-Z][a-z]{2})(\d+)([A-Z][a-z]{2})", "p.Arg175His")
short = f"{seq1(m[1])}{m[2]}{seq1(m[3])}"        # 'R175H'

nt_search expands IUPAC ambiguity in the pattern, so searching for GGWCC finds both GGACC and GGTCC.

Bioinformatics examples

# GC content across a set of contigs
df = pd.DataFrame([
    {"id": r.id, "length": len(r), "gc": gc_fraction(r.seq),
     "n_frac": str(r.seq).upper().count("N") / len(r)}
    for r in SeqIO.parse("contigs.fasta", "fasta")
])
sns.scatterplot(data=df, x="length", y="gc", size="n_frac", alpha=0.5)
plt.xscale("log")

A length-vs-GC scatter of an assembly is a standard contamination check: a distinct cloud at a different GC is usually a different organism.

# primer design sanity checks
for name, p in primers.items():
    tm = mt.Tm_NN(p, Na=50, Mg=1.5, dnac1=250)
    gc = gc_fraction(p)
    print(f"{name}: Tm {tm:.1f}°C, GC {gc:.1%}, len {len(p)}")
    assert 55 <= tm <= 65, f"{name} Tm out of range"
    assert 0.4 <= gc <= 0.6, f"{name} GC out of range"

# GC skew to find a bacterial origin of replication
skew = gc_skew(genome.seq, window=10_000)
cum = np.cumsum(skew)
origin = int(np.argmin(cum) * 10_000)

# codon bias per gene
df["gc3"] = [GC123(r.seq)[3] / 100 for r in cds_records]
sns.histplot(df, x="gc3", bins=50)

# protein properties for a construct
pa = ProteinAnalysis(str(protein).replace("*", ""))
print(f"MW {pa.molecular_weight()/1000:.1f} kDa, pI {pa.isoelectric_point():.2f}, "
      f"ε280 {pa.molar_extinction_coefficient()[1]} M⁻¹cm⁻¹")

# hydropathy profile: find transmembrane segments
from Bio.SeqUtils import ProtParamData
kd = pa.protein_scale(ProtParamData.kd, window=19)
fig, ax = plt.subplots(figsize=(9, 2.5))
ax.plot(kd, lw=0.9)
ax.axhline(1.6, ls="--", c="crimson")           # a common TM threshold
ax.set(xlabel="Residue", ylabel="Kyte-Doolittle hydropathy")

Common mistakes

• SeqUtils.GC(): removed, and 100× different from gc_fraction.
• Tm_Wallace for real primer design.
• Ignoring salt and Mg²⁺ in Tm calculation.
• ProteinAnalysis on a Seq: needs str.
• ProteinAnalysis on a sequence containing * or X → error.
• Counting N in the GC denominator without deciding to. Use ambiguous=.
• Average vs monoisotopic mass confusion in proteomics.
• Not case-normalising soft-masked sequence before composition analysis.

See also

Seq · SeqIO · Bio.motifs · Biopython Deprecations · Applied - Sequence Composition Analysis

Exercise

gc_fraction, GC123 and Tm_Wallace are three small computations, and Biopython is not installed here, so write them. Doing so also makes the removed GC() trap concrete: the fraction and the percentage differ by exactly a factor of 100, and only one of them is what a threshold in old code expects.

Write profile(text), taking a FASTA file as one string and returning one report line per record. Take the record id as the header text up to the first whitespace, and join wrapped sequence lines before analysing. Upper-case the sequence, then count only A, C, G and T in both the numerator and the denominator of the GC calculations, so an ambiguous N is excluded rather than counted as a non-GC base. GC3 uses the third position of every complete codon, so a trailing partial codon is ignored. The Wallace melting temperature is 2 * (A + T) + 4 * (G + C) over the whole upper-cased sequence. Print each record as <id> len=<n> gc=<x> gc3=<x> tm=<n> with both fractions to three decimals and the length counted before any filtering. A record with no unambiguous base has gc=0.000, and one with no complete codon has gc3=0.000. If the input holds no records at all, print no records.