# Python Bio Sets

> Python set ops (union/intersection/difference) and collections.Counter for gene-list comparisons and k-mer/codon counting. Use when comparing gene lists, finding shared orthologs, computing k-mer Jaccard similarity, or tallying GC/codon usage.

- Skill: `pavel-kravchenko/python-bio-sets` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/python-bio-sets`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/python-bio-sets/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

| Operation | Syntax | Meaning |
|---|---|---|
| Union | `A \| B` | All elements in either set |
| Intersection | `A & B` | Elements in both sets |
| Difference | `A - B` | Elements in A but not B |
| Symmetric diff | `A ^ B` | Elements in exactly one set |
| Subset test | `A <= B` | All of A is in B |

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

```python
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.

```python
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.

```python
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.