diff --git a/.github/workflows/apptainer.yml b/.github/workflows/apptainer.yml index 7cf8cab..3542cf2 100644 --- a/.github/workflows/apptainer.yml +++ b/.github/workflows/apptainer.yml @@ -5,10 +5,11 @@ on: branches: - main - dev + workflow_dispatch: pull_request: branches-ignore: [] schedule: - - cron: '0 0 * * SUN' + - cron: '0 0 1 * *' jobs: build-linux: @@ -16,46 +17,79 @@ jobs: strategy: max-parallel: 5 matrix: - python: [3.8, '3.10'] + python: ['3.11'] fail-fast: false - steps: - - name: precleanup + # Clean up unnecessary preinstalled packages to free disk space + - name: Pre-cleanup run: | sudo rm -rf /usr/share/dotnet sudo rm -rf "$AGENT_TOOLSDIRECTORY" - - name: install graphviz + + # Cache APT .deb packages + - name: Cache APT archives + uses: actions/cache@v3 + with: + path: /var/cache/apt/archives + key: ${{ runner.os }}-apt-cache-v1 + + # Cache Apptainer installation + - name: Cache Apptainer install + id: cache-apptainer + uses: actions/cache@v3 + with: + path: | + /usr/bin/apptainer + /usr/lib/apptainer + /etc/apptainer + key: ${{ runner.os }}-apptainer-v1 + + # Install Apptainer only if not cached + - name: Install Apptainer + if: steps.cache-apptainer.outputs.cache-hit != 'true' run: | - sudo apt update - sudo apt-get install -y graphviz software-properties-common + sudo apt-get update + sudo apt-get install -y software-properties-common graphviz sudo add-apt-repository -y ppa:apptainer/ppa - sudo apt update - sudo apt install -y apptainer + sudo apt-get update + sudo apt-get install -y apptainer + + # Cache Apptainer image cache (~/.apptainer/cache) + - name: Cache Apptainer images + uses: actions/cache@v3 + with: + path: ~/.apptainer/cache + key: ${{ runner.os }}-apptainer-images-v1 - - name: checkout git repo - uses: actions/checkout@v3 + # Checkout repository + - name: Checkout repo + uses: actions/checkout@v4 - - name: Set up Python 3.X - uses: actions/setup-python@v3 + # Set up Python + - name: Set up Python ${{ matrix.python }} + uses: actions/setup-python@v5 with: python-version: ${{ matrix.python }} + # Install dependencies - name: Install dependencies run: | - pip install .[testing] + python -m pip install --upgrade pip + pip install . - - name: install package itself + # Install pinned pulp + - name: Install pinned pulp run: | - pip install . + pip install "pulp>=2.8" --no-deps - - name: testing + # Run tests using Apptainer + - name: Run Apptainer tests run: | + sequana_ribofinder --input-directory test/data \ + --reference-file test/data/Lepto.fa \ + --gff-file test/data/Lepto.gff \ + --use-apptainer --apptainer-prefix ~/.apptainer/cache - sequana_ribofinder --input-directory test/data --reference-file test/data/Lepto.fa --gff-file test/data/Lepto.gff --use-apptainer && cd ribofinder && sh ribofinder.sh - - - - - + cd ribofinder && bash ribofinder.sh diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 773e125..37aaf67 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -17,32 +17,33 @@ jobs: strategy: max-parallel: 5 matrix: - python: [3.8, 3.9, '3.10'] + python: ['3.11', '3.12'] fail-fast: false - steps: - name: install graphviz run: | + sudo apt-get update -qq sudo apt-get install -y graphviz - name: checkout git repo - uses: actions/checkout@v2 + uses: actions/checkout@v4 - name: conda/mamba - uses: mamba-org/provision-with-micromamba@main + uses: mamba-org/setup-micromamba@v2 with: environment-file: environment.yml - extra-specs: | + create-args: >- python=${{ matrix.python }} + cache-environment: true + cache-downloads: true - name: install package itself shell: bash -l {0} run: | pip install . - - name: Install dependencies shell: bash -l {0} run: | @@ -54,6 +55,7 @@ jobs: pytest -v --cov-report term-missing --cov=sequana_pipelines.ribofinder - name: coveralls + continue-on-error: true shell: bash -l {0} run: | echo $COVERALLS_REPO_TOKEN diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml index 435f810..5b7ca38 100644 --- a/.github/workflows/pypi.yml +++ b/.github/workflows/pypi.yml @@ -1,5 +1,5 @@ name: Publish to PyPI -on: +on: workflow_dispatch: push: tags: @@ -8,15 +8,17 @@ on: jobs: build-n-publish: name: Build and publish to PyPI and TestPyPI - runs-on: ubuntu-20.04 + runs-on: ubuntu-latest + permissions: + id-token: write steps: - - uses: actions/checkout@main - - name: Set up Python 3.8 - uses: actions/setup-python@v2 + - uses: actions/checkout@v4 + - name: Set up Python 3.11 + uses: actions/setup-python@v5 with: - python-version: 3.8 + python-version: 3.11 - - name: Install package + - name: Install package run: | pip install build poetry @@ -24,16 +26,15 @@ jobs: run: | rm -rf dist; poetry build - + - name: Publish distribution to Test PyPI + if: startsWith(github.ref, 'refs/tags') uses: pypa/gh-action-pypi-publish@release/v1 with: - user: __token__ password: ${{ secrets.TEST_PYPI_API_TOKEN }} - repository_url: https://test.pypi.org/legacy/ + repository-url: https://test.pypi.org/legacy/ - name: Publish distribution to PyPI if: startsWith(github.ref, 'refs/tags') uses: pypa/gh-action-pypi-publish@release/v1 with: - user: __token__ password: ${{ secrets.PYPI_API_TOKEN }} diff --git a/README.rst b/README.rst index 6c691f6..d63125b 100644 --- a/README.rst +++ b/README.rst @@ -9,9 +9,9 @@ .. image:: https://github.com/sequana/ribofinder/actions/workflows/main.yml/badge.svg :target: https://github.com/sequana/ribofinder/actions/workflows/main.yml -.. image:: https://img.shields.io/badge/python-3.8%20%7C%203.9%20%7C3.10-blue.svg +.. image:: https://img.shields.io/badge/python-3.9%20%7C%203.10%20%7C%203.11%20%7C%203.12-blue.svg :target: https://pypi.python.org/pypi/sequana - :alt: Python 3.8 | 3.9 | 3.10 + :alt: Python 3.9 | 3.10 | 3.11 | 3.12 @@ -67,12 +67,16 @@ Requirements This pipelines requires the following executable(s): -- bowtie1 >= 1.3.0 +- bowtie2 >= 2.4.0 +- bwa +- sambamba - bedtools - samtools -- bamtools - pigz +The aligner is selectable at the command line with ``--aligner bowtie2`` (default) +or ``--aligner bwa``. You only need to install the one(s) you intend to use. + .. image:: https://raw.githubusercontent.com/sequana/ribofinder/master/sequana_pipelines/ribofinder/dag.png Details @@ -109,6 +113,23 @@ Changelog ========= ==================================================================== Version Description ========= ==================================================================== +1.2.0 * Replace bowtie1 with bowtie2 (default) and bwa as aligner options + for rRNA mapping; the aligner is now user-selectable via + ``--aligner bowtie2|bwa``. + * Drop the bowtie1-specific ``fix_bowtie1_log`` workaround and the + standalone ``bam_indexing`` / ``samtools_faidx`` rules; the new + sequana-wrappers already produce sorted+indexed BAMs and a + faidx'ed reference. + * For bwa, feed multiqc with ``samtools stats`` output (bwa has no + native multiqc module); the mapping-rate plot in the HTML report + is built from ``multiqc_samtools_stats.txt``. + * For bowtie2, the HTML summary uses ``sequana.multiqc.plots.Bowtie2`` + which handles both single-end and paired-end outputs. + * Rewrite the HTML summary introduction to explain what the + pipeline does, how it works step-by-step, and how to interpret + each of the three plots (mapping rate, per-sequence proportions, + per-sequence RPKM). + (``bowtie2``, ``bwa``, ``sambamba`` replace ``bowtie``/``bamtools``). 1.1.1 * hotfix for running on HPC (slurm) 1.1.0 * Uses click (refactoring of sequana_pipetools) 1.0.1 * add sequana_wrappers in the config/pipeline diff --git a/environment.yml b/environment.yml index e30f4d3..1027dff 100644 --- a/environment.yml +++ b/environment.yml @@ -7,8 +7,10 @@ channels: - r dependencies: -- 'bowtie>=1.3.0' +- 'bowtie2>=2.4.0' +- bwa +- sambamba - samtools>1.7 -- bamtools - bedtools - pigz +- graphviz diff --git a/pyproject.toml b/pyproject.toml index 6e39ead..3bb6e58 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "poetry.core.masonry.api" [tool.poetry] name = "sequana-ribofinder" -version = "1.1.1" +version = "1.2.0" description = "A multi-sample identification of ribosomal content" authors = ["Sequana Team"] license = "BSD-3" @@ -19,9 +19,9 @@ classifiers = [ "Intended Audience :: Science/Research", "License :: OSI Approved :: BSD License", "Operating System :: POSIX :: Linux", - "Programming Language :: Python :: 3.8", - "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", "Topic :: Software Development :: Libraries :: Python Modules", "Topic :: Scientific/Engineering :: Bio-Informatics", "Topic :: Scientific/Engineering :: Information Analysis", @@ -33,11 +33,13 @@ packages = [ [tool.poetry.dependencies] -python = ">=3.8,<4.0" -sequana = ">=0.15.0" -sequana_pipetools = ">=0.16.0" +python = ">=3.9,<4.0" +sequana = ">=0.18.0" +sequana_pipetools = ">=1.5.4" +sequana-wrappers = ">=26.3.19" click-completion = "^0.5.2" - +snakemake = ">=7.32" +pulp = ">=2.8" [tool.poetry.scripts] sequana_ribofinder = "sequana_pipelines.ribofinder.main:main" @@ -45,8 +47,15 @@ sequana_ribofinder = "sequana_pipelines.ribofinder.main:main" [tool.poetry.group.dev.dependencies] black = "^23.7.0" -pytest = "^7.4.0" +pytest = "^8.0.0" mock = "^5.1.0" pytest-mock = "^3.11.1" pytest-cov = "^4.1.0" +[tool.pytest.ini_options] +filterwarnings = [ + "ignore::pytest.PytestDeprecationWarning", + "ignore::DeprecationWarning:click_completion", + "ignore::DeprecationWarning:importlib", +] +asyncio_mode = "auto" diff --git a/sequana_pipelines/ribofinder/config.yaml b/sequana_pipelines/ribofinder/config.yaml index 4c7abe8..9603d80 100644 --- a/sequana_pipelines/ribofinder/config.yaml +++ b/sequana_pipelines/ribofinder/config.yaml @@ -1,5 +1,5 @@ # ============================================================================ -# Config file for RNA-seq +# Config file for Ribofinder # # ==================[ Sections for the users ]================================ # @@ -7,7 +7,6 @@ # If input_directory provided, use it otherwise if input_pattern provided, # use it, otherwise use input_samples. # ============================================================================ -sequana_wrappers: "v23.11.18" input_directory: input_readtag: _R[12]_ @@ -19,23 +18,16 @@ input_pattern: '*fastq.gz' # # :Parameters: # -# - aligner: bowtie (no other mapper for now) -# - rRNA_file: path to an existing fasta file for ribosomal RNA (to be found in -# genome_directory) -# - rRNA_feature: if rRNA_file is not provided, ribosomal RNA will be extract -# from GFF using this feature name. It must be found. -# - indexing: if indexing is set to True, the index for bowtie1 will be done as -# well as the indexing of the aligner provided. If the files exists already, not -# indexing is performed. If you want to force the index building despite the -# presence of the index files, then, use the force_indexing parameter and set -# it to True. Indexing is followed by force_indexing to make sure we do not -# erase the index files, which may be large. +# - aligner: bowtie2 or bwa +# - rRNA_file: path to an existing fasta file for ribosomal RNA +# - rRNA_feature: if rRNA_file is not provided, ribosomal RNA will be extracted +# from GFF using this feature name. It must be found. # # Only one annotation file must be provided. If so, we extract the feature -# rRNA from the annotation, identifiy start/end and create fasta file on the +# rRNA from the annotation, identify start/end and create fasta file on the # fly. You may simply provide an input fasta file with rRNA if you have one. general: - aligner: bowtie1 + aligner: bowtie2 rRNA_file: '' rRNA_feature: rRNA genbank_file: '' @@ -43,28 +35,39 @@ general: reference_file: apptainers: - sequana_ribofinder: "https://zenodo.org/record/7348115/files/sequana_ribofinder_0.12.0.img" + sequana_tools: "https://zenodo.org/record/18257162/files/sequana_tools_26.1.14.img" pigz: "https://zenodo.org/record/7346805/files/pigz_2.4.0.img" - bedtools: "https://zenodo.org/record/7346774/files/bedtools_2.30.0.img" - samtools: "https://zenodo.org/record/7215278/files/samtools_1.15.0.img" - graphviz: "https://zenodo.org/record/7928262/files/graphviz_7.0.5.img" + graphviz: "https://zenodo.org/record/7928262/files/graphviz_7.0.5.img" multiqc: "https://zenodo.org/record/10205070/files/multiqc_1.16.0.img" ############################################################################# -# bowtie1_mapping_rna used to align reads against ribosomal RNA +# bowtie2_mapping_rna used to align reads against ribosomal RNA # # :Parameters: # -# - do: if unchecked, this rule is ignored -# - options: any options recognised by bowtie1 tool +# - options: any options recognised by bowtie2 # - threads: number of threads to be used # -bowtie1_mapping_rna: - do: true +bowtie2_mapping_rna: options: '' threads: 4 +############################################################################# +# bwa_mapping_rna used to align reads against ribosomal RNA +# +# :Parameters: +# +# - index_algorithm: 'is' (default, small genomes) or 'bwtsw' (large genomes) +# - options: any options recognised by bwa mem +# - threads: number of threads to be used +# - tmp_directory: scratch directory for sambamba sort +# +bwa_mapping_rna: + index_algorithm: is + options: '' + threads: 4 + tmp_directory: ./tmp ############################################################################# @@ -74,13 +77,12 @@ bowtie1_mapping_rna: # :Parameters: # # - options: any options recognised by multiqc -# - output-directory: Create report in the specified output directory -# - config_file: by default, we use sequana RNA-seq multiqc_config file. -# If you want your own multiqc, fill this entry +# - modules: space-separated multiqc modules to enable +# ('bowtie2' for bowtie2 aligner, 'samtools' for bwa aligner) +# - config_file: multiqc config file multiqc: options: -f input_directory: . - modules: "" + modules: bowtie2 config_file: multiqc_config.yaml - diff --git a/sequana_pipelines/ribofinder/main.py b/sequana_pipelines/ribofinder/main.py index 4f19366..dac1bee 100755 --- a/sequana_pipelines/ribofinder/main.py +++ b/sequana_pipelines/ribofinder/main.py @@ -48,7 +48,11 @@ @include_options_from(ClickInputOptions) @include_options_from(ClickGeneralOptions) @click.option( - "--aligner", "aligner", type=click.Choice(["bowtie1"]), default="bowtie1", help="the alignement tool (bowtie1)" + "--aligner", + "aligner", + type=click.Choice(["bowtie2", "bwa"]), + default="bowtie2", + help="the alignment tool (bowtie2 or bwa). Default: bowtie2.", ) @click.option( "--rRNA-feature", @@ -87,6 +91,12 @@ def main(**options): # --------------------------------------------------------- general def fill_aligner(): cfg.general.aligner = options.aligner + # Keep multiqc modules in sync with the selected aligner. bwa has no + # dedicated multiqc module, so we rely on samtools stats instead. + if options.aligner == "bwa": + cfg.multiqc.modules = "samtools" + else: + cfg.multiqc.modules = "bowtie2" def fill_rRNA_feature(): cfg.general.rRNA_feature = options.rRNA_feature @@ -107,7 +117,7 @@ def fill_reference_file(): if options.reference_file: cfg.general.reference_file = os.path.abspath(options.reference_file) - if options["from_project"]: + if options.from_project: if "--aligner" in sys.argv: fill_aligner() if "--rRNA-feature" in sys.argv: @@ -118,7 +128,7 @@ def fill_reference_file(): fill_genbank_file() if "--gff-file" in sys.argv: fill_gff_file() - if "--reference-fle" in sys.argv: + if "--reference-file" in sys.argv: fill_reference_file() else: @@ -141,7 +151,7 @@ def fill_reference_file(): gbk = GenBank(options.genbank_file) if options.genbank_file is None and options.gff_file is None: click.echo("Most probably you want to provide an annotation (genbank-file or gff-file)", err=True) - sys.exit(0) + sys.exit(1) # finalise the command and save it; copy the snakemake. update the config # file and save it. diff --git a/sequana_pipelines/ribofinder/ribofinder.rules b/sequana_pipelines/ribofinder/ribofinder.rules index 170ba19..886cb7a 100644 --- a/sequana_pipelines/ribofinder/ribofinder.rules +++ b/sequana_pipelines/ribofinder/ribofinder.rules @@ -26,6 +26,13 @@ configfile: "config.yaml" manager = PipelineManager("ribofinder", config) +aligner = config["general"]["aligner"] +if aligner not in ("bowtie2", "bwa"): + raise ValueError( + f"config['general']['aligner'] must be 'bowtie2' or 'bwa', got {aligner!r}" + ) + + rule pipeline: input: ".sequana/rulegraph.svg", @@ -89,123 +96,152 @@ elif get_annot_file().endswith(".gff") or get_annot_file().endswith(".gff3"): output: fasta = "build_feature_fasta/feature.fasta", gff = "build_feature_fasta/feature_rRNA.gff" + log: + "logs/build_feature_fasta/build_feature_fasta.log" params: feature = config['general']['rRNA_feature'] container: - config["apptainers"]["bedtools"] + config["apptainers"]["sequana_tools"] shell: """ # used to be gawk but awk is more generic and available in the container. - awk '{{ if ($3=="{params.feature}") print }}' {input.gff} > {output.gff} + awk '{{ if ($3=="{params.feature}") print }}' {input.gff} > {output.gff} 2> {log} if [ -s {output.gff} ] then - bedtools getfasta -fi {input.fasta} -bed {output.gff} -fo {output.fasta} + bedtools getfasta -fi {input.fasta} -bed {output.gff} -fo {output.fasta} >> {log} 2>&1 else : echo -e ">empty\\nAAAAAAAAAAAAAA" > {output.fasta} fi """ - rule samtools_faidx: - input: - fasta = "build_feature_fasta/feature.fasta" - output: - fai = "build_feature_fasta/feature.fasta.fai" - container: - config["apptainers"]["samtools"] - shell: - """ - samtools faidx {input.fasta} - """ - - -rule bowtie1_indexing: - input: - reference= "build_feature_fasta/feature.fasta" - output: - "build_feature_fasta/feature_rRNA.1.ebwt" - log: - "indexing/bowtie_rRNA.log" - params: - options="" - threads: - 2 - container: - config['apptainers']['sequana_ribofinder'] - wrapper: - f"{manager.wrappers}/wrappers/bowtie1/build" + # The bowtie2/build and bwa/build wrappers both call `samtools faidx` on + # the reference, so a dedicated samtools_faidx rule is no longer needed. +# ====================================================== Indexing and mapping +# +# Reads are decompressed, then mapped onto the rRNA feature fasta using the +# aligner selected by the user (bowtie2 or bwa). Both wrappers also index the +# sorted BAM file, so no separate bam_indexing rule is required. -"""With paired data, alignement on rRNA leads to 0% alignment if we use R1 and +"""With paired data, alignment on rRNA leads to 0% alignment if we use R1 and R2. If we use R1 only, the percentage is >0. First reason is that reads are not -trimmed properly. In truth, bowtie2 supports local alignments which means it can -soft-clip non-matching (=adapter) content while still align the local part of -the read that matches the reference. With Bowtie1 the read will probably go -unaligned due to the many mismatches. So we do not include R2. +trimmed properly. bowtie2 supports local alignments which means it can +soft-clip non-matching (=adapter) content while still aligning the local part +of the read that matches the reference. For this reason we only feed R1 to the +rRNA mapping step. """ rule unpigz: input: manager.getrawdata() output: temp("{sample}/unpigz/{sample}.fastq") + log: + "{sample}/unpigz/unpigz.log" threads: 4 container: config['apptainers']['pigz'] shell: """ - unpigz -p {threads} -fk --stdout {input[0]} > {output} + unpigz -p {threads} -fk --stdout {input[0]} > {output} 2> {log} """ -# ========================================================== Bowtie1 indexing -# -rule bowtie1_mapping_rna: - input: - fastq="{sample}/unpigz/{sample}.fastq", - index="build_feature_fasta/feature_rRNA.1.ebwt" - output: - bam="{sample}/bowtie1/{sample}.bam", - sorted="{sample}/bowtie1/{sample}.sorted.bam" - log: - "{sample}/bowtie1/{sample}.log" - params: - options="" - threads: - config['bowtie1_mapping_rna']['threads'] - container: - config['apptainers']['sequana_ribofinder'] - wrapper: - f"{manager.wrappers}/wrappers/bowtie1/align" +if aligner == "bowtie2": + + rule bowtie2_indexing: + input: + reference="build_feature_fasta/feature.fasta" + output: + multiext( + "build_feature_fasta/feature", + ".1.bt2", ".2.bt2", ".3.bt2", ".4.bt2", ".rev.1.bt2", ".rev.2.bt2", + ) + log: + "indexing/bowtie2_rRNA.log" + params: + options="" + threads: 2 + container: + config['apptainers']['sequana_tools'] + shell: + manager.get_shell("bowtie2/build", "v1") -# ======================================================== Fix bowtie log -# This fix_bowtie1_log is used to fix a bug in multiqc showing 0% match in all samples. -# Note: same rule as in rnaseq pipeline. -rule fix_bowtie1_log: - input: - expand("{sample}/bowtie1/{sample}.log", sample=manager.samples) - output: - "logs/fix_bowtie1_log/done" - run: - for filename in input: - # we read the file - with open(filename) as fin: - data = fin.readlines() - # we update the file - with open(filename, "w") as fout: - for line in data: - if "least one alignment" in line: - fout.write(line) - fout.write(line.replace("least one alignment", "least one reported alignment")) - else: - fout.write(line) - with open(output[0], "w") as fout: - fout.write("") + rule bowtie2_mapping_rna: + input: + fastq=["{sample}/unpigz/{sample}.fastq"], + idx=multiext( + "build_feature_fasta/feature", + ".1.bt2", ".2.bt2", ".3.bt2", ".4.bt2", ".rev.1.bt2", ".rev.2.bt2", + ) + output: + bam="{sample}/bowtie2/{sample}.sorted.bam" + log: + "{sample}/bowtie2/{sample}.log" + params: + options=config["bowtie2_mapping_rna"]["options"] + threads: + config["bowtie2_mapping_rna"]["threads"] + container: + config['apptainers']['sequana_tools'] + shell: + manager.get_shell("bowtie2/align", "v1") + + + # Alias variables used downstream by plotting and multiqc rules. + mapper_bam_pattern = "{sample}/bowtie2/{sample}.sorted.bam" + mapper_log_pattern = "{sample}/bowtie2/{sample}.log" + + +elif aligner == "bwa": + + rule bwa_indexing: + input: + reference="build_feature_fasta/feature.fasta" + output: + bwt="build_feature_fasta/feature.fasta.bwt", + fai="build_feature_fasta/feature.fasta.fai" + log: + "indexing/bwa_rRNA.log" + params: + options="", + index_algorithm=config["bwa_mapping_rna"]["index_algorithm"] + threads: 2 + container: + config['apptainers']['sequana_tools'] + shell: + manager.get_shell("bwa/build", "v1") + + + rule bwa_mapping_rna: + input: + fastq="{sample}/unpigz/{sample}.fastq", + reference="build_feature_fasta/feature.fasta", + bwt="build_feature_fasta/feature.fasta.bwt", + fai="build_feature_fasta/feature.fasta.fai" + output: + sorted="{sample}/bwa/{sample}.sorted.bam" + log: + "{sample}/bwa/{sample}.log" + params: + options=config["bwa_mapping_rna"]["options"], + tmp_directory=config["bwa_mapping_rna"]["tmp_directory"] + threads: + config["bwa_mapping_rna"]["threads"] + container: + config['apptainers']['sequana_tools'] + shell: + manager.get_shell("bwa/align", "v1") + + + mapper_bam_pattern = "{sample}/bwa/{sample}.sorted.bam" + mapper_log_pattern = "{sample}/bwa/{sample}.log" rule plotting: input: - bam_files=expand("{sample}/bowtie1/{sample}.sorted.bam", sample=manager.samples), - fasta_file= "build_feature_fasta/feature.fasta" + bam_files=expand(mapper_bam_pattern, sample=manager.samples), + fasta_file="build_feature_fasta/feature.fasta" output: png="outputs/proportions.png", rpkm="outputs/RPKM.png" @@ -227,12 +263,12 @@ rule plotting: results.append(c) df = pd.DataFrame(results) - # a copy to preserve correcte sequence name for later + # a copy to preserve correct sequence name for later df1 = df.copy() df1.columns = [x[0:30] for x in df1.columns ] - # first plot shown percentage of hits on each sequence + # first plot shown percentage of hits on each sequence import seaborn as sns sns.violinplot(data=df1, orient='h') xlabel("Ribosomal content (%)", fontsize=16) @@ -240,7 +276,7 @@ rule plotting: tight_layout() savefig(output['png']) - # now, we shown the RPKM (read count normalised by gene length) + # now, we show the RPKM (read count normalised by gene length) f = FastA(input['fasta_file']) df = df.T @@ -254,26 +290,34 @@ rule plotting: savefig(output['rpkm']) - -# ========================================================== BAM indexing +# ========================================== samtools stats (for bwa multiqc) # -rule bam_indexing: - input: - "{sample}/bowtie1/{sample}.sorted.bam" - output: - "{sample}/bowtie1/{sample}.sorted.bam.bai" - container: - config['apptainers']['sequana_ribofinder'] - shell: - """ - bamtools index -in {input} - """ +# bwa itself produces no structured stats usable by multiqc. To expose +# mapped/unmapped percentages in both multiqc and the HTML summary, we feed +# multiqc the output of `samtools stats` on the sorted BAM file. +if aligner == "bwa": + rule samtools_stats: + input: + bam=mapper_bam_pattern + output: + "{sample}/samtools_stats/{sample}.stats.txt" + log: + "{sample}/samtools_stats/{sample}.log" + container: + config['apptainers']['sequana_tools'] + shell: + """ + samtools stats {input.bam} > {output} 2> {log} + """ # ========================================================== multiqc # -sequana_multiqc_input = expand("{sample}/bowtie1/{sample}.sorted.bam.bai", sample=sorted(manager.samples)) -sequana_multiqc_input += ["logs/fix_bowtie1_log/done"] +sequana_multiqc_input = expand(mapper_log_pattern, sample=sorted(manager.samples)) +if aligner == "bwa": + sequana_multiqc_input += expand( + "{sample}/samtools_stats/{sample}.stats.txt", sample=sorted(manager.samples) + ) rule multiqc: @@ -290,20 +334,20 @@ rule multiqc: "multiqc/multiqc.log" container: config["apptainers"]["multiqc"] - wrapper: - f"{manager.wrappers}/wrappers/multiqc" + shell: + manager.get_shell("multiqc/run", "v1") # ========================================================== rulegraph rule rulegraph: input: str(manager.snakefile) output: - svg = "rulegraph/rulegraph.dot" + temp("rulegraph/rulegraph.dot") params: mapper = {"multiqc": "../multiqc/multiqc_report.html"}, configname = "config.yaml" - wrapper: - f"{manager.wrappers}/wrappers/rulegraph" + run: + manager.get_run("rulegraph/run", "v1")(input, output, params) rule dot2svg: @@ -314,13 +358,13 @@ rule dot2svg: container: config['apptainers']['graphviz'] shell: - """dot -Tsvg {input} -o {output}""" + manager.get_shell("graphviz/dot2svg", "v1") # Those rules takes a couple of seconds so no need for a cluster -localrules: rulegraph +localrules: rulegraph, dot2svg onsuccess: @@ -338,18 +382,128 @@ onsuccess: "pipeline_version": ribofinder.version } - HTML = """

