Merge de novo pipeline

CHANGELOG.md

@@ -4,40 +4,15 @@ 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).
 
-## [v0.0.7]
+## [0.1.0]
+## Added 
+- Added the denovo pipeline 
+## Changed
+- Updates to the report plots
+
+## [0.0.1]
 ### Added
-- Fastqingress module for common handling of (possibly
+- First release
-  multiplexed) inputs.
+- Initial port of Snakemake WF from https://github.com/nanoporetech/pipeline-nanopore-ref-isoforms
-- Optimized container size through removal of various
-  conda cruft.
-### Changed
-- Use mamba by default for building conda environments.
-- Cut down README to items specific to workflow.
-### Fixed
-- Incorrect specification of conda environment file in Nextflow config.
 
-## [v0.0.6]
-### Changed
-- Explicitely install into base conda env
 
-## [v0.0.5]
-### Added
-- Software versioning report example.
-
-## [v0.0.4]
-### Changed
-- Version bump to test CI.
-
-## [v0.0.3]
-### Changed
-- Moved all CI to templates.
-- Use canned aplanat report components.
-
-## [v0.0.2]
-### Added
-- CI release checks.
-- Create pre-releases in CI from dev branch.
-
-## [v0.0.1]
-
-First release.

Dockerfile

@@ -17,3 +17,4 @@ RUN \
 
 USER $WF_UID
 WORKDIR $HOME
+

README.md

@@ -1,16 +1,40 @@
-# Pipeline for annotating genomes using long read transcriptomics
+# wf-isoforms
+
 
 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. <br>
-
+It has been adapted from two existing Snakemake pipelines:
+* https://github.com/nanoporetech/pipeline-nanopore-ref-isoforms
+* https://github.com/nanoporetech/pipeline-nanopore-denovo-isoforms
+---
 ## Overview
