diff --git a/bin/read_lengths.py b/bin/read_lengths.py deleted file mode 100755 index 669ed47..0000000 --- a/bin/read_lengths.py +++ /dev/null @@ -1,50 +0,0 @@ -#!/usr/bin/env python -"""Create a simple summary of a fastq file.""" - -import argparse -import glob -import itertools -import os - -import numpy as np -import pysam - - -def mean_qual(quals): - """Calculate mean quality of a read.""" - 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(): - """Run entry point.""" - parser = argparse.ArgumentParser() - parser.add_argument( - "fastq", help="Directory containing .fastq(.gz) files, or single fastq") - parser.add_argument( - "output", help="Output file") - args = parser.parse_args() - - if os.path.isfile(args.fastq): - fastqs = [args.fastq] - elif os.path.isdir(args.fastq): - fastqs = glob.glob(os.path.join(args.directory, "*.fastq*")) - else: - raise IOError("fastq argument should be directory of file.") - reads = itertools.chain.from_iterable( - pysam.FastxFile(fname) for fname in fastqs) - - 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 index 14fae2a..c8d85c5 100755 --- a/bin/report.py +++ b/bin/report.py @@ -3,18 +3,8 @@ import argparse -from aplanat import annot, hist, report -from bokeh.layouts import gridplot -import numpy as np -import pandas as pd - - -def read_files(summaries): - """Combine a list of files into a single dataframe.""" - dfs = list() - for fname in sorted(summaries): - dfs.append(pd.read_csv(fname, sep="\t")) - return pd.concat(dfs) +from aplanat.components import fastcat +from aplanat.report import HTMLReport def main(): @@ -24,52 +14,15 @@ def main(): parser.add_argument("summaries", nargs='+', help="Read summary file.") args = parser.parse_args() - report_doc = report.HTMLReport( + 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. -''') + report.add_section( + section=fastcat.full_report(args.summaries)) - np_blue = '#0084A9' - - # 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(''' + report.markdown(''' ### About **Oxford Nanopore Technologies products are not intended for use for health @@ -84,7 +37,7 @@ workflow can be run using `nextflow epi2me-labs/wf-template --help` ''') # write report - report_doc.write(args.report) + report.write(args.report) if __name__ == "__main__": diff --git a/environment.yaml b/environment.yaml index 4860b82..873af7e 100644 --- a/environment.yaml +++ b/environment.yaml @@ -1,9 +1,11 @@ name: epi2melabs-nf-template-workflow channels: + - epi2melabs - bioconda - conda-forge - defaults dependencies: - python==3.6.* - - aplanat + - aplanat >=0.3.5 - pysam + - fastcat diff --git a/main.nf b/main.nf index fc3fec8..dac286e 100644 --- a/main.nf +++ b/main.nf @@ -27,50 +27,18 @@ Script Options: } -process concatFastq { +process summariseReads { // concatenate fastq and fastq.gz in a dir label "pysam" cpus 1 input: file "input" - output: - file "reads.fastq.gz" - - shell: - ''' -#!/usr/bin/env python -from glob import glob -import gzip -import itertools -import os -import pysam - -# we use pysam just because it will read both fastq and fastq.gz -# and we don't have to worry about having a combination or not -with gzip.open("reads.fastq.gz", "wt") as fh: - files = itertools.chain( - glob("input/*.fastq"), glob("input/*.fastq.gz")) - records = itertools.chain.from_iterable( - pysam.FastxFile(fn) for fn in files) - for rec in records: - annot = " {}".format(rec.comment) if rec.comment else "" - qual = rec.quality if rec.quality else "+"*len(rec.sequence) - fh.write("@{}{}\\n{}\\n+\\n{}\\n".format(rec.name, annot, rec.sequence, qual)) - ''' -} - - -process readSeqs { - // Just write a file with sequence lengths - label "pysam" - input: - file reads output: file "seqs.txt" - + shell: """ - read_lengths.py $reads seqs.txt + fastcat -r seqs.txt input/*.fastq* > /dev/null """ } @@ -80,9 +48,9 @@ process makeReport { input: file "seqs.txt" output: - file "report.html" + file "wf-template-report.html" """ - report.py report.html seqs.txt + report.py wf-template-report.html seqs.txt """ } @@ -109,8 +77,7 @@ workflow pipeline { take: reads main: - reads = concatFastq(reads) - summary = readSeqs(reads) + summary = summariseReads(reads) report = makeReport(summary) emit: summary.concat(report)