Merge branch 'CW-6125' into 'dev'

Filter out unstranded annotations and log warning [CW-6125]

Closes CW-6125

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

.gitlab-ci.yml

@@ -54,7 +54,7 @@ docker-run:
                 "only_differential_expression", "differential_expression_gff3",
                 "ncbi_gzip", "ncbi_no_gene_id", "ensembl_with_versions",
                 "differential_expression_mouse", "no_ref_annotation",
-                "unstranded_annotation_error", "igv", "igv_fai_gz",
+                "unstranded_annotation_warn", "igv", "igv_fai_gz",
                 "mismatch_seqid_error","mismatch_seqid_warn"
             ]
     rules:
@@ -68,21 +68,21 @@ docker-run:
               NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
               NF_WORKFLOW_OPTS: "--fastq  ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq --transcriptome_source reference-guided \
                   --ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf --pychopper_backend phmm"
-              NF_IGNORE_PROCESSES: preprocess_reads,validate_ref_annotation,faidx,gz_faidx,check_annotation_strand,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
+              NF_IGNORE_PROCESSES: preprocess_reads,filter_unstranded_annotation,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
         - if: $MATRIX_NAME == "isoforms_bam"
           variables:
             NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
             NF_WORKFLOW_OPTS: "--bam  ${CI_PROJECT_NAME}/data/ERR6053095_chr20.bam --transcriptome_source reference-guided \
                   --ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf --pychopper_backend phmm \
                   -c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
-            NF_IGNORE_PROCESSES: preprocess_reads,validate_ref_annotation,faidx,gz_faidx,check_annotation_strand,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
+            NF_IGNORE_PROCESSES: preprocess_reads,filter_unstranded_annotation,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
         - if: $MATRIX_NAME == "no_ref_annotation"
           variables:
               NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf  ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
               NF_WORKFLOW_OPTS: "--fastq  ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq --transcriptome_source reference-guided \
                   --ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa \
                   -c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
-              NF_IGNORE_PROCESSES: run_gffcompare,validate_ref_annotation,check_annotation_strand,preprocess_reads,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
+              NF_IGNORE_PROCESSES: run_gffcompare,filter_unstranded_annotation,validate_ref_annotation,preprocess_reads,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
         - if: $MATRIX_NAME == "differential_expression"
           variables:
               NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
@@ -179,7 +179,7 @@ docker-run:
               NF_IGNORE_PROCESSES: >
                 preprocess_reads,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,
                 build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam    
-        - if: $MATRIX_NAME == "unstranded_annotation_error"
+        - if: $MATRIX_NAME == "unstranded_annotation_warn"
           variables:
               NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/  &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
               NF_WORKFLOW_OPTS: "--fastq  ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
@@ -189,8 +189,7 @@ docker-run:
                 --direct_rna --minimap2_index_opts '-k 15'  --sample_sheet ${CI_PROJECT_NAME}/data/differential_expression/sample_sheet.csv \
                 -c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
               NF_IGNORE_PROCESSES: preprocess_reads,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
-              ASSERT_NEXTFLOW_FAILURE: "test_fail" # set to any non-zero length str to allow the nextflow CMD to fail
+              AFTER_NEXTFLOW_CMD: "grep 'Unstranded entries found and excluded from differential expression analysis.' '.nextflow.log'"
-              ASSERT_NEXTFLOW_FAILURE_REXP: "In ref_annotation, transcript features must have a strand of either '+' or '-'"
         - if: $MATRIX_NAME == "igv"
           variables:
               NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/  &&  wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"

CHANGELOG.md

@@ -18,6 +18,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
     - Sequence summary read length N50 incorrectly displayed minimum read length, it now correctly shows the N50.
     - Sequence summary component alignment and coverage plots failed to plot under some conditions.
 - Error in `deAnalysis` process - `mode(counts) %in% "numeric" is not TRUE` - caused by hyphens in sample sheet aliases.
+- Error in `deAnalysis` process - `values in 'transcripts$tx_strand' must be "+" or "-"`.
+  - The workflow will now filter out any unstranded annotations from downstream analysis and log a warning.
 
 ## [v1.7.0]
 ### Changed

main.nf

@@ -540,16 +540,37 @@ process publish_results {
 
 
 // Check ref_annotation transcript strand column for "." if in de_analysis mode
-process check_annotation_strand {
+process filter_unstranded_annotation {
     label "isoforms"
     cpus 1
     memory "2 GB"
     input:
         path "ref_annotation.gtf"
     output:
-        tuple stdout, path("ref_annotation.gtf")   
+        tuple stdout, path("ref_annotation_stranded.gtf")   
+
+    script:
     """
-    awk '{if (\$3=="transcript" && \$7 != "+" && \$7 != "-") print \$3, \$7}' "ref_annotation.gtf"
+    awk '
+        BEGIN { OFS = "\\t" }
+        /^#/ { print; next }
+    # Valid strands
+    (\$7 == "+" || \$7 == "-") { print; next }
+
+    # Invalid strand
+    {
+        print \$0 >> "unstranded.gtf"
+    }
+    ' ref_annotation.gtf > ref_annotation_stranded.gtf
+
+    # Check and log warning if unstranded entries exist
+    if [ -s unstranded.gtf ]; then
+        echo "Warning: Unstranded entries found and excluded from differential expression analysis.
+        If running with reference-guided transcriptome source, consider increasing read depth to reduce unstranded annotations.
+        If running with precomputed transcriptome source, ensure your ref_annotation gtf file contains only '+' or '-' strand entries."
+        echo "A sample of unstranded entries:"
+        head -n 20 unstranded.gtf
+    fi
     """
 }
 
@@ -727,12 +748,6 @@ workflow pipeline {
         }
         if (params.de_analysis){
             sample_sheet = file(params.sample_sheet, type:"file")
-            // check ref annotation contains only + or - strand as DE analysis will error on .
-            check_annotation_strand(ref_annotation).map { stdoutput, annotation ->
-            // check if there was an error message
-            if (stdoutput) error "In ref_annotation, transcript features must have a strand of either '+' or '-'."
-                    stdoutput
-                }
             if (!params.ref_transcriptome){
                 validate_ref_annotation(ref_annotation, ref_genome).map { stdoutput ->
                     if (stdoutput) {
@@ -751,7 +766,19 @@ workflow pipeline {
                 transcriptome = preprocess_ref_transcriptome(transcriptome)
                 gtf = ref_annotation
             }
-            de = differential_expression(transcriptome, full_len_reads.map{ sample_id, reads -> [[alias:sample_id], reads]}, sample_sheet, gtf)
+            // Filter out any .(unstranded) records
+            // only + or - strand allowed for DE analysis.
+            stranded_annotation = filter_unstranded_annotation(gtf)
+                .map{ stdoutput, annotation -> annotation
+                    if (stdoutput) {
+                        log.warn(stdoutput)
+                    }
+                annotation
+                }
+            de = differential_expression(
+                transcriptome,
+                full_len_reads.map{ sample_id, fq_reads -> [[alias:sample_id], fq_reads]},
+                sample_sheet, stranded_annotation)
             de_report = de.all_de
             de_outputs = de.de_outputs
             de_alignment_stats = de.de_alignment_stats