pythonforbio.
[WASM idle]
Biopython07/15

Bio.SeqUtils

Starting sandbox…
Beginnerlesson

Bio.SeqUtils

In one line: small utility functions for sequence composition, molecular weight, melting temperature and codon usage — with one notorious removed function.

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

scratch

No output yet — run the code to populate this drawer.