Bio Interdomain Hgt

---
name: bio-interdomain-hgt
description: Detect and polarize interdomain horizontal gene transfer with reciprocal hits, genomic context, and gene trees. Use for virus-host gene exchange, endogenous viral elements, or donor direction.
---

# Bio Interdomain HGT

Detect, polarize, and confirm horizontal gene transfer between a query genome
(virus, MAG, isolate, or bin) and other domains of life. Built for the common
asymmetric case where the query is well annotated but the comparison set is
genome-only (proteins missing). Pairs with `/bio-annotation` (homology/taxonomy),
`/bio-phylogenomics` (trees), `/bio-viromics` (viral classification), and
`/bio-fasta-database-curator` (building the arbiter database).

## Instructions

Run the steps in order; capture outputs and provenance at each step. Steps 0
(database gate) and 5 (frame-aware context guard on eukaryotic DNA) are the ones
most often skipped and most often responsible for wrong conclusions.

Run every search and tree step on compute nodes with an explicit thread count
(`diamond --threads`, `iqtree3 -T`, `mafft --thread`); on HPC, submit through
`sbatch` with CPUs matched to those counts.

After Steps 0-7 have produced their normalized TSVs, apply the evidence gates with
the driver:

```bash
uv run --script skills/bio-interdomain-hgt/scripts/run_hgt_evidence.py \
  forward_hits.tsv --arbiter-hits arbiter_hits.tsv --reciprocal reciprocal.tsv \
  --context context.tsv --trees trees.tsv --sampling-depth sampling_depth.tsv \
  --databases databases.json --hypotheses hypotheses.tsv --reflections reflections.tsv \
  --query-domain ncldv --out results/bio-interdomain-hgt
```

The driver checksum-verifies the arbiter, lineage labels, and comparison collection
listed in `databases.json`, applies the gates below, normalizes confirmed candidates
by lineage sampling depth, and writes `run_manifest.json` under
[schemas/hgt-evidence.schema.json](schemas/hgt-evidence.schema.json). Input
contracts (tab-separated, exact column order) are in `fixtures/`; the gates are:

| Gate | Pass rule in the driver |
|------|-------------------------|
| Homology | query and subject coverage >= 0.5 and e-value <= 1e-5 |
| Reciprocal best hit | `forward_rank` and `reverse_rank` both 1 |
| Direction | query's best arbiter domain equals the recipient domain: `recipient_to_query`; recipient locus best-matches `--query-domain`: `query_to_recipient`; otherwise `ambiguous` |
| Context | `recipient_domain_fraction` >= 0.6; eukaryotic recipients also need `method` `diamond_blastx` |
| Phylogeny | `nesting_clade` equals `expected_clade` and `support` >= 0.9 on a 0-1 scale |

A row is `confirmed` only when all gates pass and the direction is not ambiguous,
`candidate` when homology passes but another gate fails, and `rejected` otherwise. Domain labels use
the arbiter vocabulary `eukaryota`, `bacteria`, `archaea`, `ncldv`, `phage`,
`organelle`, and `--query-domain` must be one of them. Hypotheses need five distinct
IDs with one `technical` or `null` row; reflections must cover the gates `database`,
`forward`, `reciprocal`, `context`, `phylogeny`, `final` in that order.

### Step 0: Database availability gate (DO THIS FIRST; never hardcode paths)

HGT calls are only as good as the reference. Resolve the site/project DB root from
`$BIO_DB_ROOT` (or ask); never bake absolute paths into the analysis. Verify that
BOTH of the following exist before any search; if one is missing, build it or STOP.

1. A **comprehensive multi-domain reciprocal-arbiter proteome**: a single protein
   search database (DIAMOND `.dmnd` or MMseqs2) that contains eukaryotes + bacteria
   + archaea + viruses (including NCLDV/giant viruses and phages) + organelles,
   with a parallel `genome_id -> lineage` labels table. This one database is what
   makes "best-hit taxon", and therefore transfer direction, meaningful.
   - Building blocks: EukProt, GTDB, NCBI nr/RefSeq, IMG/VR, a giant-virus proteome
     (GVDB / gvclass-style), organelle RefSeq.
   - Check: list `$BIO_DB_ROOT` for an existing combined-proteome `.dmnd` + labels.
   - If absent: build it with `/bio-fasta-database-curator` (prefix every header by
     domain, e.g. `EUK__`, `BAC__`, `ARC__`, `NCLDV__`, `PHAGE__`, then
     `diamond makedb`). A clustered build (clusterednr / MMseqs2-reduced) is much
     faster at comparable sensitivity; prefer it.
   - A euk-only or virus-only arbiter CANNOT polarize transfer. Confirm it spans
     every candidate donor domain.