Summary

Please see the ribosomal rates here below - (a similar plot from multiqc is here )

""" + aligner = config["general"]["aligner"] + + HTML = f""" +

Overview

+

+The ribofinder pipeline estimates the proportion of ribosomal RNA +(rRNA) contamination in a set of sequencing samples. This information is +useful as a QC step for RNA-seq experiments: after rRNA depletion (e.g. +Ribo-Zero) or poly-A selection, the residual rRNA content is expected to +be low; a high value typically indicates that the library preparation +failed or that the depletion protocol did not work as intended. +

+

How it works

+
    +
  1. A reference fasta of rRNA sequences is prepared: either provided + directly via --rRNA-file, or extracted on-the-fly from a + GFF/GenBank annotation by selecting the + {config['general']['rRNA_feature']} feature.
    2. +
  2. Each sample is decompressed (only R1 is kept for paired-end data, + to avoid low mapping rates on short, adapter-contaminated mates).
    3. +
  3. Reads are mapped onto the rRNA reference with + {aligner}, producing a sorted and indexed BAM file per + sample.
    4. +
  4. Mapping statistics and per-sequence hit counts are aggregated in a + MultiQC report and in the + summary plots below.
    5. +
+ +

Mapping rate: how much of each sample is rRNA?

+

+This first plot shows, for every sample, the fraction of reads that +aligned to the rRNA reference. Higher values mean more ribosomal +contamination. As a rule of thumb, values < 10 % are expected +after a successful rRNA depletion, whereas values > 20 % suggest that +the depletion did not work well and that a large fraction of the +sequencing budget was spent on ribosomal reads. The equivalent plot from +multiqc is also available in the +multiqc report. +

