bio-sequence-properties

What this skill does

## Version Compatibility

Reference examples tested with: BioPython 1.83+, samtools 1.19+

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 biological sequences using Biopython.

**"Calculate GC content"** -> Compute the fraction of G+C bases in a nucleotide sequence.
- Python: `gc_fraction(seq)` (BioPython SeqUtils)

**"Analyze protein properties"** -> Compute MW, pI, stability, hydrophobicity from an amino acid sequence.
- Python: `ProteinAnalysis(str_seq)` (BioPython ProtParam)

## Required Imports

```python
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
```

## DNA/RNA Properties

### GC Content

```python
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:
```python
gc = gc_fraction(seq, ambiguous='ignore')  # Ignore N bases in calculation
gc = gc_fraction(seq, ambiguous='weighted')  # Weight by probability
```

### GC at Codon Positions (GC123)

Analyze GC content at each codon position (useful for codon bias analysis):

```python
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)
```

### GC Skew

Calculate (G-C)/(G+C) in sliding windows to identify replication origins:

```python
from Bio.SeqUtils import GC_skew

seq = Seq('ATGCGATCGATCGATCGATCG' * 10)
skew_values = GC_skew(seq, window=100)  # Returns list of skew values
```

### Molecular Weight

```python
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')
```

### Melting Temperature

```python
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)
```

### Sequence Search with IUPAC Codes (nt_search)

Built-in search that handles IUPAC ambiguity codes:

```python
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)
```

### Base Composition

```python
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'}
```

## Amino Acid Code Conversion

Convert between 1-letter and 3-letter amino acid codes:

```python
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'
```

## Protein Properties

**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.

### Using ProteinAnalysis

```python
from Bio.SeqUtils.ProtParam import ProteinAnalysis

protein = ProteinAnalysis('MAEGEITTFTALTEKFNLPPGNYKKPKLLYCSNG')
```

### Basic Properties

```python
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
```

### Stability and Hydropathy

```python
instability = protein.instability_index()  # < 40 = stable, > 40 = unstable
gravy = protein.gravy()                    # Negative = hydrophilic, Positive = hydrophobic
arom = protein.aromaticity()               # Fraction of Phe+Trp+Tyr
```

### Charge at Specific pH

```python
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 Profile

```python
flexibility = protein.flexibility()  # List of flexibility values per residue
# Based on Vihinen, 1994 scale
```

### Secondary Structure

```python
helix, turn, sheet = protein.secondary_structure_fraction()
```

### Extinction Coefficient

```python
eps_reduced, eps_cystine = protein.molar_extinction_coefficient()
# eps_reduced: all Cys are reduced
# eps_cystine: all Cys form disulfide bonds
```

### Protein Scale Profiles

Calculate profiles using any amino acid scale:

```python
# Kyte-Doolittle hydropathy profile
kd_scale = {'A': 1.8, 'R': -4.5, 'N': -3.5, ...}
profile = protein.protein_scale(kd_scale, window=7)
```

## Code Patterns

### Analyze Multiple Sequences

```python
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
```

### GC Content Distribution (Sliding Windows)

```python
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
```

### GC Skew Plot Data

```python
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
```

### Full Protein Analysis Report

**Goal:** Generate a comprehensive summary of protein biophysical properties.

**Approach:** Compute all key metrics from a single `ProteinAnalysis` object and return as a dictionary.

```python
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,
    }
```

### Dinucleotide Frequencies

```python
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: