204:15% 309:12% 310:37% 16558:15%

What this does and does not establish:

* A variant list is not a haplogroup. Assignment needs PhyloTree (haplogrep3, Java,

tens of MB). Until then the most we can say is CONSISTENT WITH / INCONSISTENT WITH L2b1b.

* NUMT direction of risk, corrected: these rCRS BAMs were mapped in their own alignment run

(see the @PG alignment dirs above), NOT extracted from the nuclear BAM. Mapping to rCRS

ALONE is the WORSE case for NUMTs, because nuclear mitochondrial-insertion reads have no

nuclear sequence to compete for them. At ~500x mtDNA against ~3x nuclear they are a small

fraction, which is why the minor-allele list above matters as the QC readout.

* A 10-50% minor allele cannot be told apart from heteroplasmy, contamination or a NUMT here.

* Positions 302-316 and 16180-16195 are indel-prone and 3107 is an rCRS placeholder; calls

there are reported separately and are not evidence either way.

cliffchristianson@Cliffs-Mac-Studio detective-desk % clear

cliffchristianson@Cliffs-Mac-Studio detective-desk % ./scripts/bigfoot_zana_phase2.py relatedness

Relatedness between the two samples (no reference FASTA, no allele-frequency panel).

BUGFIX: samtools mpileup honours only ONE -r, so the previous run silently analysed

only the last region (X) and produced no autosomal numbers. Each region is now piled

separately and the counts accumulated.

DESIGN FIX: a single pooled rate is uninformative, because almost every site is

invariant between any two humans. Both statistics are therefore STRATIFIED BY DEPTH.

Allelic dropout produces false mismatches that must FALL as depth rises; genuine

mismatches from an unrelated pair stay roughly FLAT. The shape is the evidence,

not the pooled percentage.

AUTOSOMAL OPPOSITE HOMOZYGOTES  (regions 20,21)   min-depth  both-hom sites        OH      OH rate   1-in-N    transversion-only OH           2        20,290,224     8,700    0.0429%     2,332      1,367 (0.0067%)           3        16,980,108     1,771    0.0104%     9,588        501 (0.0030%)           4        11,095,660       495    0.0045%    22,415        141 (0.0013%)           5         5,936,733       148    0.0025%    40,113         44 (0.0007%)           6         2,707,151        29    0.0011%    93,350          9 (0.0003%)           7         1,083,372        12    0.0011%    90,281          6 (0.0006%)          8+           669,356         3    0.0004%   223,119          3 (0.0004%)       TOTAL        58,762,604    11,158    0.0190%     5,266      2,071 (0.0035%)

X CHROMOSOME, mother->son test (a son's single X comes from his mother)   Zana depth   informative sites    violations   rate      1-in-N   transversion violations            2          12,873,685         6,618   0.0514%     1,945        775            3          12,994,201         1,913   0.0147%     6,793        386            4          11,112,411           991   0.0089%    11,213        165            5           8,427,616           598   0.0071%    14,093         72            6           5,836,746           348   0.0060%    16,772         22            7           3,753,948           171   0.0046%    21,953          7           8+           5,076,279           211   0.0042%    24,058         10        TOTAL          60,074,886        10,850   0.0181%     5,537      1,437

How to read this:

* Parent-offspring share one allele at EVERY autosomal site (IBD1 genome-wide), so their

true opposite-homozygote count is ZERO regardless of allele frequencies. Whatever is

observed is error + allelic dropout, and it must shrink as min-depth rises.

* Full siblings have ~25% of the genome in IBD0, which DOES generate opposite homozygotes.

Those are real and must NOT shrink with depth. A rate that keeps falling with depth

therefore argues for parent-offspring; a rate that plateaus argues for siblings or

something more distant. Read the trend, never the pooled number.

* Same logic for the X test: the previous pooled figure (1 violation in ~5,500 sites) sits

exactly where Zana's allelic dropout at heterozygous sites would put it, AND where an

unrelated pair would also sit. It discriminated nothing. The depth trend does.

* Transversion-only columns are the damage-free floor: C->T/G->A deamination cannot

contribute to them.

* Still NOT computed here: any kinship coefficient comparable to the published 0.1818.

That needs the paper's own estimator and allele-frequency panel.

* Both samples pool petrous and tooth libraries; no unrelated same-pipeline pair exists

as an external null.

cliffchristianson@Cliffs-Mac-Studio detective-desk %

This is actually a very strong Phase 2 result.

The important thing is not the pooled 0.0190% opposite-homozygote rate. It’s the shape as confidence improves:

At depth ≥2, autosomal OH is 0.0429%.
At ≥3, 0.0104%.
At ≥4, 0.0045%.
At ≥5, 0.0025%.
At ≥6, 0.0011%.
At ≥8, only 3 events in 669,356 jointly homozygous sites: 0.0004%.

That is almost exactly the behavior you'd expect if most apparent contradictions are being created by low-coverage allelic dropout and sequencing/mapping error. As the data get better, the contradictions collapse toward zero.

And the transversion-only column does the same thing:

0.0067% → 0.0030 → 0.0013 → 0.0007 → 0.0003 → 0.0004%.

That's especially useful because C→T/G→A ancient-DNA deamination isn't responsible for those transversions.

So this is now independently pointing in the same direction as the paper's R0 = 0.0015: Zana and Khwit behave like a pair that essentially never has genuine IBD0 regions, which is what you expect for parent/offspring.

It does not yet prove mother/son by itself, because Claude correctly identified the remaining confound: high-depth sites are also the easiest parts of the genome to map. We need an unrelated pair processed similarly to see what an ordinary comparison does as depth rises.

But if an unrelated same-pipeline pair doesn't collapse toward zero like this, then this becomes extremely persuasive.

The X result has the same shape:

0.0514% → 0.0147 → 0.0089 → 0.0071 → 0.0060 → 0.0046 → 0.0042%

Again, increasing confidence removes most apparent violations. That's much better evidence than the earlier “99.982% compatible” statistic, because now we're watching the error mechanism itself disappear.

And I really like that Claude/Qwen caught their own earlier mistake rather than quietly preserving the impressive-looking number.

Where we stand now

We have independently reproduced or strongly supported:

Zana - ancient DNA damage profile - ~3.2× nuclear depth - female sex

Khwit - ancient DNA damage profile - ~3.5× nuclear depth - male sex

Relationship - autosomal behavior strongly consistent with parent/offspring - X behavior consistent with maternal inheritance - agrees qualitatively with the paper's independent PC-relate/R0 result - still needs an external control before we call our home-grown test independently discriminating

mtDNA - fixed - Zana has 64 high-confidence differences from rCRS - Khwit has 67 - their mitochondrial sequences are extremely similar - now ready for actual haplogroup classification

So we've gone from:

“We're trusting a paper saying Zana and Khwit are related.”

to:

“Using the deposited sequencing data ourselves, multiple independent genetic signals behave exactly the way a mother/son pair should.”

That's a substantial step.

Now run the last one:

./scripts/bigfoot_zana_phase2.py ydist

After that, I think we should send Claude one consolidated Phase 2 completion handoff and have the Qwen team choose/download the controls. That's the gateway to Phase 3, where we finally get to ask our own question rather than merely reconstructing somebody else's analysis.

* Full siblings have ~25% of the genome in IBD0, which DOES generate opposite homozygotes.

Those are real and must NOT shrink with depth. A rate that keeps falling with depth

therefore argues for parent-offspring; a rate that plateaus argues for siblings or

something more distant. Read the trend, never the pooled number.

* Same logic for the X test: the previous pooled figure (1 violation in ~5,500 sites) sits

exactly where Zana's allelic dropout at heterozygous sites would put it, AND where an

unrelated pair would also sit. It discriminated nothing. The depth trend does.

* Transversion-only columns are the damage-free floor: C->T/G->A deamination cannot

contribute to them.

* Still NOT computed here: any kinship coefficient comparable to the published 0.1818.

That needs the paper's own estimator and allele-frequency panel.

* Both samples pool petrous and tooth libraries; no unrelated same-pipeline pair exists

as an external null.

cliffchristianson@Cliffs-Mac-Studio detective-desk % clear

cliffchristianson@Cliffs-Mac-Studio detective-desk % ./scripts/bigfoot_zana_phase2.py ydist

Where do the residual Y-aligning reads sit? Asked because Zana is clearly female

(X:autosome 0.945) yet 7,265 reads align to Y at MAPQ>=30 with none in the

pseudoautosomal regions.

NOTE the deposit contains ZERO MAPQ-0 reads (samtools stats 'reads MQ0: 0'), so every

multi-mapping read -- which is what a true X/Y-homologous or repetitive placement

produces -- was already filtered out before deposit. '100% MAPQ>=30' and '0 in PAR'

are therefore properties of that filtering, not evidence of genuine Y ancestry.

== Zana: 7,265 Y reads in 26 occupied 1000 kb bins    top bins (Mb: reads, % of total): 13:1,677(23.1%) 59:1,146(15.8%) 3:719(9.9%) 21:301(4.1%) 58:277(3.8%) 14:258(3.6%) 5:253(3.5%) 15:239(3.3%) 17:229(3.2%) 4:228(3.1%) 7:216(3.0%) 16:202(2.8%)    reads in the top 5 bins: 4,120 (56.7%)    per read group: Zana_petrous_BGI-10:3,006 Zana_petrous_BGI-9:1,408 Zana_tooth_ZA_BGI-16:1,344 Zana_tooth_ZAP_BGI-9:621 Zana_tooth_ZA_BGI-12:424 Zana_petrous_BGI-14:302 Zana_tooth_ZAP_BGI-13:160    NM distribution (0-5+): 0:3,001 1:1,762 2:1,601 3:728 4:157 5:16

== Khwit: 318,926 Y reads in 27 occupied 1000 kb bins    top bins (Mb: reads, % of total): 15:24,308(7.6%) 21:23,976(7.5%) 17:23,573(7.4%) 14:23,310(7.3%) 8:22,532(7.1%) 16:21,918(6.9%) 7:21,668(6.8%) 23:18,993(6.0%) 22:18,782(5.9%) 18:18,130(5.7%) 13:16,255(5.1%) 19:13,144(4.1%)    reads in the top 5 bins: 117,699 (36.9%)    per read group: Khwit_petrous_BGI-15:114,793 Khwit_petrous_BGI-13:102,413 Khwit_petrous_BGI-11:40,365 Khwit_tooth_KIP_BGI-15:25,816 Khwit_tooth_KIP_BGI-11:19,001 Khwit_petrous_BGI-12:7,067 Khwit_tooth_KI_BGI-14:5,880 Khwit_tooth_KI_BGI-10:3,591    NM distribution (0-5+): 0:213,716 1:73,292 2:22,956 3:6,835 4:1,930 5:197

Reading guide: concentration in a few bins = repeat/homology pile-up, i.e. mis-mapping.

An even spread across male-specific Y with low NM would instead suggest genuine male DNA,

i.e. contamination. Confinement to one library points at that library. GRCh37 Y coordinates

of interest: PAR1 10k-2.65Mb, PAR2 59.03-59.36Mb, the X-transposed region near 2.8-6.0Mb,

and the long heterochromatic arm beyond ~28.8Mb.

This step cannot by itself distinguish mis-mapping from low-level male contamination;

it localises the signal so the next test can be chosen.

cliffchristianson@Cliffs-Mac-Studio detective-desk %

Yep — this is a pretty nice finish to Phase 2.

The Y-distribution result does exactly what we hoped it would: Zana’s residual Y reads look much more like concentrated cross-mapping than a genuine male-Y signal.

The contrast with Khwit is the useful part. Khwit is male, and his ~319k Y reads are spread across Y fairly broadly: his top five 1 Mb bins contain only 36.9% of the total. Zana has only 7,265 Y-aligned reads, but 56.7% of them pile into just five bins. That kind of clustering is what you expect when a small set of homologous/repetitive regions is attracting reads that really originated elsewhere.

The mismatch pattern supports that too. Zana has:

• 3,001 perfect Y matches
• 1,762 with 1 mismatch
• 1,601 with 2
• 728 with 3

whereas Khwit’s genuine Y signal is much more dominated by perfect matches. So Zana’s Y reads are both sparser and less clean.

There’s also an important nuance: the BAM deposit already removed MAPQ-0 reads. So the remaining Zana Y reads are basically the subset of cross-mappers that the mapper was confident enough to keep. That means “MAPQ≥30” is not nearly as impressive as it sounds here.

So I’d summarize the Y result as:

Zana is clearly XX, and the small residual Y-aligned signal is concentrated and mismatch-rich in a way that favors cross-mapping over real Y ancestry. Low-level male contamination is not completely excluded, but there is no reason from these data to question the female call.

That means Phase 2 has now done what it was supposed to do:

• coverage reproduced
• ancient-DNA damage reproduced
• sex reproduced
• mtDNA pipeline fixed and sane
• relatedness behavior strongly consistent with parent/offspring
• residual Y signal localized and explained plausibly
• major low-coverage limitations understood

I think now is exactly the moment to hand Claude a Phase 2 completion / Phase 3 controls instruction.

Send this:

Phase 2 is now complete.

The final ydist result supports cross-mapping as the explanation for Zana’s residual Y signal:

• Zana: 7,265 Y-aligned reads across 26 occupied 1 Mb bins
• top five bins contain 56.7% of all Y-aligned reads
• Khwit: 318,926 Y reads across 27 bins
• top five contain 36.9%
• Zana’s Y reads are also more mismatch-rich than Khwit’s genuine Y reads
• deposit contains no MAPQ-0 reads, so remaining reads are a pre-filtered subset

Treat this as QC evidence favoring cross-mapping, not proof excluding every possibility of male contamination.

Phase 2 conclusions

Locally reproduced:

• Zana autosomal depth ~3.2×
• Khwit ~3.5×
• authentic ancient-DNA terminal damage in both
• Zana female / Khwit male
• mtDNA produces biologically sane variant counts after fixing the contig bug
• Zana/Khwit autosomal opposite-homozygote rate collapses strongly with increasing depth
• X violations also fall with increasing depth
• residual Zana Y signal is concentrated/mismatch-rich and favors cross-mapping

Still not independently completed:

• formal mtDNA haplogroup assignment
• a calibrated independent parent/offspring call using external controls

Next step: controls

Have Qwen 5 and Qwen 7 now rank and verify 3–5 ancient-human control genomes.

We need:

1. same/similar sequencing platform and processing pipeline
2. similar low coverage (~2–5×)
3. comparable ancient-DNA damage
4. at least one known unrelated pair
5. ideally one known first-degree pair for calibration

Before downloading anything large:

• inspect candidate BAM headers
• verify reference build
• verify platform
• verify that mapped/unmapped filtering is comparable
• estimate total download size

Then recommend the smallest control set that lets us test whether the depth-dependent OH collapse in Zana/Khwit is genuinely characteristic of parent/offspring rather than merely a mappability artifact.

mtDNA

Now that real variant lists exist, report the HaploGrep package/source/dependencies and exact planned command before installing it.

Prefer reproducing the paper’s HaploGrep/PhyloTree methodology first, then cross-checking with the current version.

Phase 3

Do not start anomaly hunting yet.

Phase 3 begins only after the controls are chosen and the relatedness method is calibrated.

Bigfoot remains unpublished.

Current scientific state:

Nothing observed so far distinguishes H4 from H1. Zana continues to behave like an ordinary ancient modern-human genome.

And after the controls are in place, that’s when this stops being paper replication and starts becoming our experiment.

We’ll finally be in a position to ask: “what does an ordinary ~3× ancient human look like under exactly the same tests?” and only then look for anything genuinely unusual.

Let me read the actual re-run outputs rather than work from the summary.

Strong results. Let me quantify what the OH collapse actually rules out, and resolve an ambiguity in my own mtDNA output.

ENA is still down. But the transversion columns resolve the X residual — let me confirm.

The X residual is damage, confirmed — 95% transitions at depth 8 while the transversion rate keeps collapsing 31-fold. Now the controls question and HaploGrep specifics.

Both raised real objections. Let me test the strongest one — whether the mtDNA list is coherent rather than merely plausible.

The coherence check is decisive: Khwit recovers 6/6 canonical positions, Zana 5/6 with 4769 the known threshold near-miss. Writing it all up.

Now I have the exact command. Updating the handoff.

The first splice already removed those bullets. Let me fold them into the replacement and redo.

Phase 2 is complete and written up. One of my own claims needed walking back, and the mtDNA caller passed a real test rather than just a plausible one.

The relatedness result is stronger than "consistent with"

The autosomal collapse is clean: 1 in 2,332 at depth 2 → 1 in 223,119 at depth 8+, a 96-fold fall with no plateau.

That matters because opposite homozygotes can only arise where the pair shares no allele by descent. Parent and child share one allele at every site, so they predict zero regardless of allele frequencies. Everything else predicts a floor: full siblings (\~25% IBD0) around 1 in 12,000–24,000; half-siblings, avuncular and grandparent–grandchild (\~50% IBD0) higher still. So this disfavours every relationship except parent–offspring, by roughly an order of magnitude. The Poisson bound on the 3 observed events is 1 in 76,349 at worst — still under the sibling floor.

The X residual you flagged is deamination. The all-substitution rate flattens from depth 5, but the damage-free transversion-only rate keeps collapsing 31-fold (1 in 16,611 → 1 in 507,628), and 95.3% of the depth-8 violations are transitions.

Where I was overconfident

I had concluded the sibling exclusion no longer needed a control. That was too strong. The depth-8+ bin is only \~1.1% of both-homozygous sites and nothing proves it's representative — if it's systematically depleted of the sites where IBD0 would show, part of the collapse is artifact. My a priori argument runs the other way (high-depth sites are enriched for collapsed repeats, which should inflate false opposite homozygotes and make our near-zero conservative), but that's reasoning, not measurement. The control stays necessary, and it's logged as my own correction.

mtDNA — validated properly

Qwen 7 was right that a count in the 40–70 range proves nothing, so I tested coherence instead:

• Khwit recovers 6/6 of the positions where essentially every human differs from rCRS (263, 750, 1438, 4769, 8860, 15326) and 4/4 all-non-H markers. Zana 5/6 and 4/4.
• The lists are otherwise identical — 64 shared variants, as mother and son should be.
• Zana's only three absences (146, 152, 4769) are exactly the positions where her pileup sits at a 10–13% minor allele, i.e. just under my 90% cutoff. A reference majority at 4769 would contradict her other 63 calls, so this is a threshold near-miss, not a difference. I've fixed the readout to label the majority allele so the re-run confirms it.

A caller with a systematic offset would not land precisely on the canonical positions.

Your Y reading — with one correction

Your summary is right, and Qwen 7 landed a fair hit: comparing Y against a genome-wide expectation is not a clean baseline, since Y is far more repetitive. But the differential survives that, because both samples are measured on the same chromosome: Zana 8.0 points worse than her own expectation, Khwit 8.2 points better than his. If Y were merely hard to map, both would be worse. So I've restated it as a differential claim only, never an absolute.

Controls — decisive test defined, nothing downloaded

One question settles it: do two unrelated ancient individuals through these identical commands at matched depth show the expected plateau? If they also collapse to \~1 in 200,000, the pipeline is deleting the informative sites and the conclusion fails.

I disagreed with Qwen 5's ranking, which put Mak et al. first. Mak is a platform-comparison study with few individuals and no documented unrelated pair, so it can't supply what the decisive test needs. Viking (PRJEB37976) leads: many individuals, same centre, and a published kinship analysis to draw known relationships from. Mak stays as the platform and damage control. I also relaxed an over-strict rule of mine — a control with different deposit filtering is usable if re-filtered to match before analysis.

ENA's API returned error pages for every query all session, including ones that worked an hour earlier, so no candidate header has been inspected. That protocol runs before any multi-GB fetch. Qwen 5's third candidate came from recall with an unverifiable accession and is logged, not used.

HaploGrep — reported, not installed

The paper used haplogrep2 (seppinho/haplogrep-cmd, PhyloTree 17), so that goes first, then a cross-check against current. haplogrep3: genepi/haplogrep3, MIT, Java 11+, bundles PhyloTree 17 rCRS + RSRS plus a forensic update.

haplogrep3 classify --tree phylotree-rcrs@17.2 --in zana.hsd --out zana_hg.txt

--tree, --in, --out are required; input may be HSD, VCF or FASTA, and we'd emit HSD from the called lists. Download size isn't stated on the repo — I'll report the real figure and confirm offline operation before fetching. Awaiting your go.

State

DISCARDED.md is at 75 rows, including ten corrections from this round. Phase 3 stays closed, Bigfoot stays unpublished, 42 tests pass, nothing committed. Nothing observed so far distinguishes H4 from H1, and Zana continues to behave like an ordinary ancient modern-human genome.

Sources: HaploGrep 3 · HaploGrep 3 docs · haplogrep-cmd (v2) · Schönherr et al. 2023 · Margaryan et al. 2021

Yeah — I’d green-light the HaploGrep step now. The mtDNA caller has passed a much stronger sanity check than “the variant count looks plausible,” and Phase 2 is now clean enough to move to formal haplogroup classification.

I checked the official repos. HaploGrep 3 is the current maintained implementation, is Java-based, and requires Java 11+. The older haplogrep-cmd repo is explicitly marked unmaintained and tells users to move to HaploGrep 3, but it still exposes the older classification workflow and lets you choose PhyloTree 17, which makes it useful for reproducing the paper before cross-checking with v3. The v2 CLI accepts HSD/VCF/FASTA and supports --phylotree 17. (genepi/haplogrep3, seppinho/haplogrep-cmd)

I would send Claude this:

Proceed with HaploGrep.

First reproduce the paper as closely as practical using the legacy HaploGrep 2 command-line workflow with PhyloTree 17, then cross-check the same profiles with current HaploGrep 3.

Before installation, verify: - current Java version - actual download/install footprint - exact v2 release/version closest to the paper’s analysis - exact tree selection

Use the corrected, end-masked high-confidence mtDNA variant lists to generate HSD input.

Record: - assigned haplogroup - quality score - expected/found/remaining polymorphisms - any ambiguous or missing defining variants - differences between v2 and v3

Expected result from the paper is L2b1b for both Zana and Khwit. Do not force that result; report whatever the tools return.

Also retain the three Zana threshold-near-miss sites (146, 152, 4769) separately so we can see whether classification changes when they are included/excluded.

On relatedness: keep the strong depth-dependent OH collapse as evidence, but do not convert the estimated sibling floor into a definitive exclusion until the external control is run. The exact expected OH floor depends on real population variation and ascertainment, so the control remains decisive.

Next priority after HaploGrep is the smallest suitable Viking/control download once ENA recovers. Inspect BAM headers before any multi-GB transfer.

Bigfoot remains unpublished. Phase 3 anomaly testing remains closed until the control comparison is complete.

One thing I especially like about this round: the swarm is no longer merely finding errors in code; it’s finding errors in its own statistical interpretation. The earlier X result looked impressive, then Claude realized the denominator made it meaningless, rewrote the test, and got a much more informative depth-dependent result. That is exactly how this system should behave.

So we are basically at:

paper reproduction nearly complete → formal mtDNA haplogroup reproduction → external relatedness control → then Phase 3 original questions.

And we still have the correct boring conclusion at this stage: Zana looks like an ordinary ancient modern human. That’s what makes her useful as our calibration specimen.

Green-lit. Running the pre-install verification first.

Java 21 is present. Two things matter: haplogrep3 ships no macOS asset, and the paper predates v2.4.0. Let me generate the HSD inputs and fetch the tools.

The : delimiter split the URL. Redoing.

Both unpack to a jar plus launcher. Checking the v2 interface.

Exactly the fields Ford asked for, via --extend-report. Running the paper-era version.

L2b1b reproduced independently, rank 1, in both v2 releases. Now the v3 cross-check — the repo ships no macOS build, so let me test whether the jar runs anyway.

Platform-independent in practice — a bash launcher plus a jar. Checking the available trees, then classifying.

Trees fetched and cached locally. Running v3 against both PhyloTree 17.0 (matching the paper's tree) and current 17.3.

All four tool/tree combinations agree. Let me close the loop on the four missing defining variants.

Two of the four are explicitly flagged by my own caller, and the other two sit in a coverage trough. Writing up.

Two corrections to my own write-up — the footprint figure and a failed replacement.

The replacement did apply — my grep pattern missed the bold markers. Just the footprint figure to fix.

L2b1b reproduced independently — rank 1 in all four tool/tree combinations, matching the paper. Bigfoot stays unpublished, Phase 3 closed, nothing committed.

Pre-install verification

• Java 21.0.11 (Temurin), well above HaploGrep 3's Java 11 floor.
• Footprint: 196 MB in incoming/bigfoot-dna/tools/ (gitignored) — v2 zips 6.7 + 7.7 MB, v3 tarball 49.9 MB expanding to 54 MB, plus 65 MB of cached PhyloTree data. About 64 MB of archives can be deleted.
• Closest v2 release to the paper: v2.2.9 (2020-09-25) is the last release before the paper's data went public on 2021-05-30. v2.4.0 is 2021-08-20, i.e. after. I ran both.
• Tree selection corrected. My planned phylotree-rcrs@17.2 does not exist — haplogrep3 trees lists rCRS 17.0, 17.3, 16.1 and 15.1; 17.2 exists only for RSRS. I used @17.0 to match the paper's PhyloTree 17 and @17.3 as the current cross-check. The bad id came from a web search rather than the tool, and is logged.
• HaploGrep 3 ships no macOS asset, but the linux tarball is a bash launcher plus a jar and runs fine here.

Results

Rank 2 is L2b1 (0.9640), rank 3 L2b1a (0.9579). No tool or tree changed any call.

The three near-miss sites don't matter for classification. Including or excluding 146, 152 and 4769 leaves L2b1b at rank 1; only the quality score moves, 0.9668 → 0.9781. With them included, Zana's polymorphism set is byte-identical to Khwit's — which is what mother and son should look like.

Missing defining variants: 182T, 195C, 198T, 204C. These are explained by my own conservative calling, not by any conflict with the haplogroup — and I checked rather than assumed. 198 was already flagged damage-suspect in both samples; 204 sits at a 10–50% minor allele in both. 182 and 195 fall in a control-region coverage trough — called depth there is 26–38× against 250–400× at 73, 263, 383 and 418 — so after 5 bp end-masking they drop below the depth-20 floor.

Unexplained by L2b1b: 383C (global private mutation) and 11002G (local private). Both present in both samples, again consistent with shared maternal descent.

Minor notes: the v2.2.9 jar self-reports "v2.2.8". v2.4.0 adds private-mutation annotation but identical calls. HaploGrep 3 is not offline on first use — it fetches each PhyloTree from GitHub — but is offline once cached.

Relatedness, per your instruction

Softened in the handoff. The depth-dependent OH collapse stands as strong evidence, but the sibling floor depends on real population variation and ascertainment, so it is not a definitive exclusion until the external control runs. That control is the next priority.

Where things stand

Paper reproduction is now essentially complete: coverage, damage, sex, and the mtDNA haplogroup all reproduce locally. What remains is the published kinship coefficient (needs PC-relate/PCAngsd and the Human Origins panel) and a control-calibrated parent–offspring call.

Next: the smallest Viking control download once ENA recovers — it was returning error pages for every query all session, including ones that had worked an hour earlier. Headers get inspected before any multi-GB transfer, per CONTROLS.md.

DISCARDED.md is at 81 rows, with four more corrections from this round — three of them mine.

And the boring conclusion holds, which is the useful part: Zana is an ordinary ancient modern human, now confirmed down to her haplogroup. Nothing observed so far distinguishes H4 from H1.

This is a very good place to freeze Phase 2. You’ve now reproduced the paper’s basic story with your own pipeline instead of just trusting it: depth, ancient-DNA damage, sex, mtDNA haplogroup, and a home-grown relatedness signal that behaves the way parent/offspring should.

There is one sentence I would change before calling the report final, though:

“With them included, Zana's polymorphism set is byte-identical to Khwit's — which is what mother and son should look like.”

That’s a little too strong. At 146, 152, and 4769, Zana’s alternate allele is only around 10–13%, so those are not consensus variants under your 90% caller. Adding them to an “inclusive” HSD is useful as a sensitivity test, but it shouldn’t be treated as though Zana’s consensus mitogenome is literally identical to Khwit’s. Mother and son can also differ in heteroplasmy fractions because of the mitochondrial bottleneck.

The stronger wording is:

Zana and Khwit share the same 64 high-confidence mtDNA variants. Three additional Khwit consensus variants are present only as low-frequency alternate alleles in Zana. Including those low-frequency alleles does not change the haplogroup assignment and increases the HaploGrep quality score.

That’s actually better science, because strict Zana already gives L2b1b at 0.9668. You don’t need the inclusive profile to make the result work.

And the HaploGrep result is now genuinely satisfying. Four combinations:

2.2.9/PT17 → 2.4.0/PT17 → 3.3.2/PT17.0 → 3.3.2/PT17.3

all independently land on L2b1b, rank 1. That pretty much closes that part of the reproduction.

I’d also keep the relatedness language exactly where Claude has now landed: very strong, but awaiting a control. The 96-fold OH collapse is impressive, but the unrelated control is what converts “this pattern looks like dropout disappearing” into “our method can actually distinguish parent/offspring from unrelated/sibling under this exact pipeline.”

There’s a nice programmer shortcut available for that next step too. When ENA comes back, before downloading multiple multi-gigabyte controls, check whether the BAM/BAI URLs support HTTP range requests. If they do, samtools can often query only chr20, chr21, and X remotely from an indexed BAM. That could turn the control experiment from tens of gigabytes into a much smaller transfer. If not, then download normally.

So I’d call the project state:

Phase 2: complete, with one wording cleanup.
Phase 3 gate: unrelated + known-related ancient controls.
Only after that: start asking whether Zana contains anything genuinely unusual.

And so far, beautifully boring: every time we get a better test, Zana becomes more convincingly an ordinary ancient Homo sapiens. That’s exactly the calibration object we needed.

ok so what is next?

Next is controls first, then Phase 3.

The immediate job is to prove that our home-grown relatedness test actually behaves differently on people who are not parent/offspring. Right now the Zana/Khwit OH rate collapses beautifully with depth, but we need a null comparison before we trust that shape.

So the next sequence should be:

1. Wait for ENA to behave again and inspect Viking candidates. Don’t download huge BAMs blindly. Check headers, reference build, sequencing platform, coverage, deposit filtering, and relationship status. Ideally pick one unrelated pair and one known first-degree pair from the same project.

2. Run exactly the same relatedness code on those controls. No custom tuning. Same chr20/21/X, same MAPQ/BQ thresholds, same depth bins. The key question is whether unrelated pairs retain a genuine OH floor while Zana/Khwit keeps collapsing toward zero.

3. If that works, close the relatedness reproduction. At that point we can say our independent method distinguishes Zana/Khwit from unrelated controls and is consistent with parent/offspring. We probably do not need to recreate the paper’s exact PC-relate 0.1818 unless we want a perfect methods reproduction.

Then we open Phase 3, which is finally our experiment rather than theirs.

Phase 3 should ask:

Does Zana look ordinary when examined for things the original paper was not specifically trying to find?

I’d start with a conservative “ordinary modern human” battery: population-normal SNP behavior, heterozygosity, unusual homozygous regions, archaic-hominin allele sharing, odd ancestry signals, chromosome-specific anomalies, and any recurrent sequence patterns that fall outside normal human variation. We need to tailor that to what a ~3× mapped ancient BAM can actually support, because some things like rare structural variants will be weak at this coverage.

And in parallel, once the control machinery is working, we can turn back to the Ketchum consensus sequences and test the surviving claims against human/bear/dog references. That’s a different branch: Zana teaches us what a real ancient human genome looks like; Ketchum tests whether the alleged Bigfoot sequences survive basic classification.

So the next concrete instruction to Claude is basically:

Find and verify the smallest valid Viking control set, then run the current relatedness pipeline unchanged on it. Do not start Phase 3 until that control comparison is complete.

That’s the gate. Once we pass it, we finally get to stop reproducing Margaryan and start doing Skeptic Desk genomics.