feat: fix threads and new IO options in Bowtie2 (#1324)

<!-- Ensure that the PR title follows conventional commit style (<type>:
<description>)-->
<!-- Possible types are here:
https://github.com/commitizen/conventional-commit-types/blob/master/index.json
-->

### Description

Fixes the following issues in the Bowtie2 wrapper:

* For `x` threads reserved in the snakemake rule, `2x` threads were
used. (Once in Bowtie2, once again in the Samtools pipe). This is now
fixed. For `x` threads reserved in the snakemake rule, `x-1` is used in
Bowtie2, and 1 in through the Samtools pipe.
* Optional metrics and alignment are now available through `output`
interface.
* Samtools compressed BAM/CRAM file does now explicitly keep header.
* Samtool indexing and compression is now documented in `meta.yaml` and
tested.
* Input file format is now automatically inferred when it is possible.

Open questions:

* `snakemake-wrappers-utils.samtools.get_samtools_opts()` tests
`params.extra` parameters. How do I let it test `params.extra_samtools`
? Then, we could let user set-up extra parameters for samtools view.
* Bowtie2 can concatenate multiple input files, as long as they are
separated by commas. This could be implemented in the wrapper. However,
when user inputs exactly two files in `input.sample`, then I can't find
how to guess if they need to be concatenated through `-U` option, or
considered as paired reads through `-1` and `-2`. My solutions always
include breaking retro-compatibility.

Both open questions can be ignored and the wrapper would still follow
its current behaviour. They are not breaking/blocking.

### QC
<!-- Make sure that you can tick the boxes below. -->

* [X] 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](https://en.wikipedia.org/wiki/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](https://stackoverflow.com/a/64594513/2352071),
* 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 `Snakefile`s 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](https://docs.python.org/3/library/tempfile.html#tempfile.gettempdir);
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,
* `Snakefile`s pass the linting (`snakemake --lint`),
* `Snakefile`s are formatted with
[snakefmt](https://github.com/snakemake/snakefmt),
* Python wrapper scripts are formatted with
[black](https://black.readthedocs.io).
* 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).

---------

Co-authored-by: tdayris <tdayris@gustaveroussy.fr>
Co-authored-by: tdayris <thibault.dayris@gustaveroussy.fr>
Co-authored-by: Johannes Köster <johannes.koester@uni-due.de>
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: snakedeploy-bot[bot] <115615832+snakedeploy-bot[bot]@users.noreply.github.com>
Co-authored-by: Felix Mölder <felix.moelder@uni-due.de>
Co-authored-by: Christopher Schröder <christopher.schroeder@tu-dortmund.de>

bio/bowtie2/align/meta.yaml

@@ -4,11 +4,22 @@ url: http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml
 authors:
   - Johannes Köster
   - Filipe G. Vieira
+  - Thibault Dayris
 input:
-  - FASTQ file(s)
-  - Bowtie2 indexed reference index
+  - sample: FASTQ file(s)
+  - idx: Bowtie2 indexed reference index
+  - ref: Optional path to genome sequence (FASTA)
+  - ref_fai: Optional path to reference genome sequence index (FAI)
 output:
-  - SAM/BAM/CRAM file
+  - SAM/BAM/CRAM file. This must be the first output file in the output file list.
+  - idx: Optional path to bam index.
+  - metrics: Optional path to metrics file.
+  - unaligned: Optional path to unaligned unpaired reads.
+  - unpaired: Optional path to unpaired reads that aligned at least once.
+  - unconcordant: Optional path to pairs that didn't align concordantly.
+  - concordant: Optional path to pairs that aligned concordantly at least once.
+params:
+  - extra: additional program arguments (except for `-x`, `-U`, `-1`, `-2`, `--interleaved`, `-b`, `--met-file`, `--un`, `--al`, `--un-conc`, `--al-conc`, `-f`, `--tab6`, `--tab5`, `-q`, or `-p/--threads`)
+  - interleaved: Input `sample` contains interleaved paired-end FASTQ/FASTA reads. `False`(default) or `True`.
 notes: |
-  * The `extra` param allows for additional program arguments.
   * This wrapper uses an inner pipe. Make sure to use at least two threads in your Snakefile.

bio/bowtie2/align/test/Snakefile

@@ -1,4 +1,4 @@
-rule bowtie2:
+rule test_bowtie2:
     input:
         sample=["reads/{sample}.1.fastq", "reads/{sample}.2.fastq"],
         idx=multiext(
@@ -10,8 +10,74 @@ rule bowtie2:
             ".rev.1.bt2",
             ".rev.2.bt2",
         ),
+        # ref="genome.fasta", #Required for CRAM output
     output:
         "mapped/{sample}.bam",
+        # idx="",
+        # metrics="",
+        # unaligned="",
+        # unpaired="",
+        # unconcordant="",
+        # concordant="",
+    log:
+        "logs/bowtie2/{sample}.log",
+    params:
+        extra="",  # optional parameters
+    threads: 8  # Use at least two threads
+    wrapper:
+        "master/bio/bowtie2/align"
+
+
+rule test_bowtie2_index:
+    input:
+        sample=["reads/{sample}.1.fastq", "reads/{sample}.2.fastq"],
+        idx=multiext(
+            "index/genome",
+            ".1.bt2",
+            ".2.bt2",
+            ".3.bt2",
+            ".4.bt2",
+            ".rev.1.bt2",
+            ".rev.2.bt2",
+        ),
+    output:
+        "mapped_idx/{sample}.bam",
+        idx="mapped_idx/{sample}.bam.bai",
+        metrics="mapped_idx/{sample}.metrics.txt",
+        unaligned="mapped_idx/{sample}.unaligned.sam",
+        unpaired="mapped_idx/{sample}.unpaired.sam",
+        # unconcordant="",
+        # concordant="",
+    log:
+        "logs/bowtie2/{sample}.log",
+    params:
+        extra="",  # optional parameters
+    threads: 8  # Use at least two threads
+    wrapper:
+        "master/bio/bowtie2/align"
+
+
+rule test_bowtie2_cram:
+    input:
+        sample=["reads/{sample}.1.fastq", "reads/{sample}.2.fastq"],
+        idx=multiext(
+            "index/genome",
+            ".1.bt2",
+            ".2.bt2",
+            ".3.bt2",
+            ".4.bt2",
+            ".rev.1.bt2",
+            ".rev.2.bt2",
+        ),
+        ref="genome.fasta",
+    output:
+        "mapped_idx/{sample}.cram",
+        # idx="",
+        # metrics="",
+        # unaligned="",
+        # unpaired="",
+        # unconcordant="",
+        # concordant="",
     log:
         "logs/bowtie2/{sample}.log",
     params:

bio/bowtie2/align/wrapper.py

@@ -9,7 +9,27 @@
 from snakemake_wrapper_utils.samtools import get_samtools_opts
 
 
-samtools_opts = get_samtools_opts(snakemake)
+def get_format(path: str) -> str:
+    """
+    Return file format since Bowtie2 reads files that
+    could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
+    """
+    if path.endswith((".gz", ".bz2")):
+        return path.split(".")[-2].lower()
+    return path.split(".")[-1].lower()
+
+
+bowtie2_threads = snakemake.threads - 1
+if bowtie2_threads < 1:
+    raise ValueError(
+        f"This wrapper expected at least two threads, got {snakemake.threads}"
+    )
+
+# Setting parse_threads to false since samtools performs only
+# bam compression. Thus the wrapper would use *twice* the amount
+# of threads reserved by user otherwise.
+samtools_opts = get_samtools_opts(snakemake, parse_threads=False)
+
 extra = snakemake.params.get("extra", "")
 log = snakemake.log_fmt_shell(stdout=True, stderr=True)
 
@@ -19,22 +39,62 @@
     n == 1 or n == 2
 ), "input->sample must have 1 (single-end) or 2 (paired-end) elements."
 
+reads = ""
 if n == 1:
-    reads = "-U {}".format(*snakemake.input.sample)
+    if get_format(snakemake.input.sample) in ("bam", "sam"):
+        reads = f"-b {snakemake.input.sample}"
+    else:
+        if snakemake.params.get("interleaved", False):
+            reads = f"--interleaved {snakemake.input.sample}"
+        else:
+            reads = f"-U {snakemake.input.sample}"
 else:
     reads = "-1 {} -2 {}".format(*snakemake.input.sample)
 
 
+if all(get_format(sample) in ("fastq", "fq") for sample in snakemake.input.sample):
+    extra += " -q "
+elif all(get_format(sample) == "tab5" for sample in snakemake.input.sample):
+    extra += " --tab5 "
+elif all(get_format(sample) == "tab6" for sample in snakemake.input.sample):
+    extra += " --tab6 "
+elif all(
+    get_format(sample) in ("fa", "mfa", "fasta") for sample in snakemake.input.sample
+):
+    extra += " -f "
+
+
+metrics = snakemake.output.get("metrics")
+if metrics:
+    extra += f" --met-file {metrics} "
+
+unaligned = snakemake.output.get("unaligned")
+if unaligned:
+    extra += f" --un {unaligned} "
+
+unpaired = snakemake.output.get("unpaired")
+if unpaired:
+    extra += f" --al {unpaired} "
+
+unconcordant = snakemake.output.get("unconcordant")
+if unconcordant:
+    extra += f" --un-conc {unconcordant} "
+
+concordant = snakemake.output.get("concordant")
+if concordant:
+    extra += f" --al-conc {concordant} "
+
+
 index = os.path.commonprefix(snakemake.input.idx).rstrip(".")
 
 
 shell(
     "(bowtie2"
-    " --threads {snakemake.threads}"
+    " --threads {bowtie2_threads}"
     " {reads} "
     " -x {index}"
     " {extra}"
-    "| samtools view"
+    "| samtools view --with-header "
     " {samtools_opts}"
     " -"
     ") {log}"

test.py

@@ -1560,7 +1560,17 @@ def test_blast_blastn():
 def test_bowtie2_align():
     run(
         "bio/bowtie2/align",
-        ["snakemake", "--cores", "1", "mapped/a.bam", "--use-conda", "-F"],
+        ["snakemake", "--cores", "2", "mapped_idx/a.cram", "--use-conda", "-F"],
+    )
+
+    run(
+        "bio/bowtie2/align",
+        ["snakemake", "--cores", "2", "mapped_idx/a.bam", "--use-conda", "-F"],
+    )
+
+    run(
+        "bio/bowtie2/align",
+        ["snakemake", "--cores", "2", "mapped/a.bam", "--use-conda", "-F"],
     )