# Python Core Bio

> Write pure-Python DNA/RNA sequence code (reverse complement, GC%, translation, ORF finding) and parse FASTA/FASTQ with generators. Use when writing sequence utilities without Biopython or parsing bio files from scratch.

- Skill: `pavel-kravchenko/python-core-bio` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/python-core-bio`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/python-core-bio/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/python-core-bio

---


# Python Core for Bioinformatics

## When to Use
- Writing DNA/RNA/protein sequence manipulation code from scratch (no Biopython dependency)
- Parsing FASTA, FASTQ, CSV, TSV, or JSON files line-by-line or as generators
- Implementing GC content, codon tables, ORF finding, reverse complement, or melting temperature
- Designing reusable, memory-efficient functions for sequence pipelines
- Handling files too large to load entirely into memory

## Version Compatibility
Standard library only — Python ≥3.9 (uses `str.maketrans`, f-strings, `pathlib`). No third-party packages required; contrast with `biopython` which wraps the same operations behind `Seq`/`SeqRecord` objects.

## Prerequisites
- Core Python: strings, slicing, functions, `dict`/`set`, `with` context managers
- Concepts: 0-based indexing, string immutability, generators (`yield`)
- No installs needed (`csv`, `pathlib` are stdlib)

## DNA/RNA Alphabet and Validation

**Goal:** classify a sequence and validate its alphabet before processing.
**Approach:** use `set` membership (`<=`) against reference alphabets; always `.upper()` first since databases mix case.

```python
VALID_DNA = set("ATGC")
VALID_RNA = set("AUGC")

def detect_seq_type(seq: str) -> str:
    """Classify a sequence as DNA, RNA, Protein, or Unknown by its alphabet."""
    chars = set(seq.upper())
    if chars <= VALID_DNA: return "DNA"
    if chars <= VALID_RNA: return "RNA"
    if chars <= set("ACDEFGHIKLMNPQRSTVWY"): return "Protein"
    return "Unknown"
```

## Reverse Complement, GC Content, Transcription

**Goal:** compute the reverse complement and composition stats of a DNA string.
**Approach:** `replace()` cannot swap A<->T and G<->C in one pass (chaining `.replace('A','T').replace('T','A')` turns everything into `A`). Use `str.maketrans` + `str.translate` instead, then reverse with `[::-1]`.

```python
RC_TABLE = str.maketrans("ATGC", "TACG")

def reverse_complement(seq: str) -> str:
    """Return the reverse complement of a DNA sequence (5'->3')."""
    return seq.upper().translate(RC_TABLE)[::-1]

def complement(seq: str) -> str:
    """Return the complement only, no reversal (3'->5' strand read forward)."""
    return seq.upper().translate(RC_TABLE)

def gc_content(seq: str) -> float:
    """Return GC content as a percentage (0-100)."""
    s = seq.upper()
    return (s.count("G") + s.count("C")) / len(s) * 100

def transcribe(dna: str) -> str:
    """Transcribe a DNA coding (sense) strand to mRNA."""
    return dna.upper().replace("T", "U")

def reverse_transcribe(rna: str) -> str:
    """Convert mRNA back to a DNA coding strand."""
    return rna.upper().replace("U", "T")
```

## Translation and ORF Finding

**Goal:** translate DNA to protein and locate open reading frames.
**Approach:** walk the sequence 3 nucleotides at a time using the standard genetic code; stop at the first in-frame stop codon. For ORFs, scan all 3 forward frames for ATG...stop spans, and run the same function on the reverse complement for the other strand.

```python
CODON_TABLE = {
    'TTT':'F','TTC':'F','TTA':'L','TTG':'L','CTT':'L','CTC':'L','CTA':'L','CTG':'L',
    'ATT':'I','ATC':'I','ATA':'I','ATG':'M','GTT':'V','GTC':'V','GTA':'V','GTG':'V',
    'TCT':'S','TCC':'S','TCA':'S','TCG':'S','CCT':'P','CCC':'P','CCA':'P','CCG':'P',
    'ACT':'T','ACC':'T','ACA':'T','ACG':'T','GCT':'A','GCC':'A','GCA':'A','GCG':'A',
    'TAT':'Y','TAC':'Y','TAA':'*','TAG':'*','CAT':'H','CAC':'H','CAA':'Q','CAG':'Q',
    'AAT':'N','AAC':'N','AAA':'K','AAG':'K','GAT':'D','GAC':'D','GAA':'E','GAG':'E',
    'TGT':'C','TGC':'C','TGA':'*','TGG':'W','CGT':'R','CGC':'R','CGA':'R','CGG':'R',
    'AGT':'S','AGC':'S','AGA':'R','AGG':'R','GGT':'G','GGC':'G','GGA':'G','GGG':'G',
}
STOP_CODONS = {'TAA', 'TAG', 'TGA'}

def translate(dna: str) -> str:
    """Translate a DNA coding sequence to protein, stopping at the first stop codon."""
    dna = dna.upper()
    protein = []
    for i in range(0, len(dna) - 2, 3):
        aa = CODON_TABLE.get(dna[i:i+3], '?')
        if aa == '*':
            break
        protein.append(aa)
    return ''.join(protein)

def find_orfs(seq: str, min_len: int = 30) -> list[dict]:
    """Find all ATG...stop ORFs >= min_len bp in all 3 forward reading frames.

    Run again on reverse_complement(seq) and tag strand='-' to cover both strands.
    """
    seq = seq.upper()
    orfs = []
    for frame in range(3):
        i = frame
        while i <= len(seq) - 3:
            if seq[i:i+3] == 'ATG':
                j = i + 3
                while j <= len(seq) - 3:
                    if seq[j:j+3] in STOP_CODONS:
                        length = j + 3 - i
                        if length >= min_len:
                            orfs.append({'start': i, 'end': j + 3,
                                         'length': length, 'frame': frame + 1,
                                         'seq': seq[i:j+3]})
                        break
                    j += 3
            i += 3
    return orfs
```

## Melting Temperature, Motif Search, Sliding-Window GC

**Goal:** primer Tm estimation and motif/palindrome detection for restriction sites.
**Approach:** Wallace rule for short primers (<14 nt), salt-adjusted formula otherwise; scan motifs with repeated `str.find()`.

```python
def tm(primer: str) -> float:
    """Estimate primer melting temperature (Wallace rule below 14 nt, salt-adjusted above)."""
    s = primer.upper()
    a, t, g, c = s.count('A'), s.count('T'), s.count('G'), s.count('C')
    if len(s) < 14:
        return 2 * (a + t) + 4 * (g + c)
    return 64.9 + 41 * (g + c - 16.4) / len(s)

def find_motif(seq: str, motif: str) -> list[int]:
    """Return all 0-based start positions of motif in seq (overlaps included)."""
    positions, pos = [], seq.find(motif)
    while pos != -1:
        positions.append(pos)
        pos = seq.find(motif, pos + 1)
    return positions

def is_palindrome(seq: str) -> bool:
    """Check if a DNA sequence is a reverse-complement palindrome (e.g. EcoRI site)."""
    s = seq.upper()
    return s == s.translate(RC_TABLE)[::-1]

def sliding_gc(seq: str, window: int = 100, step: int = 1) -> list[tuple]:
    """Compute GC% in a sliding window across seq; returns (start, gc_pct) pairs."""
    s = seq.upper()
    return [(i, (s[i:i+window].count('G') + s[i:i+window].count('C')) / window * 100)
            for i in range(0, len(s) - window + 1, step)]
```

## FASTA/FASTQ Parsing and File I/O

**Goal:** read/write FASTA and FASTQ without loading the whole file into memory.
**Approach:** generator functions that `yield` one record at a time; state machine on lines starting with `>` (FASTA) or 4-line cycles (FASTQ, Phred+33).

```python
def parse_fasta(filename: str):
    """Yield (header, sequence) tuples, one record at a time (constant memory)."""
    header, parts = None, []
    with open(filename) as f:
        for line in f:
            line = line.strip()
            if not line:
                continue
            if line.startswith('>'):
                if header is not None:
                    yield header, ''.join(parts)
                header, parts = line[1:], []
            else:
                parts.append(line)
    if header is not None:
        yield header, ''.join(parts)

def write_fasta(seqs: dict, filename: str, width: int = 60) -> None:
    """Write a dict of {header: sequence} to a FASTA file, wrapped at `width` chars."""
    with open(filename, 'w') as f:
        for header, seq in seqs.items():
            f.write(f">{header}\n")
            for i in range(0, len(seq), width):
                f.write(seq[i:i+width] + '\n')

def parse_fastq(filename: str):
    """Yield (id, sequence, quality_scores) tuples; quality_scores are Phred+33 ints."""
    with open(filename) as f:
        while True:
            header = f.readline().strip()
            if not header:
                break
            seq = f.readline().strip()
            f.readline()  # '+' separator line
            qual = f.readline().strip()
            scores = [ord(c) - 33 for c in qual]
            yield header[1:], seq, scores
```

CSV/TSV round-trips use the stdlib `csv` module; build paths with `pathlib.Path`:

```python
import csv
from pathlib import Path

with open('genes.csv') as f:
    for row in csv.DictReader(f):
        gc = float(row['gc_content'])

with open('results.csv', 'w', newline='') as f:
    writer = csv.DictWriter(f, fieldnames=['gene', 'fold_change', 'p_value'])
    writer.writeheader()
    writer.writerows(results)

output = Path('results') / 'analysis' / 'output.csv'
output.parent.mkdir(parents=True, exist_ok=True)
for fasta_file in Path('data').glob('*.fasta'):
    ...
```

## Pitfalls

| Pitfall | Fix |
|---------|-----|
| `seq[0] = 'G'` raises TypeError | Strings are immutable: `seq = 'G' + seq[1:]` |
| Chained `.replace('A','T').replace('T','A')` for complement | All A's become T then all T's (including original A's) become A — use `str.maketrans`/`translate` |
| Mixed case breaks counting | Always `.upper()` before processing |
| `0.1 + 0.2 == 0.3` is False | Use `abs(a - b) < 1e-9` for float comparison |
| Mutable default arg `def f(items=[])` | Use `items=None`, then `if items is None: items = []` |
| `int()` on file values | All values from file reads are `str`; cast explicitly |
| `seq.count('ATG')` counts non-overlapping only; `find` loop counts overlaps | Pick the one matching biological intent |
| `is` vs `==` for None | Always `if x is None`, never `if x == None` |
| `f.close()` not called on error | Always use `with open(...) as f:` |
| Loading entire genome into RAM | Use generator-based parsers (`yield`) for large files |
| GC formula without parentheses | `(g + c) / total * 100`, not `g + c / total * 100` |
| Reading frame confusion | Frame 0 starts at index 0; `frame = pos % 3` |

## See Also
- `biopython` — SeqIO, Entrez, BLAST wrappers, SeqRecord objects (use once Biopython is available)
- `python-bio-file-operations` — deeper file I/O patterns (chunked reads, `pathlib`, error handling)
- `python-bio-sequences` — string/slicing fundamentals this skill builds on
- `python-bio-numpy` — vectorised sequence stats, expression data, DataFrames