Merge branch 'direct_rna' into 'dev'

Re-add ability to do direct RNA

Closes CW-566

See merge request epi2melabs/workflow-containers/wf-isoforms!46

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

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) 
+* 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. <br>
+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) <br>
@@ -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
+* Full-length reads are then mapped to these polished CDS.
-* Transcripts are assembled by stringtie as for the reference-based approach
+* 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:<br>
 
 - 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` 
 <br>
@@ -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'`
-<br>Note: edlib is used in the configs as it's quite a lot faster. However it may be less sensitive than nhmmscan. 
+<br>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
@@ -130,8 +133,6 @@ Each sample will also have it's own directory containing the following (dependen
   - 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
   
 
 ## Useful links

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 \
+        if not set(['sample_id', 'barcode']).intersection(samples.columns):
-                or 'sample_id' not in samples.columns:
             raise IOError()
     except Exception:
         raise IOError(

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,\
+from sklearn.metrics.cluster import (
-    homogeneity_score, v_measure_score
+    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),
+    return (
-            min_class_size, max_class_size, mean_class_size, median_class_size,
+        total_nr_classes - len(singleton_classes),
-            total_nr_clusters, len(singleton_clusters)
+        len(singleton_classes),
-            + len(omitted_from_output_singletons),
+        min_class_size,
-            min_cluster_size, max_cluster_size, mean_cluster_size,
+        max_class_size,
-            median_cluster_size, len(unaligned_but_nontrivially_clustered),
+        mean_class_size,
-            upper_75_class_size, upper_75_cluster_size, e_class_size,
+        median_class_size,
-            n50_class_size, e_cluster_size, n50_cluster_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,\
+        nr_non_singleton_classes, singleton_classes, min_class_size,
-        singleton_clusters, min_cluster_size, max_cluster_size, \
+        max_class_size, mean_class_size, median_class_size, total_nr_clusters,
-        mean_cluster_size, median_cluster_size, \
+        singleton_clusters, min_cluster_size, max_cluster_size,
-        unaligned_but_nontrivially_clustered, \
+        mean_cluster_size, median_cluster_size,
-        upper_75_class_size, upper_75_cluster_size, e_class_size, \
+        unaligned_but_nontrivially_clustered,
-        n50_class_size, e_cluster_size, n50_cluster_size = \
+        upper_75_class_size, upper_75_cluster_size, e_class_size,
-        get_cluster_information(clusters, classes)
+        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'],
+    dfs2 = pd.DataFrame(
-                         'Value': [
+        {'Statistic': ['N50ClsSize', 'N50ClassSize'],
-                        n50_cluster_size,
+            'Value': [n50_cluster_size, n50_class_size]
-                        n50_class_size]}).set_index('Statistic')
+         }).set_index('Statistic')
 
     rdo = Path(args.raw_data_out)
     dfc.to_csv(rdo / 'v_ari_com_hom.csv')

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,
+    generate_tracking_summary(
-                              annotations=args.annotation)
+        args.tracking, output_dir=args.output_dir, annotations=args.annotation)
 
 
 if __name__ == '__main__':

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,7 +242,8 @@ 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,
+    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
@@ -304,7 +259,8 @@ 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)
+        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]
 
@@ -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,
+        p.vbar(
-               width=width, source=source, color=color,
+            x=dodge('x_groups', current_dodge, range=p.x_range), top=col,
-               legend_label=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
+    gffcompare_md = ('''
-    section.markdown('''
     ### 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')
 
+        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",
+            bar_performance = grouped_bar(
-                                      tilted_xlabs=True)
+                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_))
 
             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:=',
@@ -421,9 +383,10 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
     print(gffcompare_outdirs)
     track_files = [x / 'str_merged.tracking' for x in gffcompare_outdirs]
 
-    df_tracking = load_sample_data(track_files, sample_ids,
+    df_tracking = load_sample_data(
-                                   read_func=lambda x:
+        track_files, sample_ids,
-                                   pd.read_csv(x, sep="\t", header=None,
+        read_func=lambda x: pd.read_csv(
+            x, sep="\t", header=None,
             usecols=[0, 3],
             names=['Count', 'Overlaps']))
 
@@ -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)
+        cols = [TableColumn(
-                for Ci in tracking.columns]
+            field=Ci, title=Ci, width=100) for Ci in tracking.columns]
 