+""" - from sequana.multiqc.plots import Bowtie1Reader - if os.path.exists("multiqc/multiqc_report_data/multiqc_bowtie1.txt"): - br = Bowtie1Reader("multiqc/multiqc_report_data/multiqc_bowtie1.txt") - fig = br.plot_bar(html_code=True) - from plotly import offline - HTML += offline.plot(fig, output_type="div", include_plotlyjs=True) + from plotly import offline - HTML += '

Here below, you can see how the ribosomal hits are distributed amongst the sequences (note that this is not normalised against the lenth of each sequence)

' - HTML += '

Here below, you can see how the reads are distributed amongst the sequences normalised by their lengths (RPKM=reads per kilo base)

' + try: + if aligner == "bowtie2": + from sequana.multiqc.plots import Bowtie2 + + # multiqc writes either a combined table or a per-mode plot table + # depending on the version and on whether samples are SE or PE. + candidates = [ + "multiqc/multiqc_data/mqc_bowtie2_se_plot_1.txt", + "multiqc/multiqc_data/mqc_bowtie2_pe_plot_1.txt", + "multiqc/multiqc_data/multiqc_bowtie2.txt", + "multiqc/multiqc_report_data/mqc_bowtie2_se_plot_1.txt", + "multiqc/multiqc_report_data/mqc_bowtie2_pe_plot_1.txt", + "multiqc/multiqc_report_data/multiqc_bowtie2.txt", + ] + filename = next((f for f in candidates if os.path.exists(f)), None) + if filename: + br = Bowtie2(filename) + fig = br.plot(html_code=True) + HTML += offline.plot(fig, output_type="div", include_plotlyjs=True) + elif aligner == "bwa": + # bwa has no dedicated multiqc module, so build the mapping-rate + # bar plot from the samtools_stats multiqc output. + import pandas as pd + import plotly.graph_objects as go + + candidates = [ + "multiqc/multiqc_data/multiqc_samtools_stats.txt", + "multiqc/multiqc_report_data/multiqc_samtools_stats.txt", + ] + filename = next((f for f in candidates if os.path.exists(f)), None) + if filename: + df = pd.read_csv(filename, sep="\t") + mapped_pct = df["reads_mapped"] / df["raw_total_sequences"] * 100 + unmapped_pct = 100 - mapped_pct + fig = go.Figure() + fig.add_trace(go.Bar( + y=df["Sample"], x=mapped_pct, name="Aligned (rRNA)", + marker_color="#389f1f", orientation="h", + )) + fig.add_trace(go.Bar( + y=df["Sample"], x=unmapped_pct, name="Not aligned", + marker_color="#120946", orientation="h", + )) + fig.update_layout( + barmode="stack", xaxis_tickangle=-45, height=400, + title="Mapping on ribosomal/contaminant reference", + xaxis_title="Percentage of reads", + ) + HTML += offline.plot(fig, output_type="div", include_plotlyjs=True) + except Exception as err: + print(f"WARNING: could not add mapping rate plot to the HTML report: {err}") + + HTML += """ +