-* cDNA or direct RNA reads are optionally preprocessed by [pychopper](https://github.com/nanoporetech/pychopper) for trimming and orientation. <br>
+* cDNA or direct RNA reads are initially and optionally preprocessed by [pychopper](https://github.com/nanoporetech/pychopper) 
-* Reads are then mapped to a supplied reference genome using [minimap2](https://github.com/lh3/minimap2) <br>
+for identification of full-length reads, as well as trimming and orientation correction. <br>
-* Transcripts are assembled by[stringtie](http://ccb.jhu.edu/software/stringtie) in long read mode (with or without a guide reference annotation) to generate the GFF annotation.
-* The annotation generated by the pipeline is compared to the reference annotation (if supplied) using [gffcompare](http://ccb.jhu.edu/software/stringtie/gffcompare.shtml) 
-* An html report is generated, which contains various summary statistics and plots of the data.
 
+
+Reference-based approach
+* Full length reads are mapped to a supplied reference genome using [minimap2](https://github.com/lh3/minimap2) <br>
+* Transcripts are assembled by [stringtie](http://ccb.jhu.edu/software/stringtie) 
+in long read mode (with or without a guide reference annotation) to generate the GFF annotation.
+* 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 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 
+* 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
+
+
+For both approaches:
+* An html report is generated, which contains various summary statistics and plots of the data.
+---
 ## Quickstart
+
+
 The workflow uses [nextflow](https://www.nextflow.io/) to manage compute and 
 software resources, as such nextflow will need to be installed before attempting
 to run the workflow.
@@ -25,7 +49,8 @@ either docker, singularity or conda is installed.
 It is not required to clone or download the git repository in order to run the workflow.
 For more information on running EPI2ME Labs workflows [visit out website](https://labs.epi2me.io/wfindex).
 
-**Workflow options**
+
+### Workflow options
 
 To obtain the workflow, having installed `nextflow`, users can run:
 
@@ -36,41 +61,91 @@ nextflow run epi2me-labs/wf-isoforms --help
 to see the options for the workflow.
 
 **Workflow inputs**
-- Directory containing cDNA/direct RNA reads (or path to single file) in fastq/fastq.gz format
+- Directory containing cDNA/direct RNA reads. Or a directory containing subdirectories each with reads from different samples
-- Reference genome in fast format
+  (in fastq/fastq.gz format)
+- Reference genome in fasta format (required for reference-based assembly)
 - Optional reference annotation in GFF2/3 format 
 
-**Example workflow run**
+**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 
+--minimap2_opts '-uf --splice-flank=no'
 ```
-# To run a small and quick example using synthetic data
+OUTPUT=~/output;
-nextflow run wf-isoforms/ --fastq test_data/fastq  --ref_genome genome.fasta --ref_annotation reference.gff
+nextflow run wf-isoforms/ --fastq test_data/fastq  --ref_genome test_data/SIRV_150601a.fasta --ref_annotation test_data/SIRV_isofroms.gtf
---out_dir outdir/ -profile conda
+--minimap2_opts '-uf --splice-flank=no' --out_dir outdir -w workspace_dir -profile conda -resume
 ```
 
 ```
 # To evaluate the workflow on a larger Drosophila dataset
-chmod u+x ./run_evaluation_dmel.sh outdir
+./evaluation/run_evaluation_dmel.sh outdir
 ```
 
-**Workflow outputs**
+**Example workflow for denovo transcript assembly**
+```
+OUTPUT=~/output
+nextflow run . --fastq test_data/fastq --denovo --ref_genome test_data/SIRV_150601a.fasta  -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
+--sample sample_id -resume
+```
+A full list of options can be seen in nextflow_schema.json. Below are some commonly used ones.
+They can be set in the config like this: <br>
+&nbsp;&nbsp; opt = 50<br>
+or at the command line:<br>
+&nbsp;&nbsp;  --opt 50
+
+- Threshold for including isoforms into interactive table `transcript_table_cov_thresh = 50`
+- Run the denovo pipeline `denovo = true` (default false)
+
+Pychopper and minimap2 can take options via `minimap2_opts` and `pychopper_opts` 
+<br>
+
+For example:
+- When using the SIRV synthetic test data  
+  - `minimap2_opts = '-uf --splice-flank=no'`
+- pychopper needs to know which cDNA synthesis kit used
+  - SQK-PCS109: use `pychopper_opts = '-k PCS109'` (default)
+  - SQK-PCS110: use `pychopper_opts = '-k PCS110'`
+- 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. 
+  
+---
+## Workflow outputs
 
 * wf-isoforms-report.html
   - Summary and plots of reads, alignments and the transcript assembly and annotation
-* str_merged.gff
+  - One file per run with results from each sample merged
-  - The stringtie-generated transcript annotations
+
-* str_merged.stats
+####Transcript isoform table
-  - gffcomare file with statistics regarding the accuracy of the assembled query transcripts in relation to the reference annotation 
+Note: Currently only available for the reference-guilded assembly
-* str_merged.annotated.gtf
+The report includes a searchable filterable table of transcript isoforms<br>
-  - A gffcomapre output file with extra columns relating to comparison with the reference annotation
+Note: If using a large dataset with 10s of thousands of predicted isoforms, to speed up searching, the table can be limited
-* str_transcriptome.fas
+by read coverage, which can be set with `transcript_table_cov_thresh = 50` (default 50x)
-  - A transcriptome made from the query reads 
+
+Each sample will also have it's own directory containing the following (dependent on assembly approach and input options)
+* {sampleid}_gffcompare
+  * Directory containing all output from gffcompare
+* transcripts_{sample_id}.gff
+  * Annotation of query transcripts
+* {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
-
+* final_polished_cds.fas (de novo only)
+  * CDS sequences derived from reads
   
+
+
 ## Useful links
 
 * [nextflow](https://www.nextflow.io/)
 * [docker](https://www.docker.com/products/docker-desktop)
 * [Singularity](https://sylabs.io/singularity/)
 * [conda](https://docs.conda.io/en/latest/miniconda.html)
+* [racon](https://github.com/isovic/racon)
+* [spoa](https://github.com/rvaser/spoa)
+* [inONclust](https://github.com/ksahlin/isONclust)
+* [isONclust2](https://github.com/nanoporetech/isONclust2)
+

bin/compute_cluster_quality.py

@@ -0,0 +1,638 @@
+#!/usr/bin/env python
+
+"""Generate cluster quality data."""
+
+# Adapted form script by Kristoffer Sahlin for
+# isONclust: https://github.com/ksahlin/isONclust
+
+import argparse
+from collections import defaultdict
+import math
+from pathlib import Path
+
+import matplotlib
+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
+matplotlib.use('Agg')
+
+
+def parse_inferred_clusters_tsv(tsv_file, args):
+    """parse_inferred_clusters_tsv."""
+    infile = open(tsv_file, "r")
+    infile.readline()
+    clusters = {}
+    for line in infile:
+        cluster_id, _, read_acc = line.strip().split("\t")
+        if args.simulated:
+            read_acc = "_".join([item for item in read_acc.split("_")[:-1]])
+        elif args.ont:
+            read_acc = read_acc.split(" ")[0]
+        elif args.modified_ont:
+            read_acc = read_acc
+        else:
+            read_acc = read_acc.split("_strand")[0]
+
+        clusters[read_acc] = int(cluster_id)
+    return clusters
+
+
+def parse_true_clusters(ref_file):
+    """parse_true_clusters."""
+    classes = defaultdict(dict)
+    ref_id_to_chrom = {}
+    alignment_counter = defaultdict(int)
+
+    prev_chrom = -1
+    curr_class_id = -1
+    prev_class_stop = -1
+    prev_read_id = ""
+    unique_reads = set()
+    unclassified = 0
+    for read in ref_file.fetch(until_eof=True):
+        unique_reads.add(read.query_name)
+        if read.is_unmapped:
+            unclassified += 1
+            continue
+        # deal with supplementary alignments!!
+        if read.is_secondary or read.is_supplementary:
+            continue
+        # print(read.query_name, read.flag)
+        assert prev_read_id != read.query_name
+
+        chrom = read.reference_name
+        if chrom != prev_chrom:
+            curr_class_id += 1
+            classes[read.query_name] = curr_class_id
+            prev_chrom = chrom
+            prev_class_stop = read.reference_end
+
+        else:
+            read_ref_start = read.reference_start
+            if read_ref_start > prev_class_stop:
+                curr_class_id += 1
+                classes[read.query_name] = curr_class_id
+                prev_class_stop = read.reference_end
+            else:
+                classes[read.query_name] = curr_class_id
+
+        prev_class_stop = max(read.reference_end, prev_class_stop)
+        prev_read_id = read.query_name
+
+        # classes[read.query_name] = int(read.reference_id)  # chrom
+        ref_id_to_chrom[int(read.reference_id)] = chrom
+        alignment_counter[int(read.reference_id)] += 1
+        # if chrom not in class_ranges:
+        #     class_ranges[chrom] = {}
+
+        # print(chrom, read_ref_start, read_ref_end)
+        # for start, stop in class_ranges[chrom]:
+        #     if start <= read_ref_start and  read_ref_end <= stop:
+        #         # entirly within
+        #     elif
+
+        # classes[read.query_name] =  #read.reference_name.split("|")[0]
+        # #"|".join(read.reference_name.split("|")[:2])
+    return classes, len(unique_reads), unclassified
+
+
+def parse_true_clusters_simulated(ref_file):
+    """parse_true_clusters_simulated."""
+    classes = defaultdict(dict)
+    for read in ref_file.fetch(until_eof=True):
+        # by gene id
+        classes[read.query_name] = read.reference_name.split("|")[0]
+        # by transcript id
+        # classes[read.query_name] = read.reference_name.split("|")[1]
+
+    return classes
+
+
+def compute_V_measure(clusters, classes):
+    """compute_V_measure."""
+    class_list, cluster_list = [], []
+    # not_found_id = 1000000
+    clustered_but_unaligned = 0
+    for read in clusters:
+        if read in classes:
+            class_list.append(classes[read])
+            cluster_list.append(clusters[read])
+        else:
+            # print("Read was clustered but unaligned:", read)
+            clustered_but_unaligned += 1
+
+    # added the unprocessed reads to the measure
+    not_clustered = set(classes.keys()) - set(clusters.keys())
+    highest_cluster_id = max(clusters.values())
+    highest_cluster_id += 1
+    for read in not_clustered:
+        class_list.append(classes[read])
+        cluster_list.append(highest_cluster_id)
+        highest_cluster_id += 1
+
+    v_score = v_measure_score(class_list, cluster_list)
+    compl_score = completeness_score(class_list, cluster_list)
+    homog_score = homogeneity_score(class_list, cluster_list)
+    ari = adjusted_rand_score(class_list, cluster_list)
+
+    print("Not included in clustering but aligned:", len(not_clustered))
+    print(
+        "V:",
+        v_score,
+        "Completeness:",
+        compl_score,
+        "Homogeneity:",
+        homog_score)
+    print(
+        "Nr reads clustered but unaligned "
+        "(i.e., no class and excluded from V-measure): ",
+        clustered_but_unaligned)
+    return v_score, compl_score, homog_score, clustered_but_unaligned, ari
+
+
+def compute_V_measure_non_singleton_classes(clusters, classes):
+    """V measure for non-singleton classes."""
+    max_cluster_id = max(clusters.values())
+    new_id = max_cluster_id + 1
+    classes_dict = {}
+    for read_acc, cl_id in classes.items():
+        if cl_id not in classes_dict:
+            classes_dict[cl_id] = [read_acc]
+        else:
+            classes_dict[cl_id].append(read_acc)
+
+    nontrivial_classes_reads = []
+    for cl_id in classes_dict:
+        if len(classes_dict[cl_id]) < 5:
+            continue
+        else:
+            for read in classes_dict[cl_id]:
+                nontrivial_classes_reads.append(read)
+
+    class_list, cluster_list = [], []
+    for read in nontrivial_classes_reads:
+        if read in clusters:
+            class_list.append(classes[read])
+            cluster_list.append(clusters[read])
+        else:
+            class_list.append(classes[read])
+            cluster_list.append(new_id)
+            new_id += 1
+
+    v_score = v_measure_score(class_list, cluster_list)
+    compl_score = completeness_score(class_list, cluster_list)
+    homog_score = homogeneity_score(class_list, cluster_list)
+    nr_filtered_classes = len(
+        [1 for cl_id in classes_dict if len(classes_dict[cl_id]) >= 5])
+    print(
+        "NONTRIVIAL CLASSES: V:",
+        v_score,
+        "Completeness:",
+        compl_score,
+        "Homogeneity:",
+        homog_score)
+    print("NUMBER OF CLASSES (FILTERED):", len(
+        [1 for cl_id in classes_dict if len(classes_dict[cl_id]) >= 5]))
+    return v_score, compl_score, homog_score, nr_filtered_classes
+
+
+def compute_V_measure_non_singletons(clusters, classes):
+    """V measure for non-singletons."""
+    cluster_dict = {}
+    for read_acc, cl_id in clusters.items():
+        if cl_id not in cluster_dict:
+            cluster_dict[cl_id] = [read_acc]
+        else:
+            cluster_dict[cl_id].append(read_acc)
+
+    nontrivial_clustered_reads = []
+    for cl_id in cluster_dict:
+        if len(cluster_dict[cl_id]) <= 1:
+            continue
+        else:
+            for read in cluster_dict[cl_id]:
+                nontrivial_clustered_reads.append(read)
+
+    class_list, cluster_list = [], []
+    # not_found_id = 1000000
+    clustered_but_unaligned = 0
+    for read in nontrivial_clustered_reads:
+        if read in classes:
+            class_list.append(classes[read])
+            cluster_list.append(clusters[read])
+        else:
+            # print("Read was clustered but unaligned:", read)
+            clustered_but_unaligned += 1
+
+    v_score = v_measure_score(class_list, cluster_list)
+    compl_score = completeness_score(class_list, cluster_list)
+    homog_score = homogeneity_score(class_list, cluster_list)
+    print(
+        "NONTRIVIAL CLUSTERS: V:",
+        v_score,
+        "Completeness:",
+        compl_score,
+        "Homogeneity:",
+        homog_score)
+    print(
+        "NONTRIVIAL CLUSTERS: Nr reads clustered but unaligned "
+        "(i.e., no class and excluded from V-veasure): ",
+        clustered_but_unaligned)
+    return v_score, compl_score, homog_score, clustered_but_unaligned
+
+
+def percentile(N, percent, key=lambda x: x):
+    """
+    Find the percentile of a list of values.
+
+    @parameter N - is a list of values. Note N MUST BE already sorted.
+    @parameter percent - a float value from 0.0 to 1.0.
+    @parameter key - optional key function to compute value
+    from each element of N.
+
+    @return - the percentile of the values
+    """
+    if not N:
+        return None
+    k = (len(N) - 1) * percent
+    f = math.floor(k)
+    c = math.ceil(k)
+    if f == c:
+        return key(N[int(k)])
+    d0 = key(N[int(f)]) * (c - k)
+    d1 = key(N[int(c)]) * (k - f)
+    return d0 + d1
+
+# end of http://code.activestate.com/recipes/511478/ }}}
+
+
+def get_cluster_information(clusters, classes):
+    """Get cluster info."""
+    # class distribution
+    class_dict = {}
+    for read_acc, class_id in classes.items():
+        if class_id not in class_dict:
+            class_dict[class_id] = [read_acc]
+        else:
+            class_dict[class_id].append(read_acc)
+
+    total_nr_classes = len(class_dict)
+    class_distribution = sorted([len(cl) for cl in class_dict.values()])
+    singleton_classes = set(
+        [acc_list[0] for cl_id, acc_list in class_dict.items()
+         if len(acc_list) == 1])
+    min_class_size = min(class_distribution)
+    max_class_size = max(class_distribution)
+    mean_class_size = sum(class_distribution) / float(len(class_distribution))
+    median_class_size = class_distribution[int(len(class_distribution) / 2)] \
+        if len(class_distribution) % 2 == 1 else (
+        class_distribution[int(len(class_distribution) / 2)] +
+        class_distribution[int(len(class_distribution) / 2) - 1]) / 2.0
+
+    upper_75_class_size = percentile(class_distribution, 0.75)
+    median_class_size = percentile(class_distribution, 0.5)
+
+    tot_size = sum(class_distribution)
+    e_class_size = sum(
+        [c_s**2 for c_s in class_distribution]) / float(tot_size)
+    tot_iterated_size = 0
+    for c_s in class_distribution[::-1]:
+        tot_iterated_size += c_s
+        if tot_iterated_size >= tot_size / 2.0:
+            n50_class_size = c_s
+            break
+
+    # cluster distribution
+    cluster_dict = {}
+    for read_acc, cl_id in clusters.items():
+        if cl_id not in cluster_dict:
+            cluster_dict[cl_id] = [read_acc]
+        else:
+            cluster_dict[cl_id].append(read_acc)
+
+    cluster_distribution = sorted([len(cl) for cl in cluster_dict.values()])
+
+    # in case unclustered reads are missing from output (as for isoseq3)
+    omitted_from_output_singletons = set(classes.keys()) - set(clusters.keys())
+    cluster_distribution = [
+        1 for i in range(
+            len(omitted_from_output_singletons))] + cluster_distribution
+    total_nr_clusters = len(cluster_distribution)
+
+    singleton_clusters = set(
+        [acc_list[0] for cl_id, acc_list in cluster_dict.items()
+         if len(acc_list) == 1])
+    min_cluster_size = min(cluster_distribution)
+    max_cluster_size = max(cluster_distribution)
+    mean_cluster_size = sum(cluster_distribution) / \
+        float(len(cluster_distribution))
+
+    upper_75_cluster_size = percentile(cluster_distribution, 0.75)
+    median_cluster_size = percentile(cluster_distribution, 0.5)
+
+    tot_size = sum(cluster_distribution)
+    e_cluster_size = sum(
+        [c_s**2 for c_s in cluster_distribution]) / float(tot_size)
+    tot_iterated_size = 0
+    for c_s in cluster_distribution[::-1]:
+        tot_iterated_size += c_s
+        if tot_iterated_size >= tot_size / 2.0:
+            n50_cluster_size = c_s
+            break
+
+    unaligned_but_nontrivially_clustered = set(
+        clusters.keys()) - singleton_clusters - set(classes.keys())
+
+    # not_considered = set([read for read in classes if read not in clusters ])
+
+    not_clustered_classes = defaultdict(int)
+    clustered_classes = defaultdict(int)
+    reads_not_clustered = defaultdict(list)
+    for read in clusters:
+        if read in classes:
+            class_id = classes[read]
+        else:
+            class_id = "unaligned"
+
+        if read in singleton_clusters:
+            not_clustered_classes[class_id] += 1
+            reads_not_clustered[class_id].append(read)
+        else:
+            clustered_classes[class_id] += 1
+
+    print("UNCLUSTERED:", "Tot classes:", len(not_clustered_classes))
+    print("CLUSTERED:", "Tot classes:", len(clustered_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)
+
+
+def main(args):
+    """Entry point."""
+    clusters = parse_inferred_clusters_tsv(args.clusters, args)
+    if not clusters:
+        outfile = open(args.outfile, "w")
+        outfile.write("No clusters created\n")
+        outfile.close()
+        return
+
+    if args.simulated:
+        ref_file = pysam.AlignmentFile(args.classes, "r", check_sq=False)
+        classes = parse_true_clusters_simulated(ref_file)
+        # by simulation we know classes of all reads, they are therefore the
+        # same number.
+        tot_nr_reads = len(classes)
+    else:
+        ref_file = pysam.AlignmentFile(args.classes, "rb", check_sq=False)
+        classes, tot_nr_reads, unclassified = parse_true_clusters(ref_file)
+
+    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)
+
+    outfile = open(args.outfile, "w")
+
+    outfile.write("CLASSES\n")
+
+    # reads, unaligned, classes, singleton, min,max, mean,median
+
+    outfile.write(
+        "{0},{1},{2},{3},{4},{5},{6},{7}\n".format(
+            "tot_nr_reads",
+            "unclassified",
+            "nr_non_singleton_classes",
+            "singleton_classes",
+            "upper_75_class_size",
+            "median_class_size",
+            "e_class_size",
+            "n50_class_size"))
+    outfile.write(
+        "{0},{1},{2},{3},{4},{5},{6},{7}\n".format(
+            tot_nr_reads,
+            unclassified,
+            nr_non_singleton_classes,
+            singleton_classes,
+            upper_75_class_size,
+            median_class_size,
+            e_class_size,
+            n50_class_size))
+
+    # Reads_nontrivially_clustered_(%), Singletons_(%),
+    # Reads_Nontrivially_clustered_but_unaligned, V, c,h ,V_nt, c_nt,h_nt,
+    # non_singleton_clusters, min, max, median, mean
+
+    Reads_nontrivially_clustered_percent = round(
+        100 * (float(tot_nr_reads - singleton_clusters) / tot_nr_reads), 1)
+    # round(1.0 - Reads_nontrivially_clustered_percent, 2)
+    Reads_Nontrivially_clustered_but_unaligned = \
+        unaligned_but_nontrivially_clustered
+    V, c, h = round(v_score, 3), round(compl_score, 3), round(homog_score, 3)
+
+    non_singleton_clusters = total_nr_clusters - singleton_clusters
+
+    print("NONTRIVIAL CLUSTERS: ", (total_nr_clusters - singleton_clusters))
+
+    outfile.write("CLUSTERS\n")
+    outfile.write(
+        "{0},{1},{2},{3},{4},{5},{6},{7},{8},{9},{10},{11}\n".format(
+            "V",
+            "c",
+            "h",
+            "ARI",
+            "Reads_nontrivially_clustered_percent",
+            "Reads_Nontrivially_clustered_but_unaligned",
+            "non_singleton_clusters",
+            "singleton_clusters",
+            "upper_75_cluster_size",
+            "median",
+            "e_cluster_size",
+            "n50_cluster_size"))
+    outfile.write(
+        "{0},{1},{2},{3},{4},{5},{6},{7},{8},{9},{10},{11}\n".format(
+            V,
+            c,
+            h,
+            ari,
+            Reads_nontrivially_clustered_percent,
+            Reads_Nontrivially_clustered_but_unaligned,
+            non_singleton_clusters,
+            singleton_clusters,
+            upper_75_cluster_size,
+            median_cluster_size,
+            e_cluster_size,
+            n50_cluster_size))
+    outfile.close()
+
+    if args.ctsv is not None:
+        cfh = open(args.ctsv, "w")
+        cfh.write("Read\tCluster\n")
+        for r, c in classes.items():
+            cfh.write("{}\t{}\n".format(r, c))
+        cfh.flush()
+        cfh.close()
+
+    dfc = pd.DataFrame(
+        {
+            'Statistic': [
+                'V-measure',
+                'ARI',
+                'Completeness',
+                'Homogeneity'],
+            'Value': [
+                V,
+                ari,
+                c,
+                h]}).set_index('Statistic')
+    dfn = pd.DataFrame(
+        {'Statistic': ['NonSingleton',
+                       'Singletons'],
+         'Value': [non_singleton_clusters,
+                   singleton_clusters]}).set_index('Statistic')
+    dfs = pd.DataFrame(
+        {
+            'Statistic': [
+                'Upper75ClsSize',
+                'Upper75ClassSize',
+                'MedianClsSize',
+                'MedianClassSize'],
+            'Value': [
+                upper_75_cluster_size,
+                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')
+
+    rdo = Path(args.raw_data_out)
+    dfc.to_csv(rdo / 'v_ari_com_hom.csv')
+    dfn.to_csv(rdo / 'sing_nonsing.csv')
+    dfs.to_csv(rdo / 'class_sizes1.csv')
+    dfs2.to_csv(rdo / 'class_sizes2.csv')
+
+    pages = PdfPages(args.report)
+
+    yd = 7 * 2.5
+
+    ax = dfc.plot(kind='bar', fontsize=7, rot=0)
+    for p in ax.patches:
+        ax.annotate(
+            "{:.2f}".format(
+                p.get_height()),
+            (p.get_x() +
+             p.get_width() /
+             2.0,
+             p.get_height()),
+            ha='center')
+    pages.savefig()
+    plt.clf()
+
+    ax = dfn.plot(kind='bar', fontsize=7, rot=0)
+    for p in ax.patches:
+        ax.annotate(
+            "{:.2f}".format(
+                p.get_height()),
+            (p.get_x() +
+             p.get_width() /
+             2.0,
+             p.get_height()),
+            ha='center')
+    pages.savefig()
+    plt.clf()
+
+    ax = dfs.plot(kind='bar', fontsize=7, rot=0)
+    for p in ax.patches:
+        ax.annotate(
+            "{:.2f}".format(
+                p.get_height()),
+            (p.get_x() +
+             p.get_width() /
+             2.0,
+             p.get_height() +
+             yd),
+            ha='center')
+    pages.savefig()
+    plt.clf()
+    ax = dfs2.plot(kind='bar', fontsize=7, rot=0)
+    for p in ax.patches:
+        ax.annotate(
+            "{:.2f}".format(
+                p.get_height()),
+            (p.get_x() +
+             p.get_width() /
+             2.0,
+             p.get_height() +
+             yd),
+            ha='center')
+    pages.savefig()
+    plt.clf()
+
+    pages.close()
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(
+        description="Align predicted transcripts to transcripts in ensembl "
+                    "reference data base.")
+    parser.add_argument(
+        '--clusters',
+        type=str,
+        help='Inferred clusters (tsv file)')
+    parser.add_argument(
+        '--classes',
+        type=str,
+        help='A sorted and indexed bam file.')
+    parser.add_argument(
+        '--ctsv',
+        default=None,
+        type=str,
+        help='Write true classes in this TSV file.')
+    parser.add_argument(
+        '--simulated',
+        action="store_true",
+        help='Simulated data, we can simply read correct classes '
+             'from the ref field.')
+    parser.add_argument(
+        '--ont',
+        action="store_true",
+        help='ONT data, parsing accessions differently.')
+    parser.add_argument(
+        '--modified_ont',
+        action="store_true",
+        help='ONT data preprocessed accessions, parsing '
+             'accessions differently.')
+    parser.add_argument('--outfile', type=str, help='Output file with results')
+    parser.add_argument(
+        '--report',
+        type=str,
+        help='Output PDF file with report')
+    parser.add_argument('--sizes', type=str, help='Cluster sizes')
+    parser.add_argument(
+        '--raw_data_out',
+        type=str,
+        help='dir to save raw data for plotting')
+    args = parser.parse_args()
+
+    print("------------------------------------------------------------")
+    main(args)
+    print("------------------------------------------------------------")

bin/plot_aln_stats.py

@@ -1,186 +0,0 @@
-#!/usr/bin/env python
-"""Plot a seqkit alignment stats file."""
-
-import argparse
-import warnings
-
-import matplotlib
-from matplotlib.backends.backend_pdf import PdfPages
-import matplotlib.pyplot as plt
-import numpy as np
-import pandas as pd
-import six
-
-matplotlib.use('Agg')
-
-with warnings.catch_warnings():
-    warnings.simplefilter("ignore")
-    import seaborn as sns
-
-warnings.resetwarnings()
-_ = sns
-
-# Parse command line arguments:
-parser = argparse.ArgumentParser(
-    description="""Plot a seqkit alignment stats file.""")
-parser.add_argument(
-    '-r', metavar='report_pdf', type=str, help="Report PDF (stats.pdf).",
-    default="stats.pdf")
-parser.add_argument(
-    'input', metavar='input_tsv', type=str, help="Input TSV.")
-
-
-class Report:
-    """Class for plotting utilities on the top of matplotlib.
-
-    Plots are saved in the specified file through the PDF backend.
-    """
-
-    def __init__(self, pdf):
-        """
-        Init Report with a matplotlib PdfPages instance.
-
-        :param self: object.
-        :param pdf: Output pdf.
-        :returns: The report object.
-        :rtype: Report
-        """
-        self.pdf = pdf
-        self.plt = plt
-        self.pages = PdfPages(pdf)
-
-    def _set_properties_and_close(self, fig, title, xlab, ylab):
-        """Set title, axis labels and close the figure.
-
-        :param self: object.
-        :param fig: The current figure.
-        :param title: Figure title.
-        :param xlab: X axis label.
-        :param ylab: Y axis label.
-        :returns: None
-        :rtype: object
-        """
-        plt.xlabel(xlab)
-        plt.ylabel(ylab)
-        plt.title(title)
-        self.pages.savefig(fig)
-        plt.close(fig)
-
-    def plot_boxplots(self, data_map, title="", xlab="", ylab="",
-                      xticks_rotation=0, xticks_fontsize=5):
-        """Plot multiple pairs of data arrays.
-
-        :param self: object.
-        :param data_map:
-            A dictionary with labels as keys and lists as data values.
-        :param title: Figure title.
-        :param xlab: X axis label.
-        :param ylab: Y axis label.
-        :param xticks_rotation: Rotation value for x tick labels.
-        :param xticks_fontsize: Fontsize for x tick labels.
-        :returns: None
-        :rtype: object
-        """
-        fig = plt.figure()
-        plt.boxplot(list(data_map.values()))
-        plt.xticks(np.arange(len(data_map)) + 1, data_map.keys(),
-                   rotation=xticks_rotation, fontsize=xticks_fontsize)
-        self._set_properties_and_close(fig, title, xlab, ylab)
-
-    def plot_bars_simple(self, data_map, title="", xlab="", ylab="", alpha=0.6,
-                         xticks_rotation=0, auto_limit=False):
-        """Plot simple bar chart from input dictionary.
-
-        :param self: object.
-        :param data_map: A dictionary with labels as keys and data as values.
-        :param title: Figure title.
-        :param xlab: X axis label.
-        :param ylab: Y axis label.
-        :param alpha: Alpha value.
-        :param xticks_rotation: Rotation value for x tick labels.
-        :param auto_limit: Set y axis limits automatically.
-        :returns: None
-        :rtype: object
-        """
-        fig = plt.figure()
-
-        labels = list(data_map.keys())
-        data = list(data_map.values())
-        positions = np.arange(len(labels))
-        plt.bar(positions, data, align='center', alpha=alpha)
-        plt.xticks(positions, labels, rotation=xticks_rotation)
-
-        if auto_limit:
-            low, high = min(data), max(data)
-            plt.ylim([(low - 0.5 * (high - low)), (high + 0.5 * (high - low))])
-
-        self._set_properties_and_close(fig, title, xlab, ylab)
-
-    def plot_histograms(self, data_map, title="", xlab="", ylab="", bins=50,
-                        alpha=0.7, legend_loc='best', legend=True,
-                        vlines=None):
-        """Plot histograms of multiple data arrays.
-
-        :param self: object.
-        :param data_map:
-            A dictionary with labels as keys and data arrays as values.
-        :param title: Figure title.
-        :param xlab: X axis label.
-        :param ylab: Y axis label.
-        :param bins: Number of bins.
-        :param alpha: Transparency value for histograms.
-        :param legend_loc: Location of legend.
-        :param legend: Plot legend if True.
-        :param vlines:
-            Dictionary with labels and positions of vertical lines to draw.
-        :returns: None
-        :rtype: object
-        """
-        fig = plt.figure()
-
-        for label, data in six.iteritems(data_map):
-            if len(data) > 0:
-                plt.hist(data, bins=bins, label=label, alpha=alpha)
-        if vlines is not None:
-            for label, pos in six.iteritems(vlines):
-                plt.axvline(x=pos, label=label)
-        if legend:
-            plt.legend(loc=legend_loc)
-
-        self._set_properties_and_close(fig, title, xlab, ylab)
-
-    def close(self):
-        """Close PDF backend.
-
-        Do not forget to call this at the end of your
-        script or your output will be damaged!
-
-        :param self: object
-        :returns: None
-        :rtype: object
-        """
-        self.pages.close()
-
-
-if __name__ == '__main__':
-    args = parser.parse_args()
-
-    plotter = Report(args.r)
-
-    # Plot overview panel:
-    stats = pd.read_csv(args.input, sep="\t")
-    perc = stats[["PrimAlnPerc", "MultimapPerc"]].copy()
-    num = stats[["PrimAln", "SecAln", "SupAln", "Unmapped", "TotalReads",
-                "TotalRecords"]].copy()
-
-    perc.plot(kind='bar')
-    plt.title("Percent primary and multimapping reads")
-    plt.tight_layout()
-    plotter.pages.savefig()
-
-    num.plot(kind='bar')
-    plt.title("Number of alignment records")
-    plt.tight_layout()
-    plotter.pages.savefig()
-
-    plotter.close()

bin/plot_gffcmp_stats.py

@@ -1,471 +0,0 @@
-#!/usr/bin/env python
-
-"""Plot a gffcompare stats file."""
-
-import argparse
-from collections import OrderedDict
-import warnings
-
-import matplotlib
-from matplotlib.backends.backend_pdf import PdfPages
-import matplotlib.pyplot as plt
-import numpy as np
-import pandas as pd
-import six
-
-
-matplotlib.use('Agg')
-
-with warnings.catch_warnings():
-    warnings.simplefilter("ignore")
-    import seaborn as sns
-
-warnings.resetwarnings()
-_ = sns
-
-# Parse command line arguments:
-parser = argparse.ArgumentParser(
-    description="""Plot a gffcompare stats file.""")
-parser.add_argument(
-    '-r', metavar='report_pdf', type=str,
-    help="Report PDF (plot_gffcmp_stats.pdf).",
-    default="plot_gffcmp_stats.pdf")
-parser.add_argument(
-    '-t', metavar='tracking_tsv', type=str,
-    help="Tracking file produced by gffcompare.", default=None)
-parser.add_argument(
-    'input', metavar='input_txt', type=str,
-    help="Input gffcompare stats file.")
-
-
-class Report:
-    """Matplotlib plotting utilities."""
-
-    def __init__(self, pdf):
-        """Init class with PdfPahges instance.
-
-        Plots are saved in the specified file through the PDF backend.
-
-        :param self: object.
-        :param pdf: Output pdf.
-        :returns: The report object.
-        :rtype: Report
-
-        """
-        self.pdf = pdf
-        self.plt = plt
-        self.pages = PdfPages(pdf)
-
-    def _set_properties_and_close(self, fig, title, xlab, ylab):
-        """Set title, axis labels and close the figure.
-
-        :param self: object.
-        :param fig: The current figure.
-        :param title: Figure title.
-        :param xlab: X axis label.
-        :param ylab: Y axis label.
-        :returns: None
-        :rtype: object
-        """
-        plt.xlabel(xlab)
-        plt.ylabel(ylab)
-        plt.title(title)
-        self.pages.savefig(fig)
-        plt.close(fig)
-
-    def plot_boxplots(self, data_map, title="", xlab="", ylab="",
-                      xticks_rotation=0, xticks_fontsize=5):
-        """Plot multiple pairs of data arrays.
-
-        :param self: object.
-        :param data_map: A dictionary with labels as keys and lists as data
-            values.
-        :param title: Figure title.
-        :param xlab: X axis label.
-        :param ylab: Y axis label.
-        :param xticks_rotation: Rotation value for x tick labels.
-        :param xticks_fontsize: Fontsize for x tick labels.
-        :returns: None
-        :rtype: object
-        """
-        fig = plt.figure()
-        plt.boxplot(list(data_map.values()))
-        plt.xticks(np.arange(len(data_map)) + 1, data_map.keys(),
-                   rotation=xticks_rotation, fontsize=xticks_fontsize)
-        self._set_properties_and_close(fig, title, xlab, ylab)
-
-    def plot_bars_simple(self, data_map, title="", xlab="", ylab="", alpha=0.6,
-                         xticks_rotation=0, auto_limit=False):
-        """Plot simple bar chart from input dictionary.
-
-        :param self: object.
-        :param data_map: A dictionary with labels as keys and data as values.
-        :param title: Figure title.
-        :param xlab: X axis label.
-        :param ylab: Y axis label.
-        :param alpha: Alpha value.
-        :param xticks_rotation: Rotation value for x tick labels.
-        :param auto_limit: Set y axis limits automatically.
-        :returns: None
-        :rtype: object
-        """
-        fig = plt.figure()
-
-        labels = list(data_map.keys())
-        data = list(data_map.values())
-        positions = np.arange(len(labels))
-        plt.bar(positions, data, align='center', alpha=alpha)
-        plt.xticks(positions, labels, rotation=xticks_rotation)
-
-        if auto_limit:
-            low, high = min(data), max(data)
-            plt.ylim([(low - 0.5 * (high - low)), (high + 0.5 * (high - low))])
-
-        self._set_properties_and_close(fig, title, xlab, ylab)
-
-    def plot_histograms(self, data_map, title="", xlab="", ylab="", bins=50,
-                        alpha=0.7, legend_loc='best', legend=True,
-                        vlines=None):
-        """Plot histograms of multiple data arrays.
-
-        :param self: object.
-        :param data_map: A dictionary with labels as keys and data arrays
-            as values.
-        :param title: Figure title.
-        :param xlab: X axis label.
-        :param ylab: Y axis label.
-        :param bins: Number of bins.
-        :param alpha: Transparency value for histograms.
-        :param legend_loc: Location of legend.
-        :param legend: Plot legend if True.
-        :param vlines: Dictionary with labels and positions of vertical lines
-            to draw.
-        :returns: None
-        :rtype: object
-        """
-        fig = plt.figure()
-
-        for label, data in six.iteritems(data_map):
-            if len(data) > 0:
-                plt.hist(data, bins=bins, label=label, alpha=alpha)
-        if vlines is not None:
-            for label, pos in six.iteritems(vlines):
-                plt.axvline(x=pos, label=label)
-        if legend:
-            plt.legend(loc=legend_loc)
-
-        self._set_properties_and_close(fig, title, xlab, ylab)
-
-    def close(self):
-        """Close PDF backend.
-
-        Do not forget to call this at the end of your
-        script or your output will be damaged!
-
-        :param self: object
-        :returns: None
-        :rtype: object
-        """
-        self.pages.close()
-
-
-def _parse_stat_line(sl):
-    """Parse a stats line."""
-    res = {}
-    tmp = sl.split(':')[1]
-    tmp = tmp.split('|')
-    res['sensitivity'] = float(tmp[0].strip())
-    res['precision'] = float(tmp[1].strip())
-    return res
-
-
-def _parse_matching_line(line):
-    """Parse a metching line."""
-    tmp = line.split(':')[1].strip()
-    return int(tmp)
-
-
-def _parse_mn_line(line):
-    """Parse a miss or novel line."""
-    res = {}
-    tmp = line.split(':')[1].strip()
-    tmp = tmp.split('/')
-    res['value'] = int(tmp[0])
-    tmp = tmp[1].split('(')
-    res['value_total'] = int(tmp[0].strip())
-    res['percent'] = float(tmp[1].split('%)')[0])
-    return res
-
-
-def _parse_total_line(line):
-    """Parse a total line."""
-    res = {}
-    tmp = line.split(':')[1].strip()
-    tmp = tmp.split('in')
-    res['transcripts'] = int(tmp[0].strip())
-    tmp = tmp[1].split('loci')
-    res['loci'] = int(tmp[0].strip())
-    tmp = int(tmp[1].split('(')[1].split(' ')[0])
-    res['me_transcripts'] = tmp
-    return res
-
-
-def parse_gffcmp_stats(txt):
-    """Parse a gffcompare stats file.
-
-    :param txt: Path to the gffcompare stats file.
-    :returns: Return as tuple of dataframes containing:
-        perfromance statistics, match statistics, miss statistics,
-        novel statistics, total statistics.
-    :rtype: tuple
-    """
-    sensitivity = []
-    precision = []
-    level = []
-
-    matching = OrderedDict()
-
-    missed_level = []
-    missed = []
-    missed_total = []
-    missed_percent = []
-
-    novel_level = []
-    novel = []
-    novel_total = []
-    novel_percent = []
-
-    total_target = []
-    total_loci = []
-    total_transcripts = []
-    total_multiexonic = []
-
-    fh = open(txt, 'r')
-    for line in fh:
-        line = line.strip()
-        if len(line) == 0:
-            continue
-        # Parse totals:
-        if line.startswith('#     Query mRNAs'):
-            total_target.append('Query')
-            r = _parse_total_line(line)
-            total_loci.append(r['loci'])
-            total_transcripts.append(r['transcripts'])
-            total_multiexonic.append(r['me_transcripts'])
-
-        if line.startswith('# Reference mRNAs '):
-            total_target.append('Reference')
-            r = _parse_total_line(line)
-            total_loci.append(r['loci'])
-            total_transcripts.append(r['transcripts'])
-            total_multiexonic.append(r['me_transcripts'])
-
-        # Parse basic statistics:
-        if line.startswith('Base level'):
-            st = _parse_stat_line(line)
-            level.append('Base')
-            sensitivity.append(st['sensitivity'])
-            precision.append(st['precision'])
-        if line.startswith('Exon level'):
-            st = _parse_stat_line(line)
-            level.append('Exon')
-            sensitivity.append(st['sensitivity'])
-            precision.append(st['precision'])
-        if line.startswith('Intron level'):
-            st = _parse_stat_line(line)
-            level.append('Intron')
-            sensitivity.append(st['sensitivity'])
-            precision.append(st['precision'])
-        if line.startswith('Intron chain level'):
-            st = _parse_stat_line(line)
-            level.append('Intron chain')
-            sensitivity.append(st['sensitivity'])
-            precision.append(st['precision'])
-        if line.startswith('Transcript level'):
-            st = _parse_stat_line(line)
-            level.append('Transcript')
-            sensitivity.append(st['sensitivity'])
-            precision.append(st['precision'])
-        if line.startswith('Locus level'):
-            st = _parse_stat_line(line)
-            level.append('Locus')
-            sensitivity.append(st['sensitivity'])
-            precision.append(st['precision'])
-
-        # Parse match statistics:
-        if line.startswith('Matching intron chains'):
-            m = _parse_matching_line(line)
-            matching['Intron chains'] = [m]
-        if line.startswith('Matching transcripts'):
-            m = _parse_matching_line(line)
-            matching['Transcripts'] = [m]
-        if line.startswith('Matching loci'):
-            m = _parse_matching_line(line)
-            matching['Loci'] = [m]
-
-        # Parse missing statistics:
-        if line.startswith('Missed exons'):
-            missed_level.append('Exons')
-            r = _parse_mn_line(line)
-            missed.append(r['value'])
-            missed_total.append(r['value_total'])
-            missed_percent.append(r['percent'])
-        if line.startswith('Missed introns'):
-            missed_level.append('Introns')
-            r = _parse_mn_line(line)
-            missed.append(r['value'])
-            missed_total.append(r['value_total'])
-            missed_percent.append(r['percent'])
-        if line.startswith('Missed loci'):
-            missed_level.append('Loci')
-            r = _parse_mn_line(line)
-            missed.append(r['value'])
-            missed_total.append(r['value_total'])
-            missed_percent.append(r['percent'])
-
-        # Parse novel statistics:
-        if line.startswith('Novel exons'):
-            novel_level.append('Exons')
-            r = _parse_mn_line(line)
-            novel.append(r['value'])
-            novel_total.append(r['value_total'])
-            novel_percent.append(r['percent'])
-        if line.startswith('Novel introns'):
-            novel_level.append('Introns')
-            r = _parse_mn_line(line)
-            novel.append(r['value'])
-            novel_total.append(r['value_total'])
-            novel_percent.append(r['percent'])
-        if line.startswith('Novel loci'):
-            novel_level.append('Loci')
-            r = _parse_mn_line(line)
-            novel.append(r['value'])
-            novel_total.append(r['value_total'])
-            novel_percent.append(r['percent'])
-
-    fh.close()
-
-    df_stats = pd.DataFrame(OrderedDict(
-        [('Sensitivity', sensitivity), ('Precision', precision)]), index=level)
-    df_match = pd.DataFrame(matching, index=['Matching'])
-
-    df_miss = pd.DataFrame(
-        OrderedDict(
-                [('Total', missed_total),
-                 ('Missed', missed),
-                 ('Percent missed', missed_percent)]), index=missed_level)
-
-    df_novel = pd.DataFrame(
-        OrderedDict(
-                [('Total', novel_total),
-                 ('Novel', novel),
-                 ('Percent novel', novel_percent)]), index=novel_level)
-
-    df_total = pd.DataFrame(OrderedDict(
-        [('Loci', total_loci), ('Transcripts', total_transcripts),
-         ('Multiexonic', total_multiexonic)]), index=total_target)
-
-    return df_stats, df_match, df_miss, df_novel, df_total
-
-
-if __name__ == '__main__':
-    args = parser.parse_args()
-
-    stats, match, miss, novel, total = parse_gffcmp_stats(args.input)
-    tracking = pd.read_csv(args.t, sep="\t", header=None, usecols=[0, 3],
-                           names=['Count', 'Overlaps'])
-
-    tracking = tracking.groupby("Overlaps").count().reset_index()
-    tracking = tracking.sort_values("Overlaps")
-
-    plotter = Report(args.r)
-
-    # Plot overview panel:
-
-    plt.figure(1)
-    plt.subplot(2, 2, 1)
-    total.plot(ax=plt.gca(), kind='barh', sharex=False, title='Totals')
-    plt.tight_layout()
-
-    plt.subplot(2, 2, 2)
-    stats.plot(ax=plt.gca(), kind='barh', legend=True, sharex=False,
-               title='Performance').legend(loc='best')
-    plt.tight_layout()
-
-    plt.subplot(2, 2, 3)
-
-    miss.copy().drop(
-        'Percent missed', axis=1).plot(
-            ax=plt.gca(), kind='barh',
-            legend=True, sharex=False,
-            title='Missed')
-
-    plt.tight_layout()
-
-    plt.subplot(2, 2, 4)
-
-    novel.copy().drop(
-        'Percent novel', axis=1).plot(
-            ax=plt.gca(), kind='barh', legend=True, sharex=False,
-            title='Novel')
-
-    plt.tight_layout()
-    plotter.pages.savefig()
-
-    # Plot individual panels:
-
-    total.plot(kind='barh', subplots=True, legend=False, sharex=False)
-    plt.tight_layout()
-    plotter.pages.savefig()
-
-    stats.plot(kind='barh', subplots=True, legend=False, sharex=False)
-    plt.tight_layout()
-    plotter.pages.savefig()
-
-    match.plot(kind='barh', subplots=True, legend=False)
-    plt.tight_layout()
-    plotter.pages.savefig()
-
-    miss.plot(kind='barh', subplots=True, legend=False, sharex=False)
-    plt.tight_layout()
-    plotter.pages.savefig()
-
-    novel.plot(kind='barh', subplots=True, legend=False, sharex=False)
-    plt.tight_layout()
-    plotter.pages.savefig()
-
-    def fix_names(s):
-        """Map trancript classification codes."""
-        names = {
-            '=': 'ExactMatch:=',
-            'c': 'Contained:c',
-            'k': 'ReverseContained:k',
-            'm': 'RetainedIntron:m',
-            'n': 'PartRetainedIntron:n',
-            'j': 'PartialMatch:j',
-            'e': 'TransFragMatch:e',
-            's': 'OppositeMatch:s',
-            'o': 'OtherSameStrand:o',
-            'x': 'ExonicOpposite:o',
-            'y': 'RefInIntrons:y',
-            'p': 'PolymeraseRunon:p',
-            'r': 'Repeat:r',
-            'u': 'Intergenic:u',
-            'i': 'FullyIntronic:i',
-            }
-        return names[s]
-
-    # Plot overlaps panel:
-    tracking.Overlaps = tracking.Overlaps.apply(fix_names)
-    tracking = tracking.set_index("Overlaps")
-    tracking.plot(kind='bar', title="Overlaps detected by gffcompare",
-                  colormap='Paired')
-    plt.tight_layout()
-    plotter.pages.savefig()
-    tracking["Percent"] = tracking.Count * 100 / tracking.Count.sum()
-    tracking[["Percent"]].plot(
-        kind='bar', title="Overlaps detected by gffcompare", colormap='Paired')
-    plt.tight_layout()
-    plotter.pages.savefig()
-
-    plotter.close()

bin/report.py

@@ -2,7 +2,10 @@
 """Create workflow report."""
 
 import argparse
-from collections import OrderedDict
+from collections import Counter, defaultdict, OrderedDict
+from functools import reduce
+import math
+from pathlib import Path
 
 from aplanat import bars, lines
 from aplanat.components import fastcat
@@ -10,12 +13,65 @@ from aplanat.components import simple as scomponents
 from aplanat.report import WFReport
 from aplanat.util import Colors
 from bokeh.layouts import gridplot
-from bokeh.models import ColumnDataSource
+from bokeh.models import ColumnDataSource, Panel, Tabs
+from bokeh.models.widgets import DataTable, TableColumn
 from bokeh.palettes import Category10_10
 from bokeh.plotting import figure
 from bokeh.transform import dodge
+import gffutils
+from jinja2 import Template
 import numpy as np
 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
+
+
+class Table:
+    """A table report component.
+
+    Adapted from aplanat
+    """
+
+    def __init__(self, template, data_frame, index, table_id, **kwargs):
+        """Initialize table component.
+
+        :param dataframe: dataframe to turn in to simple table.
+        """
+        template = Template(template)
+
+        for key, val in kwargs.items():
+            if isinstance(val, bool):
+                kwargs[key] = str(val).lower()
+
+        self.div = template.render(dataframe=data_frame.to_html(
+            table_id=table_id,
+            index=index),
+            table_id=table_id,
+            kwargs=kwargs)
 
 
 def simple_hbar(df, y, right, title="", color=Colors.cerulean,
@@ -25,14 +81,11 @@ def simple_hbar(df, y, right, title="", color=Colors.cerulean,
     :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.
-
-    Move to planat when it's working?
     """
     defaults = {
         'output_backend': 'webgl',
-        'plot_height': 300, 'plot_width': 600}
+        'height': 300, 'width': 600}
     defaults.update(fig_kwargs)
-
     p = figure(y_range=df[y], height=250, title=title,
                toolbar_location=None, tools="")
 
@@ -224,15 +277,15 @@ def parse_gffcmp_stats(txt):
 
     df_miss = pd.DataFrame(
         OrderedDict(
-                [('Total', missed_total),
+            [('Total', missed_total),
-                 ('Missed', missed),
+             ('Missed', missed),
-                 ('Percent missed', missed_percent)]), index=missed_level)
+             ('Percent missed', missed_percent)]), index=missed_level)
 
     df_novel = pd.DataFrame(
         OrderedDict(
-                [('Total', novel_total),
+            [('Total', novel_total),
-                 ('Novel', novel),
+             ('Novel', novel),
-                 ('Percent novel', novel_percent)]), index=novel_level)
+             ('Percent novel', novel_percent)]), index=novel_level)
 
     df_total = pd.DataFrame(OrderedDict(
         [('Loci', total_loci), ('Transcripts', total_transcripts),
@@ -241,7 +294,7 @@ def parse_gffcmp_stats(txt):
     return df_stats, df_match, df_miss, df_novel, df_total
 
 
-def grouped_bar(df, title=""):
+def grouped_bar(df, title="", tilted_xlabs=False):
     """Create grouped bar plot from pandas dataframe.
 
     :param pandas.DataFrame
@@ -268,14 +321,17 @@ def grouped_bar(df, title=""):
     # https://docs.bokeh.org/en/latest/docs/user_guide/categorical.html
     dodge_range = (-0.25, 0.25)
     current_dodge = dodge_range[0]
-    dodge_increment = abs(dodge_range[0] - dodge_range[1]) \
+    dodge_increment = abs(dodge_range[0] - dodge_range[1]) / \
-        / (len(df.columns) - 1)
+        (len(df.columns) - 1)
+
+    if tilted_xlabs:
+        p.xaxis.major_label_orientation = math.pi / 4
 
     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)]))
+                            for x in range(num_colors)]))
         colors = [item for sublist in colors for item in sublist]
 
         if col == 'x_groups':
@@ -285,7 +341,8 @@ def grouped_bar(df, title=""):
 
         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)
+               width=width, source=source, color=color,
+               legend_label=col)
         current_dodge += dodge_increment
 
     p.x_range.range_padding = 0.1
@@ -295,53 +352,22 @@ def grouped_bar(df, title=""):
     return p
 
 
-def workflow_plots(report, df_aln_stats_file,
+def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
-                   gff_cmp_stats_file, gff_cmp_tracking_file):
     """Create various sections and plots in a WfReport.
 
     :param report: aplanat WFReport
-    :param df_aln_stats_file: alignment stats. Output of `seqkit bam -s`
+    :param gffcompare_outdirs: List of output directories from run_gffcompare
-    :param gff_cmp_stats_file: gffcompare stats file
-    :param gff_cmp_tracking_file: gffcompare tracking file
     :return: None
+
+    TODO: split this into separate functions
     """
-    df_aln_stats = pd.read_csv(df_aln_stats_file, sep='\t')
+    # If any of the gffcompare dirs are empty, skip this section
-    df_aln_stats = df_aln_stats.select_dtypes([np.number]).dropna(axis=1)
+    if not all([any(Path(x).iterdir()) for x in gffcompare_outdirs]):
-
-    section = report.add_section()
-    section.markdown('''
-    ### Read mapping summary
-
-    Output of [seqkit](https://bioinf.shenwei.me/seqkit/) bam -s''')
-
-    section.table(df_aln_stats)
-
-    # Percentage primary and secondary mapping
-    df_perc = df_aln_stats[['PrimAlnPerc', 'MultimapPerc']]
-    bar_perc = bars.simple_bar(
-        df_perc.columns.values, df_perc.iloc[0].values,
-        title='% primary and multimapping reads', colors=Colors.cerulean)
-
-    # Counts of read mapping class
-    df_counts = df_aln_stats.drop(columns=['PrimAlnPerc', 'MultimapPerc'])
-    bar_counts = bars.simple_bar(
-        df_counts.columns.values, df_counts.iloc[0].values,
-        title='Number of alignment records', colors=Colors.cerulean
-    )
-
-    grid = gridplot([bar_perc, bar_counts], ncols=2,
-                    plot_width=400, plot_height=400)
-    section.plot(grid)
-
-    # If gffcompare has not been run, finish report here
-    if not gff_cmp_stats_file or not gff_cmp_tracking_file:
         return
 
-    stats, _, miss, novel, total = \
-        parse_gffcmp_stats(gff_cmp_stats_file)
-
     # Plot overview panel:
     section = report.add_section()
+    # TODO: update this based on current version
     section.markdown('''
     ### Annotation summary
 
@@ -365,42 +391,47 @@ def workflow_plots(report, df_aln_stats_file,
     Features present in the query transcripts, but absent in the reference
     ''')
 
-    bar_totals = grouped_bar(total, title="Totals")
+    for id_, dir_ in zip(sample_ids, gffcompare_outdirs):
-    bar_performance = grouped_bar(stats, title="Performance")
+        stats, _, miss, novel, total = \
-    bar_missed = grouped_bar(miss, title="Missed")
+            parse_gffcmp_stats(dir_ / 'str_merged.stats')
-    bar_novel = grouped_bar(novel, title="Novel")
 
-    grid = gridplot([bar_totals, bar_performance, bar_missed, bar_novel],
+        bar_totals = grouped_bar(total, title="Totals")
-                    ncols=2, plot_width=400, plot_height=400)
+        bar_performance = grouped_bar(stats, title="Performance",
-    section.plot(grid)
+                                      tilted_xlabs=True)
+        bar_missed = grouped_bar(miss, title="Missed")
+        bar_novel = grouped_bar(novel, title="Novel")
 
-    def fix_names(s):
+        grid = gridplot([bar_totals, bar_performance, bar_missed, bar_novel],
-        """Map trancript classification codes."""
+                        ncols=4, width=270, height=280)
-        names = {
+        section.markdown("""
-            '=': 'ExactMatch:=',
+        #### Sample_id: {}
-            'c': 'Contained:c',
+        """.format(id_))
-            'k': 'ReverseContained:k',
+        section.plot(grid)
-            'm': 'RetainedIntron:m',
+
-            'n': 'PartRetainedIntron:n',
+    names = {
-            'j': 'PartialMatch:j',
+        '=': 'ExactMatch:=',
-            'e': 'TransFragMatch:e',
+        'c': 'Contained:c',
-            's': 'OppositeMatch:s',
+        'k': 'ReverseContained:k',
-            'o': 'OtherSameStrand:o',
+        'm': 'RetainedIntron:m',
-            'x': 'ExonicOpposite:o',
+        'n': 'PartRetainedIntron:n',
-            'y': 'RefInIntrons:y',
+        'j': 'PartialMatch:j',
-            'p': 'PolymeraseRunon:p',
+        'e': 'TransFragMatch:e',
-            'r': 'Repeat:r',
+        's': 'OppositeMatch:s',
-            'u': 'Intergenic:u',
+        'o': 'OtherSameStrand:o',
-            'i': 'FullyIntronic:i',
+        'x': 'ExonicOpposite:o',
-            }
+        'y': 'RefInIntrons:y',
-        return names[s]
+        'p': 'PolymeraseRunon:p',
+        'r': 'Repeat:r',
+        'u': 'Intergenic:u',
+        'i': 'FullyIntronic:i',
+    }
 
     # Plot overlaps panel:
     section = report.add_section()
     section.markdown('''
-    ## Query transfrag class assignments
+    ### Query transfrag class assignments
 
-    The classes that are assinged by
+    The classes that are assigned by
     [gffcompare](https://ccb.jhu.edu/software/stringtie/gffcompare.shtml),
     which describe the relationship between query transfrag and the most
     similar reference transcript.
@@ -409,112 +440,406 @@ def workflow_plots(report, df_aln_stats_file,
     gffcompare_codes.png) illustrates the different classes.
         ''')
 
-    tracking = pd.read_csv(gff_cmp_tracking_file, sep="\t", header=None,
+    tracking_dfs = []
-                           usecols=[0, 3], names=['Count', 'Overlaps'])
 
-    tracking = tracking.groupby("Overlaps").count().reset_index()
+    print(gffcompare_outdirs)
-    tracking = tracking.sort_values("Overlaps")
+    track_files = [x / 'str_merged.tracking' for x in gffcompare_outdirs]
-    tracking.Overlaps = tracking.Overlaps.apply(fix_names)
-    tracking["Percent"] = tracking.Count * 100 / tracking.Count.sum()
 
-    tracking_bar = simple_hbar(
+    df_tracking = load_sample_data(track_files, sample_ids,
-        tracking, 'Overlaps', 'Count', title="totals")
+                                   read_func=lambda x:
-    tracking_bar_perc = simple_hbar(
+                                   pd.read_csv(x, sep="\t", header=None,
-        tracking, 'Overlaps', 'Percent', title='percent'
+                                               usecols=[0, 3],
-    )
+                                               names=['Count', 'Overlaps']))
 
-    grid = gridplot([tracking_bar, tracking_bar_perc], ncols=2,
+    tabs = []
-                    plot_width=400, plot_height=400)
+    for id_, df_track in df_tracking.groupby('sample_id'):
-    section.plot(grid)
+        tracking = df_track.groupby("Overlaps").count().reset_index()
+        tracking.Overlaps = tracking.Overlaps.map(names)
+        tracking['Percent'] = tracking.Count * 100 / tracking.Count.sum()
+        tracking = tracking.sort_values("Overlaps")
+        track_bar = _hbar(
+            tracking['Overlaps'].values.tolist(),
+            tracking['Percent'].values.tolist(), title="{}".format(id_))
+
+        # Edit for creating unified table
+        tracking.rename(columns={'sample_id': id_ + ' count'}, inplace=True)
+        tracking.drop(columns=['Count'], inplace=True)
+        tracking_dfs.append(tracking)
+
+        # section.plot(grid)
+
+        # Tracking table
+        df_class_table = reduce(
+            lambda left, right: pd.merge(left, right), tracking_dfs)
+
+        desc = pd.Series(df_class_table.Overlaps.apply(
+            lambda x: x.split(':')[0]))
+        df_class_table.insert(0, 'description', desc)
+
+        code = pd.Series(df_class_table.Overlaps.apply(
+            lambda x: x.split(':')[1]))
+        df_class_table.insert(0, 'code', code)
+
+        cols = [TableColumn(field=Ci, title=Ci, width=100)
+                for Ci in df_class_table.columns]
+        track_table = DataTable(columns=cols,
+                                source=ColumnDataSource(df_class_table),
+                                index_position=None,
+                                width=500)
+
+        df_class_table.drop(columns=['Overlaps'], inplace=True)
+        tabs.append(Panel(
+            child=gridplot([track_bar, track_table], ncols=2), title=id_)
+        )
+
+    cover_panel = Tabs(tabs=tabs)
+    section.plot(cover_panel)
+
+    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)]
+        else:
+            bins = np.linspace(10, max_, 15)
+        bins = np.unique(bins)  # Low max_ can end up with duplicated bins
+
+        groups = pd.cut(df_code.TPM, bins).value_counts()
+        df_code.to_csv('dfcode.csv')
+        df_temp = pd.DataFrame.from_dict(dict(
+            x=[x.mid for x in groups.index], y=groups.values
+        ))
+        df_temp.sort_values(by='x', inplace=True)
+
+        x = [str(math.ceil(x)) for x in df_temp.x]
+        y = df_temp.y
+
+        reads_per_iso_geom = \
+            bars.simple_bar(x, y,
+                            title="{} - Num isoforms/TPM bin - "
+                                  "gffcompare class code - '{}'".format(
+                                      sample_id, class_code),
+                            colors=Colors.cerulean,
+                            x_axis_label='TPM',
+                            y_axis_label='Number of isoforms')
+
+        reads_per_iso_geom.xaxis.major_label_orientation = math.pi / 2.8
+        return reads_per_iso_geom
+
+    log_plots = defaultdict(list)
+
+    try:
+        tmap_files = [next(x.glob('*.tmap')) for x in gffcompare_outdirs]
+    except StopIteration:
+        print("Cannot find .tmap files in {}".format(gffcompare_outdirs))
+        return
+
+    df_tmap = load_sample_data(tmap_files, sample_ids)
+
+    for id_, df in df_tmap.groupby('sample_id'):
+
+        log_plots[id_].append(plot_isoforms_per_tpm_bin(
+            df, 'all', id_, geomspace=True))
+
+        for class_code, df_code in df.groupby('class_code'):
+            log_plots[id_].append(plot_isoforms_per_tpm_bin(
+                df_code, class_code, id_, geomspace=True))
+
+    tabs = []
+    for id_, sample_plots in log_plots.items():
+        tabs.append(Panel(
+            child=gridplot(sample_plots, ncols=2), title=id_))
+
+    section.markdown('''
+    ### Read coverage by gffcompare transfrag class''')
+    cover_panel = Tabs(tabs=tabs)
+    section.plot(cover_panel)
+
+    return df_tmap
 
 
 def pychopper_plots(report, df):
     """Make plots from pychopper.cdna_classifier.py.
 
     :param report: aplanat WFReport
-    :param df: result DataFrame
+    :param df: DataFrame of Pychopper stats
     """
     section = report.add_section()
     section.markdown('''
-    ### pychopper summary statisitcs
+    ### Pychopper summary statisitcs
 
     The following plots summarize the output of [cdna_classifier.py]
     (https://github.com/nanoporetech/pychopper)
 
-    * **Classification of output reads**:
+    * **Pr.found**: Reads with primers found in correct orientation at
-        * Primers_found: Reads with primers found in correct orientation at
+     both ends.
-        both ends.
+    * **Resc**: Reads 'rescued' from fused reads
-        * Rescue: Reads 'rescued' from fused reads
+    * **Unusable**: Read with missing or incorrect primer orientation
-        * Unusable: Read with missing or incorrect primer orientation
+    * **+/-**: Orientation of reads relative to the mRNA
-    * **Strand of oriented reads**:
-        * Strand of read relative to the mRNA
-    * **Strand of rescued read**:
-        * Strand of read that were rescued from fused reads
-    * **Number of primer alignment hits in unclassified reads**:
-        * Note: Need to look into what this means
-    * **Number of primer alignment hits in rescued reads**:
-        * Note: Need to look into what this means
-    * **Number of usable segments per rescued read**:
-        * Number of usable segments (primer-flanked, correctly oriented
-        regions) per fused read.
-    * **Usable bases as a function of cutoff**:
-        * The cutoff value supplied to the primer alignment tool.
-        Note: What are usabel bases in this conext
-    * ** Log10 length distribution of trimmed away sequences**:
-        * todo
     ''')
 
-    def g(df, index, title):
+    plots = []
-        df_ = df[df.index == index]
+    for id_, df in df.groupby('sample_id'):
-        groups = df_.Name.values
+        df1 = df.set_index('Name', drop=True)
-        bar_ = bars.simple_bar(
+        df1 = df1.T[['Primers_found', 'Rescue', 'Unusable']]
-            groups, df_['Value'].values,
+        df1.rename(columns={'Primers_found': 'Pr.found',
-            title=title, colors=Colors.cerulean)
+                            'Rescue': 'Resc',
-        return bar_
+                            'Unsable': 'Un'}, inplace=True)
+        df2 = df[df.index == 'Strand']
+        bar_chop = bars.simple_bar(
+            df1.columns.values.tolist() + df2.Name.values.tolist(),
+            df1.iloc[0].values.tolist() + df2.Value.values.tolist(),
+            title='{} - Pychopper stats'.format(id_),
+            colors=Colors.cerulean)
 
-    df1 = df.set_index('Name', drop=True)
+        plots.extend([bar_chop])
-    df1 = df1.T[['Primers_found', 'Rescue', 'Unusable']]
+    grid = gridplot(plots, ncols=4,
-    bar_class = bars.simple_bar(
+                    width=300, height=300)
-            df1.columns.values, df1.iloc[0].values,
-            title='Classification of output reads', colors=Colors.cerulean)
-
-    plots = [
-        bar_class,
-        g(df, 'Strand', 'Strand of oriented read'),
-        g(df, 'RescueStrand', 'Strand of rescued reads'),
-        g(df, 'UnclassHitNr', 'Number of hits in unclassified reads'),
-        g(df, 'RescueHitNr', 'Number of hits in rescued reads'),
-        g(df, 'RescueSegmentNr', 'Number of usable segments per rescued read')
-        ]
-
-    q = round(df.loc['Parameter', 'Value'], 4)
-    df_at = df[df.index == 'AutotuneSample'].astype('float')
-
-    # Add vertical line at x=q
-    ymin, ymax = df_at['Value'].min(), df_at['Value'].max()
-    plots.append(lines.line([df_at['Name'].values.tolist(), [q, q]],
-                            [df_at['Value'].values.tolist(), [ymin, ymax]],
-                            title=("Usable bases as function of cutoff(q).Best"
-                                   " q={}").format(q), colors=['blue', 'red']
-                            ))
-    df_unusable = df[df.index == 'Unusable'].astype('float')
-
-    plots.append(lines.line([np.log10(1 + df_unusable['Name'])],
-                            [df_unusable['Value']],
-                            title=("Log10 length distribution of trimmed away"
-                                   " sequences.")
-                            ))
-
-    grid = gridplot(plots, ncols=2,
-                    plot_width=400, plot_height=400)
     section.plot(grid)
 
 
+def cluster_quality(cluster_qc_dir, report, sample_ids):
+    """Make cluster quality section."""
+    section = report.add_section()
+    section.markdown('''
+    ### De novo clustering quality
+
+    This section shows plots relating to the clustering quality performed
+    by isONclust2. The full length reads are mapped to a reference genome
+    to create a ground truth of reads mapped to clusters. This is then compared
+    to the de novo-generated clusters, and the following statistics are
+    generated.
+
+    * [Homogeneity](https://scikit-learn.org/stable/modules/generated/
+    sklearn.metrics.homogeneity_score.html): Penalises over-clustering.
+
+    * [Completeness](https://scikit-learn.org/stable/modules/generated/
+    sklearn.metrics.completeness_score.html): Penalises under-clustering.
+
+    * [V-measure](https://clusteringjl.readthedocs.io/en/latest/vmeasure.html):
+    The harmonic mean of the homogeneity and completeness
+
+    * [Adjusted Rand Index](https://scikit-learn.org/stable/modules/generated/
+    sklearn.metrics.adjusted_rand_score.html): Intuitively,  measures the
+    percentage of read pairs correctly clustered, normalized so that a perfect
+    clustering = 1 and a random cluster assignment achieves = 0
+
+    * NonSingleton: Clusters with multiple reads
+    * Singleton: Clusters consisting of a single read (These do not contribute
+    to the final transcript calling - I need to check this!)
+
+    ''')
+
+    tabs = []
+    for id_, cluster_dir in zip(sample_ids, cluster_qc_dir):
+        plots = []
+        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)
+            bar = bars.simple_bar(
+                df.Statistic.values.tolist(), df.Value.values.tolist(),
+                colors=Colors.cerulean
+            )
+            bar.xaxis.major_label_orientation = math.pi / 2.8
+            plots.append(bar)
+        tabs.append(Panel(
+            child=gridplot(plots, ncols=4,
+                    width=300, height=300), title=id_))
+
+    cover_panel = Tabs(tabs=tabs)
+    section.plot(cover_panel)
+
+
+def transcript_table(report, df_tmaps, covr_threshold, table_template):
+    """Create searchable table of transcripts."""
+    section = report.add_section()
+
+    # Should we put data from each sample into it's own table or have it
+    # 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.sort_values('cov', ascending=True, inplace=True)
+    counts = list(range(len(df)))
+
+    # Keep Isoforms with coverage > threshold
+    vline_x = np.argmax(df['cov'] > covr_threshold)
+    vline_y = [0, df['cov'].max()]
+
+    cov_plt = lines.line(
+        [counts, [vline_x, vline_x]],  # x-values
+        [df['cov'].values.tolist(), vline_y],  # y-values
+        title=(
+            "Read Coverage. Threshold = {}x coverage".format(
+                covr_threshold)
+        ), x_axis_label='Num Isoforms',
+        y_axis_label='Coverage',
+        colors=['blue', 'red'])
+
+    section.markdown('''
+    ### Query transcript table
+
+    Low coverage transcripts are removed to speed up the table viewing. <br>
+    This can be set with the parameter `args.min_isoform_cov` in the config.
+      ''')
+    section.plot(cov_plt)
+    # Filter on converge threshold
+    df = df[df['cov'] >= covr_threshold]
+
+    # Make a column of number of isoforms in parent gene
+    gb = df.groupby(['ref_gene_id', 'sample_id']).count()
+    # gb = gb.set_index(['ref_gene_id', 'sample_id'])
+    gb.rename(columns={'ref_id': 'num_isoforms'}, inplace=True)
+
+    df['parent gene iso num'] = df.apply(
+        lambda x: gb.loc[(x.ref_gene_id, x.sample_id), 'num_isoforms'], axis=1)
+
+    df.sort_values('parent gene iso num', inplace=True, ascending=True)
+
+    with open(table_template, 'r') as fh:
+        tabletempl = fh.read()
+
+    bigtable = Table(tabletempl, df, index=False, table_id='bigtable')
+    section._add_item(bigtable.div)
+
+
+def tanscriptome_summary(report, gffs, sample_ids, denovo=False):
+    """
+    Plot transcriptome summaries.
+
+    Some of this data is available via gffcompare output, but the de novo
+    pipeline skips that, so we do it al 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()
+    section.markdown('''
+    ### Transcriptome summary
+    ''')
+
+    tabs = []
+    for id_, gff in zip(sample_ids, gffs):
+
+        plots = []
+
+        db = gffutils.create_db(
+            gff, dbfn=':memory:', force=True, keep_order=True,
+            merge_strategy='merge', sort_attribute_values=True
+        )
+
+        num_transcripts = db.count_features_of_type('transcript')
+        num_genes = db.count_features_of_type('gene')
+
+        transcript_lens = []
+        exons_per_transcript = Counter()
+        isoforms_per_gene = []
+
+        for g in db.features_of_type('gene'):
+
+            n_isos = len(list(db.children(g, featuretype='transcript')))
+            isoforms_per_gene.append(n_isos)
+
+            for t in db.children(
+                    g, featuretype='transcript', order_by='start'):
+                tr_len = 0
+                exons = list(enumerate(db.children(t, featuretype='exon')))
+                if len(exons) == 0:
+                    continue
+                for nx, ex in exons:
+                    tr_len += abs(ex.end - ex.start)
+
+                exons_per_transcript[nx] += 1
+
+                transcript_lens.append(tr_len)
+
+        num_bins = max(isoforms_per_gene)
+        if num_bins > 12:
+            num_bins = 12
+
+        hist_, bins = np.histogram(isoforms_per_gene, bins=num_bins)
+
+        bar_isos = bars.simple_bar([str(round(x_, 1)) for x_ in bins], hist_,
+                                   title="Isoforms per gene",
+                                   colors=Colors.cerulean,
+                                   x_axis_label='Num. isoforms',
+                                   y_axis_label='Num. genes')
+
+        bar_isos.xaxis.major_label_orientation = math.pi / 2.8
+        plots.append(bar_isos)
+
+        box = bars.boxplot_series(
+            [id_] * len(transcript_lens), transcript_lens,
+            width=70, ylim=(min(transcript_lens), max(transcript_lens)),
+            title='transcript lengths')
+        plots.append(box)
+
+        if not denovo:
+            x, y = zip(*sorted(exons_per_transcript.items()))
+
+            ept_bar = _vbar(x, list(y),
+                            title="Exons per transcript",
+                            color=Colors.cerulean)
+
+            ept_bar.xaxis.major_label_orientation = math.pi / 2.8
+            plots.append(ept_bar)
+
+        df_sum = pd.DataFrame.from_dict(
+            {'Total genes': [num_genes],
+             'Total transcripts': [num_transcripts],
+             'Max trans. len': max(transcript_lens),
+             'Min trans. len': min(transcript_lens)}).T
+        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)
+        plots.append(data_table)
+
+        tabs.append(Panel(
+            child=gridplot(plots, ncols=4,
+                           width=300, height=300), title=id_))
+
+    cover_panel = Tabs(tabs=tabs)
+    section.plot(cover_panel)
+
+
+def load_sample_data(files, sample_ids, read_func=None):
+    """Load CSVs into dataframe, and assign sample_id column."""
+    df_ = pd.DataFrame()
+    if not files:
+        return None
+    for id_, x in zip(sample_ids, files):
+        if read_func:
+            d = read_func(x)
+        else:
+            d = pd.read_csv(x, sep='\t+')
+        d['sample_id'] = id_
+        df_ = pd.concat([df_, d])
+    return df_
+
+
 def main():
     """Run the entry point."""
     parser = argparse.ArgumentParser()
-    parser.add_argument("report", help="Report output file")
+    parser.add_argument("--report", help="Report output file")
-    parser.add_argument("summaries", nargs='+', help="Read summary file.")
+    parser.add_argument("--summaries", nargs='+', help="Read summary file.")
     parser.add_argument(
         "--versions", required=True,
         help="directory containing CSVs containing name,version.")
@@ -529,36 +854,95 @@ def main():
         help="git commit of the executed workflow")
 
     parser.add_argument(
-        "--alignment_stats", required=True,
+        "--alignment_stats", required=False, default=None, nargs='*',
-        help="TSV summary file of alignment statistics")
-
-    parser.add_argument(
-        "--gffcompare_tracking", required=False, default=None,
         help="TSV summary file of alignment statistics")
     parser.add_argument(
-        "--gffcompare_stats", required=False, default=None,
+        "--gff_annotation", required=True, nargs='+',
-        help="TSV summary file of alignment statistics")
+        help="transcriptome annotation gff file")
-
     parser.add_argument(
-        "--pychop_report", required=True,
+        "--gffcompare_dir", required=False, default=None, nargs='*',
+        help="gffcompare outout dir")
+    parser.add_argument(
+        "--pychop_report", required=True, nargs='+',
         help="TSV summary file of pychopper statistics")
+    parser.add_argument(
+        "--sample_ids", required=True, nargs='+',
+        help="List of sample ids")
+    parser.add_argument(
+        "--report_template", required=True,
+        help="Jinja template")
+    parser.add_argument(
+        "--table_template", required=True,
+        help="Template for big transcript table")
+    parser.add_argument(
+        "--transcript_table_cov_thresh", required=False, type=int, default=50,
+        help="Isoforms without this support will be excluded from the table")
+    parser.add_argument(
+        "--cluster_qc_dirs", required=False, type=str, default=None, nargs='*',
+        help="Directory with various cluster quality csvs")
+    parser.add_argument('--denovo', dest='denovo', action='store_true')
+
     args = parser.parse_args()
+    print('denovo', args.denovo)
+
+    sample_ids = args.sample_ids
 
     report = WFReport(
-        "Workflow for assembling transcript isoforms", "wf-isoforms",
+        "Transcript isoform report", "wf-isoforms",
         revision=args.revision, commit=args.commit)
 
     # Add reads summary section
-    report.add_section(
+    for id_, summ in zip(sample_ids, args.summaries):
-        section=fastcat.full_report(args.summaries))
+        report.add_section(
+            section=fastcat.full_report(
+                [summ],
+                header='#### Read stats: {}'.format(id_)
+            ))
+
+    if args.alignment_stats is not None:
+        df_aln_stats = load_sample_data(args.alignment_stats, sample_ids)
+        section = report.add_section()
+        section.markdown('''
+          ### Read mapping summary
+
+          Output of [seqkit](https://bioinf.shenwei.me/seqkit/) bam -s''')
+
+        section.table(df_aln_stats)
 
     # workflow-specific plotting
-    workflow_plots(report, args.alignment_stats,  args.gffcompare_stats,
+    tanscriptome_summary(report, args.gff_annotation, sample_ids,
-                   args.gffcompare_tracking)
+                         denovo=args.denovo)
 
-    if args.pychop_report:
+    df_tmaps = gff_compare_plots(
-        df_chop_stats = pd.read_csv(args.pychop_report, sep='\t', index_col=0)
+        report,
-        pychopper_plots(report, df_chop_stats)
+        [Path(x) for x in args.gffcompare_dir],
+        sample_ids)
+
+    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)
+
+    with open(args.report_template, "r") as fh:
+        reptempl = fh.read()
+
+    report.template = Template(reptempl)
+
+    if df_tmaps is not None:
+        transcript_table(report, df_tmaps, args.transcript_table_cov_thresh,
+                         args.table_template)
+
+    if args.cluster_qc_dirs is not None:
+        cluster_quality(args.cluster_qc_dirs, report, sample_ids)
 
     # Arguments and software versions
     report.add_section(
@@ -566,7 +950,6 @@ def main():
     report.add_section(
         section=scomponents.params_table(args.params))
 
-    # write report
     report.write(args.report)
 
 

bin/report_template.html

@@ -0,0 +1,42 @@
+<!doctype html>
+<html lang="en" xmlns="http://www.w3.org/1999/html">
+<head>
+
+    <meta charset="utf-8">
+    <meta http-equiv="X-UA-Compatible" content="IE=edge">
+    <meta name="viewport" content="width=device-width, initial-scale=1">
+    <title>wf-isoforms report</title>
+    <link rel="stylesheet"
+          href="https://stackpath.bootstrapcdn.com/bootstrap/3.4.1/css/bootstrap.min.css"
+          integrity="sha384-HSMxcRTRxnN+Bdg0JdbxYKrThecOKuH5zCYotlSAcp1+c8xmyTe9GYg1l9a69psu"
+          crossorigin="anonymous">
+    <link href="https://cdn.jsdelivr.net/npm/simple-datatables@latest/dist/style.css"
+          rel="stylesheet" type="text/css">
+    <script src="https://cdn.jsdelivr.net/npm/simple-datatables@latest"
+            type="text/javascript"></script>
+<link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/4.7.0/css/font-awesome.min.css" rel="stylesheet"/>
+    <script
+            src="https://code.jquery.com/jquery-3.6.0.slim.min.js"
+            integrity="sha256-u7e5khyithlIdTpu22PHhENmPcRdFiHRjhAuHcs05RI="
+            crossorigin="anonymous"></script>
+    <script src="https://cdn.datatables.net/1.11.3/js/jquery.dataTables.min.js"
+            type="text/javascript"></script>
+    <link rel="stylesheet"
+          href="https://cdn.datatables.net/1.11.3/css/jquery.dataTables.min.css"
+          type="text/css">
+
+    {{ resources }}
+    {{ script }}
+<!--    delete?-->
+    {{ bigtable_js }}
+</head>
+<body>
+<div class="container">
+    <h1>{{ title }}</h1>
+    <p class="lead">{{ lead }}
+        {{ div }}
+    <div class="container" id="bigtable">
+        {{ big_table }}
+    </div>
+</body>
+</html>

bin/run_fastq_qc.py

@@ -1,59 +0,0 @@
-#!/usr/bin/env python
-"""Get fastq QC reports."""
-
-# -*- coding: utf-8 -*-
-
-import argparse
-import os
-import sys
-
-import numpy as np
-from pysam import FastxFile
-
-
-def parse_args(argv=sys.argv[1:]):
-    """Parse args."""
-    description = """Script to run the isoform workflow """
-    parser = argparse.ArgumentParser(description=description)
-    parser.add_argument("--fastq", required=True, help="")
-    parser.add_argument("--output_dir", required=True, help="")
-    return parser.parse_args(argv)
-
-
-def run_fastq_qc(fastq_path, output):
-    """Write QC info to files."""
-    qualities = list()
-    mean_qualities = list()
-    lengths = list()
-
-    with FastxFile(fastq_path) as fq:
-        for rec in fq:
-            # ONT calculation for "mean Q score"
-            quals = np.fromiter(
-                (ord(x) - 33 for x in rec.quality),
-                dtype=int, count=len(rec.quality))
-            mean_p = np.mean(np.power(10, quals / -10))
-            mean_qualities.append(-10 * np.log10(mean_p))
-            # all qualities
-            qualities.extend(quals)
-            lengths.append(len(quals))
-
-    with open(os.path.join(output, "base_qual.txt"), 'w') as f:
-        f.write("\n".join((str(q) for q in qualities)))
-
-    with open(os.path.join(output, "read_qual.txt"), 'w') as f:
-        f.write("\n".join((str(q) for q in mean_qualities)))
-
-    with open(os.path.join(output, "lengths.txt"), 'w') as f:
-        f.write("\n".join((str(_l) for _l in lengths)))
-
-
-def main(args):
-    """Run entry point."""
-    assert os.path.isfile(args.fastq)
-    assert os.path.isdir(args.output_dir)
-    run_fastq_qc(fastq_path=args.fastq, output=args.output_dir)
-
-
-if __name__ == '__main__':
-    main(args=parse_args())

bin/run_isonclust2.py

@@ -0,0 +1,124 @@
+#!/usr/bin/env python
+"""Dynamically generate isONclust2 processes."""
+from collections import OrderedDict
+from glob import glob
+from itertools import zip_longest
+from pathlib import Path
+import re
+import subprocess as sub
+
+
+class Node:
+    """Node."""
+
+    def __init__(self, Id, File, Left, Right, Parent, Level):
+        """Set node attaributes."""
+        self.Id = Id
+        self.File = File
+        self.Left = Left
+        self.Right = Right
+        self.Parent = Parent
+        self.Level = Level
+        self.Done = False
+        self.RightSide = False
+
+    def __repr__(self):
+        """Get string repr of a node."""
+        return "Node:{} Level: {} File: {} Done: {} Left: {} Right: " \
+               "{} Parent: {}".format(
+                self.Id, self.Level,
+                self.File, self.Done, self.Left.Id if
+                self.Left is not None else None,
+                self.Right.Id if self.Right is not None else None,
+                self.Parent.Id if self.Parent is not None else None)
+
+
+def grouper(n, iterable, fillvalue=None):
+    """
+    Group adjacent nodes.
+
+    grouper(3, 'ABCDEFG', 'x') --> ABC DEF Gxx.
+    """
+    args = [iter(iterable)] * n
+    return zip_longest(fillvalue=fillvalue, *args)
+
+
+def build_job_tree():
+    """Build a job tree of nodes."""
+    JOB_TREE = OrderedDict()
+    batches = glob("batches/isONbatch_*.cer")
+    batch_ids = [int(re.search('batches/isONbatch_(.*)\\.cer$', x).group(1))
+                 for x in batches]
+    LEVELS = OrderedDict()
+    LEVELS[0] = []
+    for Id, bf in sorted(zip(batch_ids, batches), key=lambda x: x[0]):
+        n = Node(
+            Id,
+            "clusters/isONcluster_{}.cer".format(Id),
+            None,
+            None,
+            None,
+            0)
+        n.Done = True
+        JOB_TREE[Id] = n
+        LEVELS[0].append(n)
+    level = 0
+    max_id = LEVELS[0][-1].Id
+    while len(LEVELS[level]) != 1:  # Final level will be link
+        next_level = level + 1
+        LEVELS[next_level] = []
+        for l_, r in grouper(2, LEVELS[level]):
+            if r is None:  # End of a level
+                LEVELS[level].pop()  # remove last node?
+                l_.Level += 1  # ncrement level
+                LEVELS[next_level].append(l_)  # Add the left to the next level
+                continue
+            max_id += 1
+            new_batch = "clusters/isONcluster_{}.cer".format(max_id)
+            new_node = Node(max_id, new_batch, l_, r, None, next_level)
+            l_.Parent = new_node
+            r.Parent = new_node
+            r.RightSide = True
+            LEVELS[next_level].append(new_node)
+            JOB_TREE[max_id] = new_node
+        level = next_level
+    ROOT = JOB_TREE[len(JOB_TREE) - 1].Id
+    JOB_TREE[ROOT].RightSide = True
+
+    return JOB_TREE, LEVELS
+
+
+def main():
+    """Entry point."""
+    Path('clusters').mkdir(exist_ok=True)
+    job_tree, levels = build_job_tree()
+
+    init_template = 'isONclust2 cluster -x {} -v -Q -l batches/' \
+                    'isONbatch_{}.cer -o clusters/isONcluster_{}.cer {}; ' \
+                    'sync;\n'
+    template = 'isONclust2 cluster -x {} -v -Q -l clusters/isONcluster_{}' \
+               '.cer -r clusters/isONcluster_{}.cer -o clusters/isONcluster' \
+               '_{}.cer {}; sync\n'
+
+    for nr, l in levels.items():
+        jobs_out = 'jobs_level_{}.sh'.format(nr)
+        with open(jobs_out, 'w') as fh:
+            for n in l:
+                purge = "-z" if n.RightSide else ""
+                if nr == 0 or n.Left is None or n.Right is None:
+                    jr = init_template.format('sahlin', n.Id, n.Id, purge)
+                    fh.write(jr)
+                else:
+                    jr = template.format('sahlin', n.Left.Id,
+                                         n.Right.Id, n.Id, purge)
+                    fh.write(jr)
+        # Run a level in parallel
+        cmd = "parallel < {}".format(jobs_out)
+        sub.call(cmd, shell=True)
+    sub.call("ln -s `realpath clusters/isONcluster_{}"
+             ".cer` isONcluster_ROOT.cer".format(n.Id), shell=True)
+
+
+if __name__ == '__main__':
+    # The cwd should be the process dir that contains 'batches/'
+    main()

bin/table_template.html

@@ -0,0 +1,106 @@
+<script type="text/javascript">
+$(document).ready(function(){
+    $("#load_msg").hide();
+    $("#{{ table_id }}").show();
+});
+</script>
+
+<style>
+table{
+table-layout: fixed;
+word-wrap: break-word;
+}
+#{{ table_id }}{
+<!-- To make the column widths smaller -->
+
+font-family: Arial, Helvetica, sans-serif;
+border-collapse: collapse;
+width: 100%;
+display: none;
+}
+#{{ table_id }} td, #{{ table_id }} th {
+border: 1px solid #ddd;
+padding: 4px;
+}
+#{{ table_id }} tr:nth-child(even){background-color: #f2f2f2;}
+#{{ table_id }} tr:hover {background-color: #90C5E7;}
+#{{ table_id }} th {
+padding-top: 4px;
+padding-bottom: 4px;
+text-align: left;
+background-color: #0084A9;
+}
+</style>
+
+<body>
+<div id='load_msg'>Table loading</div>
+{{ dataframe }}
+</body>
+<script type="text/javascript">
+    $(document).ready(function () {
+    // Setup - add a text input to each footer cell
+    $('#{{ table_id }} thead tr')
+        .clone(true)
+        .addClass('filters')
+        .appendTo('#{{ table_id }} thead');
+
+    var table = $('#{{ table_id }}').DataTable({
+        "columnDefs": [
+                { "width": "5%", "targets": [2, 4] }],
+        sDom: 'lrtip', // removes search box while still allowing search
+        searching: true,
+        pageLength: 30,
+        orderCellsTop: true,
+        fixedHeader: true,
+        initComplete: function () {
+            var api = this.api();
+
+            // For each column
+            api
+                .columns()
+                .eq(0)
+                .each(function (colIdx) {
+                    // Set the header cell to contain the input element
+                    var cell = $('.filters th').eq(
+                        $(api.column(colIdx).header()).index()
+                    );
+                    var title = $(cell).text();
+                    $(cell).html('<input type="text" placeholder="&#xF002; Search" style="font-family:Arial, FontAwesome; width:100%" />');
+
+                    // On every keypress in this input
+                    $(
+                        'input',
+                        $('.filters th').eq($(api.column(colIdx).header()).index())
+                    )
+                        .off('keyup change')
+                        .on('keyup change', function (e) {
+                            e.stopPropagation();
+
+                            // Get the search value
+                            $(this).attr('title', $(this).val());
+                            var regexr = '({search})'; //$(this).parents('th').find('select').val();
+
+                            var cursorPosition = this.selectionStart;
+                            // Search the column for that value
+                            api
+                                .column(colIdx)
+                                .search(
+                                    this.value != ''
+                                        ? regexr.replace('{search}', '(((' + this.value + ')))')
+                                        : '',
+                                    this.value != '',
+                                    this.value == ''
+                                )
+                                .draw();
+
+                            $(this)
+                                .focus()[0]
+                                .setSelectionRange(cursorPosition, cursorPosition);
+                        });
+                });
+        },
+    });
+});
+
+</script>
+

data/OPTIONAL_FILE_1

@@ -1,11 +0,0 @@
-# Nothing to see here. A sentinel file to replace real data.
-# e.g.:
-#
-# input:
-#     file some_data
-#     file extra_data
-# script:
-# def extra = extra_data.name != 'OPTIONAL_FILE' ? "--extra-data $opt" : ''
-# """
-# command ${some_data} ${extra}
-# """

denovo.nf

@@ -1,205 +0,0 @@
-import sys
-import time
-import os
-from os import path
-import pandas as pd
-from glob import glob
-import re
-from threading import Lock
-from itertools import zip_longest
-from collections import OrderedDict, namedtuple
-
-if not workflow.overwrite_configfiles:
-    configfile: "config.yml"
-
-WORKDIR = path.abspath(path.join(config["workdir_top"], config["pipeline"]))
-workdir: WORKDIR
-SNAKEDIR = path.dirname(workflow.snakefile)
-
-include: "snakelib/utils.snake"
-
-in_fastq = config["reads_fastq"]
-if not path.isabs(in_fastq):
-	in_fastq = path.join(SNAKEDIR, in_fastq)
-
-
-class Node:
-    def __init__(self, Id, File, Left, Right, Parent, Level):
-        self.Id = Id
-        self.File = File
-        self.Left = Left
-        self.Right = Right
-        self.Parent = Parent
-        self.Level = Level
-        self.Done = False
-        self.RightSide = False
-    def __repr__(self):
-        return "Node:{} Level: {} File: {} Done: {} Left: {} Right: {} Parent: {}".format(self.Id, self.Level, self.File, self.Done, self.Left.Id if self.Left is not None else None, self.Right.Id if self.Right is not None else None, self.Parent.Id if self.Parent is not None else None)
-
-def grouper(n, iterable, fillvalue=None):
-    "grouper(3, 'ABCDEFG', 'x') --> ABC DEF Gxx"
-    args = [iter(iterable)] * n
-    return zip_longest(fillvalue=fillvalue, *args)
-
-def build_job_tree(batch_dir):
-    batches = glob("{}/isONbatch_*.cer".format(batch_dir))
-    batch_ids = [int(re.search('/isONbatch_(.*)\.cer$', x).group(1)) for x in batches]
-    global LEVELS
-    global JOB_TREE
-    LEVELS = OrderedDict()
-    LEVELS[0] = []
-    for Id, bf in sorted(zip(batch_ids, batches), key=lambda x: x[0]):
-        n = Node(Id, "clusters/isONcluster_{}.cer".format(Id), None, None, None, 0)
-        n.Done = True
-        JOB_TREE[Id] = n
-        LEVELS[0].append(n)
-    level = 0
-    max_id = LEVELS[0][-1].Id
-    while len(LEVELS[level]) != 1:
-        next_level = level + 1
-        LEVELS[next_level] = []
-        for l, r in grouper(2, LEVELS[level]):
-            if r is None:
-                LEVELS[level].pop()
-                l.Level += 1
-                LEVELS[next_level].append(l)
-                continue
-            max_id += 1
-            new_batch = "clusters/isONcluster_{}.cer".format(max_id)
-            new_node = Node(max_id, new_batch, l, r, None, next_level)
-            l.Parent = new_node
-            r.Parent =  new_node
-            r.RightSide = True
-            LEVELS[next_level].append(new_node)
-            JOB_TREE[max_id] = new_node
-        level = next_level
-    global ROOT
-    ROOT = JOB_TREE[len(JOB_TREE)-1].Id
-    JOB_TREE[ROOT].RightSide = True
-    print("Merge clustering job tree nodes:",file=sys.stderr)
-    for n in JOB_TREE.values():
-        print("\t{}".format(n),file=sys.stderr)
-
-
-def generate_rules(levels, snk):
-    init_template = """
-rule cluster_job_%d:
-    input:
-        left = "sorted/batches/isONbatch_%d.cer",
-    output: "clusters/isONcluster_%d.cer"
-    shell: "isONclust2 cluster -x %s -v -Q -l %s -o %s %s; sync"
-
-    """
-    template = """
-rule cluster_job_%d:
-    input:
-        left = "clusters/isONcluster_%d.cer",
-        right = "clusters/isONcluster_%d.cer",
-    output: "clusters/isONcluster_%d.cer"
-    shell: "isONclust2 cluster -x %s -v -Q -l %s -r %s -o %s %s; sync"
-
-    """
-
-    link_template="""
-rule link_root:
-    input: "clusters/isONcluster_%d.cer",
-    output: "clusters/isONcluster_ROOT.cer"
-    shell: "ln -s `realpath {input}` {output}"
-
-    """
-
-    fh = open(snk, "w")
-    for nr, l in levels.items():
-        for n in l:
-            purge = "-z" if n.RightSide else ""
-            if nr == 0 or n.Left is None or n.Right is None:
-                jr = init_template % (n.Id, n.Id, n.Id, config["cls_mode"], "{input.left}", "{output}", purge)
-                fh.write(jr)
-            else:
-                jr = template % (n.Id, n.Left.Id, n.Right.Id, n.Id, config["cls_mode"], "{input.left}", "{input.right}", "{output}", purge)
-                fh.write(jr)
-    global ROOT
-    fh.write(link_template % ROOT)
-    fh.write("\nROOT = %d" % ROOT)
-    fh.flush()
-    fh.close()
-
-def count_fastq_bases(fname, size=128000000):
-    fh = open(fname, "r")
-    count = 0
-    while True:
-        b = fh.read(size)
-        if not b:
-            break
-        count += b.count("A")
-        count += b.count("T")
-        count += b.count("G")
-        count += b.count("C")
-        count += b.count("U")
-    fh.close()
-    return count
-
-def preprocess_reads(fq):
-    pc_opts = config["pychopper_opts"]
-    concat = config["concatenate"]
-    thr = config["cores"]
-
-    out_fq = "processed_reads/full_length_reads.fq"
-    if os.path.isdir("processed_reads"):
-        return out_fq
-
-    shell("mkdir -p processed_reads")
-    if concat:
-        print("Concatenating reads under directory: " + fq)
-        shell("find %s  -regextype posix-extended -regex '.*\.(fastq|fq)$' -exec cat {{}} \\; > processed_reads/input_reads.fq" % fq)
-    else:
-        shell("ln -s `realpath %s` processed_reads/input_reads.fq" % fq)
-
-    if config["run_pychopper"]:
-        print("Running pychopper of fastq file: processed_reads/input_reads.fq")
-        shell("(cd processed_reads; cdna_classifier.py -t %d %s input_reads.fq full_length_reads.fq)" % (thr, pc_opts))
-    else:
-        shell("ln -s `realpath processed_reads/input_reads.fq` processed_reads/full_length_reads.fq")
-
-    return out_fq
-
-
-ROOT = None
-DYNAMIC_RULES="job_rules.snk"
-if ((not os.path.isfile(os.path.join(WORKDIR,"sorted","sorted_reads.fastq"))) or (not os.path.isfile(os.path.join(SNAKEDIR, DYNAMIC_RULES)))):
-    print("Preprocessing read in fastq file:", in_fastq)
-    proc_fastq = preprocess_reads(in_fastq)
-    print("Counting records in input fastq:", proc_fastq)
-    nr_bases = count_fastq_bases(proc_fastq)
-    print("Bases in input: {} megabases".format(int(nr_bases/10**6)))
-    if config['batch_size'] < 0:
-        config['batch_size'] = int(nr_bases/1000/config["cores"])
-
-    print("Batch size is: {}".format(config['batch_size']))
-
-    init_cls_options = """ --batch-size {} --kmer-size {} --window-size {} --min-shared {} --min-qual {}\
-                         --mapped-threshold {} --aligned-threshold {} --min-fraction {} --min-prob-no-hits {} -M {} -P {} -g {} -c {} -F {} """
-    init_cls_options = init_cls_options.format(config["batch_size"], config["kmer_size"], config["window_size"], config["min_shared"], config["min_qual"], \
-                        config["mapped_threshold"], config["aligned_threshold"], config["min_fraction"], config["min_prob_no_hits"], config["batch_max_seq"], config["consensus_period"],
-    config["consensus_minimum"], config["consensus_maximum"], config["min_left_cls"])
-
-    shell("""
-            rm -fr clusters sorted
-            mkdir -p sorted; isONclust2 sort {} -v -o sorted {};
-            mkdir -p clusters;
-        """.format(init_cls_options, proc_fastq))
-    JOB_TREE = OrderedDict()
-    LEVELS = None
-
-    build_job_tree("sorted/batches")
-    generate_rules(LEVELS, "{}/job_rules.snk".format(SNAKEDIR))
-
-include: DYNAMIC_RULES
-
-
-
-rule all:
-    input: rules.link_root.output
-    output: directory("final_clusters")
-    shell:
-        """ isONclust2 dump -v -i sorted/sorted_reads_idx.cer -o final_clusters {input}; sync """

denovo_assembly.nf

@@ -0,0 +1,295 @@
+import groovy.json.JsonSlurper
+import nextflow.util.BlankSeparatedList;
+
+
+map_sample_ids_cls = {it -> 
+/* Harmonize tuples
+output:
+    tuple val(sample_id), path('*.gff')
+When there are multiple paths, will emit:
+    [sample_id, [path, path ..]]
+when there's a single path, this:
+    [sample_id, path]
+This closure makes both cases:
+    [[sample_id, path][sample_id, path]].
+*/
+    if (it[1].getClass() != java.util.ArrayList){
+        // If only one path, `it` will be [sample_id, path]
+        return [it]
+    }
+    l = [];
+    for (x in it[1]){
+        l.add(tuple(it[0], x))
+    }
+    return l
+}
+
+
+process dump_clusters {
+    label "isoforms"
+    input:
+        tuple val(sample_id), path(root_cluster), path(sorted_reads_dir)
+    output:
+        tuple val(sample_id), path("final_clusters"), emit: final_clusters_dir
+        tuple val(sample_id), path("final_clusters/cluster_fastq/*.fq"), emit: final_clusters
+    shell:
+    """ isONclust2 dump -v -i $sorted_reads_dir/sorted_reads_idx.cer -o final_clusters $root_cluster; sync """
+
+}
+
+process build_backbones {
+    /*
+    This step can fail in what seems to at the racon stage giving an 'empty overlap error set' message
+    As a temporary fix, do racon_cmd || true to prevent it crashing the pipeline, and then move on to next cluster
+
+    This process needs some work
+    */
+    label "isoforms"
+	input:
+		tuple val(sample_id), path(cluster_fq)
+	output:
+		tuple val(sample_id), path("*final_polished_cds.fa"), emit: polished_cds
+	script:
+	def cluster_fq_bl = new BlankSeparatedList(cluster_fq)
+	"""
+    # Get one of the cluster ids to give the output a unique name
+    UNID=\$(echo ${cluster_fq_bl[1]} | grep -o -E '[0-9]+')
+
+	for cluster in $cluster_fq_bl
+      do
+        clfq=`basename \$cluster`
+        cln=\${clfq%.*}
+
+        echo Building backbone for cluster: \$cln
+        echo "\tSampling input reads for backbone construction."
+        sample=\${cln}_sample.fq
+        seqkit head --quiet -n 100 \$cluster    > \$sample
+        seqkit sample --quiet -n 500 -2 -s 100 \$cluster    >> \$sample
+        echo "\tConstructing spoa consensus."
+        spoa_cons=\${cln}_spoa.fa
+        spoa -m 5 -n -4 -g -8 -e -6 -q -10 -c -15 -l 1 -r 0 \$sample > \$spoa_cons
+
+        echo "\tPolishing the consensus using racon."
+
+        # polish 1
+        samgz=\${cln}_aln.sam.gz
+        racon_cons=\${cln}_racon.fa
+        racon_cons1=\${cln}_racon1.fa
+
+        minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$spoa_cons \$sample | gzip - > \$samgz
+        racon -t ${params.threads} --no-trimming -u -w 2000 \$sample \$samgz \$spoa_cons > \$racon_cons || true
+        if [ -f \$racon_cons ];
+        then
+            cat \$racon_cons | seqkit replace -p ".*" -r cluster_\${cln} > \$racon_cons1
+        else
+            continue
+        fi
+
+        # polish 2
+        samgz1=\${cln}_aln_1.sam.gz
+        racon_cons2=\${cln}_racon2.fa
+        minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons1 \$sample  | gzip - > \$samgz1
+        racon -t ${params.threads}  -u --no-trimming \$sample \$samgz1 \$racon_cons1 > \$racon_cons2 || true
+        if [ -f \$racon_cons2 ];
+          then
+           echo "success polish 2"
+        else
+            continue
+        fi
+
+         # polish 3
+        samgz2=\${cln}_aln_2.sam.gz
+        minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} \$racon_cons2 \$sample | gzip - > \$samgz2
+        EXITCODE=0
+        racon -t ${params.threads}  -u \$sample \$samgz2 \$racon_cons2 >> ${sample_id}_\${UNID}_final_polished_cds.fa || EXITCODE=\$?
+        if [ \$EXITCODE -eq 0 ];
+          then
+            echo "finished"
+        fi
+
+      done
+    echo "Finished backbones"
+	"""
+
+
+}
+
+process merge_cds {
+    label "isoforms"
+    input:
+        tuple val(sample_id), path(cds)
+    output:
+        tuple val(sample_id), path("${sample_id}_cds.fa"), emit: final_polished_cds
+    script:
+    def merge_list
+    """
+    for FILE in *final_polished_cds.fa
+    do
+        cat \$FILE  >> "${sample_id}_cds.fa"
+    done
+    """
+}
+
+process cds_align {
+    label "isoforms"
+    input:
+        tuple val(sample_id), path(polished_cds), path(sorted_reads_dir)
+    output:
+        tuple val(sample_id), path("${sample_id}_reads_aln_sorted.bam"), emit: bam
+        tuple val(sample_id), path("${sample_id}_read_aln_stats.tsv"), emit: stats
+    script:
+    def ab = "${sample_id}_reads_aln_sorted.bam"
+    """
+		minimap2 -t ${params.threads} \
+		-ax splice ${params.minimap2_opts} $polished_cds $sorted_reads_dir/sorted_reads.fastq |\
+		 samtools view -b - |\
+		samtools sort -o $ab;
+		samtools index $ab;
+		((seqkit bam -s -j ${params.threads} ${ab} 2>&1)  | tee ${sample_id}_read_aln_stats.tsv ) || true
+	"""
+}
+
+
+process make_batches {
+    /*
+    Take a fasta file and creates batches for isONclust2 to work with
+
+    */
+    label "isoforms"
+    input:
+        tuple val(sample_id), path(fastq)
+    output:
+        tuple val(sample_id), path('sorted/batches'), emit: sorted_batches
+        tuple val(sample_id), path('sorted'), emit: sorted_reads_dir
+    
+    script:
+
+    maxcpus = Runtime.runtime.availableProcessors()
+
+    minimum_batch_size = 2000
+    """
+    nr_bases=\$(seqkit stats -T $fastq|cut -f 5| sed '2q;d')
+
+    b=0
+    if [ ${params.batch_size} -lt \$b ];
+    then
+        nr_bases=\$(seqkit stats -T $fastq|cut -f 5| sed '2q;d')
+        let batch_size=\$nr_bases/1000/$maxcpus
+        if [ \$batch_size -lt $minimum_batch_size ];
+        then
+            batch_size=$minimum_batch_size
+        fi
+    else
+        batch_size=${params.batch_size}
+    fi
+
+    echo "Batch size:\$batch_size";
+    echo "Num bases: \$nr_bases";
+
+    init_cls_options="--batch-size \$batch_size --kmer-size ${params.kmer_size} \
+    --window-size ${params.window_size} --min-shared ${params.min_shared} --min-qual ${params.min_qual} \
+                         --mapped-threshold ${params.mapped_threshold} --aligned-threshold ${params.aligned_threshold} \
+                          --min-fraction ${params.min_fraction} --min-prob-no-hits ${params.min_prob_no_hits} \
+                          -M ${params.batch_max_seq} -P ${params.consensus_period} -g ${params.consensus_minimum} -c ${params.consensus_maximum} -F ${params.min_left_cls} "
+    mkdir -p sorted; isONclust2 sort \$init_cls_options -v -o sorted $fastq;
+    """
+
+}
+
+process clustering() {
+    label "isoforms"
+
+    input:
+        tuple val(sample_id), path(sorted_batches)
+    output:
+        tuple val(sample_id), path('isONcluster_ROOT.cer'), emit: root_cluster
+    script:
+    """
+    run_isonclust2.py $sorted_batches
+
+    """
+}
+
+process cluster_quality() {
+    // Run Kristoffer Sahlin's QC code.
+    // For now just write out the PDF and CSV results. Move these to the report at some point
+    label "isoforms"
+
+    input:
+        tuple val(sample_id), path(reads_fl), path(final_clusters_dir)
+        path reference
+    output:
+        tuple val(sample_id),
+            path("${sample_id}_cluster_qc"), emit: cluster_qc_dir
+        tuple val(sample_id),
+            path("${sample_id}_cluster_qc_raw"), emit: cluster_qc_raw
+    script:
+    def qc_dir = "${sample_id}_cluster_qc"
+    def qc_dir_raw = "${sample_id}_cluster_qc_raw" // To generate plots in report
+    def bam = "${qc_dir}/ref_aln.bam"
+
+    """
+    echo $qc_dir_raw
+    mkdir $qc_dir
+    mkdir $qc_dir_raw
+    minimap2 -ax splice -t 2 $reference $reads_fl |\
+     samtools view -q 2 -F 2304 -b - |\
+     samtools sort - -o $bam;
+    samtools index $bam;
+    compute_cluster_quality.py --sizes $final_clusters_dir/clusters_info.tsv \
+     --outfile ${qc_dir}/cluster_quality.csv --ont --clusters $final_clusters_dir/clusters.tsv \
+     --classes $bam --report ${qc_dir}/cluster_quality.pdf --raw_data_out $qc_dir_raw
+    """
+
+}
+
+workflow denovo_assembly {
+    take:
+       fastq_reads_fl
+       reference
+    main:
+
+        make_batches(fastq_reads_fl)
+
+        clustering(make_batches.output.sorted_batches)
+
+        dump_clusters(clustering.output.root_cluster
+                      .join(make_batches.output.sorted_reads_dir))
+
+        build_backbones(dump_clusters.output.final_clusters
+                                        .flatMap(map_sample_ids_cls)
+                                        .groupTuple(size: 10, remainder: true)
+                        )
+
+        merge_cds(build_backbones.output.polished_cds
+                                            .flatMap(map_sample_ids_cls)
+                                            .groupTuple()
+                  )
+
+        cds_align(merge_cds.out.final_polished_cds.view()
+                  .join(make_batches.output.sorted_reads_dir))
+
+        if (!reference.name.startsWith('OPTIONAL_FILE')){
+
+            cluster_quality(fastq_reads_fl
+                .join(dump_clusters.output.final_clusters_dir), reference)
+
+            cluster_quality.out.cluster_qc_dir
+                .set { opt_qual_ch }
+
+            cluster_quality.output.cluster_qc_raw
+                .set { opt_qual_raw_ch }
+
+        } else{
+            Channel.empty().set { opt_qual_ch }
+            Channel.empty().set { opt_qual_raw_ch }
+        }
+
+    emit:
+       bam = cds_align.output.bam
+       cds = merge_cds.out.final_polished_cds
+       stats = cds_align.output.stats
+       opt_qual_ch
+       opt_qual_raw_ch
+}
+

environment.yaml

@@ -1,4 +1,4 @@
-name: epi2melabs-wf-isoforms
+name: epi2melabs-wf-isoforms-test
 channels:
     - epi2melabs
     - bioconda
@@ -8,20 +8,21 @@ dependencies:
     - python==3.8.*
     - aplanat >=0.5.0
     - epi2melabs
-    - minimap2
+    - minimap2 ==2.24
-    - samtools
+    - samtools ==1.14
-    - bedtools
+    - bedtools ==2.30.0
-    - pychopper
+    - pychopper==2.5.0
-    - pandas
+    - pandas==1.3.5
-    - seaborn
+    - gffread==0.12.7
-    - requests
+    - gffcompare==0.11.2
-    - gffread
+    - gffutils=0.10.1
-    - seqkit
+    - seqkit==2.1.0
-    - csvtk
     - stringtie==2.1.1
-    - gffcompare
     - curl
-    - pysam
+    - pysam==0.17.0
+    - racon==1.4.20
     - spoa==3.4.0
-    - fastcat
+    - fastcat==0.4.10
-    - isonclust2
+    - isonclust2==2.3
+    - parallel
+    - scikit-learn==1.0.2

evaluation/run_evaluation_dmel.sh

@@ -0,0 +1,59 @@
+#!/usr/bin/env bash
+
+# Usage: ./run_evaluation_dmel.sh pathto/outputdir
+
+# See the isONcorrect paper https://www.nature.com/articles/s41467-020-20340-8 where this dataset is described
+
+
+if [[ "$#" -lt 1 ]]; then
+    echo "usage: run_evaluation_dmel.sh <outdir> [nextflow.config]"
+    exit 1
+fi
+
+if [[ "$#" -eq 1 ]]; then
+    config=''
+fi
+
+if [[ "$#" -eq 2 ]]; then
+  config="-c $2";
+fi
+
+OUTDIR=$1;
+
+FASTQ_URL="http://ftp.sra.ebi.ac.uk/vol1/fastq/ERR358/005/ERR3588905/ERR3588905_1.fastq.gz"
+REF_URL="http://ftp.ensembl.org/pub/release-99/fasta/drosophila_melanogaster/dna/Drosophila_melanogaster.BDGP6.28.dna.toplevel.fa.gz"
+GFF_URL="http://ftp.ensembl.org/pub/release-99/gff3/drosophila_melanogaster/Drosophila_melanogaster.BDGP6.28.99.gff3.gz"
+
+DATA_DIR="$OUTDIR/data"
+READS_DIR="$DATA_DIR/reads"
+FASTQ="$READS_DIR/ERR3588905_1.fastq.gz"
+REF="$DATA_DIR/Drosophila_melanogaster.BDGP6.28.dna.toplevel.fa"
+GFF="$DATA_DIR/Drosophila_melanogaster.BDGP6.28.99.gff3"
+
+mkdir -p $READS_DIR
+
+if [ ! -f $REF ];
+then (echo "downloading reference genome"; cd $DATA_DIR; curl -L -C - -O $REF_URL); gzip -d ${REF}.gz
+fi
+
+if [ ! -f $GFF ];
+then
+	(echo "downloading reference annotation"; cd $DATA_DIR; curl -L -C - -O $GFF_URL); gzip -d ${GFF}.gz
+fi
+
+if [ ! -f $FASTQ ];
+then (echo "downloading reads"; cd $READS_DIR; curl -L -C - -O $FASTQ_URL); gzip -d ${FASTQ}.gz
+fi
+
+
+OUT_REF="$OUTDIR/ref"
+OUT_DENOVO="$OUTDIR/denovo"
+
+
+nextflow run ../ --fastq $READS_DIR  $config \
+--ref_genome $REF --ref_annotation $GFF -profile local --out_dir $OUT_REF --minimap2_opts '-uf --splice-flank=no' \
+-w $OUT_REF/workspace -resume;
+
+echo "Doing de novo evaluation"
+nextflow run ../ --fastq $READS_DIR  $config --denovo -profile local --out_dir $OUT_DENOVO \
+-w $OUT_DENOVO/workspace -resume;

evaluation/tests.sh

@@ -0,0 +1,73 @@
+#!/usr/bin/env bash
+
+# A few simple tests with different combinations of CLI options
+# Run from within an appropriate active conda environment
+
+if [[ "$#" -lt 1 ]]; then
+    echo "usage: tests.sh <outdir> [nextflow.config]"
+    exit 1
+fi
+
+if [[ "$#" -eq 1 ]]; then
+    config=''
+fi
+
+if [[ "$#" -eq 2 ]]; then
+  config="-c $2";
+fi
+
+SCRIPT_DIR="$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )"
+cd $SCRIPT_DIR/../;
+
+singledir="test_data/fastq"
+multisampledir="test_data/demultiplexed_fastq"
+
+# This is for when using SIRV dataset with non-canonical spice junctions
+#"--minimap2_opts '-uf --splice-flank=no'"
+results=()
+
+
+OUTPUT=$1/denovo_multi_sample_no_ref_genome;
+nextflow run . --fastq $multisampledir $config --denovo --ref_genome test_data/SIRV_150601a.fasta  -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
+--sample_sheet test_data/sample_sheet -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+OUTPUT=$1/denovo_single;
+nextflow run . --fastq $singledir $config --denovo --ref_genome test_data/SIRV_150601a.fasta -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
+--sample_sheet test_data/sample_sheet -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+# Reference based tests
+OUTPUT=$1/reference_single_dir;
+nextflow run . --fastq $singledir $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no' \
+--ref_annotation test_data/SIRV_isofroms.gtf -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+OUTPUT=$1/multiple_samples;
+nextflow run . --fastq $multisampledir $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no'\
+--ref_annotation test_data/SIRV_isofroms.gtf -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
+--sample_sheet test_data/sample_sheet -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+OUTPUT=$1/reference_no_ref_annotation;
+nextflow run . --fastq $singledir $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no'\
+ -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+# Force split_bam to make multiple alignment bundles
+OUTPUT=$1/reference_frce_split_bam;
+nextflow run . --fastq $singledir  $config --ref_genome test_data/SIRV_150601a.fasta --minimap2_opts '-uf --splice-flank=no'\
+--ref_annotation test_data/SIRV_isofroms.gtf -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
+--bundle_min_reads 5 -resume;
+r=$?
+results+=("$(basename $OUTPUT): $r")
+
+echo "Exit status codes for each test"
+for value in "${results[@]}"; do
+     echo "${value}"
+done

main.nf

@@ -1,39 +1,37 @@
 #!/usr/bin/env nextflow
 
-// Developer notes
+/* This workflow is a adapted from two previous pipeline written in Snakemake:
-// 
+- https://github.com/nanoporetech/pipeline-nanopore-ref-isoforms
-// This template workflow provides a basic structure to copy in order
+- https://github.com/nanoporetech/pipeline-nanopore-denovo-isoforms
-// to create a new workflow. Current recommended pratices are:
+*/
-//     i) create a simple command-line interface.
-//    ii) include an abstract workflow scope named "pipeline" to be used
-//        in a module fashion.
-//   iii) a second concreate, but anonymous, workflow scope to be used
-//        as an entry point when using this workflow in isolation.
 
-import groovy.json.JsonBuilder
+import groovy.json.JsonBuilder;
+import nextflow.util.BlankSeparatedList;
 import java.util.ArrayList;
 nextflow.enable.dsl = 2
 
 include { fastq_ingress } from './lib/fastqingress' 
 include { start_ping; end_ping } from './lib/ping'
+include { reference_assembly } from './reference_assembly'
+include { denovo_assembly } from './denovo_assembly'
 
 
-process summariseReads {
+process summariseConcatReads {
-    // concatenate fastq and fastq.gz in a dir
+    // concatenate fastq and fastq.gz in a dir write stats
 
     label "isoforms"
     cpus 1
     input:
-        tuple path(directory), val(sample_name), val(type)
+        tuple path(directory), val(sample_id), val(type)
     output:
-        tuple val(sample_name), path('*'), emit: summary
+        tuple val(sample_id), path("${sample_id}.fastq"), emit: input_reads
+        tuple val(sample_id), path('*.stats'), emit: summary
     script:
     """
-    fastcat -s ${sample_name} -r ${sample_name}.stats -x ${directory} > /dev/null
+    fastcat -s ${sample_id} -r ${sample_id}.stats -x ${directory} >  ${sample_id}.fastq
     """
 }
 
-
 process getVersions {
     label "isoforms"
     cpus 1
@@ -44,7 +42,7 @@ process getVersions {
     python -c "import pysam; print(f'pysam,{pysam.__version__}')" >> versions.txt
     python -c "import aplanat; print(f'aplanat,{aplanat.__version__}')" >> versions.txt
     python -c "import pandas; print(f'pandas,{pandas.__version__}')" >> versions.txt
-    python -c "import seaborn; print(f'seaborn,{seaborn.__version__}')" >> versions.txt
+    python -c "import sklearn; print(f'scikit-learn,{sklearn.__version__}')" >> versions.txt
     fastcat --version | sed 's/^/fastcat,/' >> versions.txt
     minimap2 --version | sed 's/^/minimap2,/' >> versions.txt
     samtools --version | head -n 1 | sed 's/ /,/' >> versions.txt
@@ -52,10 +50,10 @@ process getVersions {
     python -c "import pychopper; print(f'pychopper,{pychopper.__version__}')" >> versions.txt
     gffread --version | sed 's/^/gffread,/' >> versions.txt
     seqkit version | head -n 1 | sed 's/ /,/' >> versions.txt
-    csvtk version | head -n 1 | sed 's/ /,/' >> versions.txt
     stringtie --version | sed 's/^/stringtie,/' >> versions.txt
     gffcompare --version | head -n 1 | sed 's/ /,/' >> versions.txt
     spoa --version | sed 's/^/spoa,/' >> versions.txt
+    isONclust2 version | sed 's/ version: /,/' >> versions.txt
     """
 }
 
@@ -80,44 +78,27 @@ process preprocess_reads {
     */
  
     label "isoforms"
-    cpus params.threads
+    cpus 4
 
     input:
-        tuple path(directory), val(sample_id), val(type) 
+        tuple val(sample_id), path(input_reads)
     output:
-         tuple val(sample_id), path("full_length_reads.fq"), emit: full_len_reads
+         tuple val(sample_id), path("${sample_id}_full_length_reads.fq"), emit: full_len_reads
-         tuple val(sample_id), path('cdna_classifier_report.tsv'),  emit: report
+         tuple val(sample_id), path('*.tsv'),  emit: report
-//          val "${sample_id}", emit: sample_id
-//          path 'cdna_classifier_report.tsv', optional: true, emit: cdna_class_report
-         // Not sure if this is an antipattern, but it's function is to publish all the files to the publishDir dir
     script:
         """
-        fastcat -s ${sample_id} -r ${sample_id}.stats -x ${directory} > input_reads.fq
-    
         if [[ ${params.use_pychopper} == true ]];
         then
-            cdna_classifier.py -t ${params.threads} ${params.pychopper_opts} input_reads.fq full_length_reads.fq
+            cdna_classifier.py -t ${params.threads} ${params.pychopper_opts} ${input_reads} ${sample_id}_full_length_reads.fq
-            generate_pychopper_stats.py --data cdna_classifier_report.tsv --output .
+            mv cdna_classifier_report.tsv ${sample_id}_cdna_classifier_report.tsv
+            generate_pychopper_stats.py --data ${sample_id}_cdna_classifier_report.tsv --output .
         else
-            ln -s `realpath input_reads.fq` full_length_reads.fq
+            ln -s `realpath $input_reads` "${sample_id}_full_length_reads.fq"
+            touch $sample_id}_cdna_classifier_report.tsv
         fi
         """
 }
 
-process generate_fq_stats {
-    label "isoforms"
-
-    input:
-        tuple(sample_id), path(fastq), val(unused)
-
-    output:
-        path "*"
-    script:
-    """
-        run_fastq_qc.py --fastq ${fastq} --output .
-    """
-}
-
 process build_minimap_index{
     /*
     Build minimap index from reference genome
@@ -126,7 +107,7 @@ process build_minimap_index{
     cpus params.threads
 
     input:
-        file reference
+        path reference
     output:
         path "genome_index.mmi", emit: index
     script:
@@ -135,66 +116,6 @@ process build_minimap_index{
     """
 }
 
-
-process map_reads{
-    /*
-    Map reads to reference using minimap2.
-    Filter reads by mapping quality.
-    Filter reads where length of poly(A) > max_poly_run at either ends of the read (defined by poly_context)
-    */
-    label "isoforms"
-    cpus params.threads
-
-    input:
-       file index
-       file reference
-       tuple val(sample_id), file(fastq_reads)
-
-    output:
-       tuple val(sample_id), path("reads_aln_sorted.bam"), emit: bam
-       tuple val(sample_id), path("read_aln_stats.tsv"), emit: stats
-    script:
-        def ab = "reads_aln_sorted.bam"
-        def af = "internal_priming_fail.tsv"
-        def fs = "context_internal_priming_fail_start.fasta"
-        def fe = "context_internal_priming_fail_end.fasta"
-        def fasta_reads = "reads.fa"
-        def ContextFilter = """AlnContext: { Ref: "${reference}", LeftShift: -${params.poly_context}, RightShift: ${params.poly_context},
-           RegexEnd: "[Aa]{${params.max_poly_run},}",
-           Stranded: True,Invert: True, Tsv: "internal_priming_fail.tsv"} """
-       """
-       seqkit fq2fa ${fastq_reads} -o ${fasta_reads};
-       minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} ${index} ${fasta_reads}\
-       | samtools view -q ${params.minimum_mapping_quality} -F 2304 -Sb -\
-       | seqkit bam -j ${params.threads} -x -T '${ContextFilter}' -\
-       | samtools sort -@ ${params.threads} -o ${ab} -;
-
-       ((seqkit bam -s -j ${params.threads} ${ab} 2>&1)  | tee read_aln_stats.tsv ) || true
-
-       if [[ -s ${af} ]];
-       then
-           tail -n +2 ${af} | awk '{{print ">" \$1 "\\n" \$4 }}' - > ${fs}
-           tail -n +2 ${af} | awk '{{print ">" \$1 "\\n" \$6 }}' - > ${fe}
-       fi
-       """
-}
-
-process plot_aln_stats{
-    /*
-    Create a pdf of alignemnt statistis. 
-    */
-    label 'isoforms'
-
-    input:
-        tuple val(sample_id), path (aln_stats)
-    output:
-        path "*"
-    script:
-    """
-    plot_aln_stats.py ${aln_stats} -r read_aln_stats.pdf
-    """
-}
-
 process split_bam{
     /*
     Partition BAM file into loci or bundles with `params.bundle_min_reads` minimum size
@@ -215,19 +136,21 @@ process split_bam{
         """
         seqkit bam -j ${params.threads} -N ${params.bundle_min_reads} ${bam} -o  bam_bundles/
         mv bam_bundles/* .
+        for f in *:*; do mv -v "\$f" \$(echo "\$f" | tr ':' '-'); done
         """
     else
         """
         mkdir -p ./${sample_id}_bam_bundles
-        ln -s ${bam} ${sample_id}_bam_bundles/000000000_ALL:0:1_bundle.bam
+        ln -s ${bam} ${sample_id}_bam_bundles-000000000-ALL-0-1_bundle.bam
         """
 }
 
 
-process stringtie{
+process assemble_transcripts{
     /*
-    Takes in aligned reads in bam format that may be a chunk of a larger alignment file.
+    Assemble transcripts using stringtie.
-    $G_FLAG specifies whether or not to use reference annotation as a guide in transcript assembly. 
+    Take aligned reads in bam format that may be a chunk of a larger alignment file.
+    Optionally use reference annotation to guide assembly.
 
     Output gff annotation files in a tuple with `sample_id` for combining into samples late rin the pipeline.
     */
@@ -236,15 +159,14 @@ process stringtie{
 
     input:
         tuple val(sample_id), path(bam)
-        file ref_annotation
+        path ref_annotation
     output:
         tuple val(sample_id), path('*.gff'), emit: gff_bundles
     script:
-        def out_filename = bam.name.replaceFirst(~/\.[^\.]+$/, '') + '.gff'
+        def out_filename = bam.name.replaceFirst(~/\.[^\.]+$/, '') + "_${sample_id}.gff"
-        def label = "STR.${bam.name.split('_')[0].toInteger()}."
         def G_FLAG = ref_annotation.name.startsWith('OPTIONAL_FILE') ? '' : "-G ${ref_annotation}"
     """
-    stringtie --rf ${G_FLAG} -L -v -p ${params.threads} ${params.stringtie_opts} -o  ${out_filename} \
+    stringtie --rf ${G_FLAG} -L -v -A gene_abund.tab -p ${params.threads} ${params.stringtie_opts} -o  ${out_filename} \
     ${bam} 2>/dev/null
     """
 }
@@ -261,9 +183,9 @@ process merge_gff_bundles{
     output:
         tuple val(sample_id), path('*.gff'), emit: gff
     script:
-    def merged_gff = "str_merged_${sample_id}.gff"
+    def merged_gff = "transcripts_${sample_id}.gff"
     """
-    echo '#gff-version 2' >> $merged_gff;
+    echo '##gff-version 2' >> $merged_gff;
     echo '#pipeline-nanopore-isoforms: stringtie' >> $merged_gff;
 
     for fn in ${gff_bundle};
@@ -274,9 +196,11 @@ process merge_gff_bundles{
     """
 }
 
-process run_gff_compare{
+process run_gffcompare{
     /*
     Compare query and reference annotations.
+    If ref_annotation is an optional file, just make an empty directory to satisfy
+    the requirements of the downstream processes.
     */
 
     label 'isoforms'
@@ -285,45 +209,55 @@ process run_gff_compare{
        tuple val(sample_id), path(query_annotation)
        path ref_annotation
     output:
-        tuple val(sample_id), path('str_merged.annotated.gtf'), emit: merged_annotated
+        tuple val(sample_id), path("${sample_id}_gffcompare"), emit: gffcmp_dir
-        tuple val(sample_id), path('str_merged.stats'), emit: stats
-        tuple val(sample_id), path('str_merged.tracking'), emit: tracking
     script:
-        """
+    def out_dir = "${sample_id}_gffcompare"
-        echo "Doing comparison of reference annotation: ${ref_annotation} and the current annotation"
-        gffcompare -o str_merged -r ${ref_annotation} ${params.gffcompare_opts} ${query_annotation}
-        generate_tracking_summary.py --tracking str_merged.tracking --output_dir . --annotation ${ref_annotation}
 
-        if [[ ${params.plot_gffcmp_stats} == true ]];
+    if ( ref_annotation.name.startsWith('OPTIONAL_FILE') ){
-        then
-            plot_gffcmp_stats.py -r str_gffcmp_report.pdf -t str_merged.tracking str_merged.stats;
-        fi
         """
+        mkdir $out_dir
+        """
+    } else {
+        """
+        mkdir $out_dir
+        echo "Doing comparison of reference annotation: ${ref_annotation} and the query annotation"
+
+       gffcompare -o ${out_dir}/str_merged -r ${ref_annotation} \
+        ${params.gffcompare_opts} ${query_annotation}
+
+        generate_tracking_summary.py --tracking $out_dir/str_merged.tracking \
+        --output_dir ${out_dir} --annotation ${ref_annotation}
+
+       mv *.tmap $out_dir
+       mv *.refmap $out_dir
+        """
+    }
 }
 
 
-process run_gffread{
+process get_transcriptome{
         /*
-        Write out a transctiptome file based on the gff annotations.
+        Write out a transcriptome file based on the query gff annotations.
+        TODO: Do we need to touch merged_transcriptome or can we just pass it?
         */
         label 'isoforms'
 
         input:
-            tuple val(sample_id), path(gff_merged), path(merged_ann_gff)
+            tuple val(sample_id), path(transcripts_gff), path(gffcmp_dir), path(reference_seq)
-            path reference_seq
         output:
-            tuple val(sample_id), path('*.fas'), emit: transcriptome
+            tuple val(sample_id), path("*.fas"), emit: transcriptome
+
         script:
-        def str_transcriptome = "${sample_id}_str_transcriptome.fas"
+        def transcriptome = "${sample_id}_transcriptome.fas"
         def merged_transcriptome = "${sample_id}_merged_transcriptome.fas"
         """
 
-        gffread -g ${reference_seq} -w ${str_transcriptome} ${gff_merged}
+        gffread -g ${reference_seq} -w ${transcriptome} ${transcripts_gff}
-        if [ -f ${merged_ann_gff} ]
+        if  [ "\$(ls -A $gffcmp_dir)" ];
         then
-            gffread -F -g ${reference_seq} -w ${merged_transcriptome} ${merged_ann_gff}
+            echo "Yes"
-        else
+            gffread -F -g ${reference_seq} -w ${merged_transcriptome} \
-            touch ${merged_transcriptome}
+            $gffcmp_dir/str_merged.annotated.gtf
         fi
         """
 }
@@ -333,26 +267,45 @@ process makeReport {
     label "isoforms"
 
     input:
+        path report_template
+        path table_template
         path versions
-        path params
+        path "params.json"
-        tuple val(sample_id), path(seqs), path(aln_stats), path(gffcompare_tracking), path(gffcompare_stats),
+        val denovo
-                path(cdna_class_report)
+        tuple val(sample_ids),
+                path(seq_summaries),
+                path(aln_stats),
+                path(gffcmp_dir),
+                path(cdna_class_report),
+                path(gff_annotation)
     output:
-        tuple val(sample_id), path("wf-isoforms-*.html"), emit: report
+        path("wf-isoforms-*.html"), emit: report
     script:
-        def report_name = "wf-isoforms-${sample_id}_report.html"
+        // Convert the sample_id arrayList.
-        def opt_gff_track = gffcompare_tracking.name.startsWith('OPTIONAL_FILE') ? '' : "--gffcompare_tracking ${gffcompare_tracking}"
+        sids = new BlankSeparatedList(sample_ids)
-        def opt_gff_stats = gffcompare_stats.name.startsWith('OPTIONAL_FILE') ? '' : "--gffcompare_stats ${gffcompare_stats}"
+
+        def report_name = "wf-isoforms-report.html"
+        def OPT_ALN = denovo ?  '' : "--alignment_stats ${aln_stats}"
+        def OPT_DENOVO = denovo ? "--denovo" : ''
     """
-    report.py ${report_name} --versions ${versions} ${seqs} --params params.json \
+    report.py --report $report_name \
-    --alignment_stats ${aln_stats} \
+    --report_template $report_template \
-    ${opt_gff_track} \
+    --table_template $table_template \
-    ${opt_gff_stats} \
+    --versions $versions \
-    --pychop_report ${cdna_class_report}
+    --params params.json \
+    $OPT_ALN \
+    --pychop_report $cdna_class_report \
+    --sample_ids $sids \
+    --summaries $seq_summaries \
+    --gffcompare_dir $gffcmp_dir \
+    --gff_annotation $gff_annotation \
+    --transcript_table_cov_thresh $params.transcript_table_cov_thresh \
+    $OPT_DENOVO
+
+
     """
 }
 
-
 // See https://github.com/nextflow-io/nextflow/issues/1636
 // This is the only way to publish files from a workflow whilst
 // decoupling the publish from the process steps.
@@ -362,11 +315,24 @@ process output {
     publishDir "${params.results_dir}/${sample_id}", mode: 'copy', pattern: "*"
 
     input:
-        tuple val(sample_id), file(fname)
+        tuple val(sample_id), path(fname)
     output:
-        file fname
+        path fname
     """
     echo "Writing output files"
+    echo $fname
+    """
+}
+
+process output_report {
+    publishDir "${params.results_dir}", mode: 'copy', pattern: "*report.html"
+
+    input:
+        path fname
+    output:
+        path fname
+    """
+    echo "Copying report"
     """
 }
 
@@ -377,6 +343,8 @@ workflow pipeline {
         reads
         ref_genome
         ref_annotation
+        report_template
+        table_template
     main:
 
         map_sample_ids_cls = {it -> 
@@ -401,60 +369,92 @@ workflow pipeline {
             return l
         }
 
-        summariseReads(reads)
+        summariseConcatReads(reads)
-        sample_ids = summariseReads.out.summary.collect({it -> it[0]})
+        sample_ids = summariseConcatReads.out.summary.flatMap({it -> it[0]})
+
         software_versions = getVersions()
         workflow_params = getParams()
-        preprocess_reads(reads)
+        preprocess_reads(summariseConcatReads.out.input_reads)
-//         generate_fq_stats(preprocess_reads.out.sample)  skip for now
 
-        build_minimap_index(ref_genome)
+        if (params.denovo){
-        map_reads(build_minimap_index.out.index, ref_genome, preprocess_reads.out.full_len_reads)
+            println("Doing de novo assembly")
-//         plot_aln_stats(map_reads.out.bam)
+             m = denovo_assembly(preprocess_reads.out.full_len_reads, ref_genome)
-        split_bam(map_reads.out.bam)
+
-        stringtie(split_bam.out.bundles.flatMap(map_sample_ids_cls), ref_annotation)
+        } else {
-        merge_gff_bundles(stringtie.out.gff_bundles.groupTuple())
+            build_minimap_index(ref_genome)
+            println("Doing reference based transcript analysis")
+            m = reference_assembly(build_minimap_index.out.index, ref_genome, preprocess_reads.out.full_len_reads)
+        }
+
+        split_bam(m.bam)
+
+        assemble_transcripts(split_bam.out.bundles.flatMap(map_sample_ids_cls), ref_annotation)
+
+        merge_gff_bundles(assemble_transcripts.out.gff_bundles.groupTuple())
 
         use_ref_ann = !ref_annotation.name.startsWith('OPTIONAL_FILE')
 
-        if (use_ref_ann){
+        run_gffcompare(merge_gff_bundles.out.gff, ref_annotation)
-            run_gff_compare(merge_gff_bundles.out.gff, ref_annotation)
+
-            run_gffread(merge_gff_bundles.out.gff.join(run_gff_compare.out.merged_annotated),
+        if (params.denovo){
-                        ref_genome)
+            // Use the perd-sample, de novo-assembled CDS
-            gff_tracking = run_gff_compare.out.tracking
+            seq_for_transcriptome_build = m.cds
-            gff_stats = run_gff_compare.out.stats
+        }else {
-        }else{
+            // If doing reference based assembly, there is only one reference
-            // Create dummy file paths to satisfy required path inputs of processes
+            // So map this reference to all sample_ids
-            // Add to tuple with sample_id to guide to correct sample
+            seq_for_transcriptome_build = sample_ids.flatten().combine(Channel.fromPath(params.ref_genome))
-            gff_tracking = sample_ids.combine(Channel.fromPath("$projectDir/data/OPTIONAL_FILE")).view()
-            gff_stats = sample_ids.combine(Channel.fromPath("$projectDir/data/OPTIONAL_FILE_1"))
         }
 
-        makeReport(
+        makeReport(report_template,
+                    table_template,
                     software_versions,
                     workflow_params,
-                    summariseReads.out.summary
+                    params.denovo,
-                    .join(map_reads.out.stats)
+                    summariseConcatReads.out.summary
-                    .join(gff_tracking)
+                    .join(m.stats)
-                    .join(gff_stats)
+                    .join(run_gffcompare.out.gffcmp_dir)
                     .join(preprocess_reads.out.report)
+                    .join(merge_gff_bundles.out.gff)
+                    .toList().transpose().toList()
                     )
+
+        report = makeReport.out.report
+
+       get_transcriptome(merge_gff_bundles.out.gff
+            .join(run_gffcompare.out.gffcmp_dir)
+            .join(seq_for_transcriptome_build))
+
         if (use_ref_ann){
             results = preprocess_reads.out.report
-                        .concat(merge_gff_bundles.out.gff,
+                        .concat(run_gffcompare.output.gffcmp_dir,
-                        run_gff_compare.out.stats,
+                        m.stats,
-                        run_gff_compare.out.merged_annotated,
+                        get_transcriptome.out.flatMap(map_sample_ids_cls),
-                        run_gffread.out.transcriptome,
-                        makeReport.out.report
                         )
-        }else{  // Write a minimal report if no reference annotation is given
+        }
+        if (!use_ref_ann && !params.denovo){
             results = preprocess_reads.out.report
-                        .concat(merge_gff_bundles.out.gff,
+                        .concat(m.stats,
-                                makeReport.out.report
+                        get_transcriptome.out.flatMap(map_sample_ids_cls),
-                                )
+                        )
-            }
+        }
+        if (params.denovo){
+            results = m.cds
+                      .concat(m.stats,
+                      seq_for_transcriptome_build,
+                      get_transcriptome.out.flatMap(map_sample_ids_cls),
+                      merge_gff_bundles.out.gff,
+                      m.opt_qual_ch.flatMap {it -> 
+                      l = []
+                       for (x in it[1..-1]){
+                            l.add(tuple(it[0], x))
+                            }
+                        return l
+                      })
+        }
+
     emit:
         results
+        report
         telemetry = workflow_params
 }
 
@@ -465,6 +465,9 @@ workflow {
     start_ping()
     params.results_dir = "${params.out_dir}/output"
 
+    report_template = file("$projectDir/bin/report_template.html")
+    table_template = file("$projectDir/bin/table_template.html")
+
     fastq = file(params.fastq, type: "file")
 
      if (!fastq.exists()) {
@@ -472,15 +475,27 @@ workflow {
         exit 1
     }
 
+    if (!params.denovo && !params.ref_genome){
+        println("--ref_genome must be supplied unless doing de novo assembly (--denovo)")
+        exit 1
+    }
+
+
     if (params.ref_genome){
         ref_genome = file(params.ref_genome, type: "file")
         if (!ref_genome.exists()) {
-            println("--reference: File doesn't exist, check path.")
+            println("--ref_genome: File doesn't exist, check path.")
             exit 1
         }
+    }else {
+        ref_genome = file("$projectDir/data/OPTIONAL_FILE")
+    }
+
+    if (params.denovo && params.ref_annotation) {
+        println("Reference annotation with de denovo assembly is not supported")
+        exit 1
     }
 
-    ref_annotation = null
     if (params.ref_annotation){
         ref_annotation = file(params.ref_annotation, type: "file")
         if (!ref_annotation.exists()) {
@@ -495,11 +510,10 @@ workflow {
         params.fastq, params.out_dir, params.sample, params.sample_sheet, params.sanitize_fastq
     )
 
-    pipeline(reads, ref_genome, ref_annotation)
+    pipeline(reads, ref_genome, ref_annotation, report_template, table_template)
 
-    output(
+    output(pipeline.out.results)
-        pipeline.out.results
+    output_report(pipeline.out.report)
-        )
 
     end_ping(pipeline.out.telemetry)
 }

nextflow.config

@@ -13,17 +13,19 @@
 params {
     help = false
     fastq = null
-    ref_genome = null
+    ref_genome = false
     ref_annotation = null
     // Process cDNA reads using pychopper, turn off for direct RNA:
     use_pychopper = true
     threads = 4
+     // Thresholds for viewing isoforms in report table
+    transcript_table_cov_thresh = 50
 
     out_dir = null
     sample = null
     sample_sheet = null
     sanitize_fastq = false
-    wfversion = "v0.0.1"
+    wfversion = "v0.1.0"
     aws_image_prefix = null
     aws_queue = null
     report_name = "report"
@@ -35,15 +37,16 @@ params {
 
 
     // Options passed to pychopper:
-    pychopper_opts = ""
+    pychopper_opts = "-m edlib"
 
     // Extra option passed to minimap2 when generating index
     minimap_index_opts = "-k14"
 
     // Extra options passed to minimap2
+    // For SIRV data
+    //minimap2_opts = "-uf --splice-flank=no"
+    // AFor non-SIRV data:
     minimap2_opts = "-uf"
-    // Add this for SIRV data:
-    // "--splice-flank=no"
 
     // Minmum mapping quality
     minimum_mapping_quality = 40
@@ -67,6 +70,56 @@ params {
     plot_gffcmp_stats = true
 
     disable_ping = false
+
+    //// Denovo-specific parameters
+    denovo = false
+
+    // Batch size in kilobases (if -1 then it is calculated based on the number of cores and bases):
+    batch_size = -1
+
+    // Maximum sequences per input batch (-1 means no limit):
+    batch_max_seq = -1
+
+    // Clustering mode:
+    cls_mode = "sahlin"
+
+    // Kmer size:
+    kmer_size = 11
+
+    // Window size:
+    window_size = 15
+
+    // Minimum cluser size in the left batch:
+    min_left_cls = 2
+
+    // Consensus period (-1 means no consensus):
+    consensus_period = 500
+
+    // Minimum consensus sample size:
+    consensus_minimum = 50
+
+    // Maximum consensus sample size:
+    consensus_maximum = -150
+
+    // Minimum number of minimizers shared between read and cluster:
+    min_shared = 5
+
+    // Minimum average quality value:
+    min_qual = 7.0
+
+    // Minmum mapped fraction of read to be included in cluster:
+    mapped_threshold = 0.65
+
+    // Minimum aligned fraction of read to be included in cluster:
+    aligned_threshold = 0.2
+
+    // Minimum fraction of minimizers shared compared to best hit, in order to continue mapping:
+    min_fraction = 0.8
+
+    // Minimum probability for i consecutive minimizers to be different between read and representative:
+    min_prob_no_hits = 0.1
+
+
 }
 
 manifest {
@@ -76,7 +129,7 @@ manifest {
     description     = 'RNA/cDNA isoform analysis workflow'
     mainScript      = 'main.nf'
     nextflowVersion = '>=20.10.0'
-    //version       = 'v0.0.7'  // TODO: do switch to this?
+    //version       = 'v0.0.1'
 }
 
 executor {

nextflow_schema.json

@@ -63,7 +63,8 @@
                 },
                 "pychopper_opts": {
                     "type": "string",
-                    "description": "Extra pychopper opts"
+                    "description": "Extra pychopper opts",
+                    "default": "-m edlib"
                 },
                 "threads": {
                     "type": "integer",
@@ -110,11 +111,95 @@
                 },
                 "disable_ping": {
                     "type": "boolean"
+                },
+                "transcript_table_cov_thresh": {
+                    "type": "integer",
+                    "description": "Minimum coverage for a transcript to appear in the report table",
+                    "default": 50
+                },
+                "denovo": {
+                    "type": "boolean",
+                    "description": "Use denovo transcript assembly rather than reference guided",
+                    "default": false
+                },
+                "batch_size": {
+                    "type": "integer",
+                    "description": "Maximum sequences per input batch (-1 means no limit)",
+                    "default": -1
+                },
+                "batch_max_seq": {
+                    "type": "integer",
+                    "description": "Maximum sequences per input batch (-1 means no limit)",
+                    "default": -1
+                },
+                "cls_mode": {
+                    "type": "string",
+                    "description": "Clustering mode",
+                    "default": "sahlin"
+                },
+                "kmer_size": {
+                    "type": "integer",
+                    "description": "Kmer size",
+                    "default": 11
+                },
+                "window_size": {
+                    "type": "integer",
+                    "description": "Window size",
+                    "default": 15
+                },
+                "min_left_cls": {
+                    "type": "integer",
+                    "description": "Minimum cluser size in the left batch",
+                    "default": 2
+                },
+                "consensus_period": {
+                    "type": "integer",
+                    "description": "Consensus period (-1 means no consensus)",
+                    "default": 500
+                },
+                "consensus_minimum": {
+                    "type": "integer",
+                    "description": "Minimum consensus sample size:",
+                    "default": 50
+                },
+                "consensus_maximum": {
+                    "type": "integer",
+                    "description": "Maximum consensus sample size",
+                    "default": -150
+                },
+                "min_shared": {
+                    "type": "integer",
+                    "description": "Minimum number of minimizers shared between read and cluster",
+                    "default": 5
+                },
+                "min_qual": {
+                    "description": "Minimum average quality value",
+                    "type": "number",
+                    "default": 7.0
+                },
+                "mapped_threshold": {
+                    "description": "Minimum mapped fraction of read to be included in cluster",
+                    "type": "number",
+                    "default": 0.65
+                },
+                "aligned_threshold": {
+                    "tpye": "number",
+                    "description": "Minimum aligned fraction of read to be included in cluster",
+                    "default": 0.2
+                },
+                "min_fraction": {
+                    "type": "number",
+                    "description": "Minimum fraction of minimizers shared compared to best hit, in order to continue mapping",
+                    "default": 0.8
+                },
+                "min_prob_no_hits" : {
+                    "type": "number",
+                    "description": "Minimum probability for i consecutive minimizers to be different between read and representative",
+                    "default": 0.2
                 }
             },
             "required": [
-                "fastq",
+                "fastq"
-                "ref_genome"
             ]
         },
         "meta_data": {

reference_assembly.nf

@@ -0,0 +1,48 @@
+process map_reads{
+    /*
+    Map reads to reference using minimap2.
+    Filter reads by mapping quality.
+    Filter reads where length of poly(A) > max_poly_run at either ends of the read (defined by poly_context)
+    */
+    label "isoforms"
+    cpus params.threads
+
+    input:
+       path index
+       path reference
+       tuple val(sample_id), path (fastq_reads)
+
+    output:
+       tuple val(sample_id), path("${sample_id}_reads_aln_sorted.bam"), emit: bam
+       tuple val(sample_id), path("${sample_id}_read_aln_stats.tsv"), emit: stats
+    script:
+        def ContextFilter = """AlnContext: { Ref: "${reference}", LeftShift: -${params.poly_context}, RightShift: ${params.poly_context},
+           RegexEnd: "[Aa]{${params.max_poly_run},}",
+           Stranded: True,Invert: True, Tsv: "internal_priming_fail.tsv"} """
+       """
+       seqkit fq2fa ${fastq_reads} -o "reads.fa";
+       minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} ${index} "reads.fa"\
+       | samtools view -q ${params.minimum_mapping_quality} -F 2304 -Sb -\
+       | seqkit bam -j ${params.threads} -x -T '${ContextFilter}' -\
+       | samtools sort -@ ${params.threads} -o "${sample_id}_reads_aln_sorted.bam" - \
+       | ((seqkit bam -s -j ${params.threads} - 2>&1)  | tee ${sample_id}_read_aln_stats.tsv ) || true
+
+       if [[ -s "internal_priming_fail.tsv" ]];
+       then
+           tail -n +2 "internal_priming_fail.tsv" | awk '{{print ">" \$1 "\\n" \$4 }}' - > "context_internal_priming_fail_start.fasta"
+           tail -n +2 "internal_priming_fail.tsv" | awk '{{print ">" \$1 "\\n" \$6 }}' - > "context_internal_priming_fail_end.fasta"
+       fi
+       """
+}
+
+workflow reference_assembly {
+    take:
+       index
+       reference
+       fastq_reads
+    main:
+        map_reads(index, reference, fastq_reads)
+    emit:
+       bam = map_reads.out.bam
+       stats = map_reads.out.stats
+}

run_evaluation_dmel.sh

@@ -1,40 +0,0 @@
-#!/bin/bash
-
-# Usage: ./run_evaluation_dmel.sh pathto/outputdir
-
-# See the isONcorrect paper https://www.nature.com/articles/s41467-020-20340-8 where this dataset is described
-
-OUTDIR=$1;
-
-FASTQ_URL="http://ftp.sra.ebi.ac.uk/vol1/fastq/ERR358/005/ERR3588905/ERR3588905_1.fastq.gz"
-REF_URL="http://ftp.ensembl.org/pub/release-99/fasta/drosophila_melanogaster/dna/Drosophila_melanogaster.BDGP6.28.dna.toplevel.fa.gz"
-GFF_URL="http://ftp.ensembl.org/pub/release-99/gff3/drosophila_melanogaster/Drosophila_melanogaster.BDGP6.28.99.gff3.gz"
-
-RESULTS_DIR="$OUTDIR/results"
-DATA_DIR="$OUTDIR/data"
-READS_DIR="$DATA_DIR/reads"
-FASTQ="$READS_DIR/ERR3588905_1.fastq"
-REF="$DATA_DIR/Drosophila_melanogaster.BDGP6.28.dna.toplevel.fa"
-GFF="$DATA_DIR/Drosophila_melanogaster.BDGP6.28.99.gff3"
-
-echo $READS_DIR;
-
-rm -fr $OUT_DIR/results
-mkdir -p $OUTDIR/data
-
-if [ ! -f $REF ];
-then (cd $DATA_DIR; curl -L -C - -O $REF_URL); gzip -d ${REF}.gz
-fi
-
-if [ ! -f $GFF ]
-then
-	(cd $DATA_DIR; curl -L -C - -O $GFF_URL); gzip -d ${GFF}.gz
-fi
-
-if [ ! -f $FASTQ ];
-then (cd $READS_DIR; curl -L -C - -O $FASTQ_URL); gzip -d ${FASTQ}.gz
-fi
-
-nextflow run ../wf-isoforms --fastq $READS_DIR  \
---reference_genome $REF --annotation $GFF -profile conda --out_dir $OUTDIR \
--w $OUTDIR/workspace -resume

tests.sh

@@ -1,53 +0,0 @@
-#!/usr/bin/env bash
-
-# A couple of simple tests with different combinations of CLI options
-
-if [[ "$#" -ne 1 ]]; then
-    echo "Please supply path to out_dir"
-    exit 1
-fi
-
-SCRIPT_DIR="$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )"
-cd $SCRIPT_DIR;
-
-singledir="test_data/fastq"
-multisampledir="test_data/demultiplexed_fastq"
-
-echo $singledir
-echo $multisampledir
-results=()
-
-# Test1 single sample
-OUTPUT=$1/test1;
-nextflow run . --fastq $singledir --ref_genome test_data/SIRV_150601a.fasta \
---ref_annotation test_data/SIRV_isofroms.gtf -profile conda --out_dir ${OUTPUT} -w ${OUTPUT}/workspace -resume;
-
-results+=("test1: $?")
-
-# Test2 multiple samples demultiplexed
-OUTPUT=$1/test2;
-nextflow run . --fastq $multisampledir --ref_genome test_data/SIRV_150601a.fasta \
---ref_annotation test_data/SIRV_isofroms.gtf -profile conda --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
---sample_sheet test_data/sample_sheet -resume;
-
-results+=("test2: $?")
-
-# Test3 single sample. No reference annotation
-OUTPUT=$1/test3;
-nextflow run . --fastq $singledir --ref_genome test_data/SIRV_150601a.fasta \
- -profile conda --out_dir ${OUTPUT} -w ${OUTPUT}/workspace -resume;
-
-results+=("test3: $?")
-
-# Test4 single sample. force split_bam to make multiple alignment bundles
-OUTPUT=$1/test4;
-nextflow run . --fastq $singledir  --ref_genome test_data/SIRV_150601a.fasta \
---ref_annotation test_data/SIRV_isofroms.gtf -profile conda --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \
---bundle_min_reads 5 -resume;
-
-results+=("test4: $?")
-
-echo "Exit status codes for each test"
-for value in "${results[@]}"; do
-     echo "${value}"
-done