cactus-graphmap-join: a `--reference` sample name that isn't in the graph fails with an opaque `VPKG::load_one` error

[diekhans]

cactus-graphmap-join: a --reference sample name that isn't in the graph fails with an opaque VPKG::load_one error

Diagnosis and bug report by Dr. Claude Shannon

Related: see #1970 (cactus-graphmap-join v3.2.1: multi-reference VCF is broken) for a separate, more serious
issue in the deconstruct job hit by the same workflow.

Versions

• Cactus v3.2.1 (quay.io/comparative-genomics-toolkit/cactus:v3.2.1, --binariesMode local)
• Toil 9.4.1-c6a39f2c827899f9a1dd0018ca2b72c675120096
• Source line numbers below are from cactus-bin-v3.0.1, where the relevant code is
  unchanged; the traceback line numbers are from v3.2.1.

Summary

If a sample name passed to --reference / --vcfReference does not exist in the input
--vg graphs, cactus-graphmap-join does not report that. Instead the clip_vg job
fails several tools downstream with:

error[VPKG::load_one]: Correct input type not found in standard input while loading handlegraph::MutablePathMutableHandleGraph

Nothing in the message names the reference, the sample names that are present, or even
hints that the reference is the problem. The only clue is a non-fatal warning[vg paths]
line earlier in the same stderr blob.

This bites in two distinct situations, and both produce the identical unhelpful message:

1. The reference name is simply wrong (typo, or a name that never existed).
2. The reference name is correct, but that sample has no path in this chromosome's
   graph — e.g. a female sample when joining chrY, or any sample whose mitochondrion
   is not in the graph when joining chrM. This is a legitimate and expected data
   situation, and it deserves a real diagnostic rather than a parse error.

Observed failure

Command (one chromosome per Slurm array task):

cactus-graphmap-join jobstores/jobstore-cactus-chrM \
    --vg pangenome/chrom-alignments/chrM.vg \
    --hal pangenome/chrom-alignments/chrM.hal \
    --outDir multiref-vcfs.tmp/cactus/chrM \
    --outName t2t_h9_hs.chrM \
    --reference T2T H9 HS REF4 \
    --vcfReference T2T H9 HS REF4 \
    --vcf \
    --maxCores 32 \
    --binariesMode local

There is no sample named T2T in the graph. The actual sample names are CHM13,
GRCh38, HG002_MAT, HG002_PAT, RPE1_HAP1, RPE1_HAP2, H9_HAP1, H9_HAP2.

Relevant log, abridged:

[I] [toil-rt] Running the command: "gfaffix .../chrM.vg.gfa --output_refined .../chrM.vg.gfaffixed.gfa --check_transformation --threads 1 --dont_collapse T2T#[.]*"
[I] [toil-rt] Successfully ran: "gfaffix ..." in 0.0124 seconds
[I] [toil-rt] Running the command: "bash -c set -eo pipefail && clip-vg .../chrM.vg.gfaffixed -f -e T2T -d _MINIGRAPH_ -L | vg mod -X 1024 - | vg paths -x - -n -Q T2T -t 1 | vg clip -d 1 - -P T2T | vg ids -s -"
[C] [toil.worker] Worker crashed with traceback:
  File ".../cactus/refmap/cactus_graphmap_join.py", line 1084, in clip_vg
    cactus_call(parameters=cmd, outfile=clipped_path, job_memory=job.memory)
  File ".../cactus/shared/common.py", line 916, in cactus_call
    raise RuntimeError("{}Command {} exited {}: {}".format(sigill_msg, call, process.returncode, out))
RuntimeError: Command ... exited 1: stderr=warning[vg paths] no matching paths found in graph
error[VPKG::load_one]: Correct input type not found in standard input while loading handlegraph::MutablePathMutableHandleGraph
error[VPKG::load_one]: Correct input type not found in standard input while loading handlegraph::MutablePathMutableHandleGraph

Why the error is unhelpful

Every stage in the chain tolerates the unmatched name, so the failure surfaces two tools
away from its cause:

• gfaffix --dont_collapse T2T#[.]* — an unmatched regex is not an error, exits 0.
• clip-vg -e T2T — the "don't clip the reference" prefix matches nothing, so every
  path is dropped. Exits 0, emits a graph with no paths.
• vg paths -x - -n -Q T2T — prints warning[vg paths] no matching paths found in graph
  (a warning, not an error) and writes nothing to stdout.
• vg clip -P T2T and vg ids -s - — get an empty stream and fail to deserialize a
  handlegraph. This is the exit status that pipefail propagates, and its message is
  what the user sees.

So the reported error describes a serialization problem in the last tool of the pipe,
while the actual cause is a command-line argument validated nowhere.

The cost is amplified by when it happens. This is a per-chromosome Slurm array over 25
tasks; each one queues, starts a container, imports its .vg/.hal, and only then dies
in clip_vg. A bad reference name wastes the whole array. Everything needed to reject the
name is available before the first job is scheduled.

Root cause in the source

graphmap_join_validate_options() does call check_sample_names()
(cactus/refmap/cactus_graphmap_join.py:186 in v3.0.1):

check_sample_names(options.reference, options.reference[0])

but it passes options.reference as the sample_names argument. So the existence test in
check_sample_names() (cactus/refmap/cactus_minigraph.py:237) —

if references[0] not in sample_names:
    raise RuntimeError("Specified reference, \"{}\" not in seqfile".format(references[0]))

— is checking the reference list against itself and is trivially satisfied. The suffix and
.-in-name checks still do useful work, but the "does this sample exist" check is
effectively a no-op here.

That is understandable: cactus-graphmap-join has no seqfile, so unlike
cactus-minigraph / cactus-pangenome it has no list of sample names to validate
against. But it does have the input graphs, and the sample names are recoverable from
them. Note also that --vcfReference is validated, against --reference
(cactus_graphmap_join.py:315-317), which makes the missing check on --reference itself
more surprising.

Suggested fix

Two independent improvements; either alone would have made this self-diagnosing.

1. Validate reference names against the input graphs up front. After the graphs are
imported and before any clip_vg job is scheduled, collect the sample names from each
input graph (vg paths -x <vg> -L, or vg paths -M, taking the field before the first
#) and check the --reference list against the union:

• A reference present in no input graph is unambiguously a mistake — fail immediately
  with the name and the sample names that were found:

  Specified --reference "T2T" has no paths in any input graph.
  Sample names found: CHM13, GRCh38, HG002_MAT, HG002_PAT, RPE1_HAP1, RPE1_HAP2, H9_HAP1, H9_HAP2
  

• A reference present in some graphs but absent from others is the chrY / chrM case.
  Failing is defensible, but the message must say which graph and which reference, e.g.:

  Specified --reference "H9_HAP1" has no paths in pangenome/chrom-alignments/chrY.vg
  (samples in that graph: CHM13, GRCh38, HG002_PAT). Either omit this reference or omit this graph.
  

  Better still would be to skip that reference for that graph and warn, since a sample
  genuinely has no chrY/chrM coordinate system to build a VCF against and there is nothing
  the user can do about it other than partition their inputs.

2. Make clip_vg fail with its own diagnostic. Independent of (1), when the pipeline
in clip_vg produces a graph with no reference paths, that is worth detecting directly
rather than letting an empty stream reach vg clip. Checking whether
clip-vg -e <ref> retained any paths, and raising

No paths matching reference "<ref>" in <graph>; found samples: ...

would keep the diagnostic attached to the graph and the reference that caused it, rather
than to the last tool in the pipe.

Reproduction

Any minigraph-cactus per-chromosome .vg will do:

cactus-graphmap-join ./js \
    --vg chrom-alignments/chr1.vg \
    --hal chrom-alignments/chr1.hal \
    --outDir out --outName test \
    --reference NO_SUCH_SAMPLE \
    --vcf --maxCores 4

Expected: an immediate error naming NO_SUCH_SAMPLE and listing the samples in the graph.
Actual: the clip_vg job runs and fails with error[VPKG::load_one]: Correct input type not found in standard input.