Add 5' truncation support and overlap file export

DESCRIPTION

@@ -6,15 +6,21 @@ Description: Various mRNA sequencing library preparation methods generate
     truncated versions of transcriptome annotations, both at the alignment or 
     pseudo-alignment stage, as well as in downstream analysis. This package 
     implements some convenience methods for readily generating such truncated 
-    annotations and their corresponding sequences.
-Version: 1.15.2
-Date: 2025-09-07
-Authors@R: 
+    annotations from either their 5' or 3' transcript ends and their 
+    corresponding sequences.
+Version: 1.17.0
+Date: 2026-03-10
+Authors@R: c(
     person(given = "Mervin",
            family = "Fansler",
            role = c("aut", "cre"),
            email = "mervin.fansler@bric.ku.dk",
-           comment = c(ORCID = "0000-0002-4108-4218"))
+           comment = c(ORCID = "0000-0002-4108-4218")),
+    person(given = "Guillermo",
+           family = "Rocamora Pérez",
+           role = c("ctb"),
+           email = "guillermorocamora@gmail.com",
+           comment = c(ORCID = "0000-0002-4795-3648")))
 License: GPL-3
 URL: https://github.com/mfansler/txcutr
 BugReports: https://github.com/mfansler/txcutr/issues
@@ -36,7 +42,7 @@ Imports:
     methods,
     utils
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.3
 biocViews: 
     Alignment,
     Annotation,
@@ -52,6 +58,7 @@ Suggests:
     testthat (>= 3.0.0),
     TxDb.Scerevisiae.UCSC.sacCer3.sgdGene,
     BSgenome.Scerevisiae.UCSC.sacCer3,
-    GenomeInfoDbData
+    GenomeInfoDbData,
+    withr
 VignetteBuilder: knitr
 Config/testthat/edition: 3

NAMESPACE

@@ -4,10 +4,15 @@ export(exportFASTA)
 export(exportGTF)
 export(exportMergeTable)
 export(generateMergeTable)
+export(truncate3primeTxome)
+export(truncate5primeTxome)
 export(truncateTxome)
 export(txdbToGRangesList)
 exportMethods(generateMergeTable)
+exportMethods(truncate3primeTxome)
+exportMethods(truncate5primeTxome)
 exportMethods(truncateTxome)
+importFrom(AnnotationDbi,metadata)
 importFrom(AnnotationDbi,select)
 importFrom(AnnotationDbi,taxonomyId)
 importFrom(BiocGenerics,paste)
@@ -27,6 +32,7 @@ importFrom(GenomicRanges,GRangesList)
 importFrom(GenomicRanges,end)
 importFrom(GenomicRanges,findOverlaps)
 importFrom(GenomicRanges,intersect)
+importFrom(GenomicRanges,invertStrand)
 importFrom(GenomicRanges,mcols)
 importFrom(GenomicRanges,reduce)
 importFrom(GenomicRanges,resize)
@@ -49,4 +55,5 @@ importFrom(rtracklayer,export)
 importFrom(stats,setNames)
 importFrom(txdbmaker,makeTxDbFromGRanges)
 importFrom(utils,capture.output)
+importFrom(utils,write.csv)
 importFrom(utils,write.table)

R/merge.R

@@ -10,8 +10,9 @@
 #'
 #' @importFrom methods setGeneric
 #' @export
-setGeneric("generateMergeTable", signature=c("txdb", "minDistance"),
-           function(txdb, minDistance=200) standardGeneric("generateMergeTable")
+setGeneric("generateMergeTable",
+  signature = c("txdb", "minDistance"),
+  function(txdb, minDistance = 200) standardGeneric("generateMergeTable")
 )
 
 #' Generate Merge Table
@@ -40,74 +41,96 @@ setGeneric("generateMergeTable", signature=c("txdb", "minDistance"),
 #' txdb_w500
 #'
 #' ## last 100 nts per tx
-#' txdb_w100 <- truncateTxome(txdb, maxTxLength=100)
+#' txdb_w100 <- truncateTxome(txdb, maxTxLength = 100)
 #' txdb_w100
 #'
-#' @importFrom GenomicRanges mcols resize findOverlaps
+#' @importFrom GenomicRanges mcols resize findOverlaps invertStrand
 #' @importFrom GenomicFeatures transcripts
+#' @importFrom AnnotationDbi metadata
 #' @importFrom methods setMethod
 #' @importFrom stats setNames
 #' @export
-setMethod("generateMergeTable", "TxDb", function(txdb, minDistance=200L) {
-    grTxs <- transcripts(txdb,
-                         columns=c("gene_id", "tx_id", "tx_name"))
-
-    ## use only ends
-    grTxs <- resize(grTxs, width=minDistance, fix="end", ignore.strand=FALSE)
-
-    ## generate raw overlaps
-    overlaps <- as.data.frame(findOverlaps(grTxs, ignore.strand=FALSE,
-                                           drop.self=TRUE))
-    overlaps["tx_in"] <- unlist(grTxs$tx_name[overlaps$queryHits])
-    overlaps["tx_out"] <- unlist(grTxs$tx_name[overlaps$subjectHits])
-    overlaps["gene_in"] <- unlist(grTxs$gene_id[overlaps$queryHits])
-    overlaps["gene_out"] <- unlist(grTxs$gene_id[overlaps$subjectHits])
-
-    ## filter unmatched genes
-    overlaps <- overlaps[overlaps$gene_in == overlaps$gene_out,]
-
-    ## keep downstream
-    overlaps["strand"] <- as.character(strand(grTxs))[overlaps$queryHits]
-    overlaps["end_in"] <- ifelse(overlaps$strand == "+",
-                                 end(grTxs[overlaps$queryHits]),
-                                 -start(grTxs[overlaps$queryHits]))
-    overlaps["end_out"] <- ifelse(overlaps$strand == "+",
-                                  end(grTxs[overlaps$subjectHits]),
-                                  -start(grTxs[overlaps$subjectHits]))
-    overlaps <- overlaps[overlaps$end_in <= overlaps$end_out,]
-
-    ## if equal, prioritize early transcript
-    ## e.g., a transcript ENSMUST000000001 will replace ENSMUST000001000
-    idxEqual <- which(overlaps$end_in == overlaps$end_out)
-    if (length(idxEqual) > 0) {
-        isLaterTx <- overlaps[idxEqual, "tx_in"] < overlaps[idxEqual, "tx_out"]
-        overlaps <- overlaps[-idxEqual[isLaterTx], ]
+setMethod("generateMergeTable", "TxDb", function(txdb, minDistance = 200L) {
+  ############################################################################
+  # Extract txEnd from txdb metadata
+  if (!"Truncation End" %in% metadata(txdb)$name) {
+    warning("'Truncation End' parameter not found in the TxDb metadata. Defaults to '3prime'. Was this object created with txcutr::truncateTxome?")
+    txEnd <- "3prime"
+  } else {
+    txEnd <- metadata(txdb)[metadata(txdb)$name == "Truncation End", "value"]
+  }
+  if (!txEnd %in% c("3prime", "5prime")) stop("'Truncation End' parameter is not valid - only '3prime' or '5prime' are accepted.")
+
+  ############################################################################
+  # Merge pipeline
+  grTxs <- transcripts(txdb,
+    columns = c("gene_id", "tx_id", "tx_name")
+  )
+
+  ## Invert strand if 5' truncation
+  if (txEnd == "5prime") grTxs <- invertStrand(grTxs)
+
+  ## use only ends
+  grTxs <- resize(grTxs, width = minDistance, fix = "end", ignore.strand = FALSE)
+
+  ## generate raw overlaps
+  overlaps <- as.data.frame(findOverlaps(grTxs,
+    ignore.strand = FALSE,
+    drop.self = TRUE
+  ))
+  overlaps["tx_in"] <- unlist(grTxs$tx_name[overlaps$queryHits])
+  overlaps["tx_out"] <- unlist(grTxs$tx_name[overlaps$subjectHits])
+  overlaps["gene_in"] <- unlist(grTxs$gene_id[overlaps$queryHits])
+  overlaps["gene_out"] <- unlist(grTxs$gene_id[overlaps$subjectHits])
+
+  ## filter unmatched genes
+  overlaps <- overlaps[overlaps$gene_in == overlaps$gene_out, ]
+
+  ## keep downstream
+  overlaps["strand"] <- as.character(strand(grTxs))[overlaps$queryHits]
+  overlaps["end_in"] <- ifelse(overlaps$strand == "+",
+    end(grTxs[overlaps$queryHits]),
+    -start(grTxs[overlaps$queryHits])
+  )
+  overlaps["end_out"] <- ifelse(overlaps$strand == "+",
+    end(grTxs[overlaps$subjectHits]),
+    -start(grTxs[overlaps$subjectHits])
+  )
+  overlaps <- overlaps[overlaps$end_in <= overlaps$end_out, ]
+
+  ## if equal, prioritize early transcript
+  ## e.g., a transcript ENSMUST000000001 will replace ENSMUST000001000
+  idxEqual <- which(overlaps$end_in == overlaps$end_out)
+  if (length(idxEqual) > 0) {
+    isLaterTx <- overlaps[idxEqual, "tx_in"] < overlaps[idxEqual, "tx_out"]
+    overlaps <- overlaps[-idxEqual[isLaterTx], ]
+  }
+
+  ## pick unique out
+  groupedOverlaps <- split(overlaps, overlaps$queryHits)
+  singleOverlaps <- lapply(
+    groupedOverlaps,
+    function(df) {
+      df[with(df, order(-end_out, tx_out)), ][1, ]
     }
+  )
+  dfMerge <- do.call(rbind, singleOverlaps)[, c("tx_in", "tx_out")]
+  dfMerge <- .propagateMap(dfMerge)
 
-    ## pick unique out
-    groupedOverlaps <- split(overlaps, overlaps$queryHits)
-    singleOverlaps <- lapply(groupedOverlaps,
-                             function (df) {
-                                 df[with(df, order(-end_out, tx_out)),][1,]
-                             })
-    dfMerge <- do.call(rbind, singleOverlaps)[, c("tx_in", "tx_out")]
-    dfMerge <- .propagateMap(dfMerge)
+  ## include self-maps
+  unmergedTxs <- grTxs$tx_name[!(grTxs$tx_name %in% dfMerge$tx_in)]
+  dfMerge <- rbind(dfMerge, data.frame(tx_in = unmergedTxs, tx_out = unmergedTxs))
 
-    ## include self-maps
-    unmergedTxs <- grTxs$tx_name[!(grTxs$tx_name %in% dfMerge$tx_in)]
-    dfMerge <- rbind(dfMerge, data.frame(tx_in=unmergedTxs, tx_out=unmergedTxs))
+  ## append gene
+  txToGeneMap <- setNames(object = unlist(grTxs$gene_id), nm = grTxs$tx_name)
+  dfMerge["gene_out"] <- txToGeneMap[as.character(dfMerge$tx_out)]
 
-    ## append gene
-    txToGeneMap <- setNames(object=unlist(grTxs$gene_id), nm=grTxs$tx_name)
-    dfMerge['gene_out'] <- txToGeneMap[as.character(dfMerge$tx_out)]
+  ## reorder and drop rownames
+  dfMerge <- dfMerge[with(dfMerge, order(tx_in)), ]
+  rownames(dfMerge) <- NULL
 
-    ## reorder and drop rownames
-    dfMerge <- dfMerge[with(dfMerge, order(tx_in)),]
-    rownames(dfMerge) <- NULL
-
-    dfMerge
-}
-)
+  dfMerge
+})
 
 
 #' Propagate Transcript Merge Map
@@ -119,32 +142,32 @@ setMethod("generateMergeTable", "TxDb", function(txdb, minDistance=200L) {
 #'     present in any \code{tx_in}
 #'
 #' @importFrom stats setNames
-.propagateMap <- function (df, MAXITERS=1000) {
-    ## NOTE: This method assumes no loops. To provide execution safety, we
-    ## enforce a maximum number of iterations.
-    iteration <- 0
+.propagateMap <- function(df, MAXITERS = 1000) {
+  ## NOTE: This method assumes no loops. To provide execution safety, we
+  ## enforce a maximum number of iterations.
+  iteration <- 0
 
-    ## while there are any outputs in the inputs,
-    while(any(df$tx_out %in% df$tx_in) && iteration < MAXITERS) {
-        iteration <- iteration + 1
+  ## while there are any outputs in the inputs,
+  while (any(df$tx_out %in% df$tx_in) && iteration < MAXITERS) {
+    iteration <- iteration + 1
 
-        ## generate merge map for simple translation
-        mergeMap <- setNames(object=df$tx_out, nm=df$tx_in)
+    ## generate merge map for simple translation
+    mergeMap <- setNames(object = df$tx_out, nm = df$tx_in)
 
-        ## these will be propagated
-        idxMappable <- df$tx_out %in% df$tx_in
+    ## these will be propagated
+    idxMappable <- df$tx_out %in% df$tx_in
 
-        ## outputs to translate
-        oldOut <- as.character(df$tx_out[idxMappable])
+    ## outputs to translate
+    oldOut <- as.character(df$tx_out[idxMappable])
 
-        ## translations
-        newOut <- mergeMap[oldOut]
+    ## translations
+    newOut <- mergeMap[oldOut]
 
-        ## replace old with new
-        df[idxMappable, "tx_out"] <- newOut
-    }
-    ## did we fail to converge?
-    stopifnot(iteration <= MAXITERS && !any(df$tx_out %in% df$tx_in))
+    ## replace old with new
+    df[idxMappable, "tx_out"] <- newOut
+  }
+  ## did we fail to converge?
+  stopifnot(iteration <= MAXITERS && !any(df$tx_out %in% df$tx_in))
 
-    df
+  df
 }