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.

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.

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.

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

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.