kfukuta/wf-transcriptomes-v202
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

CW-2769 remove denovo

Browse Source
This commit is contained in:
Sarah Griffiths 2023-09-27 08:02:00 +00:00
parent 5d16332a2c
commit f683d44411
14 changed files with 36 additions and 598 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
.gitlab-ci.yml
View File

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

3
CHANGELOG.md
View File

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

28
README.md
View File

@ -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).
@ -298,6 +275,3 @@ in `${out_dir}/jaffal_output_${sample_id}` you will find:
* [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)

11
bin/workflow_glue/report.py
View File

@ -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,7 +765,6 @@ 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()))
fig = figure(title="Exons per transcript")
@ -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(

138
bin/workflow_glue/run_isonclust2.py
View File

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

1
docs/header.md
View File

@ -2,4 +2,3 @@
This repository contains a [nextflow](https://www.nextflow.io/) workflow
for assembly and annotation of transcripts from Oxford Nanopore cDNA or direct RNA reads.

10
docs/intro.md
View File

@ -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.

3
docs/links.md
View File

@ -4,6 +4,3 @@
* [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)

13
docs/quickstart.md
View File

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

12
evaluation/tests.sh
View File

@ -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' \

50
main.nf
View File

@ -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,12 +281,6 @@ 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"
@ -306,7 +296,7 @@ process run_gffcompare{
cp ${out_dir}/str_merged.annotated.gtf ${sample_id}_annotated.gtf
"""
}
}
process get_transcriptome{
@ -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"){
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)
}
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){

9
nextflow.config
View File

@ -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,

31
nextflow_schema.json
View File

File diff suppressed because one or more lines are too long

285
subworkflows/denovo_assembly.nf
View File

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