Merge branch 'CW-2498' into 'dev'

tidy up DE outputs

Closes CW-2498

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

CHANGELOG.md

@@ -4,11 +4,19 @@ 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.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [updated]
-### Changed 
+## [unreleased]
+### Changed
+- Improve differential expression outputs.
+- Include transcript and gene count tables in DE_final folder.
+- If differential expression subworkflow is used a non redundant transcriptome will be output which includes novel transcripts.
+- Added wording to the report about how to identify novel transcripts in the DE tables.
 - Nextflow minimum required version to 23.04.2
 - `--minimap_index_opts` parameter has been changed to `minimap2_index_opts` for consistency.
 
+### Added
+- An additional gene name column to the differential gene expression results. This is especially handy for transcriptomes where the gene ID is not the same as gene name (e.g. Ensembl).
+- Wording to the report about how to identify novel transcripts in the DE tables.
+
 ## [v0.2.1]
 ### Changed
 - Any sample aliases that contain spaces will be replaced with underscores.

README.md

@@ -276,6 +276,11 @@ in `${out_dir}/jaffal_output_${sample_id}` you will find:
 * `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.
+* `de_analysis/all_gene_counts.tsv` - Gene counts generated by the `salmon` tool before filtering.
+* `de_analysis/de_transcript_counts.tsv` - Transcript counts generated by the `salmon` tool before filtering. 
+* `de_analysis/de_tpm_transcript_counts.tsv` - To facilitate comparisons across samples this file shows transcript per million (TPM) of the raw counts. 
+* `de_analysis/all_counts_filtered.tsv` - Transcript counts filtered with input criteria. Used for DE analysis.
+* `final_non_redundant_transcriptome.fasta` - Transcripts that were used for differential expression analysis including novel transcripts with the identifiers used for DE analysis. 
 
 ### References
 

bin/workflow_glue/de_plots.py

@@ -237,6 +237,13 @@ def dtu_section(dtu_file, section, gt_dic, ge_dic):
     section.table(dtu_results.loc[dtu_pvals.index])
 
 
+def dge_names(dge_file, geid_gname):
+    """Add gene name column to DGE tsv."""
+    dge_results = pd.read_csv(dge_file, sep='\t')
+    dge_results["gene_name"] = dge_results.index.map(lambda x: geid_gname.get(x))
+    dge_results.to_csv('results_dge.tsv', index=True, index_label="gene_id")
+
+
 def dge_section(dge_file, section, ids_dic):
     """Create DGE table and plot."""
     section.markdown('### Differential gene expression')
@@ -312,6 +319,7 @@ def get_translations(gtf):
     fn = open(gtf).readlines()
     gene_txid = {}
     gene_geid = {}
+    geid_gname = {}
 
     def get_feature(row, feature):
         return row.split(feature)[1].split(
@@ -347,7 +355,8 @@ def get_translations(gtf):
             gene_id = gene_name
         gene_txid[transcript_id] = gene_name
         gene_geid[gene_id] = gene_reference
-    return gene_txid, gene_geid
+        geid_gname[gene_reference] = gene_name
+    return gene_txid, gene_geid, geid_gname
 
 
 def de_section(
@@ -364,6 +373,10 @@ the GTF-format annotation.
 These counts were used to perform a statistical analysis to identify
 the genes and isoforms that show differences in abundance between
 the experimental conditions.
+Any novel genes or transcripts that do not have relevant gene or transcript IDs
+are prefixed with MSTRG for use in differential expression analysis.
+Find the full sequences of any transcripts in the
+`final_non_redundant_transcriptome.fasta` file.
     """)
     section.markdown("### Alignment summary stats")
     alignment_stats = pool_csvs("seqkit")
@@ -371,8 +384,9 @@ the experimental conditions.
     alignment_summary_df = alignment_summary_df.fillna(0).applymap(np.int64)
     section.table(alignment_summary_df, key='alignment-stats', index=True)
     salmon_table(tpm, section)
-    gene_txid, gene_name = get_translations(stringtie)
+    gene_txid, gene_name, geid_gname = get_translations(stringtie)
     dge_section(dge, section, gene_name)
+    dge_names(dge, geid_gname)
     dexseq_section(dexseq, section, gene_name)
     dtu_section(dtu, section, gene_txid, gene_name)
     # missing dtu plots at the moment as too many

bin/workflow_glue/report.py

@@ -888,7 +888,8 @@ def de_section(report):
     dge = os.path.join("de_report", "results_dge.tsv")
     dtu = os.path.join("de_report", "results_dtu_stageR.tsv")
     stringtie = os.path.join("de_report", "stringtie_merged.gtf")
-    tpm = os.path.join("de_report", "tpm_counts.tsv")
+    tpm = os.path.join("de_report", "de_tpm_transcript_counts.tsv")
+    # This will also add a gene name column to the "results_dge.tsv"
     de_plots.de_section(
         stringtie=stringtie,
         dexseq=dexseq,

docs/quickstart.md

@@ -186,6 +186,11 @@ in `${out_dir}/jaffal_output_${sample_id}` you will find:
 * `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.