Per-sequence distribution of ribosomal hits

+

+The rRNA reference may contain several sequences (e.g. 5S, 16S, 23S for +bacteria, or 18S/28S for eukaryotes). The violin plot below shows, for +every rRNA sequence (y-axis), the distribution across all samples of the +percentage of total reads that aligned to it (x-axis). Each violin +therefore summarises how much that specific ribosomal sequence +contributes to the overall rRNA signal. This is a raw count and is +not normalised against the length of each sequence, so longer +sequences naturally tend to collect more reads. +

+ + +

Per-sequence RPKM distribution

+

+The same data is shown below, but normalised by the length of each +reference sequence (RPKM = Reads Per Kilobase). This makes the +contributions of short and long rRNA sequences directly comparable: for +instance, a short 5S rRNA and a long 23S rRNA that are equally abundant +in the sample will produce similar RPKM values, while the raw count plot +above would make the 23S look dominant. Use this plot to identify which +specific rRNA sequences are over- or under-represented relative to what +you would expect. +

+ +""" # Now the final report. add the original command in the HTML report try: diff --git a/sequana_pipelines/ribofinder/schema.yaml b/sequana_pipelines/ribofinder/schema.yaml index fb2ebb7..7ec91bb 100644 --- a/sequana_pipelines/ribofinder/schema.yaml +++ b/sequana_pipelines/ribofinder/schema.yaml @@ -3,9 +3,6 @@ type: map mapping: - "sequana_wrappers": - type: str - required: true "input_directory": type: str required: True @@ -15,6 +12,9 @@ mapping: "input_pattern": type: str required: True + "exclude_pattern": + type: str + required: False "apptainers": type: any @@ -24,7 +24,7 @@ mapping: "aligner": type: str required: True - enum: ["bowtie1"] + enum: ["bowtie2", "bwa"] "rRNA_feature": type: str "rRNA_file": @@ -36,7 +36,7 @@ mapping: "gff_file": type: str - 'bowtie1_mapping_ref': + 'bowtie2_mapping_rna': type: map mapping: "options": @@ -44,20 +44,22 @@ mapping: "threads": type: int required: True - range: { min: 1 } - + range: { min: 1} - 'bowtie1_mapping_rna': + 'bwa_mapping_rna': type: map mapping: - "do": - type: bool + "index_algorithm": + type: str + enum: ["is", "bwtsw", "rb2"] "options": type: str "threads": type: int required: True range: { min: 1} + "tmp_directory": + type: str "multiqc": type: map @@ -74,6 +76,3 @@ mapping: type: str "modules": type: str - - - diff --git a/sequana_pipelines/ribofinder/tools.txt b/sequana_pipelines/ribofinder/tools.txt index e5c02cc..8eb13b3 100644 --- a/sequana_pipelines/ribofinder/tools.txt +++ b/sequana_pipelines/ribofinder/tools.txt @@ -1,5 +1,8 @@ -bowtie +bowtie2 +bwa +sambamba samtools bedtools -bamtools pigz +graphviz +multiqc diff --git a/test/test_main.py b/test/test_main.py index f2b3ffa..a560ee9 100644 --- a/test/test_main.py +++ b/test/test_main.py @@ -1,60 +1,229 @@ -import easydev import os -import tempfile import subprocess -import sys +import tempfile + +import pytest from click.testing import CliRunner -from sequana_pipelines.ribofinder.main import main +from sequana_pipelines.ribofinder.main import main from . import test_dir sharedir = f"{test_dir}/data" + +# ----------------------------------------------------------------------------- +# Setup-only tests (no snakemake execution). Fast and safe to run on any host +# because they don't need bowtie2/bwa/samtools to be available on PATH. +# ----------------------------------------------------------------------------- + + def test_standalone_subprocess(): - directory = tempfile.TemporaryDirectory() - cmd = """sequana_ribofinder --input-directory {} - --working-directory {} --force""".format(sharedir, directory.name) - subprocess.call(cmd.split()) + with tempfile.TemporaryDirectory() as directory: + cmd = f"sequana_ribofinder --input-directory {sharedir} " + cmd += f"--working-directory {directory} --force " + cmd += f"--rRNA-file {sharedir}/feature.fasta" + subprocess.call(cmd.split()) def test_standalone_script(): - directory = tempfile.TemporaryDirectory() + with tempfile.TemporaryDirectory() as directory: + runner = CliRunner() + results = runner.invoke( + main, + [ + "--input-directory", sharedir, + "--working-directory", directory, + "--force", + "--rRNA-file", f"{sharedir}/feature.fasta", + ], + ) + assert results.exit_code == 0 + + +def test_standalone_script_bowtie2_explicit(): + """Explicit --aligner bowtie2 should write bowtie2 in the generated config.""" + with tempfile.TemporaryDirectory() as directory: + runner = CliRunner() + results = runner.invoke( + main, + [ + "--input-directory", sharedir, + "--working-directory", directory, + "--force", + "--rRNA-file", f"{sharedir}/feature.fasta", + "--aligner", "bowtie2", + ], + ) + assert results.exit_code == 0 + with open(os.path.join(directory, "config.yaml")) as fin: + content = fin.read() + assert "aligner: bowtie2" in content + # multiqc module should auto-track the aligner selection + assert "modules: bowtie2" in content + + +def test_standalone_script_bwa(): + """--aligner bwa must switch both the aligner and the multiqc modules.""" + with tempfile.TemporaryDirectory() as directory: + runner = CliRunner() + results = runner.invoke( + main, + [ + "--input-directory", sharedir, + "--working-directory", directory, + "--force", + "--rRNA-file", f"{sharedir}/feature.fasta", + "--aligner", "bwa", + ], + ) + assert results.exit_code == 0 + with open(os.path.join(directory, "config.yaml")) as fin: + content = fin.read() + assert "aligner: bwa" in content + assert "modules: samtools" in content + + +def test_standalone_script_invalid_aligner(): + """click should reject any aligner outside the declared choice.""" + with tempfile.TemporaryDirectory() as directory: + runner = CliRunner() + results = runner.invoke( + main, + [ + "--input-directory", sharedir, + "--working-directory", directory, + "--force", + "--rRNA-file", f"{sharedir}/feature.fasta", + "--aligner", "bowtie1", + ], + ) + assert results.exit_code != 0 + + +def test_standalone_script_gff(): + """Setup pipeline from reference fasta + GFF annotation.""" + with tempfile.TemporaryDirectory() as directory: + runner = CliRunner() + results = runner.invoke( + main, + [ + "--input-directory", sharedir, + "--working-directory", directory, + "--force", + "--reference-file", f"{sharedir}/Lepto.fa", + "--gff-file", f"{sharedir}/Lepto.gff", + ], + ) + assert results.exit_code == 0 + + +def test_standalone_script_missing_inputs(): + """Neither --rRNA-file nor --reference-file should fail with a clear exit.""" + with tempfile.TemporaryDirectory() as directory: + runner = CliRunner() + results = runner.invoke( + main, + [ + "--input-directory", sharedir, + "--working-directory", directory, + "--force", + ], + ) + assert results.exit_code != 0 + + +def test_version(): + cmd = "sequana_ribofinder --version" + subprocess.call(cmd.split()) + + +# ----------------------------------------------------------------------------- +# End-to-end pipeline tests. These require bowtie2 / bwa / samtools / sambamba +# / bedtools / pigz to be available on PATH (CI provides them via +# environment.yml). +# ----------------------------------------------------------------------------- - runner = CliRunner() - results = runner.invoke(main, [ "--input-directory", sharedir, - "--working-directory", directory.name, "--force", "--rRNA-file", sharedir+"feature.fasta"]) - assert results.exit_code == 0 +def _assert_pipeline_outputs(wk): + """Shared post-run assertions for end-to-end tests.""" + assert os.path.exists(os.path.join(wk, "summary.html")), "summary.html missing" + assert os.path.exists(os.path.join(wk, "outputs/proportions.png")) + assert os.path.exists(os.path.join(wk, "outputs/RPKM.png")) + assert os.path.exists(os.path.join(wk, "multiqc/multiqc_report.html")) + # New explanatory sections from v1.2.0 must be embedded in the summary. + with open(os.path.join(wk, "summary.html")) as fin: + html = fin.read() + assert "How it works" in html + assert "Per-sequence distribution of ribosomal hits" in html + assert "Per-sequence RPKM distribution" in html def test_full_rRNA_file(): - with tempfile.TemporaryDirectory() as directory: - wk = directory + """Default aligner (bowtie2) + pre-built rRNA fasta.""" + with tempfile.TemporaryDirectory() as wk: cmd = f"sequana_ribofinder --input-directory {sharedir} " cmd += f"--working-directory {wk} --force --rRNA-file {sharedir}/feature.fasta" subprocess.call(cmd.split()) - stat = subprocess.call("sh ribofinder.sh".split(), cwd=wk) - assert os.path.exists(wk + "/summary.html") + rc = subprocess.call("bash ribofinder.sh".split(), cwd=wk) + assert rc == 0 + _assert_pipeline_outputs(wk) + # bowtie2 mapper directory should exist for each sample + assert os.path.isdir(os.path.join(wk, "data/bowtie2")) + def test_full_rRNA_extract(): - with tempfile.TemporaryDirectory() as directory: - wk = directory + """Default aligner (bowtie2) + on-the-fly rRNA extraction from GFF.""" + with tempfile.TemporaryDirectory() as wk: cmd = f"sequana_ribofinder --input-directory {sharedir} " - cmd += f"--working-directory {wk} --force --reference-file {sharedir}/Lepto.fa --gff-file {sharedir}/Lepto.gff" + cmd += ( + f"--working-directory {wk} --force " + f"--reference-file {sharedir}/Lepto.fa --gff-file {sharedir}/Lepto.gff" + ) subprocess.call(cmd.split()) - stat = subprocess.call("sh ribofinder.sh".split(), cwd=wk) - - - if os.path.exists(wk + "/summary.html"): - pass - else: - with open(f"{wk}/indexing/bowtie_rRNA.log", "r") as fout: - print(fout.read()) - raise IOError - -def test_version(): - cmd = "sequana_ribofinder --version" - subprocess.call(cmd.split()) - - + rc = subprocess.call("bash ribofinder.sh".split(), cwd=wk) + if rc != 0: + log = os.path.join(wk, "indexing/bowtie2_rRNA.log") + if os.path.exists(log): + with open(log) as fout: + print(fout.read()) + raise IOError("pipeline failed") + _assert_pipeline_outputs(wk) + + +def test_full_bwa(): + """--aligner bwa must produce the same set of outputs as the bowtie2 run.""" + with tempfile.TemporaryDirectory() as wk: + cmd = f"sequana_ribofinder --input-directory {sharedir} " + cmd += ( + f"--working-directory {wk} --force " + f"--rRNA-file {sharedir}/feature.fasta --aligner bwa" + ) + subprocess.call(cmd.split()) + rc = subprocess.call("bash ribofinder.sh".split(), cwd=wk) + assert rc == 0 + _assert_pipeline_outputs(wk) + # bwa-specific outputs: per-sample sorted BAM and samtools stats file. + assert os.path.isdir(os.path.join(wk, "data/bwa")) + assert os.path.exists( + os.path.join(wk, "data/samtools_stats/data.stats.txt") + ) + # multiqc should pick up samtools-stats data when aligner == bwa. + assert os.path.exists( + os.path.join(wk, "multiqc/multiqc_data/multiqc_samtools_stats.txt") + ) + + +def test_full_bwa_extract(): + """--aligner bwa combined with reference+gff extraction.""" + with tempfile.TemporaryDirectory() as wk: + cmd = f"sequana_ribofinder --input-directory {sharedir} " + cmd += ( + f"--working-directory {wk} --force " + f"--reference-file {sharedir}/Lepto.fa " + f"--gff-file {sharedir}/Lepto.gff --aligner bwa" + ) + subprocess.call(cmd.split()) + rc = subprocess.call("bash ribofinder.sh".split(), cwd=wk) + assert rc == 0 + _assert_pipeline_outputs(wk)