Python Bio Sets

Sets and Counter in Bioinformatics

When to Use

• Comparing gene lists from two or more experiments/studies (shared hits, study-only genes, all genes mentioned).
• Finding conserved orthologs across organisms (human vs. mouse vs. zebrafish gene symbols).
• Deduplicating k-mers or checking a sequence's alphabet (validating DNA vs. RNA vs. ambiguous codes).
• Computing k-mer or amino-acid composition and comparing profiles (Jaccard similarity, codon usage bias).
• Building an O(1) membership filter against a large reference gene/oncogene list instead of scanning a list repeatedly.

Version Compatibility

Pure standard library — no version dependency. Applies to Python ≥3.8 (relies only on built-in set/frozenset and collections.Counter/defaultdict).

Prerequisites

• No external packages required (stdlib only: collections.Counter, collections.defaultdict).
• Familiarity with python-bio-dictionaries (Counter and defaultdict are dict subclasses) helps but is not required.

Set Operations Reference

Gene Set Comparison

Goal: find shared/unique genes across two or more gene lists (e.g. cancer genes vs. DNA-repair genes, or orthologs across organisms). Approach: cast each list to a set, then combine with &, |, -, ^. For cross-organism comparisons, normalize case first since orthologous symbols differ in capitalization (human TP53 vs. mouse Trp53).

def compare_gene_sets(*named_gene_lists):
    """
    Compare an arbitrary number of gene sets.

    Parameters:
        named_gene_lists: dict mapping a label (e.g. organism/study name)
            to an iterable of gene symbols.
    Returns:
        dict with 'core' (in every set), 'union' (in any set), and
        'unique' (dict of label -> genes found only in that set).
    """
    sets = {name: set(genes) for name, genes in named_gene_lists[0].items()}
    core = set.intersection(*sets.values())
    union = set.union(*sets.values())
    unique = {}
    for name, genes in sets.items():
        others = set.union(*(g for n, g in sets.items() if n != name))
        unique[name] = genes - others
    return {"core": core, "union": union, "unique": unique}


cancer_genes = {"BRCA1", "TP53", "EGFR", "MYC", "KRAS"}
dna_repair_genes = {"BRCA1", "BRCA2", "ATM", "MLH1", "TP53"}

shared = cancer_genes & dna_repair_genes          # {'BRCA1', 'TP53'}
cancer_only = cancer_genes - dna_repair_genes      # genes not in DNA repair
exclusive = cancer_genes ^ dna_repair_genes        # in exactly one set

## Cross-organism comparison needs case normalization
human_genes = {"TP53", "BRCA1", "EGFR", "MYC", "KRAS", "RB1"}
mouse_genes = {"Trp53", "Brca1", "Egfr", "Myc", "Kras", "Rb1", "Pax6"}
orthologs = {g.upper() for g in human_genes} & {g.upper() for g in mouse_genes}

K-mer Sets, Jaccard Similarity, and Sequence Validation

Goal: dedupe k-mers, measure sequence similarity, or validate a sequence's alphabet. Approach: build a set of substrings with a sliding window; use set arithmetic for shared/unique k-mers and len(intersection) / len(union) for Jaccard similarity; use set difference against {'A','T','G','C'} to flag invalid characters.

def unique_kmers(sequence, k=3):
    """Return the set of distinct k-mers in a sequence."""
    return {sequence[i:i + k] for i in range(len(sequence) - k + 1)}


def jaccard_similarity(seq1, seq2, k=3):
    """Jaccard similarity of two sequences' k-mer sets (0.0-1.0)."""
    kmers1, kmers2 = unique_kmers(seq1, k), unique_kmers(seq2, k)
    if not kmers1 and not kmers2:
        return 0.0
    return len(kmers1 & kmers2) / len(kmers1 | kmers2)


def validate_dna(sequence):
    """Check whether a sequence contains only A/T/G/C (case-insensitive)."""
    valid = {'A', 'T', 'G', 'C'}
    found = set(sequence.upper())
    invalid = found - valid
    if invalid:
        return False, f"Invalid characters: {invalid}"
    return True, "Valid DNA"


seq1, seq2 = "ATGCGATCGATCG", "ATGCGGGGATCCC"
print(jaccard_similarity(seq1, seq2, k=3))
print(validate_dna("ATGCUGATC"))   # False -- 'U' is RNA, not DNA

Counter for Nucleotide/Codon Frequencies

Goal: tally nucleotide, k-mer, or amino-acid composition and compare frequency profiles between sequences. Approach: Counter is a dict subclass that counts hashable items; .most_common(n) ranks them; & (min counts) and + (sum counts) combine two Counters for profile comparison.

from collections import Counter, defaultdict

def gc_content(sequence):
    """Return GC percentage of a sequence using Counter."""
    counts = Counter(sequence.upper())
    return (counts['G'] + counts['C']) / len(sequence) * 100


def codon_usage_bias(dna_sequence, genetic_code):
    """
    Group codons by the amino acid they encode and count each codon.

    Parameters:
        dna_sequence: in-frame coding sequence (length multiple of 3 ideally)
        genetic_code: dict mapping codon -> amino acid letter
    Returns:
        dict: {amino_acid: Counter({codon: count, ...}), ...}
    """
    usable_len = len(dna_sequence) - len(dna_sequence) % 3
    codons = [dna_sequence[i:i + 3] for i in range(0, usable_len, 3)]
    usage = defaultdict(Counter)
    for codon in codons:
        aa = genetic_code.get(codon, 'X')
        usage[aa][codon] += 1
    return dict(usage)


sequence = "ATGCGATCGATCGTAGCGATCGATCGATGCGA"
nt_counts = Counter(sequence)
print(nt_counts.most_common(2))
print(f"GC%: {gc_content(sequence):.1f}")

kmers_a = Counter(seq1[i:i + 3] for i in range(len(seq1) - 2))
kmers_b = Counter(seq2[i:i + 3] for i in range(len(seq2) - 2))
shared_min = kmers_a & kmers_b     # per-kmer minimum count
combined = kmers_a + kmers_b       # per-kmer summed count

Pitfalls

• Sets are unordered: my_set[0] raises TypeError. Convert to sorted(my_set) when you need ordering.
• Set elements must be hashable: lists and dicts cannot be set members. Use frozenset for a set of sets, or as a dict key (e.g. mapping a group of synonymous codons to one amino acid).
• Empty set: {} creates an empty dict, not an empty set. Use set().
• in on a set is O(1): converting a reference gene list (e.g. known oncogenes) to a set once, before repeated membership tests, is a major speedup over scanning a list (O(n) per lookup).
• Case mismatches silently break intersections: cross-organism or cross-database gene symbols often differ only in case (TP53 vs. Trp53) — normalize with .upper()/.lower() before set operations, or you'll get an empty intersection with no error.
• Counter subtraction discards negatives: Counter(a) - Counter(b) drops any key whose result would be ≤ 0, unlike plain dict subtraction — use + and & deliberately, and check Counter docs if you need signed differences.

See Also

• python-bio-dictionaries — Counter/defaultdict are dict subclasses; covers core dict patterns.
• python-bio-strings — sequence slicing and string methods used to build k-mers here.
• python-bio-comprehensions — set/dict comprehension syntax used throughout ({seq[i:i+k] for i in ...}).
• bio-population-genetics-scikit-allel-analysis — set-like allele/variant comparisons at population scale.