From 367d264bd3341575b571d1d263bfd33455c6828a Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 30 Aug 2026 23:16:58 -0400 Subject: [PATCH] Trim the tangled fringe at satellite boundaries Both aligners work by extending anchors, so an anchor in unique sequence beside a satellite array gets extended into it, and inside the array the orthologous alignment is undefined. Each haplotype takes some arbitrary but locally high-scoring path for a while past the boundary, and they disagree with each other. clip-vg's -u removes the runs longer than its threshold; what survives is a tangle at every boundary, in runs too short for -u to reach. clip-vg v1.2.2 adds -k/--flank for this, and this wires it up: --clipFlank alongside --clip, how far to trim outward. This is the control that matters -- inside a repeat the gate seldom closes on its own, so the cap decides how much goes. clipFlank config, default 50000, so it is on by default. clipFlankThreshold config, negative to calibrate per graph, which is the default and what you want. One command-line option, the rest config, because only the cap is worth reaching for from a command line. --clipFlank rejects a negative value and rejects being combined with --clip 0; both would otherwise fall through the same test as 0 and disable the trim silently, which next to clipFlankThreshold -- where negative is the meaningful "calibrate" value -- is an easy thing to type. The calibrated threshold decides how much sequence goes and clip-vg reports it on stderr, which cactus_call only surfaces when a command fails, so clip_vg() captures it into a flank-stats.tsv per chromosome. The threshold gets its own column because that is the number worth reading down; the whole message sits beside it so a rewording upstream costs the column rather than the record. Validated whole-genome on HPRC v2.1 against CHM13, on the same chrom-alignments as the release so only the clipping differs. Graph: nodes -1.53%, edges -1.98%, length -1.49%, concentrated where it should be (chrY -9.4% of edges, the acrocentrics next, chr6 -0.1%). VCF: raw -0.92%, final -0.61%. Accuracy against dipcall truth over 231 HPRC samples, aardvark: F1 improves in 231 of 231, driven by recall. Needs hal2vg v1.2.2, pinned here: v1.2.1 predates -k, so with clipFlank on by default every run using downloaded binaries died at the clip phase on an unknown option. All six hal2vg binaries move together, since they ship from one release. Co-Authored-By: Claude Opus 5 (1M context) --- build-tools/downloadPangenomeTools | 12 +++--- src/cactus/cactus_progressive_config.xml | 6 ++- src/cactus/refmap/cactus_graphmap_join.py | 48 ++++++++++++++++++++++- src/cactus/refmap/cactus_pangenome.py | 4 +- 4 files changed, 60 insertions(+), 10 deletions(-) diff --git a/build-tools/downloadPangenomeTools b/build-tools/downloadPangenomeTools index 802741bde..8fdd2e8da 100755 --- a/build-tools/downloadPangenomeTools +++ b/build-tools/downloadPangenomeTools @@ -252,7 +252,7 @@ fi # hal2vg cd ${pangenomeBuildDir} -wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.1/hal2vg +wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.2/hal2vg chmod +x hal2vg if [[ $STATIC_CHECK -ne 1 || $(ldd hal2vg | grep so | wc -l) -eq 0 ]] then @@ -262,7 +262,7 @@ else fi # clip-vg cd ${pangenomeBuildDir} -wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.1/clip-vg +wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.2/clip-vg chmod +x clip-vg if [[ $STATIC_CHECK -ne 1 || $(ldd clip-vg | grep so | wc -l) -eq 0 ]] then @@ -272,7 +272,7 @@ else fi # halRemoveDupes cd ${pangenomeBuildDir} -wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.1/halRemoveDupes +wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.2/halRemoveDupes chmod +x halRemoveDupes if [[ $STATIC_CHECK -ne 1 || $(ldd halRemoveDupes | grep so | wc -l) -eq 0 ]] then @@ -282,7 +282,7 @@ else fi # halMergeChroms cd ${pangenomeBuildDir} -wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.1/halMergeChroms +wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.2/halMergeChroms chmod +x halMergeChroms if [[ $STATIC_CHECK -ne 1 || $(ldd halMergeChroms | grep so | wc -l) -eq 0 ]] then @@ -293,7 +293,7 @@ fi # halUnclip cd ${pangenomeBuildDir} -wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.1/halUnclip +wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.2/halUnclip chmod +x halUnclip if [[ $STATIC_CHECK -ne 1 || $(ldd halUnclip | grep so | wc -l) -eq 0 ]] then @@ -304,7 +304,7 @@ fi # filter-paf-deletions cd ${pangenomeBuildDir} -wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.1/filter-paf-deletions +wget -nv --tries=5 --waitretry=3 --timeout=60 --retry-connrefused https://github.com/ComparativeGenomicsToolkit/hal2vg/releases/download/v1.2.2/filter-paf-deletions chmod +x filter-paf-deletions if [[ $STATIC_CHECK -ne 1 || $(ldd filter-paf-deletions | grep so | wc -l) -eq 0 ]] then diff --git a/src/cactus/cactus_progressive_config.xml b/src/cactus/cactus_progressive_config.xml index a21576163..175ed4c5a 100644 --- a/src/cactus/cactus_progressive_config.xml +++ b/src/cactus/cactus_progressive_config.xml @@ -457,6 +457,8 @@ + + @@ -477,7 +479,9 @@ maxNodeLength="1024" gfaffix="1" pathNormalize="1" - clipNonMinigraph="1" + clipNonMinigraph="1" + clipFlank="50000" + clipFlankThreshold="-1" minimizerOptions="-k 29 -w 11" lrMinimizerOptions="-k 31 -w 50 -W" minFilterFragment="1000" diff --git a/src/cactus/refmap/cactus_graphmap_join.py b/src/cactus/refmap/cactus_graphmap_join.py index 6dd3aef31..40ace0993 100644 --- a/src/cactus/refmap/cactus_graphmap_join.py +++ b/src/cactus/refmap/cactus_graphmap_join.py @@ -162,6 +162,8 @@ def graphmap_join_options(parser): """ we share these options with cactus-pangenome """ parser.add_argument("--clip", type=int, default=10000, help = "Generate clipped graph by removing anything longer than this amount that is unaligned to the underlying minigraph. Set to 0 to disable (must also set --filter 0 as well). [default=10000]") + parser.add_argument("--clipFlank", type=int, help = "Trim the tangled fringe outward from each clipped interval by up to this many bp. 0 to disable (overrides config option of same name)") + parser.add_argument("--filter", type=int, default=2, help = "Generate a frequency filtered graph (from the clipped graph) by removing any sequence present in fewer than this many sequences. Set to 0 to disable. [default=2]") parser.add_argument("--gfa", nargs='*', default=None, help = "Produce a GFA for given graph type(s) if specified. Valid types are 'full', 'clip', and 'filter'. If no type specified 'clip' will be used ('full' used if clipping disabled). Multiple types can be provided separated by a space. [--gfa clip assumed by default]") @@ -254,6 +256,12 @@ def graphmap_join_defaults(options): setattr(options, name, value) return options +def graphmap_join_config_overrides(options, config_node): + """ apply the command-line options that override config attributes. + shared by cactus-graphmap-join and cactus-pangenome, which load the config separately """ + if getattr(options, 'clipFlank', None) is not None: + findRequiredNode(config_node, "graphmap_join").attrib["clipFlank"] = str(options.clipFlank) + def graphmap_join_validate_options(options): """ make sure the options make sense and fill in sensible defaults """ @@ -355,6 +363,15 @@ def graphmap_join_validate_options(options): if options.filter and not options.clip and not options.bypass: raise RuntimeError('--filter cannot be used without also disabling --clip.') + # 0 disables the flank trim; a negative would too, but only by falling through the same test, + # so say so rather than let it look like the negative clipFlankThreshold that means "calibrate" + if getattr(options, 'clipFlank', None) is not None: + if options.clipFlank < 0: + raise RuntimeError('--clipFlank cannot be negative (0 disables it).') + if options.clipFlank and not options.clip and not options.bypass: + raise RuntimeError('--clipFlank cannot be used with --clip disabled: it trims outward ' + 'from the intervals --clip removes.') + # check the reference name suffix check_sample_names(options.reference, options.reference[0]) @@ -633,6 +650,8 @@ def graphmap_join(options): if options.collapse: findRequiredNode(configNode, "graphmap").attrib["collapse"] = 'all' + graphmap_join_config_overrides(options, configNode) + # load up the hals hal_ids = [] for hal_path in options.hal: @@ -1273,6 +1292,12 @@ def clip_vg(job, options, config, vg_path, vg_id, phase): clip_vg_cmd += ['-u', str(options.clip)] if getOptionalAttrib(join_xml_node, "clipNonMinigraph", typeFn=bool, default=True): clip_vg_cmd += ['-a', graph_event] + # trim the tangled fringe left where an aligner extended anchors into a repeat + flank = getOptionalAttrib(join_xml_node, "clipFlank", typeFn=int, default=0) + if flank > 0: + clip_vg_cmd += ['-k', str(flank), + '-T', str(getOptionalAttrib(join_xml_node, "clipFlankThreshold", + typeFn=float, default=-1.0))] # disable reference cycle check if desiired if getOptionalAttrib(findRequiredNode(config.xmlRoot, "graphmap"), "collapse", typeFn=str, default="none") in ["all", "reference"]: @@ -1328,7 +1353,8 @@ def clip_vg(job, options, config, vg_path, vg_id, phase): if phase == 'full': cmd.append(['vg', 'ids', '-s', '-']) - cactus_call(parameters=cmd, outfile=clipped_path, job_memory=job.memory) + clip_stderr = cactus_call(parameters=cmd, outfile=clipped_path, job_memory=job.memory, + returnStdErr=True) # worth it cactus_call(parameters=['vg', 'validate', clipped_path]) @@ -1359,8 +1385,25 @@ def clip_vg(job, options, config, vg_path, vg_id, phase): stat_of.get('length', 'NA'))) # sample-stats and contig-stats are written by the exclusion report instead: they need the # input contig lengths, which are not available here + # -T calibrates against each graph, and the threshold it lands on decides how much + # sequence -k removes. It is reported on clip-vg's stderr, which cactus_call only surfaces + # when a command fails, so record it here or it is lost. The threshold gets its own column + # because that is the number worth reading down; the whole message is kept beside it so a + # rewording upstream costs the column rather than the record. + flank_stats_path = vg_path + '.flank-stats.tsv' + calib_line, threshold = 'flank trimming not run', 'NA' + for err_line in (clip_stderr or '').split('\n'): + if 'Flank calibration' in err_line: + calib_line = err_line.split(']:', 1)[-1].strip() + match = re.search(r'using threshold ([0-9.eE+-]+)', calib_line) + if match: + threshold = match.group(1) + break + with open(flank_stats_path, 'w') as flank_stats_file: + flank_stats_file.write('{}\t{}\t{}\n'.format(chr_name, threshold, calib_line)) out_stats = { 'path-stats.tsv' : job.fileStore.writeGlobalFile(path_stats_path), - 'graph-stats.tsv' : job.fileStore.writeGlobalFile(graph_stats_path) } + 'graph-stats.tsv' : job.fileStore.writeGlobalFile(graph_stats_path), + 'flank-stats.tsv' : job.fileStore.writeGlobalFile(flank_stats_path) } else: out_stats = None return job.fileStore.writeGlobalFile(clipped_path), out_stats @@ -2757,6 +2800,7 @@ def cat_bed_files(job, bed_ids): STATS_HEADERS = { 'path-stats.tsv': '#ref_chrom\tpath\tlength', 'graph-stats.tsv': '#ref_chrom\tnodes\tedges\tlength', + 'flank-stats.tsv': '#ref_chrom\tflank_threshold\tcalibration', } # path-stats has a row per surviving path fragment, which runs to millions on a filtered graph. diff --git a/src/cactus/refmap/cactus_pangenome.py b/src/cactus/refmap/cactus_pangenome.py index 8fd948392..e44029301 100644 --- a/src/cactus/refmap/cactus_pangenome.py +++ b/src/cactus/refmap/cactus_pangenome.py @@ -48,7 +48,7 @@ from cactus.refmap.cactus_graphmap import apply_mgsplit_filter_overrides, add_separate_ref_contigs_job from cactus.refmap.cactus_graphmap_split import graphmap_split_workflow, export_split_data from cactus.setup.cactus_align import make_batch_align_jobs, batch_align_jobs -from cactus.refmap.cactus_graphmap_join import graphmap_join_workflow, export_join_data, graphmap_join_options, graphmap_join_validate_options, vcflib_checks +from cactus.refmap.cactus_graphmap_join import graphmap_join_workflow, export_join_data, graphmap_join_options, graphmap_join_validate_options, graphmap_join_config_overrides, vcflib_checks from cactus.refmap.pangenome_exclusions import contig_sizes_job def pangenome_options(parser): @@ -223,6 +223,8 @@ def pangenome_config_overrides(options, config_node): if options.minIdentity is not None: findRequiredNode(config_node, "graphmap").attrib["minIdentity"] = str(options.minIdentity) + graphmap_join_config_overrides(options, config_node) + def main(): parser = Job.Runner.getDefaultArgumentParser()