+* `de_analysis/all_gene_counts.tsv` - Gene counts generated by the `salmon` tool before filtering.
+* `de_analysis/de_transcript_counts.tsv` - Transcript counts generated by the `salmon` tool before filtering. 
+* `de_analysis/de_tpm_transcript_counts.tsv` - To facilitate comparisons across samples this file shows transcript per million (TPM) of the raw counts. 
+* `de_analysis/all_counts_filtered.tsv` - Transcript counts filtered with input criteria. Used for DE analysis.
+* `final_non_redundant_transcriptome.fasta` - Transcripts that were used for differential expression analysis including novel transcripts with the identifiers used for DE analysis. 
 
 ### References
 

main.nf

@@ -309,14 +309,14 @@ process merge_transcriptomes {
         path ref_annotation
         path ref_genome
     output:
-        path "non_redundant.fasta", emit: fasta
+        path "final_non_redundant_transcriptome.fasta", emit: fasta
         path "stringtie.gtf", emit: gtf
     """
     stringtie --merge -G $ref_annotation -p ${task.cpus} -o stringtie.gtf query_annotations/*
     seqkit subseq --feature "transcript" --gtf-tag "transcript_id" --gtf stringtie.gtf $ref_genome > temp_transcriptome.fasta
     seqkit rmdup -s < temp_transcriptome.fasta > temp_del_repeats.fasta
     cat temp_del_repeats.fasta | sed 's/>.* />/'  | sed -e 's/_[0-9]* \\[/ \\[/' > temp_rm_empty_seq.fasta
-    awk 'BEGIN {RS = ">" ; FS = "\\n" ; ORS = ""} \$2 {print ">"\$0}' temp_rm_empty_seq.fasta > non_redundant.fasta
+    awk 'BEGIN {RS = ">" ; FS = "\\n" ; ORS = ""} \$2 {print ">"\$0}' temp_rm_empty_seq.fasta > "final_non_redundant_transcriptome.fasta"
     rm temp_transcriptome.fasta
     rm temp_del_repeats.fasta
     rm temp_rm_empty_seq.fasta
@@ -342,6 +342,9 @@ process makeReport {
         path "seqkit/*"
     output:
         path("wf-transcriptomes-*.html"), emit: report
+        // for DE analysis, a `gene_name` column will be added to
+        // `de_report/results_dge.tsv`
+        path "results_dge.tsv", emit: de_analysis, optional: true
     script:
         // Convert the sample_id arrayList.
         sids = new BlankSeparatedList(sample_ids)
@@ -387,7 +390,6 @@ process makeReport {
     \$OPT_JAFFAL_CSV \
     $OPT_DENOVO \
     \$dereport 
-
     """
 }
 
@@ -495,7 +497,7 @@ workflow pipeline {
             return l
         }
 
-        
+     
 
         software_versions = getVersions()
         workflow_params = getParams()
