lol haplogrep, so this is strong programmer terrain we are walking on here

Oh yeah 😄 We are deep in programmer territory, just pointed at biology instead of Laravel.

A lot of modern genomics is basically:

huge structured files + weird legacy formats + command-line tools + statistical assumptions + pipelines + reference data + edge cases + “why the hell did this parser return zero?”

That is extremely familiar terrain for you.

HaploGrep itself is basically a specialized classifier: you give it mitochondrial variants, it compares them against a phylogenetic tree, scores candidate haplogroups, and tells you the best fit. That’s not conceptually that different from taking an input record and running it through a rules engine backed by a hierarchical dataset.

The hard part isn’t “can we run the software?” You absolutely can.

The hard part is understanding whether the input data and assumptions are valid.

That’s exactly where programmers are useful because genomics has tons of traps that look like software bugs:

• wrong contig name
• wrong reference build
• reference mismatch
• bad coordinate system
• missing index
• filtering everything by accident
• tool expects phased genotypes but gets low-coverage calls
• one pipeline silently drops unmapped reads
• another normalizes indels differently
• BAM header says one thing, actual data imply another
• “0 variants” because a flag was wrong

Honestly, the mtDNA issue we just hit is the perfect example. A biologist might see “0 variants” and think something biological is strange. A programmer instinct is:

“Nope. That smells like a pipeline bug.”

And that instinct is probably right.

So yeah — genomics is much more software-heavy than most people realize. Once the wet-lab side has produced the data, a huge amount of the science is basically sophisticated data engineering and statistical computing.

You are not walking into alien territory nearly as much as it feels like. You’re walking into a new domain with very familiar machinery.

so explain to me where we are right now?

Right now we’re at a pretty clean transition point: we have finished proving that our local genomics pipeline can reproduce and interrogate the Zana/Khwit data, but we have not yet started the actual “is there anything unusual about Zana?” analysis.

We started with the broad Bigfoot hypothesis: maybe some alleged Bigfoot samples come back “human” because the population is a very close Homo lineage that can interbreed with us, with much of the visible difference coming from developmental/regulatory biology rather than a wildly different genome. Zana is our control case for learning how to test that idea correctly, not evidence for Bigfoot by herself.

So far we have done real independent work rather than merely trusting the paper. We downloaded roughly 20 GB of public Zana/Khwit sequencing data from ENA, verified every file against its MD5, installed samtools/bcftools, and reproduced the basic coverage results. Zana comes out around 3.20× nuclear coverage, Khwit around 3.48×, very close to the published ~3.1×/~3.3× numbers.

We then moved into Phase 2 and found several things that are working extremely well:

• Ancient-DNA authenticity signal: Zana has the expected strong 5′ C→T and 3′ G→A damage that decays away from the read ends. Khwit shows the same pattern more weakly. That is exactly the sort of degradation pattern ancient DNA should show.
• Sex: independently reproduced very cleanly. Zana’s X/autosome ratio is about 0.945, strongly XX/female; Khwit’s is about 0.498, essentially textbook XY/male.
• Coverage distribution: we confirmed that “3× genome” is deceptive. A lot of positions still have only 0, 1, or 2 reads, which means allelic dropout is a serious limitation. We now understand how cautious we have to be with rare variants.
• Maternal relationship signal: on Khwit’s X chromosome, about 99.982% of usable alleles are compatible with being present in Zana. That is a very strong sanity check for the claimed maternal relationship.
• Ancient-data limitations: the deposited BAM files contain no unmapped reads, so anything discarded before alignment is simply unavailable to us. That is important because it means we cannot use these BAMs to search for exotic sequence that failed to map to the human reference.

We also caught bugs, which is exactly what we wanted the swarm to do. The biggest one is the mitochondrial analysis: the script currently says zero variants versus the human mitochondrial reference for both Zana and Khwit. That makes no sense alongside an L2b1b haplogroup claim, so we are treating that as a software/pipeline problem, not a biological result. We also caught a documentation error where the script said -q/-Q trimmed damaged read ends; they actually control mapping/base quality, not terminal trimming.

