Merge branch 'CW-1044' into 'dev'

demo data de

See merge request epi2melabs/workflows/wf-transcriptomes!76

.gitlab-ci.yml

@@ -46,10 +46,10 @@ docker-run:
               NF_IGNORE_PROCESSES: preprocess_reads,merge_transcriptomes
         - if: $MATRIX_NAME == "differential_expression"
           variables:
-              NF_BEFORE_SCRIPT: wget -O differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_differential_expression.tar.gz && tar -xzvf  differential_expression.tar.gz
+              NF_BEFORE_SCRIPT: tar -xzvf test_data/differential_expression.tar.gz
-              NF_WORKFLOW_OPTS: "--fastq  differential_expression_dataset/fastq \
+              NF_WORKFLOW_OPTS: "--fastq  differential_expression/differential_expression_fastq \
                   --de_analysis \
-                  --ref_genome differential_expression_dataset/hg38_chr20.fa \
+                  --ref_genome differential_expression/hg38_chr20.fa \
-                --ref_annotation differential_expression_dataset/gencode.v22.annotation.chr20.gtf \
+                --ref_annotation differential_expression/gencode.v22.annotation.chr20.gtf \
                 --direct_rna"
               NF_IGNORE_PROCESSES: preprocess_reads,merge_transcriptomes

.pre-commit-config.yaml

@@ -22,8 +22,13 @@ repos:
     hooks:
       - id: flake8
         additional_dependencies:
-          - flake8-import-order==0.18.1
+          - flake8-rst-docstrings
-          - flake8-docstrings==1.6.0
+          - flake8-docstrings
-          - flake8-rst-docstrings==0.2.5
+          - flake8-import-order
-          - flake8-forbid-visual-indent==0.0.2
+          - flake8-forbid-visual-indent
+          - pep8-naming
+          - flake8-no-types
+          - flake8-builtins
+          - flake8-absolute-import
+          - flake8-print
         entry: flake8 bin --import-order-style google --statistics

CHANGELOG.md

@@ -4,13 +4,20 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [unreleased]
+### Updated
+- Removed sanitize option
+### Added
+- Demo differential expression data in repository.
+- Improved DE explanation in docs
+
 ## [v0.1.5]
 ### Added
 - Differential transcript and gene expression subworkflow
 
 ## [v0.1.4]
 ### Added
-- JAFFAL fusion detectoion subworkflow
+- JAFFAL fusion detection subworkflow
 ### Changed
 - Args parser for fastqingress
 - Set out_dir option type to ensure output is written to correct directory on Windows

README.md

