GPTomics/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

| Method | Model | Use when | Caveat | |--------|-------|----------|--------| | Tm_Wallace | 4(G+C) + 2(A+T) "2+4 rule" | Oligos <=14 nt only | Ignores order/salt; WRONG for primers (off 5-10 C+) | | Tm_GC | GC-content empirical equation | Longer sequences, rough estimate | Composition-only, no nearest-neighbor info | | Tm_NN | Nearest-neighbor thermodynamics | PCR primers, probes, accurate work | Needs realistic salt/strand conc for absolute values |

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.

Amino-Acid Code Conversion

from Bio.SeqUtils import seq1, seq3

seq1('MetAlaGlyTrp')           # 'MAGW'  (3-letter -> 1-letter)
seq3('MAGW')                   # 'MetAlaGlyTrp'  (1-letter -> 3-letter, no separator)

Code Patterns

Per-Record GC Across a FASTA

Goal: Report length and GC percent for every record in a file.

Approach: Stream records with SeqIO.parse, compute GC per record, multiply the fraction by 100.

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

def analyze_fasta(filename):
    return [{'id': r.id, 'length': len(r.seq), 'gc': gc_fraction(r.seq) * 100} for r in SeqIO.parse(filename, 'fasta')]

Cumulative GC Skew (Headless)

Goal: Locate a candidate replication origin without the Tkinter graphics routine.

Approach: Take per-window skew from GC_skew, accumulate it, and read off the minimum (oriC) and maximum (terminus).

from Bio.SeqUtils import GC_skew

def cumulative_skew(seq, window=10000):
    skew = GC_skew(seq, window=window)
    positions, cumulative, total = [], [], 0
    for i, s in enumerate(skew):
        total += s
        positions.append(i * window)
        cumulative.append(total)
    ori = positions[cumulative.index(min(cumulative))]
    return positions, cumulative, ori

Full Protein Report

Goal: Summarize the key biophysical metrics of a protein in one pass.

Approach: Sanitize non-standard residues, build one ProteinAnalysis object, and collect each metric into a dict.

from Bio.SeqUtils.ProtParam import ProteinAnalysis

def protein_report(sequence):
    clean = str(sequence).upper().replace('*', '').replace('X', '')
    protein = ProteinAnalysis(clean)
    helix, turn, sheet = protein.secondary_structure_fraction()
    return {
        'length': len(clean),
        'molecular_weight': protein.molecular_weight(),
        'isoelectric_point': protein.isoelectric_point(),
        'charge_at_pH7': protein.charge_at_pH(7.0),
        'instability_index': protein.instability_index(),
        'gravy': protein.gravy(),
        'aromaticity': protein.aromaticity(),
        'helix_fraction': helix, 'turn_fraction': turn, 'sheet_fraction': sheet,
    }

CpG Observed/Expected Ratio

def cpg_ratio(seq):
    s = str(seq).upper()
    expected = (s.count('C') * s.count('G')) / len(s) if s else 0
    return s.count('CG') / expected if expected > 0 else 0

Property Reference

| Property | Function | Units / Notes | |----------|----------|---------------| | GC content | gc_fraction() | FRACTION 0-1 (multiply by 100 for percent) | | GC by codon position | GC123() | FOUR PERCENTAGES 0-100 (total + pos 1/2/3) | | GC skew | GC_skew() | (G-C)/(G+C) per window; 0 = no G/C in window | | Molecular weight | molecular_weight() | Daltons; single-strand by DEFAULT | | Melting temp | MeltingTemp.Tm_NN() | Celsius; accurate for primers | | pI / charge | isoelectric_point() / charge_at_pH() | Linear pKa only, ignores 3D/PTMs | | Instability | instability_index() | > 40 => predicted unstable | | Hydropathy | gravy() | neg = hydrophilic; default scale misspelled 'KyteDoolitle' |

Common Errors

| Symptom | Cause | Fix | |---------|-------|-----| | GC values look 100x too small; filters reject all | Used gc_fraction() (fraction) where percent expected | Multiply by 100 | | GC differs from legacy GC() on real data | Default ambiguous='remove' vs legacy 'ignore' | Use gc_fraction(seq, ambiguous='ignore') * 100 for a faithful drop-in | | GC_skew returns 0 across a region | Window had no G or C (guarded division), not true zero skew | Treat 0 as "no data"; widen the window | | xGC_skew raises TclError/ImportError | It is a Tkinter graphics routine, fails headless | Sum GC_skew() yourself for cumulative skew | | Primer Tm off by 5-10 C | Used Tm_Wallace/Tm_GC (composition-only) | Use Tm_NN with realistic salt and strand concentration | | Mg/dNTPs change nothing in Tm_NN | Mg/dNTPs ignored unless saltcorr is 6 or 7 | Set saltcorr=7 (Owczarzy 2008) for divalent correction | | MW ~half of expected for genomic dsDNA | molecular_weight default double_stranded=False | Pass double_stranded=True | | KeyError from ProteinAnalysis | Non-standard residue (B, Z, X, U, *, -) | Strip or replace before analysis | | KeyError from gravy('KyteDoolittle') | Default scale literal is misspelled 'KyteDoolitle' (one 't') | Omit the argument or pass 'KyteDoolitle' |

References

Lobry JR (1996) Asymmetric substitution patterns in the two DNA strands of bacteria. Mol Biol Evol 13(5):660-665.

Lobry JR, Gautier C (1994) Hydrophobicity, expressivity and aromaticity are the major trends of amino-acid usage in 999 Escherichia coli chromosome-encoded genes. Nucleic Acids Res 22(15):3174-3180.

SantaLucia J Jr (1998) A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics. PNAS 95(4):1460-1465.

Kyte J, Doolittle RF (1982) A simple method for displaying the hydropathic character of a protein. J Mol Biol 157(1):105-132.

Guruprasad K, Reddy BVB, Pandit MW (1990) Correlation between stability of a protein and its dipeptide composition: a novel approach for predicting in vivo stability of a protein from its primary sequence. Protein Eng 4(2):155-161.

Related Skills

• seq-objects - Create and modify Seq objects before property calculation
• codon-usage - GC123 and codon-bias indices for coding-sequence analysis
• transcription-translation - Translate a CDS before protein property analysis
• sequence-io/sequence-statistics - File-level statistics (N50, totals, dataset GC)
• primer-design/primer-basics - Design primers where Tm_NN and GC content drive the choices
• restriction-analysis/restriction-sites - Locate enzyme recognition sites in the same sequence