There are three other open technical questions. The script advertised autosomal opposite-homozygote relatedness numbers, but those numbers didn’t show up in the completed output, so we need to find out whether that calculation failed or was simply omitted. The paper’s 0.1818 kinship coefficient also needs to be traced back to the exact estimator before we compare it with a textbook 0.25 parent/child expectation. And Zana has a small number of Y-aligning reads that look high-quality; her X ratio still makes her clearly female, but those reads are worth mapping to see whether they’re cross-mapping, contamination, or something else mundane.

Parallel to Zana, we recovered the surviving Ketchum Bigfoot DNA archive from the Wayback Machine: paper, methods, tables, primers, correspondence, and three nuclear consensus sequences. Crucially, those are not raw sequencing reads even though the paper apparently calls them “raw sequence data.” Without the reads, we cannot reconstruct contamination, mixtures, heterozygosity, or read support. But we can later test the surviving consensus sequences against human/bear/dog references and see whether Hart’s assignments reproduce.

On the software side, Skeptic Desk itself is also much stronger now. The eight Qwens are functioning as a team, Claude is reviewing them, bad ideas go into DISCARDED.md, draft investigations are hidden, automated updates can be held as proposals, weekly synthesis has an empty-output guard, and the Bigfoot investigation is still unpublished. Nothing observed so far distinguishes our interesting H4 relic-Homo/introgression idea from the ordinary H1 explanation.

So the immediate next step is not “look for Bigfoot genes.” It is:

debug mtDNA → recover the missing autosomal relatedness output → understand Zana’s residual Y reads → identify the paper’s kinship estimator → select proper ancient-human controls.

Only after those are clean do we start the first genuinely original question:

Does Zana’s nuclear genome contain anything that falls meaningfully outside ordinary modern-human variation in a way the original ancestry-focused paper would not have noticed?

That’s where we are: we’ve built and validated the microscope; we haven’t started claiming we see monsters through it.

Let me look at the actual outputs — the @SQ line in the mtDNA report should identify the zero-variant bug immediately.

Both bugs are confirmed, and I can see a third problem in my own X test. Let me verify the key statistical point before rewriting.

Confirmed — my X test as written can't discriminate. Rewriting with depth stratification, plus the two real bugs.

Both bugs are fixed, and the paper's kinship method is now pinned down — which resolves item 6 entirely. Dispatching Qwen 5 and Qwen 7 on the corrected state.

Both converge, and Qwen 7 made one sharp observation I can test quantitatively right now.

The asymmetry is clean and runs in opposite directions. Now checking a candidate control's pipeline before ranking.

ENA's API is down for every query, including ones that worked earlier — so the Viking reference check has to wait. Writing up the corrected state and the control ranking.

The assert fired before any write, so the file is untouched. Let me redo it with index-based edits.

Corrected Phase 2 report. Three bugs found, all fixed, and one of them was in my own test design. Bigfoot stays unpublished; nothing committed.

Damage — reproduced, and it's genuine

Both samples show the real double-stranded-library signature: 5′ C→T and the complementary 3′ G→A, both decaying over about ten bases, against a flat transversion floor of 0.05–0.07%. Zana runs 5.70% → 0.75% at the 5′ end; Khwit 1.89%.

In the length-matched 80–100 bp band Zana's terminal C→T is about 2.8× Khwit's, so the difference survives the read-length confound and is real. But I've rejected Qwen 5's reading that this means Zana is older or more thermally cycled — burial microenvironment, skeletal element (both samples pool petrous and tooth libraries) and library chemistry all move deamination, and they died decades apart. The measurement stands; the cause doesn't follow.

The three bugs

