Have clip-vg trim flanks of big unaligned sequences

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

src/cactus/cactus_progressive_config.xml

@@ -457,6 +457,8 @@
 	<!-- gfaffix: toggle gfaffix normalization on/off -->
 	<!-- pathNormalize: toggle path normalization (merge equivalent paths through same site) on/off -->
 	<!-- clipNonMinigraph: clip out regions if the don't align to minigraph (if disabled, clip if they don't align to anything) -->
+	<!-- clipFlank: trim the fringe an aligner leaves where it extended anchors into a repeat, outward from each clipped interval by up to this many bp, for as long as unaligned sequence stays dense.  0 disables.  This value is the main control: inside a repeat the gate seldom closes on its own, so the cap decides how much goes.  Measured on a 124-haplotype chrY, 3.4% of edges at 20000, 9.4% at 50000, 29.4% at 200000.  50000 is the value benchmarked whole-genome (231/231 samples improved on F1); lowering it does not sharpen the trim, it just removes less -->
+	<!-- clipFlankThreshold: advanced.  keep trimming while more than this fraction of bases are unaligned, by the same test the clipping itself uses.  Negative calibrates per graph, which is almost always what you want: over 24 human chromosomes the calibrated value ran from 0.054 to 0.230, a 4.3x spread that tracks neither haplotype count nor chromosome size -->
 	<!-- minimizerOptions: options for vg minimizer -->
 	<!-- minFilterFragment: discard fragments that result from the vg clip frequency filter if they are smaller than this -->
 	<!-- clipContext: use this many context steps when doing edge clipping (activated via delEdgeFilter cli option) -->	
@@ -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"

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 <graphmap_join> 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.

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()