Mask the blind spots a lift between distant anchors used to drop

CHANGELOG.md

@@ -84,6 +84,16 @@ based on [Keep a Changelog](https://keepachangelog.com/), and the project follow
 
 ### Fixed
 
+- **The blind-spot mask missed a third of the blind spots on any other reference.** Lifted onto
+  the run's reference by interpolating between shared unique k-mers, a blind spot was dropped
+  wherever the anchors around it were far apart, and blind spots are repeats, where anchors are
+  scarce. On two MTBC assemblies of a 185-sample cohort, the mask covered 63% and 65% of the blind
+  spots minimap2 places on each, leaving 46,700 to 49,800 bp of problematic sites to the variant
+  caller. The lift now aligns the stretch between two anchors where it cannot interpolate
+  (`--align-gaps`): 97% and 98% are masked at minimap2's place, with 700 to 1,100 bp masked
+  elsewhere. About 15% of the H37Rv blind spots are stretches those assemblies do not have at all,
+  with nothing to mask.
+
 - **The genome page drew one reference's genes under every sample.** The report took the GFF and the
   masked regions of the first reference of the run, and its length, for the whole landscape, so on a
   cohort mapped against two assemblies the second one's samples were drawn over the first one's genes,

bin/pathotypr_liftover.README.md

@@ -132,8 +132,10 @@ strand` for `lift --global-chain`. `lifted.bed` = the positions collapsed into i
   `pos_h37rv`, looked up by contig and position, and `bin/lift_vcf.py`, which moves each SNP to its H37Rv
   position so the canonical snpEff annotation reads H37Rv's gene and codon.
 - **(B) Blind-spots mask on any reference.** Lift `assets/H37Rv_blindspots.bed` **from H37Rv to the run's
-  reference**: `A = H37Rv`, `B = mapping reference`. The resulting BED (target coords) goes into the
-  per-reference exclusion — so the mask works even when the reference is not H37Rv.
+  reference**: `A = H37Rv`, `B = mapping reference`, with `--sample 10 --align-gaps` (PREPARE_REFERENCE). The
+  resulting BED (target coords, every contig) goes into the per-reference exclusion — so the mask works even
+  when the reference is not H37Rv. The blind spots are repeats, where anchors are scarce: interpolation alone
+  masked 63–65% of those minimap2 places on two MTBC assemblies, `--align-gaps` 97–98%, ~1 GB and ~15 s.
 
 ## Calibration and limits (be honest)
 

docs/pathotypr.md

@@ -130,7 +130,12 @@ lifts it onto the run reference, every contig of it, with the **same k-mer lifto
 `--canonical_ref` as the FASTA the blind-spots are defined on (default: the
 H37Rv-colinear MTBC ancestor bundled at `/opt/pathotypr/reference.fasta`; point it at a
 true H37Rv FASTA for an exact lift), so the mask is correct even when the reference is
-**not** H37Rv. Set `--blindspot_liftover false` — which appends the BED directly and
-ignores `--canonical_ref` — only when the reference already shares H37Rv coordinates.
+**not** H37Rv. The blind spots sit in repeats, where shared unique k-mers are scarce, so
+the stretch between two anchors is aligned rather than interpolated (`--align-gaps`): on
+two MTBC assemblies of a 185-sample cohort, 97–98% of the blind spots minimap2 places on
+each are masked at the same place, where interpolation alone masked 63–65%; about 15% of
+the H37Rv blind spots are stretches those assemblies do not have, with nothing to mask. Set
+`--blindspot_liftover false` — which appends the BED directly and ignores
+`--canonical_ref` — only when the reference already shares H37Rv coordinates.
 The repetitive fraction is already covered by the pipeline's own repeat /
 mappability masking — the blind-spots add the non-repetitive problematic sites.

modules/reference.nf

@@ -107,11 +107,13 @@ process PREPARE_REFERENCE {
         CONTIG=$(head -1 reference.fa | sed 's/^>//; s/[[:space:]].*//')
         if [[ "!{params.blindspot_liftover}" == "true" && -s "!{params.canonical_ref}" ]]; then
             # Whole-genome anchor-chain liftover of the H37Rv blind-spots onto THIS reference, so the mask is
-            # correct without assuming shared coordinates: every position is placed by interpolation between
-            # flanking unique anchors (handles SNP sites + indels; anchor-desert positions drop, never mis-map),
-            # on whichever contig of the reference it lies, and the BED names that contig.
+            # correct without assuming shared coordinates: every position is placed between flanking unique
+            # anchors, on whichever contig of the reference it lies, and the BED names that contig. The blind
+            # spots sit in repeats, where anchors are scarce: interpolating only across short colinear gaps
+            # left a third of them unmasked on another MTBC assembly, so the stretch between two anchors is
+            # aligned instead (--align-gaps), which masks 97-98% of them where minimap2 puts them.
             python3 !{projectDir}/bin/pathotypr_liftover.py lift "!{params.blindspot_bed}" "!{params.canonical_ref}" reference.fa \
-                --out-bed bs_lifted.bed --kmer-size 21 --global-chain --sample 10
+                --out-bed bs_lifted.bed --kmer-size 21 --global-chain --sample 10 --align-gaps
             awk 'BEGIN{OFS="\t"} $2 ~ /^[0-9]+$/ && $3 ~ /^[0-9]+$/ {print $1, $2+1, $3, "blindspot", ""}' bs_lifted.bed >> "$out"
         else
             # reference already shares H37Rv coordinates: append the blind-spots directly (contig rewrite + 1-based)