mtDNA zero variants — contig-name mismatch. The deposited rCRS BAMs name their contig gi|251831106|ref|NC_012920.1|; my FASTA was headed NC_012920.1. bcftools couldn't match it and failed outright — and my wrapper had sent stderr to /dev/null, so a hard error surfaced as "0 variants" against 628× coverage. That suppression was my mistake and the real lesson. Replaced with a pileup built directly from each read's CIGAR, with the terminal 5 bp genuinely masked. Your -q/-Q correction was right and the false comment is gone.

Autosomal numbers missing — dropped regions. samtools mpileup honours only one -r, so passing three regions analysed X alone. The counters stayed at zero and their report lines were never emitted. Not a silent numerical failure — a dropped-region bug. Now one pileup per region.

The X result was weaker than it looked — my design flaw. Before interpreting it I checked whether that rate could discriminate at all, and it can't. 1 violation in 5,535 sites is exactly where Zana's allelic dropout at heterozygous sites lands (X heterozygosity around 1/1000–1/2000 × dropout of 25–50% at 3.17× gives 1 in 2,000–4,000), and also where an unrelated pair would land. The denominator is dominated by sites invariant between any two humans. So please drop "very strong result" — it demonstrated nothing.

The fix is self-calibrating: stratify both statistics by depth. Dropout-driven mismatches must fall as depth rises; genuine mismatches from an unrelated pair stay flat. The trend is the evidence, never the pooled number.

Kinship — closed

From the paper's methods: PC-relate in PCAngsd plus realSFS on the Affymetrix Human Origins panel. Quoting it: a coefficient of "0.1818 ... falling within the first degree relationship interval (0.177‐0.354)", and separately "R0 = 0.0015 ... suggesting a parent‐offspring relationship."

So 0.1818-against-0.25 was the wrong comparison, and my earlier flag was wrong — it sits inside that estimator's own interval, and the parent-offspring call rests on R0, not the coefficient. R0 counts sites where the pair shares no allele, which a parent and child cannot produce except through error. My opposite-homozygote statistic is the same family, which is a useful convergence. Contamination was contamMix plus ANGSD; haplogroup was haplogrep2.

Zana's Y reads — mis-mapping, and now quantified

Qwen 7 spotted that Zana's zero-mismatch fraction on Y (41.3%) is lower than Khwit's (67.0%). I tested that against each sample's own genome-wide expectation, which is the fair comparison since Zana's overall mismatch rate is 1.8× Khwit's:

Opposite directions. Reads forced onto a paralogous target carry excess mismatches; genuine Y sequence carries fewer. That favours cross-mapping over male contamination. Her sex is untouched — X:autosome 0.945 can't be male however badly Y maps.

Also: the deposit holds zero MAPQ-0 reads, so "100% MAPQ≥30" and "0 reads in the pseudoautosomal regions" are artefacts of pre-deposit filtering, not evidence about ancestry. ydist will localise the reads along Y.

What changed in the plan

Red Team's objection I'm accepting: depth correlates with mappability, so the depth trend could fall for any pair, related or not. A control is therefore now a prerequisite for the parent-offspring-versus-sibling claim, not an optional extra. Ranked candidates are in data/bigfoot-dna/CONTROLS.md, led by an unrelated ancient pair through the same pipeline (Viking, PRJEB37976) plus known-relationship pairs from the same project for calibration. Nothing downloaded — and ENA's API went down mid-session for every query, so each candidate's reference build still has to be read from its own BAM header first.

DISCARDED.md is now 64 rows, including ten corrections from this round.

To re-run, in order

/Users/Shared/AI/detective-desk/scripts/bigfoot_zana_phase2.py mtdna

Then relatedness, then ydist. Don't re-run damage, sexdetail or coverage — those worked.

Sanity check before trusting the mtDNA output: an African macrohaplogroup-L mitogenome should differ from rCRS at roughly 40–70 homoplasmic positions. Around 0 or 10 means still broken; around 200 means NUMTs or mis-mapping. Haplogrep stays uninstalled until real variants come back.

Current state is unchanged: nothing observed so far distinguishes H4 from H1, and Zana continues to behave as an ordinary ancient modern-human genome.

