From 0654fa91c1e464e97cca74b7ed49aa5b6f795b14 Mon Sep 17 00:00:00 2001 From: Thomas Cokelaer Date: Sat, 22 Aug 2026 23:57:12 +0200 Subject: [PATCH 1/2] Add a multiqc report and containerize the trimming MultiQC aggregates the FastQC and cutadapt results in multiqc/multiqc_report.html and each sample report links it. The rule uses the multiqc container and the sequana_wrappers shell, as in the other pipelines. The trimming used to build its command line inside a python 'run' directive, which cannot be executed in a container. It is now a plain shell rule using the cutadapt container (cutadapt 4.6). The atropos variant is kept as a separate rule; atropos has no container so it must be available locally. Its own report is no longer overwritten by the standard output. The rules that remain in python (the fastq statistics and the extraction of the mapped reads) now have a log file, which they write as they go: the pipeline monitor follows the modification times to report the elapsed time. Those rules use the sequana library directly, which the report generation requires on the host anyway, hence no container for them. Also require sequana >=0.24.1, declare multiqc in tools.txt and in the conda environment, and run the summary rule locally. --- README.rst | 6 + environment.yml | 1 + pyproject.toml | 6 +- sequana_pipelines/quality_control/config.yaml | 20 ++ .../quality_control/quality_control.rules | 252 +++++++++++------- sequana_pipelines/quality_control/schema.yaml | 19 ++ sequana_pipelines/quality_control/tools.txt | 1 + 7 files changed, 207 insertions(+), 98 deletions(-) diff --git a/README.rst b/README.rst index 6e13faf..9be3d2e 100644 --- a/README.rst +++ b/README.rst @@ -93,6 +93,12 @@ ChangeLog ========= ==================================================================== Version Description ========= ==================================================================== +1.3.0 * add a multiqc report aggregating the FastQC and cutadapt results + * the trimming is now a shell rule running in the cutadapt + container (it was a python rule using the local installation) + * add the missing log files so that the pipeline monitor reports + the correct elapsed times + * require sequana >=0.24.1 1.2.0 * fix --skip-phix-removal and --disable-trimming, which both ended in a NameError * compute the FastQ statistics of R2 as well (R1 only before) diff --git a/environment.yml b/environment.yml index 260b8ce..695091a 100644 --- a/environment.yml +++ b/environment.yml @@ -10,6 +10,7 @@ dependencies: - fastqc - sambamba - graphviz +- multiqc - pulp>=2.8 - pip: - sequana diff --git a/pyproject.toml b/pyproject.toml index e7d1b37..16d51a4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,13 +4,13 @@ build-backend = "poetry.core.masonry.api" [tool.poetry] name = "sequana-quality-control" -version = "1.2.0" +version = "1.3.0" description = "Quality control pipeline for NGS data (phix removal, adapter trimming, FastQC)" authors = ["Sequana Team"] license = "BSD-3" repository = "https://github.com/sequana/quality_control" readme = "README.rst" -keywords = ["snakemake", "NGS", "sequana", "fastqc", "cutadapt"] +keywords = ["snakemake", "NGS", "sequana", "fastqc", "cutadapt", "multiqc"] classifiers = [ "Development Status :: 5 - Production/Stable", "Intended Audience :: Education", @@ -35,7 +35,7 @@ packages = [ [tool.poetry.dependencies] python = ">=3.9,<4.0" -sequana = ">=0.18.0" +sequana = ">=0.24.1" sequana_pipetools = ">=1.1.1" sequana-wrappers = ">=26.4.16" snakemake = ">=7.32" diff --git a/sequana_pipelines/quality_control/config.yaml b/sequana_pipelines/quality_control/config.yaml index 6897d79..eabb211 100644 --- a/sequana_pipelines/quality_control/config.yaml +++ b/sequana_pipelines/quality_control/config.yaml @@ -13,6 +13,8 @@ reference_file: phiX174.fa apptainers: + cutadapt: https://zenodo.org/record/18834062/files/cutadapt_4.6.0.img + multiqc: https://zenodo.org/record/10155626/files/multiqc_1.17.0.img pigz: https://zenodo.org/record/7346805/files/pigz_2.4.0.img fastqc: https://zenodo.org/record/7015004/files/fastqc_0.11.9-py3.img sequana_minimal: https://zenodo.org/record/19127803/files/sequana_minimal_26.3.19.img @@ -109,3 +111,21 @@ fastqc: fastq_stats: max_reads: 500000 + + +############################################################################# +# MultiQC aggregates the FastQC and cutadapt reports in a single HTML file +# +# :Parameters: +# +# - do: if unchecked, this rule is ignored +# - options: any options recognised by multiqc +# - input_directory: where to look for the reports +# - modules: the multiqc modules to use (space separated) +# - config_file: an optional multiqc configuration file +multiqc: + do: true + options: -p -f + input_directory: '.' + modules: fastqc cutadapt + config_file: diff --git a/sequana_pipelines/quality_control/quality_control.rules b/sequana_pipelines/quality_control/quality_control.rules index ef1e804..a44a6fe 100644 --- a/sequana_pipelines/quality_control/quality_control.rules +++ b/sequana_pipelines/quality_control/quality_control.rules @@ -36,8 +36,12 @@ def fastq_stats_output(directory): } -def compute_fastq_stats(fastqs, output, max_reads): - """Compute the sequana statistics and plots of each input FastQ file""" +def compute_fastq_stats(fastqs, output, max_reads, logfile): + """Compute the sequana statistics and plots of each input FastQ file + + The progress is written in *logfile* as we go: the pipeline monitor follows + the modification time of the log files to report the elapsed time. + """ import shutil from pathlib import Path @@ -46,6 +50,8 @@ def compute_fastq_stats(fastqs, output, max_reads): pylab.ioff() for filename, json_file, gc_file, boxplot_file in zip(fastqs, output.json, output.gc, output.boxplot): + with open(logfile, "a") as fout: + fout.write(f"Computing the statistics of {filename}\n") fastq = FastQC(filename, max_sample=max_reads) if len(fastq.fastq) != 0: pylab.clf() @@ -86,6 +92,14 @@ if manager.config.trimming.do: if manager.config.fastqc.do_after_adapter_removal: expected_output += expand("{sample}/fastqc_trimmed/fastqc.done", sample=manager.samples) +# multiqc aggregates all the reports, hence it must be the last rule to run +__multiqc__input = list(expected_output) +if manager.config.trimming.do: + __multiqc__input += expand("{sample}/logs/cutadapt/cutadapt.txt", sample=manager.samples) + +if config['multiqc']['do']: + expected_output += ["multiqc/multiqc_report.html"] + rule pipeline: input: ".sequana/rulegraph.svg", @@ -134,8 +148,10 @@ rule fastq_stat_samples: boxplot=__fastq_stats_samples__output["boxplot"] params: max_reads=config['fastq_stats']['max_reads'] + log: + "logs/{sample}/fastq_stats_samples.log" run: - compute_fastq_stats(input.fastq, output, params.max_reads) + compute_fastq_stats(input.fastq, output, params.max_reads, log[0]) @@ -225,9 +241,14 @@ if manager.config.bwa_mem_phix.do: stats=__bwa_bam_to_fastq__stats_output params: wkdir = __bwa_bam_to_fastq__wkdir + log: + "logs/{sample}/bwa_bam_to_fastq.log" run: # save some stats for later from sequana.tools import StatsBAM2Mapped + + with open(log[0], "w") as fout: + fout.write(f"Extracting the mapped and unmapped reads of {input.bam}\n") stats = StatsBAM2Mapped(input["bam"], wkdir=params.wkdir) stats.to_json(output["stats"]) @@ -279,8 +300,10 @@ if manager.config.bwa_mem_phix.do: boxplot=__fastq_stats_phix__output["boxplot"] params: max_reads=config['fastq_stats']['max_reads'] + log: + "logs/{sample}/fastq_stats_phix.log" run: - compute_fastq_stats(input.fastq, output, params.max_reads) + compute_fastq_stats(input.fastq, output, params.max_reads, log[0]) @@ -296,7 +319,7 @@ if manager.config['trimming']['do']: rev = manager.config.cutadapt.rev if adapter_tool in ["cutadapt", "atropos"]: - adapter_tool = "cutadapt" + executable = config['cutadapt']['tool_choice'] if manager.config.bwa_mem_phix.do: __cutadapt__input_fastq = [x for x in __bwa_bam_to_fastq__fastq_output_gz if "unmapped" in x] @@ -315,93 +338,105 @@ if manager.config['trimming']['do']: __cutadapt__output = [ "{sample}/cutadapt/{sample}_R1_.cutadapt.fastq.gz"] - rule cutadapt: - input: - fastq = __cutadapt__input_fastq - output: - fastq = __cutadapt__output - params: - fwd = manager.config.cutadapt.fwd, - rev = manager.config.cutadapt.rev, - m = config['cutadapt']['m'], # cutoff - quality = config['cutadapt']['quality'], - kwargs = manager.config.cutadapt.options, - mode = manager.config.cutadapt.mode, - wkdir = "{sample}/cutadapt", - sample = "{sample}", - threads: - config['cutadapt']['threads'] - log: "{sample}/logs/cutadapt/cutadapt.txt" - run: - executable = config['cutadapt']['tool_choice'] # could be cutadapt or atropos - - # Fill the fwd and revcomp variables depending on the config file - fwd = params.fwd - revcomp = params.rev - - # For atropos, output must be txt AND json for the module to work. - # - # Paired data - if len(input.fastq) == 2: - if executable == "atropos": - cmd = "%s trim" % executable - cmd += " -pe1 {input.fastq[0]} -pe2 {input.fastq[1]} " - # atropos requires at least 2 threads - if int(threads) > 1: - cmd += " -T {threads} " - cmd += " --process-timeout 600 " - logpath = log[0].replace(".txt", "") - cmd += " --report-formats json txt --report-file %s" % logpath - else: - cmd = "{}".format(executable) - mode2 = params.mode.upper() - - if fwd: - cmd += " -{params.mode} %s " % fwd - else: - cmd += " -{params.mode} XXXX " - - if revcomp: - cmd += " -%s %s " % (mode2, revcomp) - else: - cmd += " -%s XXXX " % mode2 - - cmd += " -m {params.m} -q {params.quality} " - # -o is common to atropos and cutadapt for the output R1 and -p for - # R2 - cmd += " {params.kwargs} -o {output.fastq[0]} -p {output.fastq[1]}" - - # input of cutdapt must be at the end. earlier version handled - # mixing of optional arguments and positional argument but not - # latest version. In any case, it is safer to do it correctly - # (changes sequana 0.8) - if executable == "cutadapt": - cmd += " {input.fastq[0]} {input.fastq[1]} " - cmd += "> {log}" - shell(cmd) - else: - if executable == "atropos": - cmd = "%s trim -se {input.fastq[0]} " % executable - # atropos requires at least 2 threads - - if int(threads) > 1: - cmd += " -T {threads} " - cmd += " --process-timeout 600 " - logpath = log[0].replace(".txt", "") - cmd += " --report-formats json txt --report-file %s" % logpath - else: - cmd = "%s" % executable - - if fwd: - cmd += " -{params.mode} %s " % fwd - else: - cmd += " -{params.mode} XXXX " - cmd += " -m {params.m} -q {params.quality} " - cmd += " {params.kwargs} -o {output.fastq[0]}" - if executable == "cutadapt": - cmd += " {input.fastq[0]} " - cmd += " > {log}" - shell(cmd) + # the adapters are optional. A dummy one is used when none is provided so + # that the quality and length trimming is still performed + mode = manager.config.cutadapt.mode + adapter_options = f" -{mode} {fwd if fwd else 'XXXX'} " + if manager.paired: + adapter_options += f" -{mode.upper()} {rev if rev else 'XXXX'} " + + trimming_options = "{} -m {} -q {} {}".format( + adapter_options, + config['cutadapt']['m'], + config['cutadapt']['quality'], + manager.config.cutadapt.options) + + if executable == "cutadapt": + if manager.paired: + rule cutadapt: + input: + fastq=__cutadapt__input_fastq + output: + fastq=__cutadapt__output + params: + options=trimming_options + threads: + config['cutadapt']['threads'] + log: + "{sample}/logs/cutadapt/cutadapt.txt" + container: + config['apptainers']['cutadapt'] + shell: + """ + cutadapt -j {threads} {params.options} \ + -o {output.fastq[0]} -p {output.fastq[1]} \ + {input.fastq[0]} {input.fastq[1]} > {log} 2>&1 + """ + else: + rule cutadapt: + input: + fastq=__cutadapt__input_fastq + output: + fastq=__cutadapt__output + params: + options=trimming_options + threads: + config['cutadapt']['threads'] + log: + "{sample}/logs/cutadapt/cutadapt.txt" + container: + config['apptainers']['cutadapt'] + shell: + """ + cutadapt -j {threads} {params.options} \ + -o {output.fastq[0]} {input.fastq[0]} > {log} 2>&1 + """ + else: + # atropos has no container; it must be available in the environment. + # Its own report (cutadapt.txt) is the one read by the HTML report, + # hence the separate log file for the standard output + if manager.paired: + rule atropos: + input: + fastq=__cutadapt__input_fastq + output: + fastq=__cutadapt__output + params: + options=trimming_options, + report="{sample}/logs/cutadapt/cutadapt" + threads: + config['cutadapt']['threads'] + log: + "{sample}/logs/cutadapt/atropos.log" + shell: + """ + atropos trim -pe1 {input.fastq[0]} -pe2 {input.fastq[1]} \ + -T {threads} --process-timeout 600 \ + --report-formats json txt --report-file {params.report} \ + {params.options} \ + -o {output.fastq[0]} -p {output.fastq[1]} > {log} 2>&1 + """ + else: + rule atropos: + input: + fastq=__cutadapt__input_fastq + output: + fastq=__cutadapt__output + params: + options=trimming_options, + report="{sample}/logs/cutadapt/cutadapt" + threads: + config['cutadapt']['threads'] + log: + "{sample}/logs/cutadapt/atropos.log" + shell: + """ + atropos trim -se {input.fastq[0]} \ + -T {threads} --process-timeout 600 \ + --report-formats json txt --report-file {params.report} \ + {params.options} \ + -o {output.fastq[0]} > {log} 2>&1 + """ else: raise ValueError(f"trimming must be one of {valid_trimmer}") @@ -415,8 +450,10 @@ if manager.config['trimming']['do']: boxplot=__fastq_stats_trimmed__output["boxplot"] params: max_reads=config['fastq_stats']['max_reads'] + log: + "logs/{sample}/fastq_stats_trimmed.log" run: - compute_fastq_stats(input.fastq, output, params.max_reads) + compute_fastq_stats(input.fastq, output, params.max_reads, log[0]) # Now we can perform again a FastQC and FastQ stats @@ -439,6 +476,27 @@ if manager.config['trimming']['do']: +# ============================================================ multiqc +if config['multiqc']['do']: + + rule multiqc: + input: + __multiqc__input + output: + "multiqc/multiqc_report.html" + params: + options=config['multiqc']['options'], + input_directory=config['multiqc']['input_directory'], + config_file=config['multiqc']['config_file'], + modules=config['multiqc']['modules'] + log: + "multiqc/multiqc.log" + container: + config['apptainers']['multiqc'] + shell: + manager.get_shell("multiqc/run", "v1") + + # create a json file that summarise information of your pipeline __summary_pipeline__inputs = manager.getrawdata() if manager.config.trimming.do: @@ -527,7 +585,7 @@ rule dot2svg: # Those rules takes a couple of seconds so no need for a cluster -localrules: rulegraph, dot2svg +localrules: rulegraph, dot2svg, summary onsuccess: @@ -576,9 +634,13 @@ onsuccess: fqmod = FastQStatsModule(indir, "fastqc_raw") sample_summary["fastq_stats_samples_json"] = json.loads(fqmod.get_stats().to_json()) + multiqc_link = '' + if config['multiqc']['do']: + multiqc_link = ' | MultiQC' + conf.summary_sections.append({ "name": "Stats (input data)", - "title_links": 'FastQC', + "title_links": 'FastQC' + multiqc_link, "anchor": "stats", "content": fqmod._get_stats_section() }) diff --git a/sequana_pipelines/quality_control/schema.yaml b/sequana_pipelines/quality_control/schema.yaml index 9a19a28..eac2243 100644 --- a/sequana_pipelines/quality_control/schema.yaml +++ b/sequana_pipelines/quality_control/schema.yaml @@ -54,6 +54,25 @@ mapping: type: bool + "multiqc": + type: map + mapping: + "do": + type: bool + default: True + "options": + type: str + required: False + "input_directory": + type: str + required: False + "modules": + type: str + required: False + "config_file": + type: str + required: False + "fastq_stats": type: map mapping: diff --git a/sequana_pipelines/quality_control/tools.txt b/sequana_pipelines/quality_control/tools.txt index b98550c..f570fb5 100644 --- a/sequana_pipelines/quality_control/tools.txt +++ b/sequana_pipelines/quality_control/tools.txt @@ -6,3 +6,4 @@ - cutadapt - atropos - graphviz +- multiqc From 8177299d750f6b7178bd764396f1b9e778e32e8d Mon Sep 17 00:00:00 2001 From: Thomas Cokelaer Date: Mon, 24 Aug 2026 11:32:22 +0200 Subject: [PATCH 2/2] Declare the atropos report and reset the statistics logs With atropos, the cutadapt.txt report was created as a side effect of --report-file. Nothing declared it, so multiqc, which reads that file, could not be resolved and the DAG could not be built: MissingInputException in rule multiqc: affected files: data/logs/cutadapt/cutadapt.txt The report is now declared in the log section of the atropos rules, next to the standard output. A test builds the workflow with atropos. The statistics rules also truncate their log file before writing, so a rerun no longer appends to the messages of the previous run. --- .../quality_control/quality_control.rules | 15 +++++++++++---- test/test_main.py | 5 +++++ 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/sequana_pipelines/quality_control/quality_control.rules b/sequana_pipelines/quality_control/quality_control.rules index a44a6fe..ce086e9 100644 --- a/sequana_pipelines/quality_control/quality_control.rules +++ b/sequana_pipelines/quality_control/quality_control.rules @@ -49,6 +49,11 @@ def compute_fastq_stats(fastqs, output, max_reads, logfile): from sequana import FastQC, sequana_data pylab.ioff() + + # start from a fresh log, otherwise a rerun appends to the previous messages + with open(logfile, "w"): + pass + for filename, json_file, gc_file, boxplot_file in zip(fastqs, output.json, output.gc, output.boxplot): with open(logfile, "a") as fout: fout.write(f"Computing the statistics of {filename}\n") @@ -407,14 +412,15 @@ if manager.config['trimming']['do']: threads: config['cutadapt']['threads'] log: - "{sample}/logs/cutadapt/atropos.log" + stdout="{sample}/logs/cutadapt/atropos.log", + report="{sample}/logs/cutadapt/cutadapt.txt" shell: """ atropos trim -pe1 {input.fastq[0]} -pe2 {input.fastq[1]} \ -T {threads} --process-timeout 600 \ --report-formats json txt --report-file {params.report} \ {params.options} \ - -o {output.fastq[0]} -p {output.fastq[1]} > {log} 2>&1 + -o {output.fastq[0]} -p {output.fastq[1]} > {log.stdout} 2>&1 """ else: rule atropos: @@ -428,14 +434,15 @@ if manager.config['trimming']['do']: threads: config['cutadapt']['threads'] log: - "{sample}/logs/cutadapt/atropos.log" + stdout="{sample}/logs/cutadapt/atropos.log", + report="{sample}/logs/cutadapt/cutadapt.txt" shell: """ atropos trim -se {input.fastq[0]} \ -T {threads} --process-timeout 600 \ --report-formats json txt --report-file {params.report} \ {params.options} \ - -o {output.fastq[0]} > {log} 2>&1 + -o {output.fastq[0]} > {log.stdout} 2>&1 """ else: diff --git a/test/test_main.py b/test/test_main.py index abdaae1..19ce15e 100644 --- a/test/test_main.py +++ b/test/test_main.py @@ -74,3 +74,8 @@ def test_skip_phix_removal_and_trimming(): def test_skip_fastqc(): dryrun("--skip-fastqc-raw", "--skip-fastqc-cleaned") + + +def test_atropos(): + # atropos writes the report read by multiqc and by the HTML report + dryrun("--software-choice", "atropos")