diff --git a/Makefile b/Makefile index 58bfd58ce..8993e78eb 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,21 @@ 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_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 + 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..69c5f4c12 100644 --- a/doc/pangenome.md +++ b/doc/pangenome.md @@ -267,6 +267,161 @@ 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`. +### Starting From an Existing Graph or Mappings + +`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 \ + --inGFA pg1/pg.sv.gfa.gz --inGAF pg1/pg.gaf.gz +``` + +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 \ + --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 --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 is an ordinary +pangenome that can itself be used as an `--inGFA`/`--inGAF` input again. + +#### 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 +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. `--inGAF` is the option that does this, and it turns the mapping stage from hours +of `minigraph` into minutes of file conversion. + +`--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 +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 `--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 +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 + +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 +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 `--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 --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 + `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 a4b19823a..f5ef7fcec 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("--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 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 --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", 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.inGAF: + if options.batch: + raise RuntimeError("--inGAF cannot be used with --batch") + if options.collapse or options.collapseRefPAF: + 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 --inGAF, 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 + 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)) + output_dict = toil.start(Job.wrapJobFn(minigraph_batch_separate_workflow, options, config_wrapper, input_dict, graph_event, True, + in_gaf_id=in_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, in_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, in_gaf_id=in_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, 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) + 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): +def minigraph_workflow(job, options, config, seq_id_map, gfa_id, graph_event, sanitize, ref_collapse_paf_id, pansn_gfa_input=True, + 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 @@ -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 + 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: @@ -381,7 +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] - paf_job = Job.wrapJobFn(minigraph_map_all, options, config, gfa_id, seq_id_map, graph_event) + 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, 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, in_gaf_map) root_job.addFollowOn(paf_job) collapse_paf_id = ref_collapse_paf_id @@ -469,8 +511,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, in_gaf_map=None): + """ top-level job to run the minigraph mapping in parallel, returns paf. + + 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() job.addChild(top_job) @@ -488,6 +533,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 +546,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 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 = 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)) + 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 +578,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 --inGFA 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 +744,175 @@ 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) + +# 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 --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)) + + 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. + + 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. + + Nothing in cactus-gfa-tools does this today. gaffilter's rebase_path() is the same algorithm in + C++, but it runs only on the records -t actually trims, and only when trimming is on at all, so + it never sees the reuse case. gaf2unstable, which is the thing changing the granularity, is + where this belongs if it ever moves. """ + 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' @@ -608,7 +950,23 @@ def minigraph_map_one(job, config, event_name, fa_file_id, gfa_file_id): overlap_cmd += ['-t', '-e', str(trim_edge), '-Q', str(trim_min_mapq), '-l', mg_lengths_path] cmd = [cmd, overlap_cmd] - 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 reused. 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..3fae34f09 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("--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; 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") @@ -121,11 +126,22 @@ def main(): if '://' not in options.outputGFA: options.outputGFA = os.path.abspath(options.outputGFA) + in_gfa_id = None + if options.inGFA: + if options.batch: + raise RuntimeError('--inGFA cannot be used with --batch') + if options.refOnly: + 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)) + output_dict = toil.start(Job.wrapJobFn(minigraph_construct_batch_workflow, options, config_node, input_dict, options.outputGFA, + in_gfa_id=in_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, 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 @@ -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, 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): + 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 @@ -293,12 +309,73 @@ 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 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.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), 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, + 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()) + # 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 {}, 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]} + 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.inGFA, gfa_path, + disk=12 * (seed_gfa_id.size if hasattr(seed_gfa_id, 'size') else 0)) + # 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: refonly_seq_id_map = {} refonly_seq_order = [] @@ -309,7 +386,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 +396,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 +405,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 +420,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 +488,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 +541,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 +583,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 +592,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 +638,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.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] @@ -529,8 +656,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 +671,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 +720,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 +780,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 + --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 + 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 +797,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 +827,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 --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 + # 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 +846,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..513a24b1e 100644 --- a/src/cactus/refmap/cactus_pangenome.py +++ b/src/cactus/refmap/cactus_pangenome.py @@ -70,6 +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("--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 --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") @@ -198,6 +212,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.inGFA: + if options.mgSplit: + 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('--inGFA cannot be used with --collapse or --collapseRefPAF: collapse PAFs are minimap2 ' + 'self-alignments and are not derived from the GAF') + else: + 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 --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 # (but we want to do as much error-checking upfront as possible) @@ -299,6 +327,17 @@ def main(): if options.scoresFile: last_scores_id = toil.importFile(makeURL(options.scoresFile)) + #import the pangenome being extended + 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 = {} input_path_map = {} @@ -315,7 +354,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, + in_gfa_id=in_gfa_id, in_gaf_id=in_gaf_id)) end_time = timeit.default_timer() run_time = end_time - start_time @@ -567,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): + 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) @@ -598,7 +638,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 @@ -617,7 +664,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, + 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: @@ -626,15 +674,16 @@ 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_") 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, + 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 d3e1762b8..472b79804 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, '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('--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): """ 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..98ec8f82c 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,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): - """ run the primates start in using high-level cactus-pangenome interface """ + 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 --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') @@ -491,26 +522,26 @@ 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 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 += ['--inGFA', base['gfa'], '--inGAF', base['gaf']] if mgSplit: cactus_pangenome_cmd += ['--mgSplit'] else: @@ -525,6 +556,121 @@ 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 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 --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'}) + 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 """ + 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 --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 --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] + + 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 += ['--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 + + 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', '--inGAF', 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 +795,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 --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 """ orig_seq_file_path = './examples/yeastPangenome.txt' @@ -677,12 +829,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 += ['--inGFA', base['gfa'], '--inGAF', 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 +1871,134 @@ 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', + '--inGFA', 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 + --inGFA 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 --inGFA, + 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 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 """ @@ -1734,6 +2036,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 --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 + 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'