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