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..ce086e9 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
@@ -45,7 +49,14 @@ def compute_fastq_stats(fastqs, output, max_reads):
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")
fastq = FastQC(filename, max_sample=max_reads)
if len(fastq.fastq) != 0:
pylab.clf()
@@ -86,6 +97,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 +153,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 +246,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 +305,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 +324,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 +343,107 @@ 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:
+ 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.stdout} 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:
+ 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.stdout} 2>&1
+ """
else:
raise ValueError(f"trimming must be one of {valid_trimmer}")
@@ -415,8 +457,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 +483,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 +592,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 +641,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
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")