+
{{ title }}
+
{{ lead }}
+ {{ div }}
+
+ {{ big_table }}
+
+
+
diff --git a/bin/run_fastq_qc.py b/bin/run_fastq_qc.py
deleted file mode 100755
index 2512a06..0000000
--- a/bin/run_fastq_qc.py
+++ /dev/null
@@ -1,59 +0,0 @@
-#!/usr/bin/env python
-"""Get fastq QC reports."""
-
-# -*- coding: utf-8 -*-
-
-import argparse
-import os
-import sys
-
-import numpy as np
-from pysam import FastxFile
-
-
-def parse_args(argv=sys.argv[1:]):
- """Parse args."""
- description = """Script to run the isoform workflow """
- parser = argparse.ArgumentParser(description=description)
- parser.add_argument("--fastq", required=True, help="")
- parser.add_argument("--output_dir", required=True, help="")
- return parser.parse_args(argv)
-
-
-def run_fastq_qc(fastq_path, output):
- """Write QC info to files."""
- qualities = list()
- mean_qualities = list()
- lengths = list()
-
- with FastxFile(fastq_path) as fq:
- for rec in fq:
- # ONT calculation for "mean Q score"
- quals = np.fromiter(
- (ord(x) - 33 for x in rec.quality),
- dtype=int, count=len(rec.quality))
- mean_p = np.mean(np.power(10, quals / -10))
- mean_qualities.append(-10 * np.log10(mean_p))
- # all qualities
- qualities.extend(quals)
- lengths.append(len(quals))
-
- with open(os.path.join(output, "base_qual.txt"), 'w') as f:
- f.write("\n".join((str(q) for q in qualities)))
-
- with open(os.path.join(output, "read_qual.txt"), 'w') as f:
- f.write("\n".join((str(q) for q in mean_qualities)))
-
- with open(os.path.join(output, "lengths.txt"), 'w') as f:
- f.write("\n".join((str(_l) for _l in lengths)))
-
-
-def main(args):
- """Run entry point."""
- assert os.path.isfile(args.fastq)
- assert os.path.isdir(args.output_dir)
- run_fastq_qc(fastq_path=args.fastq, output=args.output_dir)
-
-
-if __name__ == '__main__':
- main(args=parse_args())
diff --git a/bin/run_isonclust2.py b/bin/run_isonclust2.py
new file mode 100755
index 0000000..7560dbc
--- /dev/null
+++ b/bin/run_isonclust2.py
@@ -0,0 +1,124 @@
+#!/usr/bin/env python
+"""Dynamically generate isONclust2 processes."""
+from collections import OrderedDict
+from glob import glob
+from itertools import zip_longest
+from pathlib import Path
+import re
+import subprocess as sub
+
+
+class Node:
+ """Node."""
+
+ def __init__(self, Id, File, Left, Right, Parent, Level):
+ """Set node attaributes."""
+ self.Id = Id
+ self.File = File
+ self.Left = Left
+ self.Right = Right
+ self.Parent = Parent
+ self.Level = Level
+ self.Done = False
+ self.RightSide = False
+
+ def __repr__(self):
+ """Get string repr of a node."""
+ return "Node:{} Level: {} File: {} Done: {} Left: {} Right: " \
+ "{} Parent: {}".format(
+ self.Id, self.Level,
+ self.File, self.Done, self.Left.Id if
+ self.Left is not None else None,
+ self.Right.Id if self.Right is not None else None,
+ self.Parent.Id if self.Parent is not None else None)
+
+
+def grouper(n, iterable, fillvalue=None):
+ """
+ Group adjacent nodes.
+
+ grouper(3, 'ABCDEFG', 'x') --> ABC DEF Gxx.
+ """
+ args = [iter(iterable)] * n
+ return zip_longest(fillvalue=fillvalue, *args)
+
+
+def build_job_tree():
+ """Build a job tree of nodes."""
+ JOB_TREE = OrderedDict()
+ batches = glob("batches/isONbatch_*.cer")
+ batch_ids = [int(re.search('batches/isONbatch_(.*)\\.cer$', x).group(1))
+ for x in batches]
+ LEVELS = OrderedDict()
+ LEVELS[0] = []
+ for Id, bf in sorted(zip(batch_ids, batches), key=lambda x: x[0]):
+ n = Node(
+ Id,
+ "clusters/isONcluster_{}.cer".format(Id),
+ None,
+ None,
+ None,
+ 0)
+ n.Done = True
+ JOB_TREE[Id] = n
+ LEVELS[0].append(n)
+ level = 0
+ max_id = LEVELS[0][-1].Id
+ while len(LEVELS[level]) != 1: # Final level will be link
+ next_level = level + 1
+ LEVELS[next_level] = []
+ for l_, r in grouper(2, LEVELS[level]):
+ if r is None: # End of a level
+ LEVELS[level].pop() # remove last node?
+ l_.Level += 1 # ncrement level
+ LEVELS[next_level].append(l_) # Add the left to the next level
+ continue
+ max_id += 1
+ new_batch = "clusters/isONcluster_{}.cer".format(max_id)
+ new_node = Node(max_id, new_batch, l_, r, None, next_level)
+ l_.Parent = new_node
+ r.Parent = new_node
+ r.RightSide = True
+ LEVELS[next_level].append(new_node)
+ JOB_TREE[max_id] = new_node
+ level = next_level
+ ROOT = JOB_TREE[len(JOB_TREE) - 1].Id
+ JOB_TREE[ROOT].RightSide = True
+
+ return JOB_TREE, LEVELS
+
+
+def main():
+ """Entry point."""
+ Path('clusters').mkdir(exist_ok=True)
+ job_tree, levels = build_job_tree()
+
+ init_template = 'isONclust2 cluster -x {} -v -Q -l batches/' \
+ 'isONbatch_{}.cer -o clusters/isONcluster_{}.cer {}; ' \
+ 'sync;\n'
+ template = 'isONclust2 cluster -x {} -v -Q -l clusters/isONcluster_{}' \
+ '.cer -r clusters/isONcluster_{}.cer -o clusters/isONcluster' \
+ '_{}.cer {}; sync\n'
+
+ for nr, l in levels.items():
+ jobs_out = 'jobs_level_{}.sh'.format(nr)
+ with open(jobs_out, 'w') as fh:
+ for n in l:
+ purge = "-z" if n.RightSide else ""
+ if nr == 0 or n.Left is None or n.Right is None:
+ jr = init_template.format('sahlin', n.Id, n.Id, purge)
+ fh.write(jr)
+ else:
+ jr = template.format('sahlin', n.Left.Id,
+ n.Right.Id, n.Id, purge)
+ fh.write(jr)
+ # Run a level in parallel
+ cmd = "parallel < {}".format(jobs_out)
+ sub.call(cmd, shell=True)
+ sub.call("ln -s `realpath clusters/isONcluster_{}"
+ ".cer` isONcluster_ROOT.cer".format(n.Id), shell=True)
+
+
+if __name__ == '__main__':
+ # The cwd should be the process dir that contains 'batches/'
+ main()
diff --git a/bin/table_template.html b/bin/table_template.html
new file mode 100755
index 0000000..1faf6d9
--- /dev/null
+++ b/bin/table_template.html
@@ -0,0 +1,106 @@
+
+
+
+
+
+
Table loading
+{{ dataframe }}
+
+
+
diff --git a/data/OPTIONAL_FILE_1 b/data/OPTIONAL_FILE_1
deleted file mode 100644
index 65aa70a..0000000
--- a/data/OPTIONAL_FILE_1
+++ /dev/null
@@ -1,11 +0,0 @@
-# Nothing to see here. A sentinel file to replace real data.
-# e.g.:
-#
-# input:
-# file some_data
-# file extra_data
-# script:
-# def extra = extra_data.name != 'OPTIONAL_FILE' ? "--extra-data $opt" : ''
-# """
-# command ${some_data} ${extra}
-# """
diff --git a/denovo.nf b/denovo.nf
deleted file mode 100644
index 3d08111..0000000
--- a/denovo.nf
+++ /dev/null
@@ -1,205 +0,0 @@
-import sys
-import time
-import os
-from os import path
-import pandas as pd
-from glob import glob
-import re
-from threading import Lock
-from itertools import zip_longest
-from collections import OrderedDict, namedtuple
-
-if not workflow.overwrite_configfiles:
- configfile: "config.yml"
-
-WORKDIR = path.abspath(path.join(config["workdir_top"], config["pipeline"]))
-workdir: WORKDIR
-SNAKEDIR = path.dirname(workflow.snakefile)
-
-include: "snakelib/utils.snake"
-
-in_fastq = config["reads_fastq"]
-if not path.isabs(in_fastq):
- in_fastq = path.join(SNAKEDIR, in_fastq)
-
-
-class Node:
- def __init__(self, Id, File, Left, Right, Parent, Level):
- self.Id = Id
- self.File = File
- self.Left = Left
- self.Right = Right
- self.Parent = Parent
- self.Level = Level
- self.Done = False
- self.RightSide = False
- def __repr__(self):
- return "Node:{} Level: {} File: {} Done: {} Left: {} Right: {} Parent: {}".format(self.Id, self.Level, self.File, self.Done, self.Left.Id if self.Left is not None else None, self.Right.Id if self.Right is not None else None, self.Parent.Id if self.Parent is not None else None)
-
-def grouper(n, iterable, fillvalue=None):
- "grouper(3, 'ABCDEFG', 'x') --> ABC DEF Gxx"
- args = [iter(iterable)] * n
- return zip_longest(fillvalue=fillvalue, *args)
-
-def build_job_tree(batch_dir):
- batches = glob("{}/isONbatch_*.cer".format(batch_dir))
- batch_ids = [int(re.search('/isONbatch_(.*)\.cer$', x).group(1)) for x in batches]
- global LEVELS
- global JOB_TREE
- LEVELS = OrderedDict()
- LEVELS[0] = []
- for Id, bf in sorted(zip(batch_ids, batches), key=lambda x: x[0]):
- n = Node(Id, "clusters/isONcluster_{}.cer".format(Id), None, None, None, 0)
- n.Done = True
- JOB_TREE[Id] = n
- LEVELS[0].append(n)
- level = 0
- max_id = LEVELS[0][-1].Id
- while len(LEVELS[level]) != 1:
- next_level = level + 1
- LEVELS[next_level] = []
- for l, r in grouper(2, LEVELS[level]):
- if r is None:
- LEVELS[level].pop()
- l.Level += 1
- LEVELS[next_level].append(l)
- continue
- max_id += 1
- new_batch = "clusters/isONcluster_{}.cer".format(max_id)
- new_node = Node(max_id, new_batch, l, r, None, next_level)
- l.Parent = new_node
- r.Parent = new_node
- r.RightSide = True
- LEVELS[next_level].append(new_node)
- JOB_TREE[max_id] = new_node
- level = next_level
- global ROOT
- ROOT = JOB_TREE[len(JOB_TREE)-1].Id
- JOB_TREE[ROOT].RightSide = True
- print("Merge clustering job tree nodes:",file=sys.stderr)
- for n in JOB_TREE.values():
- print("\t{}".format(n),file=sys.stderr)
-
-
-def generate_rules(levels, snk):
- init_template = """
-rule cluster_job_%d:
- input:
- left = "sorted/batches/isONbatch_%d.cer",
- output: "clusters/isONcluster_%d.cer"
- shell: "isONclust2 cluster -x %s -v -Q -l %s -o %s %s; sync"
-
- """
- template = """
-rule cluster_job_%d:
- input:
- left = "clusters/isONcluster_%d.cer",
- right = "clusters/isONcluster_%d.cer",
- output: "clusters/isONcluster_%d.cer"
- shell: "isONclust2 cluster -x %s -v -Q -l %s -r %s -o %s %s; sync"
-
- """
-
- link_template="""
-rule link_root:
- input: "clusters/isONcluster_%d.cer",
- output: "clusters/isONcluster_ROOT.cer"
- shell: "ln -s `realpath {input}` {output}"
-
- """
-
- fh = open(snk, "w")
- for nr, l in levels.items():
- for n in l:
- purge = "-z" if n.RightSide else ""
- if nr == 0 or n.Left is None or n.Right is None:
- jr = init_template % (n.Id, n.Id, n.Id, config["cls_mode"], "{input.left}", "{output}", purge)
- fh.write(jr)
- else:
- jr = template % (n.Id, n.Left.Id, n.Right.Id, n.Id, config["cls_mode"], "{input.left}", "{input.right}", "{output}", purge)
- fh.write(jr)
- global ROOT
- fh.write(link_template % ROOT)
- fh.write("\nROOT = %d" % ROOT)
- fh.flush()
- fh.close()
-
-def count_fastq_bases(fname, size=128000000):
- fh = open(fname, "r")
- count = 0
- while True:
- b = fh.read(size)
- if not b:
- break
- count += b.count("A")
- count += b.count("T")
- count += b.count("G")
- count += b.count("C")
- count += b.count("U")
- fh.close()
- return count
-
-def preprocess_reads(fq):
- pc_opts = config["pychopper_opts"]
- concat = config["concatenate"]
- thr = config["cores"]
-
- out_fq = "processed_reads/full_length_reads.fq"
- if os.path.isdir("processed_reads"):
- return out_fq
-
- shell("mkdir -p processed_reads")
- if concat:
- print("Concatenating reads under directory: " + fq)
- shell("find %s -regextype posix-extended -regex '.*\.(fastq|fq)$' -exec cat {{}} \\; > processed_reads/input_reads.fq" % fq)
- else:
- shell("ln -s `realpath %s` processed_reads/input_reads.fq" % fq)
-
- if config["run_pychopper"]:
- print("Running pychopper of fastq file: processed_reads/input_reads.fq")
- shell("(cd processed_reads; cdna_classifier.py -t %d %s input_reads.fq full_length_reads.fq)" % (thr, pc_opts))
- else:
- shell("ln -s `realpath processed_reads/input_reads.fq` processed_reads/full_length_reads.fq")
-
- return out_fq
-
-
-ROOT = None
-DYNAMIC_RULES="job_rules.snk"
-if ((not os.path.isfile(os.path.join(WORKDIR,"sorted","sorted_reads.fastq"))) or (not os.path.isfile(os.path.join(SNAKEDIR, DYNAMIC_RULES)))):
- print("Preprocessing read in fastq file:", in_fastq)
- proc_fastq = preprocess_reads(in_fastq)
- print("Counting records in input fastq:", proc_fastq)
- nr_bases = count_fastq_bases(proc_fastq)
- print("Bases in input: {} megabases".format(int(nr_bases/10**6)))
- if config['batch_size'] < 0:
- config['batch_size'] = int(nr_bases/1000/config["cores"])
-
- print("Batch size is: {}".format(config['batch_size']))
-
- init_cls_options = """ --batch-size {} --kmer-size {} --window-size {} --min-shared {} --min-qual {}\
- --mapped-threshold {} --aligned-threshold {} --min-fraction {} --min-prob-no-hits {} -M {} -P {} -g {} -c {} -F {} """
- init_cls_options = init_cls_options.format(config["batch_size"], config["kmer_size"], config["window_size"], config["min_shared"], config["min_qual"], \
- config["mapped_threshold"], config["aligned_threshold"], config["min_fraction"], config["min_prob_no_hits"], config["batch_max_seq"], config["consensus_period"],
- config["consensus_minimum"], config["consensus_maximum"], config["min_left_cls"])
-
- shell("""
- rm -fr clusters sorted
- mkdir -p sorted; isONclust2 sort {} -v -o sorted {};
- mkdir -p clusters;
- """.format(init_cls_options, proc_fastq))
- JOB_TREE = OrderedDict()
- LEVELS = None
-
- build_job_tree("sorted/batches")
- generate_rules(LEVELS, "{}/job_rules.snk".format(SNAKEDIR))
-
-include: DYNAMIC_RULES
-
-
-
-rule all:
- input: rules.link_root.output
- output: directory("final_clusters")
- shell:
- """ isONclust2 dump -v -i sorted/sorted_reads_idx.cer -o final_clusters {input}; sync """
\ No newline at end of file
diff --git a/denovo_assembly.nf b/denovo_assembly.nf
new file mode 100644
index 0000000..03abae1
--- /dev/null
+++ b/denovo_assembly.nf
@@ -0,0 +1,295 @@
+import groovy.json.JsonSlurper
+import nextflow.util.BlankSeparatedList;
+
+
+map_sample_ids_cls = {it ->
+/* Harmonize tuples
+output:
+ tuple val(sample_id), path('*.gff')
+When there are multiple paths, will emit:
+ [sample_id, [path, path ..]]
+when there's a single path, this:
+ [sample_id, path]
+This closure makes both cases:
+ [[sample_id, path][sample_id, path]].
+*/
+ if (it[1].getClass() != java.util.ArrayList){
+ // If only one path, `it` will be [sample_id, path]
+ return [it]
+ }
+ l = [];
+ for (x in it[1]){
+ l.add(tuple(it[0], x))
+ }
+ return l
+}
+
+
+process dump_clusters {
+ label "isoforms"
+ input:
+ tuple val(sample_id), path(root_cluster), path(sorted_reads_dir)
+ output:
+ tuple val(sample_id), path("final_clusters"), emit: final_clusters_dir
+ tuple val(sample_id), path("final_clusters/cluster_fastq/*.fq"), emit: final_clusters
+ shell:
+ """ isONclust2 dump -v -i $sorted_reads_dir/sorted_reads_idx.cer -o final_clusters $root_cluster; sync """
+
+}
+
+process build_backbones {
+ /*
+ This step can fail in what seems to at the racon stage giving an 'empty overlap error set' message
+ As a temporary fix, do racon_cmd || true to prevent it crashing the pipeline, and then move on to next cluster
+
+ This process needs some work
+ */
+ label "isoforms"
+ input:
+ tuple val(sample_id), path(cluster_fq)
+ output:
+ tuple val(sample_id), path("*final_polished_cds.fa"), emit: polished_cds
+ script:
+ def cluster_fq_bl = new BlankSeparatedList(cluster_fq)
+ """
+ # Get one of the cluster ids to give the output a unique name
+ UNID=\$(echo ${cluster_fq_bl[1]} | grep -o -E '[0-9]+')
+
+ for cluster in $cluster_fq_bl
+ do
+ clfq=`basename \$cluster`
+ cln=\${clfq%.*}
+
+ echo Building backbone for cluster: \$cln
+ echo "\tSampling input reads for backbone construction."
+ sample=\${cln}_sample.fq
+ seqkit head --quiet -n 100 \$cluster > \$sample
+ seqkit sample --quiet -n 500 -2 -s 100 \$cluster >> \$sample
+ echo "\tConstructing spoa consensus."
+ spoa_cons=\${cln}_spoa.fa
+ spoa -m 5 -n -4 -g -8 -e -6 -q -10 -c -15 -l 1 -r 0 \$sample > \$spoa_cons
+
+ echo "\tPolishing the consensus using racon."
+
+ # polish 1
+ samgz=\${cln}_aln.sam.gz
+ racon_cons=\${cln}_racon.fa
+ racon_cons1=\${cln}_racon1.fa
+
+ minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$spoa_cons \$sample | gzip - > \$samgz
+ racon -t ${params.threads} --no-trimming -u -w 2000 \$sample \$samgz \$spoa_cons > \$racon_cons || true
+ if [ -f \$racon_cons ];
+ then
+ cat \$racon_cons | seqkit replace -p ".*" -r cluster_\${cln} > \$racon_cons1
+ else
+ continue
+ fi
+
+ # polish 2
+ samgz1=\${cln}_aln_1.sam.gz
+ racon_cons2=\${cln}_racon2.fa
+ minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons1 \$sample | gzip - > \$samgz1
+ racon -t ${params.threads} -u --no-trimming \$sample \$samgz1 \$racon_cons1 > \$racon_cons2 || true
+ if [ -f \$racon_cons2 ];
+ then
+ echo "success polish 2"
+ else
+ continue
+ fi
+
+ # polish 3
+ samgz2=\${cln}_aln_2.sam.gz
+ minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons2 \$sample | gzip - > \$samgz2
+ EXITCODE=0
+ racon -t ${params.threads} -u \$sample \$samgz2 \$racon_cons2 >> ${sample_id}_\${UNID}_final_polished_cds.fa || EXITCODE=\$?
+ if [ \$EXITCODE -eq 0 ];
+ then
+ echo "finished"
+ fi
+
+ done
+ echo "Finished backbones"
+ """
+
+
+}
+
+process merge_cds {
+ label "isoforms"
+ input:
+ tuple val(sample_id), path(cds)
+ output:
+ tuple val(sample_id), path("${sample_id}_cds.fa"), emit: final_polished_cds
+ script:
+ def merge_list
+ """
+ for FILE in *final_polished_cds.fa
+ do
+ cat \$FILE >> "${sample_id}_cds.fa"
+ done
+ """
+}
+
+process cds_align {
+ label "isoforms"
+ input:
+ tuple val(sample_id), path(polished_cds), path(sorted_reads_dir)
+ output:
+ tuple val(sample_id), path("${sample_id}_reads_aln_sorted.bam"), emit: bam
+ tuple val(sample_id), path("${sample_id}_read_aln_stats.tsv"), emit: stats
+ script:
+ def ab = "${sample_id}_reads_aln_sorted.bam"
+ """
+ minimap2 -t ${params.threads} \
+ -ax splice ${params.minimap2_opts} $polished_cds $sorted_reads_dir/sorted_reads.fastq |\
+ samtools view -b - |\
+ samtools sort -o $ab;
+ samtools index $ab;
+ ((seqkit bam -s -j ${params.threads} ${ab} 2>&1) | tee ${sample_id}_read_aln_stats.tsv ) || true
+ """
+}
+
+
+process make_batches {
+ /*
+ Take a fasta file and creates batches for isONclust2 to work with
+
+ */
+ label "isoforms"
+ input:
+ tuple val(sample_id), path(fastq)
+ output:
+ tuple val(sample_id), path('sorted/batches'), emit: sorted_batches
+ tuple val(sample_id), path('sorted'), emit: sorted_reads_dir
+
+ script:
+
+ maxcpus = Runtime.runtime.availableProcessors()
+
+ minimum_batch_size = 2000
+ """
+ nr_bases=\$(seqkit stats -T $fastq|cut -f 5| sed '2q;d')
+
+ b=0
+ if [ ${params.batch_size} -lt \$b ];
+ then
+ nr_bases=\$(seqkit stats -T $fastq|cut -f 5| sed '2q;d')
+ let batch_size=\$nr_bases/1000/$maxcpus
+ if [ \$batch_size -lt $minimum_batch_size ];
+ then
+ batch_size=$minimum_batch_size
+ fi
+ else
+ batch_size=${params.batch_size}
+ fi
+
+ echo "Batch size:\$batch_size";
+ echo "Num bases: \$nr_bases";
+
+ init_cls_options="--batch-size \$batch_size --kmer-size ${params.kmer_size} \
+ --window-size ${params.window_size} --min-shared ${params.min_shared} --min-qual ${params.min_qual} \
+ --mapped-threshold ${params.mapped_threshold} --aligned-threshold ${params.aligned_threshold} \
+ --min-fraction ${params.min_fraction} --min-prob-no-hits ${params.min_prob_no_hits} \
+ -M ${params.batch_max_seq} -P ${params.consensus_period} -g ${params.consensus_minimum} -c ${params.consensus_maximum} -F ${params.min_left_cls} "
+ mkdir -p sorted; isONclust2 sort \$init_cls_options -v -o sorted $fastq;
+ """
+
+}
+
+process clustering() {
+ label "isoforms"
+
+ input:
+ tuple val(sample_id), path(sorted_batches)
+ output:
+ tuple val(sample_id), path('isONcluster_ROOT.cer'), emit: root_cluster
+ script:
+ """
+ run_isonclust2.py $sorted_batches
+
+ """
+}
+
+process cluster_quality() {
+ // Run Kristoffer Sahlin's QC code.
+ // For now just write out the PDF and CSV results. Move these to the report at some point
+ label "isoforms"
+
+ input:
+ tuple val(sample_id), path(reads_fl), path(final_clusters_dir)
+ path reference
+ output:
+ tuple val(sample_id),
+ path("${sample_id}_cluster_qc"), emit: cluster_qc_dir
+ tuple val(sample_id),
+ path("${sample_id}_cluster_qc_raw"), emit: cluster_qc_raw
+ script:
+ def qc_dir = "${sample_id}_cluster_qc"
+ def qc_dir_raw = "${sample_id}_cluster_qc_raw" // To generate plots in report
+ def bam = "${qc_dir}/ref_aln.bam"
+
+ """
+ echo $qc_dir_raw
+ mkdir $qc_dir
+ mkdir $qc_dir_raw
+ minimap2 -ax splice -t 2 $reference $reads_fl |\
+ samtools view -q 2 -F 2304 -b - |\
+ samtools sort - -o $bam;
+ samtools index $bam;
+ compute_cluster_quality.py --sizes $final_clusters_dir/clusters_info.tsv \
+ --outfile ${qc_dir}/cluster_quality.csv --ont --clusters $final_clusters_dir/clusters.tsv \
+ --classes $bam --report ${qc_dir}/cluster_quality.pdf --raw_data_out $qc_dir_raw
+ """
+
+}
+
+workflow denovo_assembly {
+ take:
+ fastq_reads_fl
+ reference
+ main:
+
+ make_batches(fastq_reads_fl)
+
+ clustering(make_batches.output.sorted_batches)
+
+ dump_clusters(clustering.output.root_cluster
+ .join(make_batches.output.sorted_reads_dir))
+
+ build_backbones(dump_clusters.output.final_clusters
+ .flatMap(map_sample_ids_cls)
+ .groupTuple(size: 10, remainder: true)
+ )
+
+ merge_cds(build_backbones.output.polished_cds
+ .flatMap(map_sample_ids_cls)
+ .groupTuple()
+ )
+
+ cds_align(merge_cds.out.final_polished_cds.view()
+ .join(make_batches.output.sorted_reads_dir))
+
+ if (!reference.name.startsWith('OPTIONAL_FILE')){
+
+ cluster_quality(fastq_reads_fl
+ .join(dump_clusters.output.final_clusters_dir), reference)
+
+ cluster_quality.out.cluster_qc_dir
+ .set { opt_qual_ch }
+
+ cluster_quality.output.cluster_qc_raw
+ .set { opt_qual_raw_ch }
+
+ } else{
+ Channel.empty().set { opt_qual_ch }
+ Channel.empty().set { opt_qual_raw_ch }
+ }
+
+ emit:
+ bam = cds_align.output.bam
+ cds = merge_cds.out.final_polished_cds
+ stats = cds_align.output.stats
+ opt_qual_ch
+ opt_qual_raw_ch
+}
+
diff --git a/environment.yaml b/environment.yaml
index 75b28c9..d771bf9 100644
--- a/environment.yaml
+++ b/environment.yaml
@@ -1,4 +1,4 @@
-name: epi2melabs-wf-isoforms
+name: epi2melabs-wf-isoforms-test
channels:
- epi2melabs
- bioconda
@@ -8,20 +8,21 @@ dependencies:
- python==3.8.*
- aplanat >=0.5.0
- epi2melabs
- - minimap2
- - samtools
- - bedtools
- - pychopper
- - pandas
- - seaborn
- - requests
- - gffread
- - seqkit
- - csvtk
+ - minimap2 ==2.24
+ - samtools ==1.14
+ - bedtools ==2.30.0
+ - pychopper==2.5.0
+ - pandas==1.3.5
+ - gffread==0.12.7
+ - gffcompare==0.11.2
+ - gffutils=0.10.1
+ - seqkit==2.1.0
- stringtie==2.1.1
- - gffcompare
- curl
- - pysam
+ - pysam==0.17.0
+ - racon==1.4.20
- spoa==3.4.0
- - fastcat
- - isonclust2
+ - fastcat==0.4.10
+ - isonclust2==2.3
+ - parallel
+ - scikit-learn==1.0.2
diff --git a/evaluation/run_evaluation_dmel.sh b/evaluation/run_evaluation_dmel.sh
new file mode 100644
index 0000000..5da39ea
--- /dev/null
+++ b/evaluation/run_evaluation_dmel.sh
@@ -0,0 +1,59 @@
+#!/usr/bin/env bash
+
+# Usage: ./run_evaluation_dmel.sh pathto/outputdir
+
+# See the isONcorrect paper https://www.nature.com/articles/s41467-020-20340-8 where this dataset is described
+
+
+if [[ "$#" -lt 1 ]]; then
+ echo "usage: run_evaluation_dmel.sh
[nextflow.config]"
+ exit 1
+fi
+
+if [[ "$#" -eq 1 ]]; then
+ config=''
+fi
+
+if [[ "$#" -eq 2 ]]; then
+ config="-c $2";
+fi
+
+OUTDIR=$1;
+
+FASTQ_URL="http://ftp.sra.ebi.ac.uk/vol1/fastq/ERR358/005/ERR3588905/ERR3588905_1.fastq.gz"
+REF_URL="http://ftp.ensembl.org/pub/release-99/fasta/drosophila_melanogaster/dna/Drosophila_melanogaster.BDGP6.28.dna.toplevel.fa.gz"
+GFF_URL="http://ftp.ensembl.org/pub/release-99/gff3/drosophila_melanogaster/Drosophila_melanogaster.BDGP6.28.99.gff3.gz"
+
+DATA_DIR="$OUTDIR/data"
+READS_DIR="$DATA_DIR/reads"
+FASTQ="$READS_DIR/ERR3588905_1.fastq.gz"
+REF="$DATA_DIR/Drosophila_melanogaster.BDGP6.28.dna.toplevel.fa"
+GFF="$DATA_DIR/Drosophila_melanogaster.BDGP6.28.99.gff3"
+
+mkdir -p $READS_DIR
+
+if [ ! -f $REF ];
+then (echo "downloading reference genome"; cd $DATA_DIR; curl -L -C - -O $REF_URL); gzip -d ${REF}.gz
+fi
+
+if [ ! -f $GFF ];
+then
+ (echo "downloading reference annotation"; cd $DATA_DIR; curl -L -C - -O $GFF_URL); gzip -d ${GFF}.gz
+fi
+
+if [ ! -f $FASTQ ];
+then (echo "downloading reads"; cd $READS_DIR; curl -L -C - -O $FASTQ_URL); gzip -d ${FASTQ}.gz
+fi
+
+
+OUT_REF="$OUTDIR/ref"
+OUT_DENOVO="$OUTDIR/denovo"
+
+
+nextflow run ../ --fastq $READS_DIR $config \
+--ref_genome $REF --ref_annotation $GFF -profile local --out_dir $OUT_REF --minimap2_opts '-uf --splice-flank=no' \
+-w $OUT_REF/workspace -resume;
+
+echo "Doing de novo evaluation"
+nextflow run ../ --fastq $READS_DIR $config --denovo -profile local --out_dir $OUT_DENOVO \
+-w $OUT_DENOVO/workspace -resume;
diff --git a/evaluation/tests.sh b/evaluation/tests.sh
new file mode 100644
index 0000000..515016f
--- /dev/null
+++ b/evaluation/tests.sh
@@ -0,0 +1,73 @@
+#!/usr/bin/env bash
+
+# A few simple tests with different combinations of CLI options
+# Run from within an appropriate active conda environment
+
+if [[ "$#" -lt 1 ]]; then
+ echo "usage: tests.sh [nextflow.config]"
+ exit 1
+fi
+
+if [[ "$#" -eq 1 ]]; then
+ config=''
+fi
+
+if [[ "$#" -eq 2 ]]; then
+ config="-c $2";
+fi
+
+SCRIPT_DIR="$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )"
+cd $SCRIPT_DIR/../;
+
+singledir="test_data/fastq"
+multisampledir="test_data/demultiplexed_fastq"
+
+# This is for when using SIRV dataset with non-canonical spice junctions
+#"--minimap2_opts '-uf --splice-flank=no'"
+results=()
+
+
+OUTPUT=$1/denovo_multi_sample_no_ref_genome;
+nextflow run . --fastq $multisampledir $config --denovo --ref_genome test_data/SIRV_150601a.fasta -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
+--sample_sheet test_data/sample_sheet -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+OUTPUT=$1/denovo_single;
+nextflow run . --fastq $singledir $config --denovo --ref_genome test_data/SIRV_150601a.fasta -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
+--sample_sheet test_data/sample_sheet -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+# Reference based tests
+OUTPUT=$1/reference_single_dir;
+nextflow run . --fastq $singledir $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no' \
+--ref_annotation test_data/SIRV_isofroms.gtf -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+OUTPUT=$1/multiple_samples;
+nextflow run . --fastq $multisampledir $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no'\
+--ref_annotation test_data/SIRV_isofroms.gtf -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
+--sample_sheet test_data/sample_sheet -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+OUTPUT=$1/reference_no_ref_annotation;
+nextflow run . --fastq $singledir $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no'\
+ -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+# Force split_bam to make multiple alignment bundles
+OUTPUT=$1/reference_frce_split_bam;
+nextflow run . --fastq $singledir $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no'\
+--ref_annotation test_data/SIRV_isofroms.gtf -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
+--bundle_min_reads 5 -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+echo "Exit status codes for each test"
+for value in "${results[@]}"; do
+ echo "${value}"
+done
\ No newline at end of file
diff --git a/main.nf b/main.nf
index e833789..6cd0250 100644
--- a/main.nf
+++ b/main.nf
@@ -1,39 +1,37 @@
#!/usr/bin/env nextflow
-// Developer notes
-//
-// This template workflow provides a basic structure to copy in order
-// to create a new workflow. Current recommended pratices are:
-// i) create a simple command-line interface.
-// ii) include an abstract workflow scope named "pipeline" to be used
-// in a module fashion.
-// iii) a second concreate, but anonymous, workflow scope to be used
-// as an entry point when using this workflow in isolation.
+/* This workflow is a adapted from two previous pipeline written in Snakemake:
+- https://github.com/nanoporetech/pipeline-nanopore-ref-isoforms
+- https://github.com/nanoporetech/pipeline-nanopore-denovo-isoforms
+*/
-import groovy.json.JsonBuilder
+import groovy.json.JsonBuilder;
+import nextflow.util.BlankSeparatedList;
import java.util.ArrayList;
nextflow.enable.dsl = 2
include { fastq_ingress } from './lib/fastqingress'
include { start_ping; end_ping } from './lib/ping'
+include { reference_assembly } from './reference_assembly'
+include { denovo_assembly } from './denovo_assembly'
-process summariseReads {
- // concatenate fastq and fastq.gz in a dir
+process summariseConcatReads {
+ // concatenate fastq and fastq.gz in a dir write stats
label "isoforms"
cpus 1
input:
- tuple path(directory), val(sample_name), val(type)
+ tuple path(directory), val(sample_id), val(type)
output:
- tuple val(sample_name), path('*'), emit: summary
+ tuple val(sample_id), path("${sample_id}.fastq"), emit: input_reads
+ tuple val(sample_id), path('*.stats'), emit: summary
script:
"""
- fastcat -s ${sample_name} -r ${sample_name}.stats -x ${directory} > /dev/null
+ fastcat -s ${sample_id} -r ${sample_id}.stats -x ${directory} > ${sample_id}.fastq
"""
}
-
process getVersions {
label "isoforms"
cpus 1
@@ -44,7 +42,7 @@ process getVersions {
python -c "import pysam; print(f'pysam,{pysam.__version__}')" >> versions.txt
python -c "import aplanat; print(f'aplanat,{aplanat.__version__}')" >> versions.txt
python -c "import pandas; print(f'pandas,{pandas.__version__}')" >> versions.txt
- python -c "import seaborn; print(f'seaborn,{seaborn.__version__}')" >> versions.txt
+ python -c "import sklearn; print(f'scikit-learn,{sklearn.__version__}')" >> versions.txt
fastcat --version | sed 's/^/fastcat,/' >> versions.txt
minimap2 --version | sed 's/^/minimap2,/' >> versions.txt
samtools --version | head -n 1 | sed 's/ /,/' >> versions.txt
@@ -52,10 +50,10 @@ process getVersions {
python -c "import pychopper; print(f'pychopper,{pychopper.__version__}')" >> versions.txt
gffread --version | sed 's/^/gffread,/' >> versions.txt
seqkit version | head -n 1 | sed 's/ /,/' >> versions.txt
- csvtk version | head -n 1 | sed 's/ /,/' >> versions.txt
stringtie --version | sed 's/^/stringtie,/' >> versions.txt
gffcompare --version | head -n 1 | sed 's/ /,/' >> versions.txt
spoa --version | sed 's/^/spoa,/' >> versions.txt
+ isONclust2 version | sed 's/ version: /,/' >> versions.txt
"""
}
@@ -80,44 +78,27 @@ process preprocess_reads {
*/
label "isoforms"
- cpus params.threads
+ cpus 4
input:
- tuple path(directory), val(sample_id), val(type)
+ tuple val(sample_id), path(input_reads)
output:
- tuple val(sample_id), path("full_length_reads.fq"), emit: full_len_reads
- tuple val(sample_id), path('cdna_classifier_report.tsv'), emit: report
-// val "${sample_id}", emit: sample_id
-// path 'cdna_classifier_report.tsv', optional: true, emit: cdna_class_report
- // Not sure if this is an antipattern, but it's function is to publish all the files to the publishDir dir
+ tuple val(sample_id), path("${sample_id}_full_length_reads.fq"), emit: full_len_reads
+ tuple val(sample_id), path('*.tsv'), emit: report
script:
"""
- fastcat -s ${sample_id} -r ${sample_id}.stats -x ${directory} > input_reads.fq
-
if [[ ${params.use_pychopper} == true ]];
then
- cdna_classifier.py -t ${params.threads} ${params.pychopper_opts} input_reads.fq full_length_reads.fq
- generate_pychopper_stats.py --data cdna_classifier_report.tsv --output .
+ cdna_classifier.py -t ${params.threads} ${params.pychopper_opts} ${input_reads} ${sample_id}_full_length_reads.fq
+ mv cdna_classifier_report.tsv ${sample_id}_cdna_classifier_report.tsv
+ generate_pychopper_stats.py --data ${sample_id}_cdna_classifier_report.tsv --output .
else
- ln -s `realpath input_reads.fq` full_length_reads.fq
+ ln -s `realpath $input_reads` "${sample_id}_full_length_reads.fq"
+ touch $sample_id}_cdna_classifier_report.tsv
fi
"""
}
-process generate_fq_stats {
- label "isoforms"
-
- input:
- tuple(sample_id), path(fastq), val(unused)
-
- output:
- path "*"
- script:
- """
- run_fastq_qc.py --fastq ${fastq} --output .
- """
-}
-
process build_minimap_index{
/*
Build minimap index from reference genome
@@ -126,7 +107,7 @@ process build_minimap_index{
cpus params.threads
input:
- file reference
+ path reference
output:
path "genome_index.mmi", emit: index
script:
@@ -135,66 +116,6 @@ process build_minimap_index{
"""
}
-
-process map_reads{
- /*
- Map reads to reference using minimap2.
- Filter reads by mapping quality.
- Filter reads where length of poly(A) > max_poly_run at either ends of the read (defined by poly_context)
- */
- label "isoforms"
- cpus params.threads
-
- input:
- file index
- file reference
- tuple val(sample_id), file(fastq_reads)
-
- output:
- tuple val(sample_id), path("reads_aln_sorted.bam"), emit: bam
- tuple val(sample_id), path("read_aln_stats.tsv"), emit: stats
- script:
- def ab = "reads_aln_sorted.bam"
- def af = "internal_priming_fail.tsv"
- def fs = "context_internal_priming_fail_start.fasta"
- def fe = "context_internal_priming_fail_end.fasta"
- def fasta_reads = "reads.fa"
- def ContextFilter = """AlnContext: { Ref: "${reference}", LeftShift: -${params.poly_context}, RightShift: ${params.poly_context},
- RegexEnd: "[Aa]{${params.max_poly_run},}",
- Stranded: True,Invert: True, Tsv: "internal_priming_fail.tsv"} """
- """
- seqkit fq2fa ${fastq_reads} -o ${fasta_reads};
- minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} ${index} ${fasta_reads}\
- | samtools view -q ${params.minimum_mapping_quality} -F 2304 -Sb -\
- | seqkit bam -j ${params.threads} -x -T '${ContextFilter}' -\
- | samtools sort -@ ${params.threads} -o ${ab} -;
-
- ((seqkit bam -s -j ${params.threads} ${ab} 2>&1) | tee read_aln_stats.tsv ) || true
-
- if [[ -s ${af} ]];
- then
- tail -n +2 ${af} | awk '{{print ">" \$1 "\\n" \$4 }}' - > ${fs}
- tail -n +2 ${af} | awk '{{print ">" \$1 "\\n" \$6 }}' - > ${fe}
- fi
- """
-}
-
-process plot_aln_stats{
- /*
- Create a pdf of alignemnt statistis.
- */
- label 'isoforms'
-
- input:
- tuple val(sample_id), path (aln_stats)
- output:
- path "*"
- script:
- """
- plot_aln_stats.py ${aln_stats} -r read_aln_stats.pdf
- """
-}
-
process split_bam{
/*
Partition BAM file into loci or bundles with `params.bundle_min_reads` minimum size
@@ -215,19 +136,21 @@ process split_bam{
"""
seqkit bam -j ${params.threads} -N ${params.bundle_min_reads} ${bam} -o bam_bundles/
mv bam_bundles/* .
+ for f in *:*; do mv -v "\$f" \$(echo "\$f" | tr ':' '-'); done
"""
else
"""
mkdir -p ./${sample_id}_bam_bundles
- ln -s ${bam} ${sample_id}_bam_bundles/000000000_ALL:0:1_bundle.bam
+ ln -s ${bam} ${sample_id}_bam_bundles-000000000-ALL-0-1_bundle.bam
"""
}
-process stringtie{
+process assemble_transcripts{
/*
- Takes in aligned reads in bam format that may be a chunk of a larger alignment file.
- $G_FLAG specifies whether or not to use reference annotation as a guide in transcript assembly.
+ Assemble transcripts using stringtie.
+ Take aligned reads in bam format that may be a chunk of a larger alignment file.
+ Optionally use reference annotation to guide assembly.
Output gff annotation files in a tuple with `sample_id` for combining into samples late rin the pipeline.
*/
@@ -236,15 +159,14 @@ process stringtie{
input:
tuple val(sample_id), path(bam)
- file ref_annotation
+ path ref_annotation
output:
tuple val(sample_id), path('*.gff'), emit: gff_bundles
script:
- def out_filename = bam.name.replaceFirst(~/\.[^\.]+$/, '') + '.gff'
- def label = "STR.${bam.name.split('_')[0].toInteger()}."
+ def out_filename = bam.name.replaceFirst(~/\.[^\.]+$/, '') + "_${sample_id}.gff"
def G_FLAG = ref_annotation.name.startsWith('OPTIONAL_FILE') ? '' : "-G ${ref_annotation}"
"""
- stringtie --rf ${G_FLAG} -L -v -p ${params.threads} ${params.stringtie_opts} -o ${out_filename} \
+ stringtie --rf ${G_FLAG} -L -v -A gene_abund.tab -p ${params.threads} ${params.stringtie_opts} -o ${out_filename} \
${bam} 2>/dev/null
"""
}
@@ -261,9 +183,9 @@ process merge_gff_bundles{
output:
tuple val(sample_id), path('*.gff'), emit: gff
script:
- def merged_gff = "str_merged_${sample_id}.gff"
+ def merged_gff = "transcripts_${sample_id}.gff"
"""
- echo '#gff-version 2' >> $merged_gff;
+ echo '##gff-version 2' >> $merged_gff;
echo '#pipeline-nanopore-isoforms: stringtie' >> $merged_gff;
for fn in ${gff_bundle};
@@ -274,9 +196,11 @@ process merge_gff_bundles{
"""
}
-process run_gff_compare{
+process run_gffcompare{
/*
Compare query and reference annotations.
+ If ref_annotation is an optional file, just make an empty directory to satisfy
+ the requirements of the downstream processes.
*/
label 'isoforms'
@@ -285,45 +209,55 @@ process run_gff_compare{
tuple val(sample_id), path(query_annotation)
path ref_annotation
output:
- tuple val(sample_id), path('str_merged.annotated.gtf'), emit: merged_annotated
- tuple val(sample_id), path('str_merged.stats'), emit: stats
- tuple val(sample_id), path('str_merged.tracking'), emit: tracking
+ tuple val(sample_id), path("${sample_id}_gffcompare"), emit: gffcmp_dir
script:
- """
- echo "Doing comparison of reference annotation: ${ref_annotation} and the current annotation"
- gffcompare -o str_merged -r ${ref_annotation} ${params.gffcompare_opts} ${query_annotation}
- generate_tracking_summary.py --tracking str_merged.tracking --output_dir . --annotation ${ref_annotation}
+ def out_dir = "${sample_id}_gffcompare"
- if [[ ${params.plot_gffcmp_stats} == true ]];
- then
- plot_gffcmp_stats.py -r str_gffcmp_report.pdf -t str_merged.tracking str_merged.stats;
- fi
+ if ( ref_annotation.name.startsWith('OPTIONAL_FILE') ){
"""
+ mkdir $out_dir
+ """
+ } else {
+ """
+ mkdir $out_dir
+ echo "Doing comparison of reference annotation: ${ref_annotation} and the query annotation"
+
+ gffcompare -o ${out_dir}/str_merged -r ${ref_annotation} \
+ ${params.gffcompare_opts} ${query_annotation}
+
+ generate_tracking_summary.py --tracking $out_dir/str_merged.tracking \
+ --output_dir ${out_dir} --annotation ${ref_annotation}
+
+ mv *.tmap $out_dir
+ mv *.refmap $out_dir
+ """
+ }
}
-process run_gffread{
+process get_transcriptome{
/*
- Write out a transctiptome file based on the gff annotations.
+ Write out a transcriptome file based on the query gff annotations.
+ TODO: Do we need to touch merged_transcriptome or can we just pass it?
*/
label 'isoforms'
input:
- tuple val(sample_id), path(gff_merged), path(merged_ann_gff)
- path reference_seq
+ tuple val(sample_id), path(transcripts_gff), path(gffcmp_dir), path(reference_seq)
output:
- tuple val(sample_id), path('*.fas'), emit: transcriptome
+ tuple val(sample_id), path("*.fas"), emit: transcriptome
+
script:
- def str_transcriptome = "${sample_id}_str_transcriptome.fas"
+ def transcriptome = "${sample_id}_transcriptome.fas"
def merged_transcriptome = "${sample_id}_merged_transcriptome.fas"
"""
- gffread -g ${reference_seq} -w ${str_transcriptome} ${gff_merged}
- if [ -f ${merged_ann_gff} ]
+ gffread -g ${reference_seq} -w ${transcriptome} ${transcripts_gff}
+ if [ "\$(ls -A $gffcmp_dir)" ];
then
- gffread -F -g ${reference_seq} -w ${merged_transcriptome} ${merged_ann_gff}
- else
- touch ${merged_transcriptome}
+ echo "Yes"
+ gffread -F -g ${reference_seq} -w ${merged_transcriptome} \
+ $gffcmp_dir/str_merged.annotated.gtf
fi
"""
}
@@ -333,26 +267,45 @@ process makeReport {
label "isoforms"
input:
+ path report_template
+ path table_template
path versions
- path params
- tuple val(sample_id), path(seqs), path(aln_stats), path(gffcompare_tracking), path(gffcompare_stats),
- path(cdna_class_report)
+ path "params.json"
+ val denovo
+ tuple val(sample_ids),
+ path(seq_summaries),
+ path(aln_stats),
+ path(gffcmp_dir),
+ path(cdna_class_report),
+ path(gff_annotation)
output:
- tuple val(sample_id), path("wf-isoforms-*.html"), emit: report
+ path("wf-isoforms-*.html"), emit: report
script:
- def report_name = "wf-isoforms-${sample_id}_report.html"
- def opt_gff_track = gffcompare_tracking.name.startsWith('OPTIONAL_FILE') ? '' : "--gffcompare_tracking ${gffcompare_tracking}"
- def opt_gff_stats = gffcompare_stats.name.startsWith('OPTIONAL_FILE') ? '' : "--gffcompare_stats ${gffcompare_stats}"
+ // Convert the sample_id arrayList.
+ sids = new BlankSeparatedList(sample_ids)
+
+ def report_name = "wf-isoforms-report.html"
+ def OPT_ALN = denovo ? '' : "--alignment_stats ${aln_stats}"
+ def OPT_DENOVO = denovo ? "--denovo" : ''
"""
- report.py ${report_name} --versions ${versions} ${seqs} --params params.json \
- --alignment_stats ${aln_stats} \
- ${opt_gff_track} \
- ${opt_gff_stats} \
- --pychop_report ${cdna_class_report}
+ report.py --report $report_name \
+ --report_template $report_template \
+ --table_template $table_template \
+ --versions $versions \
+ --params params.json \
+ $OPT_ALN \
+ --pychop_report $cdna_class_report \
+ --sample_ids $sids \
+ --summaries $seq_summaries \
+ --gffcompare_dir $gffcmp_dir \
+ --gff_annotation $gff_annotation \
+ --transcript_table_cov_thresh $params.transcript_table_cov_thresh \
+ $OPT_DENOVO
+
+
"""
}
-
// See https://github.com/nextflow-io/nextflow/issues/1636
// This is the only way to publish files from a workflow whilst
// decoupling the publish from the process steps.
@@ -362,11 +315,24 @@ process output {
publishDir "${params.results_dir}/${sample_id}", mode: 'copy', pattern: "*"
input:
- tuple val(sample_id), file(fname)
+ tuple val(sample_id), path(fname)
output:
- file fname
+ path fname
"""
echo "Writing output files"
+ echo $fname
+ """
+}
+
+process output_report {
+ publishDir "${params.results_dir}", mode: 'copy', pattern: "*report.html"
+
+ input:
+ path fname
+ output:
+ path fname
+ """
+ echo "Copying report"
"""
}
@@ -377,6 +343,8 @@ workflow pipeline {
reads
ref_genome
ref_annotation
+ report_template
+ table_template
main:
map_sample_ids_cls = {it ->
@@ -401,60 +369,92 @@ workflow pipeline {
return l
}
- summariseReads(reads)
- sample_ids = summariseReads.out.summary.collect({it -> it[0]})
+ summariseConcatReads(reads)
+ sample_ids = summariseConcatReads.out.summary.flatMap({it -> it[0]})
+
software_versions = getVersions()
workflow_params = getParams()
- preprocess_reads(reads)
-// generate_fq_stats(preprocess_reads.out.sample) skip for now
+ preprocess_reads(summariseConcatReads.out.input_reads)
- build_minimap_index(ref_genome)
- map_reads(build_minimap_index.out.index, ref_genome, preprocess_reads.out.full_len_reads)
-// plot_aln_stats(map_reads.out.bam)
- split_bam(map_reads.out.bam)
- stringtie(split_bam.out.bundles.flatMap(map_sample_ids_cls), ref_annotation)
- merge_gff_bundles(stringtie.out.gff_bundles.groupTuple())
+ if (params.denovo){
+ println("Doing de novo assembly")
+ m = denovo_assembly(preprocess_reads.out.full_len_reads, ref_genome)
+
+ } else {
+ build_minimap_index(ref_genome)
+ println("Doing reference based transcript analysis")
+ m = reference_assembly(build_minimap_index.out.index, ref_genome, preprocess_reads.out.full_len_reads)
+ }
+
+ split_bam(m.bam)
+
+ assemble_transcripts(split_bam.out.bundles.flatMap(map_sample_ids_cls), ref_annotation)
+
+ merge_gff_bundles(assemble_transcripts.out.gff_bundles.groupTuple())
use_ref_ann = !ref_annotation.name.startsWith('OPTIONAL_FILE')
- if (use_ref_ann){
- run_gff_compare(merge_gff_bundles.out.gff, ref_annotation)
- run_gffread(merge_gff_bundles.out.gff.join(run_gff_compare.out.merged_annotated),
- ref_genome)
- gff_tracking = run_gff_compare.out.tracking
- gff_stats = run_gff_compare.out.stats
- }else{
- // Create dummy file paths to satisfy required path inputs of processes
- // Add to tuple with sample_id to guide to correct sample
- gff_tracking = sample_ids.combine(Channel.fromPath("$projectDir/data/OPTIONAL_FILE")).view()
- gff_stats = sample_ids.combine(Channel.fromPath("$projectDir/data/OPTIONAL_FILE_1"))
+ run_gffcompare(merge_gff_bundles.out.gff, ref_annotation)
+
+ if (params.denovo){
+ // Use the perd-sample, de novo-assembled CDS
+ seq_for_transcriptome_build = m.cds
+ }else {
+ // If doing reference based assembly, there is only one reference
+ // So map this reference to all sample_ids
+ seq_for_transcriptome_build = sample_ids.flatten().combine(Channel.fromPath(params.ref_genome))
}
- makeReport(
+ makeReport(report_template,
+ table_template,
software_versions,
workflow_params,
- summariseReads.out.summary
- .join(map_reads.out.stats)
- .join(gff_tracking)
- .join(gff_stats)
+ params.denovo,
+ summariseConcatReads.out.summary
+ .join(m.stats)
+ .join(run_gffcompare.out.gffcmp_dir)
.join(preprocess_reads.out.report)
+ .join(merge_gff_bundles.out.gff)
+ .toList().transpose().toList()
)
+
+ report = makeReport.out.report
+
+ get_transcriptome(merge_gff_bundles.out.gff
+ .join(run_gffcompare.out.gffcmp_dir)
+ .join(seq_for_transcriptome_build))
+
if (use_ref_ann){
results = preprocess_reads.out.report
- .concat(merge_gff_bundles.out.gff,
- run_gff_compare.out.stats,
- run_gff_compare.out.merged_annotated,
- run_gffread.out.transcriptome,
- makeReport.out.report
+ .concat(run_gffcompare.output.gffcmp_dir,
+ m.stats,
+ get_transcriptome.out.flatMap(map_sample_ids_cls),
)
- }else{ // Write a minimal report if no reference annotation is given
+ }
+ if (!use_ref_ann && !params.denovo){
results = preprocess_reads.out.report
- .concat(merge_gff_bundles.out.gff,
- makeReport.out.report
- )
- }
+ .concat(m.stats,
+ get_transcriptome.out.flatMap(map_sample_ids_cls),
+ )
+ }
+ if (params.denovo){
+ results = m.cds
+ .concat(m.stats,
+ seq_for_transcriptome_build,
+ get_transcriptome.out.flatMap(map_sample_ids_cls),
+ merge_gff_bundles.out.gff,
+ m.opt_qual_ch.flatMap {it ->
+ l = []
+ for (x in it[1..-1]){
+ l.add(tuple(it[0], x))
+ }
+ return l
+ })
+ }
+
emit:
results
+ report
telemetry = workflow_params
}
@@ -465,6 +465,9 @@ workflow {
start_ping()
params.results_dir = "${params.out_dir}/output"
+ report_template = file("$projectDir/bin/report_template.html")
+ table_template = file("$projectDir/bin/table_template.html")
+
fastq = file(params.fastq, type: "file")
if (!fastq.exists()) {
@@ -472,15 +475,27 @@ workflow {
exit 1
}
+ if (!params.denovo && !params.ref_genome){
+ println("--ref_genome must be supplied unless doing de novo assembly (--denovo)")
+ exit 1
+ }
+
+
if (params.ref_genome){
ref_genome = file(params.ref_genome, type: "file")
if (!ref_genome.exists()) {
- println("--reference: File doesn't exist, check path.")
+ println("--ref_genome: File doesn't exist, check path.")
exit 1
}
+ }else {
+ ref_genome = file("$projectDir/data/OPTIONAL_FILE")
+ }
+
+ if (params.denovo && params.ref_annotation) {
+ println("Reference annotation with de denovo assembly is not supported")
+ exit 1
}
- ref_annotation = null
if (params.ref_annotation){
ref_annotation = file(params.ref_annotation, type: "file")
if (!ref_annotation.exists()) {
@@ -495,11 +510,10 @@ workflow {
params.fastq, params.out_dir, params.sample, params.sample_sheet, params.sanitize_fastq
)
- pipeline(reads, ref_genome, ref_annotation)
+ pipeline(reads, ref_genome, ref_annotation, report_template, table_template)
- output(
- pipeline.out.results
- )
+ output(pipeline.out.results)
+ output_report(pipeline.out.report)
end_ping(pipeline.out.telemetry)
}
diff --git a/nextflow.config b/nextflow.config
index 1754941..afa7082 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -13,17 +13,19 @@
params {
help = false
fastq = null
- ref_genome = null
+ ref_genome = false
ref_annotation = null
// Process cDNA reads using pychopper, turn off for direct RNA:
use_pychopper = true
threads = 4
+ // Thresholds for viewing isoforms in report table
+ transcript_table_cov_thresh = 50
out_dir = null
sample = null
sample_sheet = null
sanitize_fastq = false
- wfversion = "v0.0.1"
+ wfversion = "v0.1.0"
aws_image_prefix = null
aws_queue = null
report_name = "report"
@@ -35,15 +37,16 @@ params {
// Options passed to pychopper:
- pychopper_opts = ""
+ pychopper_opts = "-m edlib"
// Extra option passed to minimap2 when generating index
minimap_index_opts = "-k14"
// Extra options passed to minimap2
+ // For SIRV data
+ //minimap2_opts = "-uf --splice-flank=no"
+ // AFor non-SIRV data:
minimap2_opts = "-uf"
- // Add this for SIRV data:
- // "--splice-flank=no"
// Minmum mapping quality
minimum_mapping_quality = 40
@@ -67,6 +70,56 @@ params {
plot_gffcmp_stats = true
disable_ping = false
+
+ //// Denovo-specific parameters
+ denovo = false
+
+ // Batch size in kilobases (if -1 then it is calculated based on the number of cores and bases):
+ batch_size = -1
+
+ // Maximum sequences per input batch (-1 means no limit):
+ batch_max_seq = -1
+
+ // Clustering mode:
+ cls_mode = "sahlin"
+
+ // Kmer size:
+ kmer_size = 11
+
+ // Window size:
+ window_size = 15
+
+ // Minimum cluser size in the left batch:
+ min_left_cls = 2
+
+ // Consensus period (-1 means no consensus):
+ consensus_period = 500
+
+ // Minimum consensus sample size:
+ consensus_minimum = 50
+
+ // Maximum consensus sample size:
+ consensus_maximum = -150
+
+ // Minimum number of minimizers shared between read and cluster:
+ min_shared = 5
+
+ // Minimum average quality value:
+ min_qual = 7.0
+
+ // Minmum mapped fraction of read to be included in cluster:
+ mapped_threshold = 0.65
+
+ // Minimum aligned fraction of read to be included in cluster:
+ aligned_threshold = 0.2
+
+ // Minimum fraction of minimizers shared compared to best hit, in order to continue mapping:
+ min_fraction = 0.8
+
+ // Minimum probability for i consecutive minimizers to be different between read and representative:
+ min_prob_no_hits = 0.1
+
+
}
manifest {
@@ -76,7 +129,7 @@ manifest {
description = 'RNA/cDNA isoform analysis workflow'
mainScript = 'main.nf'
nextflowVersion = '>=20.10.0'
- //version = 'v0.0.7' // TODO: do switch to this?
+ //version = 'v0.0.1'
}
executor {
diff --git a/nextflow_schema.json b/nextflow_schema.json
index 499c94e..cbf0488 100644
--- a/nextflow_schema.json
+++ b/nextflow_schema.json
@@ -63,7 +63,8 @@
},
"pychopper_opts": {
"type": "string",
- "description": "Extra pychopper opts"
+ "description": "Extra pychopper opts",
+ "default": "-m edlib"
},
"threads": {
"type": "integer",
@@ -110,11 +111,95 @@
},
"disable_ping": {
"type": "boolean"
+ },
+ "transcript_table_cov_thresh": {
+ "type": "integer",
+ "description": "Minimum coverage for a transcript to appear in the report table",
+ "default": 50
+ },
+ "denovo": {
+ "type": "boolean",
+ "description": "Use denovo transcript assembly rather than reference guided",
+ "default": false
+ },
+ "batch_size": {
+ "type": "integer",
+ "description": "Maximum sequences per input batch (-1 means no limit)",
+ "default": -1
+ },
+ "batch_max_seq": {
+ "type": "integer",
+ "description": "Maximum sequences per input batch (-1 means no limit)",
+ "default": -1
+ },
+ "cls_mode": {
+ "type": "string",
+ "description": "Clustering mode",
+ "default": "sahlin"
+ },
+ "kmer_size": {
+ "type": "integer",
+ "description": "Kmer size",
+ "default": 11
+ },
+ "window_size": {
+ "type": "integer",
+ "description": "Window size",
+ "default": 15
+ },
+ "min_left_cls": {
+ "type": "integer",
+ "description": "Minimum cluser size in the left batch",
+ "default": 2
+ },
+ "consensus_period": {
+ "type": "integer",
+ "description": "Consensus period (-1 means no consensus)",
+ "default": 500
+ },
+ "consensus_minimum": {
+ "type": "integer",
+ "description": "Minimum consensus sample size:",
+ "default": 50
+ },
+ "consensus_maximum": {
+ "type": "integer",
+ "description": "Maximum consensus sample size",
+ "default": -150
+ },
+ "min_shared": {
+ "type": "integer",
+ "description": "Minimum number of minimizers shared between read and cluster",
+ "default": 5
+ },
+ "min_qual": {
+ "description": "Minimum average quality value",
+ "type": "number",
+ "default": 7.0
+ },
+ "mapped_threshold": {
+ "description": "Minimum mapped fraction of read to be included in cluster",
+ "type": "number",
+ "default": 0.65
+ },
+ "aligned_threshold": {
+ "tpye": "number",
+ "description": "Minimum aligned fraction of read to be included in cluster",
+ "default": 0.2
+ },
+ "min_fraction": {
+ "type": "number",
+ "description": "Minimum fraction of minimizers shared compared to best hit, in order to continue mapping",
+ "default": 0.8
+ },
+ "min_prob_no_hits" : {
+ "type": "number",
+ "description": "Minimum probability for i consecutive minimizers to be different between read and representative",
+ "default": 0.2
}
},
"required": [
- "fastq",
- "ref_genome"
+ "fastq"
]
},
"meta_data": {
diff --git a/reference_assembly.nf b/reference_assembly.nf
new file mode 100644
index 0000000..e0da46a
--- /dev/null
+++ b/reference_assembly.nf
@@ -0,0 +1,48 @@
+process map_reads{
+ /*
+ Map reads to reference using minimap2.
+ Filter reads by mapping quality.
+ Filter reads where length of poly(A) > max_poly_run at either ends of the read (defined by poly_context)
+ */
+ label "isoforms"
+ cpus params.threads
+
+ input:
+ path index
+ path reference
+ tuple val(sample_id), path (fastq_reads)
+
+ output:
+ tuple val(sample_id), path("${sample_id}_reads_aln_sorted.bam"), emit: bam
+ tuple val(sample_id), path("${sample_id}_read_aln_stats.tsv"), emit: stats
+ script:
+ def ContextFilter = """AlnContext: { Ref: "${reference}", LeftShift: -${params.poly_context}, RightShift: ${params.poly_context},
+ RegexEnd: "[Aa]{${params.max_poly_run},}",
+ Stranded: True,Invert: True, Tsv: "internal_priming_fail.tsv"} """
+ """
+ seqkit fq2fa ${fastq_reads} -o "reads.fa";
+ minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} ${index} "reads.fa"\
+ | samtools view -q ${params.minimum_mapping_quality} -F 2304 -Sb -\
+ | seqkit bam -j ${params.threads} -x -T '${ContextFilter}' -\
+ | samtools sort -@ ${params.threads} -o "${sample_id}_reads_aln_sorted.bam" - \
+ | ((seqkit bam -s -j ${params.threads} - 2>&1) | tee ${sample_id}_read_aln_stats.tsv ) || true
+
+ if [[ -s "internal_priming_fail.tsv" ]];
+ then
+ tail -n +2 "internal_priming_fail.tsv" | awk '{{print ">" \$1 "\\n" \$4 }}' - > "context_internal_priming_fail_start.fasta"
+ tail -n +2 "internal_priming_fail.tsv" | awk '{{print ">" \$1 "\\n" \$6 }}' - > "context_internal_priming_fail_end.fasta"
+ fi
+ """
+}
+
+workflow reference_assembly {
+ take:
+ index
+ reference
+ fastq_reads
+ main:
+ map_reads(index, reference, fastq_reads)
+ emit:
+ bam = map_reads.out.bam
+ stats = map_reads.out.stats
+}
\ No newline at end of file
diff --git a/run_evaluation_dmel.sh b/run_evaluation_dmel.sh
deleted file mode 100755
index bc7e63b..0000000
--- a/run_evaluation_dmel.sh
+++ /dev/null
@@ -1,40 +0,0 @@
-#!/bin/bash
-
-# Usage: ./run_evaluation_dmel.sh pathto/outputdir
-
-# See the isONcorrect paper https://www.nature.com/articles/s41467-020-20340-8 where this dataset is described
-
-OUTDIR=$1;
-
-FASTQ_URL="http://ftp.sra.ebi.ac.uk/vol1/fastq/ERR358/005/ERR3588905/ERR3588905_1.fastq.gz"
-REF_URL="http://ftp.ensembl.org/pub/release-99/fasta/drosophila_melanogaster/dna/Drosophila_melanogaster.BDGP6.28.dna.toplevel.fa.gz"
-GFF_URL="http://ftp.ensembl.org/pub/release-99/gff3/drosophila_melanogaster/Drosophila_melanogaster.BDGP6.28.99.gff3.gz"
-
-RESULTS_DIR="$OUTDIR/results"
-DATA_DIR="$OUTDIR/data"
-READS_DIR="$DATA_DIR/reads"
-FASTQ="$READS_DIR/ERR3588905_1.fastq"
-REF="$DATA_DIR/Drosophila_melanogaster.BDGP6.28.dna.toplevel.fa"
-GFF="$DATA_DIR/Drosophila_melanogaster.BDGP6.28.99.gff3"
-
-echo $READS_DIR;
-
-rm -fr $OUT_DIR/results
-mkdir -p $OUTDIR/data
-
-if [ ! -f $REF ];
-then (cd $DATA_DIR; curl -L -C - -O $REF_URL); gzip -d ${REF}.gz
-fi
-
-if [ ! -f $GFF ]
-then
- (cd $DATA_DIR; curl -L -C - -O $GFF_URL); gzip -d ${GFF}.gz
-fi
-
-if [ ! -f $FASTQ ];
-then (cd $READS_DIR; curl -L -C - -O $FASTQ_URL); gzip -d ${FASTQ}.gz
-fi
-
-nextflow run ../wf-isoforms --fastq $READS_DIR \
---reference_genome $REF --annotation $GFF -profile conda --out_dir $OUTDIR \
--w $OUTDIR/workspace -resume
diff --git a/test_data/demultiplexed_fastq/barcode02/SIRV_E0_PCS109_25.fq.gz b/test_data/demultiplexed_fastq/barcode02/SIRV_E0_PCS109_25.fq.gz
new file mode 100644
index 0000000..6497611
Binary files /dev/null and b/test_data/demultiplexed_fastq/barcode02/SIRV_E0_PCS109_25.fq.gz differ
diff --git a/test_data/demultiplexed_fastq/barcode02/SIRV_E0_PCS109_50.fq.gz b/test_data/demultiplexed_fastq/barcode02/SIRV_E0_PCS109_50.fq.gz
deleted file mode 100644
index 54a22d3..0000000
Binary files a/test_data/demultiplexed_fastq/barcode02/SIRV_E0_PCS109_50.fq.gz and /dev/null differ
diff --git a/tests.sh b/tests.sh
deleted file mode 100755
index 109bdac..0000000
--- a/tests.sh
+++ /dev/null
@@ -1,53 +0,0 @@
-#!/usr/bin/env bash
-
-# A couple of simple tests with different combinations of CLI options
-
-if [[ "$#" -ne 1 ]]; then
- echo "Please supply path to out_dir"
- exit 1
-fi
-
-SCRIPT_DIR="$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )"
-cd $SCRIPT_DIR;
-
-singledir="test_data/fastq"
-multisampledir="test_data/demultiplexed_fastq"
-
-echo $singledir
-echo $multisampledir
-results=()
-
-# Test1 single sample
-OUTPUT=$1/test1;
-nextflow run . --fastq $singledir --ref_genome test_data/SIRV_150601a.fasta \
---ref_annotation test_data/SIRV_isofroms.gtf -profile conda --out_dir ${OUTPUT} -w ${OUTPUT}/workspace -resume;
-
-results+=("test1: $?")
-
-# Test2 multiple samples demultiplexed
-OUTPUT=$1/test2;
-nextflow run . --fastq $multisampledir --ref_genome test_data/SIRV_150601a.fasta \
---ref_annotation test_data/SIRV_isofroms.gtf -profile conda --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
---sample_sheet test_data/sample_sheet -resume;
-
-results+=("test2: $?")
-
-# Test3 single sample. No reference annotation
-OUTPUT=$1/test3;
-nextflow run . --fastq $singledir --ref_genome test_data/SIRV_150601a.fasta \
- -profile conda --out_dir ${OUTPUT} -w ${OUTPUT}/workspace -resume;
-
-results+=("test3: $?")
-
-# Test4 single sample. force split_bam to make multiple alignment bundles
-OUTPUT=$1/test4;
-nextflow run . --fastq $singledir --ref_genome test_data/SIRV_150601a.fasta \
---ref_annotation test_data/SIRV_isofroms.gtf -profile conda --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
---bundle_min_reads 5 -resume;
-
-results+=("test4: $?")
-
-echo "Exit status codes for each test"
-for value in "${results[@]}"; do
- echo "${value}"
-done
\ No newline at end of file