add rna velocity as possible input, fixes #112

R/method_process_definition.R

@@ -40,7 +40,7 @@ definition <- function(
       input_id = c(definition$wrapper$input_required, definition$wrapper$input_optional, map_chr(definition$parameters$parameters, "id")),
       required = input_id %in% definition$wrapper$input_required,
       type = case_when(
-        input_id %in% c("counts", "expression", "rna_velocity") ~ "expression",
+        input_id %in% c("counts", "expression", "expression_projected") ~ "expression",
         input_id %in% priors$prior_id ~ "prior_information",
         TRUE ~ "parameter"
       )

R/wrap_add_attraction.R

@@ -0,0 +1,77 @@
+#' Calculate the attraction of cells to other cells using velocity
+#'
+#' @param current Current expression
+#' @param projected Projected expression based on RNA velocity
+#'
+#' @return Matrix containing the attraction ([-1, 1]) of each cell to the waypoint cells
+calculate_attraction <- function(
+  current,
+  projected,
+  cells = colnames(projected),
+  n_waypoints = 50,
+  k = ceiling(n_waypoints / 2),
+  corr_sigma = 1,
+  debug = FALSE
+) {
+  assertthat::assert_that(nrow(current) == nrow(projected))
+  assertthat::assert_that(ncol(current) == ncol(projected))
+
+  # select waypoint cells
+  n_waypoints <- min(n_waypoints, length(cells))
+  k <- min(n_waypoints, k)
+  waypoint_cells <- sample(cells, n_waypoints)
+
+  # prepare for correlation calculation
+  # em = expression
+  # nd = velocity
+
+  em <- as.matrix(current)
+  ccells <- cells
+  em <- em[, ccells]
+  nd <- as.matrix(projected[, ccells] - current[, ccells])
+  cgenes <- intersect(rownames(em), rownames(nd))
+  nd <- nd[cgenes, ]
+  em <- em[cgenes, ]
+
+  # calculate correlation
+  # this is an adapted version of colDeltaCorLog10 with waypoints
+  transfo <- function(x) (log10(abs(x) + 1) * sign(x))
+
+  nd2 <- transfo(nd)
+
+  waypoint <- waypoint_cells[[1]]
+  cell <- cells[[1]]
+
+  emw <- em[, waypoint_cells]
+
+  # cors <- map(cells, function(cell) {
+  #   print(cell)
+  #   diff <- transfo(emw - em[, cell])
+  #   cors <- cor(diff, nd2[, cell])
+  #   colnames(cors) <- cell
+  #   cors[is.na(cors)] <- 0
+  #   cors
+  # })
+  # attraction <- do.call(cbind, cors)
+  # attraction
+
+
+  cors <- map(waypoint_cells, function(waypoint) {
+    print(waypoint)
+    diff <- transfo(emw[, waypoint] - em)
+    cors <- pcor(diff, nd2)
+    rownames(cors) <- waypoint
+    cors[is.na(cors)] <- 0
+    cors
+  })
+  attraction <- do.call(rbind, cors)
+  colnames(attraction) <- colnames(em)
+  attraction
+}
+
+
+
+pcor <- function(x, y = x, method = "pearson", use = "everything") {
+  assertthat::assert_that(ncol(x) == ncol(y));
+  purrr::map_dbl(seq_len(ncol(x)), ~ cor(x[,.], y[,.], method = method, use = use))
+}

R/wrap_add_expression.R

@@ -3,7 +3,7 @@
 #' @inheritParams common_param
 #' @param counts The counts with genes in columns and cells in rows
 #' @param expression The normalised expression values with genes in columns and cells in rows
-#' @param rna_velocity RNA velocity feature WIP
+#' @param expression_projected Projected expression using RNA velocity
 #' @param feature_info Optional meta-information of the features, a data.frame with at least feature_id as column
 #' @param ... extra information to be stored in the dataset
 #' @param expression_source The source of expression, can be "counts", "expression", an expression matrix, or another dataset which contains expression
@@ -19,45 +19,17 @@ add_expression <- function(
   counts,
   expression,
   feature_info = NULL,
-  rna_velocity = NULL,
+  expression_projected = NULL,
   ...
 ) {
   testthat::expect_true(is_data_wrapper(dataset))
 
   assert_that(!(is.null(counts) && is.null(expression)), msg = "counts and expression can't both be NULL")
 
-  if (!is.null(counts)) {
-    if (is.matrix(counts)) {
-      counts <- Matrix::Matrix(counts, sparse = TRUE)
-    }
-    if (is_sparse(counts) && !"dgCMatrix" %in% class(counts)) {
-      counts <- as(counts, "dgCMatrix")
-    }
-    assert_that(
-      "dgCMatrix" %in% class(counts),
-      identical(rownames(counts), dataset$cell_ids)
-    )
-  }
-
-  if (!is.null(expression)) {
-    if (is.matrix(expression)) {
-      expression <- Matrix::Matrix(expression, sparse = TRUE)
-    }
-    assert_that(
-      "dgCMatrix" %in% class(expression),
-      identical(rownames(expression), dataset$cell_ids)
-    )
-  }
-
-  if (!is.null(rna_velocity)) {
-    if (is.matrix(rna_velocity)) {
-      rna_velocity <- Matrix::Matrix(rna_velocity, sparse = TRUE)
-    }
-    assert_that(
-      "dgCMatrix" %in% class(rna_velocity),
-      identical(rownames(rna_velocity), dataset$cell_ids)
-    )
-  }
+  # convert expression if needed
+  counts <- convert_expression(counts, dataset$cell_ids)
+  expression <- convert_expression(expression, dataset$cell_ids)
+  expression_projected <- convert_expression(expression_projected, dataset$cell_ids)
 
   if (!is.null(feature_info)) {
     assert_that(
@@ -75,7 +47,7 @@ add_expression <- function(
     "dynwrap::with_expression",
     counts = counts,
     expression = expression,
-    rna_velocity = rna_velocity,
+    expression_projected = expression_projected,
     feature_info = feature_info,
     ...
   )
@@ -129,7 +101,7 @@ wrap_expression <- function(
   counts,
   cell_info = NULL,
   feature_info = NULL,
-  rna_velocity = NULL,
+  expression_projected = NULL,
   ...
 ) {
   cell_ids <- rownames(expression) %||% rownames(counts)
@@ -147,8 +119,24 @@ wrap_expression <- function(
     add_expression(
       counts = counts,
       expression = expression,
-      rna_velocity = rna_velocity,
+      expression_projected = expression_projected,
       feature_ids = feature_ids,
       feature_info = feature_info
     )
 }
+
+
+convert_expression <- function(x, cell_ids) {
+  if (!is.null(x)) {
+    if (is.matrix(x)) {
+      x <- Matrix::Matrix(x, sparse = TRUE)
+    }
+    if (is_sparse(x) && !"dgCMatrix" %in% class(x)) {
+      x <- as(x, "dgCMatrix")
+    }
+    assert_that(
+      "dgCMatrix" %in% class(x),
+      identical(rownames(x), cell_ids)
+    )
+  }
+}

data-raw/allowed_inputs_outputs.R

@@ -53,7 +53,7 @@ allowed_inputs <- tribble(
   ~input_id, ~description,
   "expression", "Expression matrix (sparse)",
   "counts", "Raw counts matrix (sparse)",
-  "rna_velocity", "RNA velocity matrix (sparse)"
+  "expression_projected", "Projected expression matrix based on RNA velocity (sparse)"
 ) %>% bind_rows(
   priors %>% select(input_id = prior_id, description = description)
 )

tests/testthat/test-metawrap_get_output_processor.R

@@ -12,6 +12,6 @@ test_that("get_output_processor works correctly", {
   expect_is(op_expression, "list")
   expect_equal(op_expression$processor, add_expression)
   expect_equal(sort(op_expression$required_args), c("counts", "expression"))
-  expect_equal(sort(op_expression$optional_args), c("feature_info", "rna_velocity"))
+  expect_equal(sort(op_expression$optional_args), c("feature_info", "expression_projected"))
   expect_equal(sort(op_expression$args), sort(c(op_expression$required_args, op_expression$optional_args)))
 })