#' Estimate Parameters From Real Datasets by dyntoy
#'
#' This function is used to estimate useful parameters from a real dataset by
#' using \code{infer_trajectory} function in dynwrap package.
#'
#' @param ref_data A count matrix. Each row represents a gene and each column
#' represents a cell.
#' @param verbose Logical.
#' @param seed An integer of a random seed.
#' @param other_prior A list with names of certain parameters. Some methods need
#' extra parameters to execute the estimation step, so you must input them. In
#' simulation step, the number of cells, genes, groups, batches, the percent of
#' DEGs are usually customed, so before simulating a dataset you must point it out.
#' See `Details` below for more information.
#' @importFrom dynwrap infer_trajectory wrap_expression add_grouping
#' @return A list contains the estimated parameters and the results of execution
#' detection.
#' @export
#' @details
#' In dyntoy, users can input cell group information if it is available. If cell
#' group information is not provided, the procedure will detect cell groups by
#' kmeans automatically.
#' See `Examples` for more instructions.
#'
#' @references
#' Github URL: <https://github.com/dynverse/dyntoy>
#' @examples
#' \dontrun{
#' ref_data <- simmethods::data
#'
#' estimate_result <- simmethods::dyntoy_estimation(
#' ref_data = ref_data,
#' other_prior = NULL,
#' verbose = TRUE,
#' seed = 111
#' )
#'
#' ## estimation with cell group information
#' group_condition <- paste0("Group", as.numeric(simmethods::group_condition))
#' estimate_result <- simmethods::dyntoy_estimation(
#' ref_data = ref_data,
#' other_prior = list(group.condition = group_condition),
#' verbose = TRUE,
#' seed = 111
#' )
#' }
#'
dyntoy_estimation <- function(ref_data,
verbose = FALSE,
other_prior = NULL,
seed){
##############################################################################
#### Check ###
##############################################################################
if(!requireNamespace("tislingshot", quietly = TRUE)){
message("tislingshot is not installed on your device...")
message("Installing tislingshot...")
devtools::install_github("dynverse/ti_slingshot/package/")
}
if(!requireNamespace("NbClust", quietly = TRUE)){
message("NbClust is not installed on your device...")
message("Installing NbClust...")
utils::install.packages("NbClust")
}
if(!is.matrix(ref_data)){
ref_data <- as.matrix(ref_data)
}
if(is.null(other_prior[["group.condition"]])){
message("Performing k-means and determin the best number of clusters...")
clust <- NbClust::NbClust(data = t(ref_data),
distance = 'euclidean',
min.nc = 2,
max.nc = sqrt(nrow(t(ref_data))),
method = "kmeans",
index = "dunn")
other_prior[["group.condition"]] <- paste0("Group", clust[["Best.partition"]])
}
ref_data <- dynwrap::wrap_expression(counts = t(ref_data),
expression = log2(t(ref_data) + 1))
## Add group
if(verbose){
message("Add grouping to data...")
}
ref_data <- dynwrap::add_grouping(dataset = ref_data,
grouping = other_prior[["group.condition"]])
##############################################################################
#### Estimation ###
##############################################################################
if(verbose){
message("Estimating parameters using dyntoy")
}
# Estimation
estimate_detection <- peakRAM::peakRAM(
estimate_result <- dynwrap::infer_trajectory(dataset = ref_data,
method = tislingshot::ti_slingshot(),
parameters = NULL,
give_priors = NULL,
seed = seed,
verbose = verbose)
)
estimate_result <- list(estimate_result = estimate_result,
data_dim = c(length(ref_data[["feature_ids"]]),
length(ref_data[["cell_ids"]])))
##############################################################################
#### Ouput ###
##############################################################################
estimate_output <- list(estimate_result = estimate_result,
estimate_detection = estimate_detection)
return(estimate_output)
}
#' Simulate Datasets by dyntoy
#'
#' @param parameters A object generated by [dynwrap::infer_trajectory()]
#' @param other_prior A list with names of certain parameters. Some methods need
#' extra parameters to execute the estimation step, so you must input them. In
#' simulation step, the number of cells, genes, groups, batches, the percent of
#' DEGs are usually customed, so before simulating a dataset you must point it out.
#' See `Details` below for more information.
#' @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.
#' @export
#' @details
#' In dyntoy, users can only set `nCells` and `nGenes` to specify the number of cells and genes in the
#' simulated dataset. See `Examples` for instructions.
#'
#' @references
#' Github URL: <https://github.com/dynverse/dyntoy>
#' @examples
#' \dontrun{
#' ref_data <- simmethods::data
#'
#' ## estimation with cell group information
#' group_condition <- paste0("Group", as.numeric(simmethods::group_condition))
#' estimate_result <- simmethods::dyntoy_estimation(
#' ref_data = ref_data,
#' other_prior = list(group.condition = group_condition),
#' verbose = TRUE,
#' seed = 111
#' )
#'
#' # 1) Simulate with default parameters
#' simulate_result <- simmethods::dyntoy_simulation(
#' parameters = estimate_result[["estimate_result"]],
#' other_prior = NULL,
#' return_format = "list",
#' verbose = TRUE,
#' seed = 111
#' )
#' ## counts
#' counts <- simulate_result[["simulate_result"]][["count_data"]]
#' dim(counts)
#'
#' # 2) 2000 cells and 5000 genes
#' simulate_result <- simmethods::dyntoy_simulation(
#' parameters = estimate_result[["estimate_result"]],
#' other_prior = list(nCells = 2000,
#' nGenes = 5000),
#' return_format = "list",
#' verbose = TRUE,
#' seed = 111
#' )
#'
#' ## counts
#' counts <- simulate_result[["simulate_result"]][["count_data"]]
#' dim(counts)
#' }
#'
dyntoy_simulation <- function(parameters,
other_prior = NULL,
return_format,
verbose = FALSE,
seed
){
##############################################################################
#### Environment ###
##############################################################################
if(!requireNamespace("dyntoy", quietly = TRUE)){
message("dyntoy is not installed on your device")
message("Installing dyntoy...")
devtools::install_github("dynverse/dyntoy")
}
##############################################################################
#### Check ###
##############################################################################
model <- parameters[["estimate_result"]][["milestone_network"]]
# nCells
if(!is.null(other_prior[["nCells"]])){
nCells <- other_prior[["nCells"]]
}else{
nCells <- parameters[["data_dim"]][2]
}
# nGenes
if(!is.null(other_prior[["nGenes"]])){
nGenes <- other_prior[["nGenes"]]
}else{
nGenes <- parameters[["data_dim"]][1]
}
# de.prob
if(!is.null(other_prior[["de.prob"]])){
de.prob <- other_prior[["de.prob"]]
}else{
de.prob <- 0.1
}
# Return to users
message(paste0("nCells: ", nCells))
message(paste0("nGenes: ", nGenes))
message(paste0("de.prob: ", de.prob))
##############################################################################
#### Simulation ###
##############################################################################
if(verbose){
message("Simulating datasets using dyntoy")
}
# Seed
set.seed(seed)
# Estimation
simulate_detection <- peakRAM::peakRAM(
simulate_result <- dyntoy::generate_dataset(model = model,
num_cells = nCells,
num_features = nGenes,
differentially_expressed_rate = de.prob)
)
##############################################################################
#### Format Conversion ###
##############################################################################
de_gene <- simulate_result[["tde_overall"]][["differentially_expressed"]]
simulate_result <- t(as.matrix(simulate_result[["counts"]]))
colnames(simulate_result) <- paste0("Cell", 1:ncol(simulate_result))
rownames(simulate_result) <- paste0("Gene", 1:nrow(simulate_result))
## col_data
col_data <- data.frame("cell_name" = colnames(simulate_result))
## row_data
row_data <- data.frame("gene_name" = rownames(simulate_result),
"de_gene" = ifelse(de_gene, "yes", "no"))
# Establish SingleCellExperiment
simulate_result <- SingleCellExperiment::SingleCellExperiment(list(counts = simulate_result),
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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.