Review of sketchy v0.1, and what v0.2 changes

This is the review that motivated the rewrite. It was written after reading every source file, every benchmark script and every number in the v0.1 README, and after re-running the v0.1 tool and its scoring scripts to confirm that the published numbers reproduce (they do: 94.1% of HG002 fragments within 100 kb, 60.8% for alpha-satellite HOR, 42.1% for HSat2).

In summary, the empirical work in v0.1 was careful and largely correct about its own limits; the method it evaluated was poorly matched to the problem, and the benchmark needed to decide the question it posed was absent. Version 0.2 addresses both.

1. What v0.1 got right

• The density-equalizing sketch. Keeping a window's whole vocabulary when it is small (satellite arrays: ~2,300 distinct 21-mers per 10 kb) and a fixed budget when it is large (unique sequence: ~9,950) is the right allocation for localization, and v0.1 measured it fairly against bottom-k and FracMinHash at matched index bytes. This idea survives unchanged into v0.2.
• Transparent reporting. Four headline claims were retracted in the git history after adversarial review; confidence columns were audited and one (p_value) was found to be a constant and removed; the fragment benchmark's cluster structure (design effect 5.9) was measured and reported. Version 0.2 retains these practices.
• The truth-free centromere walk (cut a contig into consecutive pieces, place each independently, score the order) is a genuinely useful design for comparing localizers where no truth set exists, and its metric self-test is sound. It is kept as-is.
• The measurement that locus-private k-mers transfer between haplotypes only 2-8% of the time in satellite arrays (probes 8 and 9) is the most important scientific fact in the repository. It says that any method anchoring on rare k-mers -- VerityMap, TandemMapper, centroFlye, and v0.1 itself -- is bounded a priori when the query comes from a different haplotype than the reference. v0.2 is designed around that fact rather than against it.

2. What was wrong with the method

2.1 The score was a heuristic with no probabilistic reading

v0.1 scored a tile against a window by a tf-idf cosine with the rarity weight raised to the power 2.5 in the numerator but power 1 in the reference norm. The exponent was tuned on the benchmark, and the norm asymmetry was chosen because the symmetric form scored worse. The resulting score was, as the README stated, "not bounded by 1, not comparable between queries". Every confidence statistic (snr, margin, tile_support, support_span, rival_frac) was then a post-hoc function of that score, and the audit found that support_span covered the truth only 26-65% of the time inside satellites.

A cosine is a similarity, not a likelihood. Nothing in it says how probable it is that the query came from this window rather than that one, so nothing downstream could be calibrated except by fitting a separate model on top (the split-conformal interval in stat_audit.py, which was never wired into the tool).

Version 0.2 replaces the score with identity evidence that has a direct probabilistic reading: each non-overlapping segment of shared k-mers contributes ln(N/df) for its rarest word, the log-odds of observing that word in the true window rather than in one drawn at random. (A term penalising unexplained words, the textbook likelihood, was implemented and rejected; docs/METHODS.md §3 gives the measurement.) The profile over windows becomes a posterior through a softmax whose temperature is fitted on reads whose origin is known exactly. The posterior mass of the call is the confidence and its Phred form is mapq; the README's Section 4.4 reports its calibration on held-out error models and read lengths.

2.2 The query was sketched, and that created a problem that then needed three fixes

The reference must be sketched (that is what makes the index 2 GB rather than 30 GB). The query does not need to be: it is one sequence, its k-mers cost nothing to enumerate, and looking up ten thousand hashes takes a millisecond. v0.1 nevertheless sketched the query by a rule of its own -- first the "ramp", then a per-query blend between ramp and step, chosen from seen_frac -- and then discovered that sketching the query more densely than the reference biases the cosine (the norms are over different hash spaces), tried a per-window norm correction, found pruning worked better, and shipped the blend. Roughly a third of locate.rs and of the README's method section dealt with this.

Version 0.2 does not sketch the query. Every query k-mer is looked up; a window can share only the words whose hash lies at or below its own threshold, a prefix of the hash space that is exactly what the reference stored, so no correction is required. The ramp, the step, the blend, --query-blend, --per-window-norm and seen_frac-driven pruning are removed.

2.3 The link channel did not earn its half of the index

The order-aware "link" channel (hashed pairs of consecutive minimizers) was the only positional information in v0.1. The matched ablation put its contribution at +0.009 with p = 0.29, and it accounted for 94 M of the 222 M postings. It was an indirect way to recover order from a bag-of-words model.

