Datasets
Models
Diff of
/R/scAI_model.R
[000000]
..
[4a0329]
Switch to side-by-side view
--- a +++ b/R/scAI_model.R @@ -0,0 +1,1691 @@ +#' 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) +} + + + +