From 46015bd0cda2e16878e0b74ed30e25cbb60b015f Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Fri, 11 Sep 2026 16:27:43 -0400 Subject: [PATCH 1/4] Add genomes to an existing pangenome instead of rebuilding it A Minigraph-Cactus pangenome has been immutable: adding one genome to an HPRC-scale graph meant rebuilding all of it. minigraph construction is iterative in the input genomes and dominates the wall time of a large run -- weeks, for a release -- so that is the cost this removes. Adding 50 genomes to a graph of 450 now costs 50 genomes' worth of construction. cactus-pangenome --extendGFA seeds construction with an existing graph and --extendGAF reuses that run's mappings. cactus-minigraph --extendGFA and cactus-graphmap --extendGAF/--remap are the same thing at the step-by-step level, and the two interoperate in both directions: a graph built either way extends either way, and an extended pangenome is an ordinary one that can be extended again. Neither half needed a new algorithm. minigraph_construct_in_batches already chains "minigraph -cxggs prev.gfa a.fa b.fa" in batches of 50, so extending is seeding batch 0 with a graph somebody else wrote. The mappings did not need translating at all. cactus-graphmap runs minigraph without --vc, so the GAF it publishes as .gaf.gz is raw minigraph output in *stable* coordinates -- rGFA SN/SO names and offsets -- and it is published before any filtering. Adding genomes never moves those: nodes are appended and existing ones are only ever split, so the stable sequence a node covers stays where it was. The node-space PAF cactus consumes is derived from that GAF afterwards, by gaf2unstable | gaffilter | gaf2paf, so re-expressing 450 genomes against an augmented graph is running those three commands again with the new GFA, and gaf2unstable resolves the stable coordinates into the new, finer node ids on the way past. That chain is now stable_gaf_to_paf(), shared verbatim by mapping and by reuse. One thing does not survive the change of granularity. gaf2paf reads a record's path start as an offset into its *first* step, which holds for a GAF gaf2unstable resolved against the graph it was mapped to, where each of minigraph's stable steps is one node. Against a graph that has since been extended, the same stable step resolves into the finer nodes it was split into and the offset can reach past the first of them, which gaf2paf asserts on. On yeast that is 11 of SK1's 52 records; primates never hits it. trim_unstable_gaf() takes off the end steps that hold none of the alignment and moves the offsets with them -- a no-op, byte for byte, when the graph has not changed. It belongs in gaf2unstable, which is the thing changing the granularity; it is here so that this does not wait on a cactus-gfa-tools release. Two things make an extended pangenome differ from one built all at once, and it is worth being precise about which is which. Construction order. Building A B C D in one go sorts all four by mash distance; building A C and extending with B D gives A C B D, because the genomes already in the graph cannot be reordered around the ones being added. There is no way around that. It is also the only construction-side difference: with minigraphSortInput="none" so both sides use the same order, extending and building in one go produce a byte-identical GFA. Reused mappings. The genomes being added are mapped against the whole extended graph, so they come out exactly as a from-scratch run maps them. The ones already in the graph do not, so they never see nodes contributed by the new arrivals. Largely self-limiting -- a node is only there because some genome carries that allele, and an existing genome carrying it would have contributed it when the graph was built -- but not nothing. Measured against a from-scratch mapping on the same graph: primates (2 extended by 2) puts all 1183 alignments at identical coordinates with 5 cigars differing by a 2bp indel shift; yeast (3 extended by 3) agrees on 97.8% of alignments, with the reused genomes' aligned bases differing by at most 0.09%. --remap removes this difference entirely and reproduces a from-scratch PAF and GAF byte for byte, at the cost of the mapping stage -- which, unlike construction, is embarrassingly parallel. --mgSplit and --collapse are rejected with explicit errors. --mgSplit has per-chromosome graphs and mappings that would need extending as well; --collapse self-alignments come from minimap2 rather than from the GAF, so there is nothing in the GAF to reuse. A graph whose reference is not already in it is refused too: it would otherwise be constructed in last, at the highest rGFA rank rather than rank 0. Verification, all run locally: 24 offline unit tests covering the PanSN round trip on a published GAF, splitting a merged GAF back into the per-genome pieces it was concatenated from, and trim_unstable_gaf's invariants extending by nothing reproduces the GFA, PAF, GAF and minigraph fasta byte for byte -- what pins the whole reuse path -- plus an uncompressed output path, since an unchanged graph is the one graph construction hands back without writing it and so the one that can arrive at the wrong compression for where it is going extending and building in one go produce the same graph, the reused mappings land at from-scratch coordinates, and --remap reproduces the from-scratch PAF exactly the primates pipeline end to end, against the same MAF accuracy baseline as the from-scratch runs the yeast pipeline end to end, which is the one that puts a translated PAF through cactus-graphmap-split, checked against the same pinned graph statistics as the from-scratch six-strain run Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_0189H4mFGx16gLTMmqhLc2ku --- Makefile | 13 + doc/pangenome.md | 107 ++++++ src/cactus/paf/last_scoring.py | 15 +- src/cactus/refmap/cactus_graphmap.py | 323 ++++++++++++++++- .../refmap/cactus_graphmap_extendTest.py | 304 ++++++++++++++++ src/cactus/refmap/cactus_graphmap_split.py | 2 +- src/cactus/refmap/cactus_minigraph.py | 209 +++++++++-- src/cactus/refmap/cactus_pangenome.py | 49 ++- src/cactus/refmap/cactus_panpatch.py | 8 + test/evolverTest.py | 330 ++++++++++++++++-- 10 files changed, 1283 insertions(+), 77 deletions(-) create mode 100644 src/cactus/refmap/cactus_graphmap_extendTest.py diff --git a/Makefile b/Makefile index 58bfd58ce..8fc779e66 100644 --- a/Makefile +++ b/Makefile @@ -93,6 +93,7 @@ testModules = \ preprocessor/checkPreprocessedSequenceTest.py \ preprocessor/lastzRepeatMasking/cactus_lastzRepeatMaskTest.py \ progressive/multiCactusTreeTest.py \ + refmap/cactus_graphmap_extendTest.py \ refmap/cactus_panpatchTest.py \ refmap/pangenome_exclusionsTest.py \ update/cactus_update_prepareTest.py @@ -230,6 +231,18 @@ evolver_test_primates_pangenome_steps_mgsplit_docker: all ${CWD}/test/primates-t PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=docker CACTUS_DOCKER_MODE=1 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testEvolverPrimatesPangenomeStepByStepSplitDocker +pangenome_extend_null_test_local: + PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testPangenomeExtendNullLocal + +pangenome_extend_construction_test_local: + PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testPangenomeExtendConstructionLocal + +yeast_test_extend_local: + PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testYeastPangenomeExtendLocal + +evolver_test_primates_pangenome_extend_local: all ${CWD}/test/primates-truth.maf + PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testEvolverPrimatesPangenomeExtendLocal + evolver_test_all_local: evolver_test_local evolver_test_prepare_toil evolver_test_decomposed_local evolver_test_prepare_no_outgroup_local evolver_test_poa_local evolver_test_refmap_local evolver_test_minigraph_local yeast_test_local: diff --git a/doc/pangenome.md b/doc/pangenome.md index 1d68e8571..a8701818f 100644 --- a/doc/pangenome.md +++ b/doc/pangenome.md @@ -267,6 +267,113 @@ For `--vgFilter`, the filter threshold is inferred from the `.dX.vg` filename pa Note: per-chromosome output options (`--chrom-vg`, `--chrom-og`, `--viz`, `--draw`) cannot be used with bypass options, as you already have those files from the previous run. Also, bypass options are not compatible with graphs that were originally built with `--collapse`. +### Adding Genomes to an Existing Pangenome + +`cactus-pangenome --extendGFA` adds genomes to a pangenome you have already built, instead of +rebuilding it from scratch: + +``` +cactus-pangenome ./js ./seqfile.txt --reference GRCh38 --outDir pg2 --outName pg \ + --extendGFA pg1/pg.sv.gfa.gz --extendGAF pg1/pg.gaf.gz +``` + +`seqfile.txt` lists **every** genome, the ones already in the graph as well as the ones being +added. Cactus works out which are which; the genomes that are already there are left exactly where +they are in the graph, and only the new ones are constructed in. A genome cannot be *removed* from +a minigraph, so leaving one out of the seqfile is an error rather than a way to drop it. + +This matters because `minigraph` construction is iterative in the input genomes and dominates the +wall time of a large run — weeks, for an HPRC-scale release. Adding 50 genomes to a graph of 450 +costs 50 genomes' worth of construction, not 500. + +The step-by-step interface has the same two options, and they interoperate with the one-shot one in +both directions: + +``` +cactus-minigraph ./js ./seqfile.txt pg2/pg.sv.gfa.gz --reference GRCh38 \ + --extendGFA pg1/pg.sv.gfa.gz +cactus-graphmap ./js ./seqfile.txt pg2/pg.sv.gfa.gz pg2/pg.paf --reference GRCh38 \ + --outputFasta pg2/pg.sv.gfa.fa.gz --extendGAF pg1/pg.gaf.gz +``` + +Everything after these two stages — `cactus-graphmap-split`, `cactus-align`, and especially +`cactus-graphmap-join`'s `vg` indexing — is recomputed in full, and the output of an extended run +is an ordinary pangenome that can itself be extended again. + +#### What `--extendGAF` does + +`cactus-graphmap` runs `minigraph` without `--vc`, so the GAF it publishes as `.gaf.gz` is +in *stable* coordinates: rGFA `SN`/`SO` sequence names and offsets. Adding genomes to a graph never +moves those. New nodes are appended and existing ones are only ever split, so the stable sequence a +node covers stays exactly where it was. + +That means an existing genome's mappings do not have to be recomputed against the extended graph — +they can be re-derived from the published GAF by the same `gaf2unstable` / `gaf2paf` conversion that +produced the PAF in the first place, which resolves the stable coordinates into the new, finer node +ids for free. `--extendGAF` is the option that does this, and it turns the mapping stage from hours +of `minigraph` into minutes of file conversion. + +`--extendGAF` is optional. Without it — extending a published release for which only the GFA is +available, say — every genome is mapped again, which is the fallback and costs a full mapping stage. + +#### Why an extended pangenome is not identical to one built all at once + +There are exactly two sources of difference, and it is worth being precise about which is which. + +**1. Construction order.** `minigraph` construction is iterative in the input genomes, so the order +they go in decides the graph. Building `A B C D` in one go sorts all four by mash distance to the +reference; building `A C` and then extending with `B D` gives the effective order `A C B D`, +because the genomes already in the graph cannot be reordered around the ones being added. There is +no way around this, and it is why an extended graph is not the graph you would have got from +scratch. + +It is the *only* construction-side difference, though: with the order pinned +(`minigraphSortInput="none"` in the config), extending and building in one go produce a **byte +identical** graph. That equivalence is what `make pangenome_extend_construction_test_local` checks, +on `A B` + `C D` versus `A B C D`. + +**2. Reused mappings.** The genomes being *added* are mapped against the whole extended graph, so +they are mapped exactly as a from-scratch run would map them. The genomes already in the graph are +not: their alignments are the ones they had against the graph before it, re-expressed, so they +never see nodes contributed by the newly added genomes. + +This is largely self-limiting — a node is only in the graph because some genome carries that +allele, and an existing genome carrying it would have contributed it when the graph was first built +— but it is not nothing, and where it bites is near-identical alleles, where an existing genome +might now prefer a new genome's node. Measured against a from-scratch mapping on the very same +graph: + +| | agreement | +|---|---| +| evolver primates, 2 extended by 2 | all 1183 alignments at identical coordinates; 5 CIGARs differ by a 2bp indel shift | +| yeast, 3 strains extended by 3 | 97.8% of alignments at identical coordinates; the reused genomes' aligned bases differ by at most 0.09% | + +In both, the genomes that were *added* come out identical to a from-scratch mapping, because they +are mapped against the whole extended graph. All the divergence is in the reused ones. + +`--remap` removes this second difference entirely: every genome is mapped against the extended +graph, which reproduces a from-scratch mapping byte for byte, leaving construction order as the +only thing that differs. It costs the full mapping stage — which, unlike construction, is +embarrassingly parallel. + +#### Other caveats + +* The graph must have been built with a compatible configuration (`minigraphConstructOptions`, the + `` `assemblyName`) and with the same `--reference`. +* `--mgSplit` and `--collapse` are not supported with `--extendGFA`. `--mgSplit` has per-chromosome + graphs and mappings that would need extending as well; `--collapse` self-alignments come from + `minimap2` rather than from the GAF, so there is nothing in the GAF to reuse. +* Standalone `cactus-graphmap --extendGAF` still imports and sanitizes every genome's FASTA even + though the reused ones are not mapped. On the `cactus-pangenome` path that work is not wasted — + `cactus-align` needs those FASTAs anyway. +* A reused GAF needs one adjustment beyond re-running the conversion, handled by + `trim_unstable_gaf` in `cactus_graphmap.py`: `gaf2paf` reads a record's path start as an offset + into its *first* step, and splitting a node makes the steps finer while leaving the offset where + it was, so the offset can end up past the first of the nodes that replaced it. The steps it + reaches past hold none of the alignment, so they are taken off the front and back. This belongs + in `gaf2unstable`, which is the thing changing the granularity; it is done in cactus so that + reusing mappings does not wait on a `cactus-gfa-tools` release. + ### VCF Output The `--vcf` option runs `vg deconstruct` to represent the graph as sites of variation along a reference. A single run can write several VCFs, all prefixed with `--outName`: diff --git a/src/cactus/paf/last_scoring.py b/src/cactus/paf/last_scoring.py index d63c155da..b770d3f87 100644 --- a/src/cactus/paf/last_scoring.py +++ b/src/cactus/paf/last_scoring.py @@ -137,12 +137,19 @@ def apply_scores_to_config(score_dict, config_xml): poa_node.attrib['partialOrderAlignmentGapExtensionPenalty2'], poa_node.attrib['partialOrderAlignmentSubMatrix'])) -def last_train(job, config, seq_order, seq_id_map): - """ run last_train on a pair of fasta files, using the first as the database """ +def last_train(job, config, seq_order, seq_id_map, ref_name=None): + """ run last_train on a pair of fasta files, using the first as the database. - assert len(seq_order) > 1 + ref_name names the database genome when it is not seq_order[0], as when extending an existing + minigraph where the order holds only the genomes being added """ - name1 = seq_order[0] + if ref_name is None: + assert len(seq_order) > 1 + ref_name = seq_order[0] + elif ref_name not in seq_order: + seq_order = [ref_name] + list(seq_order) + + name1 = ref_name name2 = None # short circuit if ref sequence is too small to have a hope of training diff --git a/src/cactus/refmap/cactus_graphmap.py b/src/cactus/refmap/cactus_graphmap.py index 1e87e1250..3cf2c1784 100644 --- a/src/cactus/refmap/cactus_graphmap.py +++ b/src/cactus/refmap/cactus_graphmap.py @@ -8,6 +8,7 @@ """ import os, sys, re +import gzip from argparse import ArgumentParser import xml.etree.ElementTree as ET import copy @@ -57,6 +58,14 @@ def main(): parser.add_argument("--mapCores", type=int, help = "Number of cores for minigraph. Overrides graphmap cpu in configuration") parser.add_argument("--collapse", help = "Incorporate minimap2 self-alignments.", action='store_true', default=False) parser.add_argument("--collapseRefPAF", help ="Incorporate given (reference-only) self-alignments in PAF format [Experimental]") + parser.add_argument("--extendGAF", type=str, default=None, + help = "Reuse the mappings in this GAF (as published by a previous cactus-graphmap or cactus-pangenome run) " + "instead of re-running minigraph for the genomes it covers. Minigraph GAF is in stable coordinates, which " + "node splitting does not change, so these mappings are simply re-derived against the given (extended) graph. " + "Only genomes that are not in the GAF are mapped. Intended for use with cactus-minigraph --extendGFA") + parser.add_argument("--remap", action="store_true", default=False, + help = "Map every genome with minigraph even if --extendGAF already covers it. Slower, but the existing " + "genomes then see the nodes contributed by the newly added ones, as they would in a from-scratch run") parser.add_argument("--batch", action="store_true", help="Run independently on set of chromosomea inputs (chromfile as from cactus-minigraph --batch). Note that the output will be a directory and not a PAF") @@ -104,6 +113,15 @@ def main(): if options.mgSplit and options.batch: raise RuntimeError("--mgSplit is for the whole-genome splitting pass and cannot be used with --batch") + + if options.extendGAF: + if options.batch: + raise RuntimeError("--extendGAF cannot be used with --batch") + if options.collapse or options.collapseRefPAF: + raise RuntimeError("--extendGAF cannot be used with --collapse or --collapseRefPAF: collapse PAFs are minimap2 " + "self-alignments and are not derived from the GAF") + elif options.remap: + raise RuntimeError("--remap only means something with --extendGAF, which is what it overrides") # Mess with some toil options to create useful defaults. cactus_override_toil_options(options) @@ -230,9 +248,13 @@ def graph_map(options): input_dict[chrom] = seq_id_map, gfa_id, ref_collapse_paf_id, input_map[chrom][0], input_map[chrom][1] + #import the mappings to reuse + extend_gaf_id = toil.importFile(makeURL(options.extendGAF)) if options.extendGAF and not options.remap else None + # run the workflow # output_dict is chrom -> paf_id, gfa_fa_id, gaf_id, unfiltered_paf_id, paf_filter_log, paf_was_filtered - output_dict = toil.start(Job.wrapJobFn(minigraph_batch_separate_workflow, options, config_wrapper, input_dict, graph_event, True)) + output_dict = toil.start(Job.wrapJobFn(minigraph_batch_separate_workflow, options, config_wrapper, input_dict, graph_event, True, + extend_gaf_id=extend_gaf_id)) export_graphmap_output(options, config_node, input_map, output_dict, toil) @@ -293,7 +315,7 @@ def export_graphmap_output(options, config_node, input_map, output_dict, toil): if chrom_file_path.startswith('s3://'): write_s3(chrom_file_temp_path, chrom_file_path) -def minigraph_batch_workflow(job, options, config, input_dict, graph_event, sanitize, pansn_gfa_input=True): +def minigraph_batch_workflow(job, options, config, input_dict, graph_event, sanitize, pansn_gfa_input=True, extend_gaf_id=None): """ Batch wrapper to run grpahmap independently at the chromosome level.""" output_dict = {} options.mg_chrom_name = None @@ -308,7 +330,7 @@ def minigraph_batch_workflow(job, options, config, input_dict, graph_event, sani else: chrom_options = options mgwf_job = job.addChildJobFn(minigraph_workflow, chrom_options, config, seq_id_map, gfa_id, graph_event, - sanitize, ref_collapse_paf_id, pansn_gfa_input) + sanitize, ref_collapse_paf_id, pansn_gfa_input, extend_gaf_id=extend_gaf_id) output_dict[chrom] = mgwf_job.rv() return output_dict @@ -324,14 +346,15 @@ def add_separate_ref_contigs_job(batch_job, options, config, input_dict): getattr(options, 'permissiveContigFilter', None), whole_genome_ref=getattr(options, 'mgSplitWholeGenomeRef', False)) -def minigraph_batch_separate_workflow(job, options, config, input_dict, graph_event, sanitize, pansn_gfa_input=True): +def minigraph_batch_separate_workflow(job, options, config, input_dict, graph_event, sanitize, pansn_gfa_input=True, extend_gaf_id=None): """ minigraph_batch_workflow followed by the separation pass, for callers that just want the final result and add nothing after it """ batch_job = job.addChildJobFn(minigraph_batch_workflow, options, config, input_dict, graph_event, sanitize, - pansn_gfa_input) + pansn_gfa_input, extend_gaf_id=extend_gaf_id) return add_separate_ref_contigs_job(batch_job, options, config, input_dict).rv() -def minigraph_workflow(job, options, config, seq_id_map, gfa_id, graph_event, sanitize, ref_collapse_paf_id, pansn_gfa_input=True): +def minigraph_workflow(job, options, config, seq_id_map, gfa_id, graph_event, sanitize, ref_collapse_paf_id, pansn_gfa_input=True, + extend_gaf_id=None): """ Overall workflow takes command line options and returns (paf-id, (optional) fa-id) """ fa_id = None gfa_id_size = gfa_id.size @@ -365,7 +388,15 @@ def minigraph_workflow(job, options, config, seq_id_map, gfa_id, graph_event, sa new_root_job = Job() root_job.addFollowOn(new_root_job) root_job = new_root_job - gfa_id = rename_gfa_job.rv() + gfa_id = rename_gfa_job.rv(0) + + # split up any mappings we've been given to reuse, so each genome's re-derivation is its own + # job just as its mapping would have been + extend_gaf_map = None + if extend_gaf_id: + split_gaf_job = root_job.addChildJobFn(split_gaf_by_event, extend_gaf_id, genome_names, options.extendGAF, + disk=12*extend_gaf_id.size) + extend_gaf_map = split_gaf_job.rv() zipped_gfa = options.minigraphGFA.endswith('.gz') if options.outputFasta: @@ -381,7 +412,7 @@ def minigraph_workflow(job, options, config, seq_id_map, gfa_id, graph_event, sa gfa_id = gfa_unzip_job.rv() gfa_id_size *= 10 options.minigraphGFA = options.minigraphGFA[:-3] - paf_job = Job.wrapJobFn(minigraph_map_all, options, config, gfa_id, seq_id_map, graph_event) + paf_job = Job.wrapJobFn(minigraph_map_all, options, config, gfa_id, seq_id_map, graph_event, extend_gaf_map) root_job.addFollowOn(paf_job) collapse_paf_id = ref_collapse_paf_id @@ -469,8 +500,11 @@ def make_minigraph_fasta(job, gfa_file_id, gfa_file_path, name): return job.fileStore.writeGlobalFile(fa_path) -def minigraph_map_all(job, options, config, gfa_id, fa_id_map, graph_event): - """ top-level job to run the minigraph mapping in parallel, returns paf """ +def minigraph_map_all(job, options, config, gfa_id, fa_id_map, graph_event, extend_gaf_map=None): + """ top-level job to run the minigraph mapping in parallel, returns paf. + + a genome that extend_gaf_map already has mappings for has its PAF re-derived from them rather + than being mapped again -- see translate_gaf_one() """ # hang everything on this job, to self-contain workflow top_job = Job() job.addChild(top_job) @@ -488,6 +522,12 @@ def minigraph_map_all(job, options, config, gfa_id, fa_id_map, graph_event): # it must be gated on batch as well as the option: the option's own first pass maps whole-genome # queries against a whole-genome graph, where the query anchor still holds and 2x is right gfa_coefficient = 6 if options.batch and getattr(options, 'mgSplitWholeGenomeRef', False) else 2 + + # every genome whose contigs can name a step in the GAF's paths, which is more than the + # genomes being mapped: --refFromGFA takes the reference out of the sequence map + genome_names = set(fa_id_map.keys()) + if options.reference: + genome_names.add(options.reference if type(options.reference) is str else options.reference[0]) for event, fa_id in fa_id_map.items(): mem = 72*fa_id.size + gfa_coefficient*gfa_id.size event_name = event @@ -495,11 +535,19 @@ def minigraph_map_all(job, options, config, gfa_id, fa_id_map, graph_event): # the memory heuristc seems to drastically underestimate some chromosomes in batch mode... mem *= 2 event_name = '{}.{}'.format(event, options.mg_chrom_name) - minigraph_map_job = top_job.addChildJobFn(minigraph_map_one, config, event_name, fa_id, gfa_id, - cores=mg_cores, disk=5*fa_id.size + gfa_id.size, - memory=cactus_clamp_memory(mem)) - gaf_id_map[event] = minigraph_map_job.rv(0) - paf_id_map[event] = minigraph_map_job.rv(1) + if extend_gaf_map and event in extend_gaf_map: + # no minigraph, and no input fasta: gaf2unstable/gaffilter/gaf2paf against the new graph + # is the whole job. gaffilter reads its input into memory, as it does when mapping + gaf_shard_id = extend_gaf_map[event] + map_job = top_job.addChildJobFn(translate_gaf_one, config, event_name, gaf_shard_id, gfa_id, genome_names, + disk=12*gaf_shard_id.size + 2*gfa_id.size, + memory=cactus_clamp_memory(24*gaf_shard_id.size + 4*gfa_id.size)) + else: + map_job = top_job.addChildJobFn(minigraph_map_one, config, event_name, fa_id, gfa_id, + cores=mg_cores, disk=5*fa_id.size + gfa_id.size, + memory=cactus_clamp_memory(mem)) + gaf_id_map[event] = map_job.rv(0) + paf_id_map[event] = map_job.rv(1) # merge up. these two are the merges whose inputs scale with the number of genomes, so they get # sized off them rather than taking the default; the GAF one also bgzips, so give it the mapping @@ -519,6 +567,120 @@ def minigraph_map_all(job, options, config, gfa_id, fa_id_map, graph_event): # name that happens to contain id=...|, rewriting a name that exists in no graph and no input gaf_pansn_re = re.compile(r'(^|[\t><])id=([^|\t\n<>]+)\|') +def pansn_to_event_map(names): + """ SAMPLE#HAP -> seqfile event, for the genomes in names. event_to_pansn_prefix is lossy + (both S288C and S288C.0 give S288C#0), so going back needs the event names to hand """ + prefix_map = {} + for event in names: + prefix_map.setdefault(event_to_pansn_prefix(event), event) + return prefix_map + +# SAMPLE#HAP, as it appears in a published (PanSN) GAF's query column and in each of its path +# segments. anchored like gaf_pansn_re above, and stopping at the second '#' so that a PanSN +# phase block (SAMPLE#HAP#CONTIG#PHASEBLOCK) is left in the contig part where it belongs +pansn_gaf_re = re.compile(r'(^|[\t><])([^|\t\n<>#]+#[^|\t\n<>#]+)#') + +def gaf_from_pansn(names, gaf_path, out_path): + """ the inverse of gaf_to_pansn(): rewrite a published GAF's PanSN SAMPLE#HAP#CONTIG names back + to cactus's id=EVENT|CONTIG, so it can be resolved against a cactus-named GFA again. + + a prefix that is not one of the seqfile's genomes is left alone rather than mangled, which + covers both an already-cactus-named GAF (from a cactus old enough to have published one) and + any contig name that happens to look like a PanSN prefix """ + prefix_map = pansn_to_event_map(names) + + def replace(m): + event = prefix_map.get(m.group(2)) + if event is None: + return m.group(0) + return '{}id={}|'.format(m.group(1), event) + + with open(gaf_path, 'r') as in_file, open(out_path, 'w') as out_file: + for line in in_file: + out_file.write(pansn_gaf_re.sub(replace, line)) + +def split_gaf_by_event(job, gaf_id, names, gaf_path): + """ split a published (merged) GAF into one file per genome, returning {event: file id} """ + work_dir = job.fileStore.getLocalTempDir() + local_gaf_path = os.path.join(work_dir, 'extend.gaf.gz' if gaf_path.endswith('.gz') else 'extend.gaf') + job.fileStore.readGlobalFile(gaf_id, local_gaf_path) + shard_dir = os.path.join(work_dir, 'shards') + os.makedirs(shard_dir) + + shard_paths, dropped = split_gaf_file_by_event(local_gaf_path, names, shard_dir) + + if dropped: + RealtimeLogger.info('Ignoring mappings in {} for {}name(s) not in the seqfile: {}'.format( + gaf_path, 'at least ' if len(dropped) >= MAX_DROPPED_NAMES_REPORTED else '{} '.format(len(dropped)), + ' '.join(sorted(dropped)))) + RealtimeLogger.info('Reusing mappings for {} genome(s) from {}'.format(len(shard_paths), gaf_path)) + + return {event: job.fileStore.writeGlobalFile(shard_path) for event, shard_path in shard_paths.items()} + +# how many unrecognised names a split reports before it stops collecting them. the genome part of +# a name is normally one of a handful, but a name in neither naming falls back to its contig +MAX_DROPPED_NAMES_REPORTED = 20 + +def split_gaf_file_by_event(gaf_path, names, shard_dir): + """ split a published (merged) GAF into one file per genome, returning ({event: path}, dropped names). + + the merged GAF is a concatenation of the per-genome files, so this puts each genome's mappings + back exactly as minigraph_map_one() left them -- which is what lets the reused mappings be + re-derived by the very same code that produced them in the first place. + + genomes in the GAF that are not in names are dropped, not an error: --refFromGFA legitimately + takes the reference out of the sequence map. cactus-minigraph --extendGFA is where a genome + missing from the seqfile is caught, because there it is unrecoverable """ + prefix_map = pansn_to_event_map(names) + + def event_of(query_name): + """ the seqfile event a GAF query column belongs to, or None """ + if query_name.startswith('id='): + barpos = query_name.find('|') + event = query_name[3:barpos] if barpos > 3 else None + return event if event in names else None + hashpos = query_name.find('#') + if hashpos < 0: + return None + hashpos2 = query_name.find('#', hashpos + 1) + if hashpos2 < 0: + return None + return prefix_map.get(query_name[:hashpos2]) + + shard_paths = {} + dropped = set() + # the merged GAF groups each genome's records together, so one handle at a time is enough. + # append mode keeps it correct even if some other producer interleaved them + cur_event, cur_file = None, None + opener = gzip.open if gaf_path.endswith('.gz') else open + with opener(gaf_path, 'rt') as gaf_file: + for line in gaf_file: + tab = line.find('\t') + if tab < 0: + continue + query_name = line[:tab] + event = event_of(query_name) + if event is None: + # report the genome part of the name, in whichever naming it is in. a name that + # carries no genome at all falls back to the whole contig, so this is bounded by + # contig count rather than genome count and needs a cap of its own + if len(dropped) < MAX_DROPPED_NAMES_REPORTED: + dropped.add(query_name[3:query_name.find('|')] if query_name.startswith('id=') and '|' in query_name + else query_name.split('#')[0]) + continue + if event != cur_event: + if cur_file: + cur_file.close() + if event not in shard_paths: + shard_paths[event] = os.path.join(shard_dir, '{}.gaf'.format(event)) + cur_file = open(shard_paths[event], 'a') + cur_event = event + cur_file.write(line) + if cur_file: + cur_file.close() + + return shard_paths, dropped + def gaf_to_pansn(gaf_path, out_path): """ rewrite cactus's internal id=EVENT|CONTIG names as PanSN SAMPLE#HAP#CONTIG @@ -571,6 +733,119 @@ def minigraph_map_one(job, config, event_name, fa_file_id, gfa_file_id): cactus_call(parameters=cmd, job_memory=job.memory) + return stable_gaf_to_paf(job, config, gaf_path, gfa_path) + +def translate_gaf_one(job, config, event_name, gaf_file_id, gfa_file_id, genome_names): + """ Re-derive one genome's PAF from mappings it already has, against a (possibly extended) graph. + + minigraph GAF is in stable coordinates -- rGFA SN/SO names and offsets -- which adding genomes + to a graph does not change: new nodes are appended and existing ones are only ever split, so + the stable sequence a node covers stays exactly where it was. That makes re-deriving the PAF a + matter of running the same gaf2unstable/gaf2paf chain minigraph_map_one() runs, against the new + graph, which gaf2unstable resolves into the new (finer) node ids for free. + + the graph the genome was originally mapped to is not needed, and neither is minigraph """ + + work_dir = job.fileStore.getLocalTempDir() + gfa_path = os.path.join(work_dir, "mg.gfa") + gaf_path = os.path.join(work_dir, "{}.gaf".format(event_name)) + job.fileStore.readGlobalFile(gfa_file_id, gfa_path) + + # the published GAF is PanSN, but gaf2unstable resolves it against the cactus-named GFA. + # gaf_from_pansn passes through anything already in cactus naming. note the input cannot be + # named .pansn: that is where stable_gaf_to_paf() writes the copy it publishes + in_gaf_path = os.path.join(work_dir, "{}.in.gaf".format(event_name)) + job.fileStore.readGlobalFile(gaf_file_id, in_gaf_path) + gaf_from_pansn(genome_names, in_gaf_path, gaf_path) + + # the reused GAF is published back (PanSN in, PanSN out, unchanged), so a run that extends a + # pangenome can itself be extended + return stable_gaf_to_paf(job, config, gaf_path, gfa_path, regranulated=True) + +# a GAF path step: an orientation mark followed by a name that runs to the next mark +gaf_step_re = re.compile(r'[<>][^<>]+') + +def trim_unstable_gaf(gaf_path, out_path, node_lengths_path): + """ drop the path steps of each record that carry none of its alignment, moving the path + offsets along with them. + + gaf2paf reads a record's path start as an offset into its *first* step. That holds for a GAF + gaf2unstable resolved against the graph it was mapped to, where each of minigraph's stable + steps is one node. Against a graph that has since been extended, the same stable step resolves + into the several finer nodes it was split into, and the offset can now reach past the first of + them -- which gaf2paf asserts on rather than handles. + + The steps it reaches past hold no aligned bases, and gaf2paf emits nothing for them even when it + does cope, so taking them off the front and back restores the shape gaf2paf expects without + changing the alignment at all. When nothing needs trimming -- every mapping that was made + against the graph it is being resolved against -- every line is passed through untouched. + + This belongs in gaf2unstable, which is the thing changing the granularity; it lives here so + that reusing mappings does not depend on a new cactus-gfa-tools release. """ + node_len = {} + with open(node_lengths_path) as lengths_file: + for line in lengths_file: + toks = line.split() + if len(toks) >= 2: + node_len[toks[0]] = int(toks[1]) + + def first_step_len(path): + end = path.find('>', 1) + alt = path.find('<', 1) + if alt != -1 and (end == -1 or alt < end): + end = alt + return node_len[path[1:] if end == -1 else path[1:end]] + + def last_step_len(path): + start = max(path.rfind('>'), path.rfind('<')) + return node_len[path[start + 1:]] + + trimmed_records = 0 + with open(gaf_path) as in_file, open(out_path, 'w') as out_file: + for line in in_file: + toks = line.rstrip('\n').split('\t') + if len(toks) < 12 or not toks[5] or toks[5][0] not in '<>': + out_file.write(line) + continue + path, path_len, path_start, path_end = toks[5], int(toks[6]), int(toks[7]), int(toks[8]) + # the overwhelming majority of records need nothing done, and deciding that needs only + # the two end steps: with the offsets inside them, nothing in between can be outside + if path_start < first_step_len(path) and path_end > path_len - last_step_len(path): + out_file.write(line) + continue + steps = gaf_step_re.findall(path) + lo, hi = 0, len(steps) + while lo < hi - 1 and node_len[steps[lo][1:]] <= path_start: + dropped = node_len[steps[lo][1:]] + path_start -= dropped + path_end -= dropped + path_len -= dropped + lo += 1 + while hi - 1 > lo and path_len - node_len[steps[hi - 1][1:]] >= path_end: + path_len -= node_len[steps[hi - 1][1:]] + hi -= 1 + if lo == 0 and hi == len(steps): + # nothing was outside the alignment after all, so the record stands as it is + out_file.write(line) + continue + toks[5], toks[6], toks[7], toks[8] = ''.join(steps[lo:hi]), str(path_len), str(path_start), str(path_end) + out_file.write('\t'.join(toks) + '\n') + trimmed_records += 1 + + return trimmed_records + +def stable_gaf_to_paf(job, config, gaf_path, gfa_path, regranulated=False): + """ Turn a stable-coordinate (ie minigraph output) GAF into the node-coordinate PAF cactus + consumes, returning (published PanSN gaf id, paf id). Shared by mapping and by reuse of an + existing mapping, so that the two produce identical output for identical input. + + regranulated says the GAF was made against a coarser version of this graph, so its path offsets + have to be brought back inside their first and last steps -- see trim_unstable_gaf(). It is a + no-op when the graph has not changed, but it is only asked for on the reuse path so that + mapping keeps running exactly the commands it always has """ + + xml_node = findRequiredNode(config.xmlRoot, "graphmap") + # convert the gaf into unstable gaf (targets are node sequences) # note: the gfa needs to be uncompressed for this tool to work mg_lengths_path = gfa_path + '.node_lengths.tsv' @@ -586,7 +861,23 @@ def minigraph_map_one(job, config, event_name, fa_file_id, gfa_file_id): if overlap_ratio: cmd = [cmd, ['gaffilter', '-', '-r', str(overlap_ratio), '-m', str(length_ratio), '-q', str(min_mapq), '-b', str(min_block), '-i', str(min_ident)]] - cactus_call(parameters=cmd, outfile=unstable_gaf_path, job_memory=job.memory) + try: + cactus_call(parameters=cmd, outfile=unstable_gaf_path, job_memory=job.memory) + except RuntimeError as e: + if not regranulated: + raise + # gaf2unstable asserts, rather than reporting, when a record names stable sequence the + # graph does not have. Mapping cannot reach that -- the GAF came from this graph -- but + # reuse can, and the assertion on its own says nothing about which input is wrong + raise RuntimeError('Failed to resolve reused mappings against this graph. If the GAF names sequence the graph ' + 'does not have, it was made against a different pangenome: the GAF must come from the run ' + 'that produced the graph being extended. Underlying error: {}'.format(e)) + + if regranulated: + trimmed_path = unstable_gaf_path + '.trimmed' + trimmed = trim_unstable_gaf(unstable_gaf_path, trimmed_path, mg_lengths_path) + RealtimeLogger.info('Moved the path offsets of {} reused GAF record(s) back inside their end steps'.format(trimmed)) + os.replace(trimmed_path, unstable_gaf_path) # convert the unstable gaf into unstable paf, which is what cactus expects # also tack on the unique id to the target column diff --git a/src/cactus/refmap/cactus_graphmap_extendTest.py b/src/cactus/refmap/cactus_graphmap_extendTest.py new file mode 100644 index 000000000..a925192dc --- /dev/null +++ b/src/cactus/refmap/cactus_graphmap_extendTest.py @@ -0,0 +1,304 @@ +#!/usr/bin/env python3 + +""" +Unit tests for the pure (non-Toil) logic behind extending a pangenome: the PanSN <-> cactus name +round trip on a published GAF, and splitting a merged GAF back into the per-genome pieces it was +concatenated from. + +Both exist because minigraph GAF is in stable coordinates, which adding genomes to a graph does not +change, so a genome's existing mappings can be re-derived against the extended graph instead of +being recomputed. That only holds if the published GAF round trips exactly, which is what these +pin down. + +These are fast and offline. The end-to-end extension is covered by evolverTest.py. +""" + +import gzip +import os +import re +import tempfile +import unittest + +from cactus.refmap.cactus_graphmap import ( + gaf_from_pansn, gaf_to_pansn, pansn_to_event_map, split_gaf_file_by_event, trim_unstable_gaf) + + +# a stable GAF line as minigraph writes it: query name, then a path of stable segments with +# orientation marks, then the tags. only columns 1 and 6 carry sequence names +def gaf_line(query, path, tags='60\ttp:A:P\tcm:i:100'): + return '{}\t1000\t0\t900\t+\t{}\t2000\t10\t910\t880\t900\t{}\n'.format(query, path, tags) + + +class TestPansnToEventMap(unittest.TestCase): + + def test_haploid_and_diploid(self): + self.assertEqual(pansn_to_event_map({'S288C', 'HG002.1', 'HG002.2'}), + {'S288C#0': 'S288C', 'HG002#1': 'HG002.1', 'HG002#2': 'HG002.2'}) + + def test_explicit_hap_zero_beats_bare_name(self): + # event_to_pansn_prefix maps both S288C and S288C.0 onto S288C#0. whichever wins, the map + # has to be single valued or the round trip would be ambiguous + prefix_map = pansn_to_event_map({'S288C', 'S288C.0'}) + self.assertEqual(len(prefix_map), 1) + self.assertIn(prefix_map['S288C#0'], ('S288C', 'S288C.0')) + + def test_non_numeric_suffix_is_part_of_the_sample(self): + # HG002.pat is not a haplotype suffix, so the whole thing is the sample + self.assertEqual(pansn_to_event_map({'HG002.pat'}), {'HG002.pat#0': 'HG002.pat'}) + + +class TestGafPansnRoundTrip(unittest.TestCase): + + def round_trip(self, names, cactus_gaf): + """ cactus -> PanSN -> cactus, which is the path a reused GAF actually takes """ + with tempfile.TemporaryDirectory() as work_dir: + in_path = os.path.join(work_dir, 'in.gaf') + pansn_path = os.path.join(work_dir, 'pansn.gaf') + back_path = os.path.join(work_dir, 'back.gaf') + with open(in_path, 'w') as in_file: + in_file.write(cactus_gaf) + gaf_to_pansn(in_path, pansn_path) + gaf_from_pansn(names, pansn_path, back_path) + with open(pansn_path) as pansn_file, open(back_path) as back_file: + return pansn_file.read(), back_file.read() + + def test_forward_and_back(self): + names = {'S288C', 'SK1'} + gaf = gaf_line('id=SK1|chrI', '>id=S288C|chrI>id=S288C|chrII') + pansn, back = self.round_trip(names, gaf) + self.assertIn('SK1#0#chrI', pansn) + self.assertIn('>S288C#0#chrI>S288C#0#chrII', pansn) + self.assertEqual(back, gaf) + + def test_reverse_steps(self): + names = {'HG002.1', 'HG002.2', 'GRCh38'} + gaf = gaf_line('id=HG002.2|chr1', 'id=HG002.1|chr1HG002#1#chr1id=S288C|chrI:100-2000') + pansn, back = self.round_trip(names, gaf) + self.assertIn('>S288C#0#chrI:100-2000', pansn) + self.assertEqual(back, gaf) + + def test_contig_name_containing_a_hash(self): + # PanSN phase blocks put a fourth '#' field on a path name. going back, only the first two + # '#' belong to the prefix + names = {'HG002.1', 'GRCh38'} + gaf = gaf_line('id=HG002.1|chr1#0', '>id=GRCh38|chr1#0') + pansn, back = self.round_trip(names, gaf) + self.assertIn('HG002#1#chr1#0', pansn) + self.assertEqual(back, gaf) + + def test_contig_name_containing_a_pipe(self): + names = {'S288C', 'SK1'} + gaf = gaf_line('id=SK1|ctg|1', '>id=S288C|chrI') + pansn, back = self.round_trip(names, gaf) + self.assertEqual(back, gaf) + + def test_already_cactus_named_gaf_passes_through(self): + # an older cactus published the GAF in its own naming; reading one back must not mangle it + names = {'S288C', 'SK1'} + gaf = gaf_line('id=SK1|chrI', '>id=S288C|chrI') + with tempfile.TemporaryDirectory() as work_dir: + in_path = os.path.join(work_dir, 'in.gaf') + out_path = os.path.join(work_dir, 'out.gaf') + with open(in_path, 'w') as in_file: + in_file.write(gaf) + gaf_from_pansn(names, in_path, out_path) + with open(out_path) as out_file: + self.assertEqual(out_file.read(), gaf) + + def test_unknown_prefix_is_left_alone(self): + # a contig that merely looks like a PanSN prefix belongs to no genome and must not be + # rewritten into a name that exists in no graph + gaf = gaf_line('weird#name#ctg', '>other#thing#ctg') + with tempfile.TemporaryDirectory() as work_dir: + in_path = os.path.join(work_dir, 'in.gaf') + out_path = os.path.join(work_dir, 'out.gaf') + with open(in_path, 'w') as in_file: + in_file.write(gaf) + gaf_from_pansn({'S288C'}, in_path, out_path) + with open(out_path) as out_file: + self.assertEqual(out_file.read(), gaf) + + +class TestSplitGafByEvent(unittest.TestCase): + + def split(self, gaf_text, names, gzipped=False): + work_dir = tempfile.mkdtemp() + gaf_path = os.path.join(work_dir, 'merged.gaf.gz' if gzipped else 'merged.gaf') + if gzipped: + with gzip.open(gaf_path, 'wt') as gaf_file: + gaf_file.write(gaf_text) + else: + with open(gaf_path, 'w') as gaf_file: + gaf_file.write(gaf_text) + shard_dir = os.path.join(work_dir, 'shards') + os.makedirs(shard_dir) + shard_paths, dropped = split_gaf_file_by_event(gaf_path, names, shard_dir) + contents = {} + for event, shard_path in shard_paths.items(): + with open(shard_path) as shard_file: + contents[event] = shard_file.read() + return contents, dropped + + def test_partition_is_exact(self): + # every line lands in exactly one shard, in order: this is what makes the re-derived PAF + # identical to the original when nothing has been added + names = {'S288C', 'SK1', 'Y12'} + lines = [gaf_line('S288C#0#chrI', '>S288C#0#chrI'), + gaf_line('S288C#0#chrII', '>S288C#0#chrII'), + gaf_line('SK1#0#chrI', '>S288C#0#chrI'), + gaf_line('Y12#0#chrI', '>S288C#0#chrI')] + contents, dropped = self.split(''.join(lines), names) + self.assertEqual(set(contents), names) + self.assertEqual(dropped, set()) + self.assertEqual(contents['S288C'], lines[0] + lines[1]) + self.assertEqual(contents['SK1'], lines[2]) + self.assertEqual(''.join(contents[e] for e in ['S288C', 'SK1', 'Y12']), ''.join(lines)) + + def test_gzipped_input(self): + names = {'S288C'} + line = gaf_line('S288C#0#chrI', '>S288C#0#chrI') + contents, dropped = self.split(line, names, gzipped=True) + self.assertEqual(contents, {'S288C': line}) + + def test_interleaved_records_still_group(self): + # append mode means the grouping does not depend on the concatenation order + names = {'S288C', 'SK1'} + a1 = gaf_line('S288C#0#chrI', '>S288C#0#chrI') + b1 = gaf_line('SK1#0#chrI', '>S288C#0#chrI') + a2 = gaf_line('S288C#0#chrII', '>S288C#0#chrII') + contents, _ = self.split(a1 + b1 + a2, names) + self.assertEqual(contents['S288C'], a1 + a2) + self.assertEqual(contents['SK1'], b1) + + def test_diploid_haplotypes_are_separate_shards(self): + names = {'HG002.1', 'HG002.2'} + h1 = gaf_line('HG002#1#chr1', '>GRCh38#0#chr1') + h2 = gaf_line('HG002#2#chr1', '>GRCh38#0#chr1') + contents, _ = self.split(h1 + h2, names) + self.assertEqual(contents, {'HG002.1': h1, 'HG002.2': h2}) + + def test_genome_not_in_seqfile_is_dropped_and_reported(self): + names = {'S288C'} + keep = gaf_line('S288C#0#chrI', '>S288C#0#chrI') + drop = gaf_line('GONE#0#chrI', '>S288C#0#chrI') + contents, dropped = self.split(keep + drop, names) + self.assertEqual(contents, {'S288C': keep}) + self.assertEqual(dropped, {'GONE'}) + + def test_cactus_named_gaf_splits_too(self): + names = {'S288C', 'SK1'} + a = gaf_line('id=S288C|chrI', '>id=S288C|chrI') + b = gaf_line('id=SK1|chrI', '>id=S288C|chrI') + contents, dropped = self.split(a + b, names) + self.assertEqual(contents, {'S288C': a, 'SK1': b}) + self.assertEqual(dropped, set()) + + +class TestTrimUnstableGaf(unittest.TestCase): + """ gaf2paf reads a record's path start as an offset into its first step. Extending a graph + splits nodes, so a reused mapping's offset can end up past the first of the finer nodes that + replaced the one it was made against -- which gaf2paf asserts on. Trimming the steps that hold + none of the alignment puts the offsets back where gaf2paf expects them. """ + + NODE_LENS = {'s1': 100, 's2': 50, 's3': 200, 's4': 30, 's5': 80} + + def trim(self, records): + """ run trim_unstable_gaf over (path, path_len, path_start, path_end) tuples, returning the + same tuples back """ + with tempfile.TemporaryDirectory() as work_dir: + lengths_path = os.path.join(work_dir, 'lens.tsv') + with open(lengths_path, 'w') as lengths_file: + for node, length in self.NODE_LENS.items(): + lengths_file.write('{}\t{}\n'.format(node, length)) + in_path = os.path.join(work_dir, 'in.gaf') + out_path = os.path.join(work_dir, 'out.gaf') + with open(in_path, 'w') as in_file: + for path, path_len, path_start, path_end in records: + in_file.write('q\t1000\t0\t900\t+\t{}\t{}\t{}\t{}\t880\t900\t60\tcg:Z:900M\n'.format( + path, path_len, path_start, path_end)) + trimmed = trim_unstable_gaf(in_path, out_path, lengths_path) + out = [] + with open(out_path) as out_file: + for line in out_file: + toks = line.rstrip('\n').split('\t') + out.append((toks[5], int(toks[6]), int(toks[7]), int(toks[8]))) + return out, trimmed + + def assert_gaf2paf_invariant(self, record): + """ what gaf2paf assumes: the start offset is inside the first step and the end offset is + inside the last, and the steps still add up to the stated path length """ + path, path_len, path_start, path_end = record + steps = [s[1:] for s in re.findall(r'[<>][^<>]+', path)] + node_lens = [self.NODE_LENS[s] for s in steps] + self.assertEqual(sum(node_lens), path_len) + self.assertLess(path_start, node_lens[0]) + self.assertGreater(path_end, path_len - node_lens[-1]) + + def test_nothing_to_trim_passes_through(self): + # the mapping was made against this very graph: every offset is already in its end step + record = ('>s1>s2>s3', 350, 40, 300) + out, trimmed = self.trim([record]) + self.assertEqual(out, [record]) + self.assertEqual(trimmed, 0) + + def test_offset_past_the_first_step(self): + # s1+s2 replaced one 150bp node, so an offset of 120 into it now lands in s2 + out, trimmed = self.trim([('>s1>s2>s3', 350, 120, 300)]) + self.assertEqual(trimmed, 1) + self.assertEqual(out, [('>s2>s3', 250, 20, 200)]) + self.assert_gaf2paf_invariant(out[0]) + + def test_alignment_ends_before_the_last_steps(self): + out, trimmed = self.trim([('>s3>s2>s1', 350, 10, 200)]) + self.assertEqual(trimmed, 1) + self.assertEqual(out, [('>s3', 200, 10, 200)]) + self.assert_gaf2paf_invariant(out[0]) + + def test_trims_both_ends(self): + out, trimmed = self.trim([('>s1>s2>s3>s4>s5', 460, 150, 350)]) + self.assertEqual(trimmed, 1) + # s1 and s2 are wholly before the start, s4 and s5 wholly after the end + self.assertEqual(out, [('>s3', 200, 0, 200)]) + self.assert_gaf2paf_invariant(out[0]) + + def test_offsets_keep_their_span(self): + # trimming moves the window, it must never resize it + for record in [('>s1>s2>s3', 350, 120, 300), ('>s1>s2>s3>s4>s5', 460, 150, 350)]: + out, _ = self.trim([record]) + self.assertEqual(out[0][3] - out[0][2], record[3] - record[2]) + + def test_reverse_steps_trim_the_same_way(self): + # the offsets run along the path as written, so orientation does not enter into it + out, trimmed = self.trim([('s1', 100, 100, 100)]) + self.assertEqual(out, [('>s1', 100, 100, 100)]) + self.assertEqual(trimmed, 0) + + def test_non_path_lines_pass_through(self): + # a stable path (no orientation marks) is not ours to touch + with tempfile.TemporaryDirectory() as work_dir: + lengths_path = os.path.join(work_dir, 'lens.tsv') + open(lengths_path, 'w').write('s1\t100\n') + in_path = os.path.join(work_dir, 'in.gaf') + out_path = os.path.join(work_dir, 'out.gaf') + text = gaf_line('q', 'chr1') + open(in_path, 'w').write(text) + self.assertEqual(trim_unstable_gaf(in_path, out_path, lengths_path), 0) + with open(out_path) as out_file: + self.assertEqual(out_file.read(), text) + + +if __name__ == '__main__': + unittest.main() diff --git a/src/cactus/refmap/cactus_graphmap_split.py b/src/cactus/refmap/cactus_graphmap_split.py index de0ce17b0..b65a1c392 100644 --- a/src/cactus/refmap/cactus_graphmap_split.py +++ b/src/cactus/refmap/cactus_graphmap_split.py @@ -230,7 +230,7 @@ def graphmap_split_workflow(job, options, config, seq_id_map, seq_name_map, gfa_ new_root_job = Job() root_job.addFollowOn(new_root_job) root_job = new_root_job - gfa_id = rename_gfa_job.rv() + gfa_id = rename_gfa_job.rv(0) # use file extension to sniff out compressed input if gfa_path.endswith(".gz"): diff --git a/src/cactus/refmap/cactus_minigraph.py b/src/cactus/refmap/cactus_minigraph.py index b5062d089..e82d23a56 100644 --- a/src/cactus/refmap/cactus_minigraph.py +++ b/src/cactus/refmap/cactus_minigraph.py @@ -57,6 +57,11 @@ def main(): help="Use last-train to estimate scoring matrix from input data", default=False) parser.add_argument("--refOnly", action="store_true", help="Only build the graph out of reference genome(s). Can be used when it will only be used for chromosome-splitting, for example") + parser.add_argument("--extendGFA", type=str, default=None, + help="Extend this existing minigraph GFA (as made by a previous cactus-minigraph or cactus-pangenome run) " + "instead of building from scratch. Only the seqFile genomes that are not already in the graph get added, in " + "mash-distance order among themselves. The seqFile must still contain every genome in the graph: genomes " + "cannot be removed from a minigraph") parser.add_argument("--batch", action="store_true", help="Run independently on set of chromosomea inputs (chromfile as from cactus-graphmap-split). Note that the output will be a directory and not a GFA") @@ -121,11 +126,22 @@ def main(): if '://' not in options.outputGFA: options.outputGFA = os.path.abspath(options.outputGFA) + extend_gfa_id = None + if options.extendGFA: + if options.batch: + raise RuntimeError('--extendGFA cannot be used with --batch') + if options.refOnly: + raise RuntimeError('--extendGFA cannot be used with --refOnly') + if '://' not in options.extendGFA: + options.extendGFA = os.path.abspath(options.extendGFA) + extend_gfa_id = toil.importFile(makeURL(options.extendGFA)) + # maps name -> input_seq_id_map, input_seq_order input_dict = minigraph_construct_import_sequences(options, config_wrapper, input_seqfiles, toil) # output_dict: chrom-> (gfa_id, pansn_gfa_id, train_id) - output_dict = toil.start(Job.wrapJobFn(minigraph_construct_batch_workflow, options, config_node, input_dict, options.outputGFA)) + output_dict = toil.start(Job.wrapJobFn(minigraph_construct_batch_workflow, options, config_node, input_dict, options.outputGFA, + extend_gfa_id=extend_gfa_id)) export_minigraph_construct_output(options, input_seqfiles, output_dict, toil) @@ -261,7 +277,7 @@ def check_sample_names(sample_names, references): raise RuntimeError("Sample name {} with \"{}\" suffix is not supported. You must either remove this suffix or use .N where N is an integer to specify haplotype".format(sample, sample_ext)) def minigraph_construct_batch_workflow(job, options, config_node, input_dict, gfa_path, sanitize=True, - construct_ref_id_map=None): + construct_ref_id_map=None, extend_gfa_id=None): """ run the construction workflow on individual chromosomes. construct_ref_id_map, if given, swaps the whole-genome reference fastas in for the chromosome's own slice of them (--mgSplit --mgSplitWholeGenomeRef). the merge happens here rather than inside minigraph_construct_workflow @@ -279,12 +295,12 @@ def minigraph_construct_batch_workflow(job, options, config_node, input_dict, gf else: gfa_path = options.outputGFA mgwf_job = job.addChildJobFn(minigraph_construct_workflow, options, config_node, seq_id_map, seq_order, gfa_path, sanitize, - construct_seq_id_map=construct_seq_id_map) + construct_seq_id_map=construct_seq_id_map, extend_gfa_id=extend_gfa_id) output_dict[chrom] = mgwf_job.rv() return output_dict def minigraph_construct_workflow(job, options, config_node, seq_id_map, seq_order, gfa_path, sanitize=True, - construct_seq_id_map=None): + construct_seq_id_map=None, extend_gfa_id=None): """ minigraph can handle bgzipped files but not gzipped; so unzip everything in case before running construct_seq_id_map, when given, replaces seq_id_map for the graph construction alone. it is how @@ -293,12 +309,58 @@ def minigraph_construct_workflow(job, options, config_node, seq_id_map, seq_orde reference slice, since training against a whole genome would be both expensive and cross-chromosome contaminated -- and would silently produce no model at all, as last_train() requires its partner sequence to be at least half the size of the database it trains against. - The construction job itself is sized off the substituted map, since that is what it runs on. """ + The construction job itself is sized off the substituted map, since that is what it runs on. + + with a graph to extend, which genomes still need constructing is not known until that graph's + SN tags have been read, so the rest of the workflow is deferred behind the job that reads them """ + if not extend_gfa_id: + return minigraph_construct_run(job, options, config_node, seq_id_map, seq_order, gfa_path, sanitize, + construct_seq_id_map=construct_seq_id_map) + + # the renaming pass decompresses the GFA before bgzipping it back up, so it needs room for + # the raw copy (reckoned at 10x, as elsewhere) on top of the compressed input and output + rename_job = job.addChildJobFn(minigraph_gfa_from_pansn, set(seq_id_map.keys()), options.extendGFA, extend_gfa_id, + disk=extend_gfa_id.size*12) + run_job = rename_job.addFollowOnJobFn(minigraph_construct_run, options, config_node, seq_id_map, seq_order, gfa_path, + sanitize, construct_seq_id_map, + rename_job.rv(0), rename_job.rv(1), extend_gfa_id) + return run_job.rv(0), run_job.rv(1), run_job.rv(2) + +def minigraph_construct_run(job, options, config_node, seq_id_map, seq_order, gfa_path, sanitize=True, + construct_seq_id_map=None, seed_gfa_id=None, seed_events=None, seed_pansn_gfa_id=None): assert type(options.reference) is list - assert options.reference[0] == seq_order[0] # the substituted map is a plain dict here, but sanitized_seq_id_map below is a promise when # sanitize is on, so the two can't be reconciled in this job assert not (construct_seq_id_map and sanitize) + ref_size = seq_id_map[options.reference[0]].size + # the PanSN rename at the end of construction has to resolve every SN tag in the finished + # graph, which on the extend path is more genomes than minigraph is being given + graph_names = set(seq_id_map.keys()) + if seed_events is not None: + if options.reference[0] not in seed_events: + # it would otherwise be constructed in last, at the highest rGFA rank rather than rank 0, + # and every rank-0 assumption downstream would be reading the wrong genome + raise RuntimeError('Reference {} is not in the graph being extended, whose genomes are: {}. A graph can only be ' + 'extended with the reference it was built on'.format(options.reference[0], + ' '.join(sorted(seed_events)))) + # everything already in the seed graph is left alone: minigraph only gets the genomes that + # are new to it, appended after the ones the graph was built from. the reference is kept in + # the sequence map (but not the order) because the mash sort below still sketches against it + seq_order = [seq for seq in seq_order if seq not in seed_events] + seq_id_map = {name: fa_id for name, fa_id in seq_id_map.items() + if name not in seed_events or name == options.reference[0]} + RealtimeLogger.info('Extending a graph of {} genomes with {}: {}'.format( + len(seed_events), len(seq_order), ' '.join(seq_order) if seq_order else '(nothing)')) + if not seq_order: + # nothing to add, so the graph handed back is the one --extendGFA was given. its + # compression follows the input name, and everything downstream reads the *output* + # name to decide whether to unzip, so it has to be re-emitted to match + match_job = job.addChildJobFn(match_gfa_compression, seed_gfa_id, seed_pansn_gfa_id, + options.extendGFA, gfa_path, + disk=12 * (seed_gfa_id.size if hasattr(seed_gfa_id, 'size') else 0)) + return match_job.rv(0), match_job.rv(1), None + else: + assert options.reference[0] == seq_order[0] if options.refOnly: refonly_seq_id_map = {} refonly_seq_order = [] @@ -309,7 +371,6 @@ def minigraph_construct_workflow(job, options, config_node, seq_id_map, seq_orde seq_id_map, seq_order = refonly_seq_id_map, refonly_seq_order if construct_seq_id_map: construct_seq_id_map = {seq: construct_seq_id_map[seq] for seq in refonly_seq_order} - ref_size = seq_id_map[options.reference[0]].size if sanitize: sanitize_job = job.addChildJobFn(sanitize_fasta_headers, seq_id_map, pangenome=True) sanitized_seq_id_map = sanitize_job.rv() @@ -320,7 +381,8 @@ def minigraph_construct_workflow(job, options, config_node, seq_id_map, seq_orde xml_node = findRequiredNode(config_node, "graphmap") sort_type = getOptionalAttrib(xml_node, "minigraphSortInput", str, default=None) if sort_type == "mash" and len(seq_id_map) > 2: - sort_job = sanitize_job.addFollowOnJobFn(sort_minigraph_input_with_mash, options, config_node, sanitized_seq_id_map, seq_order) + sort_job = sanitize_job.addFollowOnJobFn(sort_minigraph_input_with_mash, options, config_node, sanitized_seq_id_map, seq_order, + ref_name=options.reference[0] if seed_events is not None else None) seq_order = sort_job.rv() prev_job = sort_job else: @@ -328,12 +390,14 @@ def minigraph_construct_workflow(job, options, config_node, seq_id_map, seq_orde minigraph_job = prev_job.addFollowOnJobFn(minigraph_construct_in_batches, options, config_node, construct_seq_id_map if construct_seq_id_map else sanitized_seq_id_map, seq_order, gfa_path, - whole_genome_ref=bool(construct_seq_id_map)) + whole_genome_ref=bool(construct_seq_id_map), + seed_gfa_id=seed_gfa_id, graph_names=graph_names) train_id = None if options.lastTrain and len(seq_id_map) > 1: # note: somehow last training memory overruns don't seem to be detected by slurm so we # give 12G at least whenever possible, as --doubleMem won't help... - last_train_job = prev_job.addFollowOnJobFn(last_train, config_node, seq_order, sanitized_seq_id_map, + last_train_job = prev_job.addFollowOnJobFn(last_train, config_node, seq_order, sanitized_seq_id_map, + ref_name=options.reference[0] if seed_events is not None else None, cores=options.mgCores, disk=8*ref_size, memory=cactus_clamp_memory(max(8*ref_size, 12*10**9))) @@ -341,8 +405,40 @@ def minigraph_construct_workflow(job, options, config_node, seq_id_map, seq_orde return minigraph_job.rv(0), minigraph_job.rv(1), train_id -def sort_minigraph_input_with_mash(job, options, config_node, seq_id_map, seq_order): - """ Sort the input """ +def match_gfa_compression(job, gfa_id, pansn_gfa_id, in_path, out_path): + """ re-emit an unchanged seed graph at the compression its output path asks for. + + every path out of construction bgzips iff the output path ends in .gz, and the stages after it + read that same suffix to decide whether the file needs unzipping. Extending a graph by nothing + is the one case that returns a graph nobody wrote, so it is the one case that can arrive with + the wrong compression for where it is going. """ + want_gz = out_path.endswith('.gz') + if want_gz == in_path.endswith('.gz'): + return gfa_id, pansn_gfa_id + + work_dir = job.fileStore.getLocalTempDir() + out_ids = [] + for i, file_id in enumerate([gfa_id, pansn_gfa_id]): + in_gfa_path = os.path.join(work_dir, 'seed.{}.gfa{}'.format(i, '.gz' if not want_gz else '')) + job.fileStore.readGlobalFile(file_id, in_gfa_path) + out_gfa_path = os.path.join(work_dir, 'out.{}.gfa{}'.format(i, '.gz' if want_gz else '')) + if want_gz: + cactus_call(parameters=['bgzip', '--threads', str(job.cores), '-c', in_gfa_path], outfile=out_gfa_path) + else: + cactus_call(parameters=['gzip', '-dc', in_gfa_path], outfile=out_gfa_path) + out_ids.append(job.fileStore.writeGlobalFile(out_gfa_path)) + return out_ids[0], out_ids[1] + +def sort_minigraph_input_with_mash(job, options, config_node, seq_id_map, seq_order, ref_name=None): + """ Sort the input. + + ref_name names the genome to measure distance against when it is not seq_order[0], as when + extending an existing minigraph and the order holds only the genomes being added. It is put + back at the front for the sort and taken off again on the way out, so the genomes being added + are still ordered by their distance to the reference """ + trim_ref = ref_name is not None and ref_name not in seq_order + if trim_ref: + seq_order = [ref_name] + list(seq_order) # (dist, length) pairs which will be sorted decreasing on dist, breaking ties with increasing on length # assumption : reference is first mash_dists = [(0, sys.maxsize)] @@ -377,7 +473,7 @@ def sort_minigraph_input_with_mash(job, options, config_node, seq_id_map, seq_or disk = 2 * sum(seq_id_map[x].size for x in names) + seq_id_map[seq_order[0]].size).rv() dist_maps.append(dist_map) - return dist_root_job.addFollowOnJobFn(mash_distance_order, options, config_node, seq_order, dist_maps).rv() + return dist_root_job.addFollowOnJobFn(mash_distance_order, options, config_node, seq_order, dist_maps, trim_ref).rv() def mash_sketch(job, ref_seq, seq_id_map): """ get the sketch """ @@ -430,7 +526,7 @@ def parse_mash_output(mash_output): return output_dist_map -def mash_distance_order(job, options, config_node, seq_order, mash_output_maps): +def mash_distance_order(job, options, config_node, seq_order, mash_output_maps, trim_ref=False): """ get the sequence order from the mash distance""" # we first orient the list of dicts along seq_order @@ -472,6 +568,8 @@ def mash_distance_order(job, options, config_node, seq_order, mash_output_maps): fixed_order[empty_slots[j]] = seq j += 1 for ref in options.reference[1:]: + if ref not in mash_order: + continue mash_pos = mash_order.index(ref) fix_pos = fixed_order.index(ref) assert fix_pos == seq_order.index(ref) @@ -479,15 +577,24 @@ def mash_distance_order(job, options, config_node, seq_order, mash_output_maps): RealtimeLogger.info('Secondary reference {}, which would have mash rank {}, fixed at input rank {} because minigraphSortReference is disabled'.format(ref, mash_pos, fix_pos)) mash_order = fixed_order - return mash_order + return mash_order[1:] if trim_ref else mash_order -def minigraph_construct_in_batches(job, options, config_node, seq_id_map, seq_order, gfa_path, whole_genome_ref=False): - """ Make minigraph in sequential batches""" +def minigraph_construct_in_batches(job, options, config_node, seq_id_map, seq_order, gfa_path, whole_genome_ref=False, + seed_gfa_id=None, graph_names=None): + """ Make minigraph in sequential batches. + + seed_gfa_id is an existing graph to extend: it is fed to the first batch exactly as each batch + already feeds its output to the next, so extending a graph and constructing one in batches are + the same operation """ max_size = max([x.size for x in seq_id_map.values()]) total_size = sum([x.size for x in seq_id_map.values()]) - disk = total_size * 2 - mem = cactus_clamp_memory(60 * max_size + int(total_size / 4)) + # a seed graph can dwarf the genomes being added to it (one sample onto an HPRC-scale graph), + # so it needs to be in the estimates rather than lost in the headroom the way a batch-to-batch + # intermediate is + seed_size = seed_gfa_id.size if seed_gfa_id else 0 + disk = total_size * 2 + seed_size * 12 + mem = cactus_clamp_memory(60 * max_size + int(total_size / 4) + seed_size * 12) if whole_genome_ref: # with --mgSplitWholeGenomeRef the largest input is the whole reference, so the estimate above # is already the whole-genome one, calibrated against a graph carrying far more sample @@ -516,6 +623,11 @@ def minigraph_construct_in_batches(job, options, config_node, seq_id_map, seq_or assert num_batches > 0 and num_batches <= 991 prev_job = None prev_gfa_path = None + seed_gfa_path = None + if seed_gfa_id: + # minigraph_construct() only uses this to name its local copy, but keep the compression + # suffix honest since that is what says whether the file it reads is bgzipped + seed_gfa_path = 'extend.gfa.gz' if options.extendGFA.endswith('.gz') else 'extend.gfa' for i in range(num_batches): batch_size = len(seq_order) - i * max_batch_size if i == num_batches - 1 else max_batch_size input_seq_order = seq_order[i * max_batch_size : (i * max_batch_size) + batch_size] @@ -529,8 +641,9 @@ def minigraph_construct_in_batches(job, options, config_node, seq_id_map, seq_or out_gfa_path = '{}.{}'.format(gfa_path, i) pan_sn_output = False minigraph_job = Job.wrapJobFn(minigraph_construct, options, config_node, seq_id_map, input_seq_order, out_gfa_path, - prev_job.rv() if prev_job else None, prev_gfa_path, - pan_sn_output, + prev_job.rv() if prev_job else seed_gfa_id, + prev_gfa_path if prev_job else seed_gfa_path, + pan_sn_output, graph_names, disk=disk, memory=mem, cores=options.mgCores) if prev_job: prev_job.addFollowOn(minigraph_job) @@ -543,8 +656,13 @@ def minigraph_construct_in_batches(job, options, config_node, seq_id_map, seq_or return prev_job.rv() -def minigraph_construct(job, options, config_node, seq_id_map, seq_order, gfa_path, prev_gfa_id, prev_gfa_path, pan_sn_output): - """ Make minigraph """ +def minigraph_construct(job, options, config_node, seq_id_map, seq_order, gfa_path, prev_gfa_id, prev_gfa_path, pan_sn_output, + graph_names=None): + """ Make minigraph. + + graph_names is every genome the finished graph can contain, which is only seq_id_map's keys + when the graph is built from scratch: extending one leaves the genomes already in it out of + seq_id_map, but their SN tags still have to be renamed on the way out """ work_dir = job.fileStore.getLocalTempDir() gfa_path = os.path.join(work_dir, os.path.basename(gfa_path)) @@ -587,7 +705,7 @@ def minigraph_construct(job, options, config_node, seq_id_map, seq_order, gfa_pa if pan_sn_output: # rename to pan-sn before serializing, so it's more useful (ie for anything except cactus) pansn_gfa_path = os.path.join(work_dir, 'pan-sn.' + os.path.basename(gfa_path)) - minigraph_gfa_to_pansn(set(seq_id_map.keys()), gfa_path, pansn_gfa_path, job.cores) + minigraph_gfa_to_pansn(graph_names if graph_names else set(seq_id_map.keys()), gfa_path, pansn_gfa_path, job.cores) pansn_gfa_out_id = job.fileStore.writeGlobalFile(pansn_gfa_path) return gfa_out_id, pansn_gfa_out_id else: @@ -647,6 +765,13 @@ def minigraph_gfa_from_pansn(job, names, gfa_path, gfa_id): """ hack to convert PanSN names like simChimp#0#simpChimp.chr6 to Cactus names like id=simChimp.0|simChimp.chr6 so that a minigrpah GFA (as converted panSN by minigraph_gfa_to_pansn() above) can be read back into Cactus + returns (converted gfa id, set of genomes the graph was built from). the genome set is what + --extendGFA needs to work out which of the seqfile's genomes are new to the graph, and it comes + free with the pass that has to read every SN tag anyway. + + a GFA that is already in cactus naming -- from a cactus old enough to have published one, or + handed straight from one stage to the next -- is read for its genome set and returned as-is. + todo: Cactus should probably be changed to just use PanSN internally as well, but that's a much bigger lift """ @@ -657,12 +782,22 @@ def minigraph_gfa_from_pansn(job, names, gfa_path, gfa_id): in_file, out_file, raw_out_path = open_gfa_for_rename(gfa_path, out_gfa_path) + events = set() + unresolved = set() + already_cactus = False for line in in_file: line = line.decode() if line.startswith('S'): toks = line.strip().split('\t') for i, tok in enumerate(toks[4:]): if tok.startswith('SN:Z:'): + if tok.startswith('SN:Z:id='): + # already cactus-named: nothing to rewrite, just harvest the genome + already_cactus = True + barpos = tok.find('|') + if barpos > 8: + events.add(tok[8:barpos]) + break hashpos = tok.find('#') if hashpos < 0: # no prefix found: do nothing and hope for the best @@ -677,8 +812,17 @@ def minigraph_gfa_from_pansn(job, names, gfa_path, gfa_id): # minigraph_to_pansn() will add #0 to names without any dots # we untangle that here using the names list if name not in names: - name = '{}.{}'.format(name, hap) - assert name in names + if '{}.{}'.format(name, hap) in names: + name = '{}.{}'.format(name, hap) + else: + # collected rather than asserted on: with --extendGFA this is usually + # the user leaving a genome out of the seqfile, which deserves to be + # named. the whole tag goes in the message because the other way to + # land here is an SN tag that is not SAMPLE#HAP#CONTIG at all, and + # then the prefix alone says nothing about what is wrong + unresolved.add(tok) + break + events.add(name) toks[4+i] = 'SN:Z:id={}|{}'.format(name, tok[hashpos2+1:]) break out_file.write(('\t'.join(toks) + '\n').encode()) @@ -687,8 +831,19 @@ def minigraph_gfa_from_pansn(job, names, gfa_path, gfa_id): in_file.close() out_file.close() + + if unresolved: + raise RuntimeError('{} sequence name(s) in {} could not be matched to a seqfile genome: {}. Genomes cannot ' + 'be removed from a minigraph, so every genome in the graph must be in the seqfile -- ' + 'unless these names are not SAMPLE#HAP#CONTIG, in which case the graph was not written ' + 'by cactus and cannot be read back into it'.format( + len(unresolved), gfa_path, ' '.join(sorted(unresolved)[:10]))) + + if already_cactus: + return gfa_id, events + bgzip_gfa_rename(raw_out_path, out_gfa_path, job.cores) - return job.fileStore.writeGlobalFile(out_gfa_path) + return job.fileStore.writeGlobalFile(out_gfa_path), events diff --git a/src/cactus/refmap/cactus_pangenome.py b/src/cactus/refmap/cactus_pangenome.py index fed98be04..c62476a6f 100644 --- a/src/cactus/refmap/cactus_pangenome.py +++ b/src/cactus/refmap/cactus_pangenome.py @@ -70,6 +70,19 @@ def pangenome_options(parser): 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.") + parser.add_argument("--extendGFA", type=str, default=None, + help = "Add genomes to this existing pangenome's minigraph GFA (.sv.gfa.gz from a previous run, or a " + "published release) instead of building one from scratch. The seqFile must list every genome in the graph as " + "well as the ones being added: genomes cannot be removed from a minigraph. Only the new genomes are constructed " + "in, which is where nearly all of the minigraph cost is. Everything from cactus-graphmap-split on is recomputed") + parser.add_argument("--extendGAF", type=str, default=None, + help = "Reuse the mappings in this GAF (.gaf.gz from the same run that produced --extendGFA) rather than " + "re-running minigraph for the genomes it covers. Minigraph GAF is in stable coordinates, which adding genomes to " + "a graph does not change, so they are simply re-derived against the extended graph. Without this, every genome is " + "mapped again (see --remap)") + parser.add_argument("--remap", action="store_true", default=False, + help = "With --extendGAF, map every genome with minigraph anyway. Costs the full mapping stage, but the existing " + "genomes then see the nodes contributed by the newly added ones, as they would in a from-scratch run") # cactus-graphmap options parser.add_argument("--mapCores", type=int, help = "Number of cores for minigraph map. Overrides graphmap cpu in configuration") @@ -198,6 +211,20 @@ def pangenome_validate_options(options): if options.mgSplit and options.noSplit: raise RuntimeError('you cannot use both --mgSplit and --noSplit together: pick one') + if options.extendGFA: + if options.mgSplit: + raise RuntimeError('--extendGFA cannot (yet) be used with --mgSplit: the per-chromosome graphs and mappings would ' + 'need extending too') + if options.collapse or options.collapseRefPAF: + raise RuntimeError('--extendGFA cannot be used with --collapse or --collapseRefPAF: collapse PAFs are minimap2 ' + 'self-alignments and are not derived from the GAF') + else: + if options.extendGAF: + raise RuntimeError('--extendGAF requires --extendGFA: reusing mappings only makes sense against the graph they ' + 'were made from') + if options.remap: + raise RuntimeError('--remap only means something with --extendGAF, which is what it overrides') + # Sort out the graphmap-join options, which can be rather complex # pass in dummy values for now, they will get filled in later # (but we want to do as much error-checking upfront as possible) @@ -299,6 +326,17 @@ def main(): if options.scoresFile: last_scores_id = toil.importFile(makeURL(options.scoresFile)) + #import the pangenome being extended + extend_gfa_id, extend_gaf_id = None, None + if options.extendGFA: + if '://' not in options.extendGFA: + options.extendGFA = os.path.abspath(options.extendGFA) + extend_gfa_id = toil.importFile(makeURL(options.extendGFA)) + if options.extendGAF and not options.remap: + if '://' not in options.extendGAF: + options.extendGAF = os.path.abspath(options.extendGAF) + extend_gaf_id = toil.importFile(makeURL(options.extendGAF)) + #import the sequences input_seq_id_map = {} input_path_map = {} @@ -315,7 +353,8 @@ def main(): elif genome in input_seq_order: input_seq_order.remove(genome) - toil.start(Job.wrapJobFn(pangenome_end_to_end_workflow, options, config_wrapper, input_seq_id_map, input_path_map, input_seq_order, ref_collapse_paf_id, last_scores_id)) + toil.start(Job.wrapJobFn(pangenome_end_to_end_workflow, options, config_wrapper, input_seq_id_map, input_path_map, input_seq_order, ref_collapse_paf_id, last_scores_id, + extend_gfa_id=extend_gfa_id, extend_gaf_id=extend_gaf_id)) end_time = timeit.default_timer() run_time = end_time - start_time @@ -567,7 +606,7 @@ def export_join_wrapper(job, options, wf_output, contig_sizes_id=None): job.fileStore.exportFile(contig_sizes_id, makeURL(sizes_path)) def pangenome_end_to_end_workflow(job, options, config_wrapper, seq_id_map, seq_path_map, seq_order, ref_collapse_paf_id, - last_scores_id): + last_scores_id, extend_gfa_id=None, extend_gaf_id=None): """ chain the entire workflow together, doing exports after each step to mitigate annoyance of failures """ root_job = Job() job.addChild(root_job) @@ -617,7 +656,8 @@ def pangenome_end_to_end_workflow(job, options, config_wrapper, seq_id_map, seq_ else: split_config_node = config_node split_config_wrapper = config_wrapper - minigraph_job = prev_job.addFollowOnJobFn(minigraph_construct_workflow, mg_options, split_config_node, seq_id_map, seq_order, sv_gfa_path, sanitize=False) + minigraph_job = prev_job.addFollowOnJobFn(minigraph_construct_workflow, mg_options, split_config_node, seq_id_map, seq_order, sv_gfa_path, sanitize=False, + extend_gfa_id=extend_gfa_id) sv_gfa_id = minigraph_job.rv(0) pansn_sv_gfa_id = minigraph_job.rv(1) if not last_scores_id: @@ -634,7 +674,8 @@ def pangenome_end_to_end_workflow(job, options, config_wrapper, seq_id_map, seq_ gm_options = copy.deepcopy(options) if options.mgSplit: gm_options.collapse = False - graphmap_job = minigraph_wrapper_job.addFollowOnJobFn(minigraph_workflow, gm_options, split_config_wrapper, seq_id_map, sv_gfa_id, graph_event, False, ref_collapse_paf_id, pansn_gfa_input=False) + graphmap_job = minigraph_wrapper_job.addFollowOnJobFn(minigraph_workflow, gm_options, split_config_wrapper, seq_id_map, sv_gfa_id, graph_event, False, ref_collapse_paf_id, pansn_gfa_input=False, + extend_gaf_id=extend_gaf_id) paf_id, gfa_fa_id, gaf_id, unfiltered_paf_id, paf_filter_log = graphmap_job.rv(0), graphmap_job.rv(1), graphmap_job.rv(2), graphmap_job.rv(3), graphmap_job.rv(4) graphmap_export_job = graphmap_job.addFollowOnJobFn(export_graphmap_wrapper, options, paf_id, paf_path, gaf_id, unfiltered_paf_id, paf_filter_log) diff --git a/src/cactus/refmap/cactus_panpatch.py b/src/cactus/refmap/cactus_panpatch.py index d3e1762b8..06057c8a4 100644 --- a/src/cactus/refmap/cactus_panpatch.py +++ b/src/cactus/refmap/cactus_panpatch.py @@ -366,6 +366,14 @@ def panpatch_validate_options(options): # the unpatched input raise RuntimeError('--noSplit cannot be used with cactus-panpatch: panpatch needs one graph per reference chromosome') + # getattr because the unit tests build a minimal namespace rather than going through the parser + if getattr(options, 'extendGFA', None) or getattr(options, 'extendGAF', None) or getattr(options, 'remap', False): + # these come in via pangenome_options(), which panpatch shares. the graph panpatch builds + # is a throwaway, built per sample out of that sample and its donors, so there is nothing + # an earlier run could usefully be extended from + raise RuntimeError('--extendGFA / --extendGAF / --remap cannot be used with cactus-panpatch: it builds a fresh ' + 'graph per sample being patched, which is not an extension of anything') + def disable_pangenome_outputs(options): """ switch off every graphmap-join output except the chromosome vgs, so that we don't spend hours building indexes, VCFs and GFAs that get thrown away. must be called *after* diff --git a/test/evolverTest.py b/test/evolverTest.py index 0d4ee1b39..b825337e8 100644 --- a/test/evolverTest.py +++ b/test/evolverTest.py @@ -223,13 +223,39 @@ def _run_evolver_in_docker(self, seqFile = './examples/evolverMammals.txt'): sys.stderr.write('Running {}\n'.format(' '.format(cmd))) subprocess.check_call(' '.join(cmd), shell=True) - def _write_primates_seqfile(self, seq_file_path): - """ create the primates seqfile at given path""" + PRIMATES = ['simHuman', 'simChimp', 'simGorilla', 'simOrang'] + + def _write_primates_seqfile(self, seq_file_path, genomes=None): + """ create the primates seqfile at given path, optionally for a subset of the genomes """ + url = 'https://raw.githubusercontent.com/UCSantaCruzComputationalGenomicsLab/cactusTestData/master/evolver/primates/loci1/{}.chr6' with open(seq_file_path, 'w') as seq_file: - seq_file.write('simHuman\thttps://raw.githubusercontent.com/UCSantaCruzComputationalGenomicsLab/cactusTestData/master/evolver/primates/loci1/simHuman.chr6\n') - seq_file.write('simChimp\thttps://raw.githubusercontent.com/UCSantaCruzComputationalGenomicsLab/cactusTestData/master/evolver/primates/loci1/simChimp.chr6\n') - seq_file.write('simGorilla\thttps://raw.githubusercontent.com/UCSantaCruzComputationalGenomicsLab/cactusTestData/master/evolver/primates/loci1/simGorilla.chr6\n') - seq_file.write('simOrang\thttps://raw.githubusercontent.com/UCSantaCruzComputationalGenomicsLab/cactusTestData/master/evolver/primates/loci1/simOrang.chr6\n') + for genome in genomes if genomes is not None else self.PRIMATES: + seq_file.write('{}\t{}\n'.format(genome, url.format(genome))) + + def _write_seqfile_subset(self, in_path, out_path, genomes): + """ copy a seqfile keeping only the named genomes, in the order they appear in it """ + keep = set(genomes) + with open(in_path) as in_file, open(out_path, 'w') as out_file: + for line in in_file: + toks = line.split() + if toks and toks[0] in keep: + out_file.write('{}\t{}\n'.format(toks[0], toks[1])) + + def _write_pangenome_config(self, name, graphmap_attribs=None, graphmap_join_attribs=None): + """ write a copy of the default config with the given / + attribute overrides applied, and return its path """ + # use the same logic cactus does to get default config + config_path = 'src/cactus/cactus_progressive_config.xml' + xml_root = ET.parse(config_path).getroot() + for elem_name, attribs in [("graphmap", graphmap_attribs), ("graphmap_join", graphmap_join_attribs)]: + if attribs: + xml_root.find(elem_name).attrib.update(attribs) + out_path = os.path.join(self.tempDir, "config.{}.xml".format(name)) + with open(out_path, 'w') as config_file: + xmlString = ET.tostring(xml_root, encoding='unicode') + xmlString = minidom.parseString(xmlString).toprettyxml() + config_file.write(xmlString) + return out_path def _run_evolver_primates_star(self, binariesMode, configFile = None): """ Run cactus on the evolver primates with a star topology @@ -477,8 +503,12 @@ def _run_evolver_primates_graphmap(self, binariesMode): subprocess.check_call(['cactus-align', self._job_store(binariesMode), seq_file_fix_path, paf_path, self._out_hal(binariesMode), '--pangenome', '--outVG', '--outGFA', '--pafMaskFilter', '10000', '--barMaskFilter', '10000'] + cactus_opts) - def _run_evolver_primates_pangenome(self, binariesMode, mgSplit = False): - """ run the primates start in using high-level cactus-pangenome interface """ + def _run_evolver_primates_pangenome(self, binariesMode, mgSplit = False, extend = False): + """ run the primates start in using high-level cactus-pangenome interface. + + with extend, half the genomes are built into a graph with the step-by-step tools first and + cactus-pangenome --extendGFA adds the rest, which also checks that a graph made one way can + be extended the other """ # borrow seqfile from other primates test # todo: make a seqfile and add it to the repo seq_file_path = os.path.join(self.tempDir, 'primates.txt') @@ -491,26 +521,22 @@ def _run_evolver_primates_pangenome(self, binariesMode, mgSplit = False): # tack on a vcfwave test for docker binaries cactus_opts += ['--vcf', '--vcfReference', 'simChimp', '--clip', '1000', '--vcfwave'] - # use the same logic cactus does to get default config - config_path = 'src/cactus/cactus_progressive_config.xml' - xml_root = ET.parse(config_path).getroot() - graphmap_elem = xml_root.find("graphmap") - graphmap_join_elem = xml_root.find("graphmap_join") - # force cactus to use minigraph chunking - graphmap_elem.attrib["minigraphConstructBatchSize"] = "2" - # force cactus use vcfwave chunking - graphmap_join_elem.attrib["vcfwaveChunkLines"] = "1000" - mc_config_path = os.path.join(self.tempDir, "config.mc.xml") - with open(mc_config_path, 'w') as mc_config_file: - xmlString = ET.tostring(xml_root, encoding='unicode') - xmlString = minidom.parseString(xmlString).toprettyxml() - mc_config_file.write(xmlString) + # force cactus to use minigraph chunking, and vcfwave chunking + mc_config_path = self._write_pangenome_config('mc', + graphmap_attribs={"minigraphConstructBatchSize": "2"}, + graphmap_join_attribs={"vcfwaveChunkLines": "1000"}) cactus_opts += ['--configFile', mc_config_path] out_dir = os.path.dirname(self._out_hal(binariesMode)) out_name = os.path.splitext(os.path.basename(self._out_hal(binariesMode)))[0] cactus_pangenome_cmd = ['cactus-pangenome', self._job_store(binariesMode), seq_file_path, '--reference', 'simHuman', 'simChimp', - '--outDir', out_dir, '--outName', out_name, '--odgi', '--chrom-og', '--viz', '--draw', '--haplo', '--collapse', '--lastTrain'] + '--outDir', out_dir, '--outName', out_name, '--odgi', '--chrom-og', '--viz', '--draw', '--haplo', '--lastTrain'] + if not extend: + # collapse self-alignments are not derived from the GAF, so they cannot be reused + cactus_pangenome_cmd += ['--collapse'] + else: + base = self._build_primates_base_graph(binariesMode, mc_config_path, ['simHuman', 'simChimp']) + cactus_pangenome_cmd += ['--extendGFA', base['gfa'], '--extendGAF', base['gaf']] if mgSplit: cactus_pangenome_cmd += ['--mgSplit'] else: @@ -525,6 +551,103 @@ def _run_evolver_primates_pangenome(self, binariesMode, mgSplit = False): wave_vcf_bytes = os.path.getsize(os.path.join(out_dir, out_name + '.simChimp.wave.vcf.gz')) self.assertGreaterEqual(wave_vcf_bytes, 300000) + if extend: + # the genomes added on top of the base graph are in the graph it was extended into + self.assertEqual(self._gfa_genomes(os.path.join(out_dir, out_name + '.sv.gfa.gz')), set(self.PRIMATES)) + + def _build_primates_base_graph(self, binariesMode, config_path, genomes): + """ build the graph and mappings for a subset of the primates with the step-by-step tools, + for a later run to extend. returns the gfa/gaf paths """ + work_dir = os.path.join(self.tempDir, 'extend-base') + os.makedirs(work_dir, exist_ok=True) + seqfile = os.path.join(work_dir, 'primates.txt') + self._write_primates_seqfile(seqfile, genomes) + base = {'gfa': os.path.join(work_dir, 'base.sv.gfa.gz'), + 'gaf': os.path.join(work_dir, 'base.gaf.gz'), + 'paf': os.path.join(work_dir, 'base.paf'), + 'fa': os.path.join(work_dir, 'base.sv.gfa.fa.gz')} + cactus_opts = ['--binariesMode', binariesMode, '--logInfo', '--workDir', self.tempDir, + '--configFile', config_path] + subprocess.check_call(['cactus-minigraph', os.path.join(self.tempDir, 'js-extend-base-mg'), seqfile, base['gfa'], + '--reference', 'simHuman', 'simChimp'] + cactus_opts) + subprocess.check_call(['cactus-graphmap', os.path.join(self.tempDir, 'js-extend-base-gm'), seqfile, base['gfa'], + base['paf'], '--outputFasta', base['fa'], + '--reference', 'simHuman', 'simChimp'] + cactus_opts) + return base + + def _gfa_text(self, gfa_path): + """ the decompressed contents of a (bgzipped) GFA. compared instead of the file itself + because bgzip block boundaries depend on the thread count, so two byte-identical graphs + can have different .gz bytes """ + self.assertTrue(os.path.exists(gfa_path), '{} not found'.format(gfa_path)) + return subprocess.check_output('zcat {}'.format(gfa_path), shell=True) + + def _gfa_genomes(self, gfa_path): + """ the set of genomes an rGFA's SN tags name """ + genomes = set() + for line in self._gfa_text(gfa_path).decode().split('\n'): + if line.startswith('S'): + for tok in line.split('\t')[4:]: + if tok.startswith('SN:Z:'): + genomes.add(tok[5:].split('#')[0]) + break + return genomes + + def _run_primates_extend_steps(self, binariesMode, seqfile_genomes, extend_genomes, config_path, out_prefix): + """ cactus-minigraph + cactus-graphmap on seqfile_genomes, then the same two commands again + on extend_genomes with --extendGFA/--extendGAF pointed at the first run's output. + + returns (base outputs, extended outputs) as dicts of gfa/paf/gaf paths. only the two + stages that --extendGFA/--extendGAF touch are run: split/align/join are unchanged by them + and are covered by the end-to-end test below """ + cactus_opts = ['--binariesMode', binariesMode, '--logInfo', '--workDir', self.tempDir, + '--configFile', config_path] + + def run_stage(tag, genomes, extend_from): + out = {} + work_dir = os.path.join(self.tempDir, '{}-{}'.format(out_prefix, tag)) + os.makedirs(work_dir, exist_ok=True) + # cactus-graphmap edits the seqfile it is given in place (to add _MINIGRAPH_), so each + # stage gets its own copy + seqfile = os.path.join(work_dir, 'primates.txt') + self._write_primates_seqfile(seqfile, genomes) + out['gfa'] = os.path.join(work_dir, 'pg.sv.gfa.gz') + out['paf'] = os.path.join(work_dir, 'pg.paf') + out['gaf'] = os.path.join(work_dir, 'pg.gaf.gz') + out['fa'] = os.path.join(work_dir, 'pg.sv.gfa.fa.gz') + + mg_cmd = ['cactus-minigraph', os.path.join(self.tempDir, 'js-{}-{}-mg'.format(out_prefix, tag)), + seqfile, out['gfa'], '--reference', 'simChimp'] + gm_cmd = ['cactus-graphmap', os.path.join(self.tempDir, 'js-{}-{}-gm'.format(out_prefix, tag)), + seqfile, out['gfa'], out['paf'], '--outputFasta', out['fa'], '--reference', 'simChimp'] + if extend_from: + mg_cmd += ['--extendGFA', extend_from['gfa']] + gm_cmd += ['--extendGAF', extend_from['gaf']] + subprocess.check_call(mg_cmd + cactus_opts) + subprocess.check_call(gm_cmd + cactus_opts) + return out + + base = run_stage('base', seqfile_genomes, None) + extended = run_stage('ext', extend_genomes, base) + extended['cactus_opts'] = cactus_opts + extended['genomes'] = extend_genomes + return base, extended + + def _remap_onto_extended(self, extended, out_prefix): + """ re-run cactus-graphmap on an extended graph with --remap, which maps every genome + instead of reusing any. returns the PAF path """ + work_dir = os.path.join(self.tempDir, '{}-remap'.format(out_prefix)) + os.makedirs(work_dir, exist_ok=True) + seqfile = os.path.join(work_dir, 'primates.txt') + self._write_primates_seqfile(seqfile, extended['genomes']) + paf = os.path.join(work_dir, 'pg.paf') + subprocess.check_call(['cactus-graphmap', os.path.join(self.tempDir, 'js-{}-remap'.format(out_prefix)), + seqfile, extended['gfa'], paf, + '--outputFasta', os.path.join(work_dir, 'pg.sv.gfa.fa.gz'), + '--reference', 'simChimp', '--extendGAF', extended['gaf'], '--remap'] + + extended['cactus_opts']) + return paf + def _run_evolver_primates_step_by_step_mgsplit(self, binariesMode, train=False): """ primates star test but using graphmap pangenome pipeline with chromosome splitting (there are no chromosomes to split, but it still bangs most of the interface) @@ -649,8 +772,14 @@ def _run_yeast_pangenome_step_by_step(self, binariesMode): '--vg'] + vg_files + ['--hal'] + hal_files + ['--xg', '--vcf', '--giraffe', 'clip', 'filter', '--lrGiraffe'] + cactus_opts + ['--indexCores', '4']) - def _run_yeast_pangenome(self, binariesMode, mgSplit=False, wholeGenomeRef=False, collapse=False, gref=None, vcfL=None): - """ yeast pangenome chromosome by chromosome pipeline, as run through a single invocations + def _run_yeast_pangenome(self, binariesMode, mgSplit=False, wholeGenomeRef=False, collapse=False, gref=None, + vcfL=None, extend=None): + """ yeast pangenome chromosome by chromosome pipeline, as run through a single invocations. + + extend, if given, is the list of genomes to build a graph out of first, with the + step-by-step tools, for the run below to extend with --extendGFA. unlike the primates + tests this exercises the translated PAF through cactus-graphmap-split, which is the one + downstream stage that reads it before cactus-align """ orig_seq_file_path = './examples/yeastPangenome.txt' @@ -677,12 +806,34 @@ def _run_yeast_pangenome(self, binariesMode, mgSplit=False, wholeGenomeRef=False cactus_pangenome_cmd += ['--gref', gref] if vcfL is not None: cactus_pangenome_cmd += ['--vcfL', str(vcfL)] + if extend: + base = self._build_yeast_base_graph(binariesMode, orig_seq_file_path, extend) + cactus_pangenome_cmd += ['--extendGFA', base['gfa'], '--extendGAF', base['gaf']] subprocess.check_call(cactus_pangenome_cmd + cactus_opts) #compatibility with older test subprocess.check_call(['mkdir', '-p', os.path.join(self.tempDir, 'chroms')]) subprocess.check_call(['mv', os.path.join(join_path, 'chrom-subproblems', 'contig_sizes.tsv'), os.path.join(self.tempDir, 'chroms')]) + def _build_yeast_base_graph(self, binariesMode, seqfile_path, genomes): + """ build the graph and mappings for a subset of the yeast strains with the step-by-step + tools, for the run above to extend. returns the gfa/gaf paths """ + work_dir = os.path.join(self.tempDir, 'extend-base') + os.makedirs(work_dir, exist_ok=True) + seqfile = os.path.join(work_dir, 'yeast.txt') + self._write_seqfile_subset(seqfile_path, seqfile, genomes) + base = {'gfa': os.path.join(work_dir, 'base.sv.gfa.gz'), + 'gaf': os.path.join(work_dir, 'base.gaf.gz'), + 'paf': os.path.join(work_dir, 'base.paf'), + 'fa': os.path.join(work_dir, 'base.sv.gfa.fa.gz')} + cactus_opts = ['--binariesMode', binariesMode, '--logInfo', '--workDir', self.tempDir, '--maxCores', '4'] + subprocess.check_call(['cactus-minigraph', os.path.join(self.tempDir, 'js-yeast-base-mg'), seqfile, base['gfa'], + '--reference', 'S288C', 'DBVPG6044'] + cactus_opts) + subprocess.check_call(['cactus-graphmap', os.path.join(self.tempDir, 'js-yeast-base-gm'), seqfile, base['gfa'], + base['paf'], '--outputFasta', base['fa'], + '--reference', 'S288C', 'DBVPG6044'] + cactus_opts) + return base + def _validate_sv_gfa(self, gfa_path): """ run `zcat | vg validate -` and assert it passes. catches things like missing edges in the merged minigraph SV GFA. @@ -1697,6 +1848,121 @@ def testEvolverPrimatesPangenomeSplitLocal(self): # todo: tune config so that delta can be reduced self._check_maf_accuracy(self._out_hal("local"), delta=(0.025,0.025), dataset='primates') + def testPangenomeExtendNullLocal(self): + """ Extending a pangenome by nothing must reproduce it exactly. + + This is the load-bearing test for reusing mappings: minigraph GAF is in stable coordinates, + so a genome's existing mappings can be re-derived against the extended graph rather than + recomputed. Feeding the graph back unchanged makes "re-derived" and "as originally mapped" + the same thing, so anything the round trip through PanSN naming, the per-genome GAF split, + or gaf2unstable/gaf2paf loses shows up here as a difference. """ + config_path = self._write_pangenome_config('extend-null', + graphmap_attribs={"minigraphConstructBatchSize": "2"}) + genomes = ['simHuman', 'simChimp'] + base, extended = self._run_primates_extend_steps('local', genomes, genomes, config_path, 'null') + + # no genomes were added, so the graph is the input passed straight through + self.assertEqual(self._gfa_text(extended['gfa']), self._gfa_text(base['gfa'])) + # the mappings were re-derived from the published GAF rather than remapped, and came out + # byte for byte the same as minigraph's + with open(base['paf'], 'rb') as base_paf, open(extended['paf'], 'rb') as ext_paf: + self.assertEqual(ext_paf.read(), base_paf.read()) + # and the republished GAF is intact, so the extended run can itself be extended + self.assertEqual(subprocess.check_output('zcat {}'.format(extended['gaf']), shell=True), + subprocess.check_output('zcat {}'.format(base['gaf']), shell=True)) + + # an unchanged graph is the one graph construction hands back without writing it, so it is + # the one that can come out at the compression of the file it was read from rather than of + # the file it is being written to. everything downstream reads the output suffix + plain_dir = os.path.join(self.tempDir, 'null-plain') + os.makedirs(plain_dir, exist_ok=True) + plain_seqfile = os.path.join(plain_dir, 'primates.txt') + self._write_primates_seqfile(plain_seqfile, genomes) + plain_gfa = os.path.join(plain_dir, 'pg.sv.gfa') + subprocess.check_call(['cactus-minigraph', os.path.join(self.tempDir, 'js-null-plain'), + plain_seqfile, plain_gfa, '--reference', 'simChimp', + '--extendGFA', base['gfa'], + '--binariesMode', 'local', '--logInfo', '--workDir', self.tempDir, + '--configFile', config_path]) + with open(plain_gfa, 'rb') as plain_file: + plain_bytes = plain_file.read() + self.assertFalse(plain_bytes.startswith(b'\x1f\x8b'), + 'a .gfa output path got gzip bytes: the seed graph kept its input compression') + self.assertEqual(plain_bytes, self._gfa_text(base['gfa'])) + + def testPangenomeExtendConstructionLocal(self): + """ Extending a graph must build the same graph as constructing it in one go. + + minigraph construction is sequential -- each genome is added to the graph built from the + ones before it -- which is the whole reason extending is possible. With the input order + pinned (minigraphSortInput=none), adding simGorilla and simOrang to a simChimp+simHuman + graph is the same command sequence cactus already runs when it batches construction, so the + two graphs should be identical. If this ever stops holding, the cost model behind + --extendGFA needs revisiting. """ + config_path = self._write_pangenome_config('extend-construct', + graphmap_attribs={"minigraphConstructBatchSize": "2", + "minigraphSortInput": "none"}) + base, extended = self._run_primates_extend_steps('local', ['simHuman', 'simChimp'], self.PRIMATES, + config_path, 'construct') + + # the same four genomes, built in the same order, in one run + scratch_dir = os.path.join(self.tempDir, 'construct-scratch') + os.makedirs(scratch_dir, exist_ok=True) + scratch_seqfile = os.path.join(scratch_dir, 'primates.txt') + self._write_primates_seqfile(scratch_seqfile, self.PRIMATES) + scratch_gfa = os.path.join(scratch_dir, 'pg.sv.gfa.gz') + subprocess.check_call(['cactus-minigraph', os.path.join(self.tempDir, 'js-construct-scratch'), + scratch_seqfile, scratch_gfa, '--reference', 'simChimp', + '--binariesMode', 'local', '--logInfo', '--workDir', self.tempDir, + '--configFile', config_path]) + + self.assertEqual(self._gfa_genomes(extended['gfa']), set(self.PRIMATES)) + self.assertEqual(self._gfa_text(extended['gfa']), self._gfa_text(scratch_gfa)) + self._validate_sv_gfa(extended['gfa']) + + # the extended graph is strictly bigger than the one it came from, and every genome has + # mappings in the PAF + self.assertGreater(len(self._gfa_text(extended['gfa'])), len(self._gfa_text(base['gfa']))) + with open(extended['paf']) as paf_file: + paf_queries = set(line.split('\t')[0].split('|')[0] for line in paf_file) + self.assertEqual(paf_queries, set('id=' + genome for genome in self.PRIMATES)) + + # the reused mappings put every alignment at the same place a fresh mapping does. the + # graph is identical either way, so the only thing being compared is what reuse did to the + # coordinates -- the nodes an existing genome was mapped to have been split underneath it + scratch_paf = os.path.join(scratch_dir, 'pg.paf') + subprocess.check_call(['cactus-graphmap', os.path.join(self.tempDir, 'js-construct-scratch-gm'), + scratch_seqfile, scratch_gfa, scratch_paf, + '--outputFasta', os.path.join(scratch_dir, 'pg.sv.gfa.fa.gz'), + '--reference', 'simChimp'] + extended['cactus_opts']) + self.assertEqual(self._paf_coords(extended['paf']), self._paf_coords(scratch_paf)) + + # and --remap, which maps everything instead of reusing anything, gives back exactly the + # from-scratch PAF: same graph, same mappings, nothing carried over + self.assertEqual(self._paf_bytes(self._remap_onto_extended(extended, 'construct')), + self._paf_bytes(scratch_paf)) + + def _paf_coords(self, paf_path): + """ the set of (query, target, coordinate) tuples in a PAF, ignoring the tags. two PAFs + that agree here place every alignment identically; only the cigars can still differ, and + those differ by indel placement in ambiguous runs """ + with open(paf_path) as paf_file: + return set(tuple(line.rstrip('\n').split('\t')[:12]) for line in paf_file) + + def _paf_bytes(self, paf_path): + with open(paf_path, 'rb') as paf_file: + return paf_file.read() + + def testEvolverPrimatesPangenomeExtendLocal(self): + """ Evolver (star) primates built two genomes at a time with cactus-pangenome --extendGFA, + checked against the same accuracy baseline as the from-scratch run """ + name = "local" + self._run_evolver_primates_pangenome(name, extend=True) + + # check the output + # todo: tune config so that delta can be reduced + self._check_maf_accuracy(self._out_hal("local"), delta=(0.025,0.025), dataset='primates') + def testEvolverPrimatesPangenomeStepByStepSplitLocal(self): """ Evolver primates but using the step-by-step cactus-pangenome interface with splitting """ @@ -1734,6 +2000,20 @@ def testYeastPangenomeLocal(self): # check the output self._check_yeast_pangenome(name, other_ref='DBVPG6044', expect_odgi=True, expect_haplo=False, expect_unchopped_gfa=True) + def testYeastPangenomeExtendLocal(self): + """ Yeast pangenome built three strains at a time with cactus-pangenome --extendGFA. + + Unlike the primates extend tests, this one splits by chromosome, so the reused mappings go + through cactus-graphmap-split -- the one stage downstream of cactus-graphmap that reads the + PAF before cactus-align does. It is checked against exactly the same pinned graph + statistics as the from-scratch yeast run, because the graph covers the same six strains. """ + name = "local" + self._run_yeast_pangenome(name, extend=['S288C', 'DBVPG6044', 'SK1']) + + # check the output + self._check_yeast_pangenome(name, other_ref='DBVPG6044', expect_odgi=True, expect_haplo=False, + expect_unchopped_gfa=True) + YEAST_URL = 'https://github.com/ComparativeGenomicsToolkit/cactusTestData/raw/master/yeast/{}.fa.gz' PANPATCH_TARGET = 'SK1' PANPATCH_DONOR = 'DBVPG6044' From 8c93c4638dc9ec55ab0b8c1e3529872b46a54304 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Tue, 15 Sep 2026 08:57:30 -0400 Subject: [PATCH 2/4] Check reused mappings against the graph before fanning out Resolving a reused GAF is the same work for every genome, so a GAF that does not belong to the graph fails identically in all of them -- once per genome, after the fan-out, and again on every Toil retry. On an HPRC-scale cluster run that is hours of failures for one wrong input. Resolve a sample up front instead, in one job, and say what is wrong. The pair has to come from the same graphmap run, and there is one way to get a mismatched pair that looks like a matching one. A --mgSplit run publishes an .sv.gfa.gz and an .gaf.gz that are not a pair: the GAF is the whole-genome first pass, mapped against the reference-only graph cactus-minigraph --refOnly builds, while the GFA is the merged per-chromosome graphs. Same stable coordinates, different node boundaries, so gaf2unstable fails its tiling assert rather than on a name, and nothing in that error points at the inputs. Found on a 30-way CHM13 graph extended to 36. The construction was fine -- every pre-existing genome came through with byte-identical sequence and the six new samples were added -- and every reused genome then failed in gaf2unstable. The tell is that the GAF's paths only ever walked through CHM13 and GRCh38 while the graph it was being resolved against held all 30 genomes, so the error now reports the genomes the sampled records actually traverse. Such a pangenome can still be extended with --extendGFA alone, which is the saving that matters: construction is the expensive half. --- doc/pangenome.md | 14 ++++++ src/cactus/refmap/cactus_graphmap.py | 65 ++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+) diff --git a/doc/pangenome.md b/doc/pangenome.md index a8701818f..cd68f3c69 100644 --- a/doc/pangenome.md +++ b/doc/pangenome.md @@ -316,6 +316,20 @@ of `minigraph` into minutes of file conversion. `--extendGAF` is optional. Without it — extending a published release for which only the GFA is available, say — every genome is mapped again, which is the fallback and costs a full mapping stage. +**The GFA and the GAF must come from the same graphmap run.** They describe the same node +boundaries, and reusing mappings against a graph that has different ones fails. There is one way to +get a mismatched pair that looks like a matching one: a **`--mgSplit` run publishes an +`.sv.gfa.gz` and an `.gaf.gz` that are not a pair.** Its GAF is the whole-genome +first pass, mapped against the *reference-only* graph `cactus-minigraph --refOnly` builds, while its +GFA is the merged per-chromosome graphs. The giveaway is that such a GAF's paths only ever walk +through the reference genomes, however many samples the graph holds. + +So a pangenome built with `--mgSplit` can be extended with `--extendGFA`, but not with +`--extendGAF`. That costs the mapping stage and keeps the saving that matters, since construction +is the expensive half. A sample of the reused mappings is resolved against the graph before +anything else runs, so a mismatched pair fails in one job with a diagnosis rather than in every +per-genome job after the fan-out. + #### Why an extended pangenome is not identical to one built all at once There are exactly two sources of difference, and it is worth being precise about which is which. diff --git a/src/cactus/refmap/cactus_graphmap.py b/src/cactus/refmap/cactus_graphmap.py index 3cf2c1784..8abfe6b87 100644 --- a/src/cactus/refmap/cactus_graphmap.py +++ b/src/cactus/refmap/cactus_graphmap.py @@ -412,6 +412,17 @@ def minigraph_workflow(job, options, config, seq_id_map, gfa_id, graph_event, sa gfa_id = gfa_unzip_job.rv() gfa_id_size *= 10 options.minigraphGFA = options.minigraphGFA[:-3] + if extend_gaf_id: + # resolving a reused GAF is the same work for every genome, so a GAF that does not belong + # to this graph fails identically in all of them -- once per genome, after the fan-out, and + # again on every Toil retry. Checking a sample up front turns that into one quick failure + # with something actionable in it. chained onto the unzip (when there is one) because it + # needs the same uncompressed graph the per-genome jobs use + check_parent = gfa_unzip_job if zipped_gfa else root_job + check_parent.addFollowOnJobFn(check_reusable_gaf, config, extend_gaf_id, gfa_id, genome_names, + options.extendGAF, options.minigraphGFA, + disk=4*gfa_id_size, memory=cactus_clamp_memory(2*gfa_id_size)) + paf_job = Job.wrapJobFn(minigraph_map_all, options, config, gfa_id, seq_id_map, graph_event, extend_gaf_map) root_job.addFollowOn(paf_job) @@ -735,6 +746,60 @@ def minigraph_map_one(job, config, event_name, fa_file_id, gfa_file_id): return stable_gaf_to_paf(job, config, gaf_path, gfa_path) +# how many reused GAF records the up-front check resolves before trusting the rest +GAF_REUSE_CHECK_RECORDS = 1000 + +def check_reusable_gaf(job, config, gaf_file_id, gfa_file_id, genome_names, gaf_path, gfa_path): + """ Resolve the first few reused mappings against the graph, and fail with a diagnosis if they + do not fit. + + The pair has to come from the same graphmap run. The trap is that a --mgSplit run publishes + .sv.gfa.gz and .gaf.gz that look like a pair and are not: the GAF is the + whole-genome pass against the reference-only first-pass graph, while the GFA is the merged + per-chromosome graphs. Same stable coordinates, different node boundaries, so gaf2unstable + fails on the tiling rather than on a name and the mismatch is not obvious from the error. """ + work_dir = job.fileStore.getLocalTempDir() + gfa_local = os.path.join(work_dir, 'mg.gfa') + job.fileStore.readGlobalFile(gfa_file_id, gfa_local) + pansn_head = os.path.join(work_dir, 'check.pansn.gaf') + head = os.path.join(work_dir, 'check.gaf') + job.fileStore.readGlobalFile(gaf_file_id, pansn_head + '.full') + with open(pansn_head, 'w') as out_file: + opener = gzip.open if gaf_path.endswith('.gz') else open + with opener(pansn_head + '.full', 'rt') as in_file: + for i, line in enumerate(in_file): + if i >= GAF_REUSE_CHECK_RECORDS: break + out_file.write(line) + os.remove(pansn_head + '.full') + gaf_from_pansn(genome_names, pansn_head, head) + + try: + cactus_call(parameters=['gaf2unstable', head, '-g', gfa_local, '-o', os.path.join(work_dir, 'lens.tsv')], + outfile=os.path.join(work_dir, 'check.unstable.gaf'), job_memory=job.memory) + except RuntimeError as e: + # name the genomes each side actually talks about: the mismatch above shows up as a GAF + # whose paths only ever name the references while the graph is full of samples + step_genomes = set() + with open(head) as in_file: + for line in in_file: + toks = line.split('\t') + if len(toks) > 5: + for step in gaf_step_re.findall(toks[5]): + name = step[1:].rsplit(':', 1)[0] if ':' in step else step[1:] + step_genomes.add(name[3:name.find('|')] if name.startswith('id=') and '|' in name else name) + raise RuntimeError( + 'The mappings in {} do not resolve against {}: the two must come from the same graphmap run.\n' + 'The first {} records only ever walk through: {}.\n' + 'A --mgSplit run is the usual way to get a mismatched pair that looks like a matching one: its ' + '.gaf.gz is the whole-genome pass against the reference-only first-pass graph, while its ' + '.sv.gfa.gz is the merged per-chromosome graphs. Drop --extendGAF to map every genome ' + 'against the extended graph instead -- --extendGFA still saves the construction, which is the ' + 'expensive half.\nUnderlying error: {}'.format( + gaf_path, gfa_path, GAF_REUSE_CHECK_RECORDS, + ' '.join(sorted(step_genomes)[:12]) or '(nothing)', e)) + + RealtimeLogger.info('Reused mappings from {} resolve against the graph'.format(gaf_path)) + def translate_gaf_one(job, config, event_name, gaf_file_id, gfa_file_id, genome_names): """ Re-derive one genome's PAF from mappings it already has, against a (possibly extended) graph. From 3e826b166a594cd3541a3d4db3d99ee0551fdba8 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Tue, 15 Sep 2026 09:28:11 -0400 Subject: [PATCH 3/4] Publish the --mgSplit first pass under its own prefix A --mgSplit run builds twice. The whole-genome pass runs cactus-minigraph --refOnly, so its graph holds the references alone and the mappings taken against it walk through nothing else; the real per-chromosome graphs and mappings land in chrom-minigraph/ and chrom-graphmap/. At the end of the run cactus-graphmap-join merges the per-chromosome graphs and export_join_data writes the result under the key "sv.gfa.gz", which resolves to the same .sv.gfa.gz the first pass exported twelve hours earlier. So the first pass's graph was silently replaced, while its mappings, PAF and node fasta stayed behind under the plain name next to a graph none of them describe. Nothing said so, and nothing checked: on a 30-way CHM13 run the timestamps were the only evidence, .gaf.gz/.paf/.sv.gfa.fa.gz from 22:03 and .sv.gfa.gz from 10:39 the next morning. The first pass now publishes as .refonly.sv.gfa.gz, .refonly.gaf.gz, .refonly.paf and .refonly.sv.gfa.fa.gz. Nothing collides, the names say what they are, and the two things that are genuinely a pair stay one: the refonly graph and the refonly GAF can be extended together, while the merged .sv.gfa.gz has no GAF to go with it and takes --extendGFA alone. Only the published filenames move, and only with --mgSplit; every stage after the export passes file ids, and a run without --mgSplit is unchanged. The first-pass .train is left alone: it is a scoring matrix trained on fastas, with no graph to be inconsistent with. Found because --extendGAF was handed the plain pair from an mgSplit run and gaf2unstable failed its tiling assert on every reused genome -- the only consumer that happened to notice, and only because it asserts. --- doc/pangenome.md | 29 ++++++++++++++++----------- src/cactus/refmap/cactus_pangenome.py | 13 +++++++++--- test/evolverTest.py | 14 +++++++++++++ 3 files changed, 41 insertions(+), 15 deletions(-) diff --git a/doc/pangenome.md b/doc/pangenome.md index cd68f3c69..8eb7b47e0 100644 --- a/doc/pangenome.md +++ b/doc/pangenome.md @@ -317,18 +317,23 @@ of `minigraph` into minutes of file conversion. available, say — every genome is mapped again, which is the fallback and costs a full mapping stage. **The GFA and the GAF must come from the same graphmap run.** They describe the same node -boundaries, and reusing mappings against a graph that has different ones fails. There is one way to -get a mismatched pair that looks like a matching one: a **`--mgSplit` run publishes an -`.sv.gfa.gz` and an `.gaf.gz` that are not a pair.** Its GAF is the whole-genome -first pass, mapped against the *reference-only* graph `cactus-minigraph --refOnly` builds, while its -GFA is the merged per-chromosome graphs. The giveaway is that such a GAF's paths only ever walk -through the reference genomes, however many samples the graph holds. - -So a pangenome built with `--mgSplit` can be extended with `--extendGFA`, but not with -`--extendGAF`. That costs the mapping stage and keeps the saving that matters, since construction -is the expensive half. A sample of the reused mappings is resolved against the graph before -anything else runs, so a mismatched pair fails in one job with a diagnosis rather than in every -per-genome job after the fan-out. +boundaries, and reusing mappings against a graph that has different ones fails. + +With `--mgSplit` that means the `.refonly.` files, not the plain ones. Such a run has two +construction passes: a whole-genome one that `cactus-minigraph --refOnly` builds from the +references alone, and the per-chromosome one whose graphs `cactus-graphmap-join` merges into +`.sv.gfa.gz` at the end of the run. The first pass therefore publishes its graph, mappings +and node fasta under a `.refonly.` prefix — `refonly.sv.gfa.gz`, `refonly.gaf.gz`, +`refonly.paf`, `refonly.sv.gfa.fa.gz` — so that they are not mistaken for, or overwritten by, the +merged graph they do not describe. `.refonly.sv.gfa.gz` and `.refonly.gaf.gz` are +a genuine pair and can be extended together; the merged `.sv.gfa.gz` has no GAF to go with +it, so extend it with `--extendGFA` alone. That still keeps the saving that matters, since +construction is the expensive half. + +A sample of the reused mappings is resolved against the graph before anything else runs, so a +mismatched pair fails in one job with a diagnosis rather than in every per-genome job after the +fan-out. The tell is a GAF whose paths only ever walk through the reference genomes, however many +samples the graph holds. #### Why an extended pangenome is not identical to one built all at once diff --git a/src/cactus/refmap/cactus_pangenome.py b/src/cactus/refmap/cactus_pangenome.py index c62476a6f..f633d5c16 100644 --- a/src/cactus/refmap/cactus_pangenome.py +++ b/src/cactus/refmap/cactus_pangenome.py @@ -637,7 +637,14 @@ def pangenome_end_to_end_workflow(job, options, config_wrapper, seq_id_map, seq_ assert type(options.reference) == list # cactus_minigraph - sv_gfa_path = os.path.join(options.outDir, options.outName + '.sv.gfa.gz') + # with --mgSplit the whole-genome pass is a means to an end: cactus-minigraph runs --refOnly, so + # its graph holds the references alone and the mappings below are against that graph. The real + # per-chromosome graphs and mappings land in chrom-minigraph/ and chrom-graphmap/, and at the end + # of the run cactus-graphmap-join merges the former back onto .sv.gfa.gz. Published + # under the plain name these first-pass files would be silently replaced by, or left beside, a + # graph they do not describe -- so they get a prefix that says what they are + first_pass_name = options.outName + ('.refonly' if options.mgSplit else '') + sv_gfa_path = os.path.join(options.outDir, first_pass_name + '.sv.gfa.gz') options.batch = False options.refOnly = False @@ -666,8 +673,8 @@ def pangenome_end_to_end_workflow(job, options, config_wrapper, seq_id_map, seq_ minigraph_wrapper_job = minigraph_job.addFollowOnJobFn(export_minigraph_wrapper, options, pansn_sv_gfa_id, sv_gfa_path, last_scores_id) # cactus_graphmap - paf_path = os.path.join(options.outDir, options.outName + '.paf') - gfa_fa_path = os.path.join(options.outDir, options.outName + '.sv.gfa.fa.gz') + paf_path = os.path.join(options.outDir, first_pass_name + '.paf') + gfa_fa_path = os.path.join(options.outDir, first_pass_name + '.sv.gfa.fa.gz') options.minigraphGFA = sv_gfa_path options.outputFasta = gfa_fa_path graph_event = getOptionalAttrib(findRequiredNode(config_node, "graphmap"), "assemblyName", default="_MINIGRAPH_") diff --git a/test/evolverTest.py b/test/evolverTest.py index b825337e8..9299fae1b 100644 --- a/test/evolverTest.py +++ b/test/evolverTest.py @@ -555,6 +555,20 @@ def _run_evolver_primates_pangenome(self, binariesMode, mgSplit = False, extend # the genomes added on top of the base graph are in the graph it was extended into self.assertEqual(self._gfa_genomes(os.path.join(out_dir, out_name + '.sv.gfa.gz')), set(self.PRIMATES)) + if mgSplit: + # the whole-genome first pass is reference-only, and cactus-graphmap-join later writes + # the merged per-chromosome graph to .sv.gfa.gz. Both have to survive under + # names that say which is which, or the first pass's mappings are left beside a graph + # they do not describe -- which is exactly what breaks --extendGAF + refonly_gfa = os.path.join(out_dir, out_name + '.refonly.sv.gfa.gz') + self.assertTrue(os.path.exists(refonly_gfa), 'no {} from the --mgSplit first pass'.format(refonly_gfa)) + self.assertEqual(self._gfa_genomes(refonly_gfa), {'simHuman', 'simChimp'}) + self.assertEqual(self._gfa_genomes(os.path.join(out_dir, out_name + '.sv.gfa.gz')), set(self.PRIMATES)) + # the first pass's other artifacts travel with it rather than sitting under the plain name + for ext in ['.paf', '.gaf.gz', '.sv.gfa.fa.gz']: + self.assertTrue(os.path.exists(os.path.join(out_dir, out_name + '.refonly' + ext)), + 'first-pass {} not published under the .refonly prefix'.format(ext)) + def _build_primates_base_graph(self, binariesMode, config_path, genomes): """ build the graph and mappings for a subset of the primates with the step-by-step tools, for a later run to extend. returns the gfa/gaf paths """ From 11c680186a67bedc060f9e4417f5afe5253e23d4 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Tue, 15 Sep 2026 09:57:32 -0400 Subject: [PATCH 4/4] Generalize --extendGFA/--extendGAF into --inGFA/--inGAF Starting the pipeline from a graph somebody already built and adding genomes to one are the same operation with a different seqfile, so they get one pair of options rather than a flag that says "extend". What happens depends only on how the seqfile's genomes compare to the graph's: names genomes the graph lacks construct those in, leave the rest alone names exactly the graph's construct nothing, resume from the graph omits one the graph has error, as before The resume half of that already worked -- it is the no-new-genomes path, which the null-extend test has been covering since the code review turned up the compression bug in it -- but nothing said so and the options were named for the other half. It saves driving cactus-graphmap-split, cactus-align --batch and cactus-graphmap-join by hand: --inGFA alone picks up at cactus-graphmap, and --inGFA with --inGAF picks up at cactus-graphmap-split. Two things were missing rather than renamed. Nothing ran cactus-pangenome end to end with nothing to construct, which is exactly the case being named here, so testEvolverPrimatesPangenomeResumeLocal does: a graph of all four primates, checked to come back out unchanged and to reach the same MAF accuracy baseline as a from-scratch run. --lastTrain quietly produced no scoring model when there was nothing to construct, because the early return handed back train_id=None and cactus-align fell through to the defaults. last_train reads fastas, not the graph, so it now runs on the resume path too, over the full genome set rather than the map that has been narrowed to what minigraph still has to be given. The old spellings are gone rather than aliased: nothing has shipped, and one name per thing is worth more than the compatibility nobody needs yet. --- Makefile | 3 + doc/pangenome.md | 69 ++++++++++++++++------- src/cactus/refmap/cactus_graphmap.py | 62 ++++++++++---------- src/cactus/refmap/cactus_minigraph.py | 81 ++++++++++++++++----------- src/cactus/refmap/cactus_pangenome.py | 65 ++++++++++----------- src/cactus/refmap/cactus_panpatch.py | 4 +- test/evolverTest.py | 60 +++++++++++++------- 7 files changed, 207 insertions(+), 137 deletions(-) diff --git a/Makefile b/Makefile index 8fc779e66..8993e78eb 100644 --- a/Makefile +++ b/Makefile @@ -240,6 +240,9 @@ pangenome_extend_construction_test_local: yeast_test_extend_local: PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testYeastPangenomeExtendLocal +evolver_test_primates_pangenome_resume_local: all ${CWD}/test/primates-truth.maf + PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testEvolverPrimatesPangenomeResumeLocal + evolver_test_primates_pangenome_extend_local: all ${CWD}/test/primates-truth.maf PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testEvolverPrimatesPangenomeExtendLocal diff --git a/doc/pangenome.md b/doc/pangenome.md index 8eb7b47e0..69c5f4c12 100644 --- a/doc/pangenome.md +++ b/doc/pangenome.md @@ -267,40 +267,64 @@ For `--vgFilter`, the filter threshold is inferred from the `.dX.vg` filename pa Note: per-chromosome output options (`--chrom-vg`, `--chrom-og`, `--viz`, `--draw`) cannot be used with bypass options, as you already have those files from the previous run. Also, bypass options are not compatible with graphs that were originally built with `--collapse`. -### Adding Genomes to an Existing Pangenome +### Starting From an Existing Graph or Mappings -`cactus-pangenome --extendGFA` adds genomes to a pangenome you have already built, instead of -rebuilding it from scratch: +`cactus-pangenome --inGFA` starts the pipeline from a minigraph GFA you already have instead of +building one, and `--inGAF` does the same for that graph's mappings. What it does with them depends +only on how the seqfile's genomes compare to the graph's: + +| seqfile vs. graph | what happens | +|---|---| +| names genomes the graph lacks | those are constructed into it; everything already there is left alone | +| names exactly the graph's genomes | nothing is constructed — the run resumes from the graph | +| omits a genome the graph has | error: genomes cannot be removed from a minigraph | + +So the same two options cover adding genomes to a finished pangenome and re-running the back half of +the pipeline on one, without having to drive `cactus-graphmap-split`, `cactus-align --batch` and +`cactus-graphmap-join` by hand. + +**Adding genomes.** `seqfile.txt` lists every genome, the ones already in the graph as well as the +new ones: ``` cactus-pangenome ./js ./seqfile.txt --reference GRCh38 --outDir pg2 --outName pg \ - --extendGFA pg1/pg.sv.gfa.gz --extendGAF pg1/pg.gaf.gz + --inGFA pg1/pg.sv.gfa.gz --inGAF pg1/pg.gaf.gz ``` -`seqfile.txt` lists **every** genome, the ones already in the graph as well as the ones being -added. Cactus works out which are which; the genomes that are already there are left exactly where -they are in the graph, and only the new ones are constructed in. A genome cannot be *removed* from -a minigraph, so leaving one out of the seqfile is an error rather than a way to drop it. - This matters because `minigraph` construction is iterative in the input genomes and dominates the wall time of a large run — weeks, for an HPRC-scale release. Adding 50 genomes to a graph of 450 costs 50 genomes' worth of construction, not 500. +**Resuming.** Give the same command a seqfile whose genomes the graph already holds, and there is +nothing to construct and nothing to map: + +``` +# resume at cactus-graphmap: re-map everything against the graph, then split, align, join +cactus-pangenome ./js ./seqfile.txt --reference GRCh38 --outDir pg2 --outName pg \ + --inGFA pg1/pg.sv.gfa.gz + +# resume at cactus-graphmap-split: re-derive the PAF from the mappings, then split, align, join +cactus-pangenome ./js ./seqfile.txt --reference GRCh38 --outDir pg2 --outName pg \ + --inGFA pg1/pg.sv.gfa.gz --inGAF pg1/pg.gaf.gz +``` + +The seqfile still needs every genome's FASTA either way, since `cactus-align` reads them. + The step-by-step interface has the same two options, and they interoperate with the one-shot one in both directions: ``` cactus-minigraph ./js ./seqfile.txt pg2/pg.sv.gfa.gz --reference GRCh38 \ - --extendGFA pg1/pg.sv.gfa.gz + --inGFA pg1/pg.sv.gfa.gz cactus-graphmap ./js ./seqfile.txt pg2/pg.sv.gfa.gz pg2/pg.paf --reference GRCh38 \ - --outputFasta pg2/pg.sv.gfa.fa.gz --extendGAF pg1/pg.gaf.gz + --outputFasta pg2/pg.sv.gfa.fa.gz --inGAF pg1/pg.gaf.gz ``` Everything after these two stages — `cactus-graphmap-split`, `cactus-align`, and especially -`cactus-graphmap-join`'s `vg` indexing — is recomputed in full, and the output of an extended run -is an ordinary pangenome that can itself be extended again. +`cactus-graphmap-join`'s `vg` indexing — is recomputed in full, and the output is an ordinary +pangenome that can itself be used as an `--inGFA`/`--inGAF` input again. -#### What `--extendGAF` does +#### What `--inGAF` does `cactus-graphmap` runs `minigraph` without `--vc`, so the GAF it publishes as `.gaf.gz` is in *stable* coordinates: rGFA `SN`/`SO` sequence names and offsets. Adding genomes to a graph never @@ -310,10 +334,10 @@ node covers stays exactly where it was. That means an existing genome's mappings do not have to be recomputed against the extended graph — they can be re-derived from the published GAF by the same `gaf2unstable` / `gaf2paf` conversion that produced the PAF in the first place, which resolves the stable coordinates into the new, finer node -ids for free. `--extendGAF` is the option that does this, and it turns the mapping stage from hours +ids for free. `--inGAF` is the option that does this, and it turns the mapping stage from hours of `minigraph` into minutes of file conversion. -`--extendGAF` is optional. Without it — extending a published release for which only the GFA is +`--inGAF` is optional. Without it — extending a published release for which only the GFA is available, say — every genome is mapped again, which is the fallback and costs a full mapping stage. **The GFA and the GAF must come from the same graphmap run.** They describe the same node @@ -327,7 +351,7 @@ and node fasta under a `.refonly.` prefix — `refonly.sv.gfa.gz`, `ref `refonly.paf`, `refonly.sv.gfa.fa.gz` — so that they are not mistaken for, or overwritten by, the merged graph they do not describe. `.refonly.sv.gfa.gz` and `.refonly.gaf.gz` are a genuine pair and can be extended together; the merged `.sv.gfa.gz` has no GAF to go with -it, so extend it with `--extendGFA` alone. That still keeps the saving that matters, since +it, so extend it with `--inGFA` alone. That still keeps the saving that matters, since construction is the expensive half. A sample of the reused mappings is resolved against the graph before anything else runs, so a @@ -337,7 +361,12 @@ samples the graph holds. #### Why an extended pangenome is not identical to one built all at once -There are exactly two sources of difference, and it is worth being precise about which is which. +This section is about *adding* genomes. Resuming from a graph that needs nothing constructed does +not enter into it: the graph is handed back unchanged, and with `--remap` (or without `--inGAF`) +the mappings are made against it exactly as a from-scratch run would make them. + +For adding, there are exactly two sources of difference, and it is worth being precise about which +is which. **1. Construction order.** `minigraph` construction is iterative in the input genomes, so the order they go in decides the graph. Building `A B C D` in one go sorts all four by mash distance to the @@ -379,10 +408,10 @@ embarrassingly parallel. * The graph must have been built with a compatible configuration (`minigraphConstructOptions`, the `` `assemblyName`) and with the same `--reference`. -* `--mgSplit` and `--collapse` are not supported with `--extendGFA`. `--mgSplit` has per-chromosome +* `--mgSplit` and `--collapse` are not supported with `--inGFA`. `--mgSplit` has per-chromosome graphs and mappings that would need extending as well; `--collapse` self-alignments come from `minimap2` rather than from the GAF, so there is nothing in the GAF to reuse. -* Standalone `cactus-graphmap --extendGAF` still imports and sanitizes every genome's FASTA even +* Standalone `cactus-graphmap --inGAF` still imports and sanitizes every genome's FASTA even though the reused ones are not mapped. On the `cactus-pangenome` path that work is not wasted — `cactus-align` needs those FASTAs anyway. * A reused GAF needs one adjustment beyond re-running the conversion, handled by diff --git a/src/cactus/refmap/cactus_graphmap.py b/src/cactus/refmap/cactus_graphmap.py index 8abfe6b87..dc5718c8b 100644 --- a/src/cactus/refmap/cactus_graphmap.py +++ b/src/cactus/refmap/cactus_graphmap.py @@ -58,13 +58,13 @@ def main(): parser.add_argument("--mapCores", type=int, help = "Number of cores for minigraph. Overrides graphmap cpu in configuration") parser.add_argument("--collapse", help = "Incorporate minimap2 self-alignments.", action='store_true', default=False) parser.add_argument("--collapseRefPAF", help ="Incorporate given (reference-only) self-alignments in PAF format [Experimental]") - parser.add_argument("--extendGAF", type=str, default=None, + parser.add_argument("--inGAF", type=str, default=None, help = "Reuse the mappings in this GAF (as published by a previous cactus-graphmap or cactus-pangenome run) " "instead of re-running minigraph for the genomes it covers. Minigraph GAF is in stable coordinates, which " - "node splitting does not change, so these mappings are simply re-derived against the given (extended) graph. " - "Only genomes that are not in the GAF are mapped. Intended for use with cactus-minigraph --extendGFA") + "node splitting does not change, so these mappings are re-derived against the given graph instead. Only " + "genomes the GAF does not cover are mapped. Intended for use with cactus-minigraph --inGFA") parser.add_argument("--remap", action="store_true", default=False, - help = "Map every genome with minigraph even if --extendGAF already covers it. Slower, but the existing " + help = "Map every genome with minigraph even if --inGAF already covers it. Slower, but the existing " "genomes then see the nodes contributed by the newly added ones, as they would in a from-scratch run") parser.add_argument("--batch", action="store_true", @@ -114,14 +114,14 @@ def main(): if options.mgSplit and options.batch: raise RuntimeError("--mgSplit is for the whole-genome splitting pass and cannot be used with --batch") - if options.extendGAF: + if options.inGAF: if options.batch: - raise RuntimeError("--extendGAF cannot be used with --batch") + raise RuntimeError("--inGAF cannot be used with --batch") if options.collapse or options.collapseRefPAF: - raise RuntimeError("--extendGAF cannot be used with --collapse or --collapseRefPAF: collapse PAFs are minimap2 " + raise RuntimeError("--inGAF cannot be used with --collapse or --collapseRefPAF: collapse PAFs are minimap2 " "self-alignments and are not derived from the GAF") elif options.remap: - raise RuntimeError("--remap only means something with --extendGAF, which is what it overrides") + raise RuntimeError("--remap only means something with --inGAF, which is what it overrides") # Mess with some toil options to create useful defaults. cactus_override_toil_options(options) @@ -249,12 +249,12 @@ def graph_map(options): input_dict[chrom] = seq_id_map, gfa_id, ref_collapse_paf_id, input_map[chrom][0], input_map[chrom][1] #import the mappings to reuse - extend_gaf_id = toil.importFile(makeURL(options.extendGAF)) if options.extendGAF and not options.remap else None + in_gaf_id = toil.importFile(makeURL(options.inGAF)) if options.inGAF and not options.remap else None # run the workflow # output_dict is chrom -> paf_id, gfa_fa_id, gaf_id, unfiltered_paf_id, paf_filter_log, paf_was_filtered output_dict = toil.start(Job.wrapJobFn(minigraph_batch_separate_workflow, options, config_wrapper, input_dict, graph_event, True, - extend_gaf_id=extend_gaf_id)) + in_gaf_id=in_gaf_id)) export_graphmap_output(options, config_node, input_map, output_dict, toil) @@ -315,7 +315,7 @@ def export_graphmap_output(options, config_node, input_map, output_dict, toil): if chrom_file_path.startswith('s3://'): write_s3(chrom_file_temp_path, chrom_file_path) -def minigraph_batch_workflow(job, options, config, input_dict, graph_event, sanitize, pansn_gfa_input=True, extend_gaf_id=None): +def minigraph_batch_workflow(job, options, config, input_dict, graph_event, sanitize, pansn_gfa_input=True, in_gaf_id=None): """ Batch wrapper to run grpahmap independently at the chromosome level.""" output_dict = {} options.mg_chrom_name = None @@ -330,7 +330,7 @@ def minigraph_batch_workflow(job, options, config, input_dict, graph_event, sani else: chrom_options = options mgwf_job = job.addChildJobFn(minigraph_workflow, chrom_options, config, seq_id_map, gfa_id, graph_event, - sanitize, ref_collapse_paf_id, pansn_gfa_input, extend_gaf_id=extend_gaf_id) + sanitize, ref_collapse_paf_id, pansn_gfa_input, in_gaf_id=in_gaf_id) output_dict[chrom] = mgwf_job.rv() return output_dict @@ -346,15 +346,15 @@ def add_separate_ref_contigs_job(batch_job, options, config, input_dict): getattr(options, 'permissiveContigFilter', None), whole_genome_ref=getattr(options, 'mgSplitWholeGenomeRef', False)) -def minigraph_batch_separate_workflow(job, options, config, input_dict, graph_event, sanitize, pansn_gfa_input=True, extend_gaf_id=None): +def minigraph_batch_separate_workflow(job, options, config, input_dict, graph_event, sanitize, pansn_gfa_input=True, in_gaf_id=None): """ minigraph_batch_workflow followed by the separation pass, for callers that just want the final result and add nothing after it """ batch_job = job.addChildJobFn(minigraph_batch_workflow, options, config, input_dict, graph_event, sanitize, - pansn_gfa_input, extend_gaf_id=extend_gaf_id) + pansn_gfa_input, in_gaf_id=in_gaf_id) return add_separate_ref_contigs_job(batch_job, options, config, input_dict).rv() def minigraph_workflow(job, options, config, seq_id_map, gfa_id, graph_event, sanitize, ref_collapse_paf_id, pansn_gfa_input=True, - extend_gaf_id=None): + in_gaf_id=None): """ Overall workflow takes command line options and returns (paf-id, (optional) fa-id) """ fa_id = None gfa_id_size = gfa_id.size @@ -392,11 +392,11 @@ def minigraph_workflow(job, options, config, seq_id_map, gfa_id, graph_event, sa # split up any mappings we've been given to reuse, so each genome's re-derivation is its own # job just as its mapping would have been - extend_gaf_map = None - if extend_gaf_id: - split_gaf_job = root_job.addChildJobFn(split_gaf_by_event, extend_gaf_id, genome_names, options.extendGAF, - disk=12*extend_gaf_id.size) - extend_gaf_map = split_gaf_job.rv() + in_gaf_map = None + if in_gaf_id: + split_gaf_job = root_job.addChildJobFn(split_gaf_by_event, in_gaf_id, genome_names, options.inGAF, + disk=12*in_gaf_id.size) + in_gaf_map = split_gaf_job.rv() zipped_gfa = options.minigraphGFA.endswith('.gz') if options.outputFasta: @@ -412,18 +412,18 @@ def minigraph_workflow(job, options, config, seq_id_map, gfa_id, graph_event, sa gfa_id = gfa_unzip_job.rv() gfa_id_size *= 10 options.minigraphGFA = options.minigraphGFA[:-3] - if extend_gaf_id: + if in_gaf_id: # resolving a reused GAF is the same work for every genome, so a GAF that does not belong # to this graph fails identically in all of them -- once per genome, after the fan-out, and # again on every Toil retry. Checking a sample up front turns that into one quick failure # with something actionable in it. chained onto the unzip (when there is one) because it # needs the same uncompressed graph the per-genome jobs use check_parent = gfa_unzip_job if zipped_gfa else root_job - check_parent.addFollowOnJobFn(check_reusable_gaf, config, extend_gaf_id, gfa_id, genome_names, - options.extendGAF, options.minigraphGFA, + check_parent.addFollowOnJobFn(check_reusable_gaf, config, in_gaf_id, gfa_id, genome_names, + options.inGAF, options.minigraphGFA, disk=4*gfa_id_size, memory=cactus_clamp_memory(2*gfa_id_size)) - paf_job = Job.wrapJobFn(minigraph_map_all, options, config, gfa_id, seq_id_map, graph_event, extend_gaf_map) + paf_job = Job.wrapJobFn(minigraph_map_all, options, config, gfa_id, seq_id_map, graph_event, in_gaf_map) root_job.addFollowOn(paf_job) collapse_paf_id = ref_collapse_paf_id @@ -511,10 +511,10 @@ def make_minigraph_fasta(job, gfa_file_id, gfa_file_path, name): return job.fileStore.writeGlobalFile(fa_path) -def minigraph_map_all(job, options, config, gfa_id, fa_id_map, graph_event, extend_gaf_map=None): +def minigraph_map_all(job, options, config, gfa_id, fa_id_map, graph_event, in_gaf_map=None): """ top-level job to run the minigraph mapping in parallel, returns paf. - a genome that extend_gaf_map already has mappings for has its PAF re-derived from them rather + a genome that in_gaf_map already has mappings for has its PAF re-derived from them rather than being mapped again -- see translate_gaf_one() """ # hang everything on this job, to self-contain workflow top_job = Job() @@ -546,10 +546,10 @@ def minigraph_map_all(job, options, config, gfa_id, fa_id_map, graph_event, exte # the memory heuristc seems to drastically underestimate some chromosomes in batch mode... mem *= 2 event_name = '{}.{}'.format(event, options.mg_chrom_name) - if extend_gaf_map and event in extend_gaf_map: + if in_gaf_map and event in in_gaf_map: # no minigraph, and no input fasta: gaf2unstable/gaffilter/gaf2paf against the new graph # is the whole job. gaffilter reads its input into memory, as it does when mapping - gaf_shard_id = extend_gaf_map[event] + gaf_shard_id = in_gaf_map[event] map_job = top_job.addChildJobFn(translate_gaf_one, config, event_name, gaf_shard_id, gfa_id, genome_names, disk=12*gaf_shard_id.size + 2*gfa_id.size, memory=cactus_clamp_memory(24*gaf_shard_id.size + 4*gfa_id.size)) @@ -640,7 +640,7 @@ def split_gaf_file_by_event(gaf_path, names, shard_dir): re-derived by the very same code that produced them in the first place. genomes in the GAF that are not in names are dropped, not an error: --refFromGFA legitimately - takes the reference out of the sequence map. cactus-minigraph --extendGFA is where a genome + takes the reference out of the sequence map. cactus-minigraph --inGFA is where a genome missing from the seqfile is caught, because there it is unrecoverable """ prefix_map = pansn_to_event_map(names) @@ -792,8 +792,8 @@ def check_reusable_gaf(job, config, gaf_file_id, gfa_file_id, genome_names, gaf_ 'The first {} records only ever walk through: {}.\n' 'A --mgSplit run is the usual way to get a mismatched pair that looks like a matching one: its ' '.gaf.gz is the whole-genome pass against the reference-only first-pass graph, while its ' - '.sv.gfa.gz is the merged per-chromosome graphs. Drop --extendGAF to map every genome ' - 'against the extended graph instead -- --extendGFA still saves the construction, which is the ' + '.sv.gfa.gz is the merged per-chromosome graphs. Drop --inGAF to map every genome ' + 'against the extended graph instead -- --inGFA still saves the construction, which is the ' 'expensive half.\nUnderlying error: {}'.format( gaf_path, gfa_path, GAF_REUSE_CHECK_RECORDS, ' '.join(sorted(step_genomes)[:12]) or '(nothing)', e)) diff --git a/src/cactus/refmap/cactus_minigraph.py b/src/cactus/refmap/cactus_minigraph.py index e82d23a56..3fae34f09 100644 --- a/src/cactus/refmap/cactus_minigraph.py +++ b/src/cactus/refmap/cactus_minigraph.py @@ -57,11 +57,11 @@ def main(): help="Use last-train to estimate scoring matrix from input data", default=False) parser.add_argument("--refOnly", action="store_true", help="Only build the graph out of reference genome(s). Can be used when it will only be used for chromosome-splitting, for example") - parser.add_argument("--extendGFA", type=str, default=None, - help="Extend this existing minigraph GFA (as made by a previous cactus-minigraph or cactus-pangenome run) " + parser.add_argument("--inGFA", type=str, default=None, + help="Start from this existing minigraph GFA (as made by a previous cactus-minigraph or cactus-pangenome run) " "instead of building from scratch. Only the seqFile genomes that are not already in the graph get added, in " - "mash-distance order among themselves. The seqFile must still contain every genome in the graph: genomes " - "cannot be removed from a minigraph") + "mash-distance order among themselves; if there are none, the graph is passed through untouched. The seqFile " + "must still contain every genome in the graph: genomes cannot be removed from a minigraph") parser.add_argument("--batch", action="store_true", help="Run independently on set of chromosomea inputs (chromfile as from cactus-graphmap-split). Note that the output will be a directory and not a GFA") @@ -126,22 +126,22 @@ def main(): if '://' not in options.outputGFA: options.outputGFA = os.path.abspath(options.outputGFA) - extend_gfa_id = None - if options.extendGFA: + in_gfa_id = None + if options.inGFA: if options.batch: - raise RuntimeError('--extendGFA cannot be used with --batch') + raise RuntimeError('--inGFA cannot be used with --batch') if options.refOnly: - raise RuntimeError('--extendGFA cannot be used with --refOnly') - if '://' not in options.extendGFA: - options.extendGFA = os.path.abspath(options.extendGFA) - extend_gfa_id = toil.importFile(makeURL(options.extendGFA)) + raise RuntimeError('--inGFA cannot be used with --refOnly') + if '://' not in options.inGFA: + options.inGFA = os.path.abspath(options.inGFA) + in_gfa_id = toil.importFile(makeURL(options.inGFA)) # maps name -> input_seq_id_map, input_seq_order input_dict = minigraph_construct_import_sequences(options, config_wrapper, input_seqfiles, toil) # output_dict: chrom-> (gfa_id, pansn_gfa_id, train_id) output_dict = toil.start(Job.wrapJobFn(minigraph_construct_batch_workflow, options, config_node, input_dict, options.outputGFA, - extend_gfa_id=extend_gfa_id)) + in_gfa_id=in_gfa_id)) export_minigraph_construct_output(options, input_seqfiles, output_dict, toil) @@ -277,7 +277,7 @@ def check_sample_names(sample_names, references): raise RuntimeError("Sample name {} with \"{}\" suffix is not supported. You must either remove this suffix or use .N where N is an integer to specify haplotype".format(sample, sample_ext)) def minigraph_construct_batch_workflow(job, options, config_node, input_dict, gfa_path, sanitize=True, - construct_ref_id_map=None, extend_gfa_id=None): + construct_ref_id_map=None, in_gfa_id=None): """ run the construction workflow on individual chromosomes. construct_ref_id_map, if given, swaps the whole-genome reference fastas in for the chromosome's own slice of them (--mgSplit --mgSplitWholeGenomeRef). the merge happens here rather than inside minigraph_construct_workflow @@ -295,12 +295,12 @@ def minigraph_construct_batch_workflow(job, options, config_node, input_dict, gf else: gfa_path = options.outputGFA mgwf_job = job.addChildJobFn(minigraph_construct_workflow, options, config_node, seq_id_map, seq_order, gfa_path, sanitize, - construct_seq_id_map=construct_seq_id_map, extend_gfa_id=extend_gfa_id) + construct_seq_id_map=construct_seq_id_map, in_gfa_id=in_gfa_id) output_dict[chrom] = mgwf_job.rv() return output_dict def minigraph_construct_workflow(job, options, config_node, seq_id_map, seq_order, gfa_path, sanitize=True, - construct_seq_id_map=None, extend_gfa_id=None): + construct_seq_id_map=None, in_gfa_id=None): """ minigraph can handle bgzipped files but not gzipped; so unzip everything in case before running construct_seq_id_map, when given, replaces seq_id_map for the graph construction alone. it is how @@ -313,17 +313,17 @@ def minigraph_construct_workflow(job, options, config_node, seq_id_map, seq_orde with a graph to extend, which genomes still need constructing is not known until that graph's SN tags have been read, so the rest of the workflow is deferred behind the job that reads them """ - if not extend_gfa_id: + if not in_gfa_id: return minigraph_construct_run(job, options, config_node, seq_id_map, seq_order, gfa_path, sanitize, construct_seq_id_map=construct_seq_id_map) # the renaming pass decompresses the GFA before bgzipping it back up, so it needs room for # the raw copy (reckoned at 10x, as elsewhere) on top of the compressed input and output - rename_job = job.addChildJobFn(minigraph_gfa_from_pansn, set(seq_id_map.keys()), options.extendGFA, extend_gfa_id, - disk=extend_gfa_id.size*12) + rename_job = job.addChildJobFn(minigraph_gfa_from_pansn, set(seq_id_map.keys()), options.inGFA, in_gfa_id, + disk=in_gfa_id.size*12) run_job = rename_job.addFollowOnJobFn(minigraph_construct_run, options, config_node, seq_id_map, seq_order, gfa_path, sanitize, construct_seq_id_map, - rename_job.rv(0), rename_job.rv(1), extend_gfa_id) + rename_job.rv(0), rename_job.rv(1), in_gfa_id) return run_job.rv(0), run_job.rv(1), run_job.rv(2) def minigraph_construct_run(job, options, config_node, seq_id_map, seq_order, gfa_path, sanitize=True, @@ -336,29 +336,44 @@ def minigraph_construct_run(job, options, config_node, seq_id_map, seq_order, gf # the PanSN rename at the end of construction has to resolve every SN tag in the finished # graph, which on the extend path is more genomes than minigraph is being given graph_names = set(seq_id_map.keys()) + # last-training is over fastas, so it wants every genome, not just the ones still to construct + train_seq_id_map, train_seq_order = seq_id_map, seq_order if seed_events is not None: if options.reference[0] not in seed_events: # it would otherwise be constructed in last, at the highest rGFA rank rather than rank 0, # and every rank-0 assumption downstream would be reading the wrong genome - raise RuntimeError('Reference {} is not in the graph being extended, whose genomes are: {}. A graph can only be ' - 'extended with the reference it was built on'.format(options.reference[0], - ' '.join(sorted(seed_events)))) + raise RuntimeError('Reference {} is not in {}, whose genomes are: {}. A graph can only be reused with the ' + 'reference it was built on'.format(options.reference[0], options.inGFA, + ' '.join(sorted(seed_events)))) # everything already in the seed graph is left alone: minigraph only gets the genomes that # are new to it, appended after the ones the graph was built from. the reference is kept in # the sequence map (but not the order) because the mash sort below still sketches against it seq_order = [seq for seq in seq_order if seq not in seed_events] seq_id_map = {name: fa_id for name, fa_id in seq_id_map.items() if name not in seed_events or name == options.reference[0]} - RealtimeLogger.info('Extending a graph of {} genomes with {}: {}'.format( - len(seed_events), len(seq_order), ' '.join(seq_order) if seq_order else '(nothing)')) - if not seq_order: - # nothing to add, so the graph handed back is the one --extendGFA was given. its - # compression follows the input name, and everything downstream reads the *output* - # name to decide whether to unzip, so it has to be re-emitted to match + if seq_order: + RealtimeLogger.info('Extending the {} genomes in {} with {}: {}'.format( + len(seed_events), options.inGFA, len(seq_order), ' '.join(seq_order))) + else: + # the seqfile asks for exactly the genomes the graph already holds, so there is nothing + # to construct and the run resumes from it. the graph handed back is the one --inGFA + # was given, but its compression follows the *input* name while everything downstream + # reads the output name to decide whether to unzip, so it is re-emitted to match + RealtimeLogger.info('Resuming from the {} genomes in {}: nothing left to construct'.format( + len(seed_events), options.inGFA)) match_job = job.addChildJobFn(match_gfa_compression, seed_gfa_id, seed_pansn_gfa_id, - options.extendGFA, gfa_path, + options.inGFA, gfa_path, disk=12 * (seed_gfa_id.size if hasattr(seed_gfa_id, 'size') else 0)) - return match_job.rv(0), match_job.rv(1), None + # last_train reads fastas, not the graph, so resuming is no reason to skip it: without + # this --lastTrain would quietly fall back to the default scoring matrix + train_id = None + if options.lastTrain and len(train_seq_id_map) > 1: + train_job = job.addChildJobFn(last_train, config_node, train_seq_order, train_seq_id_map, + ref_name=options.reference[0], + cores=options.mgCores, disk=8*ref_size, + memory=cactus_clamp_memory(max(8*ref_size, 12*10**9))) + train_id = train_job.rv() + return match_job.rv(0), match_job.rv(1), train_id else: assert options.reference[0] == seq_order[0] if options.refOnly: @@ -627,7 +642,7 @@ def minigraph_construct_in_batches(job, options, config_node, seq_id_map, seq_or if seed_gfa_id: # minigraph_construct() only uses this to name its local copy, but keep the compression # suffix honest since that is what says whether the file it reads is bgzipped - seed_gfa_path = 'extend.gfa.gz' if options.extendGFA.endswith('.gz') else 'extend.gfa' + seed_gfa_path = 'extend.gfa.gz' if options.inGFA.endswith('.gz') else 'extend.gfa' for i in range(num_batches): batch_size = len(seq_order) - i * max_batch_size if i == num_batches - 1 else max_batch_size input_seq_order = seq_order[i * max_batch_size : (i * max_batch_size) + batch_size] @@ -766,7 +781,7 @@ def minigraph_gfa_from_pansn(job, names, gfa_path, gfa_id): so that a minigrpah GFA (as converted panSN by minigraph_gfa_to_pansn() above) can be read back into Cactus returns (converted gfa id, set of genomes the graph was built from). the genome set is what - --extendGFA needs to work out which of the seqfile's genomes are new to the graph, and it comes + --inGFA needs to work out which of the seqfile's genomes are new to the graph, and it comes free with the pass that has to read every SN tag anyway. a GFA that is already in cactus naming -- from a cactus old enough to have published one, or @@ -815,7 +830,7 @@ def minigraph_gfa_from_pansn(job, names, gfa_path, gfa_id): if '{}.{}'.format(name, hap) in names: name = '{}.{}'.format(name, hap) else: - # collected rather than asserted on: with --extendGFA this is usually + # collected rather than asserted on: with --inGFA this is usually # the user leaving a genome out of the seqfile, which deserves to be # named. the whole tag goes in the message because the other way to # land here is an SN tag that is not SAMPLE#HAP#CONTIG at all, and diff --git a/src/cactus/refmap/cactus_pangenome.py b/src/cactus/refmap/cactus_pangenome.py index f633d5c16..513a24b1e 100644 --- a/src/cactus/refmap/cactus_pangenome.py +++ b/src/cactus/refmap/cactus_pangenome.py @@ -70,19 +70,20 @@ def pangenome_options(parser): 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.") - parser.add_argument("--extendGFA", type=str, default=None, - help = "Add genomes to this existing pangenome's minigraph GFA (.sv.gfa.gz from a previous run, or a " - "published release) instead of building one from scratch. The seqFile must list every genome in the graph as " - "well as the ones being added: genomes cannot be removed from a minigraph. Only the new genomes are constructed " - "in, which is where nearly all of the minigraph cost is. Everything from cactus-graphmap-split on is recomputed") - parser.add_argument("--extendGAF", type=str, default=None, - help = "Reuse the mappings in this GAF (.gaf.gz from the same run that produced --extendGFA) rather than " - "re-running minigraph for the genomes it covers. Minigraph GAF is in stable coordinates, which adding genomes to " - "a graph does not change, so they are simply re-derived against the extended graph. Without this, every genome is " - "mapped again (see --remap)") + 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 " + "the rest is left alone; if it names exactly the genomes the graph already holds, nothing is constructed and the " + "run resumes from it. Either way the seqFile must list every genome in the graph, since genomes cannot be removed " + "from a minigraph. Construction is where nearly all of the minigraph cost is") + parser.add_argument("--inGAF", type=str, default=None, + help = "Reuse the mappings in this GAF (.gaf.gz from the run that produced --inGFA) rather than re-running " + "minigraph for the genomes it covers. Minigraph GAF is in stable coordinates, which adding genomes to a graph does " + "not change, so they are re-derived against the graph instead. With a graph that needs nothing constructed this " + "skips straight to cactus-graphmap-split. Without it, every genome is mapped again (see --remap)") parser.add_argument("--remap", action="store_true", default=False, - help = "With --extendGAF, map every genome with minigraph anyway. Costs the full mapping stage, but the existing " - "genomes then see the nodes contributed by the newly added ones, as they would in a from-scratch run") + help = "With --inGAF, map every genome with minigraph anyway. Costs the full mapping stage, but every genome then " + "sees the whole graph, as it would in a from-scratch run") # cactus-graphmap options parser.add_argument("--mapCores", type=int, help = "Number of cores for minigraph map. Overrides graphmap cpu in configuration") @@ -211,19 +212,19 @@ def pangenome_validate_options(options): if options.mgSplit and options.noSplit: raise RuntimeError('you cannot use both --mgSplit and --noSplit together: pick one') - if options.extendGFA: + if options.inGFA: if options.mgSplit: - raise RuntimeError('--extendGFA cannot (yet) be used with --mgSplit: the per-chromosome graphs and mappings would ' - 'need extending too') + raise RuntimeError('--inGFA cannot (yet) be used with --mgSplit: the per-chromosome graphs and mappings would ' + 'have to be carried over too') if options.collapse or options.collapseRefPAF: - raise RuntimeError('--extendGFA cannot be used with --collapse or --collapseRefPAF: collapse PAFs are minimap2 ' + raise RuntimeError('--inGFA cannot be used with --collapse or --collapseRefPAF: collapse PAFs are minimap2 ' 'self-alignments and are not derived from the GAF') else: - if options.extendGAF: - raise RuntimeError('--extendGAF requires --extendGFA: reusing mappings only makes sense against the graph they ' + if options.inGAF: + raise RuntimeError('--inGAF requires --inGFA: reusing mappings only makes sense against the graph they ' 'were made from') if options.remap: - raise RuntimeError('--remap only means something with --extendGAF, which is what it overrides') + raise RuntimeError('--remap only means something with --inGAF, which is what it overrides') # Sort out the graphmap-join options, which can be rather complex # pass in dummy values for now, they will get filled in later @@ -327,15 +328,15 @@ def main(): last_scores_id = toil.importFile(makeURL(options.scoresFile)) #import the pangenome being extended - extend_gfa_id, extend_gaf_id = None, None - if options.extendGFA: - if '://' not in options.extendGFA: - options.extendGFA = os.path.abspath(options.extendGFA) - extend_gfa_id = toil.importFile(makeURL(options.extendGFA)) - if options.extendGAF and not options.remap: - if '://' not in options.extendGAF: - options.extendGAF = os.path.abspath(options.extendGAF) - extend_gaf_id = toil.importFile(makeURL(options.extendGAF)) + in_gfa_id, in_gaf_id = None, None + if options.inGFA: + if '://' not in options.inGFA: + options.inGFA = os.path.abspath(options.inGFA) + in_gfa_id = toil.importFile(makeURL(options.inGFA)) + if options.inGAF and not options.remap: + if '://' not in options.inGAF: + options.inGAF = os.path.abspath(options.inGAF) + in_gaf_id = toil.importFile(makeURL(options.inGAF)) #import the sequences input_seq_id_map = {} @@ -354,7 +355,7 @@ def main(): input_seq_order.remove(genome) toil.start(Job.wrapJobFn(pangenome_end_to_end_workflow, options, config_wrapper, input_seq_id_map, input_path_map, input_seq_order, ref_collapse_paf_id, last_scores_id, - extend_gfa_id=extend_gfa_id, extend_gaf_id=extend_gaf_id)) + in_gfa_id=in_gfa_id, in_gaf_id=in_gaf_id)) end_time = timeit.default_timer() run_time = end_time - start_time @@ -606,7 +607,7 @@ def export_join_wrapper(job, options, wf_output, contig_sizes_id=None): job.fileStore.exportFile(contig_sizes_id, makeURL(sizes_path)) def pangenome_end_to_end_workflow(job, options, config_wrapper, seq_id_map, seq_path_map, seq_order, ref_collapse_paf_id, - last_scores_id, extend_gfa_id=None, extend_gaf_id=None): + last_scores_id, in_gfa_id=None, in_gaf_id=None): """ chain the entire workflow together, doing exports after each step to mitigate annoyance of failures """ root_job = Job() job.addChild(root_job) @@ -664,7 +665,7 @@ def pangenome_end_to_end_workflow(job, options, config_wrapper, seq_id_map, seq_ split_config_node = config_node split_config_wrapper = config_wrapper minigraph_job = prev_job.addFollowOnJobFn(minigraph_construct_workflow, mg_options, split_config_node, seq_id_map, seq_order, sv_gfa_path, sanitize=False, - extend_gfa_id=extend_gfa_id) + in_gfa_id=in_gfa_id) sv_gfa_id = minigraph_job.rv(0) pansn_sv_gfa_id = minigraph_job.rv(1) if not last_scores_id: @@ -682,7 +683,7 @@ def pangenome_end_to_end_workflow(job, options, config_wrapper, seq_id_map, seq_ if options.mgSplit: gm_options.collapse = False graphmap_job = minigraph_wrapper_job.addFollowOnJobFn(minigraph_workflow, gm_options, split_config_wrapper, seq_id_map, sv_gfa_id, graph_event, False, ref_collapse_paf_id, pansn_gfa_input=False, - extend_gaf_id=extend_gaf_id) + in_gaf_id=in_gaf_id) paf_id, gfa_fa_id, gaf_id, unfiltered_paf_id, paf_filter_log = graphmap_job.rv(0), graphmap_job.rv(1), graphmap_job.rv(2), graphmap_job.rv(3), graphmap_job.rv(4) graphmap_export_job = graphmap_job.addFollowOnJobFn(export_graphmap_wrapper, options, paf_id, paf_path, gaf_id, unfiltered_paf_id, paf_filter_log) diff --git a/src/cactus/refmap/cactus_panpatch.py b/src/cactus/refmap/cactus_panpatch.py index 06057c8a4..472b79804 100644 --- a/src/cactus/refmap/cactus_panpatch.py +++ b/src/cactus/refmap/cactus_panpatch.py @@ -367,11 +367,11 @@ def panpatch_validate_options(options): raise RuntimeError('--noSplit cannot be used with cactus-panpatch: panpatch needs one graph per reference chromosome') # getattr because the unit tests build a minimal namespace rather than going through the parser - if getattr(options, 'extendGFA', None) or getattr(options, 'extendGAF', None) or getattr(options, 'remap', False): + if getattr(options, 'inGFA', None) or getattr(options, 'inGAF', None) or getattr(options, 'remap', False): # these come in via pangenome_options(), which panpatch shares. the graph panpatch builds # is a throwaway, built per sample out of that sample and its donors, so there is nothing # an earlier run could usefully be extended from - raise RuntimeError('--extendGFA / --extendGAF / --remap cannot be used with cactus-panpatch: it builds a fresh ' + raise RuntimeError('--inGFA / --inGAF / --remap cannot be used with cactus-panpatch: it builds a fresh ' 'graph per sample being patched, which is not an extension of anything') def disable_pangenome_outputs(options): diff --git a/test/evolverTest.py b/test/evolverTest.py index 9299fae1b..98ec8f82c 100644 --- a/test/evolverTest.py +++ b/test/evolverTest.py @@ -503,12 +503,13 @@ def _run_evolver_primates_graphmap(self, binariesMode): subprocess.check_call(['cactus-align', self._job_store(binariesMode), seq_file_fix_path, paf_path, self._out_hal(binariesMode), '--pangenome', '--outVG', '--outGFA', '--pafMaskFilter', '10000', '--barMaskFilter', '10000'] + cactus_opts) - def _run_evolver_primates_pangenome(self, binariesMode, mgSplit = False, extend = False): + def _run_evolver_primates_pangenome(self, binariesMode, mgSplit = False, extend = False, resume = False): """ run the primates start in using high-level cactus-pangenome interface. with extend, half the genomes are built into a graph with the step-by-step tools first and - cactus-pangenome --extendGFA adds the rest, which also checks that a graph made one way can - be extended the other """ + cactus-pangenome --inGFA adds the rest, which also checks that a graph made one way can + be extended the other. with resume, that graph holds every genome in the seqfile, so there + is nothing to construct and nothing to map and the run picks up at cactus-graphmap-split """ # borrow seqfile from other primates test # todo: make a seqfile and add it to the repo seq_file_path = os.path.join(self.tempDir, 'primates.txt') @@ -531,12 +532,16 @@ def _run_evolver_primates_pangenome(self, binariesMode, mgSplit = False, extend out_name = os.path.splitext(os.path.basename(self._out_hal(binariesMode)))[0] cactus_pangenome_cmd = ['cactus-pangenome', self._job_store(binariesMode), seq_file_path, '--reference', 'simHuman', 'simChimp', '--outDir', out_dir, '--outName', out_name, '--odgi', '--chrom-og', '--viz', '--draw', '--haplo', '--lastTrain'] - if not extend: + if not extend and not resume: # collapse self-alignments are not derived from the GAF, so they cannot be reused cactus_pangenome_cmd += ['--collapse'] + elif resume: + # a graph of every genome in the seqfile: nothing to construct, nothing to map + base = self._build_primates_base_graph(binariesMode, mc_config_path, self.PRIMATES) + cactus_pangenome_cmd += ['--inGFA', base['gfa'], '--inGAF', base['gaf']] else: base = self._build_primates_base_graph(binariesMode, mc_config_path, ['simHuman', 'simChimp']) - cactus_pangenome_cmd += ['--extendGFA', base['gfa'], '--extendGAF', base['gaf']] + cactus_pangenome_cmd += ['--inGFA', base['gfa'], '--inGAF', base['gaf']] if mgSplit: cactus_pangenome_cmd += ['--mgSplit'] else: @@ -551,15 +556,19 @@ def _run_evolver_primates_pangenome(self, binariesMode, mgSplit = False, extend wave_vcf_bytes = os.path.getsize(os.path.join(out_dir, out_name + '.simChimp.wave.vcf.gz')) self.assertGreaterEqual(wave_vcf_bytes, 300000) - if extend: - # the genomes added on top of the base graph are in the graph it was extended into + if extend or resume: + # extending adds the missing genomes; resuming hands the same graph straight back self.assertEqual(self._gfa_genomes(os.path.join(out_dir, out_name + '.sv.gfa.gz')), set(self.PRIMATES)) + if resume: + # resuming must not rebuild the graph: what comes out is what went in + self.assertEqual(self._gfa_text(os.path.join(out_dir, out_name + '.sv.gfa.gz')), + self._gfa_text(base['gfa'])) if mgSplit: # the whole-genome first pass is reference-only, and cactus-graphmap-join later writes # the merged per-chromosome graph to .sv.gfa.gz. Both have to survive under # names that say which is which, or the first pass's mappings are left beside a graph - # they do not describe -- which is exactly what breaks --extendGAF + # they do not describe -- which is exactly what breaks --inGAF refonly_gfa = os.path.join(out_dir, out_name + '.refonly.sv.gfa.gz') self.assertTrue(os.path.exists(refonly_gfa), 'no {} from the --mgSplit first pass'.format(refonly_gfa)) self.assertEqual(self._gfa_genomes(refonly_gfa), {'simHuman', 'simChimp'}) @@ -609,10 +618,10 @@ def _gfa_genomes(self, gfa_path): def _run_primates_extend_steps(self, binariesMode, seqfile_genomes, extend_genomes, config_path, out_prefix): """ cactus-minigraph + cactus-graphmap on seqfile_genomes, then the same two commands again - on extend_genomes with --extendGFA/--extendGAF pointed at the first run's output. + on extend_genomes with --inGFA/--inGAF pointed at the first run's output. returns (base outputs, extended outputs) as dicts of gfa/paf/gaf paths. only the two - stages that --extendGFA/--extendGAF touch are run: split/align/join are unchanged by them + stages that --inGFA/--inGAF touch are run: split/align/join are unchanged by them and are covered by the end-to-end test below """ cactus_opts = ['--binariesMode', binariesMode, '--logInfo', '--workDir', self.tempDir, '--configFile', config_path] @@ -635,8 +644,8 @@ def run_stage(tag, genomes, extend_from): gm_cmd = ['cactus-graphmap', os.path.join(self.tempDir, 'js-{}-{}-gm'.format(out_prefix, tag)), seqfile, out['gfa'], out['paf'], '--outputFasta', out['fa'], '--reference', 'simChimp'] if extend_from: - mg_cmd += ['--extendGFA', extend_from['gfa']] - gm_cmd += ['--extendGAF', extend_from['gaf']] + mg_cmd += ['--inGFA', extend_from['gfa']] + gm_cmd += ['--inGAF', extend_from['gaf']] subprocess.check_call(mg_cmd + cactus_opts) subprocess.check_call(gm_cmd + cactus_opts) return out @@ -658,7 +667,7 @@ def _remap_onto_extended(self, extended, out_prefix): subprocess.check_call(['cactus-graphmap', os.path.join(self.tempDir, 'js-{}-remap'.format(out_prefix)), seqfile, extended['gfa'], paf, '--outputFasta', os.path.join(work_dir, 'pg.sv.gfa.fa.gz'), - '--reference', 'simChimp', '--extendGAF', extended['gaf'], '--remap'] + '--reference', 'simChimp', '--inGAF', extended['gaf'], '--remap'] + extended['cactus_opts']) return paf @@ -791,7 +800,7 @@ def _run_yeast_pangenome(self, binariesMode, mgSplit=False, wholeGenomeRef=False """ yeast pangenome chromosome by chromosome pipeline, as run through a single invocations. extend, if given, is the list of genomes to build a graph out of first, with the - step-by-step tools, for the run below to extend with --extendGFA. unlike the primates + step-by-step tools, for the run below to extend with --inGFA. unlike the primates tests this exercises the translated PAF through cactus-graphmap-split, which is the one downstream stage that reads it before cactus-align """ @@ -822,7 +831,7 @@ def _run_yeast_pangenome(self, binariesMode, mgSplit=False, wholeGenomeRef=False cactus_pangenome_cmd += ['--vcfL', str(vcfL)] if extend: base = self._build_yeast_base_graph(binariesMode, orig_seq_file_path, extend) - cactus_pangenome_cmd += ['--extendGFA', base['gfa'], '--extendGAF', base['gaf']] + cactus_pangenome_cmd += ['--inGFA', base['gfa'], '--inGAF', base['gaf']] subprocess.check_call(cactus_pangenome_cmd + cactus_opts) #compatibility with older test @@ -1895,7 +1904,7 @@ def testPangenomeExtendNullLocal(self): plain_gfa = os.path.join(plain_dir, 'pg.sv.gfa') subprocess.check_call(['cactus-minigraph', os.path.join(self.tempDir, 'js-null-plain'), plain_seqfile, plain_gfa, '--reference', 'simChimp', - '--extendGFA', base['gfa'], + '--inGFA', base['gfa'], '--binariesMode', 'local', '--logInfo', '--workDir', self.tempDir, '--configFile', config_path]) with open(plain_gfa, 'rb') as plain_file: @@ -1912,7 +1921,7 @@ def testPangenomeExtendConstructionLocal(self): pinned (minigraphSortInput=none), adding simGorilla and simOrang to a simChimp+simHuman graph is the same command sequence cactus already runs when it batches construction, so the two graphs should be identical. If this ever stops holding, the cost model behind - --extendGFA needs revisiting. """ + --inGFA needs revisiting. """ config_path = self._write_pangenome_config('extend-construct', graphmap_attribs={"minigraphConstructBatchSize": "2", "minigraphSortInput": "none"}) @@ -1968,7 +1977,7 @@ def _paf_bytes(self, paf_path): return paf_file.read() def testEvolverPrimatesPangenomeExtendLocal(self): - """ Evolver (star) primates built two genomes at a time with cactus-pangenome --extendGFA, + """ Evolver (star) primates built two genomes at a time with cactus-pangenome --inGFA, checked against the same accuracy baseline as the from-scratch run """ name = "local" self._run_evolver_primates_pangenome(name, extend=True) @@ -1977,6 +1986,19 @@ def testEvolverPrimatesPangenomeExtendLocal(self): # todo: tune config so that delta can be reduced self._check_maf_accuracy(self._out_hal("local"), delta=(0.025,0.025), dataset='primates') + def testEvolverPrimatesPangenomeResumeLocal(self): + """ cactus-pangenome --inGFA/--inGAF on a graph that already holds every genome in the + seqfile: nothing to construct, nothing to map, so the run picks up at + cactus-graphmap-split. This is the whole point of --inGFA being an input rather than an + "extend" flag -- it replaces driving cactus-graphmap-split / cactus-align --batch / + cactus-graphmap-join by hand. """ + name = "local" + self._run_evolver_primates_pangenome(name, resume=True) + + # check the output + # todo: tune config so that delta can be reduced + self._check_maf_accuracy(self._out_hal("local"), delta=(0.025,0.025), dataset='primates') + def testEvolverPrimatesPangenomeStepByStepSplitLocal(self): """ Evolver primates but using the step-by-step cactus-pangenome interface with splitting """ @@ -2015,7 +2037,7 @@ def testYeastPangenomeLocal(self): self._check_yeast_pangenome(name, other_ref='DBVPG6044', expect_odgi=True, expect_haplo=False, expect_unchopped_gfa=True) def testYeastPangenomeExtendLocal(self): - """ Yeast pangenome built three strains at a time with cactus-pangenome --extendGFA. + """ Yeast pangenome built three strains at a time with cactus-pangenome --inGFA. Unlike the primates extend tests, this one splits by chromosome, so the reused mappings go through cactus-graphmap-split -- the one stage downstream of cactus-graphmap that reads the