Only use exons in DE analysis [CW-5786]

CHANGELOG.md

@@ -11,8 +11,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Updated project description.
 - A common user issue is providing a ref_annotation and ref_genome parameter that have mismatched reference IDs, which causes the DE_analysis to fail. The workflow will now do an upfront check and give an error message if no overlap is found or a warning if some IDs are present in one file but not in the other.
 - Reconciled workflow with wf-template v5.5.0.
+- Sort the columns and rows of the gene and transcript count files.
 ### Fixed
 - `all_gene_counts.tsv` contained the DE counts results.
+- The merged transcriptome generated for differential expression analysis now only contains the exons and not the full genomic sequence.
+- Output the gene name annotated differential expression analysis count files only.
+- Only use full length reads in the differential expression analysis.
 
 ## [v1.6.1]
 ### Fixed

bin/workflow_glue/merge_count_tsvs.py

@@ -58,5 +58,7 @@ def main(args):
         left, right, on=args.f, how=args.j), dfs)
     if args.z:
         df_merged = df_merged.fillna(0)
+    df_merged = df_merged.sort_index(axis=1)
+    df_merged = df_merged.sort_index(axis=0)
 
     df_merged.to_csv(args.o, sep="\t", index=False)

main.nf

@@ -422,14 +422,8 @@ process merge_transcriptomes {
         path "final_non_redundant_transcriptome.fasta", emit: fasta
         path "stringtie.gtf", emit: gtf
     """
-    stringtie --rf --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 > "final_non_redundant_transcriptome.fasta"
-    rm temp_transcriptome.fasta
-    rm temp_del_repeats.fasta
-    rm temp_rm_empty_seq.fasta
+    stringtie --merge -G "${ref_annotation}" -p ${task.cpus} -o stringtie.gtf query_annotations/*
+    gffread -g "${ref_genome}" -w "final_non_redundant_transcriptome.fasta" "stringtie.gtf"
     """
 }
 
@@ -760,7 +754,7 @@ workflow pipeline {
                 transcriptome = preprocess_ref_transcriptome(transcriptome)
                 gtf = ref_annotation
             }
-            de = differential_expression(transcriptome, input_reads, sample_sheet, gtf)
+            de = differential_expression(transcriptome, full_len_reads.map{ sample_id, reads -> [[alias:sample_id], reads]}, sample_sheet, gtf)
             de_report = de.all_de
             de_outputs = de.de_outputs
             count_transcripts_file = de.count_transcripts

subworkflows/differential_expression.nf

@@ -177,8 +177,8 @@ workflow differential_expression {
             analysis.gene_counts, analysis.dge, analysis.dexseq,
             analysis.stageR, sample_sheet, merged, ref_annotation, merged_TPM, analysis.unflt_counts).collect()
         // Concat files required to be output to user without any changes
-        de_outputs_concat = analysis.cpm.concat(plotResults.out.dtu_plots, analysis.dge_pdf, analysis.dge_tsv, analysis.dexseq,
-        analysis.dtu_gene, analysis.dtu_transcript, analysis.dtu_stageR, analysis.dtu_pdf, analysis.flt_counts, analysis.gene_counts, merged_TPM).collect()
+        de_outputs_concat = analysis.cpm.concat(plotResults.out.dtu_plots, analysis.dge_pdf, analysis.dge_tsv,
+        analysis.dtu_gene, analysis.dtu_transcript, analysis.dtu_stageR, analysis.dtu_pdf, merged_TPM).collect()
         count_transcripts_file = count_transcripts.out.seqkit_stats.collect()
 emit:
        all_de = de_report