# Python Bio Sequences

> Manipulate DNA/RNA/protein sequences as raw Python strings: reverse complement via str.maketrans/translate, transcription, codon/ORF extraction, motif and restriction-site scanning with find()/re, and hand-rolled FASTA parsing without Biopython. Use when writing sequence utilities from scratch, debugging off-by-one slicing or 1-based-vs-0-based coordinate errors, or when Biopython/Seq is unavailable or overkill.

- Skill: `pavel-kravchenko/python-bio-sequences` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/python-bio-sequences`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/python-bio-sequences/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-bio-sequences

---


# Biological Sequences as Python Strings

## When to Use

- Writing a small sequence utility (reverse complement, GC%, codon split) without pulling in Biopython.
- Debugging a reverse-complement, `find()`, or slicing bug that gives subtly wrong results.
- Parsing FASTA/FASTQ-like text by hand (headers, multi-line sequences) in a script or notebook.
- Converting between 0-based (Python/BED) and 1-based (VCF/GFF/SAM) coordinates.
- Scanning for motifs, restriction sites, or IUPAC degenerate patterns in raw sequence text.

## Version Compatibility

Pure standard library — Python ≥3.8 (f-strings, `str.maketrans`/`translate`). No third-party dependencies. For anything beyond ad-hoc string ops (real FASTA/FASTQ I/O, alphabets, translation tables), switch to Biopython (see `biopython` skill).

## Prerequisites

- Basic Python: strings, slicing, list comprehensions, `re` module.
- Concepts: DNA/RNA/protein alphabets, codons, reading frames, 5'→3' orientation.

## Core Operations

**Goal:** compute reverse complement, transcribe, and split into codons correctly.
**Approach:** use `str.maketrans` + `translate` for complement (never chained `replace()`), `[::-1]` to reverse, and step-3 slicing for codons.

```python
def reverse_complement(dna: str) -> str:
    """Return the reverse complement of a DNA sequence (5'->3' input and output)."""
    complement_table = str.maketrans('ATGC', 'TACG')
    return dna.upper().translate(complement_table)[::-1]

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

def gc_content(seq: str) -> float:
    """Fraction of G/C bases in seq (0.0-1.0)."""
    seq = seq.upper()
    return (seq.count('G') + seq.count('C')) / len(seq)

def extract_codons(dna: str, reading_frame: int = 1) -> list[str]:
    """Split dna into complete codons for reading_frame 1, 2, or 3 (1-based frame)."""
    dna = dna.upper()
    start = reading_frame - 1
    return [dna[i:i + 3] for i in range(start, len(dna) - 2, 3)]
```

## Motif and Restriction-Site Scanning

**Goal:** find all (possibly overlapping) positions of a motif, including IUPAC degenerate patterns.
**Approach:** loop `str.find()` with `start = pos + 1` for overlapping matches (`count()` is non-overlapping and undercounts); use `re` for degenerate IUPAC codes.

```python
import re

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

def find_restriction_sites(sequence: str, site: str) -> list[int]:
    """1-based positions of a restriction site on the forward strand.
    Callers should also scan reverse_complement(sequence) for palindromic
    or antisense-strand sites (e.g. EcoRI GAATTC is palindromic).
    """
    return [p + 1 for p in find_all(sequence, site)]

# IUPAC degenerate pattern example: TATA box (TATAWAW, W = A or T)
tata_hits = re.findall(r'TATA[AT]A[AT]', "GGGTATAAAAGGGTATATAT".upper())
```

## FASTA Parsing (No Biopython)

**Goal:** parse a multi-record FASTA string/file into `{id, description, sequence}` records.
**Approach:** split on `>`, take the first whitespace-separated token of the header line as the id, join remaining lines as the sequence.

```python
def parse_fasta(fasta_text: str) -> list[dict]:
    """Parse a multi-record FASTA string into a list of dicts with
    'id', 'description', and 'sequence' keys."""
    records = []
    for entry in fasta_text.strip().split('>'):
        if not entry.strip():
            continue
        lines = entry.strip().split('\n')
        header_parts = lines[0].split(None, 1)
        records.append({
            'id': header_parts[0],
            'description': header_parts[1] if len(header_parts) > 1 else '',
            'sequence': ''.join(line.strip() for line in lines[1:]),
        })
    return records

def read_fasta_file(filepath: str) -> dict[str, str]:
    """Stream a FASTA file into {header: sequence} without loading it as one string."""
    records = {}
    with open(filepath) as f:
        header, seq = None, []
        for line in f:
            line = line.strip()
            if line.startswith('>'):
                if header:
                    records[header] = ''.join(seq)
                header, seq = line[1:], []
            else:
                seq.append(line)
        if header:
            records[header] = ''.join(seq)
    return records
```

## Coordinate Systems

| Format | Base | Interval type | First 3 bp |
|--------|------|--------------|------------|
| Python | 0 | Half-open `[start, stop)` | `seq[0:3]` |
| BED | 0 | Half-open | `start=0, end=3` |
| VCF/GFF | 1 | Closed `[start, stop]` | `start=1, end=3` |
| SAM | 1 | Closed | `POS=1` |

Converting GFF→Python: `python_start = gff_start - 1` (stop stays the same for half-open slicing).

## Pitfalls

- **Chained `replace()` for complement is wrong**: `dna.replace('A','T').replace('T','A')` converts everything to `A` — use `str.maketrans`/`translate` which substitutes simultaneously.
- **`find()` returns `-1` on miss, not `None`**: `if pos:` is truthy for `-1`; always check `if pos != -1:`.
- **`count()` is non-overlapping**: `"ATATATAT".count("ATAT")` is `2`, not 3; use the `find()`-loop in `find_all()` above for overlaps.
- **Case sensitivity**: always `.upper()` before scanning — lowercase denotes soft-masked repeats in UCSC/Ensembl output.
- **Strings are immutable**: `seq[0] = 'C'` raises `TypeError`; build a new string via slicing/concatenation.
- **Off-by-one from 1-based coordinates**: subtract 1 when converting VCF/GFF/SAM positions to Python indices; forgetting this shifts every downstream slice by one base.
- **Reading-frame math**: frame *N* starts at index `N - 1`, not `N` — `extract_codons(dna, 2)` starts at index 1.

## See Also

- `biopython` — real `Seq`/`SeqRecord` objects, alphabets, translation tables, and robust FASTA/FASTQ I/O for anything beyond quick scripting.
- `bio-sequence-io-read-sequences` — file-based sequence I/O patterns (FASTA/FASTQ, compressed files).
- `bio-sequence-manipulation-reverse-complement` — dedicated reverse-complement recipes and edge cases (ambiguity codes, RNA).
- `bio-restriction-analysis-restriction-sites` — enzyme recognition-site databases and mapping beyond simple string search.