2. A **per-domain genome/proteome collection** for the comparison side (e.g. a
   eukaryote genome catalog such as EukProt/MMETSP/NCBI/Mycocosm; a viral genome
   catalog such as IMG/VR/RefSeq). Prefer one with a queryable metadata table
   (per-genome taxonomy + completeness + contamination) so hits can be quality-flagged.
   - Record whether the collection ships PROTEINS or only NUCLEOTIDES; this decides
     the forward-search tool in Step 2.

Record DB name / version / date / path-relative-to-root and per-genome counts in
the run log. If a required comprehensive DB is missing and cannot be built, say so
explicitly; do not silently substitute a non-comprehensive database.

### Step 1: Frame the query and register hypotheses
- Infer the query's domain/lineage first (`/tracking-taxonomy-updates` QuickClade
  `percontig`; `/bio-viromics` GVClass for giant viruses).
- Register >=5 working hypotheses, including technical nulls:
  1. genuine donor -> recipient HGT; 2. genuine recipient -> donor HGT /
  endogenization; 3. **assembly contamination** (a donor contig co-assembled into a
  recipient genome); 4. **deep homology / convergence** (ancient shared genes, not
  transfer); 5. **reference-sampling bias** (hits track database depth); 6.
  **virus <-> virus transfer** (a frequent confounder of apparent host-derived
  viral genes).

### Step 2: Forward search (query <-> comparison collection)
- If the comparison collection has PROTEINS: `diamond blastp` (query proteins as the
  small db, or vice versa).
- If proteins are MISSING for most of the collection: `diamond blastx` of the
  comparison NUCLEOTIDE genomes (6-frame) vs the query proteins (tiny db). For
  genome-length queries use `-F 15 --range-culling --top 10` so multiple genes per
  contig are reported.
- Scale: shard the collection across a SLURM array, bin-packed by cumulative size so
  no shard is dominated by one giant genome; set `--time` to cover the largest single
  genome; write a resume-safe per-shard `.done` sentinel.
- Thresholds: e-value <=1e-5, subject coverage >=0.5, plus identity/bitscore floors.
  Record id%, query AND subject coverage, e-value, bitscore for every hit.

### Step 3: Reciprocal classification against the arbiter
- `diamond blastp` the query proteins vs the comprehensive arbiter -> for each query
  protein, the best-hit DOMAIN and lineage (donor-derived vs query-core vs ORFan).
  Use a bitscore margin (e.g. best class must beat the next by >=10%) and coalesce
  empty-class scores to 0 before comparison (a `series.max()` on an empty group is
  NaN, and `NaN or 0` stays NaN; guard with `pd.notna`).
- For candidate recipient loci, reverse-search vs the arbiter -> best-hit domain.
- A reciprocal best hit = the query protein and the recipient locus are mutual best
  hits, with the arbiter confirming the partner domain.

### Step 4: Direction inference
- recipient <- donor (e.g. host -> virus): the query gene's best arbiter hit is the
  OTHER domain (e.g. eukaryote) and it nests within that clade.
- donor -> recipient (e.g. virus -> host / endogenization): a recipient-genome locus
  best-matches the query's domain across the whole arbiter AND sits in
  recipient-dominated genomic context (Step 5).
- Leave deep-homology / tied cases as `ambiguous` for the phylogeny to polarize.

### Step 5: Genomic-context contamination guard
- Require the recipient locus to sit on a contig dominated by the RECIPIENT domain
  (flanking genes best-match the recipient); otherwise flag as contamination or a
  free donor contig (e.g. a mis-binned NCLDV contig inside a protist MAG).
- **CRITICAL on eukaryotic genome assemblies**: do NOT call genes with a prokaryotic
  caller (Prodigal/pyrodigal): introns fragment euk genes, so the locus ORF comes
  back short and unclassifiable (in one project ~94% of loci came back blank). Instead
  use frame-aware, intron-tolerant `diamond blastx` of the locus +/- flank window vs
  the arbiter (`--range-culling --top 10 -F 15`); each HSP is a gene, classified by
  subject domain, giving both the locus origin and the flanking-gene domain mix.
  Optionally cross-check with geNomad ("is this contig viral").
