diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2a80d3a..e4ff3b8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,10 @@ All notable changes to this project will be documented in this file.
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+# Unreleased
+## Added
+- --direct_rna option.
+
## [v0.1.1]
## Fixed
- Incorrect numbers and of transcripts caused by merging gff files with same gene and transcript ids
diff --git a/README.md b/README.md
index d8cb0ec..b6eff18 100644
--- a/README.md
+++ b/README.md
@@ -8,9 +8,9 @@ It has been adapted from two existing Snakemake pipelines:
* https://github.com/nanoporetech/pipeline-nanopore-denovo-isoforms
---
## Overview
-* cDNA or direct RNA reads are initially and optionally preprocessed by [pychopper](https://github.com/nanoporetech/pychopper)
-for identification of full-length reads, as well as trimming and orientation correction.
-
+* cDNA reads are initially preprocessed by [pychopper](https://github.com/nanoporetech/pychopper)
+for identification of full-length reads, as well as trimming and orientation correction (This step is omited for
+ direct RNA reads)
Reference-based approach
* Full length reads are mapped to a supplied reference genome using [minimap2](https://github.com/lh3/minimap2)
@@ -22,11 +22,12 @@ using [gffcompare](http://ccb.jhu.edu/software/stringtie/gffcompare.shtml)
de novo-based approach (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
+ 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
+* Full-length reads are then mapped to these polished CDS.
+* Transcripts are assembled by stringtie as for the reference-based approach.
+* Note: This is currently not supported with direct RNA reads.
For both approaches:
@@ -67,7 +68,7 @@ to see the options for the workflow.
- Optional reference annotation in GFF2/3 format
**Example execution of a workflow for reference-based transcript assembly**
-This uses a synthetic SIRV dataset so we need to tell minimap2 about the non-canonical spplice junctions with
+This uses a synthetic SIRV dataset, so we need to tell minimap2 about the non-canonical splice junctions with
--minimap2_opts '-uf --splice-flank=no'
```
OUTPUT=~/output;
@@ -94,6 +95,8 @@ or at the command line:
- Threshold for including isoforms into interactive table `transcript_table_cov_thresh = 50`
- Run the denovo pipeline `denovo = true` (default false)
+- To run the workflow with direct RNA reads `--direct_rna` (skips the pychopper step).
+
Pychopper and minimap2 can take options via `minimap2_opts` and `pychopper_opts`
@@ -107,7 +110,7 @@ For example:
- pychopper can use one of two available backends for identifying primers in the raw reads
- nhmmscan `pychopper opts = '-m phmm'`
- edlib `pychopper opts = '-m edlib'`
-
Note: edlib is used in the configs as it's quite a lot faster. However it may be less sensitive than nhmmscan.
+
Note: edlib is set by default in the config as it's quite a lot faster. However it may be less sensitive than nhmmscan.
---
## Workflow outputs
@@ -129,9 +132,7 @@ Each sample will also have it's own directory containing the following (dependen
* {sample_id}_transcriptome.fas
- A transcriptome derived from the query reads
* {sample_id}_merged_transcriptome.fas
- - A transcriptome derived from the query reads + reference annotation
-* merged_transcriptome.fas
- -A transcriptome made from the combined query reads and reference annotation
+ - A transcriptome derived from the query reads + reference annotation
## Useful links
diff --git a/bin/check_sample_sheet.py b/bin/check_sample_sheet.py
index b8ce989..9c6a4ff 100755
--- a/bin/check_sample_sheet.py
+++ b/bin/check_sample_sheet.py
@@ -21,8 +21,7 @@ def main():
"Warning: sample sheet contains both 'alias' and "
'sample_id, using the former.')
samples['sample_id'] = samples['alias']
- if 'barcode' not in samples.columns \
- or 'sample_id' not in samples.columns:
+ if not set(['sample_id', 'barcode']).intersection(samples.columns):
raise IOError()
except Exception:
raise IOError(
diff --git a/bin/compute_cluster_quality.py b/bin/compute_cluster_quality.py
index acd40fc..5ebb20f 100755
--- a/bin/compute_cluster_quality.py
+++ b/bin/compute_cluster_quality.py
@@ -15,8 +15,9 @@ from matplotlib import pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
import pandas as pd
import pysam
-from sklearn.metrics.cluster import adjusted_rand_score, completeness_score,\
- homogeneity_score, v_measure_score
+from sklearn.metrics.cluster import (
+ adjusted_rand_score, completeness_score,
+ homogeneity_score, v_measure_score)
matplotlib.use('Agg')
@@ -368,14 +369,26 @@ def get_cluster_information(clusters, classes):
print("MIXED:", "Tot classes containing both:", len(
set(clustered_classes.keys()) & set(not_clustered_classes.keys())))
print("Total number of classes (unique gene ID):", total_nr_classes)
- return (total_nr_classes - len(singleton_classes), len(singleton_classes),
- min_class_size, max_class_size, mean_class_size, median_class_size,
- total_nr_clusters, len(singleton_clusters)
- + len(omitted_from_output_singletons),
- min_cluster_size, max_cluster_size, mean_cluster_size,
- median_cluster_size, len(unaligned_but_nontrivially_clustered),
- upper_75_class_size, upper_75_cluster_size, e_class_size,
- n50_class_size, e_cluster_size, n50_cluster_size)
+ return (
+ total_nr_classes - len(singleton_classes),
+ len(singleton_classes),
+ min_class_size,
+ max_class_size,
+ mean_class_size,
+ median_class_size,
+ total_nr_clusters,
+ len(singleton_clusters) + len(omitted_from_output_singletons),
+ min_cluster_size,
+ max_cluster_size,
+ mean_cluster_size,
+ median_cluster_size,
+ len(unaligned_but_nontrivially_clustered),
+ upper_75_class_size,
+ upper_75_cluster_size,
+ e_class_size,
+ n50_class_size,
+ e_cluster_size,
+ n50_cluster_size)
def main(args):
@@ -400,14 +413,16 @@ def main(args):
v_score, compl_score, homog_score, clustered_but_unaligned, ari = \
compute_V_measure(clusters, classes)
- nr_non_singleton_classes, singleton_classes, min_class_size, \
- max_class_size, mean_class_size, median_class_size, total_nr_clusters,\
- singleton_clusters, min_cluster_size, max_cluster_size, \
- mean_cluster_size, median_cluster_size, \
- unaligned_but_nontrivially_clustered, \
- upper_75_class_size, upper_75_cluster_size, e_class_size, \
- n50_class_size, e_cluster_size, n50_cluster_size = \
- get_cluster_information(clusters, classes)
+ (
+ nr_non_singleton_classes, singleton_classes, min_class_size,
+ max_class_size, mean_class_size, median_class_size, total_nr_clusters,
+ singleton_clusters, min_cluster_size, max_cluster_size,
+ mean_cluster_size, median_cluster_size,
+ unaligned_but_nontrivially_clustered,
+ upper_75_class_size, upper_75_cluster_size, e_class_size,
+ n50_class_size, e_cluster_size,
+ n50_cluster_size
+ ) = get_cluster_information(clusters, classes)
outfile = open(args.outfile, "w")
@@ -519,10 +534,10 @@ def main(args):
upper_75_class_size,
median_cluster_size,
median_class_size]}).set_index('Statistic')
- dfs2 = pd.DataFrame({'Statistic': ['N50ClsSize', 'N50ClassSize'],
- 'Value': [
- n50_cluster_size,
- n50_class_size]}).set_index('Statistic')
+ dfs2 = pd.DataFrame(
+ {'Statistic': ['N50ClsSize', 'N50ClassSize'],
+ 'Value': [n50_cluster_size, n50_class_size]
+ }).set_index('Statistic')
rdo = Path(args.raw_data_out)
dfc.to_csv(rdo / 'v_ari_com_hom.csv')
diff --git a/bin/generate_tracking_summary.py b/bin/generate_tracking_summary.py
index 7bb9791..0ef09b6 100755
--- a/bin/generate_tracking_summary.py
+++ b/bin/generate_tracking_summary.py
@@ -60,8 +60,8 @@ def main(args):
assert os.path.isdir(args.output_dir)
if args.annotation:
os.path.isfile(args.annotation)
- generate_tracking_summary(args.tracking, output_dir=args.output_dir,
- annotations=args.annotation)
+ generate_tracking_summary(
+ args.tracking, output_dir=args.output_dir, annotations=args.annotation)
if __name__ == '__main__':
diff --git a/bin/report.py b/bin/report.py
index 75c235f..1c1d54f 100755
--- a/bin/report.py
+++ b/bin/report.py
@@ -23,52 +23,6 @@ import pandas as pd
import sigfig
-def _vbar(x, top, title, **kwargs):
- """Vertical bar chart."""
- fig = figure(title=title)
-
- fig.vbar(x, top=top, **kwargs)
- return fig
-
-
-def _hbar(y, right, title='', fig_height=300, fig_width=300,
- bar_height=0.1, **kwargs):
- """Horizontal bar chart."""
- fig = figure(title=title, height=fig_height, width=fig_width)
- yn = list(range(len(y)))
- fig.hbar(yn,
- right=right,
- height=bar_height,
- **kwargs)
- # Overide the numerical labels with cate
- fig.yaxis.ticker = yn
- mapper = {k: v for (k, v) in zip(yn, y)}
- fig.yaxis.major_label_overrides = mapper
-
- return fig
-
-
-def simple_hbar(df, y, right, title="", color=Colors.cerulean,
- fig_kwargs={}, plot_kwargs={}):
- """Create a simple barplot.
-
- :param groups: the grouping variable (the x-axis values).
- :param values: the data for bars are drawn (the y-axis values).
- :param kwargs: kwargs for bokeh figure.
- """
- defaults = {
- 'output_backend': 'webgl',
- 'height': 300, 'width': 600}
- defaults.update(fig_kwargs)
- p = figure(y_range=df[y], height=250, title=title,
- toolbar_location=None, tools="")
-
- plot_kwargs.update({'height': 0.2})
- p.hbar(y=df[y], right=df[right], **plot_kwargs)
-
- return p
-
-
def _parse_stat_line(sl):
"""Parse a stats line."""
res = {}
@@ -288,8 +242,9 @@ def grouped_bar(df, title="", tilted_xlabs=False):
df = df.reset_index(drop=True)
source = ColumnDataSource(data=df)
- p = figure(x_range=df['x_groups'], y_range=yrange, height=250, title=title,
- toolbar_location=None, tools="")
+ p = figure(
+ x_range=df['x_groups'], y_range=yrange, height=250, title=title,
+ toolbar_location=None, tools="")
i = 0
# Use the dodge method to plot groups of bars
# https://docs.bokeh.org/en/latest/docs/user_guide/categorical.html
@@ -304,8 +259,9 @@ def grouped_bar(df, title="", tilted_xlabs=False):
for col in df.columns:
num_colors = df.shape[1] - 1
- colors = list(zip(*[[Category10_10[x]] * (len(df.columns) - 1)
- for x in range(num_colors)]))
+ colors = list(zip(
+ *[[Category10_10[x]] * (len(df.columns) - 1)
+ for x in range(num_colors)]))
colors = [item for sublist in colors for item in sublist]
if col == 'x_groups':
@@ -314,9 +270,9 @@ def grouped_bar(df, title="", tilted_xlabs=False):
i += 1
width = df.size / 60
- p.vbar(x=dodge('x_groups', current_dodge, range=p.x_range), top=col,
- width=width, source=source, color=color,
- legend_label=col)
+ p.vbar(
+ x=dodge('x_groups', current_dodge, range=p.x_range), top=col,
+ width=width, source=source, color=color, legend_label=col)
current_dodge += dodge_increment
p.x_range.range_padding = 0.1
@@ -341,8 +297,7 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
# Plot overview panel:
section = report.add_section()
- # TODO: update this based on current version
- section.markdown('''
+ gffcompare_md = ('''
### Annotation summary
The following plots summarize some of the output from
@@ -370,19 +325,26 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
stats, _, miss, novel, total = \
parse_gffcmp_stats(dir_ / 'str_merged.stats')
- bar_totals = grouped_bar(total, title="Totals")
- bar_performance = grouped_bar(stats, title="Performance",
- tilted_xlabs=True)
- bar_missed = grouped_bar(miss, title="Missed")
- bar_novel = grouped_bar(novel, title="Novel")
+ if not any([x.empty for x in [stats, miss, novel, total]]):
+ bar_totals = grouped_bar(total, title="Totals")
+ bar_performance = grouped_bar(
+ stats, title="Performance", tilted_xlabs=True)
+ bar_missed = grouped_bar(miss, title="Missed")
+ bar_novel = grouped_bar(novel, title="Novel")
+ tabs.append(Panel(
+ child=gridplot(
+ [bar_totals, bar_performance, bar_missed, bar_novel],
+ ncols=2, width=350, height=260), title=id_))
- tabs.append(Panel(
- child=gridplot(
- [bar_totals, bar_performance, bar_missed, bar_novel],
- ncols=2, width=350, height=260), title=id_))
-
- cover_panel = Tabs(tabs=tabs)
- section.plot(cover_panel)
+ cover_panel = Tabs(tabs=tabs)
+ section.markdown(gffcompare_md)
+ section.plot(cover_panel)
+ else:
+ gffcompare_md = ('''
+ __Warning__: Gffcompare summary cannot be shown.
+ This could be due to incompatible reference fasta and gff
+ files.''')
+ section.markdown(gffcompare_md)
names = {
'=': 'ExactMatch:=',
@@ -414,18 +376,19 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
[This diagram](https://ccb.jhu.edu/software/stringtie/
gffcompare_codes.png) illustrates the different classes.
- ''')
+ ''')
tracking_dfs = []
print(gffcompare_outdirs)
track_files = [x / 'str_merged.tracking' for x in gffcompare_outdirs]
- df_tracking = load_sample_data(track_files, sample_ids,
- read_func=lambda x:
- pd.read_csv(x, sep="\t", header=None,
- usecols=[0, 3],
- names=['Count', 'Overlaps']))
+ df_tracking = load_sample_data(
+ track_files, sample_ids,
+ read_func=lambda x: pd.read_csv(
+ x, sep="\t", header=None,
+ usecols=[0, 3],
+ names=['Count', 'Overlaps']))
tabs = []
for id_, df_track in df_tracking.groupby('sample_id'):
@@ -433,9 +396,10 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
tracking.Overlaps = tracking.Overlaps.map(names)
tracking['Percent'] = tracking.Count * 100 / tracking.Count.sum()
tracking = tracking.sort_values("Overlaps")
- track_bar = _hbar(
+ track_bar = bars.simple_hbar(
tracking['Overlaps'].values.tolist(),
- tracking['Percent'].values.tolist(), title="{}".format(id_))
+ tracking['Percent'].values.tolist(),
+ colors=Colors.cerulean, title=id_)
tracking_dfs.append(tracking)
@@ -448,13 +412,12 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
tracking.drop(columns=['sample_id', 'Overlaps'], inplace=True)
tracking = tracking[['Code', 'Description', 'Count', 'Percent']]
- cols = [TableColumn(field=Ci, title=Ci, width=100)
- for Ci in tracking.columns]
+ cols = [TableColumn(
+ field=Ci, title=Ci, width=100) for Ci in tracking.columns]
- track_table = DataTable(columns=cols,
- source=ColumnDataSource(tracking),
- index_position=None,
- width=500)
+ track_table = DataTable(
+ columns=cols, source=ColumnDataSource(tracking),
+ index_position=None, width=500)
tabs.append(Panel(
child=gridplot([track_bar, track_table], ncols=2), title=id_)
)
@@ -462,13 +425,12 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
cover_panel = Tabs(tabs=tabs)
section.plot(cover_panel)
- def plot_isoforms_per_tpm_bin(df_code, class_code, sample_id,
- geomspace=False):
+ def plot_isoforms_per_tpm_bin(
+ df_code, class_code, sample_id, geomspace=False):
"""Make plots of number of isoforms per TPM coverage bin."""
max_ = int(sigfig.round(df_code.TPM.max(), 2))
if geomspace:
- bins = [math.ceil(x) for x in
- np.geomspace(10, max_, num=15)]
+ bins = [math.ceil(x) for x in np.geomspace(10, max_, num=15)]
else:
bins = np.linspace(10, max_, 15)
bins = np.unique(bins) # Low max_ can end up with duplicated bins
@@ -527,11 +489,11 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
return df_tmap
-def pychopper_plots(report, df):
+def pychopper_plots(report, pychop_report):
"""Make plots from pychopper.cdna_classifier.py.
:param report: aplanat WFReport
- :param df: DataFrame of Pychopper stats
+ :param pychop_report: path to pychopper stats file
"""
section = report.add_section()
section.markdown('''
@@ -548,12 +510,15 @@ def pychopper_plots(report, df):
''')
plots = []
+ df = pd.read_csv(pychop_report, sep='\t', index_col=0)
+
for id_, df in df.groupby('sample_id'):
df1 = df.set_index('Name', drop=True)
df1 = df1.T[['Primers_found', 'Rescue', 'Unusable']]
- df1.rename(columns={'Primers_found': 'Pr.found',
- 'Rescue': 'Resc',
- 'Unsable': 'Un'}, inplace=True)
+ df1.rename(columns={
+ 'Primers_found': 'Pr.found',
+ 'Rescue': 'Resc',
+ 'Unusable': 'Un'}, inplace=True)
df2 = df[df.index == 'Strand']
bar_chop = bars.simple_bar(
df1.columns.values.tolist() + df2.Name.values.tolist(),
@@ -562,8 +527,8 @@ def pychopper_plots(report, df):
colors=Colors.cerulean)
plots.extend([bar_chop])
- grid = gridplot(plots, ncols=4,
- width=300, height=300)
+ grid = gridplot(
+ plots, ncols=4, width=300, height=300)
section.plot(grid)
@@ -602,8 +567,7 @@ def cluster_quality(cluster_qc_dir, report, sample_ids):
tabs = []
for id_, cluster_dir in zip(sample_ids, cluster_qc_dir):
plots = []
- for fn in [
- 'v_ari_com_hom.csv', 'sing_nonsing.csv']:
+ for fn in ['v_ari_com_hom.csv', 'sing_nonsing.csv']:
# Skip the next two plots for now
# 'class_sizes1.csv', 'class_sizes2.csv']:
df = pd.read_csv(Path(cluster_dir) / fn)
@@ -629,8 +593,9 @@ def transcript_table(report, df_tmaps, covr_threshold):
# all in single table and sample_id column? Currently it's the latter
# drop some columns for the big table and do some filtering
- df = df_tmaps.drop(columns=['FPKM', 'qry_gene_id', 'major_iso_id',
- 'ref_match_len', 'TPM'])
+ df = df_tmaps.drop(
+ columns=[
+ 'FPKM', 'qry_gene_id', 'major_iso_id', 'ref_match_len', 'TPM'])
df.sort_values('cov', ascending=True, inplace=True)
counts = list(range(len(df)))
@@ -672,7 +637,7 @@ def transcript_table(report, df_tmaps, covr_threshold):
section.table(df, index=False)
-def tanscriptome_summary(report, gffs, sample_ids, denovo=False):
+def transcriptome_summary(report, gffs, sample_ids, denovo=False):
"""
Plot transcriptome summaries.
@@ -745,12 +710,12 @@ def tanscriptome_summary(report, gffs, sample_ids, denovo=False):
if not denovo:
x, y = zip(*sorted(exons_per_transcript.items()))
- ept_bar = _vbar(x, list(y),
- title="Exons per transcript",
- color=Colors.cerulean)
+ fig = figure(title="Exons per transcript")
+ fig.vbar(
+ x, top=list(y), color=Colors.cerulean)
- ept_bar.xaxis.major_label_orientation = math.pi / 2.8
- plots.append(ept_bar)
+ fig.xaxis.major_label_orientation = math.pi / 2.8
+ plots.append(fig)
df_sum = pd.DataFrame.from_dict(
{'Total genes': [num_genes],
@@ -760,12 +725,11 @@ def tanscriptome_summary(report, gffs, sample_ids, denovo=False):
df_sum.reset_index(drop=False, inplace=True)
df_sum.columns = [' ', 'count']
- cols = [TableColumn(field=Ci, title=Ci, width=80)
- for Ci in df_sum.columns]
- data_table = DataTable(columns=cols,
- source=ColumnDataSource(df_sum),
- index_position=None,
- width=180)
+ cols = [TableColumn(
+ field=Ci, title=Ci, width=80) for Ci in df_sum.columns]
+ data_table = DataTable(
+ columns=cols, source=ColumnDataSource(df_sum),
+ index_position=None, width=180)
plots.append(data_table)
tabs.append(Panel(
@@ -808,7 +772,6 @@ def main():
parser.add_argument(
"--commit", default='unknown',
help="git commit of the executed workflow")
-
parser.add_argument(
"--alignment_stats", required=False, default=None, nargs='*',
help="TSV summary file of alignment statistics")
@@ -819,7 +782,7 @@ def main():
"--gffcompare_dir", required=False, default=None, nargs='*',
help="gffcompare outout dir")
parser.add_argument(
- "--pychop_report", required=True, nargs='+',
+ "--pychop_report", required=False, default=None,
help="TSV summary file of pychopper statistics")
parser.add_argument(
"--sample_ids", required=True, nargs='+',
@@ -860,8 +823,8 @@ def main():
section.table(df_aln_stats)
# workflow-specific plotting
- tanscriptome_summary(report, args.gff_annotation, sample_ids,
- denovo=args.denovo)
+ transcriptome_summary(
+ report, args.gff_annotation, sample_ids, denovo=args.denovo)
df_tmaps = gff_compare_plots(
report,
@@ -870,17 +833,8 @@ def main():
report.write(args.report)
- pc_df = pd.DataFrame()
- for id_, pyc in zip(sample_ids, args.pychop_report):
- try:
- p = pd.read_csv(pyc, sep='\t', index_col=0)
- except pd.errors.EmptyDataError:
- continue
- p['sample_id'] = id_
- pc_df = pd.concat([pc_df, p])
-
- if len(pc_df) > 0:
- pychopper_plots(report, pc_df)
+ if args.pychop_report is not None:
+ pychopper_plots(report, args.pychop_report)
if df_tmaps is not None:
transcript_table(report, df_tmaps, args.transcript_table_cov_thresh)
diff --git a/bin/run_isonclust2.py b/bin/run_isonclust2.py
index 7560dbc..a7fcff9 100755
--- a/bin/run_isonclust2.py
+++ b/bin/run_isonclust2.py
@@ -24,8 +24,8 @@ class Node:
def __repr__(self):
"""Get string repr of a node."""
- return "Node:{} Level: {} File: {} Done: {} Left: {} Right: " \
- "{} Parent: {}".format(
+ 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,
@@ -47,8 +47,10 @@ 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]
+ 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]):
@@ -93,12 +95,14 @@ def main():
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'
+ 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)
@@ -109,14 +113,15 @@ def main():
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)
+ 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)
+ sub.call((
+ "ln -s `realpath clusters/isONcluster_{}.cer` "
+ "isONcluster_ROOT.cer".format(n.Id)), shell=True)
if __name__ == '__main__':
diff --git a/environment.yaml b/environment.yaml
index fd72183..b097162 100644
--- a/environment.yaml
+++ b/environment.yaml
@@ -6,7 +6,7 @@ channels:
- defaults
dependencies:
- python==3.8.*
- - aplanat >=0.6.2
+ - aplanat >=0.6.4
- epi2melabs
- minimap2 ==2.24
- samtools ==1.14
diff --git a/main.nf b/main.nf
index 0e559f3..bf1184e 100644
--- a/main.nf
+++ b/main.nf
@@ -84,12 +84,16 @@ process preprocess_reads {
tuple val(sample_id), path(input_reads)
output:
tuple val(sample_id), path("${sample_id}_full_length_reads.fq"), emit: full_len_reads
- tuple val(sample_id), path('*.tsv'), emit: report
+ path '*.tsv', emit: report
script:
"""
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 .
+
+ # Add sample id column
+ sed "1s/\$/\tsample_id/; 1 ! s/\$/\t${sample_id}/" ${sample_id}_cdna_classifier_report.tsv > tmp
+ mv tmp ${sample_id}_cdna_classifier_report.tsv
"""
}
@@ -266,12 +270,12 @@ process makeReport {
path versions
path "params.json"
val denovo
+ path pychopper_report
tuple val(sample_ids),
- path(seq_summaries),
- path(aln_stats),
- path(gffcmp_dir),
- path(cdna_class_report),
- path(gff_annotation)
+ path(seq_summaries),
+ path(aln_stats),
+ path(gffcmp_dir),
+ path(gff_annotation)
output:
path("wf-isoforms-*.html"), emit: report
script:
@@ -280,12 +284,13 @@ process makeReport {
def report_name = "wf-isoforms-report.html"
def OPT_ALN = denovo ? '' : "--alignment_stats ${aln_stats}"
def OPT_DENOVO = denovo ? "--denovo" : ''
+ def OPT_PC_REPORT = pychopper_report.name.startsWith('OPTIONAL_FILE') ? '' : "--pychop_report ${pychopper_report}"
"""
report.py --report $report_name \
--versions $versions \
--params params.json \
$OPT_ALN \
- --pychop_report $cdna_class_report \
+ $OPT_PC_REPORT \
--sample_ids $sids \
--summaries $seq_summaries \
--gffcompare_dir $gffcmp_dir \
@@ -345,16 +350,25 @@ workflow pipeline {
software_versions = getVersions()
workflow_params = getParams()
- preprocess_reads(summariseConcatReads.out.input_reads)
+
+ if (!params.direct_rna){
+ preprocess_reads(summariseConcatReads.out.input_reads)
+ full_len_reads = preprocess_reads.out.full_len_reads
+ pychopper_report = preprocess_reads.out.report.collectFile(keepHeader: true)
+ }
+ else{
+ full_len_reads = summariseConcatReads.out.input_reads
+ pychopper_report = file("$projectDir/data/OPTIONAL_FILE")
+ }
if (params.denovo){
println("Doing de novo assembly")
- m = denovo_assembly(preprocess_reads.out.full_len_reads, ref_genome)
+ m = denovo_assembly(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)
+ m = reference_assembly(build_minimap_index.out.index, ref_genome, full_len_reads)
}
split_bam(m.bam)
@@ -380,10 +394,10 @@ workflow pipeline {
software_versions,
workflow_params,
params.denovo,
+ pychopper_report,
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())
@@ -395,19 +409,18 @@ workflow pipeline {
.join(seq_for_transcriptome_build))
if (use_ref_ann){
- results = preprocess_reads.out.report
+ results = run_gffcompare.output.gffcmp_dir
.concat(
- run_gffcompare.output.gffcmp_dir,
- m.stats,
- get_transcriptome.out.flatMap(map_sample_ids_cls))
+ m.stats,
+ get_transcriptome.out.flatMap(map_sample_ids_cls))
.map {it -> it[1]}
.concat(makeReport.out.report)
}
if (!use_ref_ann && !params.denovo){
- results = preprocess_reads.out.report
- .concat(m.stats,
- get_transcriptome.out.flatMap(map_sample_ids_cls))
+ results = m.stats
+ .concat(
+ get_transcriptome.out.flatMap(map_sample_ids_cls))
.map {it -> it[1]}
.concat(makeReport.out.report)
@@ -415,16 +428,17 @@ workflow pipeline {
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
- })
+ 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
+ })
.map {it -> it[1]}
.concat(makeReport.out.report)
}
diff --git a/nextflow.config b/nextflow.config
index 7962535..33f4bf6 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -33,7 +33,8 @@ params {
show_hidden_params = false
schema_ignore_params = 'show_hidden_params,validate_params,monochrome_logs,aws_queue,aws_image_prefix,wfversion'
-
+ // Process cDNA reads using pychopper, turn off for direct RNA:
+ direct_rna = true
// Options passed to pychopper:
pychopper_opts = "-m edlib"
diff --git a/nextflow_schema.json b/nextflow_schema.json
index 392c5c8..5e31e34 100644
--- a/nextflow_schema.json
+++ b/nextflow_schema.json
@@ -66,6 +66,11 @@
"description": "Extra pychopper opts",
"default": "-m edlib"
},
+ "direct_rna": {
+ "type": "boolean",
+ "description": "Set to true for direct RNA sequencing. Omits the pychopper step.",
+ "default": false
+ },
"bundle_min_reads": {
"type": "integer",
"description": "Minimum size of bam bundle for parallel processing."