From e3b196787272a369478faa259b3c9e337d9e9bbb Mon Sep 17 00:00:00 2001
From: Sarah Griffiths <sarah.griffiths@nanoporetech.com>
Date: Mon, 19 Feb 2024 10:29:04 +0000
Subject: [PATCH] Improve ref annotation transcript id version handling

---
 CHANGELOG.md                            |  2 ++
 bin/de_analysis.R                       | 25 +++++++++++++++++++------
 subworkflows/differential_expression.nf | 18 ++++--------------
 3 files changed, 25 insertions(+), 20 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index e8198f4..0eca851 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,8 @@ 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).
 
 ## [unreleased]
+### Changed
+- Improved handling of annotation files that do not contain version numbers in transcript_id (such as gtf's from Ensembl).
 ### Fixed
 - Differential expression failing with 10 or more samples.
 - Regression causing the DE analysis numeric parameters to not be evaluated correctly.
diff --git a/bin/de_analysis.R b/bin/de_analysis.R
index 1927c7a..cdc8f45 100755
--- a/bin/de_analysis.R
+++ b/bin/de_analysis.R
@@ -10,7 +10,6 @@ min_samps_feature_expr <- as.numeric(args[3])
 min_gene_expr <- as.numeric(args[4]) 
 min_feature_expr <- as.numeric(args[5])
 annotation_type <- args[6]
-strip_version <- args[7]
 
 cat("Loading counts, conditions and parameters.\n")
 cts <- as.matrix(read.csv("all_counts.tsv", sep="\t", row.names="Reference", stringsAsFactors=FALSE))
@@ -26,18 +25,32 @@ if(!"control" %in% coldata$condition)
        condition - unable to set reference.")
 coldata$condition <- relevel(coldata$condition, ref = "control")
 
-cat("Loading annotation database.\n")
 
+# Transcript_id versions (eg. ENTXXX.1, eg. ENTXXX.2) represent how many times that transcript reference has been changed 
+# during its time in the database.
+# Not all annotation files include it as part of the transcript_id - notably Ensembl
+# The following handles this.
+cat("Checking annotation file for presence of transcript_id versions.\n")
+# Get the first transcript_id from the annotation file by parsing
+lines <- readLines(file(ref_annotation), n=100000)
+# Find transcript_ids in first 1000 lines and check if they contain dot (format eg. ENTXXX.1)
+check_version <- sum(grepl("transcript_id[^;]+\\.", lines))
+if (check_version != 0){
+        # we do not need to strip the count file rows if ref_annotation includes versions
+        cat("Annotation file transcript_ids include versions.\n")
+    } else {
+       # otherwise remove the versions
+        rownames(cts) <- lapply(rownames(cts),  sub, pattern = "\\.\\d+$", replacement = "")
+        cat("Annotation file transcript_ids do not include versions so also strip versions from the counts df.\n")
+    }
+
+cat("Loading annotation database.\n")
 txdb <- makeTxDbFromGFF(ref_annotation,  format = annotation_type)
 txdf <- select(txdb, keys(txdb,"GENEID"), "TXNAME", "GENEID")
 tab <- table(txdf$GENEID)
 txdf$ntx<- tab[match(txdf$GENEID, names(tab))]
 
 
-if (strip_version == "true"){
-  rownames(cts) <- lapply(rownames(cts),  sub, pattern = "\\.\\d+$", replacement = "")
-}
-
 cts <- cts[rownames(cts) %in% txdf$TXNAME, ] # FIXME: filter for transcripts which are in the annotation. Why they are not all there? 
 
 # Reorder transcript/gene database to match input counts:
diff --git a/subworkflows/differential_expression.nf b/subworkflows/differential_expression.nf
index 471bc5f..a6e835e 100644
--- a/subworkflows/differential_expression.nf
+++ b/subworkflows/differential_expression.nf
@@ -61,7 +61,9 @@ process mergeTPM {
 process deAnalysis {
     label "isoforms"
     errorStrategy "retry"
-    maxRetries 3
+    // Retry if it fails to make makeTxDbFromGFF 
+    // Because a file with .gff extension may be gff2(gtf) or gff3
+    maxRetries 1
     cpus 4
     memory "16 GB"
     input:
@@ -80,26 +82,14 @@ process deAnalysis {
     script:
     // Just try both annotation file type because a .gff extension may be gff2(gtf) or gff3
     String annotation_type = "gtf"
-    String strip_version  = "false"
     if (task.attempt == 2){
         annotation_type = "gff3"
-        strip_version  = "false"
         log.info("Retry deAnalysis with gff format setting.")
     }
-    else if (task.attempt == 3){
-        annotation_type = "gff3"
-        strip_version  = "true"
-        log.info("Retry deAnalysis with gff format setting and version removal.")
-    }
-    else if (task.attempt == 4){
-        strip_version  = "true"
-        log.info("Retry deAnalysis with gtf format setting and version removal.")
-    }
-
     """
     mkdir merged
     mkdir de_analysis
-    de_analysis.R annotation.gtf $params.min_samps_gene_expr $params.min_samps_feature_expr $params.min_gene_expr $params.min_feature_expr $annotation_type $strip_version
+    de_analysis.R annotation.gtf $params.min_samps_gene_expr $params.min_samps_feature_expr $params.min_gene_expr $params.min_feature_expr $annotation_type
     """
 }