Skip to content

Commit

Permalink
feat: deeptools bamcoverage (#1237)
Browse files Browse the repository at this point in the history
<!-- 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

Provide new wrapper for [DeepTools
bamcoverage](https://deeptools.readthedocs.io/en/develop/content/tools/bamCoverage.html).

### 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>
  • Loading branch information
8 people committed May 12, 2023
1 parent d70a806 commit ee2413b
Show file tree
Hide file tree
Showing 7 changed files with 159 additions and 1 deletion.
6 changes: 6 additions & 0 deletions bio/deeptools/bamcoverage/environment.yaml
@@ -0,0 +1,6 @@
channels:
- conda-forge
- bioconda
- nodefaults
dependencies:
- deeptools =3.5.1
17 changes: 17 additions & 0 deletions bio/deeptools/bamcoverage/meta.yaml
@@ -0,0 +1,17 @@
name: deepTools bamcoverage
url: https://deeptools.readthedocs.io/en/develop/content/tools/bamCoverage.html?highlight=bamcoverage
description: >
``deepTools bamcoverage`` takes an alignment of reads or fragments as input (BAM file) and generates a coverage track (bigWig or bedGraph) as output.
For more information about ``deepTools``, also see the `source code <https://github.com/deeptools/deepTools>`_.
authors:
- Thibault Dayris
input:
- bam: Path to alignment (BAM) file
- blacklist: Path to optional blacklist region file (BED)
output:
- Path to coverage file
params:
- effective_genome_size: Optional effective genome size value
- genome: Optional parameter used to fill effective genome size with pre-computed parameters. Can only be one of `GRCm37`, `GRCm38`, `GRCh37`, `GRCh38`, `dm3`, `dm6`, `WBcel235`, or `GRCz10`.
- read_length: Optional parameter used to fill effective genome size with pre-computed parameters. Can only be one of `50`, `75`, `100`, `150`, or `200`.
- extra: Optional parameters to be given to deepTools bamcoverage
15 changes: 15 additions & 0 deletions bio/deeptools/bamcoverage/test/Snakefile
@@ -0,0 +1,15 @@
rule test_deeptools_bamcoverage:
input:
bam="a.sorted.bam",
bai="a.sorted.bam.bai",
# Optional path to a blacklist bed file
# blacklist="",
output:
"a.coverage.bw",
params:
effective_genome_size=1000,
extra="",
log:
"logs/coverage.log",
wrapper:
"master/bio/deeptools/bamcoverage"
Binary file added bio/deeptools/bamcoverage/test/a.sorted.bam
Binary file not shown.
Binary file added bio/deeptools/bamcoverage/test/a.sorted.bam.bai
Binary file not shown.
107 changes: 107 additions & 0 deletions bio/deeptools/bamcoverage/wrapper.py
@@ -0,0 +1,107 @@
__author__ = "Thibault Dayris"
__copyright__ = "Copyright 2022, Thibault Dayris"
__email__ = "thibault.dayris@gustaveroussy.fr"
__license__ = "MIT"

from snakemake.shell import shell

log = snakemake.log_fmt_shell(stdout=True, stderr=True)
extra = snakemake.params.get("extra", "")

# See: https://deeptools.readthedocs.io/en/latest/content/feature/effectiveGenomeSize.html
default_effective_genome_size = {
"GRCz10": {
"50": 1195445591,
"75": 1251132686,
"100": 1280189044,
"150": 1312207169,
"200": 1321355241,
},
"WBcel235": {
"50": 95159452,
"75": 96945445,
"100": 98259998,
"150": 98721253,
"200": 98672758,
},
"dm3": {
"50": 130428560,
"75": 135004462,
"100": 139647232,
"150": 144307808,
"200": 148524010,
},
"dm6": {
"50": 125464728,
"75": 127324632,
"100": 129789873,
"150": 129941135,
"200": 132509163,
},
"GRCh37": {
"50": 2685511504,
"75": 2736124973,
"100": 2776919808,
"150": 2827437033,
"200": 2855464000,
},
"GRCh38": {
"50": 2701495761,
"75": 2747877777,
"100": 2805636331,
"150": 2862010578,
"200": 2887553303,
},
"GRCm37": {
"50": 2304947926,
"75": 2404646224,
"100": 2462481010,
"150": 2489384235,
"200": 2513019276,
},
"GRCm38": {
"50": 2308125349,
"75": 2407883318,
"100": 2467481108,
"150": 2494787188,
"200": 2520869189,
},
}

effective_genome_size = snakemake.params.get("effective_genome_size")
if not effective_genome_size:
genome = snakemake.params.get("genome")
read_length = snakemake.params.get("read_length")
if genome and read_length:
effective_genome_size = "--effectiveGenomeSize "
effective_genome_size += default_effective_genome_size[genome][str(read_length)]
else:
effective_genome_size = "--effectiveGenomeSize " + str(effective_genome_size)


output_format = ""
bigwig_format = ["bw", "bigwig"]
bedgraph_format = ["bg", "bedgraph"]
output_ext = str(snakemake.output[0]).split(".")[-1].lower()
if output_ext in bigwig_format:
output_format = "bigwig"
elif output_ext in bedgraph_format:
output_format = "bedgraph"
else:
raise ValueError("Output file should be either a bigwig or a bedgraph file")


blacklist = snakemake.input.get("blacklist", "")
if blacklist:
blacklist = "--blackListFileName " + blacklist

shell(
"bamCoverage "
"{blacklist} {extra} "
"--numberOfProcessors {snakemake.threads} "
"{effective_genome_size} "
"--bam {snakemake.input.bam} "
"--outFileName {snakemake.output} "
"--outFileFormat {output_format} "
"{log} "
)
15 changes: 14 additions & 1 deletion test.py
Expand Up @@ -2310,6 +2310,20 @@ def test_deeptools_computematrix():


@skip_if_not_modified
def test_deeptools_bamcoverage():
run(
"bio/deeptools/bamcoverage",
[
"snakemake",
"--cores",
"1",
"a.coverage.bw",
"--use-conda",
"-F",
],
)

@skip_if_not_modified
def test_deeptools_alignmentsieve():
run(
"bio/deeptools/alignmentsieve",
Expand All @@ -2323,7 +2337,6 @@ def test_deeptools_alignmentsieve():
],
)


@skip_if_not_modified
def test_deeptools_plotheatmap():
run(
Expand Down

0 comments on commit ee2413b

Please sign in to comment.