fix: bug with -O being used for gap opening penalty and samtools output directory

[huzuner]

Description

Gap penalty of minimap2 is done via "-O" option and the very same option also exists for samtools, but to indicate output directory. Using "-O" in the "extra" leads to the following error because of using snakemake utils for samtools (https://github.com/snakemake/snakemake-wrapper-utils/blob/f9cba17cb380dc7e6729a845ae37f26972c3d1dd/snakemake_wrapper_utils/samtools.py#LL76C23-L76C23):

You have specified output format (`-O/--output-fmt`) in params.extra; this is automatically inferred from output file extension.

The fastest and not over-engineered solution to me was to add the gap_opening variable to params, but I feel that it's also not very elegant especially when it comes to modifying the other penalty scores, but they can still be set in the extra.

QC

• I confirm that:

For all wrappers added by this PR,

• there is a test case which covers any introduced changes,
• input: and output: file paths in the resulting rule can be changed arbitrarily,
• either the wrapper can only use a single core, or the example rule contains a threads: x statement with x being a reasonable default,
• rule names in the test case are in snake_case and somehow tell what the rule is about or match the tools purpose or name (e.g., map_reads for a step that maps reads),
• all environment.yaml specifications follow the respective best practices,
• wherever possible, command line arguments are inferred and set automatically (e.g. based on file extensions in input: or output:),
• all fields of the example rules in the Snakefiles and their entries are explained via comments (input:/output:/params: etc.),
• stderr and/or stdout are logged correctly (log:), depending on the wrapped tool,
• temporary files are either written to a unique hidden folder in the working directory, or (better) stored where the Python function tempfile.gettempdir() points to (see here; this also means that using any Python tempfile default behavior works),
• the meta.yaml contains a link to the documentation of the respective tool or command,
• Snakefiles pass the linting (snakemake --lint),
• Snakefiles are formatted with snakefmt,
• Python wrapper scripts are formatted with black.
• Conda environments use a minimal amount of channels, in recommended ordering. E.g. for bioconda, use (conda-forge, bioconda, nodefaults, as conda-forge should have highest priority and defaults channels are usually not needed because most packages are in conda-forge nowadays).

[dlaehnemann]

Ah, looking at this in more detail, this makes sense. The get_samtools_opts is actually only used when samtools is used for sorting of the output. So this solution is probably the best and easiest.

One thing though. Could you document this further by:

• adding an entry for this new params field in the meta.yaml under notes::

  snakemake-wrappers/bio/minimap2/aligner/meta.yaml

  Lines 13 to 16 in 3cfaa10

  	notes: |
  	* The `extra` param allows for additional arguments for minimap2.
  	* The `sort` param allows to enable sorting (if output not PAF), and can be either 'none', 'queryname' or 'coordinate'.
  	* The `sort_extra` allows for extra arguments for samtools/picard

• adding this new params: gap_opening: "" specification to at least one of the example rules with sorted sam or bam output, so for example here:

  snakemake-wrappers/bio/minimap2/aligner/test/Snakefile

  Lines 60 to 63 in 3cfaa10

  	params:
  	extra="-x map-pb", # optional
  	sorting="coordinate", # optional: Enable sorting. Possible values: 'none', 'queryname' or 'coordinate'
  	sort_extra="", # optional: extra arguments for samtools/picard

[huzuner]

Thank you! How does it look now?

[dlaehnemann]

Looks good now, just one minor formatting suggestion -- feel free to include it or not.

[johanneskoester]

I would prefer to not have a special parameter. I think it is a bug that the wrapper utils complain here about -O, they should only check sort_extra.

[dlaehnemann]

But then you will have to change all the samtools wrappers, as get_samtools_opts() just parses params: extra:, no other params: entry. So this would maybe have to become a params: samtools_extra wherever get_samtools_opts() is used. Or we cannot use get_samtools_opts() here at all. Or we have to generalize the function to have another argument which determines which params: field to parse. Maybe this last option (with the default set to "extra") is the least breaking option.

[johanneskoester]

But then you will have to change all the samtools wrappers, as get_samtools_opts() just parses params: extra:, no other params: entry. So this would maybe have to become a params: samtools_extra wherever get_samtools_opts() is used. Or we cannot use get_samtools_opts() here at all. Or we have to generalize the function to have another argument which determines which params: field to parse. Maybe this last option (with the default set to "extra") is the least breaking option.

Exactly the latter was what I meant. Here is the PR: snakemake/snakemake-wrapper-utils#30

[johanneskoester]

snakemake-wrapper-utils 0.6 have been released with the fix. Hence, this PR should just add a keyword argument param_name="sort_extra" to the call of get_samtools_opts.

[fgvieira]

This seems to be the same issue described in #624.
Do you think you could also fix the bwa-mem, bwa sampe, bwa samse, and bwa sam(pe/se) wrappers?

PS - maybe bwa sampe and bwa samse should actually be deleted, since there is already bwa sam(pe/se)...

[huzuner]

The tests currently fail because the conda package for snakemake-wrappers-utils isn't updated yet.

[fgvieira]

Some other errors, related to memory.

[huzuner]

I think removing

        sort_extra="--write-index",  # Extra args for samtools/picard.

from Snakefile_samtools of bwa/mem would fix the

You have specified writing index (`--write-index`) in `params.extra`; this is automatically infered from `idx` output file.

error but I don't know about the memory error yet.

[huzuner]

We get the error:

Invalid maximum heap size: -Xmx0M
Error: Could not create the Java Virtual Machine.
Error: A fatal exception has occurred. Program will exit.

Regardless of the changes in the current PR, -Xmx0m returns 0 here related to a recent improvement in gathering resources in get_mem() in snakemake-wrapper-utils: https://github.com/snakemake/snakemake-wrapper-utils/blob/10ecf786c1bc3bcf6197eb8c6fb7f8615f0df2ee/snakemake_wrapper_utils/snakemake.py#L1

[dlaehnemann]

I just merged and released a non-zero default value setting for get_mem over at snakemake-wrapper-utils:
snakemake/snakemake-wrapper-utils#32

Once the new version 0.6.1 of snakemake-wrapper-utils is up on bioconda, pointing to it should fix this issue here.

[dlaehnemann]

Version 0.6.1 of snakemake-wrapper-utils is available via bioconda, now.

[dlaehnemann]

Ah, and I merged in the latest master branch, so you should probably pull the latest changes before further code changes locally on your machine.

[huzuner]

Thank you!
I now have updated the versions of snakemake-wrapper-utils for bwa/mem and bwa-mem2.