mtnoc-rs 0.1.1

De-identified estimate of the number of distinct human maternal (mtDNA) contributors to a shotgun-metagenomic sample, via k-mer prefiltering, seed-anchored alignment, and a PhyloTree mixture-EM with bootstrap-stability.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

# mtnoc-rs

**mtNoC** — a de-identified estimate of the **number of distinct human maternal (mtDNA) contributors**
to a shotgun-metagenomic sample. It reports a conservative **lower bound** on the count of distinct
maternal lineages (macro-clades) plus their composition — aggregate, never an individual haplotype.

Pure-Rust reimplementation of the reference Python pipeline. A single static binary takes raw FASTQ in
and emits one TSV row per sample.

## Pipeline

```
raw FASTQ.gz ─▶ k-mer prefilter (rCRS 21-mers) ─▶ align to chrM ─▶ per-read alleles at
               PhyloTree-informative positions ─▶ mixture-EM + bootstrap-stability ─▶ contributor
               lower bound (macro/sub-clade) + 95% CI + lineage composition
```

Two modes:

- **`fast`** (100% pure Rust): raw FASTQ pair → embedded-rCRS k-mer prefilter → a `triple_accel`
  SIMD seed-anchored aligner (against the compact reference) → per-read alleles. No external tools.
- **`accurate`**: raw FASTQ pair → `strobealign` → `samtools` filter (MAPQ≥20, NM/qlen≤0.027) +
  `mpileup -Q20` allele extraction → reproduces the reference pipeline **exactly**. Pass a pre-aligned
  chrM BAM with `--bam` to skip strobealign. Shells out to `strobealign` + `samtools`.

Both modes feed the **same** mixture-EM, which is **bit-for-bit identical** to the reference Python
implementation (including a numpy-compatible MT19937 so the bootstrap matches exactly).

## Install / build

```sh
cargo build --release          # target/release/mtnoc
# or, once published:  cargo install mtnoc-rs
```

## Reference data

The EM reference dump and the rCRS k-mer prefilter sequence are **embedded in the binary** (gzip,
decompressed at startup), so you only supply the alignment reference:

- `mtNoC_compact.fa` — chrM (`NC_012920.1`) + NUMT-decoy contigs. Pass it via `--ref`. `fast` uses it
  for the pure-Rust aligner; `accurate` uses it as the strobealign index and `mpileup` reference.

Where to get it:

- **`git clone`**: it's in `data/mtNoC_compact.fa`.
- **`cargo install mtnoc-rs`**: the crate ships source only, so download the FASTA from the release —
  `curl -LO https://github.com/necoli1822/mtnoc-rs/releases/download/v0.1.1/mtNoC_compact.fa` — then
  pass `--ref mtNoC_compact.fa`.

`data/mtref.txt` (regenerable via `scripts/dump_ref_for_rust.py`) and `data/rcrs_wrap.fa` are the
sources for the embedded tables; you do not need them at runtime.

## Usage

```sh
# FAST (pure Rust) — one sample
mtnoc fast --in1 S_R1.fastq.gz --in2 S_R2.fastq.gz --ref data/mtNoC_compact.fa

# Shorthand: give the sample stem, mtnoc resolves _R1/_R2 (or .R1/.R2) + extension
mtnoc fast --in S --ref data/mtNoC_compact.fa

# Multiple samples (one TSV row each): comma-separate, or comma-separate stems with --in
mtnoc fast --in1 A_R1.fq.gz,B_R1.fq.gz --in2 A_R2.fq.gz,B_R2.fq.gz --ref data/mtNoC_compact.fa
mtnoc fast --in A,B --ref data/mtNoC_compact.fa

# Optionally save the chrM-candidate reads (omit to only emit the estimate)
mtnoc fast --in1 S_R1.fastq.gz --in2 S_R2.fastq.gz --ref data/mtNoC_compact.fa \
           --out1 S.cand_R1.fq --out2 S.cand_R2.fq

# ACCURATE (needs strobealign + samtools on PATH)
mtnoc accurate --in1 S_R1.fastq.gz --in2 S_R2.fastq.gz --ref data/mtNoC_compact.fa
# ...or from a pre-aligned chrM BAM (skips strobealign)
mtnoc accurate --bam S.chrM.bam --ref data/mtNoC_compact.fa
```

Common options: `--thread N` (0 = all cores), `--header` (prepend a column-name line), `--n-boot`,
`--seed`, `--min-prop`, `--min-identity` (fast). Progress goes to stderr; only the TSV goes to stdout.

## Output

Tab-separated, **one row per sample**, no header unless `--header` is given (so multiple runs
concatenate cleanly). Columns:

| column | meaning |
|---|---|
| `sample` | sample name (derived from the R1 filename) |
| `informative_reads` | reads scored at PhyloTree-informative positions |
| `contributors_ci95` | `ci95_low;min_contributors_subclade;ci95_high` (`;`-separated) |
| `composition` | `haplogroup:proportion:reads` components, `;`-separated |

Example:

```
sample           informative_reads   contributors_ci95   composition
mix_EAS85_AFR15  4000                1;1;3               M7c2b:0.9770:3736r;Z2:0.0230:88r
```

`min_contributors_subclade` is a floor limited by the bootstrap-stability criterion (≥8 supporting
reads, stable in ≥80% of resamples); the true contributor count is generally higher.

## Dependencies

- **Build**: Rust ≥ 1.74. All crates are pure Rust — no C toolchain, no `*-sys` build steps.
  gzip is `flate2` (zlib-rs backend); FASTX parsing is `needletail` (`default-features = false`).
- **Runtime**:
  - `fast` — none. A single self-contained binary.
  - `accurate` — external `strobealign` and `samtools` on `PATH` (override with `--strobealign` /
    `--samtools`). Only this mode uses external tools.

## Ethics

Reports only aggregate counts and maternal-lineage composition — never an individual's haplotype or a
re-identifiable profile.

## Attribution

Haplogroup definitions derive from PhyloTree (van Oven & Kayser); rCRS = NC_012920.1. Aligner kernels:
`triple_accel`. Please cite PhyloTree and this tool if used.

## License

MIT — see `LICENSE`. (Bundled reference data derived from PhyloTree/rCRS retains its original terms.)