vdjdb-harmonize

name: vdjdb-harmonize description: Harmonize antigen.gene and antigen.species fields in a VDJdb chunk to canonical VDJdb naming. Detects spurious gene/species names, resolves inconsistencies (same epitope → multiple names), and warns about epitopes that are exact substrings of longer epitopes. Invoked standalone or from /proofread when spurious values are detected.

/harmonize — VDJdb Antigen Gene/Species Harmonization Skill

Purpose

Normalize antigen.gene and antigen.species values in a chunk to the VDJdb canonical vocabulary. This skill uses a unified, layered lookup system drawing from two patch sources:

• patches/antigen_epitope_species_gene.dict — epitope-keyed authority: maps known epitopes directly to their canonical species + gene. This is the highest-confidence source.
• proofreading/gene_aliases.tsv — free-text gene name aliases → VDJdb gene name.
• proofreading/species_aliases.tsv — species substring fragments → VDJdb species name. Checked in order; first match wins.

Invocation

/harmonize [path-to-tsv]

Or called automatically from /proofread (Step 6a) when spurious antigen fields are detected.

Step 1 — Load Lookup Tables

Load all three sources at the start of the session:

import csv, re

# 1a. Epitope → (species, gene) from patches/
epitope_dict = {}  # {epitope: (species, gene)}
with open('patches/antigen_epitope_species_gene.dict') as f:
    reader = csv.DictReader(f, delimiter='\t')
    for row in reader:
        ep = row['antigen.epitope'].strip()
        epitope_dict[ep] = (row['antigen.species'].strip(), row['antigen.gene'].strip())

# 1b. Gene alias table
gene_map = {}  # {raw_name: canonical}
with open('proofreading/gene_aliases.tsv') as f:
    for line in f:
        line = line.strip()
        if not line or line.startswith('#'): continue
        parts = line.split('\t')
        if len(parts) >= 2 and parts[0] not in ('source_name', 'vdjdb_name'):
            gene_map[parts[0].strip()] = parts[1].strip()

# 1c. Species fragment list (order-sensitive)
species_fragments = []  # [(fragment_lower, canonical), ...]
with open('proofreading/species_aliases.tsv') as f:
    for line in f:
        line = line.strip()
        if not line or line.startswith('#'): continue
        parts = line.split('\t')
        if len(parts) >= 2 and parts[0] != 'fragment':
            species_fragments.append((parts[0].strip().lower(), parts[1].strip()))

Step 2 — Harmonization Functions

Use this unified pipeline for every row:

GENE_STRIP_SUFFIXES = (' protein', ' glycoprotein', ' polyprotein', ' precursor')

def harmonize_gene(raw: str) -> str:
    """Map raw antigen.gene to VDJdb canonical name."""
    if not raw:
        return raw
    # Strip embedded [species] annotation (e.g. "pp65 [CMV]")
    s = re.sub(r'\s*\[[^\]]+\]\s*$', '', raw).strip()
    # Exact match
    if s in gene_map:
        return gene_map[s]
    # Strip common trailing suffixes and retry
    for suffix in GENE_STRIP_SUFFIXES:
        if s.lower().endswith(suffix) and len(s) > len(suffix):
            stripped = s[:-len(suffix)].strip()
            if stripped in gene_map:
                return gene_map[stripped]
            if len(stripped) >= 2:
                return stripped   # cleaned fallback
    return s

def harmonize_species(raw: str) -> str:
    """Map raw antigen.species to VDJdb canonical name via substring fragments."""
    if not raw:
        return raw
    low = raw.lower()
    for fragment, canonical in species_fragments:
        if fragment in low:
            return canonical
    return raw  # unchanged if no match

def harmonize_row(row: dict) -> dict:
    """Apply full harmonization pipeline to a single row."""
    ep = row.get('antigen.epitope', '').strip()

    # Priority 1: epitope dict (most authoritative)
    if ep in epitope_dict:
        canon_species, canon_gene = epitope_dict[ep]
        row['antigen.species'] = canon_species
        row['antigen.gene'] = canon_gene
        return row

    # Priority 2: alias tables
    raw_gene = row.get('antigen.gene', '').strip()
    raw_species = row.get('antigen.species', '').strip()
    new_gene = harmonize_gene(raw_gene)
    new_species = harmonize_species(raw_species)

    if new_gene != raw_gene:
        row['antigen.gene'] = new_gene
    if new_species != raw_species:
        row['antigen.species'] = new_species

    return row

Step 3 — Detect Spurious Values

Before harmonizing, scan the chunk and flag any rows where antigen.gene or antigen.species looks suspicious. These are the patterns that trigger automatic /harmonize invocation from /proofread:

3a. Spurious antigen.gene indicators

VDJdb gene naming conventions:

• Human/mouse genes: HGNC uppercase symbol (e.g., PABPC1, SMC1A, COL18A1)
• Viral genes: use established short names from literature (e.g., pp65, BMLF1, Gag, Pol, Tax)
• Bacterial/parasitic genes: use standard gene symbol or protein abbreviation (e.g., glnA, GRA6, yeiH)
• Avoid: full UniProt protein descriptions, parenthetical qualifiers, Chain [X], prefixes from PDB entries
• Special cases: P protein (HBV) → Pol; LCMV Gp33(variant) → GPC; MULTISPECIES: prefix → strip prefix then look up

KNOWN_MULTIWORD_GENES = {'PB1-F2', 'non-structural', 'Large T antigen', 'HLA-DRB1',
                          'HLA-DQB1', 'HLA-DPB1', 'NY-ESO-1', 'CORT_0A05310'}  # extend as needed

def is_spurious_gene(gene: str) -> bool:
    if not gene: return True
    if re.search(r'\[.+\]', gene): return True
    if any(gene.lower().endswith(s) for s in GENE_STRIP_SUFFIXES): return True
    if len(gene) > 30 and gene not in KNOWN_MULTIWORD_GENES: return True
    if gene.lower() in ('polyprotein', 'unknown', 'na', 'n/a'): return True
    return False

3b. Spurious antigen.species indicators

VDJdb species naming conventions:

• Two-word binomial names → CamelCase with no space: Homo sapiens → HomoSapiens, Bacillus subtilis → BacillusSubtilis
• Common abbreviations for well-known pathogens: EBV, CMV, HIV-1, SARS-CoV-2, InfluenzaA, etc.
• Strip parenthetical common names and former names: Columba livia (carrier pigeon) → ColumbaLivia; Schinkia azotoformans (Bacillus azotiformans) → SchinkiaAzotoformans
• Genus-only entries (when species is unknown): keep as single CamelCase word: Bacillus [genus] → Bacillus
• PseudomonasFluorescens, PseudomonasAeruginosa — already CamelCase but missing space between genus and species (both are acceptable as-is if already in database)

Known canonical antigen.species values (derive from existing chunks):

EBV, CMV, MCMV, HSV-1, HSV-2, VZV, InfluenzaA, InfluenzaB, SARS-CoV-2, SARS-CoV,
HCoV-OC43, HCoV-HKU1, HIV-1, HCV, HBV, YFV, DENV, DENV1, DENV2, DENV3, DENV3/4,
LCMV, HTLV-1, MLV, RSV, MCPyV, HomoSapiens, MusMusculus, RattusNorvegicus,
MacacaMulatta, GallusGallus, M.tuberculosis, PlasmodiumFalciparum, PlasmodiumBerghei,
Trypanosoma cruzi, SIV, CoxsackievirusB, Wheat, ManducaSexta

CANONICAL_SPECIES = {
    'EBV', 'CMV', 'MCMV', 'HSV-1', 'HSV-2', 'VZV', 'InfluenzaA', 'InfluenzaB',
    'SARS-CoV-2', 'SARS-CoV', 'HCoV-OC43', 'HCoV-HKU1', 'HIV-1', 'HCV', 'HBV',
    'YFV', 'DENV', 'DENV1', 'DENV2', 'DENV3', 'DENV3/4', 'LCMV', 'HTLV-1',
    'MLV', 'RSV', 'MCPyV', 'HomoSapiens', 'MusMusculus', 'RattusNorvegicus',
    'MacacaMulatta', 'GallusGallus', 'M.tuberculosis', 'PlasmodiumFalciparum',
    'PlasmodiumBerghei', 'Trypanosoma cruzi', 'SIV', 'CoxsackievirusB',
    'Wheat', 'ManducaSexta',
}

def is_spurious_species(species: str) -> bool:
    if not species: return True
    if species in CANONICAL_SPECIES: return False
    if ' ' in species: return True   # multi-word → likely not normalised
    return True  # unknown single token → flag

Step 4 — Cross-Chunk Consistency Check

After harmonizing individual rows, check for same-epitope inconsistencies across the chunk (and optionally across all of chunks/):

from collections import defaultdict

def check_consistency(rows: list[dict]) -> list[str]:
    warnings = []
    ep_genes = defaultdict(set)
    ep_species = defaultdict(set)
    for r in rows:
        ep = r.get('antigen.epitope', '').strip()
        g = r.get('antigen.gene', '').strip()
        s = r.get('antigen.species', '').strip()
        if ep:
            if g: ep_genes[ep].add(g)
            if s: ep_species[ep].add(s)
    for ep, genes in ep_genes.items():
        if len(genes) > 1:
            warnings.append(f'INCONSISTENT antigen.gene for epitope {ep!r}: {genes} — check patches/antigen_epitope_species_gene.dict')
    for ep, sps in ep_species.items():
        if len(sps) > 1:
            warnings.append(f'INCONSISTENT antigen.species for epitope {ep!r}: {sps} — check patches/antigen_epitope_species_gene.dict')
    return warnings