Yep — this is exactly where I’d go next. Claude has cleaned up the three bad assumptions, and the remaining work is now nicely scoped.

Run these one at a time, in this order:

cd /Users/Shared/AI/detective-desk

./scripts/bigfoot_zana_phase2.py mtdna

For that first one, the sanity check is important:

• roughly 40–70 mtDNA differences vs rCRS = plausible
• around 0–10 = still broken
• around 200 = likely mapping/NUMT trouble

If that looks sane, then:

./scripts/bigfoot_zana_phase2.py relatedness

Then:

./scripts/bigfoot_zana_phase2.py ydist

I would not rerun damage, sexdetail, or coverage; those are already good.

What changed conceptually is important too:

• The ancient-DNA damage result stands.
• Sex stands.
• The original X “99.982% match” should be discarded as evidence of close relationship; Claude correctly recognized that invariant human sites dominate the denominator.
• The relatedness question now needs a depth-stratified comparison plus controls, because dropout/mappability can mimic the trend.
• The Y reads now look more like cross-mapping than male contamination, because they perform worse than Zana’s own genome-wide mismatch expectation while Khwit’s true Y reads perform better.

And the kinship coefficient issue is now basically resolved: 0.1818 is not anomalous under the estimator they actually used, and the parent-offspring interpretation came from the very low R0 = 0.0015, not from comparing the coefficient to a naive 0.25 expectation.

So the immediate goal is now:

fix mtDNA → recover autosomal relatedness stats → map the Y-read distribution → then bring in controls.

Once you paste the mtdna output, I can tell you immediately whether the fix worked or whether there’s still a pipeline problem.

rror + dropout.

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

So a near-zero OH rate argues against full siblings -- but only as far as the error floor

(the transversion OH rate) is below what IBD0 would contribute. State both numbers.

* Second-degree (avuncular, half-sib, grandparent) also produce IBD0 and so OH > 0.

* This CANNOT be converted into a kinship coefficient comparable to the paper's 0.1818 without

population allele frequencies, which we have not downloaded. No coefficient is reported here.

* A son's single X comes from his mother, so a mother-son pair should show almost no violations.

Failure modes: a violation can also be error, damage, dropout in Zana, or mis-mapping; and the

test does not distinguish mother from maternal grandmother/aunt on its own.

* Both samples pool petrous and tooth libraries; neither is a single-library observation.

* No unrelated same-population pair processed through this pipeline is available as a null.

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

cliffchristianson@Cliffs-Mac-Studio detective-desk % cd /Users/Shared/AI/detective-desk

./scripts/bigfoot_zana_phase2.py mtdna

mtDNA vs rCRS (NC_012920.1, 16569 bp). Pileup built here from CIGAR; the deposited

rCRS BAMs name their contig 'gi|...|ref|NC_012920.1|', which is why a bcftools run

against a FASTA headed 'NC_012920.1' silently produced zero variants (contig mismatch).

Calling: base quality >= 30, terminal 5 bp of each read masked,

depth(masked) >= 20, alt fraction >= 0.9. Both masked and unmasked

call sets are reported so the effect of terminal damage is visible.