Version 0.2 stores every position of every retained k-mer (u32 with a strand bit) and runs a positional stage inside the candidate regions of the coarse posterior: every query k-mer votes for the diagonal ref - query of every stored occurrence, weighted by ln(N/df). This is the FASTA diagonal method (1985) with rarity weights, which v0.1's prior-art section itself named, and it gives base-level resolution, the strand (including for single-window queries, which v0.1 could not orient) and a second calibrated posterior.

2.4 Resolution was bounded by the window

v0.1's median error in unique sequence was 2.2-2.5 kb, a consequence of snapping to a 10 kb grid and interpolating rather than a property of the sequence. The fine stage removes that floor; on clean sequence v0.2 places to the base.

3. What was wrong with the evaluation

3.1 The benchmark could not answer the question posed

Every real-data number in v0.1 scores HG002 fragments against CHM13, with truth from minimap2. Inside a satellite array that conflates three quantities: the method's error, the aligner's error (78% of alpha-satellite truth calls carried mapq ≤ 5), and the genuine difference between two haplotypes' arrays, which Logsdon et al. (2024) quantify as 45.8% of centromeric sequence being unalignable between haplotypes. v0.1 acknowledged this, reported satellite accuracy as a range rather than a number, and ran one small within-sample test (59 satellite reads on the HG002 HPRC draft assembly, whose centromeres are incomplete).

Version 0.2 uses the T2T HG002 v1.1 diploid assembly, in which both haplotypes are complete through every centromere and which carries its own cenSat annotation, and simulates reads from it with stated error models. The origin of every base is then known exactly, inside the centromere and out, for both haplotypes. This benchmark can determine whether composition places a read within its array, at what resolution, and whether an aligner does. The cross-haplotype benchmarks (HG002 and HG02723 fragments against CHM13) are retained as the second question: what remains recoverable when the reference is not the sample's haplotype.

3.2 The failure of alignment was inferred from one contig walk

The claim that alignment fails in satellites rested in v0.1 on the walk, where sketchy's order statistic exceeded every aligner's inside the HOR array, and on the mapq of the truth calls. Neither measures an aligner's accuracy in a satellite, because neither had a truth. With exact-truth simulated reads, v0.2 measures minimap2, Winnowmap2 and MashMap3 on the same reads, by region class, with wall time and peak memory recorded, and reports the aligners' advantage where it exists.

3.3 Simulation was dismissed too early

v0.1 concluded that simulated divergence is a structurally different and easier problem after finding that tiling helped on simulated data and not on real cross-haplotype data. That holds for cross-haplotype placement, where the difference is the haplotype rather than the error model. It is not a reason to avoid simulation for the same-haplotype question, where simulation is the only route to exact truth inside a satellite. v0.2 uses both and states which question each answers.

4. Compute

An intermediate v0.2 index with one position per word was 2.05 GB for CHM13 and built in 34 s; storing every position (+50%) is what bought within-array resolution, and the repeat-aware second pass adds another 17% for a five-point gain in accuracy. Query cost is measured in the write-up alongside the aligners on identical inputs.

5. What v0.2 does not change

The density-equalizing keep rule; the 10 kb window; the inverted index with a 22-bit bucket table over sorted keys; the block-parallel BGZF reader; the benchmark harness, samplers and the truth-free sweep. The README's prior-art section was accurate and is carried forward.

6. Claims from v0.1 that v0.2 retires

• "margin > 3 is the column to filter on" -- replaced by a calibrated posterior; filter on mapq or post like any other mapper.
• "snr, tile_support, support_span, rival_frac" -- gone, along with the audit that showed which of them meant what. cred_span / cred_loci (the 90% credible set) carry the shape information with a coverage guarantee the old columns lacked.
• "the query keep rule is chosen per query" -- there is no query keep rule.
• "orientation needs more than one tile" -- the fine stage orients single-tile queries.

Review and rewrite by Claude Fable 5.1 (Anthropic) for @jlanej, September 2026.

Methods (v0.2)

sketchy places a query sequence on a reference genome using its k-mer composition, and reports a posterior probability for the call. Alignment appears in this repository only offline, to establish ground truth for the cross-haplotype benchmarks, and as a competitor in the benchmarks themselves. It is never available to the method.