- Transcriptome assemblies are ~one spliced transcript per contig, so the flanking
  context signal is weak; rely more on reciprocity + phylogeny there.

### Step 6: Deep homology vs recent transfer
- Ancient shared genes sit at LOW identity; recent HGT sits HIGH. Bound the expensive
  context + phylogeny steps to high-identity candidates (state the cutoff and log how
  many were dropped). Do not treat every conserved-core hit as HGT.

### Step 7: Per-gene phylogenetic confirmation (required for confirmed calls)
- For each top candidate, gather homologs ACROSS ALL DOMAINS from the arbiter (one
  search returning subject sequences, e.g. DIAMOND `full_sseq`), taxon-balanced and
  dereplicated; align (MAFFT) -> trim (trimAl) -> tree (IQ-TREE with ultrafast
  bootstrap, fixed seed), for example
  `iqtree3 -s aln.faa -m MFP -B 1000 --seed 1729 -keep-ident -T 4`. Pass
  `-keep-ident` so the focal tip is not collapsed; make tip names unique to avoid
  duplicate-taxon failures. Give the driver the support of the nesting node on a
  0-1 scale (UFBoot 95 becomes 0.95); `/bio-phylogenomics` `--normalize-only`
  converts a tree's labels.
- Confirmed when the focal sequence nests inside the EXPECTED donor/recipient clade
  with support. Including donor + other-virus + recipient homologs is exactly what
  separates genuine host <-> virus transfer from virus <-> virus transfer.

### Step 8: Integrate, contextualize, report
- Lineage x function matrix; transfer-direction tallies; **normalize per-lineage
  counts by collection sampling depth** (control for reference bias before claiming a
  lineage is enriched).
- Literature context (`/polars-dovmed`, `/biorxiv-search`) for the inferred group.
- Produce an interesting-findings table (evidence, confidence, comparison baseline,
  follow-up test) ordered deterministically from `hgt_candidates.tsv`: `status` first
  (`confirmed`, then `candidate`, then `rejected`), then `bitscore` descending with
  blank or non-numeric values last, then `query_protein` and `recipient_locus` ascending.

## Quick Reference

| Task | Action |
|------|--------|
| Check DBs | Confirm a comprehensive multi-domain arbiter + per-domain collection under `$BIO_DB_ROOT` (Step 0). |
| Forward search | blastp if comparison has proteins; blastx (6-frame) if genome-only. |
| Polarize | Reciprocal best hit + arbiter best-hit domain -> direction. |
| Guard | Frame-aware blastx context on euk DNA; geNomad cross-check. |
| Confirm | All-domain homolog tree; focal must nest in expected clade. |
| Tool docs | [docs/README.md](docs/README.md); DB recipe in [docs/database-availability.md](docs/database-availability.md). |

## Input Requirements
- `$BIO_DB_ROOT` set; comprehensive multi-domain arbiter `.dmnd` + labels; a
  per-domain comparison collection (proteins or nucleotides) with metadata.
- Query proteins (`.faa`); query contigs (`.fna`); optional query domain annotations.
- Tools pinned in the project's Pixi environment: diamond, mafft, trimal, iqtree, geNomad, taxonkit, seqkit (see [docs/README.md](docs/README.md)).

## Output
Driver outputs under `results/bio-interdomain-hgt/`:
- `hgt_candidates.tsv`: per hit, best arbiter domains, RBH, direction, context, phylogeny, and status
- `query_protein_origin.tsv`: donor-derived vs query-core per query protein
- `lineage_sampling_normalization.tsv`: confirmed candidates per 100 sampled genomes
- `phylogeny_evidence.tsv`, `hypothesis_register.tsv`, `gate_reflections.tsv`
- `run_manifest.json`, validated against [schemas/hgt-evidence.schema.json](schemas/hgt-evidence.schema.json)
- stdout: the last line is one JSON envelope `{ok, skill, out, manifest, warnings}` (driver stdout contract in AGENTS.md)

Workflow outputs, written by the agent:
- (agent-authored) `forward_hits.tsv`, `lineage_function_matrix.tsv`
- (agent-authored) `phylogeny/<gene>/`: alignment, tree, and nesting call
- (agent-authored) `hgt_report.md`, `logs/`

## Examples

