diff --git a/R/cluster_gene.R b/R/cluster_gene.R
index d50ca7f..2a88920 100644
--- a/R/cluster_gene.R
+++ b/R/cluster_gene.R
@@ -2,7 +2,7 @@
 #'
 #' This function is used to support the external function clusterGene(). cluster genes based on their temporal patterns  of gene expression (constant time test) or the temporal patterns of gene expression group difference (sample covariate test).
 #'
-#' @import ggplot2 RColorBrewer splines gridExtra viridis
+#' @import ggplot2 RColorBrewer splines gridExtra viridis parallel
 #' @return a plot
 #' @author Wenpin Hou <whou10@jhu.edu>
 #' @param testobj object returned from lamian.test().
@@ -19,7 +19,8 @@ cluster_gene <- function(testobj,
                          type = 'Time',
                          method = 'kmeans', 
                          scale.difference = F,
-                         seed = 12345){
+                         seed = 12345,
+                         ncores = detectCores()){
   if (toupper(type) == 'TIME'){ 
     if ('populationFit' %in% names(testobj)) {
       fit <- testobj$populationFit
@@ -44,7 +45,7 @@ cluster_gene <- function(testobj,
     set.seed(seed)
     # 
     if (k.auto){
-      clu <- mykmeans(mat.scale, maxclunum = 20)$cluster
+      clu <- mykmeans(mat.scale, maxclunum = 20, ncores = ncores)$cluster
     } else {
       clu <- kmeans(mat.scale, k, iter.max = 1000)$cluster
     }
diff --git a/R/evaluate_uncertainty.R b/R/evaluate_uncertainty.R
index 1bba75c..a031f6a 100644
--- a/R/evaluate_uncertainty.R
+++ b/R/evaluate_uncertainty.R
@@ -2,6 +2,7 @@
 #'
 #' This function is designed to evaluate the pseudotime tree uncertainty, as one of the main functions in lamain module 1.
 #'
+#' @import parallel
 #' @param inferobj the output object from function infer_tree_structure().
 #' @param n.permute: a numeric number of permutation in the permutation test.
 #' @param subset.cell a character vector of the names of the selected cells where boostrap will happen on. If NULL, then boostrap from all the cells.
@@ -18,7 +19,8 @@ evaluate_uncertainty <-
            n.permute,
            subset.cell = NULL,
            design = NULL,
-           return.ctcomp = FALSE 
+           return.ctcomp = FALSE,
+           ncores = detectCores()
            # branchPropTest.method = 'ttest', ## this is a quick way to call the t-test in branchPropTest(); however, if want to use the multinom test, call branchPropTest() separately. 
            # branchPropTest.value.log = FALSE
   ) {
@@ -44,7 +46,7 @@ evaluate_uncertainty <-
       
       ## cluster cells
       invisible(capture.output(clu <-
-        mykmeans(pr.pm, number.cluster = max(inferobj$clusterid))$cluster))
+        mykmeans(pr.pm, number.cluster = max(inferobj$clusterid))$cluster), ncores = ncores)
       
       
       ## build pseudotime
diff --git a/R/infer_tree_structure.R b/R/infer_tree_structure.R
index dd06e93..3f3e307 100644
--- a/R/infer_tree_structure.R
+++ b/R/infer_tree_structure.R
@@ -5,7 +5,7 @@
 #' @author Wenpin Hou <whou10@jhu.edu>
 #' @return a list
 #' @export
-#' @import TSCAN scattermore RColorBrewer grDevices
+#' @import TSCAN scattermore RColorBrewer grDevices parallel
 #' @importFrom grDevices pdf
 #' @importFrom grDevices dev.off
 #' @param  pca cell by principal component (pc) matrix. Principal components reduction of the cells.
@@ -35,7 +35,8 @@ infer_tree_structure <-
            xlab = 'PC1',
            ylab = 'PC2',
            max.clunum = 50,
-           kmeans.seed = 12345) {
+           kmeans.seed = 12345,
+           ncores = detectCores()) {
     alls <- cellanno[, 2]
     names(alls) <- cellanno[, 1]
     ## set.seed(12345)
@@ -50,7 +51,7 @@ infer_tree_structure <-
     
     ## clustering
     clu <-
-      mykmeans(pr, maxclunum = 50, number.cluster = number.cluster, seed = kmeans.seed)$cluster
+      mykmeans(pr, maxclunum = 50, number.cluster = number.cluster, seed = kmeans.seed, ncores = ncores)$cluster
     table(clu)
     pd = data.frame(x = pr[, 1],
                     y = pr[, 2],
diff --git a/R/mykmeans.R b/R/mykmeans.R
index 09a7063..5b088eb 100644
--- a/R/mykmeans.R
+++ b/R/mykmeans.R
@@ -16,14 +16,15 @@ mykmeans <-
   function(matrix,
            number.cluster = NA,
            maxclunum = 30,
-           seed = 12345) {
+           seed = 12345,
+           ncores = detectCores()) {
     library(parallel)
     if (is.na(number.cluster)) {
       rss <- mclapply(seq_len(maxclunum), function(clunum) {
         ## set.seed(seed)
         tmp <- kmeans(matrix, clunum, iter.max = 1000)
         tmp$betweenss / tmp$totss
-      }, mc.cores = 30)
+      }, mc.cores = ncores)
       rss <- unlist(rss)
       # number.cluster <- which(diff(rss) < 1e-2)[1]
       x <- 2:maxclunum