|
| 1 | +# Benchmarking genotype inference |
| 2 | + |
| 3 | + |
| 4 | +The point of simulating from a *known* genotype is that you can run a |
| 5 | +genotype-inference tool on the resulting repertoire and score it against the |
| 6 | +planted truth — with no real-data uncertainty about what the right answer is. |
| 7 | + |
| 8 | +The recipe is the same for any tool: |
| 9 | + |
| 10 | +1. Build a `Genotype`, simulate a repertoire, write the AIRR table |
| 11 | + (`result.to_tsv(...)`) and/or reads FASTA, and the ground truth |
| 12 | + (`genotype.to_tsv(...)`). |
| 13 | +2. Run the inference tool to recover the per-individual allele set. |
| 14 | +3. Compare recovered vs planted: presence/absence, zygosity, and (for |
| 15 | + discovery tools) any novel alleles. |
| 16 | + |
| 17 | +### Worked example: TIgGER and IgDiscover recover a planted genotype |
| 18 | + |
| 19 | +To show this end to end we planted a diploid IGH genotype in `human_igh` — |
| 20 | +**3 heterozygous** V genes (two alleles each), **3 homozygous** (one allele), |
| 21 | +and **3 fully deleted** genes — and filled the rest from the reference. We |
| 22 | +simulated 4,000 reads with light SHM, then ran two independent AIRR |
| 23 | +genotype-inference tools on the result: |
| 24 | +[**TIgGER**](https://tigger.readthedocs.io) (Immcantation; consumes the AIRR |
| 25 | +table) and [**IgDiscover**](https://igdiscover.se) (germline discovery from the |
| 26 | +raw reads, with its own IgBLAST). |
| 27 | + |
| 28 | +Because GenAIRR already emits AIRR records with `v_call` **and** |
| 29 | +`sequence_alignment`, TIgGER's `inferGenotype` consumes the rearrangement table |
| 30 | +**directly — no separate IgBLAST step is needed**: |
| 31 | + |
| 32 | +```r |
| 33 | +library(tigger); library(airr) |
| 34 | +rep <- read_rearrangement("repertoire.tsv") # GenAIRR's AIRR output |
| 35 | +germ_v <- readIgFasta("germline_V.fasta") # cartridge V germline (names match v_call) |
| 36 | +# find_unmutated=TRUE asks inferGenotype to base calls on unmutated reads per |
| 37 | +# allele; the light-SHM simulation below provides them. |
| 38 | +geno <- inferGenotype(rep, germline_db = germ_v, find_unmutated = TRUE) |
| 39 | +plotGenotype(geno) |
| 40 | +``` |
| 41 | + |
| 42 | +The `germline_V.fasta` written below is **ungapped** — fine for `inferGenotype`; |
| 43 | +if you go on to TIgGER's IMGT-gapped steps (`reassignAlleles`) supply a gapped V |
| 44 | +germline instead. |
| 45 | + |
| 46 | +TIgGER recovered the planted genotype **exactly**: every heterozygous gene → two |
| 47 | +alleles, every homozygous gene → one, every deleted gene → **absent**. Across all |
| 48 | +52 V genes, allele-presence **precision = 1.00**, **recall = 1.00**, and the |
| 49 | +per-gene allele count matched the truth for **52/52** genes. |
| 50 | + |
| 51 | +**IgDiscover**, run on the raw reads with the cartridge as its starting database, |
| 52 | +independently agreed: **precision = 1.00** (zero false-positive alleles), |
| 53 | +**recall = 0.96** (50/52 carried alleles), with **all three deletions correct** |
| 54 | +and **all heterozygous genes fully resolved** (both alleles recovered). The two |
| 55 | +missed alleles were low-expression single-copy genes below IgDiscover's default |
| 56 | +expression threshold — a tool-tuning matter, not a simulation artefact. |
| 57 | + |
| 58 | +!!! note "These are upper-bound numbers on idealised data" |
| 59 | + Near-perfect recovery is expected here: the simulated germline names match the |
| 60 | + scoring database exactly, SHM is light and substitution-only, and there are no |
| 61 | + indels, chimeras, or contamination. Real data is harder — which is the point |
| 62 | + of being able to **dial difficulty up** (heavier SHM via `mutate`, sequencing |
| 63 | + artefacts via the corruption passes, lower per-allele depth) and re-measure |
| 64 | + how each tool degrades against the same known truth. |
| 65 | + |
| 66 | + |
| 67 | + |
| 68 | +*(A) The nine study genes: both tools' inferred allele counts match the planted |
| 69 | +zygosity for each (heterozygous → 2, homozygous → 1, deleted → 0). (B) All 52 V |
| 70 | +genes fall on the agreement diagonal; presence precision = 1.00 for both tools, |
| 71 | +recall 1.00 (TIgGER) / 0.96 (IgDiscover), zero false-positive alleles, all |
| 72 | +deletions correct.* |
| 73 | + |
| 74 | +### Reproduce it |
| 75 | + |
| 76 | +This builds the **exact** genotype behind the figure — 3 heterozygous, 3 |
| 77 | +homozygous, and 3 deleted study V genes, the rest filled from the reference — |
| 78 | +simulates 4,000 reads with light SHM at `seed=7`, and writes every input the two |
| 79 | +tools need plus the ground truth to score against: |
| 80 | + |
| 81 | +```python |
| 82 | +import GenAIRR as ga |
| 83 | +import GenAIRR.data as gdata |
| 84 | +from GenAIRR.genotype import Genotype |
| 85 | + |
| 86 | +cfg = gdata.HUMAN_IGH_OGRDB |
| 87 | +HET = ["IGHVF1-G1", "IGHVF1-G2", "IGHVF1-G3"] # 2 alleles each |
| 88 | +HOM = ["IGHVF2-G4", "IGHVF3-G5", "IGHVF3-G6"] # 1 allele |
| 89 | +DEL = ["IGHVF3-G7", "IGHVF3-G8", "IGHVF3-G9"] # deleted (both chromosomes) |
| 90 | + |
| 91 | +g = Genotype.from_dataconfig(cfg).complete_from_reference("homozygous_first_reference") |
| 92 | +for gene in HET: |
| 93 | + a0, a1 = (a.name for a in cfg.v_alleles[gene][:2]) |
| 94 | + g.heterozygous(gene, a0, a1) |
| 95 | +for gene in HOM: |
| 96 | + g.homozygous(gene, cfg.v_alleles[gene][0].name) |
| 97 | +for gene in DEL: |
| 98 | + g.delete_gene(gene, haplotype="both") |
| 99 | +g.with_subject("DONOR01") |
| 100 | + |
| 101 | +res = ( |
| 102 | + ga.Experiment.on(cfg).with_genotype(g).recombine() |
| 103 | + .mutate(rate=0.004) # light SHM, as in real data |
| 104 | + .run_records(n=4000, seed=7, expose_provenance=True) |
| 105 | +) |
| 106 | + |
| 107 | +res.to_tsv("repertoire.tsv") # AIRR table → TIgGER |
| 108 | +g.to_tsv("truth_genotype.tsv") # ground truth to score against |
| 109 | + |
| 110 | +with open("reads.fasta", "w") as fh: # raw reads → IgDiscover / partis |
| 111 | + for r in res: |
| 112 | + fh.write(f">{r['sequence_id']}\n{r['sequence'].upper()}\n") |
| 113 | + |
| 114 | +with open("germline_V.fasta", "w") as fh: # cartridge V germline (names match v_call) |
| 115 | + for gene, alleles in cfg.v_alleles.items(): |
| 116 | + for a in alleles: |
| 117 | + fh.write(f">{a.name}\n{a.ungapped_seq.upper()}\n") |
| 118 | +``` |
| 119 | + |
| 120 | +**Score it.** Run TIgGER (R snippet above) on `repertoire.tsv`, or IgDiscover on |
| 121 | +`reads.fasta` with the cartridge as its starting database |
| 122 | +(`igdiscover init --database db/ --single-reads reads.fasta project/ && cd project |
| 123 | +&& igdiscover run`). Then compare each tool's per-gene allele set against |
| 124 | +`g.to_table()` (the planted truth): allele-presence precision/recall, zygosity, |
| 125 | +and deletion calls. With the genotype above this yields TIgGER precision/recall |
| 126 | +1.00 (52/52 genes) and IgDiscover precision 1.00 / recall 0.96 — the figure. |
| 127 | + |
| 128 | +### Running other tools on the same data |
| 129 | + |
| 130 | +The only difference between tools is whether they consume the **AIRR table** |
| 131 | +(TIgGER) or the **raw reads** (`reads.fasta`, which you can write from `result`), |
| 132 | +running their own aligner: |
| 133 | + |
| 134 | +- **[IgDiscover](https://igdiscover.se)** — germline *discovery* from reads (its |
| 135 | + own IgBLAST + iterative filtering). Initialise with the cartridge germline as |
| 136 | + the starting database and the simulated reads, then run the pipeline; the |
| 137 | + `final/database/V.fasta` expressed-allele set is the recovered genotype: |
| 138 | + |
| 139 | + ```bash |
| 140 | + igdiscover init --database db/ --single-reads reads.fasta project/ |
| 141 | + cd project && igdiscover run |
| 142 | + ``` |
| 143 | + |
| 144 | +- **[partis](https://github.com/psathyrella/partis)** — HMM annotation with |
| 145 | + per-sample germline inference. Start from the cartridge germline and cache |
| 146 | + parameters into a parameter directory: |
| 147 | + |
| 148 | + ```bash |
| 149 | + partis cache-parameters --infname reads.fa \ |
| 150 | + --initial-germline-dir db/ --parameter-dir out/ |
| 151 | + ``` |
| 152 | + |
| 153 | + partis infers a per-sample germline set during this step (starting from |
| 154 | + `--initial-germline-dir`) and writes it under `--parameter-dir`; see partis's |
| 155 | + [germline-inference docs](https://github.com/psathyrella/partis/blob/main/docs/germline-inference.md) |
| 156 | + for the exact output location. Score that recovered allele set against |
| 157 | + `g.to_table()` the same way. |
| 158 | + |
| 159 | +Because the genotype is planted, every tool is scored identically: recovered |
| 160 | +allele set vs `genotype.to_table()` — presence precision/recall, zygosity, and |
| 161 | +deletion calls. |
| 162 | + |
0 commit comments