Python Collections Regex

Python Collections & Regex for Bioinformatics

When to Use

• Counting k-mer or codon frequencies, or computing k-mer Jaccard similarity between two sequences.
• Parsing FASTA/FASTQ files as a memory-efficient generator instead of loading the whole file.
• Deduplicating or set-comparing gene lists (union/intersection/difference of gene sets).
• Grouping records by key (e.g., genes by chromosome) without manual if key not in dict checks.
• Searching sequences for regex patterns: restriction sites, stop/start codons, homopolymer runs, tandem repeats, or PROSITE motifs; parsing FASTA/GenBank/UniProt headers.

Version Compatibility

Pure standard library — collections, re, itertools. Applies to any Python ≥3.8 (uses f-strings and dict-ordering guarantees from 3.7+). No third-party packages required.

Prerequisites

• Core Python: dicts, list/dict/set comprehensions, generators (yield).
• Familiarity with basic regex syntax (., *, +, [], (), {}) is assumed, not taught from scratch.
• For actual sequence I/O beyond raw strings (quality trimming, alignment), see bio-sequence-io-read-sequences and biopython.

Data Structure Cheat Sheet

Core Patterns

Goal: count and compare k-mer composition between two sequences. Approach: slide a window with a generator expression into Counter; use Counter set-like ops (& = min per key, + = sum) for similarity metrics.

from collections import Counter

def kmer_counts(seq: str, k: int) -> Counter:
    """Return a Counter of all overlapping k-mers of length k in seq."""
    return Counter(seq[i:i + k] for i in range(len(seq) - k + 1))

def kmer_jaccard(seq_a: str, seq_b: str, k: int = 3) -> float:
    """Jaccard similarity of the two sequences' k-mer sets (ignores counts)."""
    set_a, set_b = set(kmer_counts(seq_a, k)), set(kmer_counts(seq_b, k))
    return len(set_a & set_b) / len(set_a | set_b)

Goal: stream FASTA/FASTQ files record-by-record without holding the whole file in memory. Approach: a generator that yields as soon as a full record is assembled; chain generators to build lazy pipelines (nothing is computed until consumed).

def read_fasta(filename):
    """Yield (header, sequence) tuples. Memory: O(single record), not O(file)."""
    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 read_fastq(filename):
    """Yield {'id', 'sequence', 'quality'} dicts, 4 lines per record."""
    with open(filename) as f:
        while True:
            header = f.readline().strip()
            if not header:
                break
            seq = f.readline().strip()
            f.readline()  # '+' separator line, discarded
            qual = f.readline().strip()
            yield {'id': header[1:].split()[0], 'sequence': seq, 'quality': qual}

def quality_filter(records, min_q=30):
    """Keep only records whose mean Phred quality (Phred+33) is >= min_q."""
    for r in records:
        mean_q = sum(ord(c) - 33 for c in r['quality']) / len(r['quality'])
        if mean_q >= min_q:
            yield r

# Nothing runs until you iterate — lazy, constant-memory pipeline:
pipeline = quality_filter(read_fastq('reads.fastq'))

Goal: find restriction sites and predict digest fragment sizes. Approach: precompile enzyme recognition sequences, re.finditer for positions, re.split for fragments.

import re

RESTRICTION_ENZYMES = {
    'EcoRI': 'GAATTC', 'BamHI': 'GGATCC', 'HindIII': 'AAGCTT',
    'NotI': 'GCGGCCGC', 'XhoI': 'CTCGAG', 'NdeI': 'CATATG',
}

def restriction_map(seq: str, enzymes: dict) -> dict:
    """Map each enzyme name -> list of cut-site start positions found in seq."""
    return {name: [m.start() for m in re.finditer(site, seq)]
            for name, site in enzymes.items() if re.search(site, seq)}

def predict_fragments(seq: str, site: str) -> list:
    """Return fragment lengths after cutting seq at every occurrence of site."""
    return [len(frag) for frag in re.split(site, seq) if frag]

Goal: convert a PROSITE pattern (e.g., N-{P}-[ST]-{P}) into a Python regex. Approach: walk the pattern token by token; [...] stays a char class, {X} (excluded residues) becomes [^X], (n) becomes a {n} repeat count, x becomes ..

def prosite_to_regex(pattern: str) -> str:
    """Translate a PROSITE motif, e.g. 'N-{P}-[ST]-{P}', to 'N[^P][ST][^P]'."""
    result, parts = [], pattern.replace('-', '')
    i = 0
    while i < len(parts):
        if parts[i] == '[':
            j = parts.index(']', i); result.append(parts[i:j + 1]); i = j + 1
        elif parts[i] == '{':
            j = parts.index('}', i); result.append(f'[^{parts[i+1:j]}]'); i = j + 1
        elif parts[i] == '(':
            j = parts.index(')', i); result.append(f'{{{parts[i+1:j]}}}'); i = j + 1
        elif parts[i] == 'x':
            result.append('.'); i += 1
        else:
            result.append(parts[i]); i += 1
    return ''.join(result)

PROSITE = {
    'N-glycosylation':     'N-{P}-[ST]-{P}',
    'PKC_phosphorylation': '[ST]-x-[RK]',
    'CK2_phosphorylation': '[ST]-x(2)-[DE]',
}

Set operations on gene lists

cancer = {"BRCA1", "TP53", "EGFR", "MYC", "KRAS"}
repair = {"BRCA1", "BRCA2", "ATM", "MLH1", "TP53"}
cancer | repair   # union
cancer & repair   # intersection -> {"BRCA1", "TP53"}
cancer - repair   # cancer-only
cancer ^ repair   # symmetric difference

Regex quick reference for sequence patterns

import re

re.findall(r'TAA|TAG|TGA', dna)           # stop codons
re.finditer(r'(.)\1{4,}', dna)            # homopolymer runs >= 5 bases
re.finditer(r'([ATGC]{3})\1+', dna)       # tandem repeats (backreference)
re.findall(r'[^ATGCatgcNn]', seq)         # non-standard bases
re.findall(r'(?=(ATG))', dna)             # overlapping matches via lookahead
re.findall(r'ATG.*?TAG', dna)             # lazy ORF match (vs greedy ATG.*TAG)

# UniProt SwissProt header: >sp|P04637|P53_HUMAN ...
header_pat = re.compile(
    r'>(?:sp|tr)\|(?P<accession>[^|]+)\|(?P<entry>\S+)\s+'
    r'(?P<desc>.+?)\s+OS=(?P<organism>[^=]+?)(?:\s+\w+=|$)'
)
# GenBank feature table
genes    = re.findall(r'/gene="([^"]+)"', gb_text)
prot_ids = re.findall(r'/protein_id="([^"]+)"', gb_text)
coords   = re.findall(r'(\d+)\.\.(\d+)', gb_text)

Pitfalls

See Also

• bio-sequence-io-read-sequences — Biopython-based FASTA/FASTQ parsing for production pipelines.
• bio-sequence-manipulation-motif-search — higher-level motif/domain search (PROSITE, Pfam) beyond raw regex.
• bio-restriction-analysis-restriction-sites — full restriction-enzyme site catalogs and digest simulation.
• biopython — general Biopython skill for Seq, SeqIO, and sequence objects.