If the epitope dict has a definitive entry, the inconsistency will be resolved by Step 2. Report remaining inconsistencies (those the dict does not cover) for manual curation.

Step 5 — Epitope Substring Warning

Check whether any epitope in the chunk is an exact substring of a longer epitope in the same chunk, or in the global VDJdb database. This detects truncation artefacts (e.g., EPLPQGQLTAY being a substring of GPEPLPQGQLTAY).

def check_epitope_substrings(rows: list[dict], global_epitopes: set[str] = None) -> list[str]:
    """
    Warn when epitope A is an exact substring of longer epitope B.
    Checks within the chunk and against global_epitopes if provided.
    """
    warnings = []
    chunk_epitopes = {r.get('antigen.epitope', '').strip() for r in rows if r.get('antigen.epitope')}
    all_epitopes = chunk_epitopes | (global_epitopes or set())

    for ep in sorted(chunk_epitopes):
        if len(ep) < 4: continue
        for longer in all_epitopes:
            if longer != ep and ep in longer:
                src = 'chunk' if longer in chunk_epitopes else 'global VDJdb'
                warnings.append(
                    f'SUBSTRING WARNING: {ep!r} is contained in longer epitope {longer!r} ({src}) — '
                    f'possible truncation artefact; verify correct epitope boundaries'
                )
    return warnings

To load global VDJdb epitopes for the cross-chunk check:

import glob

def load_global_epitopes() -> set[str]:
    eps = set()
    for f in glob.glob('chunks/*.txt'):
        with open(f) as fh:
            for row in csv.DictReader(fh, delimiter='\t'):
                ep = row.get('antigen.epitope', '').strip()
                if ep: eps.add(ep)
    return eps

Step 6 — Apply and Report

Run the full harmonization on the chunk, report all changes and warnings:

def harmonize_chunk(path: str, check_global: bool = True) -> None:
    with open(path) as f:
        rows = list(csv.DictReader(f, delimiter='\t'))
    fieldnames = list(rows[0].keys()) if rows else []

    changes = []
    for i, row in enumerate(rows):
        orig_gene = row.get('antigen.gene', '')
        orig_species = row.get('antigen.species', '')
        row = harmonize_row(row)
        if row['antigen.gene'] != orig_gene:
            changes.append(f'  ROW {row["chunk.id"]}: antigen.gene {orig_gene!r} → {row["antigen.gene"]!r}')
        if row['antigen.species'] != orig_species:
            changes.append(f'  ROW {row["chunk.id"]}: antigen.species {orig_species!r} → {row["antigen.species"]!r}')
        rows[i] = row

    # Consistency check
    consistency_warnings = check_consistency(rows)

    # Substring check
    global_eps = load_global_epitopes() if check_global else None
    substring_warnings = check_epitope_substrings(rows, global_eps)

    # Write back
    with open(path, 'w', newline='') as f:
        writer = csv.DictWriter(f, fieldnames=fieldnames, delimiter='\t')
        writer.writeheader()
        writer.writerows(rows)

    # Report
    print(f'=== HARMONIZATION REPORT: {path} ===')
    print(f'Rows processed: {len(rows)}')
    print(f'Fields changed: {len(changes)}')
    for c in changes: print(c)
    if consistency_warnings:
        print(f'\nConsistency warnings ({len(consistency_warnings)}):')
        for w in consistency_warnings: print(f'  {w}')
    if substring_warnings:
        print(f'\nSubstring warnings ({len(substring_warnings)}):')
        for w in substring_warnings: print(f'  {w}')
    if not changes and not consistency_warnings and not substring_warnings:
        print('No issues found.')

Integration with /proofread

When /proofread is running Step 6 (MHC Consistency) or scanning antigen.gene/antigen.species fields, it should invoke harmonization automatically if any of these are detected:

• is_spurious_gene(row['antigen.gene']) returns True for ≥1 row
• is_spurious_species(row['antigen.species']) returns True for ≥1 row
• check_consistency(rows) returns any warnings

From /proofread Step 6, add:

Step 6a — Antigen Harmonization Trigger

After running ChunkQC, scan all antigen.gene and antigen.species values using the spurious-value detectors from /harmonize. If any are flagged:

1. Report the count and examples to the user.
2. Ask: "Run /harmonize to fix these automatically? [y/n]"
3. If yes: run the full harmonization pipeline, then re-run ChunkQC to verify no regressions.

Patch File Maintenance

When harmonization fails for a value (no alias match, no epitope dict entry), add the mapping to the appropriate patch file rather than leaving it unfixed:

• New epitope → species/gene mapping → append to patches/antigen_epitope_species_gene.dict
• New gene alias (free-text → canonical) → append to proofreading/gene_aliases.tsv
• New species fragment (substring → canonical) → append to proofreading/species_aliases.tsv (put more-specific fragments before less-specific ones)

Reference Files