@@ -578,7 +580,7 @@ workflow pipeline {
                 gtf = merge_transcriptomes.out.gtf
             }
             else {
-                transcriptome = ref_transcriptome
+                transcriptome =  Channel.fromPath(ref_transcriptome)
                 gtf = ref_annotation
             }
             de = differential_expression(transcriptome, input_reads, sample_sheet, gtf)
@@ -586,6 +588,7 @@ workflow pipeline {
             count_transcripts_file = de.count_transcripts
             dtu_plots = de.dtu_plots
             de_outputs = de.de_outputs
+            counts = de.counts
         } else{
             de_report = file("$projectDir/data/OPTIONAL_FILE")
             count_transcripts_file = file("$projectDir/data/OPTIONAL_FILE")
@@ -605,6 +608,8 @@ workflow pipeline {
             count_transcripts_file)
 
         report = makeReport.out.report
+        
+     
        
         results = results.concat(makeReport.out.report)
        
@@ -645,15 +650,19 @@ workflow pipeline {
                 .concat(gene_fusions.out.results
                 .map {it -> it[1]})
         }
-
+        
+       results = results.map{ [it, null] }.concat(fastq_ingress_results.map { [it, "fastq_ingress_results"] })
+        
         if (params.de_analysis){
-            results = results.concat(de.dtu_plots, de_outputs)
+           de_update = makeReport.out.de_analysis
+           de_results = report.concat(transcriptome, de_outputs.flatten(), counts.flatten(), de_update)
+           results = results.concat(de_results.map{ [it, "de_analysis"] })
         }
 
-        results  = fastq_ingress_results.map { [it, "fastq_ingress_results"] }.concat(results.map{ [it, null]})
+        results.concat(workflow_params.map{ [it, null]})
+       
     emit:
         results
-        telemetry = workflow_params
 }
 
 // entrypoint workflow

subworkflows/differential_expression.nf

@@ -32,9 +32,9 @@ process mergeCounts {
     input:
         path counts
     output:
-        path "all_counts.tsv"
+        path "de_transcript_counts.tsv"
     """
-    workflow-glue merge_count_tsvs -z -o all_counts.tsv -tsvs ${counts}
+    workflow-glue merge_count_tsvs -z -o de_transcript_counts.tsv -tsvs ${counts}
     """
 }
 
@@ -43,9 +43,9 @@ process mergeTPM {
     input:
         path counts
     output:
-        path "tpm_counts.tsv"
+        path "de_tpm_transcript_counts.tsv"
     """
-    workflow-glue merge_count_tsvs -o tpm_counts.tsv -z -tpm True -tsvs $counts
+    workflow-glue merge_count_tsvs -o de_tpm_transcript_counts.tsv -z -tpm True -tsvs $counts
     """
 }
 
@@ -104,7 +104,7 @@ process plotResults {
     output:
         path "de_analysis/dtu_plots.pdf", emit: dtu_plots
         path "sample_sheet.tsv", emit: sample_sheet_csv
-        path "de_analysis", emit: stageR
+        path "de_analysis/*", emit: stageR
     """
     mkdir merged
     mv $sample_sheet de_analysis/coldata.tsv
@@ -173,9 +173,11 @@ workflow differential_expression {
                     analysis.stageR).combine(plotResults.out.sample_sheet_csv).combine(merged).combine(
                     ref_annotation).combine(merged_TPM)
         count_transcripts_file = count_transcripts.out.seqkit_stats.collect()
+        all_counts = merged_TPM.concat(merged, analysis.flt_counts, analysis.gene_counts)
 emit:
        all_de = de_report
        count_transcripts = count_transcripts_file
        dtu_plots = plotResults.out.dtu_plots
        de_outputs = plotResults.out.stageR
+       counts = all_counts
 }

test_data/sample_sheet.csv

@@ -4,4 +4,4 @@ barcode02,sample02,sample02,untreated
 barcode03,sample03,sample03,untreated
 barcode04,sample04,sample04,treated
 barcode05,sample05,sample05,treated
-barcode06,sample06,sample06,treated
+barcode06,sample06,sample06,treated