From 1e161b05cf617a99d287963d9a4abe5ce08d05d0 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Wed, 16 Sep 2026 11:11:35 -0400 Subject: [PATCH 1/2] fix --mgSplitWholeGenomeRef option bug --- doc/pangenome.md | 2 +- src/cactus/refmap/cactus_pangenome.py | 20 +++++++++++++------- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/doc/pangenome.md b/doc/pangenome.md index 69c5f4c12..cd932c088 100644 --- a/doc/pangenome.md +++ b/doc/pangenome.md @@ -176,7 +176,7 @@ cactus-align js ./out-map/chromfile.gm.txt out-align --reference simChimp --outV cactus-graphmap-join js --vg out-align/*.vg --hal out-align/*.hal --sv-gfa out-construct/*.gfa.gz --reference simChimp --outDir out-join --outName ep --gb ``` -Splitting this way costs some inter-chromosome context: a region homologous to several chromosomes (the acrocentric short arms, say) has nothing to compete against in a chromosome-level graph, where the whole-genome pipeline would have filtered it out as ambiguous. `cactus-pangenome --mgSplitWholeGenomeRef` (which implies `--mgSplit`, and has no step-by-step equivalent) builds each chromosome's second-pass minigraph against the whole reference genome(s) instead, restoring that competition and pruning the off-chromosome material back out before `cactus-align`, at the cost of indexing a whole reference per chromosome. +Splitting this way costs some inter-chromosome context: a region homologous to several chromosomes (the acrocentric short arms, say) has nothing to compete against in a chromosome-level graph, where the whole-genome pipeline would have filtered it out as ambiguous. `cactus-pangenome --mgSplitWholeGenomeRef` (which implies `--mgSplit`, and has no step-by-step equivalent) builds each chromosome's second-pass minigraph against the whole primary reference genome instead, restoring that competition and pruning the off-chromosome material back out before `cactus-align`, at the cost of indexing a whole reference per chromosome. ### Pipeline diff --git a/src/cactus/refmap/cactus_pangenome.py b/src/cactus/refmap/cactus_pangenome.py index 513a24b1e..c171ab8dc 100644 --- a/src/cactus/refmap/cactus_pangenome.py +++ b/src/cactus/refmap/cactus_pangenome.py @@ -69,7 +69,7 @@ def pangenome_options(parser): parser.add_argument("--mgSplit", action="store_true", default=False, help = "Run minigraph construction and mapping independently on each chromosome") parser.add_argument("--mgSplitWholeGenomeRef", action="store_true", default=False, - help = "Implies --mgSplit, and builds each chromosome's second-pass minigraph against the whole reference genome(s) rather than just that chromosome, so off-chromosome mappings can compete and be filtered the way they are in the whole-genome pipeline. The off-chromosome material is pruned back out before cactus-align.") + help = "Implies --mgSplit, and builds each chromosome's second-pass minigraph against the whole primary reference genome rather than just that chromosome, so off-chromosome mappings can compete and be filtered the way they are in the whole-genome pipeline. The off-chromosome material is pruned back out before cactus-align.") parser.add_argument("--inGFA", type=str, default=None, help = "Start from this existing minigraph GFA (.sv.gfa.gz from a previous run, or a published release) " "rather than building one. If the seqFile names genomes the graph does not have, they are constructed into it and " @@ -413,11 +413,17 @@ def split_reference_ids(job, seq_id_map, references): """ separate the whole-genome reference fastas from everything else, so --mgSplitWholeGenomeRef can hold them back from the post-split cleanup and build the second-pass minigraphs against them. - every --reference goes in, matching the first pass: --refOnly builds that graph from all of them, - and it is the graph the chromosome bins were decided against. only reference[0] is rank-0, but a - secondary reference still carries sequence the primary lacks, and minigraph maps against the whole - graph -- so it contributes competition for exactly the off-chromosome material this is here to - catch, and leaving it out would let the two passes disagree about what a chromosome contains. + only --reference[0] belongs here, and putting a secondary reference in is a decomposition bug. + everything the split produces is binned to exactly one chromosome; reference[0] is exempt only + because it *defines* the chromosomes, entering minigraph first and so carrying rank 0 everywhere. + a whole-genome input at rank > 0 is not chromosome-assigned, and rgfa-split bins its segments by + the reference contig they align to -- so CHM13's chr13 sequence lands in chr21's subproblem on + acrocentric homology while its own chr13 subproblem keeps the rest. the same contig then exists + in two per-chromosome graphs, at whichever rank the mash ordering gave that construct, and + merge_sv_gfa (which renumbers segment ids but not ranks) emits it with conflicting SR:i:. + observed on a 459-genome GRCh38 run: 5 CHM13 contigs, e.g. CHM13#0#chr13 at rank 192 in chr13 + and 335 in chr21. samples never do this -- they come in pre-split by chromosome, and the + off-chromosome part of a sample contig is dropped rather than kept in a second subproblem. they are already sanitized, and they must also bypass sanitize_fasta_headers_batch below: that deletes its input (sanitize_fasta_header in checkUniqueHeaders.py), and every chromosome shares @@ -714,7 +720,7 @@ def pangenome_end_to_end_workflow(job, options, config_wrapper, seq_id_map, seq_ wg_ref_id_map = None clean_seq_id_map = seq_id_map if options.mgSplitWholeGenomeRef: - split_ref_job = split_export_job.addFollowOnJobFn(split_reference_ids, seq_id_map, options.reference) + split_ref_job = split_export_job.addFollowOnJobFn(split_reference_ids, seq_id_map, options.reference[:1]) wg_ref_id_map, clean_seq_id_map = split_ref_job.rv(0), split_ref_job.rv(1) split_export_job = split_ref_job From 2243ff190fc66432f75aae43848fa3591cb70d45 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Wed, 16 Sep 2026 14:39:56 -0400 Subject: [PATCH 2/2] merge_pafs: don't name parameters things toil eats 9b0ce2d9 passed name= and cores= to merge_pafs_sized, but toil's FunctionWrappingJob pops memory/cores/disk/accelerators/preemptible/checkpoint/ name out of the kwargs before calling the function, using "name" as the job's unitName. Both were swallowed: the merged file kept the default 'merged.paf' for the GAF merge too, and the bgzip ran on one thread, which is what that commit was meant to fix. Rename them merged_name and gzip_cores; gzip_cores is still passed on to the child as cores=, where toil eating it is the point. Visible in an HPRC log as jobs literally named "merged.gaf" and "merged.paf" either side of `bgzip .../merged.paf --threads 1`. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01XScUpHRSBsui9Y67uQxA6a --- src/cactus/refmap/cactus_graphmap.py | 29 ++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/src/cactus/refmap/cactus_graphmap.py b/src/cactus/refmap/cactus_graphmap.py index f5ef7fcec..d53e1207b 100644 --- a/src/cactus/refmap/cactus_graphmap.py +++ b/src/cactus/refmap/cactus_graphmap.py @@ -566,9 +566,10 @@ def minigraph_map_all(job, options, config, gfa_id, fa_id_map, graph_event, in_g merge_name = getattr(options, 'mg_chrom_name', None) if options.batch else None merge_name = merge_name if merge_name else 'merged' paf_merge_job = top_job.addFollowOnJobFn(merge_pafs_sized, paf_id_map, - name='{}.paf'.format(merge_name)) + merged_name='{}.paf'.format(merge_name)) gaf_merge_job = top_job.addFollowOnJobFn(merge_pafs_sized, gaf_id_map, gzip=True, - name='{}.gaf'.format(merge_name), cores=mg_cores) + merged_name='{}.gaf'.format(merge_name), + gzip_cores=mg_cores) return paf_merge_job.rv(), gaf_merge_job.rv() @@ -987,24 +988,32 @@ def stable_gaf_to_paf(job, config, gaf_path, gfa_path, regranulated=False): # return the stable gaf (minigraph output) and the unstable paf return job.fileStore.writeGlobalFile(pansn_gaf_path), job.fileStore.writeGlobalFile(unstable_paf_path) -def merge_pafs(job, paf_file_id_map, gzip=False, name=None): - """ merge up some pafs. name is what the merged file is called on disk: getLocalTempFile() would - give it an anonymous .tmp, which is all anyone reading the log of the bgzip below would see """ +def merge_pafs(job, paf_file_id_map, gzip=False, merged_name=None): + """ merge up some pafs. merged_name is what the merged file is called on disk: getLocalTempFile() + would give it an anonymous .tmp, which is all anyone reading the log of the bgzip below would see. + + it is NOT called "name": toil's FunctionWrappingJob pops memory/cores/disk/accelerators/ + preemptible/checkpoint/name out of the kwargs of addChildJobFn before the function is called, and + uses "name" as the job's unitName. a parameter with any of those names is silently swallowed """ paf_paths = [job.fileStore.readGlobalFile(paf_id) for paf_id in paf_file_id_map.values()] - merged_path = os.path.join(job.fileStore.getLocalTempDir(), name if name else 'merged.paf') + merged_path = os.path.join(job.fileStore.getLocalTempDir(), merged_name if merged_name else 'merged.paf') catFiles(paf_paths, merged_path) if gzip: cactus_call(parameters=['bgzip', merged_path, '--threads', str(job.cores)]) merged_path += '.gz' return job.fileStore.writeGlobalFile(merged_path) -def merge_pafs_sized(job, paf_file_id_map, gzip=False, name=None, cores=1): +def merge_pafs_sized(job, paf_file_id_map, gzip=False, merged_name=None, gzip_cores=1): """ merge_pafs, sized off its inputs. callers upstream of the mapping jobs hold promises and so cannot measure them; by the time this job runs they are resolved. the merge holds every input - plus the merged copy, and bgzip then writes a compressed copy alongside """ + plus the merged copy, and bgzip then writes a compressed copy alongside. + + gzip_cores is not called "cores" for the reason in merge_pafs above: toil would eat it as this + job's own resource request and the function would keep its default, leaving the bgzip on one + thread. it is passed on as the child's cores, where toil eating it is exactly what we want """ total_size = sum(paf_id.size for paf_id in paf_file_id_map.values() if paf_id) - return job.addChildJobFn(merge_pafs, paf_file_id_map, gzip=gzip, name=name, - cores=cores, disk=max(total_size * 3, 2**31)).rv() + return job.addChildJobFn(merge_pafs, paf_file_id_map, gzip=gzip, merged_name=merged_name, + cores=gzip_cores, disk=max(total_size * 3, 2**31)).rv() def extract_paf_from_gfa(job, gfa_id, gfa_path, ref_event, graph_event, ignore_paf_id): """ make a paf directly from the rGFA tags. rgfa2paf supports other ranks, but we're only