== Zana    contig ['gi|251831106|ref|NC_012920.1|'] | @PG 28 records    alignment dirs in @PG: ['Bigfoot_rCRS_aln']    bwa reference in @PG: ['rCRS.fasta']    reads used 153,441 | positions with masked depth >= 20: 16,518    VARIANTS vs rCRS (masked, reliable regions): 64  [61 transitions, 3 transversions]    73A>G(94%,303x) 150C>T(96%,26x) 263A>G(96%,393x) 383T>C(98%,252x) 385A>G(92%,243x) 418C>T(99%,353x) 750A>G(94%,512x) 769G>A(99%,452x) 1018G>A(99%,601x) 1438A>G(91%,350x) 1442G>A(97%,353x) 1706C>T(100%,542x) 2332C>T(99%,540x) 2358A>G(93%,595x) 2416T>C(95%,545x) 2706A>G(92%,513x) 3594C>T(100%,328x) 4104A>G(96%,477x) 4158A>G(91%,466x) 4370T>C(95%,527x) 4767A>G(91%,325x) 5027C>T(99%,409x) 5331C>A(100%,438x) 5814T>C(93%,512x) 6026G>A(99%,388x) 6629A>G(93%,489x) 6713C>T(100%,470x) 7028C>T(99%,425x) 7256C>T(100%,492x) 7521G>A(100%,402x) 7624T>A(100%,399x) 8080C>T(100%,532x) 8206G>A(99%,461x) 8387G>A(99%,493x) 8701A>G(90%,569x) 8860A>G(98%,368x) 9221A>G(95%,475x) 9540T>C(95%,315x) 10115T>C(95%,452x) 10398A>G(95%,629x) 10828T>C(95%,632x) 10873T>C(92%,422x) 11002A>G(93%,411x) 11719G>A(100%,493x) 11944T>C(97%,576x) 12236G>A(100%,289x) 12705C>T(99%,504x) 12948A>G(94%,548x) 13590G>A(99%,472x) 13650C>T(100%,438x) 13924C>T(99%,619x) 14059A>G(93%,645x) 14766C>T(100%,525x) 15110G>A(98%,456x) 15217G>A(100%,410x) 15301G>A(100%,409x) 15326A>G(94%,399x) 16114C>A(97%,307x) 16129G>A(95%,250x) 16213G>A(99%,271x) 16223C>T(99%,344x) 16278C>T(99%,341x) 16362T>C(96%,282x) 16390G>A(100%,252x)    excluded as unreliable regions: 0     calls present WITHOUT end-masking but absent with it (damage-suspect): 146 198    calls present WITH masking but absent without: 1438 4767 8701    positions with a 10-50% minor allele (heteroplasmy OR contamination OR NUMT): 7    146:12% 152:13% 204:31% 310:42% 4769:10% 16180:10% 16558:10%

== Khwit    contig ['gi|251831106|ref|NC_012920.1|'] | @PG 30 records    alignment dirs in @PG: ['Bigfoot_rCRS_aln']    bwa reference in @PG: ['rCRS.fasta']    reads used 121,220 | positions with masked depth >= 20: 16,517    VARIANTS vs rCRS (masked, reliable regions): 67  [64 transitions, 3 transversions]    73A>G(97%,291x) 146T>C(97%,38x) 150C>T(97%,29x) 152T>C(97%,31x) 263A>G(99%,348x) 383T>C(97%,252x) 385A>G(93%,246x) 418C>T(100%,292x) 750A>G(95%,426x) 769G>A(99%,399x) 1018G>A(100%,479x) 1438A>G(95%,322x) 1442G>A(98%,331x) 1706C>T(99%,397x) 2332C>T(99%,462x) 2358A>G(95%,506x) 2416T>C(96%,447x) 2706A>G(94%,392x) 3594C>T(99%,204x) 4104A>G(96%,383x) 4158A>G(92%,393x) 4370T>C(96%,450x) 4767A>G(94%,312x) 4769A>G(95%,306x) 5027C>T(99%,370x) 5331C>A(100%,415x) 5814T>C(94%,430x) 6026G>A(100%,375x) 6629A>G(96%,414x) 6713C>T(99%,389x) 7028C>T(100%,335x) 7256C>T(100%,414x) 7521G>A(99%,410x) 7624T>A(99%,366x) 8080C>T(99%,454x) 8206G>A(100%,413x) 8387G>A(99%,358x) 8701A>G(94%,409x) 8860A>G(98%,327x) 9221A>G(93%,384x) 9540T>C(97%,309x) 10115T>C(93%,355x) 10398A>G(96%,489x) 10828T>C(96%,457x) 10873T>C(95%,286x) 11002A>G(95%,382x) 11719G>A(99%,395x) 11944T>C(96%,428x) 12236G>A(100%,265x) 12705C>T(99%,414x) 12948A>G(93%,442x) 13590G>A(99%,397x) 13650C>T(99%,380x) 13924C>T(99%,487x) 14059A>G(95%,556x) 14766C>T(100%,462x) 15110G>A(99%,388x) 15217G>A(99%,324x) 15301G>A(100%,404x) 15326A>G(96%,419x) 16114C>A(99%,250x) 16129G>A(98%,227x) 16213G>A(99%,255x) 16223C>T(99%,257x) 16278C>T(100%,270x) 16362T>C(98%,260x) 16390G>A(99%,209x)    excluded as unreliable regions: 1 16180A>G    calls present WITHOUT end-masking but absent with it (damage-suspect): 198    calls present WITH masking but absent without: none    positions with a 10-50% minor allele (heteroplasmy OR contamination OR NUMT): 4    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 % 

