saezlab / mistyR

Multiview Intercellular SpaTial modeling framework
https://saezlab.github.io/mistyR/
GNU General Public License v3.0
43 stars 12 forks source link

Can I add z axis to the location?Consider three-dimensional coordinates. #13

Closed xiangrong7 closed 3 years ago

xiangrong7 commented 3 years ago

Thanks for you developed a great tool. I have spatial datasets.Their coordinates are three-dimensional.Can I use the three-dimensional to caluate the distance? Like this:

# in your code 
geometry <- GetTissueCoordinates(visium.slide,
    cols = c("row", "col"), scale = NULL
  )
# I want to this
geometry  = data.frame("row" = x, "col"=y,"z"=z)
jtanevski commented 3 years ago

If you plan to use the view generation functions that are implemented in mistyR then add_paraview will use all available dimensions in the positions data frame to calculate the distances. However, add_juxtaview will pass the positions as first argument to deldir::deldir(). Note that, unfortunately, the Delauney triangulation with deldir doesn't work with 3D coordinates.

xiangrong7 commented 3 years ago

If you plan to use the view generation functions that are implemented in mistyR then add_paraview will use all available dimensions in the positions data frame to calculate the distances. However, add_juxtaview will pass the positions as first argument to deldir::deldir(). Note that, unfortunately, the Delauney triangulation with deldir doesn't work with 3D coordinates.

How about this function geometry::delaunayn() ?Sorry, I don't know much about this.

jtanevski commented 3 years ago

Looks like a good alternative. I need to take a closer look and might include it in a future version. For now, you can consider creating your custom view using this function. Once you have the sum or mean expression for all markers for each spatial unit you can transform it into a misty view and add it to your workflow see add_views().

xiangrong7 commented 3 years ago

Thanks, I will try it!

xiangrong7 commented 3 years ago

I changed add_juxtaview function slightly,including geometry::delaunayn() and get_neighbors().On 2D,the neighbors are same(deldir::deldir() vs geometry::delaunayn()) 2D get_neighbors function is this:

get_neighbors <- function(ddobj, id) {
  dplyr::union_all(
    ddobj[,2][which(ddobj[,1] == id)],
    ddobj[,3][which(ddobj[,1] == id)],
    ddobj[,1][which(ddobj[,2] == id)],
    ddobj[,3][which(ddobj[,2] == id)],
    ddobj[,1][which(ddobj[,3] == id)],
    ddobj[,2][which(ddobj[,3] == id)],
  )%>%unique()
}

3D get_neighbors function is this:

get_neighbors <- function(ddobj, id) {
  dplyr::union_all(
    ddobj[,2][which(ddobj[,1] == id)],
    ddobj[,3][which(ddobj[,1] == id)],
    ddobj[,4][which(ddobj[,1] == id)],
    ddobj[,1][which(ddobj[,2] == id)],
    ddobj[,3][which(ddobj[,2] == id)],
    ddobj[,4][which(ddobj[,2] == id)],
    ddobj[,1][which(ddobj[,3] == id)],
    ddobj[,2][which(ddobj[,3] == id)],
    ddobj[,4][which(ddobj[,3] == id)],
    ddobj[,1][which(ddobj[,4] == id)],
    ddobj[,2][which(ddobj[,4] == id)],
    ddobj[,3][which(ddobj[,4] == id)],
  )%>%unique()
}
add_juxtaview_3D <- function(current.views, positions, neighbor.thr = 15,
                          cached = FALSE, verbose = TRUE) {
  expr <- current.views[["intraview"]][["data"]]

  cache.location <- R.utils::getAbsolutePath(paste0(
    ".misty.temp", .Platform$file.sep,
    current.views[["misty.uniqueid"]]
  ))

  juxta.cache.file <- paste0(
    cache.location, .Platform$file.sep,
    "juxta.view.", neighbor.thr, ".rds"
  )

  if (cached && !dir.exists(cache.location)) {
    dir.create(cache.location, recursive = TRUE, showWarnings = TRUE)
  }

  if (cached & file.exists(juxta.cache.file)) {
    juxta.view <- readr::read_rds(juxta.cache.file)
  }
  else {
    if (verbose) message("\nComputing triangulation")
    delaunay <- geometry::delaunayn(as.data.frame(positions))

    if (verbose) message("\nGenerating juxtaview")
    juxta.view <- seq(nrow(expr)) %>% furrr::future_map_dfr(function(cid) {
      alln <- get_neighbors(delaunay, cid)
      # suboptimal placement of dists, but makes conflict if out of scope
      # probably due to lazy evaluations
      dists <- distances::distances(as.data.frame(positions))
      actualn <- alln[which(dists[alln, cid] <= neighbor.thr)]
      data.frame(t(colSums(expr[actualn, ])))
    }, .progress = verbose)

    if (cached) readr::write_rds(juxta.view, juxta.cache.file)
  }

  return(current.views %>% add_views(create_view(
    paste0("juxtaview.", neighbor.thr),
    juxta.view,
    paste0("juxta.", neighbor.thr)
  )))
}