rename report

bin/read_lengths.py

@@ -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()

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__":

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

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)