R/28-SRTsim.R

In duohongrui/simmethods: A collection of 40+ simulation methods for singlce-cell RNA seqencing data

#' Estimate Parameters From Real Datasets by SRTsim
#'
#' This function is used to estimate useful parameters from a real dataset by
#' using `srtsim_fit` function in SRTsim package.
#'
#' @param ref_data A count matrix. Each row represents a gene and each column
#' represents a cell/spot.
#' @param verbose Logical.
#' @param other_prior A list with names of certain parameters. Some methods need
#' extra parameters to execute the estimation step, so you should input them.
#' See `Details` below for more information.
#' @param seed An integer of a random seed.
#' @return A list contains the estimated parameters and the results of execution
#' detection.
#' @export
#'
#' @details
#' In SRTsim, users should input the spatial coordinates of spots or cells in the
#' spatial transcriptome data.
#'
#' 1. spatial.x. The x-axis coordinates.
#' 2. spatial.y. The y-axis coordinates.
#' 3. group.condition. Users can also input a numeric vector to specify the tissues or domains that each cell/spot comes from,
#' like `other_prior = list(group.condition = the numeric vector)`.
#' 4. sim_scheme. If the labels of cells/spots represents the tissue layers or slices, `sim_scheme` should be set as `tissue`.
#' `sim_scheme` is set as `domain` by default.
#'
#' @references
#' Zhu J, Shang L, Zhou X. SRTsim: spatial pattern preserving simulations for spatially resolved transcriptomics. Genome Biology, 2023, 24(1): 39. <https://doi.org/10.1186/s13059-023-02879-z>
#'
#' URL: <https://cran.r-project.org/web/packages/SRTsim/index.html>
#'
SRTsim_estimation <- function(ref_data,
                              verbose = FALSE,
                              other_prior = NULL,
                              seed
){
  ##############################################################################
  ####                               Check                                   ###
  ##############################################################################
  if(!is.matrix(ref_data)){
    ref_data <- as.matrix(ref_data)
  }
  ##############################################################################
  ####                            Estimation                                 ###
  ##############################################################################
  if(verbose){
    message("Estimating parameters using SRTsim")
  }

  ### image col and row
  if(is.null(other_prior[["spatial.x"]]) |
     is.null(other_prior[["spatial.y"]])){
    stop("Please input the spatial coordinates of spots in other_prior")
  }
  ### scheme
  if(is.null(other_prior[["sim_schem"]])){
    other_prior[["sim_schem"]] <- "domain"
  }
  ### group.condition
  if(is.null(other_prior[["group.condition"]])){
    other_prior[["group.condition"]] <- rep("A", ncol(ref_data))
  }else{
    other_prior[["group.condition"]] <- LETTERS[other_prior[["group.condition"]]]
  }
  ## create a SRT object
  metadata <- data.frame("x" = other_prior[["spatial.x"]],
                         "y" = other_prior[["spatial.y"]],
                         "label" = other_prior[["group.condition"]])
  rownames(metadata) <- colnames(ref_data)
  simSRT <- SRTsim::createSRT(count_in = ref_data,
                              loc_in = metadata)
  other_prior[["simsrt"]] <-simSRT

  # Seed
  set.seed(seed)
  # Estimation
  estimate_detection <- peakRAM::peakRAM(
    estimate_result <- SRTsim::srtsim_fit(simsrt = simSRT,
                                          marginal = ifelse(is.null(other_prior[["auto_choose"]]),
                                                            "auto_choose",
                                                            other_prior[["auto_choose"]]),
                                          sim_scheme = other_prior[["sim_scheme"]],
                                          min_nonzero_num = ifelse(is.null(other_prior[["min_nonzero_num"]]),
                                                                   2,
                                                                   other_prior[["min_nonzero_num"]]),
                                          maxiter = ifelse(is.null(other_prior[["maxiter"]]),
                                                           500,
                                                           other_prior[["maxiter"]]))
  )
  ##############################################################################
  ####                           Ouput                                       ###
  ##############################################################################
  estimate_output <- list(estimate_result = estimate_result,
                          estimate_detection = estimate_detection)
  return(estimate_output)
}


#' Simulate Datasets by SRTsim
#'
#' This function is used to simulate datasets from learned parameters by `srtsim_count`
#' function in SRTsim package.
#'
#' @param parameters A object generated by [SRTsim::srtsim_fit()]
#' @param other_prior A list with names of certain parameters in SRTsim.
#' @param return_format A character. Alternative choices: list, SingleCellExperiment,
#' Seurat, h5ad. If you select `h5ad`, you will get a path where the .h5ad file saves to.
#' @param verbose Logical. Whether to return messages or not.
#' @param seed A random seed.
#' @details
#'
#' For more customed parameters in SRTsim, please check [SRTsim::srtsim_fit()].
#' For detailed information about SRTsim, go to <https://cran.r-project.org/web/packages/SRTsim/index.html>.
#'
#' @export
#' @references
#' Zhu J, Shang L, Zhou X. SRTsim: spatial pattern preserving simulations for spatially resolved transcriptomics. Genome Biology, 2023, 24(1): 39. <https://doi.org/10.1186/s13059-023-02879-z>
#'
#' URL: <https://cran.r-project.org/web/packages/SRTsim/index.html>
#'
SRTsim_simulation <- function(parameters,
                              other_prior = NULL,
                              return_format,
                              verbose = FALSE,
                              seed
){
  ##############################################################################
  ####                               Check                                   ###
  ##############################################################################
  assertthat::assert_that(class(parameters) == "simSRT")
  other_prior[["simsrt"]] <- parameters
  other_prior[["verbose"]] <- verbose

  # Return to users
  message(paste0("nSpots: ", ncol(parameters@refCounts)))
  message(paste0("nGenes: ", nrow(parameters@refCounts)))
  message(paste0("nGroups: ", length(unique(parameters@refcolData$label))))
  ##############################################################################
  ####                            Simulation                                 ###
  ##############################################################################
  if(verbose){
    message("Simulating datasets using SRTsim")
  }
  # Seed
  set.seed(seed)
  other_prior <- other_prior[intersect(names(other_prior), names(formals(SRTsim::srtsim_count)))]
  simulate_detection <- peakRAM::peakRAM(
    simulate_result <- BiocGenerics::do.call(SRTsim::srtsim_count, other_prior)
  )
  ##############################################################################
  ####                        Format Conversion                              ###
  ##############################################################################
  # counts
  counts <- as.matrix(simulate_result@simCounts)
  # col_data
  col_data <- as.data.frame(simulate_result@simcolData)
  colnames(col_data)[3] <- "group"
  col_data$"cell_name" <- rownames(col_data)

  # row_data
  row_data <- data.frame("gene_name" = rownames(counts))

  # Establish SingleCellExperiment
  simulate_result <- SingleCellExperiment::SingleCellExperiment(list(counts = counts),
                                                                colData = col_data,
                                                                rowData = row_data)
  simulate_result <- simutils::data_conversion(SCE_object = simulate_result,
                                               return_format = return_format)

  ##############################################################################
  ####                           Ouput                                       ###
  ##############################################################################
  simulate_output <- list(simulate_result = simulate_result,
                          simulate_detection = simulate_detection)
  return(simulate_output)
}