Markov Chain Construction and Probability Calculation

This function preprocesses data, constructs a Markov chain, and calculates transition probabilities based on pseudotime information.

Usage

markovProbability(
  milo,
  diffusionmap,
  terminal_state = NULL,
  root_cell,
  knn = 30L,
  diffusiontime = NULL,
  pseudotime_key = "pseudotime",
  scale_components = TRUE,
  num_waypoints = 500
)

Arguments

milo: A Milo or SingleCellExperiment object. This object should have pseudotime stored in colData, which will be used to calculate probabilities. If pseudotime is available in milo, it takes precedence over the value provided through the diffusiontime parameter.
diffusionmap: A DiffusionMap object corresponding to the milo object. Used for Markov chain construction.
terminal_state: Integer. The index of the terminal state in the Markov chain.
root_cell: Integer. The index of the root state in the Markov chain.
knn: Integer. The number of nearest neighbors for graph construction. Default is 30L.
diffusiontime: Numeric vector. If pseudotime is not stored in milo, this parameter can be used to provide pseudotime values to the function.
pseudotime_key: Character. The name of the column in colData that contains the inferred pseudotime.
scale_components: Logical. If TRUE, the components will be scaled before constructing the Markov chain. Default is FALSE.
num_waypoints: Integer. The number of waypoints to sample when constructing the Markov chain. Default is 500L.

Value

milo or SinglCellExperiment object with pseudotime, probabilities in its colData

Examples

data(sce_vdj)
sce_vdj <- setupVdjPseudobulk(sce_vdj,
    already.productive = FALSE,
    allowed_chain_status = c("Single pair", "Extra pair")
)
#> Checking productivity from productive_abT_VDJ, productive_abT_VJ ...
#> 7279 of cells filtered
#> checking allowed chain status...
#> 12 of cells filtered
#> VDJ data extraction begin:
#> Parameter extract_cols do not provided, automatically geneterate colnames for extraction.
#> Detect whether colData v_call_abT_VDJ, d_call_abT_VDJ, j_call_abT_VDJ, v_call_abT_VJ, j_call_abT_VJ already exist...
#> Extract main TCR from v_call_abT_VDJ, d_call_abT_VDJ, j_call_abT_VDJ, v_call_abT_VJ, j_call_abT_VJ ...
#> Complete.
#> Filtering cells from v_call_abT_VDJ_main, j_call_abT_VDJ_main, v_call_abT_VJ_main, j_call_abT_VJ_main ...
#> 63 of cells filtered
#> 2646 of cells remain.
# Build Milo Object
set.seed(100)
milo_object <- miloR::Milo(sce_vdj)
milo_object <- miloR::buildGraph(milo_object, k = 50, d = 20, reduced.dim = "X_scvi")
#> Constructing kNN graph with k:50
milo_object <- miloR::makeNhoods(milo_object, reduced_dims = "X_scvi", d = 20)
#> Checking valid object
#> Running refined sampling with reduced_dim

# Construct Pseudobulked VDJ Feature Space
pb.milo <- vdjPseudobulk(milo_object, col_to_take = "anno_lvl_2_final_clean")
pb.milo <- scater::runPCA(pb.milo, assay.type = "Feature_space")

# Define root and branch tips
pca <- t(as.matrix(SingleCellExperiment::reducedDim(pb.milo, type = "PCA")))
branch.tips <- c(189, 198) # which.min(pca[, 2]) and which.max(pca[, 2])
names(branch.tips) <- c("CD8+T", "CD4+T")
root <- 177 # which.max(pca[, 1])

# Construct Diffusion Map
dm <- destiny::DiffusionMap(t(pca), n_pcs = 50, n_eigs = 10)
#> 'as(<dsCMatrix>, "dgTMatrix")' is deprecated.
#> Use 'as(as(., "generalMatrix"), "TsparseMatrix")' instead.
#> See help("Deprecated") and help("Matrix-deprecated").
dif.pse <- destiny::DPT(dm, tips = c(root, branch.tips), w_width = 0.1)

# Markov Chain Construction
pb.milo <- markovProbability(
    milo = pb.milo,
    diffusionmap = dm,
    diffusiontime = dif.pse[[paste0("DPT", root)]],
    terminal_state = branch.tips,
    root_cell = root,
    pseudotime_key = "pseudotime"
)
#> Sampling and flocking waypoints...
#> Markov chain construction...
#> Computing fundamental matrix and absorption probabilities...
#> Project probabilites from waypoints to each pseudobulk...