-        track_table = DataTable(columns=cols,
+        track_table = DataTable(
-                                source=ColumnDataSource(tracking),
+            columns=cols, source=ColumnDataSource(tracking),
-                                index_position=None,
+            index_position=None, width=500)
-                                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,
+    def plot_isoforms_per_tpm_bin(
-                                  geomspace=False):
+            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
+            bins = [math.ceil(x) for x in np.geomspace(10, max_, num=15)]
-                    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',
+        df1.rename(columns={
+            'Primers_found': 'Pr.found',
             'Rescue': 'Resc',
-                            'Unsable': 'Un'}, inplace=True)
+            '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,
+    grid = gridplot(
-                    width=300, height=300)
+        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 [
+        for fn in ['v_ari_com_hom.csv', 'sing_nonsing.csv']:
-                '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',
+    df = df_tmaps.drop(
-                                'ref_match_len', 'TPM'])
+        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),
+            fig = figure(title="Exons per transcript")
-                            title="Exons per transcript",
+            fig.vbar(
-                            color=Colors.cerulean)
+                x, top=list(y), color=Colors.cerulean)
 
-            ept_bar.xaxis.major_label_orientation = math.pi / 2.8
+            fig.xaxis.major_label_orientation = math.pi / 2.8
-            plots.append(ept_bar)
+            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)
+        cols = [TableColumn(
-                for Ci in df_sum.columns]
+            field=Ci, title=Ci, width=80) for Ci in df_sum.columns]
-        data_table = DataTable(columns=cols,
+        data_table = DataTable(
-                               source=ColumnDataSource(df_sum),
+            columns=cols, source=ColumnDataSource(df_sum),
-                               index_position=None,
+            index_position=None, width=180)
-                               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,
+    transcriptome_summary(
-                         denovo=args.denovo)
+        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()
+    if args.pychop_report is not None:
-    for id_, pyc in zip(sample_ids, args.pychop_report):
+        pychopper_plots(report, 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 df_tmaps is not None:
         transcript_table(report, df_tmaps, args.transcript_table_cov_thresh)

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: " \
+        return "Node:{} Level: {} File: {} Done: " \
-               "{} Parent: {}".format(
+            "{} Left: {} Right: {} Parent: {}".format(
                 self.Id, self.Level,
                 self.File, self.Done, self.Left.Id if
                 self.Left is not None else None,
@@ -47,7 +47,9 @@ 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))
+    batch_ids = [
+        int(re.search(
+            'batches/isONbatch_(.*)\\.cer$', x).group(1))
         for x in batches]
     LEVELS = OrderedDict()
     LEVELS[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/' \
+    init_template = (
-                    'isONbatch_{}.cer -o clusters/isONcluster_{}.cer {}; ' \
+        'isONclust2 cluster -x {} -v -Q -l batches/isONbatch_{}.cer '
-                    'sync;\n'
+        '-o clusters/isONcluster_{}.cer {}; '
-    template = 'isONclust2 cluster -x {} -v -Q -l clusters/isONcluster_{}' \
+        'sync;\n')
-               '.cer -r clusters/isONcluster_{}.cer -o clusters/isONcluster' \
+    template = (
-               '_{}.cer {}; sync\n'
+        '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,
+                    jr = template.format(
-                                         n.Right.Id, n.Id, purge)
+                        '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_{}"
+    sub.call((
-             ".cer` isONcluster_ROOT.cer".format(n.Id), shell=True)
+        "ln -s `realpath clusters/isONcluster_{}.cer` "
+        "isONcluster_ROOT.cer".format(n.Id)), shell=True)
 
 
 if __name__ == '__main__':

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

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,11 +270,11 @@ 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)
     output:
         path("wf-isoforms-*.html"), emit: report
@@ -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()
+
+        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,9 +409,8 @@ 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))
                       .map {it -> it[1]}
@@ -405,8 +418,8 @@ workflow pipeline {
 
        }
         if (!use_ref_ann && !params.denovo){
-            results = preprocess_reads.out.report
+            results =  m.stats
-                      .concat(m.stats,
+                      .concat(
                        get_transcriptome.out.flatMap(map_sample_ids_cls))
                       .map {it -> it[1]}
                       .concat(makeReport.out.report)
@@ -418,7 +431,8 @@ workflow pipeline {
                       seq_for_transcriptome_build,
                       get_transcriptome.out.flatMap(map_sample_ids_cls),
                       merge_gff_bundles.out.gff,
-                          m.opt_qual_ch.flatMap {it ->
+                      m.opt_qual_ch.flatMap {
+                          it ->
                           l = []
                           for (x in it[1..-1]){
                               l.add(tuple(it[0], x))

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"
 

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