Add larger test data and sample report

This commit is contained in:
Chris Wright 2021-03-01 11:44:41 +00:00
parent d79c23ada4
commit ed1fc2c9fd
8 changed files with 1071 additions and 17 deletions

View File

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

88
bin/report.py Executable file
View File

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

View File

@ -4,5 +4,6 @@ channels:
- conda-forge
- defaults
dependencies:
- python==3.6
- python==3.6.*
- aplanat
- pysam

32
main.nf
View File

@ -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}`.")
}
}

View File

@ -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']
}

BIN
test_data/reads.fastq.gz Normal file

Binary file not shown.

Binary file not shown.

929
test_data/report.html Normal file

File diff suppressed because one or more lines are too long