fix --mgSplitWholeGenomeRef option bug

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
 

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

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 (<outName>.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