Python Bio Functions

Functions for Bioinformatics

When to Use

• Writing reusable def helpers (GC content, translation, ORF finding) instead of one-off scripts
• Debugging a function that returns None unexpectedly, or a list/dict that "remembers" values across calls
• Designing a signature that mixes required args, defaults, *args, and **kwargs (e.g. a FASTA header builder)
• Speeding up a recursive or repeatedly-called function (Fibonacci-style DP, k-mer counting) with memoization
• Building a small toolkit of composable functions (transcribe/translate/reverse-complement) for a pipeline

Version Compatibility

Python ≥3.9 (built-in generic type hints like list[dict]). No third-party dependencies — everything here is stdlib (functools).

Prerequisites

• Comfortable with Python control flow (for, while, if)
• Basic string/dict/list operations
• No packages to install

Goal: Avoid the classic mutable-default-argument bug.

Approach: Python evaluates default argument expressions once, at function-definition time — not on every call. A mutable default ([], {}) is therefore shared and accumulates state across calls. Fix: default to None, create the mutable object inside the body.

def bad_collect(item, items=[]):
    """BUG: the same list object is reused across every call."""
    items.append(item)
    return items

def good_collect(item, items=None):
    """Fixed: a fresh list is created each call unless one is passed in."""
    if items is None:
        items = []
    items.append(item)
    return items

assert bad_collect('A') == ['A']
assert bad_collect('T') == ['A', 'T']       # surprise: state leaked in
assert good_collect('A') == ['A']
assert good_collect('T') == ['T']            # correct: independent call

Goal: Find open reading frames (ATG → stop codon) in all six reading frames.

Approach: Scan each of the 3 forward frames for ATG, then walk forward in-frame until a stop codon; repeat on the reverse complement for the other 3 frames.

def reverse_complement(sequence: str) -> str:
    """Return the reverse complement of a DNA sequence (case-insensitive)."""
    complement = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G'}
    return ''.join(complement.get(n, 'N') for n in sequence.upper()[::-1])


def find_orfs(sequence: str, min_length: int = 30) -> list[dict]:
    """Find ORFs (ATG -> stop codon) in the 3 forward reading frames.

    Args:
        sequence: DNA sequence to search.
        min_length: minimum ORF length in nucleotides.

    Returns:
        List of dicts with 'start', 'end', 'length', 'sequence' keys.
    """
    sequence = sequence.upper()
    stop_codons = {'TAA', 'TAG', 'TGA'}
    orfs = []
    for frame in range(3):
        i = frame
        while i < len(sequence) - 2:
            if sequence[i:i + 3] == 'ATG':
                for j in range(i + 3, len(sequence) - 2, 3):
                    if sequence[j:j + 3] in stop_codons:
                        orf_seq = sequence[i:j + 3]
                        if len(orf_seq) >= min_length:
                            orfs.append({'start': i, 'end': j + 3,
                                         'length': len(orf_seq), 'sequence': orf_seq})
                        i = j + 3
                        break
                else:
                    i += 3
                    continue
                continue
            i += 3
    return sorted(orfs, key=lambda x: x['start'])


def find_all_orfs(sequence: str, min_length: int = 30) -> list[dict]:
    """Find ORFs on both strands (all 6 reading frames)."""
    all_orfs = []
    for orf in find_orfs(sequence, min_length):
        orf['strand'] = '+'
        all_orfs.append(orf)
    rev_seq = reverse_complement(sequence)
    for orf in find_orfs(rev_seq, min_length):
        orf['strand'] = '-'
        all_orfs.append(orf)
    return sorted(all_orfs, key=lambda x: x['length'], reverse=True)

Goal: Build a flexible FASTA header from an ID plus arbitrary metadata, and compute simple sequence-comparison metrics.

Approach: **kwargs collects any number of key=value pairs into a dict; use it for optional, open-ended metadata instead of a long fixed parameter list.

def create_fasta_header(sequence_id: str, **metadata) -> str:
    """Build a FASTA header from an ID and arbitrary key=value metadata.

    create_fasta_header("BRCA1", organism="Homo sapiens", length=5500)
    -> ">BRCA1 [organism=Homo sapiens] [length=5500]"
    """
    header = f">{sequence_id}"
    for key, value in metadata.items():
        header += f" [{key}={value}]"
    return header


def hamming_distance(seq1: str, seq2: str) -> int:
    """Count differing positions between two equal-length sequences.

    Raises:
        ValueError: if the sequences differ in length.
    """
    if len(seq1) != len(seq2):
        raise ValueError(f"Sequences must be the same length ({len(seq1)} != {len(seq2)})")
    return sum(a != b for a, b in zip(seq1.upper(), seq2.upper()))


assert create_fasta_header("BRCA1", organism="Homo sapiens") == ">BRCA1 [organism=Homo sapiens]"
assert hamming_distance("ATGC", "ATGG") == 1

Goal: Cache expensive/repeated calls (e.g. recursive k-mer or DP-style computations) without hand-rolling a memo dict.

Approach: functools.lru_cache memoizes by argument values automatically; arguments must be hashable.

from functools import lru_cache

@lru_cache(maxsize=None)
def fibonacci_cached(n: int) -> int:
    """Fibonacci with automatic memoization (stand-in for any expensive, pure, recursive bio calc)."""
    if n < 2:
        return n
    return fibonacci_cached(n - 1) + fibonacci_cached(n - 2)

assert fibonacci_cached(30) == 832040

Pitfalls

• Mutable default arguments: def f(x, data=[]) shares one list across all calls — use data=None and create it inside.
• return vs print(): result = my_function() is None if the function only print()s instead of returning.
• *args/**kwargs order: required, then default, then *args, then keyword-only, then **kwargs — any other order is a SyntaxError.
• Lambdas can't hold statements: only a single expression — no if/for/assignment inside; use def for anything more complex.
• Recursion depth: Python's default recursion limit (~1000) will blow up on long-sequence recursive scans; prefer iterative loops or @lru_cache for deep/repeated recursion.

See Also

• bio-sequence-manipulation-reverse-complement
• bio-sequence-manipulation-codon-usage
• bio-sequence-io-read-sequences