@@ -41,12 +41,38 @@ Fusion gene detection is performed using [JAFFA](https://github.com/Oshlack/JAFF
 with ONT long reads. 
 
 ### Differential expression analysis
-* Differential expression is done using the transcripts output by the workflow.
+
-* A non redundant transcriptome is found using the merge function in [stringtie](http://ccb.jhu.edu/software/stringtie).
+Differential gene expression (DGE) and differential transcript usage (DTU) analyses aim to identify genes and/or transcripts that show statistically altered expression patterns in a studied biological system. The results of the differential analyses are presented in a quantitative format and therefore the degree of change (up or down regulation) between experimental conditions can be calculated for each gene identified.
-* The reads are then aligned to the transcriptome using minimap2 in a splice-aware manner.
+
-* [salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification.
+These differential analyses work by taking a “snapshot” of mRNA abundance and calculating the relative levels of transcripts and isoforms. In this context, expression corresponds to the number of messenger RNAs (mRNA) measured from each gene isoform within the organism / tissue / culture being investigated. In order to determine expression levels across the whole genome, sequence data specifically targeting the mRNA molecules can be generated.
-* R packages [edgeR](https://bioconductor.org/packages/release/bioc/html/edgeR.html) and [stageR](https://bioconductor.org/packages/release/bioc/html/stageR.html) are used for differential expression analysis.
+
-* [DEXSeq](https://bioconductor.org/packages/release/bioc/html/DEXSeq.html) is then used for differential transcript usage analysis.
+Oxford Nanopore Technologies provides a number of sequencing solutions to allow users to generate the required snapshot of gene expression. This can be achieved by both sequencing the mRNA directly, or via a complementary DNA (cDNA) proxy. In contrast to short read sequencing technologies, entire mRNA transcripts can be captured as single reads. The example data provided with this tutorial is from a study based on the PCR-cDNA kit. This is a robust choice for performing differential transcript usage studies. This kit is suitable for preparation of sequence libraries from low mRNA input quantities. The cDNA population is enriched through PCR with low bias; an important prerequisite for the subsequent statistical analysis.
+
+[Workflow-transcriptomes](https://github.com/epi2me-labs/wf-transcriptomes) includes a subworkflow for DGE and DTU. The first step involves using either a reference alignment or _de novo_ assembly approach to create a set of mRNA sequences per sample. These are merged into a non-redundant transcriptome using [stringtie merge](http://ccb.jhu.edu/software/stringtie). The reads are then aligned to the transcriptome using minimap2 in a splice-aware manner. [Salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification, giving per transcript counts and then the following R packages are used for analysis.
+
+### Pre-filtering of quantitative data using DRIMSeq
+DRIMSeq (Nowicka and Robinson (2016)) is used to filter the transcript count data from the salmon analysis. The filter step will be used to select for genes and transcripts that satisfy rules for the number of samples in which a gene or transcript must be observed and minimum threshold levels for the number of observed reads. The parameters used for filtering are defined in the config.yaml file. The default parameters defined for this analysis include
+* min_samps_gene_expr = 3 - a transcript must be mapped to a gene in at least this minimum number of samples for the gene be included in the analysis
+*	min_samps_feature_expr = 1 - a transcript must be mapped to an isoform in at least this this minimum number of samples for the gene isoform to be included in the analysis
+*	min_gene_expr = 10 - the minimum number of total mapped sequence reads for a gene to be considered expressed
+*	min_feature_expr = 3 - the minimum number of total mapped sequence reads for a gene isoform to be considered
+
+### edgeR based differential expression analysis
++A statistical analysis is first performed using edgeR (Robinson, McCarthy, and Smyth (2010), McCarthy et al. (2012)) to identify the subset of differentially expressed genes. The filtered list of gene counts is used as input. A normalisation factor is calculated for each sequence library (using the default TMM method - please see McCarthy et al. (2012) for further details). The defined experimental design is used to calculate estimates of dispersion for each of the gene features. Statistical tests are calculated using the contrasts defined in the experimental design. The differentially expressed genes are corrected for false discovery (fdr) using the method of Benjamini & Hochberg (Benjamini and Hochberg (1995))
+
+### Differential transcript usage using DEXSeq
+Differential transcript usage analysis is performed using the R DEXSeq package (Reyes et al. (2013)). Similar to the edgeR package, DEXSeq estimates the variance between the biological replicates and applies generalised linear models for the statistical testing. The key difference is that the DEXSeq method looks for differences at the exon count level. DEXSeq uses the filtered transcript count data prepared earlier in this analysis. 
+
+### StageR stage-wise analysis of DGE and DTU
+The final component of this isoform analysis is a stage-wise statistical test using the R software package `stageR` (Van den Berge and Clement (2018)). stageR uses (1) the raw p-values for DTU from the DEXSeq analysis in the previous section and (2) a false-discovery corrected set of p-values from testing whether individual genes contain at least one exon showing DTU. A hierarchical two-stage statistical testing evaluates the set of genes for DTU.
+
+## Running the workflow
+For the differential expression analysis section you should have at least 3 repeats for each sample. 
+Your fastq data will need to be organised in to 6 directories that represent 3 repeats for each condition. You may also need to provide a condition sheet. 
+
+
+## Analysis 
+Differential gene expression is sensitive to the input data quantity and quality.  There should be equivalence between samples in the number of sequence reads, mapped reads and quality scores. The sequence and alignment summary plots in the report can be used to assess these metrics. There is also a table that shows the transcript per million(TPM) calculated from the salmon counts. TPM normalizes the data for gene length and then sequencing depth, and makes it easier to compare across samples compared to counts.
 
 ### Workflow inputs
 - Directory containing cDNA/direct RNA reads. Or a directory containing subdirectories each with reads from different samples
@@ -178,20 +204,37 @@ Differential Expression requires at least 2 replicates of each sample to compare
 
 **Example workflow for differential expression transcript assembly**
 
-Download differential expression data set 
+#### Condition sheet
+The condition sheet should be a .tsv with two columns.
+- The sample column will need to match the 6 directories in the input fastq directory.
+- The condition column will need to contain one of two keys to indicate the two samples being compared.
 
-`wget -O differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_differential_expression.tar.gz && tar -xzvf  differential_expression.tar.gz`
+In the default `condition_sheet.tsv` available in the test_data directory we have used the following.
 
-Run the cmd 
+eg. condition_sheet.tsv
+```
+sample,condition
+barcode01,untreated
+barcode02,untreated
+barcode03,untreated
+barcode04,treated
+barcode05,treated
+barcode06,treated
+```
+
+You will also need to provide a reference genome and a reference annotation file.
+Here is an example cmd to run the workflow using the test_data provided. 
 
 ```
 OUTPUT=~/output;
-nexflow run epi2me-labs/wf-transcriptomes --fastq  differential_expression_dataset/fastq --de_analysis \
+nexflow run epi2me-labs/wf-transcriptomes --fastq test_data/differential_expression_fastq \
---ref_genome differential_expression_dataset/hg38_chr20.fa \
+  --de_analysis \
---ref_annotation differential_expression_dataset/gencode.v22.annotation.chr20.gtf \
+  --ref_genome test_data/hg38_chr20.fa \
---direct_rna
+  --ref_annotation test_data/gencode.v22.annotation.chr20.gtf \
+  --direct_rna
 ```
 
+
 ## Workflow outputs
 * an HTML report document detailing the primary findings of the workflow.
 * for each sample:
@@ -201,15 +244,25 @@ nexflow run epi2me-labs/wf-transcriptomes --fastq  differential_expression_datas
   * merged_transcritptome.fas - annotated, assembled transcriptome
   * [jaffal](https://github.com/Oshlack/JAFFA) ooutput directories
   
-
 ### Fusion detection outputs
 in `${out_dir}/jaffal_output_${sample_id}` you will find:
 * jaffa_results.csv - the csv results summary file 
 * jaffa_results.fasta - fusion transcritpt sequences
 
 ### Differential Expression outputs
-* dtu_plots.pdf - a pdf with differntial transcript usage plots
+* `de_analysis/results_dge.tsv` and `de_analysis/results_dge.pdf`- results of `edgeR` differential gene expression analysis.
-## Useful links
+* `de_analysis/results_dtu_gene.tsv`, `de_analysis/results_dtu_transcript.tsv` and `de_analysis/results_dtu.pdf` - results of differential transcript usage by `DEXSeq`.
+* `de_analysis/results_dtu_stageR.tsv` - results of the `stageR` analysis of the `DEXSeq` output.
+* `de_analysis/dtu_plots.pdf` - DTU results plot based on the `stageR` results and filtered counts.
+
+### References
+
+* Benjamini, Yoav, and Yosef Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society. Series B (Methodological) 57 (1): 289–300. http://www.jstor.org/stable/2346101.
+* McCarthy, Davis J., Chen, Yunshun, Smyth, and Gordon K. 2012. “Differential Expression Analysis of Multifactor Rna-Seq Experiments with Respect to Biological Variation.” Nucleic Acids Research 40 (10): 4288–97.
+* Nowicka, Malgorzata, and Mark D. Robinson. 2016. “DRIMSeq: A Dirichlet-Multinomial Framework for Multivariate Count Outcomes in Genomics [Version 2; Referees: 2 Approved].” F1000Research 5 (1356). https://doi.org/10.12688/f1000research.8900.2.
+* Patro, Robert, Geet Duggal, Michael I Love, Rafael A Irizarry, and Carl Kingsford. 2017. “Salmon Provides Fast and Bias-Aware Quantification of Transcript Expression.” Nature Methods 14 (March). https://doi.org/10.1038/nmeth.4197.
+* Robinson, Mark D, Davis J McCarthy, and Gordon K Smyth. 2010. “EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1): 139–40.
+* Love, Michael I., et al. Swimming Downstream: Statistical Analysis of Differential Transcript Usage Following Salmon Quantification. 7:952, F1000Research, 14 Sept. 2018. f1000research.com, https://f1000research.com/articles/7-952## Useful links
 
 * [nextflow](https://www.nextflow.io/)
 * [docker](https://www.docker.com/products/docker-desktop)

bin/compute_cluster_quality.py

@@ -9,6 +9,7 @@ import argparse
 from collections import defaultdict
 import math
 from pathlib import Path
+import sys
 
 import matplotlib
 from matplotlib import pyplot as plt
@@ -61,7 +62,6 @@ def parse_true_clusters(ref_file):
         # 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
@@ -89,7 +89,6 @@ def parse_true_clusters(ref_file):
         # 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
@@ -112,8 +111,8 @@ def parse_true_clusters_simulated(ref_file):
     return classes
 
 
-def compute_V_measure(clusters, classes):
+def compute_v_measure(clusters, classes):
-    """compute_V_measure."""
+    """compute_v_measure."""
     class_list, cluster_list = [], []
     # not_found_id = 1000000
     clustered_but_unaligned = 0
@@ -122,7 +121,6 @@ def compute_V_measure(clusters, 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
@@ -139,22 +137,22 @@ def compute_V_measure(clusters, classes):
     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))
+    sys.stdout("Not included in clustering but aligned:", len(not_clustered))
-    print(
+    sys.stdout(
-        "V:",
+        "v:",
         v_score,
         "Completeness:",
         compl_score,
         "Homogeneity:",
         homog_score)
-    print(
+    sys.stdout(
         "Nr reads clustered but unaligned "
-        "(i.e., no class and excluded from V-measure): ",
+        "(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):
+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
@@ -188,19 +186,19 @@ def compute_V_measure_non_singleton_classes(clusters, classes):
     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(
+    sys.stdout(
-        "NONTRIVIAL CLASSES: V:",
+        "NONTRIvIAL CLASSES: v:",
         v_score,
         "Completeness:",
         compl_score,
         "Homogeneity:",
         homog_score)
-    print("NUMBER OF CLASSES (FILTERED):", len(
+    sys.stdout("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):
+def compute_v_measure_non_singletons(clusters, classes):
     """V measure for non-singletons."""
     cluster_dict = {}
     for read_acc, cl_id in clusters.items():
@@ -225,46 +223,46 @@ def compute_V_measure_non_singletons(clusters, classes):
             class_list.append(classes[read])
             cluster_list.append(clusters[read])
         else:
-            # print("Read was clustered but unaligned:", read)
+            # sys.stdout("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(
+    sys.stdout(
-        "NONTRIVIAL CLUSTERS: V:",
+        "NONTRIvIAL CLUSTERS: v:",
         v_score,
         "Completeness:",
         compl_score,
         "Homogeneity:",
         homog_score)
-    print(
+    sys.stdout(
-        "NONTRIVIAL CLUSTERS: Nr reads clustered but unaligned "
+        "NONTRIvIAL CLUSTERS: Nr reads clustered but unaligned "
-        "(i.e., no class and excluded from V-veasure): ",
+        "(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):
+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 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:
+    if not n:
         return None
-    k = (len(N) - 1) * percent
+    k = (len(n) - 1) * percent
     f = math.floor(k)
     c = math.ceil(k)
     if f == c:
-        return key(N[int(k)])
+        return key(n[int(k)])
-    d0 = key(N[int(f)]) * (c - k)
+    d0 = key(n[int(f)]) * (c - k)
-    d1 = key(N[int(c)]) * (k - f)
+    d1 = key(n[int(c)]) * (k - f)
     return d0 + d1
 
 # end of http://code.activestate.com/recipes/511478/ }}}
@@ -364,11 +362,11 @@ def get_cluster_information(clusters, classes):
         else:
             clustered_classes[class_id] += 1
 
-    print("UNCLUSTERED:", "Tot classes:", len(not_clustered_classes))
+    sys.stdout("UNCLUSTERED:", "Tot classes:", len(not_clustered_classes))
-    print("CLUSTERED:", "Tot classes:", len(clustered_classes))
+    sys.stdout("CLUSTERED:", "Tot classes:", len(clustered_classes))
-    print("MIXED:", "Tot classes containing both:", len(
+    sys.stdout("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)
+    sys.stdout("Total number of classes (unique gene ID):", total_nr_classes)
     return (
         total_nr_classes - len(singleton_classes),
         len(singleton_classes),
@@ -411,7 +409,7 @@ def main(args):
         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)
+        compute_v_measure(clusters, classes)
 
     (
         nr_non_singleton_classes, singleton_classes, min_class_size,
@@ -451,30 +449,31 @@ def main(args):
             e_class_size,
             n50_class_size))
 
-    # Reads_nontrivially_clustered_(%), Singletons_(%),
+    # reads_nontrivially_clustered_(%), Singletons_(%),
-    # Reads_Nontrivially_clustered_but_unaligned, V, c,h ,V_nt, c_nt,h_nt,
+    # 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(
+    reads_nontrivially_clustered_percent = round(
         100 * (float(tot_nr_reads - singleton_clusters) / tot_nr_reads), 1)
-    # round(1.0 - Reads_nontrivially_clustered_percent, 2)
+    # round(1.0 - reads_nontrivially_clustered_percent, 2)
-    Reads_Nontrivially_clustered_but_unaligned = \
+    reads_nontrivially_clustered_but_unaligned = \
         unaligned_but_nontrivially_clustered
-    V, c, h = round(v_score, 3), round(compl_score, 3), round(homog_score, 3)
+    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))
+    sys.stdout(
+        "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",
+            "v",
             "c",
             "h",
             "ARI",
-            "Reads_nontrivially_clustered_percent",
+            "reads_nontrivially_clustered_percent",
-            "Reads_Nontrivially_clustered_but_unaligned",
+            "reads_nontrivially_clustered_but_unaligned",
             "non_singleton_clusters",
             "singleton_clusters",
             "upper_75_cluster_size",
@@ -483,12 +482,12 @@ def main(args):
             "n50_cluster_size"))
     outfile.write(
         "{0},{1},{2},{3},{4},{5},{6},{7},{8},{9},{10},{11}\n".format(
-            V,
+            v,
             c,
             h,
             ari,
-            Reads_nontrivially_clustered_percent,
+            reads_nontrivially_clustered_percent,
-            Reads_Nontrivially_clustered_but_unaligned,
+            reads_nontrivially_clustered_but_unaligned,
             non_singleton_clusters,
             singleton_clusters,
             upper_75_cluster_size,
@@ -508,19 +507,19 @@ def main(args):
     dfc = pd.DataFrame(
         {
             'Statistic': [
-                'V-measure',
+                'v-measure',
                 'ARI',
                 'Completeness',
                 'Homogeneity'],
-            'Value': [
+            'value': [
-                V,
+                v,
                 ari,
                 c,
                 h]}).set_index('Statistic')
     dfn = pd.DataFrame(
         {'Statistic': ['NonSingleton',
                        'Singletons'],
-         'Value': [non_singleton_clusters,
+         'value': [non_singleton_clusters,
                    singleton_clusters]}).set_index('Statistic')
     dfs = pd.DataFrame(
         {
@@ -529,14 +528,14 @@ def main(args):
                 'Upper75ClassSize',
                 'MedianClsSize',
                 'MedianClassSize'],
-            'Value': [
+            '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]
+            'value': [n50_cluster_size, n50_class_size]
          }).set_index('Statistic')
 
     rdo = Path(args.raw_data_out)
@@ -648,6 +647,6 @@ if __name__ == '__main__':
         help='dir to save raw data for plotting')
     args = parser.parse_args()
 
-    print("------------------------------------------------------------")
+    sys.stdout("------------------------------------------------------------")
     main(args)
-    print("------------------------------------------------------------")
+    sys.stdout("------------------------------------------------------------")

bin/de_plots.py

@@ -36,7 +36,7 @@ def number_of_alignments(df, field_name):
 
 def create_summary_table(df):
     """Create summary table."""
-    all = number_of_alignments(df, "Read mappings")
+    all_aln = number_of_alignments(df, "Read mappings")
     primary = number_of_alignments(df.loc[df['Type'] == 'Primary'], "Primary")
     secondary = number_of_alignments(
         df.loc[df['Type'] == 'Secondary'], "Secondary")
@@ -47,7 +47,7 @@ def create_summary_table(df):
     avg_mapq = df.loc[df['Type'] == 'Primary'].groupby(
         'fname').agg(**{"Median MAPQ": ('MapQual', 'median'), }).transpose()
     return pd.concat([
-        all, primary, secondary, supplementary,
+        all_aln, primary, secondary, supplementary,
         avg_acc, avg_mapq])
 
 

bin/generate_tracking_summary.py

@@ -44,14 +44,15 @@ def generate_tracking_summary(tracking_file, output_dir, annotations=None):
         # write a separate table for each class
         for class_code, table in tracking.groupby('class'):
             if not write_empty_tsvs and table.empty:
-                print("Skipping: No transcripts found for: {}".format(
+                sys.stdout("Skipping: No transcripts found for: {}".format(
                     class_code))
                 continue
 
             path = tracking_file + ".{}.tsv".format(class_code)
             table.to_csv(path)
     else:
-        print("Skipping classification summary as no annotation provided.")
+        sys.stdout(
+            "Skipping classification summary as no annotation provided.")
 
 
 def main(args):

bin/report.py

@@ -6,6 +6,7 @@ from collections import Counter, defaultdict, OrderedDict
 import math
 import os
 from pathlib import Path
+import sys
 
 from aplanat import bars, hist
 from aplanat.components import simple as scomponents
@@ -289,7 +290,7 @@ def grouped_bar(df, title="", tilted_xlabs=False):
     return p
 
 
-def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
+def gff_compare_plots(report, gffcompare_outdirs, sample_ids):
     """Create various sections and plots in a WfReport.
 
     :param report: aplanat WFReport
@@ -387,7 +388,6 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
 
     tracking_dfs = []
 
-    print(gffcompare_outdirs)
     track_files = [x / 'str_merged.tracking' for x in gffcompare_outdirs]
 
     df_tracking = load_sample_data(
@@ -423,7 +423,7 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
         })
 
         cols = [TableColumn(
-            field=Ci, title=Ci, width=100) for Ci in tracking.columns]
+            field=ci, title=ci, width=100) for ci in tracking.columns]
 
         track_table = DataTable(
             columns=cols, source=ColumnDataSource(tracking),
@@ -473,7 +473,7 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
     try:
         tmap_files = [next(x.glob('*.tmap')) for x in gffcompare_outdirs]
     except StopIteration:
-        print("Cannot find .tmap files in {}".format(gffcompare_outdirs))
+        sys.stderr("Cannot find .tmap files in {}".format(gffcompare_outdirs))
         return
 
     df_tmap = load_sample_data(tmap_files, sample_ids)
@@ -621,7 +621,7 @@ def transcript_table(report, df_tmaps, max_rows):
             'FPKM', 'qry_gene_id', 'major_iso_id', 'ref_match_len', 'TPM'])
 
     if len(df) == 0:
-        print("No transcripts found")
+        sys.stderr("No transcripts found")
         section.markdown("No transcripts found")
         return
 
@@ -739,7 +739,7 @@ def transcriptome_summary(report, gffs, sample_ids, denovo=False):
 
         df_sum.columns = [' ', 'count']
         cols = [TableColumn(
-            field=Ci, title=Ci, width=80) for Ci in df_sum.columns]
+            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)

bin/run_isonclust2.py

@@ -11,14 +11,14 @@ import subprocess as sub
 class Node:
     """Node."""
 
-    def __init__(self, Id, File, Left, Right, Parent, Level):
+    def __init__(self, node_id, file_, left, right, parent, level):
         """Set node attaributes."""
-        self.Id = Id
+        self.Id = node_id
-        self.File = File
+        self.File = file_
-        self.Left = Left
+        self.Left = left
-        self.Right = Right
+        self.Right = right
-        self.Parent = Parent
+        self.Parent = parent
-        self.Level = Level
+        self.Level = level
         self.Done = False
         self.RightSide = False
 
@@ -45,35 +45,35 @@ def grouper(n, iterable, fillvalue=None):
 
 def build_job_tree():
     """Build a job tree of nodes."""
-    JOB_TREE = OrderedDict()
+    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 = OrderedDict()
-    LEVELS[0] = []
+    levels[0] = []
-    for Id, bf in sorted(zip(batch_ids, batches), key=lambda x: x[0]):
+    for k, bf in sorted(zip(batch_ids, batches), key=lambda x: x[0]):
         n = Node(
-            Id,
+            k,
-            "clusters/isONcluster_{}.cer".format(Id),
+            "clusters/isONcluster_{}.cer".format(k),
             None,
             None,
             None,
             0)
         n.Done = True
-        JOB_TREE[Id] = n
+        job_tree[k] = n
-        LEVELS[0].append(n)
+        levels[0].append(n)
     level = 0
-    max_id = LEVELS[0][-1].Id
+    max_id = levels[0][-1].k
-    while len(LEVELS[level]) != 1:  # Final level will be link
+    while len(levels[level]) != 1:  # Final level will be link
         next_level = level + 1
-        LEVELS[next_level] = []
+        levels[next_level] = []
-        for l_, r in grouper(2, LEVELS[level]):
+        for l_, r in grouper(2, levels[level]):
             if r is None:  # End of a level
-                LEVELS[level].pop()  # remove last node?
+                levels[level].pop()  # remove last node?
                 l_.Level += 1  # ncrement level
-                LEVELS[next_level].append(l_)  # Add the left to the next 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)
@@ -81,13 +81,13 @@ def build_job_tree():
             l_.Parent = new_node
             r.Parent = new_node
             r.RightSide = True
-            LEVELS[next_level].append(new_node)
+            levels[next_level].append(new_node)
-            JOB_TREE[max_id] = new_node
+            job_tree[max_id] = new_node
         level = next_level
-    ROOT = JOB_TREE[len(JOB_TREE) - 1].Id
+    root = job_tree[len(job_tree) - 1].k
-    JOB_TREE[ROOT].RightSide = True
+    job_tree[root].RightSide = True
 
-    return JOB_TREE, LEVELS
+    return job_tree, levels
 
 
 def main():
@@ -121,7 +121,7 @@ def main():
         sub.call(cmd, shell=True)
     sub.call((
         "ln -s `realpath clusters/isONcluster_{}.cer` "
-        "isONcluster_ROOT.cer".format(n.Id)), shell=True)
+        "isONcluster_root.cer".format(n.Id)), shell=True)
 
 
 if __name__ == '__main__':

docs/intro.md

@@ -33,12 +33,38 @@ Fusion gene detection is performed using [JAFFA](https://github.com/Oshlack/JAFF
 with ONT long reads. 
 
 ### Differential expression analysis
-* Differential expression is done using the transcripts output by the workflow.
+
-* A non redundant transcriptome is found using the merge function in [stringtie](http://ccb.jhu.edu/software/stringtie).
+Differential gene expression (DGE) and differential transcript usage (DTU) analyses aim to identify genes and/or transcripts that show statistically altered expression patterns in a studied biological system. The results of the differential analyses are presented in a quantitative format and therefore the degree of change (up or down regulation) between experimental conditions can be calculated for each gene identified.
-* The reads are then aligned to the transcriptome using minimap2 in a splice-aware manner.
+
-* [salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification.
+These differential analyses work by taking a “snapshot” of mRNA abundance and calculating the relative levels of transcripts and isoforms. In this context, expression corresponds to the number of messenger RNAs (mRNA) measured from each gene isoform within the organism / tissue / culture being investigated. In order to determine expression levels across the whole genome, sequence data specifically targeting the mRNA molecules can be generated.
-* R packages [edgeR](https://bioconductor.org/packages/release/bioc/html/edgeR.html) and [stageR](https://bioconductor.org/packages/release/bioc/html/stageR.html) are used for differential expression analysis.
+
-* [DEXSeq](https://bioconductor.org/packages/release/bioc/html/DEXSeq.html) is then used for differential transcript usage analysis.
+Oxford Nanopore Technologies provides a number of sequencing solutions to allow users to generate the required snapshot of gene expression. This can be achieved by both sequencing the mRNA directly, or via a complementary DNA (cDNA) proxy. In contrast to short read sequencing technologies, entire mRNA transcripts can be captured as single reads. The example data provided with this tutorial is from a study based on the PCR-cDNA kit. This is a robust choice for performing differential transcript usage studies. This kit is suitable for preparation of sequence libraries from low mRNA input quantities. The cDNA population is enriched through PCR with low bias; an important prerequisite for the subsequent statistical analysis.
+
+[Workflow-transcriptomes](https://github.com/epi2me-labs/wf-transcriptomes) includes a subworkflow for DGE and DTU. The first step involves using either a reference alignment or _de novo_ assembly approach to create a set of mRNA sequences per sample. These are merged into a non-redundant transcriptome using [stringtie merge](http://ccb.jhu.edu/software/stringtie). The reads are then aligned to the transcriptome using minimap2 in a splice-aware manner. [Salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification, giving per transcript counts and then the following R packages are used for analysis.
+
+### Pre-filtering of quantitative data using DRIMSeq
+DRIMSeq (Nowicka and Robinson (2016)) is used to filter the transcript count data from the salmon analysis. The filter step will be used to select for genes and transcripts that satisfy rules for the number of samples in which a gene or transcript must be observed and minimum threshold levels for the number of observed reads. The parameters used for filtering are defined in the config.yaml file. The default parameters defined for this analysis include
+* min_samps_gene_expr = 3 - a transcript must be mapped to a gene in at least this minimum number of samples for the gene be included in the analysis
+*	min_samps_feature_expr = 1 - a transcript must be mapped to an isoform in at least this this minimum number of samples for the gene isoform to be included in the analysis
+*	min_gene_expr = 10 - the minimum number of total mapped sequence reads for a gene to be considered expressed
+*	min_feature_expr = 3 - the minimum number of total mapped sequence reads for a gene isoform to be considered
+
+### edgeR based differential expression analysis
++A statistical analysis is first performed using edgeR (Robinson, McCarthy, and Smyth (2010), McCarthy et al. (2012)) to identify the subset of differentially expressed genes. The filtered list of gene counts is used as input. A normalisation factor is calculated for each sequence library (using the default TMM method - please see McCarthy et al. (2012) for further details). The defined experimental design is used to calculate estimates of dispersion for each of the gene features. Statistical tests are calculated using the contrasts defined in the experimental design. The differentially expressed genes are corrected for false discovery (fdr) using the method of Benjamini & Hochberg (Benjamini and Hochberg (1995))
+
+### Differential transcript usage using DEXSeq
+Differential transcript usage analysis is performed using the R DEXSeq package (Reyes et al. (2013)). Similar to the edgeR package, DEXSeq estimates the variance between the biological replicates and applies generalised linear models for the statistical testing. The key difference is that the DEXSeq method looks for differences at the exon count level. DEXSeq uses the filtered transcript count data prepared earlier in this analysis. 
+
+### StageR stage-wise analysis of DGE and DTU
+The final component of this isoform analysis is a stage-wise statistical test using the R software package `stageR` (Van den Berge and Clement (2018)). stageR uses (1) the raw p-values for DTU from the DEXSeq analysis in the previous section and (2) a false-discovery corrected set of p-values from testing whether individual genes contain at least one exon showing DTU. A hierarchical two-stage statistical testing evaluates the set of genes for DTU.
+
+## Running the workflow
+For the differential expression analysis section you should have at least 3 repeats for each sample. 
+Your fastq data will need to be organised in to 6 directories that represent 3 repeats for each condition. You may also need to provide a condition sheet. 
+
+
+## Analysis 
+Differential gene expression is sensitive to the input data quantity and quality.  There should be equivalence between samples in the number of sequence reads, mapped reads and quality scores. The sequence and alignment summary plots in the report can be used to assess these metrics. There is also a table that shows the transcript per million(TPM) calculated from the salmon counts. TPM normalizes the data for gene length and then sequencing depth, and makes it easier to compare across samples compared to counts.
 
 ### Workflow inputs
 - Directory containing cDNA/direct RNA reads. Or a directory containing subdirectories each with reads from different samples

docs/quickstart.md

@@ -122,20 +122,37 @@ Differential Expression requires at least 2 replicates of each sample to compare
 
 **Example workflow for differential expression transcript assembly**
 
-Download differential expression data set 
+#### Condition sheet
+The condition sheet should be a .tsv with two columns.
+- The sample column will need to match the 6 directories in the input fastq directory.
+- The condition column will need to contain one of two keys to indicate the two samples being compared.
 
-`wget -O differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_differential_expression.tar.gz && tar -xzvf  differential_expression.tar.gz`
+In the default `condition_sheet.tsv` available in the test_data directory we have used the following.
 
-Run the cmd 
+eg. condition_sheet.tsv
+```
+sample,condition
+barcode01,untreated
+barcode02,untreated
+barcode03,untreated
+barcode04,treated
+barcode05,treated
+barcode06,treated
+```
+
+You will also need to provide a reference genome and a reference annotation file.
+Here is an example cmd to run the workflow using the test_data provided. 
 
 ```
 OUTPUT=~/output;
-nexflow run epi2me-labs/wf-transcriptomes --fastq  differential_expression_dataset/fastq --de_analysis \
+nexflow run epi2me-labs/wf-transcriptomes --fastq test_data/differential_expression_fastq \
---ref_genome differential_expression_dataset/hg38_chr20.fa \
+  --de_analysis \
---ref_annotation differential_expression_dataset/gencode.v22.annotation.chr20.gtf \
+  --ref_genome test_data/hg38_chr20.fa \
---direct_rna
+  --ref_annotation test_data/gencode.v22.annotation.chr20.gtf \
+  --direct_rna
 ```
 
+
 ## Workflow outputs
 * an HTML report document detailing the primary findings of the workflow.
 * for each sample:
@@ -145,11 +162,22 @@ nexflow run epi2me-labs/wf-transcriptomes --fastq  differential_expression_datas
   * merged_transcritptome.fas - annotated, assembled transcriptome
   * [jaffal](https://github.com/Oshlack/JAFFA) ooutput directories
   
-
 ### Fusion detection outputs
 in `${out_dir}/jaffal_output_${sample_id}` you will find:
 * jaffa_results.csv - the csv results summary file 
 * jaffa_results.fasta - fusion transcritpt sequences
 
 ### Differential Expression outputs
-* dtu_plots.pdf - a pdf with differntial transcript usage plots
+* `de_analysis/results_dge.tsv` and `de_analysis/results_dge.pdf`- results of `edgeR` differential gene expression analysis.
+* `de_analysis/results_dtu_gene.tsv`, `de_analysis/results_dtu_transcript.tsv` and `de_analysis/results_dtu.pdf` - results of differential transcript usage by `DEXSeq`.
+* `de_analysis/results_dtu_stageR.tsv` - results of the `stageR` analysis of the `DEXSeq` output.
+* `de_analysis/dtu_plots.pdf` - DTU results plot based on the `stageR` results and filtered counts.
+
+### References
+
+* Benjamini, Yoav, and Yosef Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society. Series B (Methodological) 57 (1): 289–300. http://www.jstor.org/stable/2346101.
+* McCarthy, Davis J., Chen, Yunshun, Smyth, and Gordon K. 2012. “Differential Expression Analysis of Multifactor Rna-Seq Experiments with Respect to Biological Variation.” Nucleic Acids Research 40 (10): 4288–97.
+* Nowicka, Malgorzata, and Mark D. Robinson. 2016. “DRIMSeq: A Dirichlet-Multinomial Framework for Multivariate Count Outcomes in Genomics [Version 2; Referees: 2 Approved].” F1000Research 5 (1356). https://doi.org/10.12688/f1000research.8900.2.
+* Patro, Robert, Geet Duggal, Michael I Love, Rafael A Irizarry, and Carl Kingsford. 2017. “Salmon Provides Fast and Bias-Aware Quantification of Transcript Expression.” Nature Methods 14 (March). https://doi.org/10.1038/nmeth.4197.
+* Robinson, Mark D, Davis J McCarthy, and Gordon K Smyth. 2010. “EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1): 139–40.
+* Love, Michael I., et al. Swimming Downstream: Statistical Analysis of Differential Transcript Usage Following Salmon Quantification. 7:952, F1000Research, 14 Sept. 2018. f1000research.com, https://f1000research.com/articles/7-952

environment.yaml

@@ -6,6 +6,7 @@ channels:
     - defaults
 dependencies:
     - python==3.8.*
+    - bokeh==2.4.3
     - aplanat>=0.6.4
     - epi2melabs
     - minimap2==2.24

lib/fastqingress.nf

@@ -1,7 +1,15 @@
 import ArgumentParser
 
-process handleSingleFile {
+// Downstream tooling assumes FASTQ files are nicely organised into directories.
+// In the case where a single FASTQ file has been input and the parent directory
+// contains other valid FASTQ, we will create a directory in the work area to
+// hold it instead. We stage the file in with `copy` (rather than `link`)
+// to ensure that when the new dir is mounted to containers downstream it does
+// not contain a symlink that cannot be read.
+// See CW-1154
+process isolateSingleFile {
     label params.process_label
+    stageInMode 'copy'
     cpus 1
     input:
         file reads
@@ -9,7 +17,6 @@ process handleSingleFile {
         path "$reads.simpleName"
     script:
         def name = reads.simpleName
-        def reads_dir = 'reads_dir'
     """
     mkdir $name
     mv $reads $name
@@ -51,7 +58,9 @@ def compareSampleSheetFastq(int sample_sheet_count, int valid_dir_count)
 
 /**
  * Take an input file and sample name to return a channel with
- * a single named sample.
+ * a single named sample. If the input file is in a directory with other valid
+ * input files (or other directories containing valid files), a copy of it will
+ * be made to the working directory using the isolateSingleFile process.
  *
  *
  * @param input_file Single fastq file
@@ -61,9 +70,16 @@ def compareSampleSheetFastq(int sample_sheet_count, int valid_dir_count)
 def handle_single_file(input_file, sample_name)
 {
     singleFile = Channel.fromPath(input_file)
-    sample = handleSingleFile(singleFile)
+    ArrayList valid_files_in_dir = find_fastq(input_file.parent, true)
+    if (valid_files_in_dir.size() == 1) {
+        // Avoid a stageInMode copy if the parent directory contains only one valid FASTQ anyway
+        return singleFile.map { it -> tuple(it.parent, create_metamap([sample_id:sample_name ?: it.simpleName])) }
+    }
+    else {
+        // Isolate the file via copy with isolateSingleFile
+        sample = isolateSingleFile(singleFile)
         return sample.map { it -> tuple(it, create_metamap([sample_id:sample_name ?: it.simpleName])) }
-
+    }
 }
 
 
@@ -72,70 +88,45 @@ def handle_single_file(input_file, sample_name)
  * method.
  *
  * @param pattern file object corresponding to top level input folder.
- * @param maxdepth maximum depth to traverse
+ * @param search_subdirs boolean flag to search subdirectories of pattern
  * @return list of files.
  */
 
-def find_fastq(pattern, maxdepth)
+def find_fastq(pattern, search_subdirs)
 {
-    files = []
+    ArrayList files = []
-    extensions = ["fastq", "fastq.gz", "fq", "fq.gz"]
+    ArrayList extensions = ["fastq", "fastq.gz", "fq", "fq.gz"]
     for (ext in extensions) {
-        files += file(pattern.resolve("*.${ext}"), type: 'file', maxdepth: maxdepth)
+        if (search_subdirs) {
+            files += file(pattern.resolve("**.${ext}"), type: 'file')
+        }
+        else {
+            files += file(pattern.resolve("*.${ext}"), type: 'file')
+        }
     }
     return files
 }
 
 
-/**
- * Rework EPI2ME flattened directory structure into standard form
- * files are matched on barcode\d+ and moved into corresponding
- * subdirectories ready for processing.
- *
- * @param input_folder Top-level input directory.
- * @param staging Top-level output_directory.
- * @return A File object representating the staging directory created
- *     under output
- */
-def sanitize_fastq(input_folder, staging)
-{
-    // TODO: this fails if input_folder is an S3 path
-    log.info "Running sanitization."
-    log.info  " - Moving files: ${input_folder} -> ${staging}"
-    staging.mkdirs()
-    files = find_fastq(input_folder.resolve("**"), 1)
-    for (fastq in files) {
-        fname = fastq.getFileName()
-        // find barcode
-        pattern = ~/barcode\d+/
-        matcher = fname =~ pattern
-        if (!matcher.find()) {
-            // not barcoded - leave alone
-            fastq.renameTo(staging.resolve(fname))
-        } else {
-            bc_dir = file(staging.resolve(matcher[0]))
-            bc_dir.mkdirs()
-            fastq.renameTo(staging.resolve("${matcher[0]}/${fname}"))
-        }
-    }
-    log.info " - Finished sanitization."
-    return staging
-}
-
-
 /**
  * Take an input directory return the barcode and non barcode
  * sub directories contained within.
  *
  *
  * @param input_directory Top level input folder to locate sub directories
+ * @param unclassified Keep unclassified directory
+ *
  * @return A list containing sublists of barcode and non_barcode sub directories
  */
-def get_subdirectories(input_directory)
+def get_subdirectories(input_directory, unclassified)
 {
     barcode_dirs = file(input_directory.resolve("barcode*"), type: 'dir', maxdepth: 1)
     all_dirs = file(input_directory.resolve("*"), type: 'dir', maxdepth: 1)
-    non_barcoded = ( all_dirs + barcode_dirs ) - all_dirs.intersect(barcode_dirs)
+    if (!unclassified) {
+        all_dirs.removeIf(it -> it.SimpleName.toLowerCase() == "unclassified")
+    }
+    non_barcoded = (all_dirs + barcode_dirs) - all_dirs.intersect(barcode_dirs)
+    
     return [barcode_dirs, non_barcoded]
 }
 
@@ -183,7 +174,7 @@ def get_valid_directories(input_dirs)
     invalid_files_dirs = []
     for (d in input_dirs) {
         valid = true
-        fastq = find_fastq(d, 1)
+        fastq = find_fastq(d, false)
         all_files = file(d.resolve("*"), type: 'file', maxdepth: 1)
         non_fastq = ( all_files + fastq ) - all_files.intersect(fastq)
 
@@ -345,10 +336,9 @@ def create_metamap(Map arguments) {
  * @param input Top level input file or folder to locate fastq data.
  * @param sample string to name single sample data.
  * @param sample_sheet Path to sample sheet CSV file.
- * @param sanitize regularize inputs from EPI2ME platform.
- * @param output output location, required if sanitize==true
  * @param min_barcode Minimum barcode to accept.
  * @param max_barcode Maximum (inclusive) barcode to accept.
+ * @param unclassified Keep unclassified reads.
  *
  * @return Channel of tuples (path, map(sample_id, type, barcode))
  */
@@ -357,15 +347,12 @@ def fastq_ingress(Map arguments)
     def parser = new ArgumentParser(
         args:["input"],
         kwargs:[
-            "sample":null, "sample_sheet":null, "sanitize":false, "output":null,
+            "sample":null, "sample_sheet":null,
-            "min_barcode":0, "max_barcode":Integer.MAX_VALUE],
+            "min_barcode":0, "max_barcode":Integer.MAX_VALUE,
+            "unclassified":false],
         name:"fastq_ingress")
     Map margs = parser.parse_args(arguments)
 
-    if (margs.sanitize && margs.output == null) {
-        throw new Exception("Argument 'output' required if 'sanitize' is true.")
-    }
-
 
     log.info "Checking fastq input."
     input = file(margs.input)
@@ -382,14 +369,8 @@ def fastq_ingress(Map arguments)
 
     // Handle directory input
     if (input.isDirectory()) {
-        // EPI2ME harness
-        if (margs.sanitize) {
-            staging = file(margs.output).resolve("staging")
-            input = sanitize_fastq(input, staging)
-        }
-
         // Get barcoded and non barcoded subdirectories
-        (barcoded, non_barcoded) = get_subdirectories(input)
+        (barcoded, non_barcoded) = get_subdirectories(input, margs.unclassified)
 
         // Case 03: If no subdirectories, handle the single dir
         if (!barcoded && !non_barcoded) {

main.nf

@@ -475,6 +475,7 @@ workflow pipeline {
             de_report = de.all_de
             count_transcripts_file = de.count_transcripts
             dtu_plots = de.dtu_plots
+            de_outputs = de.de_outputs
         } else{
             de_report = file("$projectDir/data/OPTIONAL_FILE")
             count_transcripts_file = file("$projectDir/data/OPTIONAL_FILE")
@@ -538,7 +539,7 @@ workflow pipeline {
         }
 
         if (params.de_analysis){
-            results = results.concat(de.dtu_plots)
+            results = results.concat(de.dtu_plots, de_outputs)
         }
 
     emit:
@@ -617,9 +618,7 @@ workflow {
         reads = fastq_ingress([
             "input":params.fastq,
             "sample":params.sample,
-            "sample_sheet":params.sample_sheet,
+            "sample_sheet":params.sample_sheet])
-            "sanitize": params.sanitize_fastq,
-            "output":params.out_dir])
 
         pipeline(reads, ref_genome, ref_annotation,
             jaffal_refBase, params.jaffal_genome, params.jaffal_annotation,

nextflow.config

@@ -22,7 +22,6 @@ params {
     out_dir = "output"
     sample = null
     sample_sheet = null
-    sanitize_fastq = false
     wfversion = "v0.1.5"
     aws_image_prefix = null
     aws_queue = null
@@ -246,3 +245,7 @@ trace {
   enabled = true
   file = "${params.out_dir}/execution/trace.txt"
 }
+
+env {
+    PYTHONNOUSERSITE = 1
+}

subworkflows/differential_expression.nf

@@ -51,6 +51,7 @@ process deAnalysis {
         path "merged/all_gene_counts.tsv", emit: gene_counts
         path "de_analysis/results_dge.tsv", emit: dge
         path  "de_analysis/results_dexseq.tsv", emit: dexseq
+        path "de_analysis", emit: de_analysis
     
     """
     cp $annotation annotation.gtf
@@ -74,18 +75,17 @@ process plotResults {
         path flt_count
         path res_dtu
         path condition_sheet 
+        path de_analysis
     output:
         path "de_analysis/dtu_plots.pdf", emit: dtu_plots
         path "condition_sheet.tsv", emit: condition_sheet_tsv
+        path "de_analysis", emit: stageR
     """
     mkdir merged
-    mkdir de_analysis
-    mv $res_dtu de_analysis/results_dtu_stageR.tsv
     mv $condition_sheet de_analysis/coldata.tsv
     mv $flt_count merged/all_counts_filtered.tsv
     plot_dtu_results.R
-    cp de_analysis/coldata.tsv condition_sheet.tsv
+    mv de_analysis/coldata.tsv condition_sheet.tsv
-
     """
 }
 
@@ -144,7 +144,7 @@ workflow differential_expression {
         merged = mergeCounts(count_transcripts.out.counts.collect())
         merged_TPM = mergeTPM(count_transcripts.out.counts.collect())
         analysis = deAnalysis(condition_sheet, merged, ref_annotation)
-        plotResults(analysis.flt_counts, analysis.stageR, condition_sheet)
+        plotResults(analysis.flt_counts, analysis.stageR, condition_sheet, analysis.de_analysis)
         de_report = analysis.flt_counts.combine(analysis.gene_counts).combine(analysis.dge).combine(analysis.dexseq).combine(
                     analysis.stageR).combine(plotResults.out.condition_sheet_tsv).combine(merged).combine(
                     ref_annotation).combine(merged_TPM)
@@ -153,4 +153,5 @@ emit:
        all_de = de_report
        count_transcripts = count_transcripts_file
        dtu_plots = plotResults.out.dtu_plots
+       de_outputs = plotResults.out.stageR
 }

test_data/condition_sheet.tsv

@@ -1,7 +1,6 @@
-sample,condition,type
+sample,condition
-barcode01,untreated,single-read
+barcode01,untreated
-barcode02,untreated,single-read
+barcode02,untreated
-barcode03,untreated,single-read
+barcode04,treated
-barcode04,treated,single-read
+barcode05,treated
-barcode05,treated,single-read
+barcode06,treated
-barcode06,treated,single-read