From e7045d53cbbc10eed2d39ab2bc3790ef9cad3907 Mon Sep 17 00:00:00 2001 From: Thomas Cokelaer Date: Fri, 21 Aug 2026 23:46:37 +0200 Subject: [PATCH 1/3] Fix the anchor keys of the summary sections The key was written "anchor'" with a stray quote in the three sections, so they had no anchor and their sidebar entries pointed to nothing. The adapter section was also a copy of the phix one and reused its anchor; it now has its own. --- sequana_pipelines/quality_control/quality_control.rules | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sequana_pipelines/quality_control/quality_control.rules b/sequana_pipelines/quality_control/quality_control.rules index 5b39df6..dbb757d 100644 --- a/sequana_pipelines/quality_control/quality_control.rules +++ b/sequana_pipelines/quality_control/quality_control.rules @@ -613,7 +613,7 @@ onsuccess: conf.summary_sections.append({ "name": "Stats (input data)", "title_links": 'FastQC', - "anchor'": "stats", + "anchor": "stats", "content": fqmod._get_stats_section() }) @@ -625,7 +625,7 @@ onsuccess: conf.summary_sections.append({ "name": "Phix ", "title_links": 'FastQC', - "anchor'": "phix", + "anchor": "phix", "content": phixmod._get_html() }) @@ -640,7 +640,7 @@ onsuccess: "name": "Adapter ", "title_links": 'FastQC|' + 'Cutadapt|', - "anchor'": "phix", + "anchor": "adapter", "content": html }) From 1fec3e7a126a9be71758c54fba5d58c0baaf31b2 Mon Sep 17 00:00:00 2001 From: Thomas Cokelaer Date: Sat, 22 Aug 2026 23:12:14 +0200 Subject: [PATCH 2/3] Fix the broken options and the report of the quality_control pipeline Two of the four pipeline options could not even build the workflow: - --skip-phix-removal used an undefined __data__input variable; the raw data is now trimmed directly - --disable-trimming referenced the cutadapt output, which does not exist in that case; the summary rule now follows what was really executed. The cutadapt section of the configuration file is kept in sync by the CLI. The report was showing and linking the wrong files: - the adapter section displayed the statistics computed after the phix removal and linked fastqc_cutadapt.html, a file that is never created. A fastq_stat_trimmed rule now provides the statistics of the trimmed data and the links point to the fastqc_trimmed report. - the statistics and plots were computed for R1 only. They are now computed for each read, which is what the report module expects. - a failure of the final summary.html was caught and printed, leaving no report behind while the pipeline reported a success. - the pipeline name is used in an URL, so the space was removed. Also removed dead code: the unreachable fastp rules, the kraken entries of the rulegraph mapper, unused imports, and the duplicated statistics block. The missing sys import made an invalid trimmer name raise a NameError. The default configuration file no longer ships a hardcoded input directory. Tests cover the four options with a dry run; they would have caught the two NameErrors. --- README.rst | 7 + pyproject.toml | 2 +- sequana_pipelines/quality_control/config.yaml | 2 +- sequana_pipelines/quality_control/main.py | 2 + .../quality_control/quality_control.rules | 260 ++++++++---------- test/test_main.py | 35 ++- 6 files changed, 154 insertions(+), 154 deletions(-) diff --git a/README.rst b/README.rst index aa0d133..6e13faf 100644 --- a/README.rst +++ b/README.rst @@ -93,6 +93,13 @@ ChangeLog ========= ==================================================================== Version Description ========= ==================================================================== +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) + * the adapter section of the report shows the trimmed data (it was + showing the phix data) and links the FastQC reports that exist + * do not hide a failure of the final summary.html anymore + * remove the dead fastp rules and the kraken entries of the rulegraph 1.0.0 * switch to click, pyproject. remove kraken (see multitax pipeline instead). Uses new convention. 0.10.0 * add missing MANIFEST diff --git a/pyproject.toml b/pyproject.toml index 3e2a486..e7d1b37 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "poetry.core.masonry.api" [tool.poetry] name = "sequana-quality-control" -version = "1.1.0" +version = "1.2.0" description = "Quality control pipeline for NGS data (phix removal, adapter trimming, FastQC)" authors = ["Sequana Team"] license = "BSD-3" diff --git a/sequana_pipelines/quality_control/config.yaml b/sequana_pipelines/quality_control/config.yaml index 50f1227..6897d79 100644 --- a/sequana_pipelines/quality_control/config.yaml +++ b/sequana_pipelines/quality_control/config.yaml @@ -6,7 +6,7 @@ # If input_directory provided, use it otherwise if input_pattern provided, # use it, otherwise use input_samples. # ============================================================================ -input_directory: /home/cokelaer/Data/Hm2 +input_directory: input_readtag: _R[12]_ input_pattern: '*fastq.gz' reference_file: phiX174.fa diff --git a/sequana_pipelines/quality_control/main.py b/sequana_pipelines/quality_control/main.py index 4d48fa0..91e3b76 100755 --- a/sequana_pipelines/quality_control/main.py +++ b/sequana_pipelines/quality_control/main.py @@ -64,6 +64,8 @@ def main(**options): # --------------------------------------------------------- trimming cfg.trimming.software_choice = options.trimming_software_choice cfg.trimming.do = not options.disable_trimming + # the cutadapt section is the one used by the trimming rules + cfg.cutadapt.do = cfg.trimming.do qual = options.trimming_quality if options.trimming_software_choice in ["cutadapt", "atropos"]: diff --git a/sequana_pipelines/quality_control/quality_control.rules b/sequana_pipelines/quality_control/quality_control.rules index dbb757d..87b0bf2 100644 --- a/sequana_pipelines/quality_control/quality_control.rules +++ b/sequana_pipelines/quality_control/quality_control.rules @@ -12,16 +12,9 @@ ############################################################################## import glob import os -import shutil -import subprocess - - -# Some sequana related tools -import sequana - +import sys from sequana_pipetools import PipelineManager -from sequana_pipetools import snaketools as sm # This must be defined before the include configfile: "config.yaml" @@ -30,22 +23,68 @@ configfile: "config.yaml" manager = PipelineManager("quality_control", config) +# the statistics and plots are computed for each read (R1 and, if paired, R2) +read_tags = ["_R1_", "_R2_"] if manager.paired else ["_R1_"] + + +def fastq_stats_output(directory): + """Expected json/png files of a fastq_stat rule (one set per read)""" + return { + "json": [f"{{sample}}/{directory}/{{sample}}{tag}.json" for tag in read_tags], + "gc": [f"{{sample}}/{directory}/{{sample}}{tag}_gc.png" for tag in read_tags], + "boxplot": [f"{{sample}}/{directory}/{{sample}}{tag}_boxplot.png" for tag in read_tags], + } + + +def compute_fastq_stats(fastqs, output, max_reads): + """Compute the sequana statistics and plots of each input FastQ file""" + import shutil + from pathlib import Path + + import pylab + from sequana import FastQC, sequana_data + + pylab.ioff() + for filename, json_file, gc_file, boxplot_file in zip(fastqs, output.json, output.gc, output.boxplot): + fastq = FastQC(filename, max_sample=max_reads) + if len(fastq.fastq) != 0: + pylab.clf() + fastq.boxplot_quality() + pylab.savefig(boxplot_file) + pylab.clf() + fastq.histogram_gc_content() + pylab.savefig(gc_file) + fastq.get_stats().to_json(json_file) + else: + location = sequana_data("no_data.jpg", "images") + shutil.copy(location, gc_file) + shutil.copy(location, boxplot_file) + Path(json_file).touch() + + +__fastq_stats_samples__output = fastq_stats_output("fastq_stats_samples") +__fastq_stats_phix__output = fastq_stats_output("fastq_stats_phix") +__fastq_stats_trimmed__output = fastq_stats_output("fastq_stats_trimmed") + + expected_output = [] -# stats on raw data and optional fastac -expected_output += expand("{sample}/fastq_stats_samples/{sample}.json", sample=manager.samples) +# stats on raw data and optional fastqc +expected_output += expand(__fastq_stats_samples__output["json"], sample=manager.samples) if manager.config.fastqc.do_raw: expected_output += expand("{sample}/fastqc_raw/fastqc.done", sample=manager.samples) # if we remove the phix if manager.config.bwa_mem_phix.do: - expected_output += expand("{sample}/fastq_stats_phix/{sample}.json", sample=manager.samples) + expected_output += expand(__fastq_stats_phix__output["json"], sample=manager.samples) if manager.config.fastqc.do_after_phix_removal: expected_output += expand("{sample}/fastqc_phix/fastqc.done", sample=manager.samples) - -if manager.config.fastqc.do_after_adapter_removal: - expected_output += expand("{sample}/fastqc_trimmed/fastqc.done", sample=manager.samples) +# the trimming outputs only exist if the trimming was performed +if manager.config.trimming.do: + expected_output += expand(__fastq_stats_trimmed__output["json"], sample=manager.samples) + if manager.config.fastqc.do_after_adapter_removal: + expected_output += expand("{sample}/fastqc_trimmed/fastqc.done", sample=manager.samples) rule pipeline: @@ -87,37 +126,16 @@ if manager.config.fastqc.do_raw: # FASTQ stats on input data set rule fastq_stat_samples: - input: manager.getrawdata() + input: + fastq=manager.getrawdata() output: - json="{sample}/fastq_stats_samples/{sample}.json", - gc="{sample}/fastq_stats_samples/{sample}_gc.png", - boxplot="{sample}/fastq_stats_samples/{sample}_boxplot.png" + json=__fastq_stats_samples__output["json"], + gc=__fastq_stats_samples__output["gc"], + boxplot=__fastq_stats_samples__output["boxplot"] params: max_reads=config['fastq_stats']['max_reads'] run: - import shutil - from pathlib import Path - import pylab - from sequana import FastQC, sequana_data - from sequana_pipetools.snaketools import FileFactory - - pylab.ioff() - ff = FileFactory(input[0]) - for filename in ff.realpaths: - fastq = FastQC(filename, max_sample=params.max_reads) - if len(fastq.fastq) != 0: - pylab.clf() - fastq.boxplot_quality() - pylab.savefig(output.boxplot) - pylab.clf() - fastq.histogram_gc_content() - pylab.savefig(output.gc) - fastq.get_stats().to_json(output.json) - else: - location = sequana_data("no_data.jpg", "images") - shutil.copy(location, output.gc) - shutil.copy(location, output.boxplot) - Path(output.json).touch() + compute_fastq_stats(input.fastq, output, params.max_reads) @@ -253,44 +271,22 @@ if manager.config.bwa_mem_phix.do: rule fastq_stat_phix: input: - rules.compress_phix.output.fastq + # the unmapped reads are the ones kept after the phix removal + fastq=[x for x in __bwa_bam_to_fastq__fastq_output_gz if "unmapped" in x] output: - json="{sample}/fastq_stats_phix/{sample}.json", - gc="{sample}/fastq_stats_phix/{sample}_gc.png", - boxplot="{sample}/fastq_stats_phix/{sample}_boxplot.png" + json=__fastq_stats_phix__output["json"], + gc=__fastq_stats_phix__output["gc"], + boxplot=__fastq_stats_phix__output["boxplot"] params: max_reads=config['fastq_stats']['max_reads'] run: - import shutil - from pathlib import Path - import pylab - from sequana import FastQC, sequana_data - from sequana_pipetools.snaketools import FileFactory - - pylab.ioff() - ff = FileFactory(input[0]) - for filename in ff.realpaths: - fastq = FastQC(filename, max_sample=params.max_reads) - if len(fastq.fastq) != 0: - pylab.clf() - fastq.boxplot_quality() - pylab.savefig(output.boxplot) - pylab.clf() - fastq.histogram_gc_content() - pylab.savefig(output.gc) - fastq.get_stats().to_json(output.json) - else: - location = sequana_data("no_data.jpg", "images") - shutil.copy(location, output.gc) - shutil.copy(location, output.boxplot) - Path(output.json).touch() + compute_fastq_stats(input.fastq, output, params.max_reads) valid_trimmer = ['cutadapt', 'atropos'] if manager.config.trimming.software_choice not in valid_trimmer: - print(f"Invalid choice for trimming tool. Choose one in {valid_trimmer}") - sys.exit(1) + sys.exit(f"Invalid choice for trimming tool. Choose one in {valid_trimmer}") # Perform the adapter removal and trimming if manager.config['trimming']['do']: @@ -308,9 +304,9 @@ if manager.config['trimming']['do']: "cutadapt").replace("unmapped","cutadapt") for x in __cutadapt__input_fastq] else: - # If the fix is not yet performed, __data__input is a wildcard - # function so the output must be specified by hand - __cutadapt__input_fastq = __data__input + # the phix removal is skipped so the raw data is trimmed directly. + # getrawdata() is a wildcard function, hence the explicit output + __cutadapt__input_fastq = manager.getrawdata() if manager.paired: __cutadapt__output = [ "{sample}/cutadapt/{sample}_R1_.cutadapt.fastq.gz", @@ -407,48 +403,21 @@ if manager.config['trimming']['do']: cmd += " > {log}" shell(cmd) - elif adapter_tool in ['fastp']: - if manager.paired: - rule fastp: - input: - fastq=manager.getrawdata() - output: - r1="{sample}/fastp/{sample}_R1_.fastp.fastq.gz", - r2="{sample}/fastp/{sample}_R2_.fastp.fastq.gz", - html="{sample}/fastp/fastp_{sample}.html", - json="{sample}/fastp/fastp_{sample}.json", - log: - "logs/fastp/{sample}.log" - params: - options=config['fastp']["options"], - adapters=config["fastp"]["adapters"] - threads: - config["fastp"].get("threads", 4) - container: - config['apptainers']['fastp'] - shell: - manager.get_shell("fastp/run", "v1") - else: - rule fastp: - input: - fastq=manager.getrawdata() - output: - r1="{sample}/fastp/{sample}_R1_.fastp.fastq.gz", - html="{sample}/fastp/fastp_{sample}.html", - json="{sample}/fastp/fastp_{sample}.json", - log: - "logs/fastp/{sample}.log" - params: - options=config['fastp']["options"], - adapters=config["fastp"]["adapters"] - threads: - config["fastp"].get("threads", 4) - container: - config['apptainers']['fastp'] - shell: - manager.get_shell("fastp/run", "v1") else: - raise ValueError("trimming must be either cutadapt or atropos or fastp") + raise ValueError(f"trimming must be one of {valid_trimmer}") + + rule fastq_stat_trimmed: + input: + fastq=__cutadapt__output + output: + json=__fastq_stats_trimmed__output["json"], + gc=__fastq_stats_trimmed__output["gc"], + boxplot=__fastq_stats_trimmed__output["boxplot"] + params: + max_reads=config['fastq_stats']['max_reads'] + run: + compute_fastq_stats(input.fastq, output, params.max_reads) + # Now we can perform again a FastQC and FastQ stats if manager.config.fastqc.do_after_adapter_removal: @@ -472,8 +441,7 @@ if manager.config['trimming']['do']: # create a json file that summarise information of your pipeline __summary_pipeline__inputs = manager.getrawdata() -if manager.config['cutadapt'].do: - # todo: handle all adapter removal cases +if manager.config.trimming.do: __summary_pipeline__outputs = [ __cutadapt__output ] elif manager.config.bwa_mem_phix.do: __summary_pipeline__outputs = [x for x in __bwa_bam_to_fastq__fastq_output_gz if "unmapped" in x] @@ -490,7 +458,7 @@ __summary_pipeline__rulegraph = ".sequana/rulegraph.svg" __summary_pipeline__requirements = ".sequana/env.yaml" __summary_pipeline__snakefile = str(manager.snakefile) __summary_pipeline__config = "config.yaml" -__summary_pipeline__name = "Quality Control" +__summary_pipeline__name = "quality_control" __summary_pipeline__json_output = "{sample}/summary_pipeline/{sample}.json" rule summary: @@ -530,14 +498,11 @@ rule summary: # ========================================================== rulegraph sequana_rulegraph_mapper = { - "fastqc_raw": "../fastqc_raw.html", - "fastqc_phix": "../fastqc_phix.html", - "fastqc_trimmed": "../fastqc_trimmed.html", - "cutadapt": "../cutadapt.html", - "kraken": "../kraken/kraken/kraken.html", - "kraken_translate": "../kraken/raken/kraken.html", - "kraken_to_krona": "../kraken/kraken/kraken.html", - } + "fastqc_samples": "../fastqc_raw.html", + "fastqc_phix": "../fastqc_phix.html", + "fastqc_trimmed": "../fastqc_trimmed.html", + "cutadapt": "../cutadapt.html", +} rule rulegraph: input: str(manager.snakefile) output: @@ -596,10 +561,11 @@ onsuccess: report_dir = report_dir_format % {"proj": proj} conf.output_dir = report_dir # ensure files are stored in the correct location - # Create the 3 FastQC HTML files (independent) + # Create the FastQC HTML files (independent) for this in ["fastqc_phix", "fastqc_raw", "fastqc_trimmed"]: - FastQCModule("{}.html".format(this), - proj + "/{}/*_fastqc.html".format(this)) + if glob.glob(proj + "/{}/*_fastqc.html".format(this)): + FastQCModule("{}.html".format(this), + proj + "/{}/*_fastqc.html".format(this)) # add all sections in addition to standard summary conf.summary_sections = [] @@ -622,30 +588,33 @@ onsuccess: phixmod = PhixModule(proj) sample_summary["phix_section_json"] = json.loads(phixmod._get_stats().to_json()) sample_summary["phix_section"] = phixmod._get_summary() + html = phixmod._get_html() + html += FastQStatsModule(proj + "/fastq_stats_phix", "fastqc_phix")._get_stats_section( + tablename="phix2") conf.summary_sections.append({ - "name": "Phix ", + "name": "Phix", "title_links": 'FastQC', "anchor": "phix", - "content": phixmod._get_html() + "content": html }) # the cutadapt section ---------------------------- - if manager.config.cutadapt.do: + if manager.config.trimming.do: cutadapt_mod = CutadaptModule(f"{proj}/logs/cutadapt/cutadapt.txt", proj) sample_summary["cutadapt_json"] = json.loads(cutadapt_mod._get_stats().to_json()) html = cutadapt_mod._get_stat_section() - html += FastQStatsModule(proj + "/fastq_stats_phix", "fastqc_cutadapt")._get_stats_section( - tablename="cutadapt2") + html += FastQStatsModule(proj + "/fastq_stats_trimmed", "fastqc_trimmed")._get_stats_section( + tablename="trimmed") conf.summary_sections.append({ - "name": "Adapter ", - "title_links": 'FastQC|' + - 'Cutadapt|', + "name": "Adapter", + "title_links": 'FastQC | ' + 'Cutadapt', "anchor": "adapter", "content": html }) # The cutadapt report (independent) - if manager.config['cutadapt'].do: + if manager.config.trimming.do: from sequana.modules_report.cutadapt import CutadaptModule filename = proj + "/logs/cutadapt/cutadapt.txt" CutadaptModule(filename, proj, "cutadapt.html") @@ -666,7 +635,8 @@ onsuccess: "quality_control", __summary_pipeline__name) - SequanaReport(json.loads(open(filename).read()), intro=intro) + with open(filename, "r") as fin: + SequanaReport(json.loads(fin.read()), intro=intro) # save sample summary sample_summary['project'] = proj @@ -674,14 +644,12 @@ onsuccess: data = json.dumps(sample_summary, indent=4, sort_keys=True) fp.write(data) - try: - logger.info("Creating multi summary file") - from sequana.modules_report.multi_summary import MultiSummary - from sequana.utils import config as cfg - cfg.output_dir = "." - sms = MultiSummary(pattern="**/report_qc*/summary.json", output_filename="summary.html") - except Exception as err: - print(err) + logger.info("Creating multi summary file") + from sequana.modules_report.multi_summary import MultiSummary + from sequana.utils import config as cfg + + cfg.output_dir = "." + MultiSummary(pattern="**/report_qc*/summary.json", output_filename="summary.html") manager.teardown() onerror: diff --git a/test/test_main.py b/test/test_main.py index 2f90382..abdaae1 100644 --- a/test/test_main.py +++ b/test/test_main.py @@ -16,12 +16,9 @@ def test_standalone_subprocess(): - directory = tempfile.TemporaryDirectory() - cmd = """sequana_quality_control --input-directory {} - --working-directory --force""".format( - sharedir, directory.name - ) - subprocess.call(cmd.split()) + with tempfile.TemporaryDirectory() as directory: + cmd = f"sequana_quality_control --input-directory {sharedir} --working-directory {directory} --force" + assert subprocess.call(cmd.split()) == 0 def test_standalone_script(): @@ -51,3 +48,29 @@ def test_full(): def test_version(): cmd = "sequana_quality_control --version" subprocess.call(cmd.split()) + + +def dryrun(*args): + """Build the workflow with the given options and check that its DAG is valid""" + with tempfile.TemporaryDirectory() as directory: + cmd = ["sequana_quality_control", "--input-directory", sharedir, "--working-directory", directory, "--force"] + assert subprocess.call(cmd + list(args)) == 0 + + cmd = ["snakemake", "-s", "quality_control.rules", "--configfile", "config.yaml", "-n"] + assert subprocess.call(cmd, cwd=directory) == 0 + + +def test_skip_phix_removal(): + dryrun("--skip-phix-removal") + + +def test_disable_trimming(): + dryrun("--disable-trimming") + + +def test_skip_phix_removal_and_trimming(): + dryrun("--skip-phix-removal", "--disable-trimming") + + +def test_skip_fastqc(): + dryrun("--skip-fastqc-raw", "--skip-fastqc-cleaned") From 7e1353b290c0543d79c81c81521b30358fe0f509 Mon Sep 17 00:00:00 2001 From: Thomas Cokelaer Date: Sat, 22 Aug 2026 23:35:42 +0200 Subject: [PATCH 3/3] Do not nest the trimming outputs in the summary rule input --- sequana_pipelines/quality_control/quality_control.rules | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sequana_pipelines/quality_control/quality_control.rules b/sequana_pipelines/quality_control/quality_control.rules index 87b0bf2..ef1e804 100644 --- a/sequana_pipelines/quality_control/quality_control.rules +++ b/sequana_pipelines/quality_control/quality_control.rules @@ -442,7 +442,7 @@ if manager.config['trimming']['do']: # create a json file that summarise information of your pipeline __summary_pipeline__inputs = manager.getrawdata() if manager.config.trimming.do: - __summary_pipeline__outputs = [ __cutadapt__output ] + __summary_pipeline__outputs = __cutadapt__output elif manager.config.bwa_mem_phix.do: __summary_pipeline__outputs = [x for x in __bwa_bam_to_fastq__fastq_output_gz if "unmapped" in x] else: