Merge branch 'dev' into 'master'

prepare release v0.6.8

See merge request tron/bnt_neoants/splice2neo!126

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: splice2neo
 Title: Aberrant splice junction analysis with associated mutations
-Version: 0.6.7
+Version: 0.6.8
 Authors@R: 
     c(
     person(given = "Jonas",
@@ -41,6 +41,7 @@ Suggests:
     ggbio,
     rmarkdown,
     BSgenome,
+    liftOver,
 Config/testthat/edition: 3
 Language: en-US
 Imports: 

NAMESPACE

@@ -34,6 +34,7 @@ export(is_in_rnaseq)
 export(junc2breakpoint)
 export(junc_to_gr)
 export(leafcutter_transform)
+export(liftover_junc_id)
 export(map_requant)
 export(modify_tx)
 export(parse_cispliceai_thresh)

NEWS.md

@@ -1,3 +1,9 @@
+# splice2neo 0.6.8
+
+- Add function `liftover_junc_id() to liftOver splice junctions to other 
+ reference genomes or personalized genomes
+
+
 # splice2neo 0.6.7
 
 ## Added:

R/canonical_junctions.R

@@ -21,7 +21,7 @@
 #'
 #' @examples
 #' gtf_file <- system.file("extdata","GTF_files","Aedes_aegypti.partial.gtf",
-#'   package="GenomicFeatures")
+#'   package="splice2neo")
 #'
 #' tx <- parse_gtf(gtf_file)
 #'

R/globals.R

@@ -264,6 +264,7 @@ utils::globalVariables(
     "within_interval_3",
     "within_interval_left",
     "within_interval_right",
-    "WT_protein_length_difference"
+    "WT_protein_length_difference",
+    "junc_id_lifted_lst"
   )
 )

R/liftover_junc_id.R

@@ -0,0 +1,73 @@
+
+
+#' LiftOver junction IDs using the liftOver tool
+#'
+#' This can be used to liftOver junctions to personalized genome coordinates.
+#'
+#' @param junc_df a data.frame with at least one column `junc_id` containing junction IDs
+#' @param chain_file path to a chain file for the UCSC liftOver tool. See also \code{\link[rtracklayer]{liftOver}}
+#' @return a data.frame like the input `junc_df` with the following additional columns:
+#'   - `junc_id_lifted_is_unique` a logical vector indicating if the liftOver was successful and unique (1-to-1 correspondence).
+#'   - `junc_id_lifted_collapsed` a character vector with the lifted junction IDs.
+#'   Multiple IDs are separated by `|`.
+#'   NA represent junc_ids that could not be lifted.
+#'   - `junc_id_lifted` a character vector with a unique lifted junction IDs.
+#'   Potentially multiple lifted IDs are combined by the minmal start and maximal
+#'   end coordinate. NA represent junc_ids that could not be lifted.
+#'
+#' @examples
+#'
+#' chain_file = system.file(package="liftOver", "extdata", "hg38ToHg19.over.chain")
+#' junc_df <- toy_junc_df
+#' liftover_junc_id(junc_df, chain_file)
+#'
+#'@export
+liftover_junc_id <- function(junc_df, chain_file){
+
+
+  # convert junc_id into GenomicRanges object
+  gr <- junc_to_gr(junc_df$junc_id)
+
+  # read chain file
+  chain <- rtracklayer::import.chain(chain_file)
+
+  # Lift over junc ranges to genomic position in the personalized genome
+  junc_id_lifted_grl  <- rtracklayer::liftOver(gr, chain)
+
+  junc_df %>%
+    dplyr::mutate(
+
+      # add lifted junctions ranges as a list
+      junc_id_lifted_lst = as.list(junc_id_lifted_grl),
+
+      # check if liftOver was successful and unique
+      junc_id_lifted_is_unique = purrr::map_lgl(junc_id_lifted_lst, ~ length(.x) == 1),
+
+      # covert back to junc_id and collapse multiple IDs with `|`, and replace empty junc_id with NA
+      junc_id_lifted_collapsed = junc_id_lifted_lst %>%
+        purrr::map(as.character) %>%
+        purrr::map_chr(stringr::str_c, collapse = "|") %>%
+        dplyr::na_if(""),
+
+      # build unique lifted junctions by minimal start and maximal end
+      # coordinate of potentially multiple lifted junction ranges
+      junc_id_lifted = junc_id_lifted_lst %>%
+
+        # take min start and max end across potentially multiple junctions
+        purrr::map(base::range) %>%
+
+        # convert to junc_id and replace empty entries with NA
+        purrr::map(as.character) %>%
+        purrr::map_if(is_empty, ~ NA_character_) %>%
+        unlist()
+
+    ) %>%
+
+    # remove temporary column
+    dplyr::select(-junc_id_lifted_lst)
+
+}
+
+
+
+

R/parse_gtf.R

@@ -7,7 +7,7 @@
 #' @return A \code{\link[GenomicRanges]{GRangesList}} of exons grouped by transcripts
 #'
 #' @examples
-#' gff_file <- system.file("extdata","GFF3_files","a.gff3",package="GenomicFeatures")
+#' gff_file <- system.file("extdata","GFF3_files","a.gff3",package="splice2neo")
 #'
 #' parse_gtf(gff_file)
 #'

_pkgdown.yml

@@ -124,6 +124,7 @@ reference:
   - generate_junction_id
   - breakpoint2junc
   - junc2breakpoint
+  - liftover_junc_id
 
 - title: "Merge and filter"
   desc: "Functions annotate and filter for canonical splicing"

tests/testthat/test-canonical_junctions.R

@@ -2,7 +2,7 @@ test_that("canonical_junctions works on test gtf", {
 
 
   gtf_file <- system.file("extdata","GTF_files","Aedes_aegypti.partial.gtf",
-                          package="GenomicFeatures")
+                          package="splice2neo")
 
   tx <- parse_gtf(gtf_file, format = "gtf")
 

tests/testthat/test-liftover_junc_id.R

@@ -0,0 +1,32 @@
+test_that("liftover_junc_id works on toy example data", {
+
+  chain_file = system.file(package="liftOver", "extdata", "hg38ToHg19.over.chain")
+
+  junc_df <- toy_junc_df
+
+  junc_df_lifted <- liftover_junc_id(junc_df, chain_file)
+
+  expect_equal(nrow(junc_df_lifted), nrow(junc_df))
+  expect_true(all(c("junc_id_lifted_is_unique", "junc_id_lifted") %in% names(junc_df_lifted)))
+
+})
+
+
+
+test_that("liftover_junc_id works with non-unique mappings", {
+
+  chain_file = system.file(package="liftOver", "extdata", "hg38ToHg19.over.chain")
+
+  junc_df <- toy_junc_df %>%
+    head(3) %>%
+    add_row(junc_id = "chr2:10000000-10050000:-") %>%
+    add_row(junc_id = "chr7:123-456:+")
+
+  junc_df_lifted <- liftover_junc_id(junc_df, chain_file)
+
+  expect_equal(nrow(junc_df_lifted), nrow(junc_df))
+  expect_false(all(junc_df_lifted$junc_id_lifted_is_unique))
+  expect_true(all(junc_df_lifted$junc_id_lifted_is_unique[1:3]))
+  expect_true(any(is.na(junc_df_lifted$junc_id_lifted)))
+
+})

tests/testthat/test-parse_gtf.R

@@ -1,6 +1,6 @@
 test_that("parse_gtf works with GFF file", {
 
-  gff_file <- system.file("extdata","GFF3_files","a.gff3",package="GenomicFeatures")
+  gff_file <- system.file("extdata","GFF3_files","a.gff3",package="splice2neo")
 
   grl <- parse_gtf(gff_file, format = "gff3")
 
@@ -10,7 +10,7 @@ test_that("parse_gtf works with GFF file", {
 test_that("parse_gtf works with GTF file", {
 
   gtf_file <- system.file("extdata","GTF_files","Aedes_aegypti.partial.gtf",
-                         package="GenomicFeatures")
+                         package="splice2neo")
 
   grl <- parse_gtf(gtf_file, format = "gtf")