The version-2 method has four parts: a density-equalizing reference sketch (unchanged from v0.1), a coarse stage that scores windows by composition, a fine stage that resolves position by rare-word diagonals inside candidate regions, and a calibrated posterior that combines the two. Each is stated here in full.

1. Features

Every canonical k-mer (default k = 21) is hashed with a splitmix64 finalizer. Canonical means the lexicographically smaller of the k-mer and its reverse complement, so a sequence and its reverse complement produce identical hash multisets. Alongside each hash the index keeps whether the forward k-mer was the canonical one (one bit), which is what lets the fine stage recover strand by comparing that bit between query and reference. A non-ACGT base resets the rolling window.

2. The reference index: a density-equalizing sketch, with positions

The reference is cut into windows of W bp (default 10 kb). A window's distinct k-mers are ranked by hash and the smallest are kept -- a bottom-k sketch -- with the keep count set by the step rule:

  keep(n) = n         if n <= V^2 / (2V - budget)     (the crossover)
          = budget    otherwise

V is the vocabulary of a fully unique window of this size, estimated from the reference at build time (98th percentile of distinct counts over 4,096 sampled windows; 9,980 for CHM13 at k = 21, W = 10 kb) and budget is 256. A 10 kb alpha-satellite HOR window holds about 2,300 distinct 21-mers, an HSat3 window about 4,000, a unique window about 9,950: the windows that need every word they have are exactly those with few enough words to afford it. Windows with less than 50% A/C/G/T are not indexed.

Because the kept set is a hash prefix (every hash at or below the window's threshold), two windows sketched at different densities are comparable on the prefix they share, and an unsketched query is comparable to every window on that window's own prefix. That property is what removes any need to sketch the query.

For every kept k-mer the index stores every position at which it occurs in the window (chromosome-local, with the strand bit). v0.1 stored one position per word per window; inside a tandem repeat almost every word occurs several times per window, and a word with one stored position can vote for only one of its copies. The fine stage's within-array resolution depends on this.

On disk the index is one file of little-endian arrays, 64-byte aligned, opened by mmap: a bucket table over the top 22 bits of each hash, the sorted distinct keys, a CSR offset array, and per posting the window id (u32), multiplicity (u16) and an offset into the position array. A key lookup costs about two cache misses.

3. Coarse stage: composition over windows

The query is cut into tiles of W bp (non-overlapping by default). Every k-mer of a tile is looked up; the query is never sketched.

Evidence. For a tile and a reference window, the identity evidence is

   E(w) = sum over non-overlapping shared segments of  ln( N / df(h*) )

where df(h) is the number of windows carrying h, N the number of windows, and h* the rarest word of the segment. A segment is opened by the first shared k-mer at a query position and absorbs every shared k-mer that starts within k bases of it; the next segment opens at the next shared k-mer beyond. ln(N/df) is the log-odds that a word would be seen in a window that carries it against a window drawn at random.

Two choices here are the result of measurements, not taste:

• Segments, not k-mers. Overlapping k-mers are not independent observations. One Alu in a query shares ~300 consecutive k-mers with every window that carries an Alu, and counting each of them let a repeat-dense window that stores its whole vocabulary outscore the true unique window, whose sketch can only ever share a hash-subsample (unique-sequence accuracy 92% against 99% once segments were used). Crediting each segment with its rarest word rather than its first is what lets a private variant k-mer separate two windows of one array that share every common word.
• No penalty for unexplained words. A likelihood model that charges a window ln(1 - p) for every comparable word it fails to share is the textbook form, and it was tried. It fails on exactly the queries this tool exists for: a read from another haplotype of a satellite array carries many words that exist elsewhere in the array, and charging its true window for them let windows that share nothing win (alpha-satellite HOR right chromosome 73% against 97% without the term). Words the window does not share are therefore uninformative, and a window that shares nothing is excluded from the posterior rather than merely penalised.

Words carried by more than max_df_frac of windows (default 1%) are skipped; they carry less than 4.6 nats of information each and dominate running time.

Shift-and-sum. A tile at query offset a matching window j is evidence that the query starts at window j - round(a/W); reverse-complemented, at j - round((qlen - b)/W). Each tile's evidence is added into forward and reverse profiles over all windows. The profiles are sparse: only windows that share a segment with some tile are non-zero.

Phase correction. A query whose tiles are out of phase with the reference windows splits each tile's evidence between two consecutive starts, so its peak is about half of an in-phase query's. Two haplotypes of one genome are never in phase with each other, and on the raw profile the choice between them was decided by phase rather than evidence (haplotype accuracy 53% in unique sequence). The decision profile is therefore L(s) + L(s+1) on each strand.

Posterior. With temperature T,

   P(s) = exp(L(s)/T) / sum over starts and strands with L > 0 of exp(L/T)

The MAP start, its posterior mass within one window (coarse_post), the cred-level credible set (its bp span and number of separate runs: cred_span, cred_loci), and the evidence gap to the best start outside the guard region around the MAP (gap, nats) are reported. T is calibrated (section 6); it is not 1, because segments are still far from independent.

Candidate regions. Every maximal run of starts whose decision-profile value is at least 40% of the peak, on either strand, is a candidate region, padded by the query's length in windows plus two so the query can sit anywhere inside. The MAP's region comes first; at most six are kept. The rule is independent of T on purpose.

3a. k-mer survival

For up to eight tiles spaced along the query, the window with the most evidence is found and the fraction of the tile's k-mer occurrences it explains, divided by the fraction of the hash space that window retained, is recorded; p_hat is the median over tiles. For a read from its own reference this is (1 - e)^k: 0.69 on the 1.5%-error set (expected 0.73), 0.21 on the 7% set (expected 0.22). It sets the fine stage's smoothing width and is reported.

4. Fine stage: rare-word diagonals inside candidate regions

Inside a candidate region, every query k-mer q at query offset x that is present in the index (and not above the df cutoff) is looked up, and for every stored occurrence at chromosome position r in a window of the region it casts a vote on a diagonal:

   same strand bit:      d = r - x          (forward)
   opposite strand bit:  d = r + x          (reverse; the query start is d - qlen + k)

with weight ln(N/df(q)), into histograms binned at 20 bp. The histograms are smoothed with a box whose half-width scales with the excursion of the diagonal expected under indels, 10 * sqrt(qlen * 0.65 * e) bp (at least one bin), where e is the query's error rate inferred from its k-mer survival (section 3a) and 0.65 is the indel share of nanopore errors. The peak gives the strand and, from the weight-averaged diagonal of the anchors under the peak, the query's start to the base. The width matters, and in opposite directions on simulated and real reads. A fixed ±200 bp costs two points at 1.5% error on simulated 20 kb reads (the variant k-mers that separate one tandem copy from the next are averaged away, in HSat3 and monomeric arrays especially) and a fixed ±20 bp costs six points at 7% error (the diagonal drifts out of the window); three standard deviations of the random walk is the optimum on simulated reads (79.8% within 10 kb). On real R10.4.1 ultra-long reads, whose indels are clustered (section 7), three sigma gives the lowest agreement with minimap2 of any setting tried (75.4%) and ten sigma the highest (78.4%, against 77.2% for the fixed 1% rule), for a cost of 0.6 points on simulated reads; ten is the default. This is the FASTA diagonal method (Lipman & Pearson 1985) with rarity weights and every occurrence of every retained word; it scores no gaps and compares no bases.

Two quantities are computed from each region's histogram: the per-anchor peak, used for the position and its posterior within the region, and a segment peak (the same non-overlap rule as the coarse stage, applied along the query per diagonal bin), used to choose between regions. Per-anchor sums are what separate one HOR copy from the next -- collapsing them into segments cost 17 points of accuracy in active HOR arrays -- but on their own scale they would overrule the coarse stage wherever anchors are dense, which is exactly where the coarse stage should have the final word on which array or haplotype (section 5).

5. Choosing the region: composition updated by position

Composition and position have complementary strengths, measured on exact-truth reads from the diploid HG002 assembly. Inside HOR arrays, few k-mers are single-copy and the fine stage is starved, and composition chooses the right array and haplotype far more often (active HOR 89.5% within 10 kb against 70%). In monomeric alpha satellite, beta satellite and unique sequence, the arrangement of the words matters more than the words, and position chooses better (beta satellite 80% against 33%). Neither alone is best everywhere.

The region is therefore chosen by a Bayesian update of the coarse posterior: region i with coarse mass m_i and segment fine peak F_i has posterior

   P(i) = m_i exp(F_i / T_f)  /  ( sum_j m_j exp(F_j / T_f) + (1 - sum_j m_j) )

where the last term is the coarse mass outside every candidate, which receives no fine update. The reported post is P(i*) times the within-region posterior mass of the fine peak (within fine_tol = 2 kb of it); mapq is -10 log10(1 - post), capped at 60. When no anchor is found in any region the coarse call stands and post is coarse_post.

6. Calibration

Both temperatures are fitted on reads whose origin is known exactly -- sketchy simulate cuts reads from an assembly, records their coordinates in the header, and applies a stated error model -- over a grid of (T, T_f) pairs (sketchy calibrate), with the decision temperature held fixed. The fit set is 3,800 20 kb ONT-R10-like reads from the HG002 v1.1 diploid assembly. Temperatures scale linearly with query length from that 20 kb reference, within a factor of eight: fitted separately, the best fine temperature is about four times smaller on 5 kb reads and four times larger on 100 kb reads, because a wrong call inside a near-identical repeat accumulates evidence along the whole read exactly as a right call does. Scaling by length halves the calibration error on 5 kb reads (0.30 to 0.14). Accuracy is flat (79.8%) for any coarse T at or above 64, so the coarse prior barely matters to the decision; the reported posterior is best balanced at T = 256, T_f = 128 (expected calibration error 0.111, negative log-likelihood 0.85 per read; the best ECE on the grid was 0.090 at 512/64 and the best NLL 0.68 at 512/256, and the shipped pair sits between them). The decision temperature is 8. The evaluation sets are different error models, read lengths and, for the real-read tests, different data altogether; the write-up reports expected calibration error on each.

7. Read simulation

Per base, with probability rate, an error occurs and is a substitution, insertion or deletion in the stated proportions; deletions inside homopolymers of length >= 4 are hp_bias times more likely; an insertion extends with probability ins_ext. Presets: asm (0.2%), hifi (0.3%, indel-dominated, homopolymer-biased), ont-r10 (1.5%), ont-r9 (7%). Reads are drawn uniformly, or per annotation class from a BED, from either strand. The model does not reproduce the context dependence of real nanopore errors; the real-read benchmarks exist for that reason.

One consequence was measured directly (bench/drift.py, from minimap2's CIGAR strings): the alignment diagonal of a real HG002 R10.4.1 ultra-long read wanders far more than the model's. Real reads have a lower indel rate than the ont-r10 preset (0.0039 against 0.0103 per base) but a maximum excursion of the diagonal from its starting offset of 96 bp (median) and 433 bp (90th percentile) per read of median 77 kb, 5 to 13 times the sqrt(L * r_indel) a random walk of independent indels predicts, against 1.4 to 1.9 times for simulated reads; within any 20 kb of a real read the 90th-percentile excursion is 180 bp, against 32 bp for simulated 20 kb reads. Indels in real reads are clustered. This is why the fine stage's smoothing width, which the simulations alone would set at the random-walk scale, is set wider (section 4), and it is a caution for any benchmark that tunes a mapper's gap tolerance on simulated nanopore reads.

8. Outputs

9. Command lines behind every table

Indexes (final settings; --repeat-min-df 2 for a haploid reference, 3 for a diploid one):

sketchy index data/chm13v2.0.fa.gz        -o data/chm13v2.0.skx      --repeat-min-df 2
sketchy index data/hg002v1.1.fasta.gz     -o data/hg002v1.1.skx
sketchy index data/hg002v1.1.pat.fa.gz    -o data/hg002v1.1.pat.skx  --repeat-min-df 2 --saturation-from data/hg002v1.1.skx

Exact-truth read sets from the diploid assembly, one per error model and length:

B=data/hg002v1.1.cenSatv2.0.bed; R=data/hg002v1.1.fasta.gz
sketchy simulate $R -o sim_L20k_ontr10.fa  -n 1000 -l 20000  --model ont-r10 --bed $B --per-class 200 --seed 11
sketchy simulate $R -o sim_L20k_ontr9.fa   -n 1000 -l 20000  --model ont-r9  --bed $B --per-class 200 --seed 12
sketchy simulate $R -o sim_L20k_hifi.fa    -n 1000 -l 20000  --model hifi    --bed $B --per-class 200 --seed 13
sketchy simulate $R -o sim_L100k_ontr10.fa -n 600  -l 100000 --model ont-r10 --bed $B --per-class 120 --seed 14
sketchy simulate $R -o sim_L5k_ontr10.fa   -n 1000 -l 5000   --model ont-r10 --bed $B --per-class 200 --seed 15

Placement and competitors (5 threads each; Winnowmap2 needs a meryl repetitive-k-mer list for the reference, WM_K15):

WM_K15=meryl/hg002v11_repetitive_k15.txt bench/run_methods.sh sim_L20k_ontr10.fa $R data/hg002v1.1.skx out/bench_L20k_ontr10 5 ont
bench/summarize_sim.py --queries sim_L20k_ontr10.fa --dir out/bench_L20k_ontr10 --tol 10000 --metric w10kb
bench/score_v2.py out/bench_L20k_ontr10/sketchy.tsv --tol 10000 --strip-hap        # per-class table with ECE
sketchy calibrate -x data/hg002v1.1.skx sim_L20k_ontr10.fa --tol 10000 --fine-tol 10000   # temperature grid
bench/fig_v2.py classbars|reliability|errcdf|cost ...                                       # figures

Cross-haplotype, maternal reads on the paternal index (the haplotype FASTAs are split from the diploid assembly; data/hg002v1.1.pat.skx is built from the paternal one with --repeat-min-df 2 as above):

gzip -dc $R | awk '/^>/{p=($0 ~ /_PATERNAL/)} p' | bgzip > data/hg002v1.1.pat.fa.gz
gzip -dc $R | awk '/^>/{p=($0 ~ /_MATERNAL/)} p' | bgzip > data/hg002v1.1.mat.fa.gz
grep _MATERNAL $B > censat_mat.bed
sketchy simulate data/hg002v1.1.mat.fa.gz -o mat_L20k_ontr10.fa -n 800 -l 20000 --model ont-r10 --bed censat_mat.bed --per-class 150 --seed 21
sketchy locate -x data/hg002v1.1.pat.skx mat_L20k_ontr10.fa -o sketchy_pat.tsv
minimap2 -cx map-ont --secondary=no data/hg002v1.1.pat.fa.gz mat_L20k_ontr10.fa > minimap2_pat.paf
bench/score_xhap.py --queries mat_L20k_ontr10.fa --bed $B --sketchy sketchy_pat.tsv --paf minimap2_pat.paf

Assembly segments on the acrocentric short arms (§4.11): the short arm of each acrocentric chromosome of both HG002 haplotypes, telomere to the start of its active HOR array, is cut into non-overlapping 100 kb tiles (tiles with fewer than 90% called bases, the unresolved rDNA gaps, are dropped) and placed with every method on the diploid assembly, on the paternal haplotype (maternal tiles) and on CHM13, with the aligners' asm20 presets. The reference standard across genomes is the position lifted through a chromosome-scale alignment (blocks >= 20 kb at mapq >= 30); the identity margin of a tile is the difference, in matched fraction, between the best and the second-best CHM13 short arm among every alignment minimap2 reports for it.

D=out/v4/acro; B=data/hg002v1.1.cenSatv2.0.bed
bench/acro_tiles.py --fasta data/hg002v1.1.mat.fa.gz --bed $B --len 100000 --prefix acromat_ --out $D/tiles_mat_L100k.fa
bench/acro_tiles.py --fasta data/hg002v1.1.pat.fa.gz --bed $B --len 100000 --prefix acropat_ --out $D/tiles_pat_L100k.fa
cat $D/tiles_mat_L100k.fa $D/tiles_pat_L100k.fa > $D/tiles_L100k.fa
samtools faidx data/hg002v1.1.mat.fa.gz chr13_MATERNAL chr14_MATERNAL chr15_MATERNAL chr21_MATERNAL chr22_MATERNAL > $D/hg002_acro_chroms.fa
samtools faidx data/hg002v1.1.pat.fa.gz chr13_PATERNAL chr14_PATERNAL chr15_PATERNAL chr21_PATERNAL chr22_PATERNAL >> $D/hg002_acro_chroms.fa
samtools faidx data/chm13v2.0.fa.gz chr13:1-15547593 chr14:1-10092112 chr15:1-16678794 chr21:1-10962853 chr22:1-12788180 | sed -E 's/^>(chr[0-9]+):.*/>\1/' > $D/chm13_parms.fa
minimap2 -t 6 -cx asm20 --secondary=no data/chm13v2.0.fa.gz $D/hg002_acro_chroms.fa > $D/hg002acro_vs_chm13.paf   # reference standard
minimap2 -t 2 -cx asm20 -N 20 -p 0.5 --secondary=yes $D/chm13_parms.fa $D/tiles_L100k.fa > $D/tiles_vs_chm13parms.paf   # margins
samtools faidx data/hg002v1.1.pat.fa.gz chr13_PATERNAL chr14_PATERNAL chr15_PATERNAL chr21_PATERNAL chr22_PATERNAL > $D/hg002_acro_pat.fa; samtools faidx data/hg002v1.1.mat.fa.gz chr13_MATERNAL chr14_MATERNAL chr15_MATERNAL chr21_MATERNAL chr22_MATERNAL > $D/hg002_acro_mat.fa
minimap2 -t 8 -cx asm20 -N 20 -p 0.5 --secondary=yes $D/hg002_acro_pat.fa $D/tiles_mat_L100k.fa > $D/mat_vs_patacro.paf   # nearest copy on the other haplotype
minimap2 -t 8 -cx asm20 -N 20 -p 0.5 --secondary=yes $D/hg002_acro_mat.fa $D/tiles_pat_L100k.fa > $D/pat_vs_matacro.paf
cat $D/mat_vs_patacro.paf $D/pat_vs_matacro.paf > $D/tiles_vs_otherhap.paf
minimap2 -t 4 -cx asm20 --secondary=no data/hg002v1.1.pat.fa.gz $D/hg002_acro_mat.fa > $D/matacro_vs_pat.paf   # reference standard across haplotypes
for ref in data/hg002v1.1.pat.fa.gz data/chm13v2.0.fa.gz; do   # Winnowmap2 k=19 lists, as for the diploid assembly
  tools/winnowmap/meryl count k=19 output mdb $ref memory=8 threads=6 && tools/winnowmap/meryl print greater-than distinct=0.9998 mdb > $D/$(basename $ref .fa.gz)_repetitive_k19.txt; rm -rf mdb
done
WM_K19=out/v2/meryl/hg002v11_repetitive_k19.txt bench/run_methods.sh $D/tiles_L100k.fa     data/hg002v1.1.fasta.gz data/hg002v1.1.skx     $D/same 8 asm
WM_K19=$D/pat_repetitive_k19.txt                bench/run_methods.sh $D/tiles_mat_L100k.fa data/hg002v1.1.pat.fa.gz data/hg002v1.1.pat.skx $D/xhap 8 asm
WM_K19=$D/chm13_repetitive_k19.txt              bench/run_methods.sh $D/tiles_L100k.fa     data/chm13v2.0.fa.gz    data/chm13v2.0.skx    $D/xind 8 asm
P="--placed sketchy=SET/sketchy.tsv --placed minimap2=SET/minimap2_asm20.paf --placed winnowmap2=SET/winnowmap2_asm20.paf --placed mashmap3=SET/mashmap3.paf"
bench/score_acro.py --setting same --queries $D/tiles_L100k.fa     --target-bed $B --liftover out/v2/xhap/mat_vs_pat.paf --margin-same $D/tiles_vs_otherhap.paf ${P//SET/$D/same} --tiles-out $D/same_tiles.tsv --md $D/same.md
bench/score_acro.py --setting xhap --queries $D/tiles_mat_L100k.fa --target-bed $B --liftover $D/matacro_vs_pat.paf ${P//SET/$D/xhap} --tiles-out $D/xhap_tiles.tsv --md $D/xhap.md
bench/score_acro.py --setting xind --queries $D/tiles_L100k.fa --target-bed data/chm13v2.0_censat_v2.1.bed --liftover $D/hg002acro_vs_chm13.paf --margin $D/tiles_vs_chm13parms.paf ${P//SET/$D/xind} --tiles-out $D/xind_tiles.tsv --md $D/xind.md
bench/fig_acro.py --tiles $D/xind_tiles.tsv --out docs/v2/fig_acro_xind.png; bench/fig_acro.py --tiles $D/xhap_tiles.tsv --out docs/v2/fig_acro_xhap.png

v0.1's HPRC benchmarks (truth from minimap2 -cx asm20 -N 30 -p 0.5 against CHM13, run once and kept in out/*_truth.tsv) are re-scored with bench/score_v2.py <placed.tsv> --truth <truth.tsv> --tol 100000; the truth-free walk is bench/make_walk.py + bench/sweep.py run/score exactly as in v0.1, with this version's locate in bench/sweep_methods.json.

Written by Claude Fable 5.1 (Anthropic) with @jlanej, September 2026; v0.1 methods by Claude Opus 5 and @jlanej.