#' The scAI Class
#'
#' The scAI object is created from a paired single-cell transcriptomic and epigenomic data.
#' It takes a list of two digital data matrices as input. Genes/loci should be in rows and cells in columns. rownames and colnames should be included.
#' The class provides functions for data preprocessing, integrative analysis, and visualization.
#'
#' The key slots used in the scAI object are described below.
#'
#' @slot raw.data List of raw data matrices, one per dataset (Genes/loci should be in rows and cells in columns)
#' @slot norm.data List of normalized matrices (genes/loci by cells)
#' @slot agg.data Aggregated epigenomic data within similar cells
#' @slot scale.data List of scaled matrices
#' @slot pData data frame storing the information associated with each cell
#' @slot var.features List of informative features to be used, one giving informative genes and the other giving informative loci
#' @slot fit List of inferred low-rank matrices, including W1, W2, H, Z, R
#' @slot fit.variedK List of inferred low-rank matrices when varying the rank K
#' @slot embed List of the reduced 2D coordinates, one per method, e.g., t-SNE/FIt-SNE/umap
#' @slot identity a factor defining the cell identity
#' @slot cluster List of consensus clustering results
#' @slot options List of parameters used throughout analysis
#'
#' @exportClass scAI
#' @importFrom Rcpp evalCpp
#' @importFrom methods setClass
#' @useDynLib scAI
scAI <- methods::setClass("scAI",
slots = c(raw.data = "list",
norm.data = "list",
agg.data = "matrix",
scale.data = "list",
pData = "data.frame",
var.features = "list",
fit = "list",
fit.variedK = "list",
embed = "list",
identity = "factor",
cluster = "list",
options = "list")
)
#' show method for scAI
#'
#' @param scAI object
#' @param show show the object
#' @param object object
#' @docType methods
#'
setMethod(f = "show", signature = "scAI", definition = function(object) {
cat("An object of class", class(object), "\n", length(object@raw.data), "datasets.\n")
invisible(x = NULL)
})
#' creat a new scAI object
#'
#' @param raw.data List of raw data matrices, a paired single-cell transcriptomic and epigenomic data
#' @param do.sparse whether use sparse format
#'
#' @return
#' @export
#'
#' @examples
#' @importFrom methods as new
create_scAIobject <- function(raw.data, do.sparse = T) {
object <- methods::new(Class = "scAI",
raw.data = raw.data)
if (do.sparse) {
raw.data <- lapply(raw.data, function(x) {
as(as.matrix(x), "dgCMatrix")
})
}
object@raw.data <- raw.data
return(object)
}
#' preprocess the raw.data including quality control and normalization
#'
#' @param object scAI object
#' @param assay List of assay names to be normalized
#' @param minFeatures Filter out cells with expressed features < minimum features
#' @param minCells Filter out features expressing in less than minCells
#' @param minCounts Filter out cells with expressed count < minCounts
#' @param maxCounts Filter out cells with expressed count > minCounts
#' @param libararyflag Whether do library size normalization
#' @param logNormalize whether do log transformation
#'
#' @return
#' @export
#'
#' @examples
preprocessing <- function(object, assay = list("RNA", "ATAC"), minFeatures = 200, minCells = 3, minCounts = NULL, maxCounts = NULL, libararyflag = TRUE, logNormalize = TRUE) {
if (is.null(assay)) {
for (i in 1:length(object@raw.data)) {
object@norm.data[[i]] <- object@raw.data[[i]]
}
} else {
for (i in 1:length(assay)) {
iniData <- object@raw.data[[assay[[i]]]]
print(dim(iniData))
if (class(iniData) == "data.frame") {
iniData <- as.matrix(iniData)
}
# filter cells that have features less than #minFeatures
msum <- Matrix::colSums(iniData != 0)
proData <- iniData[, msum >= minFeatures]
# filter cells that have UMI counts less than #minCounts
if (!is.null(minCounts)) {
proData <- proData[, Matrix::colSums(proData) >= minCounts]
}
# filter cells that have expressed genes high than #maxGenes
if (!is.null(maxCounts)) {
proData <- proData[, Matrix::colSums(proData) <= maxCounts]
}
# filter genes that only express less than #minCells cells
proData <- proData[Matrix::rowSums(proData != 0) >= minCells, ]
# normalization:we employ a global-scaling normalization method that normalizes the gene expression measurements for each cell by the total expression multiplies this by a scale factor (10,000 by default)
if (libararyflag) {
proData <- sweep(proData, 2, Matrix::colSums(proData), FUN = `/`) * 10000
}
if (logNormalize) {
proData = log(proData + 1)
}
object@norm.data[[assay[[i]]]] <- proData
}
}
if (length(assay) == 2) {
X1 <- object@norm.data[[assay[[1]]]]
X2 <- object@norm.data[[assay[[2]]]]
cell.keep <- intersect(colnames(X1), colnames(X2))
object@norm.data[[assay[[1]]]] <- X1[, cell.keep]
object@norm.data[[assay[[2]]]] <- X2[, cell.keep]
} else if (length(assay) == 1) {
X1 <- object@norm.data[[assay[[1]]]]
assay2 <- setdiff(names(object@raw.data), assay[[1]])
X2 <- object@raw.data[[assay2]]
cell.keep <- intersect(colnames(X1), colnames(X2))
object@norm.data[[assay[[1]]]] <- X1[, cell.keep]
object@norm.data[[assay2]] <- X2[, cell.keep]
}
names(object@norm.data) <- names(object@raw.data)
return(object)
}
#' add the cell information into pData slot
#'
#' @param object scAi object
#' @param pdata cell information to be added
#' @param pdata.name the name of column to be assigned
#'
#' @return
#' @export
#'
#' @examples
addpData <- function(object, pdata, pdata.name = NULL) {
if (is.null(x = pdata.name) && is.atomic(x = pdata)) {
stop("'pdata.name' must be provided for atomic pdata types (eg. vectors)")
}
if (inherits(x = pdata, what = c("matrix", "Matrix"))) {
pdata <- as.data.frame(x = pdata)
}
if (is.null(x = pdata.name)) {
pdata.name <- names(pdata)
} else {
names(pdata) <- pdata.name
}
object@pData <- pdata
return(object)
}
#' run scAI model
#'
#' @param object scAI object
#' @param K An integer indicating the Rank of the inferred factor
#' @param do.fast whether use the python version for running scAI optimization model
#' @param nrun Number of times to repreat the running
#' @param hvg.use1 whether use high variable genes for RNA-seq data
#' @param hvg.use2 whether use high variable genes for ATAC-seq data
#' @param keep_all Whether keep all the results from multiple runs
#' @param s Probability of Bernoulli distribution
#' @param alpha model parameter
#' @param lambda model parameter
#' @param gamma model parameter
#' @param maxIter An integer indicating Maximum number of iteration
#' @param stop_rule Stop rule to be used
#' @param init List of the initialized low-rank matrices
#' @param rand.seed An integer indicating the seed
#'
#' @return
#' @export
#'
#' @examples
#' @importFrom foreach foreach "%dopar%"
#' @importFrom parallel makeForkCluster makeCluster detectCores
#' @importFrom doParallel registerDoParallel
#' @importFrom reticulate r_to_py source_python
run_scAI <- function(object, K, do.fast = FALSE, nrun = 5L, hvg.use1 = FALSE, hvg.use2 = FALSE, keep_all = F, s = 0.25, alpha = 1, lambda = 100000, gamma = 1, maxIter = 500L, stop_rule = 1L, init = NULL, rand.seed = 1L) {
if (!is.null(init)) {
W1.init = init$W1.init
W2.init = init$W2.init
H.init = init$H.init
Z.init = init$Z.init
R.init = init$R.init
} else {
R.init = NULL
W1.init = NULL
W2.init = NULL
H.init = NULL
Z.init = NULL
}
options(warn = -1)
# Calculate the number of cores
numCores <- min(parallel::detectCores(), nrun)
cl <- tryCatch({
parallel::makeForkCluster(numCores)
}, error = function(e) {
parallel::makeCluster(numCores)
})
doParallel::registerDoParallel(cl)
if (hvg.use1) {
X1 <- as.matrix(object@norm.data[[1]][object@var.features[[1]], ])
} else {
X1 <- as.matrix(object@norm.data[[1]])
}
if (hvg.use2) {
X2 <- as.matrix(object@norm.data[[2]][object@var.features[[2]], ])
} else {
X2 <- as.matrix(object@norm.data[[2]])
}
if (!do.fast) {
outs <- foreach(i = 1:nrun, .packages = c("Matrix")) %dopar% {
set.seed(rand.seed + i - 1)
scAImodel(X1, X2, K = K, s = s, alpha = alpha, lambda = lambda, gamma = gamma, maxIter = maxIter, stop_rule = stop_rule,
R.init = R.init, W1.init = W1.init, W2.init = W2.init, H.init = H.init, Z.init = Z.init)
}
} else {
barcodes <- colnames(X2)
feature1 <- rownames(X1)
feature2 <- rownames(X2)
names_com <- paste0("factor", seq_len(K))
path <- paste(system.file(package="scAI"), "scAImodel_py.py", sep="/")
reticulate::source_python(path)
X1py = r_to_py(X1)
X2py = r_to_py(X2)
R.initpy = r_to_py(R.init); W1.initpy = r_to_py(W1.init); W2.initpy = r_to_py(W2.init); H.initpy = r_to_py(H.init); Z.initpy = r_to_py(Z.init)
K = as.integer(K); maxIter = as.integer(maxIter)
outs <- pbapply::pbsapply(
X = 1:nrun,
FUN = function(x) {
seed = as.integer(rand.seed + x - 1)
ptm = Sys.time()
results = scAImodel_py(X1 = X1py, X2 = X2py, K = K, S = s, Alpha = alpha, Lambda = lambda, Gamma = gamma, Maxiter = maxIter, Stop_rule = stop_rule,
Seed = seed, W1 = W1.initpy, W2 = W2.initpy, H = H.initpy, Z = Z.initpy,R = R.initpy)
execution.time = Sys.time() - ptm
names(results) <- c("W1", "W2", "H", "Z", "R")
attr(results$W1, "dimnames") <- list(feature1, names_com)
attr(results$W2, "dimnames") <- list(feature2, names_com)
attr(results$H, "dimnames") <- list(names_com, barcodes)
attr(results$Z, "dimnames") <- list(barcodes, barcodes)
results$options = list(s = s, alpha = alpha, lambda = lambda, gamma = gamma, maxIter = maxIter, stop_rule = stop_rule, run.time = execution.time)
return(results)
},
simplify = FALSE
)
}
objs <- foreach(i = 1:nrun, .combine = c) %dopar% {
W1 <- outs[[i]]$W1
W2 <- outs[[i]]$W2
sum(cor(as.matrix(W1))) + sum(cor(as.matrix(W2)))
}
parallel::stopCluster(cl)
N <- ncol(X1)
C <- base::matrix(0, N, N)
for (i in seq_len(nrun)) {
H <- outs[[i]]$H
H <- sweep(H, 2, colSums(H), FUN = `/`)
clusIndex <- apply(H, 2, which.max)
# compute the consensus matrix
adjMat <- clust2Mat(clusIndex)
C <- C + adjMat
}
CM <- C/nrun
if (!keep_all) {
sprintf("The best seed is %d", which.min(objs))
outs_final <- outs[[which.min(objs)]]
#object@agg.data <- outs_final$agg.data
W = list(W1 = outs_final$W1, W2 = outs_final$W2)
names(W) <- names(object@norm.data)
object@fit <- list(W = W, H = outs_final$H, Z = outs_final$Z, R = outs_final$R)
object <- getAggregatedData(object)
object@cluster$consensus <- CM
object@options$cost <- objs
object@options$paras <- outs_final$options
object@options$paras$nrun <- nrun
object@options$best.seed <- which.min(objs)
return(object)
} else {
outs_final <- list()
outs_final$best <- outs[[which.min(objs)]]
outs_final$best$consensus <- CM
outs_final$nruns <- outs
outs_final$options$cost <- objs
return(outs_final)
}
}
#' Solving the optimization problem in scAI
#'
#' @param X1 Single-cell transcriptomic data matrix (norm.data)
#' @param X2 Single-cell epigenomic data matrix (norm.data)
#' @param K Rank of inferred factors
#' @param s Probability of Bernoulli distribution
#' @param alpha model parameter
#' @param lambda model parameter
#' @param gamma model parameter
#' @param maxIter Maximum number of iteration
#' @param stop_rule Stop rule to be used
#' @param R.init initialization of R
#' @param W1.init initialization of W1
#' @param W2.init initialization of W2
#' @param H.init initialization of H
#' @param Z.init initialization of Z
#'
#' @return
#' @export
#'
#' @examples
#' @importFrom stats rbinom runif
#' @importFrom rfunctions crossprodcpp
scAImodel <- function(X1, X2, K, s = 0.25, alpha = 1, lambda = 100000, gamma = 1, maxIter = 500, stop_rule = 1,
R.init = NULL, W1.init = NULL, W2.init = NULL, H.init = NULL, Z.init = NULL) {
# Initialization W1,W2,H and Z
p <- nrow(X1)
n <- ncol(X1)
q = nrow(X2)
if (is.null(W1.init)) {
W1.init = matrix(runif(p * K), p, K)
}
if (is.null(W2.init)) {
W2.init = matrix(runif(q * K), q, K)
}
if (is.null(H.init)) {
H.init = matrix(runif(K * n), K, n)
}
if (is.null(Z.init)) {
Z.init = matrix(runif(n), n, n)
}
if (is.null(R.init)) {
R.init = matrix(rbinom(n * n, 1, s), n, n)
}
W1 = W1.init
W2 = W2.init
H = H.init
Z = Z.init
R = R.init
# start the clock to measure the execution time
ptm = proc.time()
eps <- .Machine$double.eps
onesM_K <- matrix(1, K, K)
XtX2 <- crossprodcpp(X2)
index <- which(R == 0)
for (iter in 1:maxIter) {
# normalize H
H = H/rowSums(H)
# update W1
HHt <- tcrossprod(H)
X1Ht <- eigenMapMattcrossprod(X1, H)
W1HHt <- eigenMapMatMult(W1, HHt)
W1 <- W1 * X1Ht/(W1HHt + eps)
# update W2
ZR <- Z
ZR[index] <- 0
ZRHt <- eigenMapMattcrossprod(ZR, H)
X2ZRHt <- eigenMapMatMult(X2, ZRHt)
W2HHt <- eigenMapMatMult(W2, HHt)
W2 = W2 * X2ZRHt/(W2HHt + eps)
# update H
W1tX1 <- eigenMapMatcrossprod(W1, X1)
W2tX2 <- eigenMapMatcrossprod(W2, X2)
W2tX2ZR <- eigenMapMatMult(W2tX2, ZR)
HZZt <- eigenMapMatMult(H, Z + t(Z))
W1tW1 <- crossprodcpp(W1)
W2tW2 <- crossprodcpp(W2)
temp1 <- H * (alpha * W1tX1 + W2tX2ZR + lambda * HZZt)
temp2 <- eigenMapMatMult(alpha * W1tW1 + W2tW2 + 2 * lambda * HHt + gamma * onesM_K, H)
H <- temp1/(temp2 + eps)
# update Z
HtH <- crossprodcpp(H)
X2tW2H <- eigenMapMatcrossprod(W2tX2, H)
RX2tW2H = X2tW2H
RX2tW2H[index] = 0
XtX2ZR <- eigenMapMatMult(XtX2, ZR)
XtX2ZRR = XtX2ZR
XtX2ZRR[index] = 0
Z = Z * (RX2tW2H + lambda * HtH)/(XtX2ZRR + lambda * Z + eps)
if (stop_rule == 2) {
obj = alpha * norm(X1 - W1 %*% H, "F")^2 + norm(X2 %*% (Z * R) - W2 %*% H, "F")^2 + lambda * norm(Z - HtH, "F") + gamma * sum(colSums(H) * colSums(H))
if (iter > 1 && ((obj_old - obj)/obj_old < 10^(-6)) || iter == maxIter) {
break
}
obj_old = obj
}
}
# compute the execution time
execution.time = proc.time() - ptm
barcodes <- colnames(X2)
feature1 <- rownames(X1)
feature2 <- rownames(X2)
names_com <- paste0("factor", seq_len(K))
attr(W1, "dimnames") <- list(feature1, names_com)
attr(W2, "dimnames") <- list(feature2, names_com)
attr(H, "dimnames") <- list(names_com, barcodes)
attr(Z, "dimnames") <- list(barcodes, barcodes)
outs <- list(W1 = W1, W2 = W2, H = H, Z = Z, R = R,
options = list(s = s, alpha = alpha, lambda = lambda, gamma = gamma, maxIter = maxIter, stop_rule = stop_rule, run.time = execution.time))
return(outs)
}
#' Generate the aggregated epigenomic data
#'
#' @param object an scAI object after running run_scAI
#' @param group cell group information if available; aggregate epigenomic data based on the available cell group information instead of the learned cell-cell similarity matrix from scAI
#'
#' @return
#' @export
#'
getAggregatedData <- function(object, group = NULL) {
if (is.null(group)) {
Z <- object@fit$Z
} else {
if (is.factor(group) | is.character(group)) {
group <- as.numeric(factor(group))
}
Z <- clust2Mat(group)
diag(Z) <- 1
object@fit$Z <- Z
}
R <- object@fit$R
ZR <- Z * R
ZR <- sweep(ZR, 2, colSums(ZR), FUN = `/`)
X2 <- as.matrix(object@norm.data[[2]])
X2agg <- eigenMapMatMult(X2, ZR)
X2agg <- sweep(X2agg, 2, colSums(X2agg), FUN = `/`) * 10000
X2agg <- log(1 + X2agg)
attr(X2agg, "dimnames") <- list(rownames(object@norm.data[[2]]), colnames(object@norm.data[[2]]))
object@agg.data <- X2agg
return(object)
}
#' Perform dimensional reduction
#'
#' Dimension reduction by PCA, t-SNE or UMAP
#' @param object scAI object
#' @param return.object whether return scAI object
#' @param data.use input data
#' @param do.scale whether scale the data
#' @param do.center whether scale and center the data
#' @param method Method of dimensional reduction, one of tsne, FItsne and umap
#' @param rand.seed Set a random seed. By default, sets the seed to 42.
#' @param perplexity perplexity parameter in tsne
#' @param theta parameter in tsne
#' @param check_duplicates parameter in tsne
#' @param FItsne.path File path of FIt-SNE
#' @param dim.embed dimensions of t-SNE embedding
#' @param dim.use num of PCs used for t-SNE
#'
#' @param do.fast whether do fast PCA
#' @param dimPC the number of components to keep in PCA
#' @param weight.by.var whether use weighted pc.scores
#'
#' @param n.neighbors This determines the number of neighboring points used in
#' local approximations of manifold structure. Larger values will result in more
#' global structure being preserved at the loss of detailed local structure. In general this parameter should often be in the range 5 to 50.
#' @param n.components The dimension of the space to embed into.
#' @param distance This determines the choice of metric used to measure distance in the input space.
#' @param n.epochs the number of training epochs to be used in optimizing the low dimensional embedding. Larger values result in more accurate embeddings. If NULL is specified, a value will be selected based on the size of the input dataset (200 for large datasets, 500 for small).
#' @param learning.rate The initial learning rate for the embedding optimization.
#' @param min.dist This controls how tightly the embedding is allowed compress points together.
#' Larger values ensure embedded points are moreevenly distributed, while smaller values allow the
#' algorithm to optimise more accurately with regard to local structure. Sensible values are in the range 0.001 to 0.5.
#' @param spread he effective scale of embedded points. In combination with min.dist this determines how clustered/clumped the embedded points are.
#' @param set.op.mix.ratio Interpolate between (fuzzy) union and intersection as the set operation used to combine local fuzzy simplicial sets to obtain a global fuzzy simplicial sets.
#' @param local.connectivity The local connectivity required - i.e. the number of nearest neighbors
#' that should be assumed to be connected at a local level. The higher this value the more connected
#' the manifold becomes locally. In practice this should be not more than the local intrinsic dimension of the manifold.
#' @param repulsion.strength Weighting applied to negative samples in low dimensional embedding
#' optimization. Values higher than one will result in greater weight being given to negative samples.
#' @param negative.sample.rate The number of negative samples to select per positive sample in the
#' optimization process. Increasing this value will result in greater repulsive force being applied, greater optimization cost, but slightly more accuracy.
#' @param a More specific parameters controlling the embedding. If NULL, these values are set automatically as determined by min. dist and spread.
#' @param b More specific parameters controlling the embedding. If NULL, these values are set automatically as determined by min. dist and spread.
#'
#' @return
#' @export
#'
#' @examples
#' @importFrom Rtsne Rtsne
#' @importFrom reticulate py_module_available py_set_seed import
#'
reducedDims <- function(object, data.use = object@fit$H, do.scale = TRUE, do.center = TRUE, return.object = TRUE, method = "umap",
dim.embed = 2, dim.use = NULL, perplexity = 30, theta = 0.5, check_duplicates = F, rand.seed = 42L,
FItsne.path = NULL,
dimPC = 40,do.fast = TRUE, weight.by.var = TRUE,
n.neighbors = 30L, n.components = 2L, distance = "correlation",n.epochs = NULL,learning.rate = 1.0,min.dist = 0.3,spread = 1.0,set.op.mix.ratio = 1.0,local.connectivity = 1L,
repulsion.strength = 1,negative.sample.rate = 5,a = NULL,b = NULL) {
data.use <- as.matrix(data.use)
if (do.scale) {
data.use <- t(scale(t(data.use), center = do.center, scale = TRUE))
data.use[is.na(data.use)] <- 0
}
if (!is.null(dim.use)) {
data.use = object@embed$pca[, dim.use]
}
if (method == "pca") {
cell_coords <- runPCA(data.use, do.fast = do.fast, dimPC = dimPC, seed.use = rand.seed, weight.by.var = weight.by.var)
object@embed$pca <- cell_coords
} else if (method == "tsne") {
set.seed(rand.seed)
cell_coords <- Rtsne::Rtsne(t(data.use), pca = FALSE, dims = dim.embed, theta = theta, perplexity = perplexity, check_duplicates = F)$Y
rownames(cell_coords) <- rownames(t(data.use))
object@embed$tsne <- cell_coords
} else if (method == "FItsne") {
if (!exists("FItsne")) {
if (is.null(fitsne.path)) {
stop("Please pass in path to FIt-SNE directory as FItsne.path.")
}
source(paste0(fitsne.path, "/fast_tsne.R"), chdir = T)
}
cell_coords <- fftRtsne(t(data.use), pca = FALSE, dims = dim.embed, theta = theta, perplexity = perplexity, check_duplicates = F, rand_seed = rand.seed)
rownames(cell_coords) <- rownames(t(data.use))
object@embed$FItsne <- cell_coords
} else if (method == "umap") {
cell_coords <- runUMAP(data.use,
n.neighbors = n.neighbors,
n.components = n.components,
distance = distance,
n.epochs = n.epochs,
learning.rate = learning.rate,
min.dist = min.dist,
spread = spread,
set.op.mix.ratio = set.op.mix.ratio,
local.connectivity = local.connectivity,
repulsion.strength = repulsion.strength,
negative.sample.rate = negative.sample.rate,
a = a,
b = b,
seed.use = rand.seed
)
object@embed$umap <- cell_coords
} else {
stop("Error: unrecognized dimensionality reduction method.")
}
if (return.object) {
return(object)
} else {
return(cell_coords)
}
}
#' Dimension reduction using PCA
#'
#' @param data.use input data
#' @param do.fast whether do fast PCA
#' @param dimPC the number of components to keep
#' @param seed.use set a seed
#' @param weight.by.var whether use weighted pc.scores
#' @importFrom stats prcomp
#' @importFrom irlba irlba
#' @return
#' @export
#'
#' @examples
runPCA <- function(data.use, do.fast = T, dimPC = 50, seed.use = 42, weight.by.var = T) {
set.seed(seed = seed.use)
if (do.fast) {
dimPC <- min(dimPC, nrow(data.use) - 1)
pca.res <- irlba::irlba(t(data.use), nv = dimPC)
sdev <- pca.res$d/sqrt(max(1, ncol(data.use) - 1))
if (weight.by.var){
pc.scores <- pca.res$u %*% diag(pca.res$d)
} else {
pc.scores <- pca.res$u
}
} else {
dimPC <- min(dimPC, nrow(data.use) - 1)
pca.res <- stats::prcomp(x = t(data.use), rank. = dimPC)
sdev <- pca.res$sdev
if (weight.by.var) {
pc.scores <- pca.res$x %*% diag(pca.res$sdev[1:dimPC]^2)
} else {
pc.scores <- pca.res$x
}
}
rownames(pc.scores) <- colnames(data.use)
colnames(pc.scores) <- paste0('PC', 1:ncol(pc.scores))
cell_coords <- pc.scores
return(cell_coords)
}
#' Perform dimension reduction using UMAP
#'
#' @param data.use input data
#' @param n.neighbors This determines the number of neighboring points used in
#' local approximations of manifold structure. Larger values will result in more
#' global structure being preserved at the loss of detailed local structure. In general this parameter should often be in the range 5 to 50.
#' @param n.components The dimension of the space to embed into.
#' @param distance This determines the choice of metric used to measure distance in the input space.
#' @param n.epochs the number of training epochs to be used in optimizing the low dimensional embedding. Larger values result in more accurate embeddings. If NULL is specified, a value will be selected based on the size of the input dataset (200 for large datasets, 500 for small).
#' @param learning.rate The initial learning rate for the embedding optimization.
#' @param min.dist This controls how tightly the embedding is allowed compress points together.
#' Larger values ensure embedded points are moreevenly distributed, while smaller values allow the
#' algorithm to optimise more accurately with regard to local structure. Sensible values are in the range 0.001 to 0.5.
#' @param spread he effective scale of embedded points. In combination with min.dist this determines how clustered/clumped the embedded points are.
#' @param set.op.mix.ratio Interpolate between (fuzzy) union and intersection as the set operation used to combine local fuzzy simplicial sets to obtain a global fuzzy simplicial sets.
#' @param local.connectivity The local connectivity required - i.e. the number of nearest neighbors
#' that should be assumed to be connected at a local level. The higher this value the more connected
#' the manifold becomes locally. In practice this should be not more than the local intrinsic dimension of the manifold.
#' @param repulsion.strength Weighting applied to negative samples in low dimensional embedding
#' optimization. Values higher than one will result in greater weight being given to negative samples.
#' @param negative.sample.rate The number of negative samples to select per positive sample in the
#' optimization process. Increasing this value will result in greater repulsive force being applied, greater optimization cost, but slightly more accuracy.
#' @param a More specific parameters controlling the embedding. If NULL, these values are set automatically as determined by min. dist and spread.
#' @param b More specific parameters controlling the embedding. If NULL, these values are set automatically as determined by min. dist and spread.
#' @param seed.use Set a random seed. By default, sets the seed to 42.
#' @param metric.kwds A dictionary of arguments to pass on to the metric, such as the p value for Minkowski distance
#' @param angular.rp.forest Whether to use an angular random projection forest to initialise the
#' approximate nearest neighbor search. This can be faster, but is mostly on useful for metric that
#' use an angular style distance such as cosine, correlation etc. In the case of those metrics angular forests will be chosen automatically.
#' @param verbose Controls verbosity
#' This function is modified from Seurat package
#' @return
#' @export
#' @importFrom reticulate py_module_available py_set_seed import
#' @examples
runUMAP <- function(
data.use,
n.neighbors = 30L,
n.components = 2L,
distance = "correlation",
n.epochs = NULL,
learning.rate = 1.0,
min.dist = 0.3,
spread = 1.0,
set.op.mix.ratio = 1.0,
local.connectivity = 1L,
repulsion.strength = 1,
negative.sample.rate = 5,
a = NULL,
b = NULL,
seed.use = 42L,
metric.kwds = NULL,
angular.rp.forest = FALSE,
verbose = TRUE){
if (!reticulate::py_module_available(module = 'umap')) {
stop("Cannot find UMAP, please install through pip (e.g. pip install umap-learn or reticulate::py_install(packages = 'umap-learn')).")
}
set.seed(seed.use)
reticulate::py_set_seed(seed.use)
umap_import <- reticulate::import(module = "umap", delay_load = TRUE)
umap <- umap_import$UMAP(
n_neighbors = as.integer(n.neighbors),
n_components = as.integer(n.components),
metric = distance,
n_epochs = n.epochs,
learning_rate = learning.rate,
min_dist = min.dist,
spread = spread,
set_op_mix_ratio = set.op.mix.ratio,
local_connectivity = local.connectivity,
repulsion_strength = repulsion.strength,
negative_sample_rate = negative.sample.rate,
a = a,
b = b,
metric_kwds = metric.kwds,
angular_rp_forest = angular.rp.forest,
verbose = verbose
)
Rumap <- umap$fit_transform
umap_output <- Rumap(t(data.use))
colnames(umap_output) <- paste0('UMAP', 1:ncol(umap_output))
rownames(umap_output) <- colnames(data.use)
return(umap_output)
}
#' Identify enriched features in each factor
#'
#' Rank features in each factor by Factor loading analysis
#' @param object scAI object
#' @param assay Name of assay to be analyzed
#' @param features a vector of features
#' @param cutoff.W Threshold of feature loading values
#' @param cutoff.H Threshold of cell loading values
#' @param thresh.pc Threshold of the percent of cells enriched in one factor
#' @param thresh.fc Threshold of Fold Change
#' @param thresh.p Threshold of p-values
#' @param n.top Number of top features to be returned
#' @importFrom stats sd wilcox.test
#' @importFrom dplyr top_n slice
#' @importFrom future nbrOfWorkers
#' @importFrom future.apply future_sapply
#' @importFrom pbapply pbsapply
#' @importFrom stats p.adjust
#'
#' @return
#' @export
#'
#' @examples
identifyFactorMarkers <- function(object, assay, features = NULL, cutoff.W = 0.5, cutoff.H = 0.5,
thresh.pc = 0.05, thresh.fc = 0.25, thresh.p = 0.05, n.top = 10) {
if (assay == "RNA") {
X <- as.matrix(object@norm.data[[assay]])
} else {
X <- as.matrix(object@agg.data)
}
H <- object@fit$H
W <- object@fit$W[[assay]]
if (is.null(features)) {
features <- row.names(W)
}
H <- sweep(H, 2, colSums(H), FUN = `/`)
K = ncol(W)
lib_W <- base::rowSums(W)
lib_W[lib_W == 0] <- 1
lib_W[lib_W < mean(lib_W) - 5 * sd(lib_W)] <- 1 # omit the nearly null rows
W <- sweep(W, 1, lib_W, FUN = `/`)
MW <- base::colMeans(W)
sW <- apply(W, 2, sd)
# candidate markers for each factor
IndexW_record <- vector("list", K)
for (j in 1:K) {
IndexW_record[[j]] <- which(W[, j] > MW[j] + cutoff.W * sW[j])
}
my.sapply <- ifelse(
test = future::nbrOfWorkers() == 1,
yes = pbapply::pbsapply,
no = future.apply::future_sapply
)
# divided cells into two groups
mH <- apply(H, 1, mean)
sH <- apply(H, 1, sd)
IndexH_record1 <- vector("list", K)
IndexH_record2 <- vector("list", K)
for (i in 1:K) {
IndexH_record1[[i]] <- which(H[i, ] > mH[i] + cutoff.H * sH[i])
IndexH_record2[[i]] <- base::setdiff(c(1:ncol(H)), IndexH_record1[[i]])
}
# identify factor-specific markers
factor_markers = vector("list", K)
markersTopn = vector("list", K)
factors <- c()
Features <- c()
Pvalues <- c()
Log2FC <- c()
for (i in 1:K) {
data1 <- X[IndexW_record[[i]], IndexH_record1[[i]]]
data2 <- X[IndexW_record[[i]], IndexH_record2[[i]]]
idx1 <- which(base::rowSums(data1 != 0) > thresh.pc * ncol(data1)) # at least expressed in thresh.pc% cells in one group
FC <- log2(base::rowMeans(data1)/base::rowMeans(data2))
idx2 <- which(FC > thresh.fc)
pvalues <- my.sapply(
X = 1:nrow(x = data1),
FUN = function(x) {
return(wilcox.test(data1[x, ], data2[x, ], alternative = "greater")$p.value)
}
)
idx3 <- which(pvalues < thresh.p)
idx = intersect(intersect(idx1, idx2), idx3)
# order
FC <- FC[idx]
c = sort(FC, decreasing = T, index.return = T)$ix
ri = IndexW_record[[i]]
Pvalues <- c(Pvalues, pvalues[ri[idx[c]]])
Log2FC <- c(Log2FC, FC[c])
factors <- c(factors, rep(i, length(c)))
Features <- c(Features, features[ri[idx[c]]])
}
# stopCluster(cl)
markers.all <- cbind(factors = factors, features = Features, pvalues = Pvalues, log2FC = Log2FC)
rownames(markers.all) <- Features
markers.all <- as.data.frame(markers.all)
markers.top <- markers.all %>% dplyr::group_by(factors) %>% top_n(n.top, log2FC) %>% dplyr::slice(1:n.top)
markers = list(markers.all = markers.all, markers.top = markers.top)
return(markers)
}
#' compute the 2D coordinates of embeded cells, genes, loci and factors using VscAI visualization
#'
#' @param object scAI object
#' @param genes.embed A vector of genes to be embedded
#' @param loci.embed A vector of loci to be embedded
#' @param alpha.embed Parameter controlling the distance between the cells and the factors
#' @param snn.embed Number of neighbors
#'
#' @return
#' @export
#'
#' @examples
#' @importFrom Matrix Matrix
getEmbeddings <- function(object, genes.embed = NULL, loci.embed = NULL, alpha.embed = 1.9, snn.embed = 5) {
H <- object@fit$H
W1 <- object@fit$W[[1]]
W2 <- object@fit$W[[2]]
Z <- object@fit$Z
if (nrow(H) < 3 & length(object@fit.variedK) == 0) {
print("VscAI needs at least three factors for embedding. Now rerun scAI with rank being 3...")
outs <- run_scAI(object, K = 3, nrun = object@options$paras$nrun,
s = object@options$paras$s, alpha = object@options$paras$alpha, lambda = object@options$paras$lambda, gamma = object@options$paras$gamma,
maxIter = object@options$paras$maxIter, stop_rule = object@options$paras$stop_rule)
W <- outs@fit$W
H <- outs@fit$H
Z <- outs@fit$Z
W1 <- W[[1]]
W2 <- W[[2]]
object@fit.variedK$W <- W
object@fit.variedK$H <- H
object@fit.variedK$Z <- Z
} else if (nrow(H) < 3 & length(object@fit.variedK) != 0) {
print("Using the previous calculated factors for embedding.")
W1 <- object@fit.variedK$W[[1]]
W2 <- object@fit.variedK$W[[2]]
H <- object@fit.variedK$H
Z <- object@fit.variedK$Z
}
nmf.scores <- H
Zsym <- Z/max(Z)
Zsym[Zsym < 10^(-6)] <- 0
Zsym <- (Zsym + t(Zsym))/2
diag(Zsym) <- 1
rownames(Zsym) <- colnames(H)
colnames(Zsym) <- rownames(Zsym)
alpha.exp <- alpha.embed # Increase this > 1.0 to move the cells closer to the factors. Values > 2 start to distort the data.
snn.exp <- snn.embed # Lower this < 1.0 to move similar cells closer to each other
n_pull <- nrow(H) # The number of factors pulling on each cell. Must be at least 3.
Zsym <- Matrix(data = Zsym, sparse = TRUE)
snn <- Zsym
swne.embedding <- swne::EmbedSWNE(nmf.scores, snn, alpha.exp = alpha.exp, snn.exp = snn.exp, n_pull = n_pull, proj.method = "sammon", dist.use = "cosine")
# project cells and factors
factor_coords <- swne.embedding$H.coords
cell_coords <- swne.embedding$sample.coords
rownames(factor_coords) <- rownames(H)
rownames(cell_coords) <- colnames(H)
# project genes onto the low dimension space by VscAI
if (is.null(genes.embed)) {
genes.embed <- rownames(W1)
} else {
genes.embed <- as.character(as.vector(as.matrix((genes.embed))))
}
swne.embedding <- swne::EmbedFeatures(swne.embedding, W1, genes.embed, n_pull = n_pull)
gene_coords <- swne.embedding$feature.coords
rownames(gene_coords) <- gene_coords$name
# project loci
if (is.null(loci.embed)) {
loci.embed <- rownames(W2)
} else {
loci.embed <- as.character(as.vector(as.matrix((loci.embed))))
}
swne.embedding <- swne::EmbedFeatures(swne.embedding, W2, loci.embed, n_pull = n_pull)
loci_coords <- swne.embedding$feature.coords
rownames(loci_coords) <- loci_coords$name
object@embed$VscAI$cells <- cell_coords
object@embed$VscAI$factors <- factor_coords
object@embed$VscAI$genes <- gene_coords
object@embed$VscAI$loci <- loci_coords
return(object)
}
#' Identify cell clusters
#'
#' @param object scAI object
#' @param partition.type Type of partition to use. Defaults to RBConfigurationVertexPartition. Options include: ModularityVertexPartition, RBERVertexPartition, CPMVertexPartition, MutableVertexPartition, SignificanceVertexPartition, SurpriseVertexPartition (see the Leiden python module documentation for more details)
#' @param seed.use set seed
#' @param n.iter number of iteration
#' @param initial.membership arameters to pass to the Python leidenalg function defaults initial_membership=None
#' @param weights defaults weights=None
#' @param node.sizes Parameters to pass to the Python leidenalg function
#' @param resolution A parameter controlling the coarseness of the clusters
#' @param K Number of clusters if performing hierarchical clustering of the consensus matrix
#' @return
#' @export
#'
#' @examples
#' @importFrom stats as.dist cophenetic cor cutree dist hclust
identifyClusters <- function(object, resolution = 1, partition.type = "RBConfigurationVertexPartition",
seed.use = 42L,n.iter = 10L,
initial.membership = NULL, weights = NULL, node.sizes = NULL,
K = NULL) {
if (is.null(K) & !is.null(resolution)) {
data.use <- object@fit$H
data.use <- as.matrix(data.use)
data.use <- t(scale(t(data.use), center = TRUE, scale = TRUE))
snn <- swne::CalcSNN(data.use, k = 20, prune.SNN = 1/15)
idents <- runLeiden(SNN = snn, resolution = resolution, partition_type = partition.type,
seed.use = seed.use, n.iter = n.iter,
initial.membership = initial.membership, weights = weights, node.sizes = node.sizes)
} else {
CM <- object@cluster$consensus
d <- as.dist(1 - CM)
hc <- hclust(d, "ave")
idents<- hc %>% cutree(k = K)
coph <- cophenetic(hc)
cs <- cor(d, coph)
object@cluster$cluster.stability <- cs
}
names(idents) <- rownames(object@pData)
object@identity <- factor(idents)
object@pData$cluster <- factor(idents)
return(object)
}
#' Run Leiden clustering algorithm
#' This code is modified from Tom Kelly (https://github.com/TomKellyGenetics/leiden), where we added more parameters (seed.use and n.iter) to run the Python version. In addition, we also take care of the singleton issue after running leiden algorithm.
#' @description Implements the Leiden clustering algorithm in R using reticulate to run the Python version. Requires the python "leidenalg" and "igraph" modules to be installed. Returns a vector of partition indices.
#' @param SNN An adjacency matrix compatible with \code{\link[igraph]{igraph}} object.
#' @param seed.use set seed
#' @param n.iter number of iteration
#' @param initial.membership arameters to pass to the Python leidenalg function defaults initial_membership=None
#' @param node.sizes Parameters to pass to the Python leidenalg function
#' @param partition_type Type of partition to use. Defaults to RBConfigurationVertexPartition. Options include: ModularityVertexPartition, RBERVertexPartition, CPMVertexPartition, MutableVertexPartition, SignificanceVertexPartition, SurpriseVertexPartition (see the Leiden python module documentation for more details)
#' @param resolution A parameter controlling the coarseness of the clusters
#' @param weights defaults weights=None
#' @return A parition of clusters as a vector of integers
##'
##'
#' @keywords graph network igraph mvtnorm simulation
#' @importFrom reticulate import r_to_py
##' @export
runLeiden <- function(SNN = matrix(), resolution = 1, partition_type = c(
'RBConfigurationVertexPartition',
'ModularityVertexPartition',
'RBERVertexPartition',
'CPMVertexPartition',
'MutableVertexPartition',
'SignificanceVertexPartition',
'SurpriseVertexPartition'),
seed.use = 42L,
n.iter = 10L,
initial.membership = NULL, weights = NULL, node.sizes = NULL) {
if (!reticulate::py_module_available(module = 'leidenalg')) {
stop("Cannot find Leiden algorithm, please install through pip (e.g. pip install leidenalg).")
}
#import python modules with reticulate
leidenalg <- import("leidenalg", delay_load = TRUE)
ig <- import("igraph", delay_load = TRUE)
resolution_parameter <- resolution
initial_membership <- initial.membership
node_sizes <- node.sizes
#convert matrix input (corrects for sparse matrix input)
adj_mat <- as.matrix(SNN)
#compute weights if non-binary adjacency matrix given
is_pure_adj <- all(as.logical(adj_mat) == adj_mat)
if (is.null(weights) && !is_pure_adj) {
#assign weights to edges (without dependancy on igraph)
weights <- t(adj_mat)[t(adj_mat)!=0]
#remove zeroes from rows of matrix and return vector of length edges
}
##convert to python numpy.ndarray, then a list
adj_mat_py <- r_to_py(adj_mat)
adj_mat_py <- adj_mat_py$tolist()
#convert graph structure to a Python compatible object
GraphClass <- if (!is.null(weights) && !is_pure_adj){
ig$Graph$Weighted_Adjacency
} else {
ig$Graph$Adjacency
}
snn_graph <- GraphClass(adj_mat_py)
#compute partitions
partition_type <- match.arg(partition_type)
part <- switch(
EXPR = partition_type,
'RBConfigurationVertexPartition' = leidenalg$find_partition(
snn_graph,
leidenalg$RBConfigurationVertexPartition,
initial_membership = initial.membership, weights = weights,
resolution_parameter = resolution,
n_iterations = n.iter,
seed = seed.use
),
'ModularityVertexPartition' = leidenalg$find_partition(
snn_graph,
leidenalg$ModularityVertexPartition,
initial_membership = initial.membership, weights = weights,
n_iterations = n.iter,
seed = seed.use
),
'RBERVertexPartition' = leidenalg$find_partition(
snn_graph,
leidenalg$RBERVertexPartition,
initial_membership = initial.membership, weights = weights, node_sizes = node.sizes,
resolution_parameter = resolution_parameter,
n_iterations = n.iter,
seed = seed.use
),
'CPMVertexPartition' = leidenalg$find_partition(
snn_graph,
leidenalg$CPMVertexPartition,
initial_membership = initial.membership, weights = weights, node_sizes = node.sizes,
resolution_parameter = resolution,
n_iterations = n.iter,
seed = seed.use
),
'MutableVertexPartition' = leidenalg$find_partition(
snn_graph,
leidenalg$MutableVertexPartition,
initial_membership = initial.membership,
n_iterations = n.iter,
seed = seed.use
),
'SignificanceVertexPartition' = leidenalg$find_partition(
snn_graph,
leidenalg$SignificanceVertexPartition,
initial_membership = initial.membership, node_sizes = node.sizes,
resolution_parameter = resolution,
n_iterations = n.iter,
seed = seed.use
),
'SurpriseVertexPartition' = leidenalg$find_partition(
snn_graph,
leidenalg$SurpriseVertexPartition,
initial_membership = initial.membership, weights = weights, node_sizes = node.sizes,
n_iterations = n.iter,
seed = seed.use
),
stop("please specify a partition type as a string out of those documented")
)
partition <- part$membership+1
idents <- partition
if (min(table(idents)) == 1) {
idents <- assignSingletons(idents, SNN)
}
idents <- factor(idents)
names(idents) <- row.names(SNN)
return(idents)
}
# Group single cells that make up their own cluster in with the cluster they are most connected to.
#
# @param idents clustering result
# @param SNN SNN graph
# @return Returns scAI object with all singletons merged with most connected cluster
assignSingletons <- function(idents, SNN) {
# identify singletons
singletons <- c()
for (cluster in unique(idents)) {
if (length(which(idents %in% cluster)) == 1) {
singletons <- append(x = singletons, values = cluster)
}
}
#singletons = names(table(idents))[which(table(idents)==1)]
# calculate connectivity of singletons to other clusters, add singleton to cluster it is most connected to
cluster_names <- unique(x = idents)
cluster_names <- setdiff(x = cluster_names, y = singletons)
connectivity <- vector(mode="numeric", length = length(x = cluster_names))
names(x = connectivity) <- cluster_names
for (i in singletons) {
print(i)
for (j in cluster_names) {
subSNN = SNN[
which(idents %in% i),
which(idents %in% j)
]
if (is.object(x = subSNN)) {
connectivity[j] <- sum(subSNN) / (nrow(x = subSNN) * ncol(x = subSNN))
} else {
connectivity[j] <- mean(x = subSNN)
}
}
m <- max(connectivity, na.rm = T)
mi <- which(x = connectivity == m, arr.ind = TRUE)
closest_cluster <- sample(x = names(x = connectivity[mi]), 1)
which(idents %in% i)[which(idents %in% i)] <- closest_cluster
}
if (length(x = singletons) > 0) {
message(paste(
length(x = singletons),
"singletons identified.",
length(x = unique(idents)),
"final clusters."
))
}
return(idents)
}
#' Convert membership into an adjacent matrix
#'
#' @param memb Membership vector
#'
#' @return
#' @export
#'
#' @examples
clust2Mat <- function(memb) {
N <- length(memb)
return(as.numeric(outer(memb, memb, FUN = "==")) - outer(1:N, 1:N, "=="))
}
#' Identify markers in each cell cluster
#'
#' @param object scAI object
#' @param assay Name of assay to be analyzed
#' @param features a vector of features
#' @param use.agg whether use the aggregated data
#' @param test.use which test to use ("bimod" or "wilcox")
#' @param thresh.pc Threshold of the percent of cells enriched in one cluster
#' @param thresh.fc Threshold of Fold Change
#' @param thresh.p Threshold of p-values
#' @importFrom future nbrOfWorkers
#' @importFrom pbapply pbsapply
#' @importFrom future.apply future_sapply
#' @importFrom stats sd wilcox.test
#' @importFrom dplyr top_n slice
#' @importFrom stats p.adjust
#'
#' @return
#' @export
#'
#' @examples
identifyClusterMarkers <- function(object, assay, features = NULL, use.agg = TRUE, test.use = "bimod",
thresh.pc = 0.25, thresh.fc = 0.25, thresh.p = 0.01) {
if (assay == "RNA") {
X <- object@norm.data[[assay]]
} else {
if (use.agg) {
X <- object@agg.data
} else {
X <- object@norm.data[[assay]]
}
}
if (is.null(features)) {
features.use <- row.names(X)
} else {
features.use <- intersect(features, row.names(X))
}
data.use <- X[features.use,]
data.use <- as.matrix(data.use)
groups <- object@identity
labels <- groups
level.use <- levels(labels)
numCluster <- length(level.use)
my.sapply <- ifelse(
test = future::nbrOfWorkers() == 1,
yes = pbapply::pbsapply,
no = future.apply::future_sapply
)
mean.fxn <- function(x) {
return(log(x = mean(x = expm1(x = x)) + 1))
}
genes.de <- vector("list", length = numCluster)
for (i in 1:numCluster) {
features <- features.use
cell.use1 <- which(labels %in% level.use[i])
cell.use2 <- base::setdiff(1:length(labels), cell.use1)
# feature selection (based on percentages)
thresh.min <- 0
pct.1 <- round(
x = rowSums(x = data.use[features, cell.use1, drop = FALSE] > thresh.min) /
length(x = cell.use1),
digits = 3
)
pct.2 <- round(
x = rowSums(x = data.use[features, cell.use2, drop = FALSE] > thresh.min) /
length(x = cell.use2),
digits = 3
)
data.alpha <- cbind(pct.1, pct.2)
colnames(x = data.alpha) <- c("pct.1", "pct.2")
alpha.min <- apply(X = data.alpha, MARGIN = 1, FUN = max)
names(x = alpha.min) <- rownames(x = data.alpha)
features <- names(x = which(x = alpha.min > thresh.pc))
if (length(x = features) == 0) {
stop("No features pass thresh.pc threshold")
}
# feature selection (based on average difference)
data.1 <- apply(X = data.use[features, cell.use1, drop = FALSE],MARGIN = 1,FUN = mean.fxn)
data.2 <- apply(X = data.use[features, cell.use2, drop = FALSE],MARGIN = 1,FUN = mean.fxn)
FC <- (data.1 - data.2)
features.diff <- names(x = which(x = FC > thresh.fc))
features <- intersect(x = features, y = features.diff)
if (length(x = features) == 0) {
stop("No features pass thresh.fc threshold")
}
data1 <- data.use[features, cell.use1, drop = FALSE]
data2 <- data.use[features, cell.use2, drop = FALSE]
if (test.use == "wilcox") {
pvalues <- unlist(
x = my.sapply(
X = 1:nrow(x = data1),
FUN = function(x) {
return(wilcox.test(data1[x, ], data2[x, ], alternative = "greater")$p.value)
}
)
)
} else if (test.use == 'bimod') {
pvalues <- unlist(
x = my.sapply(
X = 1:nrow(x = data1),
FUN = function(x) {
return(DifferentialLRT(
x = as.numeric(x = data1[x,]),
y = as.numeric(x = data2[x,])
))
}
)
)
}
pval.adj = stats::p.adjust(
p = pvalues,
method = "bonferroni",
n = nrow(X)
)
genes.de[[i]] <- data.frame(clusters = level.use[i], features = as.character(rownames(data1)), pvalues = pvalues, pval_adj = pval.adj, logFC = FC[features], data.alpha[features,])
}
markers.all <- data.frame()
for (i in 1:numCluster) {
gde <- genes.de[[i]]
gde <- gde[order(gde$pvalues, -gde$logFC), ]
gde <- subset(gde, subset = pvalues < thresh.p)
if (nrow(gde) > 0) {
markers.all <- rbind(markers.all, gde)
}
}
markers.all$features <- as.character(markers.all$features)
return(markers.all)
}
# function to run mcdavid et al. DE test
#' likelood ratio test
#'
#' @param x a vector
#' @param y a vector
#' @param xmin threshold for the values in the vector
#'
#'
#' @importFrom stats pchisq
DifferentialLRT <- function(x, y, xmin = 0) {
lrtX <- bimodLikData(x = x)
lrtY <- bimodLikData(x = y)
lrtZ <- bimodLikData(x = c(x, y))
lrt_diff <- 2 * (lrtX + lrtY - lrtZ)
return(pchisq(q = lrt_diff, df = 3, lower.tail = F))
}
#' likelood ratio test
#' @importFrom stats sd dnorm
#' @param x a vector
#' @param xmin threshold for the values in the vector
bimodLikData <- function(x, xmin = 0) {
x1 <- x[x <= xmin]
x2 <- x[x > xmin]
xal <- length(x = x2) / length(x = x)
xal[xal > 1 - 1e-5] <- 1 - 1e-5
xal[xal < 1e-5] <- 1e-5
likA <- length(x = x1) * log(x = 1 - xal)
if (length(x = x2) < 2) {
mysd <- 1
} else {
mysd <- sd(x = x2)
}
likB <- length(x = x2) *
log(x = xal) +
sum(dnorm(x = x2, mean = mean(x = x2), sd = mysd, log = TRUE))
return(likA + likB)
}
#' Infer regulatory relationships
#'
#' @param object scAI object
#' @param gene.use genes to be inferred
#' @param candinate_loci cadinate loci
#' @param cutoff_H threshold of coefficient values
#' @param cutoff the difference of correlation to be considered as significant
#' @param thresh_corr threshold of correlation coefficient
#'
#' @return
#' @export
#'
#' @examples
inferRegulations <- function(object, gene.use, candinate_loci, cutoff_H = 0.5, cutoff = 0.1, thresh_corr = 0.2) {
H <- as.matrix(object@fit$H)
H <- sweep(H, 2, colSums(H), FUN = `/`)
K = nrow(H)
mH <- apply(H, 1, mean)
sH <- apply(H, 1, sd)
IndexH_record <- vector("list", K)
for (i in 1:K) {
IndexH_record[[i]] <- which(H[i, ] > mH[i] + cutoff_H * sH[i])
}
gene.use <- as.character(gene.use)
X1 <- object@norm.data[[1]]
X1 <- X1[gene.use, ]
X2a <- object@agg.data
regulations <- vector("list", length(gene.use))
names(regulations) <- gene.use
for (i in 1:length(gene.use)) {
regulations_i = vector("list", K)
names(regulations_i) <- rownames(H)
for (j in 1:K) {
loci_j <- candinate_loci$markers.all$features[candinate_loci$markers.all$factors == j]
# compute the correlation between
x1 <- X1[i, ]
x2a <- X2a[loci_j, ]
cors1 <- cor(x1, t(x2a))
# set the values of this gene and its candiate loci to be zero
X1_new <- X1
X2a_new <- X2a
X1_new[gene.use[i], IndexH_record[[j]]] <- 0
X2a_new[loci_j, IndexH_record[[j]]] <- 0
x1_new <- X1_new[gene.use[i], ]
x2a_new <- X2a_new[loci_j, ]
cors2 <- cor(x1_new, t(x2a))
cors3 <- cor(x1, t(x2a_new))
cors1[is.na(cors1)] <- 0
cors2[is.na(cors2)] <- 0
cors3[is.na(cors3)] <- 0
D <- rbind(cors1 - cors2, cors1 - cors3)
flag <- (rowSums(abs(D) > cutoff) > 0) & abs(cors1) > thresh_corr
regulations_i[[j]]$link <- as.character(loci_j[flag])
regulations_i[[j]]$intensity <- cors1[flag]
}
regulations[[i]] <- regulations_i
}
return(regulations)
}
#' select highly variable features
#'
#' @param object scAI objecy
#' @param assay Name of assay
#' @param do.plot Whether showing plot
#' @param do.text Whether adding feature names
#' @param x.low.cutoff The minimum expression level
#' @param x.high.cutoff The maximum expression level
#' @param y.cutoff The minimum fano factor values
#' @param y.high.cutoff The maximum fano factor values
#' @param num.bin Number of bins
#' @param pch.use Shape of dots in ggplot
#' @param col.use Color of dots
#' @param cex.text.use Size of text
#'
#' @return
#' @export
#'
#' @examples
#' @importFrom graphics smoothScatter text
selectFeatures <- function(object, assay = "RNA", do.plot = TRUE, do.text = TRUE,
x.low.cutoff = 0.01, x.high.cutoff = 3.5, y.cutoff = 0.5, y.high.cutoff = Inf,
num.bin = 20, pch.use = 16, col.use = "black", cex.text.use = 0.5) {
# This function is modified from Seurat Package
data <- object@norm.data[[assay]]
genes.use <- rownames(data)
gene.mean <- apply(X = data, MARGIN = 1, FUN = ExpMean)
gene.dispersion <- apply(X = data, MARGIN = 1, FUN = LogVMR)
gene.dispersion[is.na(x = gene.dispersion)] <- 0
gene.mean[is.na(x = gene.mean)] <- 0
names(x = gene.mean) <- genes.use
data_x_bin <- cut(x = gene.mean, breaks = num.bin)
names(x = data_x_bin) <- names(x = gene.mean)
mean_y <- tapply(X = gene.dispersion, INDEX = data_x_bin, FUN = mean)
sd_y <- tapply(X = gene.dispersion, INDEX = data_x_bin, FUN = sd)
gene.dispersion.scaled <- (gene.dispersion - mean_y[as.numeric(x = data_x_bin)])/sd_y[as.numeric(x = data_x_bin)]
gene.dispersion.scaled[is.na(x = gene.dispersion.scaled)] <- 0
names(x = gene.dispersion.scaled) <- names(x = gene.mean)
mv.df <- data.frame(gene.mean, gene.dispersion, gene.dispersion.scaled)
rownames(x = mv.df) <- rownames(x = data)
hvg.info <- mv.df
names(x = gene.mean) <- names(x = gene.dispersion) <- names(x = gene.dispersion.scaled) <- rownames(hvg.info)
pass.cutoff <- names(x = gene.mean)[which(x = ((gene.mean > x.low.cutoff) & (gene.mean < x.high.cutoff)) & (gene.dispersion.scaled > y.cutoff) & (gene.dispersion.scaled < y.high.cutoff))]
if (do.plot) {
smoothScatter(x = gene.mean, y = gene.dispersion.scaled, pch = pch.use, cex = 0.5, col = col.use, xlab = "Average expression", ylab = "Dispersion", nrpoints = Inf)
}
if (do.text) {
text(x = gene.mean[pass.cutoff], y = gene.dispersion.scaled[pass.cutoff], labels = pass.cutoff, cex = cex.text.use)
}
hvg.info <- hvg.info[order(hvg.info$gene.dispersion, decreasing = TRUE), ]
if (length(object@var.features) == 0) {
for (i in 1:length(object@raw.data)) {
object@var.features[[i]] <- vector("character")
}
names(object@var.features) <- names(object@raw.data)
}
object@var.features[[assay]] <- pass.cutoff
return(object)
}
#' find chromsome regions of genes
#'
#' @param genes a given set of genes
#' @param species mouse or human
#' @importFrom biomaRt useMart getBM
#' @export
#'
searchGeneRegions <- function(genes, species = "mouse") {
if (species == "mouse"){
mouse = biomaRt::useMart("ensembl", dataset = "mmusculus_gene_ensembl")
loci <- biomaRt::getBM(attributes = c( "mgi_symbol","chromosome_name",'start_position','end_position'), filters = "mgi_symbol", values = genes, mart = mouse)
} else{
human = biomaRt::useMart("ensembl", dataset = "hsapiens_gene_ensembl")
loci <- biomaRt::getBM(attributes = c( "hgnc_symbol","chromosome_name",'start_position','end_position'), filters = "hgnc_symbol", values = genes, mart = human)
}
return(loci)
}
#' find the nearby loci in a given loci set for a given set of genes
#'
#' @param genes a given set of genes
#' @param loci the background loci set
#' @param width the width from TSS
#' @param species mouse or human
#' @importFrom bedr bedr.sort.region in.region
#' @export
#'
findNearbyLoci <- function(genes, loci, width = 50000, species = "mouse") {
gene.loci <- searchGeneRegions(genes, species)
gene.loci <- gene.loci[gene.loci$chromosome_name %in% c(1:22, "X","Y"), ]
gene.loci$start_position <- gene.loci$start_position-min(width, gene.loci$start_position)
gene.loci$end_position <- gene.loci$end_position + width
gene.loci.bed <- paste0("chr",gene.loci$chromosome_name, ":", gene.loci$start_position, "-",gene.loci$end_position)
gene.loci.bed.sort <- bedr::bedr.sort.region(gene.loci.bed);
a = grep(paste0("chr", c(1:22, "X","Y"), ":"), loci)
loci <- loci[grepl(paste(paste0("chr", c(1:22, "X","Y"), ":"), collapse="|"), loci)]
loci.sort <- bedr::bedr.sort.region(loci);
is.region <- bedr::in.region(loci.sort, gene.loci.bed.sort);
loci.nearby <- loci.sort[is.region]
return(loci.nearby)
}
#' Select number of the rank
#'
#' @param object scAI object
#' @param rangeK A predefined range of K
#'
#' @return
#' @export
#'
#' @examples
selectK <- function(object, rangeK = c(2:15)) {
coph <- vector("double", length(rangeK))
i <- 0
for (k in rangeK) {
i <- i+1
outs <- run_scAI(object, K = k, nrun = object@options$paras$nrun,
s = object@options$paras$s, alpha = object@options$paras$alpha, lambda = object@options$paras$lambda, gamma = object@options$paras$gamma,
maxIter = object@options$paras$maxIter, stop_rule = object@options$paras$stop_rule)
CM <- outs@cluster$consensus
d1 <- dist(CM)
hc <- hclust(d1, method = "average")
d2 <- cophenetic(hc)
coph[i] <- cor(d1, d2)
}
df <- data.frame(k = rangeK, Coph = coph)
gg <- ggplot(df, aes(x = k, y = Coph)) + geom_line(size=1) +
geom_point() +
theme_classic() + labs(x = 'K', y='Stability score (Coph)') +
scAI_theme_opts() + theme(text = element_text(size = 10)) + labs(x = 'K', y='Stability score (Coph)') +
theme(legend.position = "right", legend.title = NULL)
gg
}
#' Reorder features according to the loading values
#'
#' @param W Basis matrix
#' @param cutoff Threshold of the loading values
#'
#' @return
#' @export
#'
#' @examples
reorderFeatures <- function(W, cutoff = 0.5) {
M <- nrow(W)
K = ncol(W)
MW <- colMeans(W)
vw <- apply(W, 2, sd)
# order features
IndexW_record <- vector("list", K)
IndexW_record[[1]] <- which(W[, 1] > MW[1] + cutoff * VW[1])
n <- length(IndexW_record[[1]])
A <- c(1:M)
c <- match(A, IndexW_record[[1]])
A <- A[-c]
for (j in 2:K) {
IndexW_record[[j]] <- which(W[, j] > MW[j] + cutoff * VW[j])
s <- 0
for (k in 1:j - 1) {
s <- s + length(IndexW_record[[k]])
}
ir2 <- match(IndexW_record[[j]], 1:s)
IndexW_record[[j]] <- IndexW_record[[j]][-ir2]
n <- length(IndexW_record[[j]])
B <- union(1:s, IndexW_record[[j]])
A <- 1:M
c <- match(A, B)
A <- A[-c]
W0 <- W
W[s + 1:s + n, ] <- W0[IndexW_record[[j]], ]
W[s + n + 1:M, ] <- W0[A, ]
}
results <- list()
list[["W"]] <- W
list[["index_record"]] <- IndexW_record
return(results)
}
#' Calculate the variance to mean ratio of logged values
#'
#' Calculate the variance to mean ratio (VMR) in non-logspace (return answer in
#' log-space)
#'
#' @param x A vector of values
#'
#' @return Returns the VMR in log-space
#'
#' @importFrom stats var
#'
#' @export
#'
#' @examples
#' LogVMR(x = c(1, 2, 3))
#'
LogVMR <- function(x) {
return(log(x = var(x = exp(x = x) - 1) / mean(x = exp(x = x) - 1)))
}
#' Calculate the mean of logged values
#'
#' Calculate mean of logged values in non-log space (return answer in log-space)
#'
#' @param x A vector of values
#'
#' @return Returns the mean in log-space
#'
#' @export
#'
#' @examples
#' ExpMean(x = c(1, 2, 3))
#'
ExpMean <- function(x) {
return(log(x = mean(x = exp(x = x) - 1) + 1))
}
#' Convert a sparse matrix to a dense matrix
#'
#' @param mat A sparse matrix
#' @export
as_matrix <- function(mat){
Rcpp::sourceCpp(code='
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
IntegerMatrix asMatrix(NumericVector rp,
NumericVector cp,
NumericVector z,
int nrows,
int ncols){
int k = z.size() ;
IntegerMatrix mat(nrows, ncols);
for (int i = 0; i < k; i++){
mat(rp[i],cp[i]) = z[i];
}
return mat;
}
' )
row_pos <- mat@i
col_pos <- findInterval(seq(mat@x)-1,mat@p[-1])
tmp <- asMatrix(rp = row_pos, cp = col_pos, z = mat@x,
nrows = mat@Dim[1], ncols = mat@Dim[2])
row.names(tmp) <- mat@Dimnames[[1]]
colnames(tmp) <- mat@Dimnames[[2]]
return(tmp)
}
Datasets
Models
