Skip to content

Commit

Permalink
Format code (#11)
Browse files Browse the repository at this point in the history
* format the indentations

and also get rid of the seed

* more styling to be consistent
  • Loading branch information
zktuong authored Dec 3, 2024
1 parent 5488930 commit 718c5e8
Show file tree
Hide file tree
Showing 22 changed files with 758 additions and 669 deletions.
18 changes: 13 additions & 5 deletions R/check.R
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,14 @@
.class.check <- function(input, must) {
requireNamespace("rlang")
requireNamespace("methods")
if (is.null(input))
if (is.null(input)) {
return()
}
if (!methods::is(input, must)) {
rlang::abort(sprintf("The '%s' must be %s, not %s", as.character(substitute(input)), must, class(input)))
rlang::abort(sprintf(
"The '%s' must be %s, not %s", as.character(substitute(input)),
must, class(input)
))
}
}

Expand All @@ -21,11 +25,15 @@
#' @param must the type we need
.type.check <- function(input, must) {
requireNamespace("methods")
if (is.null(input))
if (is.null(input)) {
return()
}
if (!methods::is(input, must)) {
requireNamespace("rlang")
requireNamespace("BiocGenerics")
rlang::abort(sprintf("The '%s' must be %s, not %s", as.character(substitute(input)), must, BiocGenerics::type(input)))
rlang::abort(sprintf(
"The '%s' must be %s, not %s", as.character(substitute(input)),
must, BiocGenerics::type(input)
))
}
}
}
127 changes: 66 additions & 61 deletions R/construct.markov.chain.R
Original file line number Diff line number Diff line change
Expand Up @@ -7,65 +7,70 @@
#' @param waypoints integer vector, index of selected waypoint used to construct markov chain
#' @return transition matrix of the markov chain
.construct.markov.chain <- function(wp_data, knn., pseudotime, waypoints) {
message("Markov chain construction...")
pseudotime <- pseudotime[waypoints]
# construct kNN graph
requireNamespace("bluster")
nbrs <- bluster::makeKNNGraph(wp_data, k = knn.)
## calculate distance of each edge
distance_m <- as.matrix(stats::dist(wp_data))
requireNamespace("igraph")
igraph::E(nbrs)$weight <- sapply(igraph::E(nbrs), function(e) {
### Get the nodes connected by each edge
nodes <- igraph::ends(nbrs, e)
node1 <- as.numeric(nodes[1])
node2 <- as.numeric(nodes[2])
### Use the distance matrix to get the distance between the two nodes
distance_m[node1, node2]
})
## generate weighted adjacent matrix of this knn graph

KNN <- igraph::as_adjacency_matrix(nbrs, attr = "weight")
## generate the index of each neighbor
idx <- KNN@p
idx_seq <- mapply(seq, (idx + 1)[-length(idx)], idx[-1])
ind <- lapply(idx_seq, function(x) {
KNN@i[x] + 1
})
## select Standard deviation allowing for 'back' edges
adaptive.k <- min(c(floor(knn./3) - 1, 30))
dist_ <- lapply(idx_seq, function(y) {
KNN@x[y]
})
dist_sort <- lapply(dist_, sort, decreasing = TRUE)
adaptive.std <- vapply(dist_sort, "[", adaptive.k, FUN.VALUE = double(1))
# Directed graph construction pseudotime position of all the neighbors
traj_nbrs <- lapply(ind, function(x) {
pseudotime[x]
})
## Remove edges that move backwards in pseudotime except for edges that are
## within the computed standard deviation
requireNamespace("purrr")
rem_edges <- purrr::pmap(list(.x = traj_nbrs, .y = (pseudotime - adaptive.std), .z = ind),
function(.x, .y, .z) {
.z[.x < .y]
message("Markov chain construction...")
pseudotime <- pseudotime[waypoints]
# construct kNN graph
requireNamespace("bluster")
nbrs <- bluster::makeKNNGraph(wp_data, k = knn.)
## calculate distance of each edge
distance_m <- as.matrix(stats::dist(wp_data))
requireNamespace("igraph")
igraph::E(nbrs)$weight <- sapply(igraph::E(nbrs), function(e) {
### Get the nodes connected by each edge
nodes <- igraph::ends(nbrs, e)
node1 <- as.numeric(nodes[1])
node2 <- as.numeric(nodes[2])
### Use the distance matrix to get the distance between the two nodes
distance_m[node1, node2]
})
for (i in 1:length(waypoints)) {
if (length(KNN[i, rem_edges[[i]]]))
KNN[i, rem_edges[[i]]] <- 0
}
# determine the indice and update adjacency matrix
cell_mapping <- 1:length(waypoints)
requireNamespace("Matrix")
ids <- Matrix::summary(KNN)
# anisotropic Diffusion Kernel
aff <- exp(-(ids$x^2)/(adaptive.std[ids$i]^2) * 0.5 - (ids$x^2)/(adaptive.std[ids$j]^2) *
0.5)
W <- Matrix::sparseMatrix(i = ids$i, j = ids$j, x = aff, dims = dim(KNN), giveCsparse = TRUE)
# Transition matrix
D <- apply(W, 1, sum)
ids <- Matrix::summary(W)
T_ <- Matrix::sparseMatrix(i = ids$i, j = ids$j, x = ids$x/D[ids$i], dims = dim(KNN),
giveCsparse = TRUE)
return(T_)
}
## generate weighted adjacent matrix of this knn graph

KNN <- igraph::as_adjacency_matrix(nbrs, attr = "weight")
## generate the index of each neighbor
idx <- KNN@p
idx_seq <- mapply(seq, (idx + 1)[-length(idx)], idx[-1])
ind <- lapply(idx_seq, function(x) {
KNN@i[x] + 1
})
## select Standard deviation allowing for 'back' edges
adaptive.k <- min(c(floor(knn. / 3) - 1, 30))
dist_ <- lapply(idx_seq, function(y) {
KNN@x[y]
})
dist_sort <- lapply(dist_, sort, decreasing = TRUE)
adaptive.std <- vapply(dist_sort, "[", adaptive.k, FUN.VALUE = double(1))
# Directed graph construction pseudotime position of all the neighbors
traj_nbrs <- lapply(ind, function(x) {
pseudotime[x]
})
## Remove edges that move backwards in pseudotime except for edges that are
## within the computed standard deviation
requireNamespace("purrr")
rem_edges <- purrr::pmap(list(
.x = traj_nbrs, .y = (pseudotime - adaptive.std),
.z = ind
), function(.x, .y, .z) {
.z[.x < .y]
})
for (i in seq_len(length(waypoints))) {
if (length(KNN[i, rem_edges[[i]]])) {
KNN[i, rem_edges[[i]]] <- 0
}
}
# determine the indice and update adjacency matrix
cell_mapping <- seq_len(length(waypoints))
requireNamespace("Matrix")
ids <- Matrix::summary(KNN)
# anisotropic Diffusion Kernel
aff <- exp(-(ids$x^2) / (adaptive.std[ids$i]^2) * 0.5 - (ids$x^2) / (adaptive.std[ids$j]^2) *
0.5)
W <- Matrix::sparseMatrix(i = ids$i, j = ids$j, x = aff, dims = dim(KNN), giveCsparse = TRUE)
# Transition matrix
D <- apply(W, 1, sum)
ids <- Matrix::summary(W)
T_ <- Matrix::sparseMatrix(
i = ids$i, j = ids$j, x = ids$x / D[ids$i], dims = dim(KNN),
giveCsparse = TRUE
)
return(T_)
}
40 changes: 20 additions & 20 deletions R/determ.multiscale.space.R
Original file line number Diff line number Diff line change
@@ -1,26 +1,26 @@
#' .determine.multiscale.space
#'
#'
#' @param diffusionmap DiffusionMap object
#' @param n_eigs integer, default is NULL. Number of eigen vectors to use.
#' - If is not specified, the number of eigen vectors will be determined using the eigen gap.
#' - If is not specified, the number of eigen vectors will be determined using the eigen gap.
#' @returns dataframe
.determine.multiscale.space <- function(diffusionmap, n_eigs = NULL) {
.class.check(diffusionmap, "DiffusionMap")
requireNamespace("destiny")
eigenvect <- destiny::eigenvectors(diffusionmap)
eigenval <- destiny::eigenvalues(diffusionmap)
# determine n_eigs
if (is.null(n_eigs)) {
val_gaps <- eigenval[2:length(eigenval) - 1] - eigenval[2:length(eigenval)]
n_eigs <- order(val_gaps)[length(val_gaps)]
if (n_eigs < 3) {
n_eigs <- order(val_gaps)[length(val_gaps) - 1]
.class.check(diffusionmap, "DiffusionMap")
requireNamespace("destiny")
eigenvect <- destiny::eigenvectors(diffusionmap)
eigenval <- destiny::eigenvalues(diffusionmap)
# determine n_eigs
if (is.null(n_eigs)) {
val_gaps <- eigenval[2:length(eigenval) - 1] - eigenval[2:length(eigenval)]
n_eigs <- order(val_gaps)[length(val_gaps)]
if (n_eigs < 3) {
n_eigs <- order(val_gaps)[length(val_gaps) - 1]
}
}
}
# Scale the data
use_eigs <- seq(1, n_eigs)
eigenvect <- eigenvect[, use_eigs]
eigenval <- eigenval[use_eigs]
multis <- apply(eigenvect, 1, `*`, (eigenval/(1 - eigenval)))
return(data.frame(t(multis)))
}
# Scale the data
use_eigs <- seq(1, n_eigs)
eigenvect <- eigenvect[, use_eigs]
eigenval <- eigenval[use_eigs]
multis <- apply(eigenvect, 1, `*`, (eigenval / (1 - eigenval)))
return(data.frame(t(multis)))
}
82 changes: 44 additions & 38 deletions R/differentiation_probabilities.R
Original file line number Diff line number Diff line change
Expand Up @@ -9,42 +9,48 @@
#' @return probabilities
#' @include construct.markov.chain.R
#' @include terminal.state.from.markov.chain.R
differentiation_probabilities <- function(wp_data, terminal_states = NULL, knn. = 30L,
pseudotime, waypoints) {
T_ <- .construct.markov.chain(wp_data, 30, pseudotime, waypoints)
# identify terminal states if not specified
if (is.null(terminal_states)) {
terminal_states <- .terminal.state.from.markov.chain(T_, wp_data, pseudotime, waypoints)
}
abs_states_idx <- which(waypoints %in% terminal_states)
T_[abs_states_idx, ] <- 0
T_ <- Reduce(function(matri, x) {
matri[x, x] <- 1
matri
}, x = abs_states_idx, init = T_)
message("Computing fundamental matrix and absorption probabilities...")
# Transition states
trans_states_idx <- which(!(waypoints %in% terminal_states))
# Q matrix
Q <- T_[-abs_states_idx, -abs_states_idx]
# Fundamental matrix
mat <- diag(dim(Q)[[1]]) - Q
requireNamespace("Matrix")
requireNamespace("MASS")
N <- tryCatch({
Matrix::solve(mat)
}, error = function(cnd) {
warning("Matrix generated is singular or nearly singular; using pseudo-inverse to construct fundamental matrix")
warning("Or you can re-run this function to reconstruct the markov chain")
MASS::ginv(as.matrix(mat))
})
R <- T_[trans_states_idx, abs_states_idx]
# absorbing probabilities:
probabilities <- N %*% R
probabilities@x[probabilities@x < 0] <- 0
# add terminal states
probabilities <- rbind(probabilities, T_[abs_states_idx, abs_states_idx])
probabilities <- probabilities[order(c(trans_states_idx, abs_states_idx)), ]
return(probabilities)
differentiation_probabilities <- function(
wp_data, terminal_states = NULL, knn. = 30L,
pseudotime, waypoints) {
T_ <- .construct.markov.chain(wp_data, 30, pseudotime, waypoints)
# identify terminal states if not specified
if (is.null(terminal_states)) {
terminal_states <- .terminal.state.from.markov.chain(
T_, wp_data, pseudotime,
waypoints
)
}
abs_states_idx <- which(waypoints %in% terminal_states)
T_[abs_states_idx, ] <- 0
T_ <- Reduce(function(matri, x) {
matri[x, x] <- 1
matri
}, x = abs_states_idx, init = T_)
message("Computing fundamental matrix and absorption probabilities...")
# Transition states
trans_states_idx <- which(!(waypoints %in% terminal_states))
# Q matrix
Q <- T_[-abs_states_idx, -abs_states_idx]
# Fundamental matrix
mat <- diag(dim(Q)[[1]]) - Q
requireNamespace("Matrix")
requireNamespace("MASS")
N <- tryCatch(
{
Matrix::solve(mat)
},
error = function(cnd) {
warning("Matrix generated is singular or nearly singular; using pseudo-inverse to construct fundamental matrix")
warning("Or you can re-run this function to reconstruct the markov chain")
MASS::ginv(as.matrix(mat))
}
)
R <- T_[trans_states_idx, abs_states_idx]
# absorbing probabilities:
probabilities <- N %*% R
probabilities@x[probabilities@x < 0] <- 0
# add terminal states
probabilities <- rbind(probabilities, T_[abs_states_idx, abs_states_idx])
probabilities <- probabilities[order(c(trans_states_idx, abs_states_idx)), ]
return(probabilities)
}

34 changes: 17 additions & 17 deletions R/filter.cells.R
Original file line number Diff line number Diff line change
Expand Up @@ -4,26 +4,26 @@
#' @param sce SingleCellExperiment object, adata in python
#' data after combineTCR, contain both vdj and seq
#' @param col_n mode for extraction the V(D)J genes.
#' @param filter_pattern character string, optional ",|None|No_contig" by default
#' @param filter_pattern character string, optional ',|None|No_contig' by default
#' @param remove_missing bool, True by default
#' - If true, will remove cells with contigs matching the filter from the object.
#' - If False, will mask them with a uniform value dependent on the column name.
#' @import SingleCellExperiment
#' @return filtered SingleCellExperiment object according to the parameter.
.filter.cells <- function(sce, col_n, filter_pattern = ",|None|No_cotig", remove_missing = TRUE) {
requireNamespace("rlang")
if (ncol(sce) < 1) {
rlang::abort("None column remains, please check whether the filtering option is correct.")
}
# find filter pattern hits in our column of interest
temp <- colData(sce)[col_n]
temp <- as.data.frame(temp@listData)
mask <- grep(filter_pattern, temp[, 1])
if (remove_missing) {
sce <- sce[, setdiff(1:ncol(sce), mask)]
} else {
# uniformly mask the filter pattern hits
colData(sce)[mask, col_n] <- sprintf("%s_missing", col_n)
}
return(sce)
}
requireNamespace("rlang")
if (ncol(sce) < 1) {
rlang::abort("None column remains, please check whether the filtering option is correct.")
}
# find filter pattern hits in our column of interest
temp <- colData(sce)[col_n]
temp <- as.data.frame(temp@listData)
mask <- grep(filter_pattern, temp[, 1])
if (remove_missing) {
sce <- sce[, setdiff(seq_len(ncol(sce)), mask)]
} else {
# uniformly mask the filter pattern hits
colData(sce)[mask, col_n] <- sprintf("%s_missing", col_n)
}
return(sce)
}
Loading

0 comments on commit 718c5e8

Please sign in to comment.