--- 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)
+}
+
+
+
+