Calculate sequence properties like GC content, molecular weight, isoelectric point, and GC skew using Biopython. Use when analyzing sequence composition, computing physical properties, or comparing sequences.
Reference examples tested with: BioPython 1.83+, samtools 1.19+
Before using code patterns, verify installed versions match. If versions differ:
pip show <package> then help(module.function) to check signaturesIf code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
Calculate physical and chemical properties of biological sequences using Biopython.
"Calculate GC content" → Compute the fraction of G+C bases in a nucleotide sequence.
gc_fraction(seq) (BioPython SeqUtils)"Analyze protein properties" → Compute MW, pI, stability, hydrophobicity from an amino acid sequence.
ProteinAnalysis(str_seq) (BioPython ProtParam)from Bio.Seq import Seq
from Bio.SeqUtils import gc_fraction, molecular_weight, GC123, GC_skew, nt_search, seq1, seq3
from Bio.SeqUtils.ProtParam import ProteinAnalysis
from Bio.SeqUtils import gc_fraction
seq = Seq('ATGCGATCGATCGATCGATCG')
gc = gc_fraction(seq) # Returns 0.476... (fraction)
gc_percent = gc * 100 # Convert to percentage
Handle ambiguous bases:
gc = gc_fraction(seq, ambiguous='ignore') # Ignore N bases in calculation
gc = gc_fraction(seq, ambiguous='weighted') # Weight by probability
Analyze GC content at each codon position (useful for codon bias analysis):
from Bio.SeqUtils import GC123
seq = Seq('ATGCGATCGATCGATCGATCG')
gc_total, gc_pos1, gc_pos2, gc_pos3 = GC123(seq)
# gc_total: overall GC%
# gc_pos1: GC% at 1st codon position
# gc_pos2: GC% at 2nd codon position
# gc_pos3: GC% at 3rd codon position (wobble)
Calculate (G-C)/(G+C) in sliding windows to identify replication origins:
from Bio.SeqUtils import GC_skew
seq = Seq('ATGCGATCGATCGATCGATCG' * 10)
skew_values = GC_skew(seq, window=100) # Returns list of skew values
from Bio.SeqUtils import molecular_weight
dna = Seq('ATGCGATCG')
mw = molecular_weight(dna) # Single-stranded DNA weight
# Double-stranded DNA
mw_ds = molecular_weight(dna, double_stranded=True)
# Circular DNA
mw_circ = molecular_weight(dna, circular=True)
# RNA
rna = Seq('AUGCGAUCG')
mw_rna = molecular_weight(rna, seq_type='RNA')
# Protein
protein = Seq('MRCRS')
mw_prot = molecular_weight(protein, seq_type='protein')
from Bio.SeqUtils import MeltingTemp
seq = Seq('ATGCGATCGATCG')
# Basic Tm (Wallace rule - short oligos)
tm_wallace = MeltingTemp.Tm_Wallace(seq)
# GC method (more accurate for longer sequences)
tm_gc = MeltingTemp.Tm_GC(seq)
# Nearest neighbor (most accurate)
tm_nn = MeltingTemp.Tm_NN(seq)
# With salt correction
tm_corrected = MeltingTemp.Tm_NN(seq, Na=50, Mg=1.5)
Built-in search that handles IUPAC ambiguity codes:
from Bio.SeqUtils import nt_search
seq = 'ATGCGATCGATCGATNGATC'
# Search for pattern with ambiguous bases
result = nt_search(seq, 'GATNGATC') # N matches any base
# Returns: ['GATNGATC', 8] - pattern and position(s)
def base_composition(seq):
seq_str = str(seq).upper()
total = len(seq_str)
return {base: seq_str.count(base) / total * 100 for base in 'ATGC'}
Convert between 1-letter and 3-letter amino acid codes:
from Bio.SeqUtils import seq1, seq3
# 3-letter to 1-letter
one_letter = seq1('MetAlaGlyTrp') # Returns 'MAGW'
# 1-letter to 3-letter
three_letter = seq3('MAGW') # Returns 'MetAlaGlyTrp'
# With custom separator
three_letter = seq3('MAGW', join='-') # Returns 'Met-Ala-Gly-Trp'
Goal: Compute physical and chemical properties of a protein from its amino acid sequence.
Approach: Create a ProteinAnalysis object, then call property methods. Remove stop codons (*) and non-standard residues (X) before analysis.
from Bio.SeqUtils.ProtParam import ProteinAnalysis
protein = ProteinAnalysis('MAEGEITTFTALTEKFNLPPGNYKKPKLLYCSNG')
mw = protein.molecular_weight() # Molecular weight in Daltons
pi = protein.isoelectric_point() # Isoelectric point
aa_comp = protein.amino_acids_percent # Amino acid composition (property)
aa_count = protein.count_amino_acids() # Raw amino acid counts
instability = protein.instability_index() # < 40 = stable, > 40 = unstable
gravy = protein.gravy() # Negative = hydrophilic, Positive = hydrophobic
arom = protein.aromaticity() # Fraction of Phe+Trp+Tyr
charge = protein.charge_at_pH(7.0) # Net charge at pH 7.0
charge_acidic = protein.charge_at_pH(4.0)
charge_basic = protein.charge_at_pH(10.0)
flexibility = protein.flexibility() # List of flexibility values per residue
# Based on Vihinen, 1994 scale
helix, turn, sheet = protein.secondary_structure_fraction()
eps_reduced, eps_cystine = protein.molar_extinction_coefficient()
# eps_reduced: all Cys are reduced
# eps_cystine: all Cys form disulfide bonds
Calculate profiles using any amino acid scale:
# Kyte-Doolittle hydropathy profile
kd_scale = {'A': 1.8, 'R': -4.5, 'N': -3.5, ...}
profile = protein.protein_scale(kd_scale, window=7)
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction
def analyze_fasta(filename):
results = []
for record in SeqIO.parse(filename, 'fasta'):
results.append({
'id': record.id,
'length': len(record.seq),
'gc': gc_fraction(record.seq) * 100
})
return results
from Bio.SeqUtils import gc_fraction
def gc_distribution(seq, window_size=100, step=50):
gc_values = []
for i in range(0, len(seq) - window_size + 1, step):
window = seq[i:i + window_size]
gc_values.append((i, gc_fraction(window) * 100))
return gc_values
from Bio.SeqUtils import GC_skew
def gc_skew_analysis(seq, window=1000):
skew = GC_skew(seq, window=window)
positions = list(range(0, len(seq) - window + 1, window))
cumulative = []
total = 0
for s in skew:
total += s
cumulative.append(total)
return positions, skew, cumulative
Goal: Generate a comprehensive summary of protein biophysical properties.
Approach: Compute all key metrics from a single ProteinAnalysis object and return as a dictionary.
from Bio.SeqUtils.ProtParam import ProteinAnalysis
def protein_report(sequence):
protein = ProteinAnalysis(str(sequence).replace('*', ''))
helix, turn, sheet = protein.secondary_structure_fraction()
return {
'length': len(sequence),
'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,
}
from collections import Counter
def dinucleotide_freq(seq):
seq_str = str(seq)
dinucs = [seq_str[i:i+2] for i in range(len(seq_str) - 1)]
counts = Counter(dinucs)
total = sum(counts.values())
return {di: count / total for di, count in counts.items()}
def cpg_ratio(seq):
seq_str = str(seq).upper()
cpg = seq_str.count('CG')
c = seq_str.count('C')
g = seq_str.count('G')
expected = (c * g) / len(seq_str) if len(seq_str) > 0 else 0
return cpg / expected if expected > 0 else 0
| Property | Function | Notes |
|---|---|---|
| GC content | gc_fraction() | Returns fraction (0-1) |
| GC by position | GC123() | Total, pos1, pos2, pos3 |
| GC skew | GC_skew() | (G-C)/(G+C) in windows |
| Molecular weight | molecular_weight() | In Daltons |
| Melting temp | MeltingTemp.Tm_*() | Multiple methods |
| IUPAC search | nt_search() | Handles ambiguity codes |
| Property | Method | Typical Range |
|---|---|---|
| MW | molecular_weight() | Daltons |
| pI | isoelectric_point() | 0-14 |
| Charge | charge_at_pH(pH) | Varies |
| Instability | instability_index() | <40 stable |
| GRAVY | gravy() | -2 to +2 |
| Aromaticity | aromaticity() | 0-0.15 |
| Flexibility | flexibility() | Per-residue list |
| Function | Input | Output |
|---|---|---|
seq1() | 'MetAlaGly' | 'MAG' |
seq3() | 'MAG' | 'MetAlaGly' |
| Error | Cause | Solution |
|---|---|---|
KeyError in ProteinAnalysis | Non-standard amino acid | Remove X, *, etc. |
| Wrong MW | Wrong seq_type | Specify 'RNA' or 'protein' |
| Invalid Tm | Sequence too short | Use >10 bp for accurate Tm |
| Empty GC_skew | Sequence shorter than window | Reduce window size |
Need sequence properties?
├── DNA/RNA sequence?
│ ├── GC content? → gc_fraction()
│ ├── GC at codon positions? → GC123()
│ ├── GC skew (replication origin)? → GC_skew()
│ ├── Molecular weight? → molecular_weight()
│ ├── Melting temperature? → MeltingTemp.Tm_NN()
│ └── Search with IUPAC codes? → nt_search()
├── Protein sequence?
│ └── Create ProteinAnalysis object
│ ├── Size → molecular_weight()
│ ├── Charge → isoelectric_point(), charge_at_pH()
│ ├── Stability → instability_index()
│ ├── Hydrophobicity → gravy()
│ ├── Flexibility → flexibility()
│ └── Structure → secondary_structure_fraction()
└── Convert amino acid codes?
└── seq1() / seq3()