vdjdb-duplicates

/vdjdb-duplicates — VDJdb Duplicate & Consistency Audit Skill

Purpose

Scan all chunks/ files to:

1. Find duplicate TCRs at two CDR3 resolution levels.
2. Classify duplicates as within-publication, cross-publication same-lab, or genuinely independent.
3. Check author overlap between publications sharing many TCR sequences.
4. Flag extremely frequent records that suggest read-count inflation rather than unique T cell clones.
5. Report epitopes presented by multiple distinct MHC molecules or with inconsistent allele resolution.

Invocation

/vdjdb-duplicates

No arguments. Run from the repo root.

---

Step 1 — Load All Chunks

import csv, glob, re
from collections import defaultdict, Counter

rows_all = []
for path in sorted(glob.glob('chunks/*.txt')):
    fname = path.split('/')[-1]
    with open(path) as f:
        for row in csv.DictReader(f, delimiter='\t'):
            row['_file'] = fname
            rows_all.append(row)

def v(row, col): return (row.get(col) or '').strip()

---

Step 2 — Beta-Only Duplicates

Key: (cdr3.beta, v.beta, antigen.epitope) — records sharing the same beta chain and epitope regardless of alpha chain or donor metadata.

key_beta = defaultdict(list)
for row in rows_all:
    cb = v(row,'cdr3.beta'); vb = v(row,'v.beta'); ep = v(row,'antigen.epitope')
    if cb and ep:
        key_beta[(cb, vb, ep)].append(row)

dups_beta = {k: vs for k, vs in key_beta.items() if len(vs) > 1}

Classify each duplicate group:

Count and report each category. List the top 30 groups by frequency.

---

Step 3 — Paired Duplicates

Key: (cdr3.alpha, v.alpha, j.alpha, cdr3.beta, v.beta, j.beta, antigen.epitope) — exact paired chain duplicates. Only applied to rows where both CDR3 chains are non-empty.

key_pair = defaultdict(list)
for row in rows_all:
    ca = v(row,'cdr3.alpha'); cb = v(row,'cdr3.beta'); ep = v(row,'antigen.epitope')
    if ca and cb and ep:
        k = (ca, v(row,'v.alpha'), v(row,'j.alpha'),
             cb, v(row,'v.beta'),  v(row,'j.beta'), ep)
        key_pair[k].append(row)

dups_pair = {k: vs for k, vs in key_pair.items() if len(vs) > 1}

Report same categories as Step 2.

---

Step 4 — Author Overlap Check

For each cross-file duplicate group involving ≥2 distinct PMIDs, retrieve author lists and compute overlap:

import urllib.request, json, time

def get_authors(pmid):
    url = f'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=pubmed&id={pmid}&retmode=json'
    with urllib.request.urlopen(url, timeout=10) as r:
        data = json.load(r)
    return [a['name'] for a in data['result'][pmid].get('authors', [])]

# Build PMID-pair → shared author count
pmid_pairs = Counter()
for k, vs in dups_beta.items():
    refs = sorted({v(r,'reference.id') for r in vs if v(r,'reference.id').startswith('PMID:')})
    for i in range(len(refs)):
        for j in range(i+1, len(refs)):
            pmid_pairs[(refs[i], refs[j])] += 1

# For top pairs (> threshold shared duplicates), check authors
OVERLAP_THRESHOLD = 20
for (r1, r2), cnt in pmid_pairs.most_common():
    if cnt < OVERLAP_THRESHOLD: break
    p1, p2 = r1.replace('PMID:',''), r2.replace('PMID:','')
    try:
        a1 = get_authors(p1); time.sleep(0.4)
        a2 = get_authors(p2); time.sleep(0.4)
        shared = set(a1) & set(a2)
        print(f"{r1} × {r2}: {cnt} shared groups, {len(shared)} common authors")
        if shared:
            print(f"  Authors: {', '.join(sorted(shared)[:6])}")
    except Exception as e:
        print(f"  ERROR: {e}")

Classification rule:

• ≥3 shared authors → "same lab, follow-up publication" (expected overlap, not spurious)
• 0–2 shared authors → "independent replication" (genuine public clonotype)

---

Step 5 — High-Frequency Record Detection

Within-file groups with frequency ≥50 from ≤3 distinct donors are suspicious. Distinguish two legitimate assay causes before flagging as data errors:

Pattern A — Single TCR tested against many epitopes (combinatorial assay): The same CDR3 appears with many different epitopes in one study. Each row is a genuine specificity claim. High n within one epitope group from few donors still indicates read-depth, not clonal abundance.

Pattern B — Pool of TCRs tested against one epitope/pattern: Many different CDR3s all assigned the same epitope from a single donor. Here n reflects different clones in the repertoire, not read inflation. This is biologically expected in large repertoire studies.

Distinguishing them: if few donors but many distinct CDR3s → Pattern B (normal). If few donors with the SAME CDR3 repeated n times → Pattern A / read inflation.

FREQ_THRESHOLD = 50
DONOR_THRESHOLD = 3

