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.

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

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.

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().

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

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:

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