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