for k, vs in sorted(dups_beta.items(), key=lambda x: -len(x[1])):
    files = {r['_file'] for r in vs}
    if len(files) > 1: continue   # only single-file
    donors = {v(r,'meta.subject.id') for r in vs}
    clones = {v(r,'meta.clone.id') for r in vs if v(r,'meta.clone.id')}
    if len(vs) >= FREQ_THRESHOLD and len(donors) <= DONOR_THRESHOLD:
        cb, vb, ep = k
        refs = {v(r,'reference.id') for r in vs}
        pattern = 'READ-INFLATION' if len(clones) <= 1 else 'DEEP-REPERTOIRE'
        print(f"{pattern} n={len(vs):4d} CDR3b={cb:22} ep={ep:15} donors={len(donors)} ref={list(refs)[0]}")

Known confirmed cases from vdjdb-db audit (2026-05-31):

Interpretation: Pattern A records represent sequencing read depth not unique T cells — flag in release notes but do not remove. Pattern B records are biologically valid and expected.

---

Step 6 — Multi-MHC Epitope Report

key_mhc = defaultdict(set)
for row in rows_all:
    ep = v(row,'antigen.epitope'); mhca = v(row,'mhc.a')
    if ep and mhca: key_mhc[ep].add(mhca)

# Flag epitopes with distinct HLA genes (not just allele sub-typing)
for ep, alleles in sorted(key_mhc.items(), key=lambda x: -len(x[1])):
    genes = set()
    for a in alleles:
        m = re.match(r'(HLA-[A-Z0-9]+|H2-\w+)', a)
        if m: genes.add(m.group(1))
    if len(genes) > 1:
        print(f"{ep:20} {len(alleles)} alleles, {len(genes)} genes: {sorted(genes)}")

# Flag allele resolution inconsistencies (same gene, coarse + fine)
for ep, alleles in key_mhc.items():
    coarse = {a for a in alleles if '*' in a and ':' not in a}
    fine   = {a for a in alleles if ':' in a}
    if coarse and fine:
        coarse_g = {a.split('*')[0] for a in coarse}
        fine_g   = {a.split('*')[0] for a in fine}
        if coarse_g & fine_g:
            print(f"RESOLUTION MIX {ep}: coarse={sorted(coarse)[:2]} fine={sorted(fine)[:2]}")

---

Step 7 — Summary Report

=== VDJDB DUPLICATE AUDIT SUMMARY ===
Total rows: N
Total chunks: N

BETA-ONLY DUPLICATES
  Unique duplicate groups: N
  Total redundant rows:    N
  Within-file:             N groups
  Cross-file:              N groups
    Same-lab (≥3 shared authors): N
    Independent replication:       N

TOP CROSS-PUBLICATION PAIRS (by shared CDR3b+Vb+epitope groups):
  [N]  PMIDX × PMIDY — [shared authors count] common authors — [lab relationship]
  ...

HIGH-FREQUENCY SUSPECTS (≥50 copies, ≤3 donors, single file):
  [N]  CDR3b / epitope — file — likely cause

MULTI-MHC EPITOPES (distinct HLA genes):
  [N epitopes] — list top cases

MHC FORMAT ISSUES (allele resolution inconsistency):
  [N epitopes] — coarse vs fine allele mix

RECOMMENDATIONS:
  - Flag high-frequency within-file duplicates in the DB release notes
  - Resolve coarse/fine allele mix per paper (check original publication)
  - Same-lab cross-publication duplicates: expected, keep all records
  - Independent cross-publication public clonotypes: expected, keep all records

---

Known Findings from 2026-05-31 Audit

Cross-publication duplicate patterns

Interpretation: The dominant source of cross-publication duplicates is same-lab follow-up publications. True independent replication (0 shared authors) reflects genuinely public clonotypes for immunodominant epitopes (GIL, GLCTLVAML, NLVPMVATV).

Most replicated epitopes (cross-publication)

Spurious high-frequency records

Extremely frequent within-file records (≥100 copies, ≤3 donors) reflect sequencing read depth, not individual T cell clones. Known cases:

• PMID_41315082: CASSEAGTGGYEQYF/VEALYLVCG n=530, 2 donors — high-throughput scTCR-seq
• PMID_39746936: 5 CDR3b/VISNDVCAQV combinations n=160–271, 1 donor
• PMID_34811538: Multiple CDR3b n=107–241, 2–3 donors — EBV/beta-cell antigen study
• 10xgenomics-2019-07-09: n=100–133, 2 donors — multiplexed 10x Genomics data

These are not data errors but should be noted in release documentation.

Genuine multi-MHC restriction

MHC notation fixes applied (2026-05-31)

• H-2Db/H-2Kb/H-2Kd/H-2Ld/H-2Dd → H2-Db/H2-Kb/H2-Kd/H2-Ld/H2-Dd (3,001 rows: strip hyphen between H and 2)
• H-2KB → H2-Kb (capitalization fix)
• H2 class I + SIINFEKL → H2-Kb (OVA/C57BL6 context)
• H2-b class I + SIINFEKL → H2-Kb
• H2-b class II + SIINFEKL → removed (MHC-I epitope mislabeled as class II)
• H2 class II + Ins2/SHLVEALYLVCGERG → H2-IAg7 (NOD mouse T1D context)
• H2-d class II + SFERFEIFPKE → H2-IEd (BALB/c HA restriction)