diff --git a/bin/read_lengths.py b/bin/read_lengths.py index d83235f..eac7a46 100755 --- a/bin/read_lengths.py +++ b/bin/read_lengths.py @@ -1,18 +1,38 @@ #!/usr/bin/env python - import argparse +import glob +import itertools +import os import pysam +import numpy as np + + +def mean_qual(quals): + qual = np.fromiter( + (ord(x) - 33 for x in quals), + dtype=int, count=len(quals)) + mean_p = np.mean(np.power(10, qual / -10)) + return -10 * np.log10(mean_p) def main(): parser = argparse.ArgumentParser() - parser.add_argument('fasta') - parser.add_argument('output') + parser.add_argument("directory", help="Directory containing .fastq(.gz) files") + parser.add_argument("output", help="Output file") args = parser.parse_args() - with open(args.output, 'w') as fh: - for rec in pysam.FastxFile(args.fasta): - fh.write("{}\t{}\n".format(rec.name, len(rec.sequence))) + fastqs = glob.glob(os.path.join(args.directory, "*.fastq*")) + reads = itertools.chain.from_iterable( + pysam.FastxFile(fname) for fname in fastqs) -if __name__ == '__main__': + with open(args.output, "w") as fh: + # names as in Guppy + fh.write("read_id\tsequence_length_template\tmean_qscore_template\n") + for read in reads: + fh.write("\t".join(str(x) for x in ( + read.name, len(read.sequence), mean_qual(read.quality)))) + fh.write("\n") + + +if __name__ == "__main__": main() diff --git a/bin/report.py b/bin/report.py new file mode 100755 index 0000000..bc355b4 --- /dev/null +++ b/bin/report.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python + +import argparse +import glob +import numpy as np +import pandas as pd + +from bokeh.layouts import gridplot, layout +import aplanat +from aplanat import annot, hist, report + + +def read_files(summaries): + dfs = list() + for fname in sorted(summaries): + dfs.append(pd.read_csv(fname, sep="\t")) + return pd.concat(dfs) + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("report", help="Report output file") + parser.add_argument("summaries", nargs='+', help="Read summary file.") + args = parser.parse_args() + + report_doc = report.HTMLReport( + "Workflow Template Sequencing report", + "Results generated through the wf-template nextflow workflow by Oxford Nanopore Technologies") + + report_doc.markdown(''' +### Read Quality control +This section displays basic QC metrics indicating read data quality. +''') + + np_blue = '#0084A9' + np_dark_grey = '#455560' + np_light_blue = '#90C6E7' + + # read length summary + seq_summary = read_files(args.summaries) + total_bases = seq_summary['sequence_length_template'].sum() + mean_length = total_bases / len(seq_summary) + median_length = np.median(seq_summary['sequence_length_template']) + datas = [seq_summary['sequence_length_template']] + length_hist = hist.histogram( + datas, colors=[np_blue], bins=100, + title="Read length distribution.", + x_axis_label='Read Length / bases', + y_axis_label='Number of reads', + xlim=(0, 2000)) + length_hist = annot.subtitle( + length_hist, + "Mean: {:.0f}. Median: {:.0f}".format( + mean_length, median_length)) + + datas = [seq_summary['mean_qscore_template']] + mean_q, median_q = np.mean(datas[0]), np.median(datas[0]) + q_hist = hist.histogram( + datas, colors=[np_blue], bins=100, + title="Read quality score", + x_axis_label="Quality score", + y_axis_label="Number of reads", + xlim=(4, 25)) + q_hist = annot.subtitle( + q_hist, + "Mean: {:.0f}. Median: {:.0f}".format( + mean_q, median_q)) + + report_doc.plot(gridplot([[length_hist, q_hist]])) + + # Footer section + report_doc.markdown(''' +### About + +**Oxford Nanopore Technologies products are not intended for use for health assessment +or to diagnose, treat, mitigate, cure or prevent any disease or condition.** + +This report was produced using the [epi2me-labs/wf-template](https://github.com/epi2me-labs/wf-template). +The workflow can be run using `nextflow epi2me-labs/wf-template --help` + +--- +''') + + # write report + report_doc.write(args.report) + +if __name__ == "__main__": + main() diff --git a/environment.yaml b/environment.yaml index 25db667..4860b82 100644 --- a/environment.yaml +++ b/environment.yaml @@ -4,5 +4,6 @@ channels: - conda-forge - defaults dependencies: - - python==3.6 + - python==3.6.* + - aplanat - pysam diff --git a/main.nf b/main.nf index d31806e..2abb3ef 100644 --- a/main.nf +++ b/main.nf @@ -18,10 +18,10 @@ def helpMessage(){ Workflow template' Usage: - nextflow run epi2melabs/workflow-template [options] + nextflow run epi2melabs/wf-template [options] Script Options: - --fastq FILE Path to FASTQ file (required) + --fastq DIR Path to directory containing FASTQ files (required) --out_dir DIR Path for output (default: $params.out_dir) """ } @@ -37,11 +37,21 @@ process readSeqs { """ read_lengths.py $reads seqs.txt - sleep 60 """ } +process makeReport { + label "pysam" + input: + file "seqs.txt" + output: + file "report.html" + """ + report.py report.html seqs.txt + """ +} + // See https://github.com/nextflow-io/nextflow/issues/1636 // This is the only way to publish files from a workflow whilst @@ -65,9 +75,10 @@ workflow pipeline { take: reads main: - seqs = readSeqs(reads) + summary = readSeqs(reads) + report = makeReport(summary) emit: - seqs + summary.concat(report) } // entrypoint workflow @@ -86,7 +97,12 @@ workflow { } - reads = channel.fromPath(params.fastq, checkIfExists:true) - results = pipeline(reads) - output(results) + reads = file("$params.fastq/*.fastq*", type: 'file', maxdepth: 1) + if (reads) { + reads = Channel.fromPath(params.fastq, type: 'dir', maxDepth: 1) + results = pipeline(reads) + output(results) + } else { + println("No .fastq(.gz) files found under `${params.fastq}`.") + } } diff --git a/nextflow.config b/nextflow.config index 705087c..0064080 100644 --- a/nextflow.config +++ b/nextflow.config @@ -36,7 +36,7 @@ profiles { } process { withLabel:pysam { - container = "ontresearch/workflow-template:${params.wfversion}" + container = "ontresearch/wf-template:${params.wfversion}" } shell = ['/bin/bash', '-euo', 'pipefail'] } diff --git a/test_data/reads.fastq.gz b/test_data/reads.fastq.gz new file mode 100644 index 0000000..b2eac34 Binary files /dev/null and b/test_data/reads.fastq.gz differ diff --git a/test_data/reads.fq.gz b/test_data/reads.fq.gz deleted file mode 100644 index f4dc3d5..0000000 Binary files a/test_data/reads.fq.gz and /dev/null differ diff --git a/test_data/report.html b/test_data/report.html new file mode 100644 index 0000000..94ecc1d --- /dev/null +++ b/test_data/report.html @@ -0,0 +1,929 @@ + + + + + + + Workflow Template Sequencing report + + + + + + + + + + + + +
+

Workflow Template Sequencing report

+

Results generated through the wf-template nextflow workflow by Oxford Nanopore Technologies +

Read Quality control

+

This section displays basic QC metrics indicating read data quality.

+ +
+

About

+

Oxford Nanopore Technologies products are not intended for use for health assessment +or to diagnose, treat, mitigate, cure or prevent any disease or condition.

+

This report was produced using the epi2me-labs/wf-artic. +The workflow can be run using nextflow epi2me-labs/wf-artic --help

+
+ + + \ No newline at end of file