### Example 1: Giant virus (NCLDV) query vs a eukaryote genome collection (proteins missing)
```text
Goal: HGT between an NCLDV MAG (524 proteins) and ~5,000 protist genomes.
Step 0: confirm a combined euk+bac+arc+viral+organelle arbiter .dmnd + labels under $BIO_DB_ROOT.
        protist collection ships NUCLEOTIDES only -> forward search = blastx.
Step 2: diamond blastx protist genomes (6-frame) vs the 524 viral proteins, sharded on SLURM.
Step 3: diamond blastp the 524 viral proteins vs the arbiter -> host-derived (best hit EUK) vs viral-core (best hit NCLDV).
Step 5: for high-id (>=70%) recipient loci, diamond blastx the +/-5kb window vs the arbiter (NOT pyrodigal) -> euk-dominated context?
Step 7: per-gene tree with EUK + NCLDV + other-virus homologs -> viral gene nests in a green-algal clade => host->virus HGT confirmed.
Outcome: a lineage x function HGT matrix + phylogeny-confirmed transfers, with virus<->virus alternatives ruled out.
```

### Example 2: Bacterium query vs archaeal + eukaryotic collections
```text
Same workflow; the arbiter must still contain ALL domains so a bacterial gene that
best-matches archaea (donor) can be polarized against eukaryotic and viral alternatives.
```

## Quality Gates
- [ ] Comprehensive multi-domain arbiter confirmed present (or built) and spans ALL candidate donor domains.
- [ ] The database manifest records versions and checksums for the arbiter, lineage labels, and comparison collection; every checksum is verified before candidate scoring.
- [ ] Each gate has a persisted reflection and candidate status is derived from the gates rather than assigned manually.
- [ ] Forward search direction chosen by protein availability (blastp vs blastx); coverage computed against the protein length.
- [ ] Every candidate carries id%, query+subject coverage, e-value, bitscore, and both reciprocal best hits.
- [ ] Recipient context guard used a frame-aware method on eukaryotic DNA (NOT prokaryotic gene-calling).
- [ ] Deep-homology vs recent-transfer cutoff stated; dropped count logged.
- [ ] Phylogeny includes all-domain homologs; the virus<->virus alternative is explicitly tested, not assumed away.
- [ ] Per-lineage counts normalized for reference sampling depth before enrichment claims.
- [ ] Contamination-prone hits (recipient genome with high assembly contamination, or donor-dominated contig) flagged, not silently kept.

## Non-Goals

- No confirmed transfer direction without both a reciprocal best hit and phylogenetic nesting in the expected clade. Tied and deep-homology cases stay `ambiguous`.
- No HGT call from single-domain hits. An arbiter that does not span every candidate donor domain cannot polarize anything.
- No dating of transfer events. Identity separates recent from ancient; it does not give an age.
- No lineage-enrichment claim before per-lineage counts are normalized by collection sampling depth.

## Performance gotchas (hard-won)
- `diamond blastx --sensitive` against a 100M+ protein arbiter hit 8-hour wall-clock
  limits with empty output in practice. Use default sensitivity for domain
  classification; reserve `--sensitive` for small or divergent focal sets.
- A clustered arbiter (clusterednr / MMseqs2-reduced) runs much faster; on a CUDA
  GPU node, MMseqs2-GPU `easy-taxonomy --gpu` is an alternative.
- SLURM: bin-pack by size; resume-safe `.done` sentinels; raise array throttle only
  into idle capacity; a watcher's "queue is empty" check must tolerate transient
  empty `squeue` (controller socket timeouts); require two consecutive empty reads
  before resubmitting, or you will fire duplicate arrays. Recover stragglers at finer
  granularity + longer `--time`, not by re-running everything.

## Troubleshooting
**Issue**: host_origin / recipient-locus class is blank for most loci on genome assemblies.
**Solution**: you are gene-calling eukaryotic DNA with a prokaryotic caller; switch to frame-aware `diamond blastx` of the locus window (Step 5).

**Issue**: context-guard / reverse search times out at the wall clock with little output.
**Solution**: drop `--sensitive` to default, shrink the flank window, and re-shard finely; the size of the arbiter is the cost driver.

**Issue**: apparent host-derived viral genes that may actually be virus-to-virus transfers.
**Solution**: include NCLDV + other-virus + cellular homologs in the per-gene tree and require nesting in the expected clade (Step 7).