From f683d4441103d53824db7e190e84e2261fc86425 Mon Sep 17 00:00:00 2001 From: Sarah Griffiths Date: Wed, 27 Sep 2023 08:02:00 +0000 Subject: [PATCH] CW-2769 remove denovo --- .gitlab-ci.yml | 6 +- CHANGELOG.md | 3 + README.md | 30 +-- bin/workflow_glue/report.py | 27 +-- bin/workflow_glue/run_isonclust2.py | 138 -------------- docs/header.md | 3 +- docs/intro.md | 10 - docs/links.md | 5 +- docs/quickstart.md | 13 +- evaluation/tests.sh | 12 -- main.nf | 62 ++---- nextflow.config | 9 +- nextflow_schema.json | 31 +-- subworkflows/denovo_assembly.nf | 285 ---------------------------- 14 files changed, 36 insertions(+), 598 deletions(-) delete mode 100755 bin/workflow_glue/run_isonclust2.py delete mode 100644 subworkflows/denovo_assembly.nf diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 13c8148..d5a391b 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -45,7 +45,7 @@ docker-run: - MATRIX_NAME: [ "fusions", "differential_expression", "isoforms", "only_differential_expression", "differential_expression_gff3", - "ncbi_gzip", "denovo", "ncbi_no_gene_id", "ensembl_with_versions", + "ncbi_gzip", "ncbi_no_gene_id", "ensembl_with_versions", "differential_expression_mouse" ] rules: @@ -60,10 +60,6 @@ docker-run: NF_WORKFLOW_OPTS: "--fastq ERR6053095_chr20.fastq --transcriptome-source reference-guided \ --ref_genome chr20/hg38_chr20.fa --ref_annotation chr20/gencode.v22.annotation.chr20.gtf" NF_IGNORE_PROCESSES: preprocess_reads,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome - - if: $MATRIX_NAME == "denovo" - variables: - NF_WORKFLOW_OPTS: "--fastq test_data/fastq/SIRV_E0_PCS109_50.fq.gz --transcriptome_source denovo" - NF_IGNORE_PROCESSES: preprocess_reads,merge_transcriptomes,decompress_annotation,decompress_ref,build_minimap_index,decompress_transcriptome,preprocess_ref_transcriptome - if: $MATRIX_NAME == "fusions" variables: NF_BEFORE_SCRIPT: wget -O test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf test_data.tar.gz diff --git a/CHANGELOG.md b/CHANGELOG.md index 627641c..33d3174 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,7 +5,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [unreleased] +### Fixed - Remove dead links from README +### Removed +- Denovo `--transcriptome_source` option.` ## [v0.3.1] ### Added diff --git a/README.md b/README.md index 904d267..35284c9 100644 --- a/README.md +++ b/README.md @@ -5,8 +5,6 @@ for assembly and annotation of transcripts from Oxford Nanopore cDNA or direct R - - ## Introduction This workflow identifies RNA isoforms using either cDNA or direct RNA (dRNA) @@ -27,16 +25,6 @@ in long read mode (with or without a guide reference annotation) to generate the * The annotation generated by the pipeline is compared to the reference annotation. using [gffcompare](http://ccb.jhu.edu/software/stringtie/gffcompare.shtml) -#### de novo-based transcript assembly (experimental!) -* Sequence clusters are generated using [isONclust2](https://github.com/nanoporetech/isONclust2) - * If a reference genome is supplied, cluster quality metrics are determined by comparing - with clusters generated from a minimap2 alignment. -* A consensus sequence for each cluster is generated using [spoa](https://github.com/rvaser/spoa) -* Three rounds of polishing using racon and minimap2 to give a final polished CDS for each gene. -* Full-length reads are then mapped to these polished CDS. -* Transcripts are assembled by stringtie as for the reference-based approach. -* __Note__: This approach is currently not supported with direct RNA reads. - ### Fusion gene detection Fusion gene detection is performed using [JAFFA](https://github.com/Oshlack/JAFFA), with the JAFFAL extension for use with ONT long reads. @@ -134,16 +122,6 @@ nextflow run epi2me-labs/wf-transcriptomes \ --out_dir outdir -w workspace_dir ``` -**Example workflow for denovo transcript assembly** -``` -OUTPUT=~/output -nextflow run epi2me-labs/wf-transcriptomes \ - --fastq test_data/fastq \ - --transcriptome_source denovo \ - --out_dir ${OUTPUT} \ - -w ${OUTPUT}/workspace \ - --sample sample_id -``` A full list of options can be seen in nextflow_schema.json. Parameters can be specified either in a config like `parameter = value` or on the command line like `--parameter value`. Below are some commonly used parameters in the format used in config files. @@ -151,8 +129,7 @@ Below are some commonly used parameters in the format used in config files. Select how the transcriptome used for analysis should be prepared: - To create a reference transcriptome using an existing reference genome `--transcriptome_source reference-guided` (default) -- Use a a supplied transcriptome `--transcriptome_source precomputed"` -- Gnerate transcriptome via the denovo pipeline `--transcriptome_source denovo"` +- Use a supplied transcriptome `--transcriptome_source precomputed"` To run the workflow with direct RNA reads `--direct_rna true` (this just skips the pychopper step). @@ -297,7 +274,4 @@ in `${out_dir}/jaffal_output_${sample_id}` you will find: * [nextflow](https://www.nextflow.io/) * [docker](https://www.docker.com/products/docker-desktop) * [Singularity](https://sylabs.io/singularity/) -* [racon](https://github.com/isovic/racon) -* [spoa](https://github.com/rvaser/spoa) -* [inONclust](https://github.com/ksahlin/isONclust) -* [isONclust2](https://github.com/nanoporetech/isONclust2) \ No newline at end of file +* [racon](https://github.com/isovic/racon) \ No newline at end of file diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index c5a45c3..5cdc433 100755 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -74,7 +74,6 @@ def argparser(): parser.add_argument( "--de_stats", required=False, type=str, default=None, nargs='*', help="Differential expression report optional") - parser.add_argument('--denovo', dest='denovo', action='store_true') return parser @@ -699,21 +698,16 @@ def transcript_table(report, df_tmaps, max_rows): section.table(df, index=False) -def transcriptome_summary(report, gffs, sample_ids, denovo=False): +def transcriptome_summary(report, gffs, sample_ids): """ Plot transcriptome summaries. Some of this data is available via gffcompare output, but the de novo pipeline skips that, so we do it all here. - We do not report exon number for the denovo assembly yet. This is because - in this case, the gff annotation is generated by aligning to the CDS not - the genome. - :param report: aplanat WFReport :param gffs: list of paths to gff transcriptome annotations :param sample_ids: list of sample ids - :param denovo: whether annotation was generated by de novo pipeline or not """ # test.db gets written to the git repo. section = report.add_section() @@ -771,17 +765,16 @@ def transcriptome_summary(report, gffs, sample_ids, denovo=False): title='transcript lengths') plots.append(box) - if not denovo: - x, y = zip(*sorted(exons_per_transcript.items())) + x, y = zip(*sorted(exons_per_transcript.items())) - fig = figure(title="Exons per transcript") - fig.vbar( - x, top=list(y), color=Colors.cerulean) - fig.xaxis.axis_label = 'Num. exons' - fig.yaxis.axis_label = 'Num. genes' + fig = figure(title="Exons per transcript") + fig.vbar( + x, top=list(y), color=Colors.cerulean) + fig.xaxis.axis_label = 'Num. exons' + fig.yaxis.axis_label = 'Num. genes' - fig.xaxis.major_label_orientation = math.pi / 2.8 - plots.append(fig) + fig.xaxis.major_label_orientation = math.pi / 2.8 + plots.append(fig) df_sum = pd.DataFrame.from_dict( {'Total genes': [num_genes], @@ -929,7 +922,7 @@ def main(args): # Results if args.gff_annotation is not None: transcriptome_summary( - report, args.gff_annotation, sample_ids, denovo=args.denovo) + report, args.gff_annotation, sample_ids) if args.gffcompare_dir is not None: df_tmaps = gff_compare_plots( diff --git a/bin/workflow_glue/run_isonclust2.py b/bin/workflow_glue/run_isonclust2.py deleted file mode 100755 index b9541f5..0000000 --- a/bin/workflow_glue/run_isonclust2.py +++ /dev/null @@ -1,138 +0,0 @@ -#!/usr/bin/env python -"""Dynamically generate isONclust2 processes.""" -from collections import OrderedDict -from glob import glob -from itertools import zip_longest -import os -from pathlib import Path -import re -import subprocess as sub - -from .util import wf_parser # noqa: ABS101 - - -def argparser(): - """Argument parser for entrypoint.""" - parser = wf_parser("report") - parser.add_argument( - "--workdir", help="directory containing batches/ dir [CWD]", - default=Path()) - - return parser - - -class Node: - """Node.""" - - def __init__(self, node_id, file_, left, right, parent, level): - """Set node attaributes.""" - self.Id = node_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(args): - """Entry point.""" - os.chdir(args.workdir) - 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) diff --git a/docs/header.md b/docs/header.md index 03d4f25..e4bbb42 100644 --- a/docs/header.md +++ b/docs/header.md @@ -1,5 +1,4 @@ # wf-transcriptomes This repository contains a [nextflow](https://www.nextflow.io/) workflow -for assembly and annotation of transcripts from Oxford Nanopore cDNA or direct RNA reads. - +for assembly and annotation of transcripts from Oxford Nanopore cDNA or direct RNA reads. \ No newline at end of file diff --git a/docs/intro.md b/docs/intro.md index 790046d..24e6152 100644 --- a/docs/intro.md +++ b/docs/intro.md @@ -18,16 +18,6 @@ in long read mode (with or without a guide reference annotation) to generate the * The annotation generated by the pipeline is compared to the reference annotation. using [gffcompare](http://ccb.jhu.edu/software/stringtie/gffcompare.shtml) -#### de novo-based transcript assembly (experimental!) -* Sequence clusters are generated using [isONclust2](https://github.com/nanoporetech/isONclust2) - * If a reference genome is supplied, cluster quality metrics are determined by comparing - with clusters generated from a minimap2 alignment. -* A consensus sequence for each cluster is generated using [spoa](https://github.com/rvaser/spoa) -* Three rounds of polishing using racon and minimap2 to give a final polished CDS for each gene. -* Full-length reads are then mapped to these polished CDS. -* Transcripts are assembled by stringtie as for the reference-based approach. -* __Note__: This approach is currently not supported with direct RNA reads. - ### Fusion gene detection Fusion gene detection is performed using [JAFFA](https://github.com/Oshlack/JAFFA), with the JAFFAL extension for use with ONT long reads. diff --git a/docs/links.md b/docs/links.md index 9d2b74c..0a71673 100644 --- a/docs/links.md +++ b/docs/links.md @@ -3,7 +3,4 @@ * [nextflow](https://www.nextflow.io/) * [docker](https://www.docker.com/products/docker-desktop) * [Singularity](https://sylabs.io/singularity/) -* [racon](https://github.com/isovic/racon) -* [spoa](https://github.com/rvaser/spoa) -* [inONclust](https://github.com/ksahlin/isONclust) -* [isONclust2](https://github.com/nanoporetech/isONclust2) \ No newline at end of file +* [racon](https://github.com/isovic/racon) \ No newline at end of file diff --git a/docs/quickstart.md b/docs/quickstart.md index c3da56d..0a33df2 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -47,16 +47,6 @@ nextflow run epi2me-labs/wf-transcriptomes \ --out_dir outdir -w workspace_dir ``` -**Example workflow for denovo transcript assembly** -``` -OUTPUT=~/output -nextflow run epi2me-labs/wf-transcriptomes \ - --fastq test_data/fastq \ - --transcriptome_source denovo \ - --out_dir ${OUTPUT} \ - -w ${OUTPUT}/workspace \ - --sample sample_id -``` A full list of options can be seen in nextflow_schema.json. Parameters can be specified either in a config like `parameter = value` or on the command line like `--parameter value`. Below are some commonly used parameters in the format used in config files. @@ -64,8 +54,7 @@ Below are some commonly used parameters in the format used in config files. Select how the transcriptome used for analysis should be prepared: - To create a reference transcriptome using an existing reference genome `--transcriptome_source reference-guided` (default) -- Use a a supplied transcriptome `--transcriptome_source precomputed"` -- Gnerate transcriptome via the denovo pipeline `--transcriptome_source denovo"` +- Use a supplied transcriptome `--transcriptome_source precomputed"` To run the workflow with direct RNA reads `--direct_rna true` (this just skips the pychopper step). diff --git a/evaluation/tests.sh b/evaluation/tests.sh index 59daaff..f6b95a7 100755 --- a/evaluation/tests.sh +++ b/evaluation/tests.sh @@ -26,18 +26,6 @@ multisampledir="test_data/demultiplexed_fastq" #"--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' \ diff --git a/main.nf b/main.nf index c7fc954..496abee 100644 --- a/main.nf +++ b/main.nf @@ -2,7 +2,6 @@ /* 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; @@ -12,7 +11,6 @@ nextflow.enable.dsl = 2 include { fastq_ingress } from './lib/fastqingress' include { reference_assembly } from './subworkflows/reference_assembly' -include { denovo_assembly } from './subworkflows/denovo_assembly' include { gene_fusions } from './subworkflows/JAFFAL/gene_fusions' include { differential_expression } from './subworkflows/differential_expression' @@ -38,8 +36,6 @@ process getVersions { seqkit 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 """ } @@ -285,13 +281,7 @@ process run_gffcompare{ path ("${sample_id}_annotated.gtf"), emit: gtf, optional: true script: def out_dir = "${sample_id}_gffcompare" - - if (params.transcriptome_source == "denovo"){ - """ - mkdir $out_dir - """ - } else { - """ + """ mkdir $out_dir echo "Doing comparison of reference annotation: ${ref_annotation} and the query annotation" @@ -305,10 +295,10 @@ process run_gffcompare{ mv *.refmap $out_dir cp ${out_dir}/str_merged.annotated.gtf ${sample_id}_annotated.gtf """ - } } + process get_transcriptome{ /* Write out a transcriptome file based on the query gff annotations. @@ -380,7 +370,6 @@ process makeReport { // Convert the sample_id arrayList. sids = new BlankSeparatedList(sample_ids) def report_name = "wf-transcriptomes-report.html" - def OPT_DENOVO = params.transcriptome_source == "denovo" ? "--denovo" : '' """ if [ -f "de_report/OPTIONAL_FILE" ]; then dereport="" @@ -419,7 +408,6 @@ process makeReport { \$OPT_GFF \ --isoform_table_nrows $params.isoform_table_nrows \ \$OPT_JAFFAL_CSV \ - $OPT_DENOVO \ \$dereport """ } @@ -550,17 +538,10 @@ workflow pipeline { pychopper_report = file("$projectDir/data/OPTIONAL_FILE") } if (params.transcriptome_source != "precomputed"){ + build_minimap_index(ref_genome) + log.info("Doing reference based transcript analysis") + assembly = reference_assembly(build_minimap_index.out.index, ref_genome, full_len_reads) - if (params.transcriptome_source == "denovo"){ - log.info("Doing de novo assembly") - log.info("WARNING: The `--transcriptome_source` denovo option may have unexpected results and errors. If possible it is preferable to use the reference-guided pipeline.") - assembly = denovo_assembly(full_len_reads, ref_genome) - - } else { - build_minimap_index(ref_genome) - log.info("Doing reference based transcript analysis") - assembly = reference_assembly(build_minimap_index.out.index, ref_genome, full_len_reads) - } assembly_stats = assembly.stats.map{ it -> it[1]}.collect() split_bam(assembly.bam) @@ -569,15 +550,10 @@ workflow pipeline { merge_gff_bundles(assemble_transcripts.out.gff_bundles.groupTuple()) run_gffcompare(merge_gff_bundles.out.gff, ref_annotation) - - if (params.transcriptome_source == "denovo"){ - // Use the per-sample, de novo-assembled CDS - seq_for_transcriptome_build = assembly.cds - }else { - // For reference based assembly, there is only one reference - // So map this reference to all sample_ids - seq_for_transcriptome_build = sample_ids.flatten().combine(ref_genome) - } + // For reference based assembly, there is only one reference + // So map this reference to all sample_ids + seq_for_transcriptome_build = sample_ids.flatten().combine(ref_genome) + get_transcriptome( merge_gff_bundles.out.gff @@ -664,22 +640,6 @@ workflow pipeline { .concat(results) } - if (params.transcriptome_source == "denovo"){ - results = assembly.cds.concat( - assembly.stats, - seq_for_transcriptome_build, - get_transcriptome.out.transcriptome.flatMap(map_sample_ids_cls), - assembly.opt_qual_ch.flatMap { - it -> - l = [] - for (x in it[1..-1]){ - l.add(tuple(it[0], x)) - } - return l - }) - .map {it -> it[1]} - .concat(results) - } if (params.jaffal_refBase){ results = results .concat(gene_fusions.out.results @@ -734,8 +694,8 @@ workflow { }else { ref_genome = file("$projectDir/data/OPTIONAL_FILE") } - if (params.transcriptome_source == "denovo" && params.ref_annotation) { - error = "Reference annotation with de denovo assembly is not supported" + if (params.containsValue("denovo")) { + error = "Denovo transcriptome source is no longer supported. Please use the reference-guided or precomputed options." } if (params.ref_annotation){ diff --git a/nextflow.config b/nextflow.config index 773f5c8..e3ba58b 100644 --- a/nextflow.config +++ b/nextflow.config @@ -71,11 +71,6 @@ params { disable_ping = false - //// Denovo-specific parameters - - isOnClust2_batch_size = -1 - isOnClust2_sort_options = "--batch-size -1 --kmer-size 11 --window-size 15 --min-shared 5 --min-qual 7.0 --mapped-threshold 0.65 --aligned-threshold 0.2 --min-fraction 0.8 --min-prob-no-hits 0.0 -M -1 -P 500 -g 50 -c 150 -F 2" - ////// Fusion detection parameters jaffal_refBase = null jaffal_genome = "hg38" @@ -103,7 +98,7 @@ params { "--ref_genome 'wf-transcriptomes-demo/hg38_chr20.fa'", ] agent = null - container_sha = "sha203915eb4b4dd444cb2e845d0b9f7814e26b7b5c" + container_sha = "shab7142eb91ba48033b690c29d2a4252c45d85a1de" common_sha = "sha0a6dc21fac17291f4acb2e0f67bcdec7bf63e6b7" } } @@ -126,7 +121,7 @@ executor { } epi2melabs { - tags = "isoforms, transcriptomics, denovo" + tags = "isoforms, transcriptomics" } // used by default for "standard" (docker) and singularity profiles, diff --git a/nextflow_schema.json b/nextflow_schema.json index 25de66a..1ceb9b3 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -26,12 +26,11 @@ "type": "string", "enum": [ "precomputed", - "reference-guided", - "denovo" + "reference-guided" ], "default": "reference-guided", "description": "Select how the transcriptome used for analysis should be prepared.", - "help_text": "To analyse only gene fusions and differential expression use of an existing transcriptome may be preferred and so 'precomputed' should be selected. In this case the 'ref_transcriptome' parameter should be specified. To create a reference transcriptome using an existing reference genome, select 'reference guided' and specify the 'ref_genome' parameter. To create a transcriptome from your sequencing data select 'denovo'." + "help_text": "To analyse only gene fusions and differential expression use of an existing transcriptome may be preferred and so 'precomputed' should be selected. In this case the 'ref_transcriptome' parameter should be specified. To create a reference transcriptome using an existing reference genome, select 'reference guided' and specify the 'ref_genome' parameter." }, "ref_genome": { "type": "string", @@ -172,25 +171,6 @@ } } }, - "denovo_wf_options": { - "title": "Options for de novo-based workflow", - "type": "object", - "description": "Parameters that are used solely for the de novo workflow", - "properties": { - "isOnClust2_batch_size": { - "type": "integer", - "description": "Number of batches to to process the data in.", - "help_text": "If set to -1 number of batches will be the same as the number of threads avaiable.", - "default": -1 - }, - "isOnClust2_sort_options": { - "type": "string", - "default": "--batch-size -1 --kmer-size 11 --window-size 15 --min-shared 5 --min-qual 7.0 --mapped-threshold 0.65 --aligned-threshold 0.2 --min-fraction 0.8 --min-prob-no-hits 0.0 -M -1 -P 500 -g 50 -c 150 -F 2", - "description": "Additional command-line options for isOnClust2 sort.", - "help_text": "isOnClust2 is used for **de novo** transcript assembly. Options for the sort command can be be supplied like so `-opt1 arg -opt2 arg`. Available arguments can be found at [isOnClust2](https://github.com/nanoporetech/isONclust2). It is recommended not to alter this parameter." - } - } - }, "gene_fusion_detection_options": { "title": "Gene Fusion Detection Options", "type": "object", @@ -341,9 +321,6 @@ { "$ref": "#/definitions/options_for_reference_based_workflow" }, - { - "$ref": "#/definitions/denovo_wf_options" - }, { "$ref": "#/definitions/gene_fusion_detection_options" }, @@ -378,7 +355,7 @@ } }, "docs": { - "intro": "## Introduction\n\nThis workflow identifies RNA isoforms using either cDNA or direct RNA (dRNA) \nOxford Nanopore reads.\n\n### Preprocesing\ncDNA reads are initially preprocessed by [pychopper](https://github.com/epi2me-labs/pychopper) \nfor the identification of full-length reads, as well as trimming and orientation correction (This step is omitted for \n direct RNA reads).\n\n\n### Transcript assembly\n\n#### Reference-aided transcript assembly approach\n* Full length reads are mapped to a supplied reference genome using [minimap2](https://github.com/lh3/minimap2)\n* Transcripts are assembled by [stringtie](http://ccb.jhu.edu/software/stringtie) \nin long read mode (with or without a guide reference annotation) to generate the GFF annotation.\n* The annotation generated by the pipeline is compared to the reference annotation. \nusing [gffcompare](http://ccb.jhu.edu/software/stringtie/gffcompare.shtml)\n\n#### de novo-based transcript assembly (experimental!)\n* Sequence clusters are generated using [isONclust2](https://github.com/nanoporetech/isONclust2)\n * If a reference genome is supplied, cluster quality metrics are determined by comparing \n with clusters generated from a minimap2 alignment.\n* A consensus sequence for each cluster is generated using [spoa](https://github.com/rvaser/spoa)\n* Three rounds of polishing using racon and minimap2 to give a final polished CDS for each gene.\n* Full-length reads are then mapped to these polished CDS.\n* Transcripts are assembled by stringtie as for the reference-based approach.\n* __Note__: This approach is currently not supported with direct RNA reads.\n\n### Fusion gene detection\nFusion gene detection is performed using [JAFFA](https://github.com/Oshlack/JAFFA), with the JAFFAL extension for use \nwith ONT long reads. \n\n### Differential expression analysis\n\nDifferential gene expression (DGE) and differential transcript usage (DTU) analyses aim to identify genes and/or transcripts that show statistically altered expression patterns in a studied biological system. The results of the differential analyses are presented in a quantitative format and therefore the degree of change (up or down regulation) between experimental conditions can be calculated for each gene identified.\n\nThese differential analyses work by taking a \u201csnapshot\u201d of mRNA abundance and calculating the relative levels of transcripts and isoforms. In this context, expression corresponds to the number of messenger RNAs (mRNA) measured from each gene isoform within the organism / tissue / culture being investigated. In order to determine expression levels across the whole genome, sequence data specifically targeting the mRNA molecules can be generated.\n\nOxford Nanopore Technologies provides a number of sequencing solutions to allow users to generate the required snapshot of gene expression. This can be achieved by both sequencing the mRNA directly, or via a complementary DNA (cDNA) proxy. In contrast to short read sequencing technologies, entire mRNA transcripts can be captured as single reads. The example data provided with this tutorial is from a study based on the PCR-cDNA kit. This is a robust choice for performing differential transcript usage studies. This kit is suitable for preparation of sequence libraries from low mRNA input quantities. The cDNA population is enriched through PCR with low bias; an important prerequisite for the subsequent statistical analysis.\n\n[Workflow-transcriptomes](https://github.com/epi2me-labs/wf-transcriptomes) includes a subworkflow for DGE and DTU. The first step involves using either a reference alignment or _de novo_ assembly approach to create a set of mRNA sequences per sample. These are merged into a non-redundant transcriptome using [stringtie merge](http://ccb.jhu.edu/software/stringtie). The reads are then aligned to the transcriptome using minimap2 in a splice-aware manner. [Salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification, giving per transcript counts and then the following R packages are used for analysis.\n\n### Pre-filtering of quantitative data using DRIMSeq\nDRIMSeq (Nowicka and Robinson (2016)) is used to filter the transcript count data from the salmon analysis. The filter step will be used to select for genes and transcripts that satisfy rules for the number of samples in which a gene or transcript must be observed and minimum threshold levels for the number of observed reads. The parameters used for filtering are defined in the config.yaml file. The default parameters defined for this analysis include\n* min_samps_gene_expr = 3 - a transcript must be mapped to a gene in at least this minimum number of samples for the gene be included in the analysis\n*\tmin_samps_feature_expr = 1 - a transcript must be mapped to an isoform in at least this this minimum number of samples for the gene isoform to be included in the analysis\n*\tmin_gene_expr = 10 - the minimum number of total mapped sequence reads for a gene to be considered expressed\n*\tmin_feature_expr = 3 - the minimum number of total mapped sequence reads for a gene isoform to be considered\n\n### edgeR based differential expression analysis\n+A statistical analysis is first performed using edgeR (Robinson, McCarthy, and Smyth (2010), McCarthy et al. (2012)) to identify the subset of differentially expressed genes. The filtered list of gene counts is used as input. A normalisation factor is calculated for each sequence library (using the default TMM method - please see McCarthy et al. (2012) for further details). The defined experimental design is used to calculate estimates of dispersion for each of the gene features. Statistical tests are calculated using the contrasts defined in the experimental design. The differentially expressed genes are corrected for false discovery (fdr) using the method of Benjamini & Hochberg (Benjamini and Hochberg (1995))\n\n### Differential transcript usage using DEXSeq\nDifferential transcript usage analysis is performed using the R DEXSeq package (Reyes et al. (2013)). Similar to the edgeR package, DEXSeq estimates the variance between the biological replicates and applies generalised linear models for the statistical testing. The key difference is that the DEXSeq method looks for differences at the exon count level. DEXSeq uses the filtered transcript count data prepared earlier in this analysis. \n\n### StageR stage-wise analysis of DGE and DTU\nThe final component of this isoform analysis is a stage-wise statistical test using the R software package `stageR` (Van den Berge and Clement (2018)). stageR uses (1) the raw p-values for DTU from the DEXSeq analysis in the previous section and (2) a false-discovery corrected set of p-values from testing whether individual genes contain at least one exon showing DTU. A hierarchical two-stage statistical testing evaluates the set of genes for DTU.\n\n## Running the workflow\nFor the differential expression analysis section you should have at least 3 repeats for each sample. \nYour FASTQ data will need to be organised in to 6 directories that represent 3 repeats for each condition. \n\n\n## Analysis \nDifferential gene expression is sensitive to the input data quantity and quality. There should be equivalence between samples in the number of sequence reads, mapped reads and quality scores. The sequence and alignment summary plots in the report can be used to assess these metrics. There is also a table that shows the transcript per million(TPM) calculated from the salmon counts. TPM normalizes the data for gene length and then sequencing depth, and makes it easier to compare across samples compared to counts.\n\n### Workflow inputs\n- Directory containing cDNA/direct RNA reads. Or a directory containing subdirectories each with reads from different samples\n (in fastq/fastq.gz format)\n- Reference genome in fasta format (required for reference-based assembly).\n- Optional reference annotation in GFF2/3 format (extensions allowed are .gtf(.gz), .gff(.gz), .gff3(.gz)) (required for differential expression analysis `--de_analysis`). Only annotation files from [Encode](https://www.encodeproject.org), [Ensembl](https://www.ensembl.org/index.html) and [NCBI](https://www.ncbi.nlm.nih.gov/) are supported.\n- For fusion detection, JAFFAL reference files (see Quickstart) \n", - "links": "## Useful links\n\n* [nextflow](https://www.nextflow.io/)\n* [docker](https://www.docker.com/products/docker-desktop)\n* [Singularity](https://sylabs.io/singularity/)\n* [racon](https://github.com/isovic/racon)\n* [spoa](https://github.com/rvaser/spoa)\n* [inONclust](https://github.com/ksahlin/isONclust)\n* [isONclust2](https://github.com/nanoporetech/isONclust2)" + "intro": "## Introduction\n\nThis workflow identifies RNA isoforms using either cDNA or direct RNA (dRNA) \nOxford Nanopore reads.\n\n### Preprocesing\ncDNA reads are initially preprocessed by [pychopper](https://github.com/epi2me-labs/pychopper) \nfor the identification of full-length reads, as well as trimming and orientation correction (This step is omitted for \n direct RNA reads).\n\n\n### Transcript assembly\n\n#### Reference-aided transcript assembly approach\n* Full length reads are mapped to a supplied reference genome using [minimap2](https://github.com/lh3/minimap2)\n* Transcripts are assembled by [stringtie](http://ccb.jhu.edu/software/stringtie) \nin long read mode (with or without a guide reference annotation) to generate the GFF annotation.\n* The annotation generated by the pipeline is compared to the reference annotation. \nusing [gffcompare](http://ccb.jhu.edu/software/stringtie/gffcompare.shtml)\n\n### Fusion gene detection\nFusion gene detection is performed using [JAFFA](https://github.com/Oshlack/JAFFA), with the JAFFAL extension for use \nwith ONT long reads. \n\n### Differential expression analysis\n\nDifferential gene expression (DGE) and differential transcript usage (DTU) analyses aim to identify genes and/or transcripts that show statistically altered expression patterns in a studied biological system. The results of the differential analyses are presented in a quantitative format and therefore the degree of change (up or down regulation) between experimental conditions can be calculated for each gene identified.\n\nThese differential analyses work by taking a \u201csnapshot\u201d of mRNA abundance and calculating the relative levels of transcripts and isoforms. In this context, expression corresponds to the number of messenger RNAs (mRNA) measured from each gene isoform within the organism / tissue / culture being investigated. In order to determine expression levels across the whole genome, sequence data specifically targeting the mRNA molecules can be generated.\n\nOxford Nanopore Technologies provides a number of sequencing solutions to allow users to generate the required snapshot of gene expression. This can be achieved by both sequencing the mRNA directly, or via a complementary DNA (cDNA) proxy. In contrast to short read sequencing technologies, entire mRNA transcripts can be captured as single reads. The example data provided with this tutorial is from a study based on the PCR-cDNA kit. This is a robust choice for performing differential transcript usage studies. This kit is suitable for preparation of sequence libraries from low mRNA input quantities. The cDNA population is enriched through PCR with low bias; an important prerequisite for the subsequent statistical analysis.\n\n[Workflow-transcriptomes](https://github.com/epi2me-labs/wf-transcriptomes) includes a subworkflow for DGE and DTU. The first step involves using either a reference alignment or _de novo_ assembly approach to create a set of mRNA sequences per sample. These are merged into a non-redundant transcriptome using [stringtie merge](http://ccb.jhu.edu/software/stringtie). The reads are then aligned to the transcriptome using minimap2 in a splice-aware manner. [Salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification, giving per transcript counts and then the following R packages are used for analysis.\n\n### Pre-filtering of quantitative data using DRIMSeq\nDRIMSeq (Nowicka and Robinson (2016)) is used to filter the transcript count data from the salmon analysis. The filter step will be used to select for genes and transcripts that satisfy rules for the number of samples in which a gene or transcript must be observed and minimum threshold levels for the number of observed reads. The parameters used for filtering are defined in the config.yaml file. The default parameters defined for this analysis include\n* min_samps_gene_expr = 3 - a transcript must be mapped to a gene in at least this minimum number of samples for the gene be included in the analysis\n*\tmin_samps_feature_expr = 1 - a transcript must be mapped to an isoform in at least this this minimum number of samples for the gene isoform to be included in the analysis\n*\tmin_gene_expr = 10 - the minimum number of total mapped sequence reads for a gene to be considered expressed\n*\tmin_feature_expr = 3 - the minimum number of total mapped sequence reads for a gene isoform to be considered\n\n### edgeR based differential expression analysis\n+A statistical analysis is first performed using edgeR (Robinson, McCarthy, and Smyth (2010), McCarthy et al. (2012)) to identify the subset of differentially expressed genes. The filtered list of gene counts is used as input. A normalisation factor is calculated for each sequence library (using the default TMM method - please see McCarthy et al. (2012) for further details). The defined experimental design is used to calculate estimates of dispersion for each of the gene features. Statistical tests are calculated using the contrasts defined in the experimental design. The differentially expressed genes are corrected for false discovery (fdr) using the method of Benjamini & Hochberg (Benjamini and Hochberg (1995))\n\n### Differential transcript usage using DEXSeq\nDifferential transcript usage analysis is performed using the R DEXSeq package (Reyes et al. (2013)). Similar to the edgeR package, DEXSeq estimates the variance between the biological replicates and applies generalised linear models for the statistical testing. The key difference is that the DEXSeq method looks for differences at the exon count level. DEXSeq uses the filtered transcript count data prepared earlier in this analysis. \n\n### StageR stage-wise analysis of DGE and DTU\nThe final component of this isoform analysis is a stage-wise statistical test using the R software package `stageR` (Van den Berge and Clement (2018)). stageR uses (1) the raw p-values for DTU from the DEXSeq analysis in the previous section and (2) a false-discovery corrected set of p-values from testing whether individual genes contain at least one exon showing DTU. A hierarchical two-stage statistical testing evaluates the set of genes for DTU.\n\n## Running the workflow\nFor the differential expression analysis section you should have at least 3 repeats for each sample. \nYour FASTQ data will need to be organised in to 6 directories that represent 3 repeats for each condition. \n\n\n## Analysis \nDifferential gene expression is sensitive to the input data quantity and quality. There should be equivalence between samples in the number of sequence reads, mapped reads and quality scores. The sequence and alignment summary plots in the report can be used to assess these metrics. There is also a table that shows the transcript per million(TPM) calculated from the salmon counts. TPM normalizes the data for gene length and then sequencing depth, and makes it easier to compare across samples compared to counts.\n\n### Workflow inputs\n- Directory containing cDNA/direct RNA reads. Or a directory containing subdirectories each with reads from different samples\n (in fastq/fastq.gz format)\n- Reference genome in fasta format (required for reference-based assembly).\n- Optional reference annotation in GFF2/3 format (extensions allowed are .gtf(.gz), .gff(.gz), .gff3(.gz)) (required for differential expression analysis `--de_analysis`). Only annotation files from [Encode](https://www.encodeproject.org), [Ensembl](https://www.ensembl.org/index.html) and [NCBI](https://www.ncbi.nlm.nih.gov/) are supported.\n- For fusion detection, JAFFAL reference files (see Quickstart) \n", + "links": "## Useful links\n\n* [nextflow](https://www.nextflow.io/)\n* [docker](https://www.docker.com/products/docker-desktop)\n* [Singularity](https://sylabs.io/singularity/)\n* [racon](https://github.com/isovic/racon)" } } \ No newline at end of file diff --git a/subworkflows/denovo_assembly.nf b/subworkflows/denovo_assembly.nf deleted file mode 100644 index e7fff94..0000000 --- a/subworkflows/denovo_assembly.nf +++ /dev/null @@ -1,285 +0,0 @@ -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 { - /* - Build coding - */ - label "isoforms" - input: - tuple val(sample_id), path(cluster_fq) - output: - tuple val(sample_id), path("*final_polished_cds.fa"), emit: polished_cds, optional: true - 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 - tmpcons=\${cln}_tmcons.fa - - minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$spoa_cons \$sample | gzip - > \$samgz - exitcode=0 - racon -t ${params.threads} --no-trimming -u -w 2000 \$sample \$samgz \$spoa_cons > \$racon_cons || exitcode=1 - if [[ \$exitcode -eq 0 ]]; - then - # Rename consensus sequence name with cluster id - cat \$racon_cons | seqkit replace -p ".*" -r cluster_\${cln} > \$tmpcons - mv \$tmpcons \$racon_cons - else - echo "Polishing failed for \${cln}" - continue - fi - - # polish 2 - minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons \$sample | gzip - > \$samgz - exitcode=0 - racon -t ${params.threads} -u --no-trimming \$sample \$samgz \$racon_cons > \$tmpcons || exitcode=1 - if [[ \$exitcode -eq 0 ]]; - then - echo "success polish 2" - mv \$tmpcons \$racon_cons - else - echo "Polishing failed for \${cln}" - continue - fi - - # polish 3 - minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons \$sample | gzip - > \$samgz - exitcode=0 - racon -t ${params.threads} -u \$sample \$samgz \$racon_cons > \$tmpcons || exitcode=1 - if [[ \$exitcode -eq 0 ]]; - then - cat \$tmpcons >> ${sample_id}_\${UNID}_final_polished_cds.fa - echo "polishing 3 success" - else - echo "Polishing failed for \${cln}" - 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: - """ - 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: - """ - minimap2 -t ${params.threads} \ - -ax splice ${params.minimap2_opts} $polished_cds $sorted_reads_dir/sorted_reads.fastq |\ - samtools view -b - |\ - samtools sort -o "${sample_id}_reads_aln_sorted.bam"; - samtools index "${sample_id}_reads_aln_sorted.bam"; - ((seqkit bam -s -j ${params.threads} "${sample_id}_reads_aln_sorted.bam" 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 - """ - b=0 - if [ ${params.isOnClust2_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 - echo "Num bases: \$nr_bases"; - else - batch_size=${params.isOnClust2_batch_size} - fi - - echo "Batch size:\$batch_size"; - - - mkdir -p sorted; isONclust2 sort $params.isOnClust2_sort_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: - """ - workflow-glue run_isonclust2 --workdir . - """ -} - -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: - path reference - tuple val(sample_id), path(reads_fl), path(final_clusters_dir) - 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" - - """ - 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; - workflow-glue compute_cluster_quality --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.out.sorted_batches) - - dump_clusters( - clustering.out.root_cluster - .join(make_batches.out.sorted_reads_dir)) - - build_backbones( - dump_clusters.out.final_clusters - .flatMap(map_sample_ids_cls) - .groupTuple(size: 10, remainder: true)) - - merge_cds( - build_backbones.out.polished_cds - .flatMap(map_sample_ids_cls) - .groupTuple()) - - cds_align( - merge_cds.out.final_polished_cds - .join(make_batches.out.sorted_reads_dir)) - if (params.ref_genome) { - cluster_quality(reference, fastq_reads_fl - .join(dump_clusters.out.final_clusters_dir)) - - cluster_quality.out.cluster_qc_dir.set { opt_qual_ch } - - cluster_quality.out.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.out.bam - cds = merge_cds.out.final_polished_cds - stats = cds_align.out.stats - opt_qual_ch - opt_qual_raw_ch -} -