R/6.2.2_ISS_cluster_spatial.R

In mashranga/MolDia: Single cell In-Situ sequencing (ISS) data analysis

######################################################################
##                         RCA spatial cluster                      ##
######################################################################
"ISS_cluster_spatial"
#' Spatial clustering for ISS by Markov random field model.
#' 
#' @description ISS spatial clustering by Markov random field with extended spatial parameters EM procedure. 
#'              This function is incomplete and needs further development to work properly. See details. 
#' 
#' @param data Input data in class MolDiaISS. Output of \link[MolDia]{readISS}.
#' @param neighbor.points Number for neighbouring points. Default is 4.
#' @param spatial.info Default is "phisical". Other value is "genetic".
#' 
#' @details This function has build on the work on the article Pettit, J at el (2014). Current function "ISS_cluster_spatial" 
#'          will only create 2 input for the actual algorithm veloped by Pettit, J at el (2014). After ruiing the function one 
#'          has to run the actual algorithm for spatial clustering. See example for application.
#' 
#' @references Pettit, J., Tomer, R., Achim, K., Richardson, S., Azizi, L., & Marioni, J. (2014). 
#'             Identifying Cell Types from Spatially Referenced Single-Cell Expression Datasets. PLoS 
#'             Computational Biology, 10(9). doi:10.1371/journal.pcbi.1003824
#' 
#' @examples 
#' ## Reading data
#' ex_data <- readISS(file = system.file("extdata", "Hypocampus_left.csv", package="MolDia"),
#'                    cellid = "CellId", centX = "centroid_x", centY = "centroid_y")
#' myres   <- ISS_cluster_spatial(ex_data, neighbor.points = 4)
#' 
#' ####################### Manual process to wun the algorithm 
#' ## Set working dirctory
#' mwd <- "c:/MRF"
#' 
#' ## Save required data in working directory 
#' main_data <- myres[[1]]
#' write.table(main_data, file = paste0(mwd,"/mdata.txt"), 
#'             col.names = FALSE, row.names = FALSE, quote = FALSE, sep = "\t")
#' neig_data <- myres[[2]]
#' write.table(neig_data, file = paste0(mwd,"/neig.txt"), 
#'             col.names = FALSE, row.names = FALSE, quote = FALSE,sep = "\t")
#'
#' ## For detail about the parameter and Download windows binary: https://github.com/jbogp/MRF_Platynereis_2014
#' ## Save the windows binary in the current working directory. Make sure the exe file has the administrator privilege            
#'           
#'  win_exe <- paste0(mwd,"/EM.exe")
#'  
#' ## Result folder name
#'  res_fold <- "myResultFolder"
#' 
#' ## Parameter set up
#'  param1   <- paste0(mwd,"/mdata.txt")
#'  param2   <- paste0(mwd,"/neig.txt")
#'  param3   <- "rand"
#'  param4   <- 0
#'  param5   <- 20
#'  param6   <- paste0(mwd,"/",res_fold)
#'  param7   <- "myResultFile"
#'  param8   <- 20
#' 
#' ## Run the function in windows CMD
#'  mycmd <- paste(win_exe, param1, param2, param3, param4, param5, param6, param7, param8)
#'  res   <- system(mycmd, intern = FALSE,invisible = FALSE)
#'
#' ## Read resulted cluster data data and convet it accordingly 
#' cell_clust <- read.csv(paste0(mwd,"/",res_fold,"/myResultFile.csv"), header = F)
#' cell_clust <- as.factor(cell_clust$V1)
#' names(cell_clust) <- rownames(ex_data@data)
#' 
#' ## Insert cluster information into main data object
#'  ex_data@cluster <- cell_clust
#' 
#' ## Map cluster data 
#' pp<- RCA_map(ex_data, what = "cluster"
#' 
#' ## Plot Spatial clustering 
#'  ex_data <- ISS_preprocess(ex_data)
#'  ex_data <- ISS_tsne(ex_data)
#'  ex_data <- ISS_marker(ex_data)
#'  
#'  ISS_map(ex_data, what = "cell", main = "All cells")
#'  ISS_map(ex_data, what = "cluster", main = "Spatial cluster: phisical distance", label.topgene = NULL)
#'  ISS_map(ex_data, what = "tsne", main = "Spatial in tsne: phisical distance", label.topgene = 2)
#'
#'@export
ISS_cluster_spatial <- function(data, neighbor.points = 4, spatial.info = "phisical")
{
  ######### Nearest neighbour cell
  ## Physical distance
  if(spatial.info == "phisical")
    {
    res1 <- fields::rdist(data@location)
    colnames(res1) <- rownames(data@location)
    rownames(res1) <- rownames(data@location)
    res1
    }
  
  ## Genetic distance
  if(spatial.info == "genetic")
    {
    res1 <- as.matrix(stats::dist(data@data))
    }
  
  ## Distance matrix
  res2 <- res1
  ## Find nearest neighbout cells
  
  res1           <- apply(res1,1,function(i)
    {
      kk<- sort(i, index.return = TRUE)$ix [2:(neighbor.points+1)]
    kk
    })
  res1 <- t(res1)
  rownames(res1) <- match(rownames(data@location),rownames(res1))
  rn <- apply(res1, 1, length)
  rownames(res1) <- rn
  res1 <- cbind(rn,res1)
  res1
  
  ## Binary gene expression
  bin_gene <- 1*(data@data > 0)
  rownames(bin_gene) <- match(rownames(data@location),rownames(bin_gene))
  rn <- rownames(bin_gene)
  bin_gene <- cbind(rn, bin_gene)
  bin_gene
  
  ## Return result
  res <- list(bin_gene,res1, res2)
  return(res)
  
}