relatedness is still running

Yep — the mtDNA bug is fixed. This output now looks biologically sane.

The big sanity check was “do we get something like 40–70 real differences from rCRS?” and we got:

• Zana: 64 variants — 61 transitions, 3 transversions
• Khwit: 67 variants — 64 transitions, 3 transversions

That is right in the expected range Claude gave you. So the previous 0 variants was absolutely just the contig-name/software error, not anything weird about the DNA.

Even better, the two mitochondrial genomes are strikingly similar, exactly what you'd expect from a mother and son. Most of the variant list is identical. The three obvious differences are interesting rather than alarming:

• 146 — Zana has ~12% minor allele; Khwit calls it ~97%
• 152 — Zana ~13% minor; Khwit ~97%
• 4769 — Zana ~10% minor; Khwit ~95%

That is actually a really nice thing to investigate. It could reflect heteroplasmy—a mother can carry a mixture of mitochondrial variants and transmit a different proportion to a child because of the mitochondrial bottleneck. But at these particular sites we also have ancient-DNA damage, NUMTs, contamination, and alignment artifacts to rule out. So don't call it maternal heteroplasmy yet.

The important part is that this is no longer:

“Why are there zero mtDNA variants?”

It's now:

“We have two very similar, high-depth mitochondrial genomes with a handful of differences worth understanding.”

And the masking is visibly doing useful work. For example, 198 disappears once read ends are masked in both samples, which is exactly the kind of thing we'd expect if terminal damage were creating a false call.

There's also a nice software lesson here: stderr suppression converted a hard bcftools failure into a biologically absurd result. That should absolutely go into DISCARDED.md/the engineering notes as a rule: genomics commands that produce scientific conclusions should never have their errors silently thrown away.

I would let relatedness finish before sending Claude another big handoff. Then we'll give him one consolidated result containing corrected mtDNA + relatedness + ydist.

After that, I'd approve haplogroup assignment—but I'd do it in two layers:

1. Reproduce the paper as closely as practical using its HaploGrep/PhyloTree methodology.
2. Then run the current HaploGrep version as a modern cross-check.

If both independently put Zana and Khwit in L2b1b, that's another piece of the original study we've reproduced ourselves.

So right now: don't interrupt anything. The mtDNA step passed. Let relatedness cook.

how long do you think this relatedness is going to run?

Probably on the order of tens of minutes, not hours, assuming it’s doing separate samtools mpileup passes over chr20, chr21, and X for both BAMs.

Given what we’ve already seen, my rough guess would be:

• best case: 5–15 minutes
• more likely: 15–40 minutes
• if it’s doing repeated full scans / Python post-processing inefficiently: maybe up to an hour

The BAMs are ~10 GB each, and a lot of these samtools steps are effectively single-core sequential scans.

You can see exactly what it’s doing right now with:

ps aux | grep '[b]igfoot_zana_phase2\|[s]amtools'

If you paste that output, I can tell you whether it’s on chr20/21/X and roughly how far along it probably is.