diff --git a/.github/workflows/check-bioc.yaml b/.github/workflows/check-bioc.yaml
index 40d726f..651dd25 100644
--- a/.github/workflows/check-bioc.yaml
+++ b/.github/workflows/check-bioc.yaml
@@ -305,7 +305,7 @@ jobs:
 
       - name: Upload check results
         if: failure()
-        uses: actions/upload-artifact@devel
+        uses: actions/upload-artifact@master
         with:
           name: ${{ runner.os }}-biocversion-RELEASE_3_18-r-4.3-results
           path: check
diff --git a/DESCRIPTION b/DESCRIPTION
index 52e4630..311cf2a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -39,7 +39,8 @@ Imports:
     mixtools,
     sparseMatrixStats,
     stats,
-    tidyr
+    tidyr,
+    utils
 Suggests: 
     BiocStyle,
     dbscan,
diff --git a/NAMESPACE b/NAMESPACE
index d44768b..9dfdada 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -20,3 +20,7 @@ import(graphics)
 import(methods)
 import(stats)
 importClassesFrom(Matrix,dgCMatrix)
+importFrom(SummarizedExperiment,SummarizedExperiment)
+importFrom(SummarizedExperiment,assay)
+importFrom(utils,globalVariables)
+importFrom(utils,stack)
diff --git a/R/AllGenerics.R b/R/AllGenerics.R
index 403ae13..d1bd51c 100644
--- a/R/AllGenerics.R
+++ b/R/AllGenerics.R
@@ -91,6 +91,8 @@ setGeneric(
 #' @param data an expression object, can be matrix or SummarizedExperiment
 #' @param slot a character, specify which slot to use when data is se object,
 #'             optional, default 'score'
+#' @param suffix a character, specify the name suffix to save score when
+#'               features is a named list
 #'
 #' @return A vector of overall score for each sample
 #'
@@ -98,7 +100,7 @@ setGeneric(
 #' @export
 #'
 #' @examples
-#' data <- matrix(rnorm(100), 10, dimnames = list(seq_len(3)))
+#' data <- matrix(rnorm(100), 10, dimnames = list(seq_len(10)))
 #' gs_score(data, features = seq_len(3))
 setGeneric(
   "gs_score",
diff --git a/R/plot.R b/R/plot.R
index 07e1d09..9c45548 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -175,3 +175,5 @@ score_barplot <- function(top_markers, column = ".dot", f_list, n = 30) {
     theme_classic()
 }
 
+utils::globalVariables(c("Group", "Score", "x", "Scores", "Comp",
+                         "..density..", "stack", "values", "ind", "Type"))
diff --git a/R/score.R b/R/score.R
index a13471a..37584ae 100644
--- a/R/score.R
+++ b/R/score.R
@@ -57,7 +57,9 @@ idf_iae_methods <- function() {
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' cal_score_init(data, label = sample(c("A", "B"), 10, replace = TRUE))
+#' label <- sample(c("A", "B"), 10, replace = TRUE)
+#' smartid:::cal_score_init(data, par.idf = list(label = label),
+#'                          par.iae = list(label = label))
 cal_score_init <- function(expr, tf = c("logtf", "tf"),
                            idf = "prob", iae = "prob",
                            par.idf = NULL, par.iae = NULL) {
diff --git a/R/smartid-package.R b/R/smartid-package.R
index d821bd7..ff6e4eb 100644
--- a/R/smartid-package.R
+++ b/R/smartid-package.R
@@ -2,6 +2,8 @@
 #' @import methods
 #' @import stats
 #' @import graphics
+#' @importFrom SummarizedExperiment SummarizedExperiment assay
+#' @importFrom utils globalVariables stack
 #' @importClassesFrom Matrix dgCMatrix
 NULL
 
diff --git a/R/tf_idf_iae_wrappers.R b/R/tf_idf_iae_wrappers.R
index 5ee44bb..c811f02 100644
--- a/R/tf_idf_iae_wrappers.R
+++ b/R/tf_idf_iae_wrappers.R
@@ -12,7 +12,7 @@
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' tf(data)
+#' smartid:::tf(data)
 tf <- function(expr, log = FALSE) {
   t.f <- sweep(expr, 2, colSums(expr, na.rm = TRUE) + 0.01, FUN = "/")
   if(log) {
@@ -40,7 +40,7 @@ tf <- function(expr, log = FALSE) {
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' idf(data)
+#' smartid:::idf(data)
 idf <- function(expr, features = NULL, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
   n_obs <- ncol(expr)  ## number of total obs
@@ -65,7 +65,7 @@ idf <- function(expr, features = NULL, thres = 0) {
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' idf_m(data)
+#' smartid:::idf_m(data)
 idf_m <- function(expr, features = NULL, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
   n_obs <- ncol(expr)  ## number of total obs
@@ -95,7 +95,7 @@ idf_m <- function(expr, features = NULL, thres = 0) {
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' idf_sd(data)
+#' smartid:::idf_sd(data)
 idf_sd <- function(expr, features = NULL, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
   n_obs <- ncol(expr)  ## number of total obs
@@ -120,8 +120,9 @@ idf_sd <- function(expr, features = NULL, thres = 0) {
 #' @return a matrix of IDF score
 #'
 #' @examples
+#' set.seed(123)
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' idf_hdb(data)
+#' smartid:::idf_hdb(data)
 idf_hdb <- function(expr, features = NULL, multi = TRUE,
                     thres = 0, minPts = 2, ...) {
   if(is.null(features)) features <- seq_len(nrow(expr))
@@ -159,7 +160,7 @@ idf_hdb <- function(expr, features = NULL, multi = TRUE,
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' idf_rf(data, label = sample(c("A", "B"), 10, replace = TRUE))
+#' smartid:::idf_rf(data, label = sample(c("A", "B"), 10, replace = TRUE))
 idf_rf <- function(expr, features = NULL, label,
                    multi = TRUE, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
@@ -202,7 +203,7 @@ idf_rf <- function(expr, features = NULL, label,
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' idf_prob(data, label = sample(c("A", "B"), 10, replace = TRUE))
+#' smartid:::idf_prob(data, label = sample(c("A", "B"), 10, replace = TRUE))
 idf_prob <- function(expr, features = NULL, label,
                      multi = TRUE, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
@@ -246,7 +247,7 @@ idf_prob <- function(expr, features = NULL, label,
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' idf_igm(data, label = sample(c("A", "B"), 10, replace = TRUE))
+#' smartid:::idf_igm(data, label = sample(c("A", "B"), 10, replace = TRUE))
 idf_igm <- function(expr, features = NULL, label, lambda = 7, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
 
@@ -284,7 +285,7 @@ idf_igm <- function(expr, features = NULL, label, lambda = 7, thres = 0) {
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' iae(data)
+#' smartid:::iae(data)
 iae <- function(expr, features = NULL, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
   n_obs <- ncol(expr)  ## number of total obs
@@ -311,7 +312,7 @@ iae <- function(expr, features = NULL, thres = 0) {
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' iae_m(data)
+#' smartid:::iae_m(data)
 iae_m <- function(expr, features = NULL, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
   n_obs <- ncol(expr)  ## number of total obs
@@ -342,7 +343,7 @@ iae_m <- function(expr, features = NULL, thres = 0) {
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' iae_sd(data)
+#' smartid:::iae_sd(data)
 iae_sd <- function(expr, features = NULL, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
   n_obs <- ncol(expr)  ## number of obs
@@ -369,8 +370,9 @@ iae_sd <- function(expr, features = NULL, thres = 0) {
 #' @return a matrix of IAE score
 #'
 #' @examples
+#' set.seed(123)
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' iae_hdb(data)
+#' smartid:::iae_hdb(data)
 iae_hdb <- function(expr, features = NULL, multi = TRUE,
                     thres = 0, minPts = 2, ...) {
   if(is.null(features)) features <- seq_len(nrow(expr))
@@ -420,7 +422,7 @@ iae_hdb <- function(expr, features = NULL, multi = TRUE,
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' iae_rf(data, label = sample(c("A", "B"), 10, replace = TRUE))
+#' smartid:::iae_rf(data, label = sample(c("A", "B"), 10, replace = TRUE))
 iae_rf <- function(expr, features = NULL, label,
                    multi = TRUE, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
@@ -463,7 +465,7 @@ iae_rf <- function(expr, features = NULL, label,
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' iae_prob(data, label = sample(c("A", "B"), 10, replace = TRUE))
+#' smartid:::iae_prob(data, label = sample(c("A", "B"), 10, replace = TRUE))
 iae_prob <- function(expr, features = NULL, label,
                      multi = TRUE, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
@@ -506,7 +508,7 @@ iae_prob <- function(expr, features = NULL, label,
 #'
 #' @examples
 #' data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-#' iae_igm(data, label = sample(c("A", "B"), 10, replace = TRUE))
+#' smartid:::iae_igm(data, label = sample(c("A", "B"), 10, replace = TRUE))
 iae_igm <- function(expr, features = NULL, label, lambda = 7, thres = 0) {
   if(is.null(features)) features <- seq_len(nrow(expr))
 
diff --git a/R/top_markers.R b/R/top_markers.R
index d09bf19..f7f55b3 100644
--- a/R/top_markers.R
+++ b/R/top_markers.R
@@ -2,7 +2,7 @@
 #'
 #' @inheritParams top_markers_abs
 #' @inheritParams top_markers_glm
-#' @param use_glm logical, if to use [stats::glm()] to compute group mean score,
+#' @param use.glm logical, if to use [stats::glm()] to compute group mean score,
 #'                if TRUE, also compute mean score difference as output
 #' @param ... params for [top_markers_abs()] or [top_markers_glm()]
 #'
@@ -190,3 +190,4 @@ softmax <- function(x, tau = 1) {
 ## tanh: [-1, 1], similar to sigmoid, no need to sum 1
 tanh <- function(x) 2 / (1 + exp(-2*x)) - 1
 
+utils::globalVariables(c(".dot", "Scores"))
diff --git a/man/cal_score_init.Rd b/man/cal_score_init.Rd
index 205ece3..71591e6 100644
--- a/man/cal_score_init.Rd
+++ b/man/cal_score_init.Rd
@@ -36,5 +36,7 @@ Calculate score for each feature in each cell
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-cal_score_init(data, label = sample(c("A", "B"), 10, replace = TRUE))
+label <- sample(c("A", "B"), 10, replace = TRUE)
+smartid:::cal_score_init(data, par.idf = list(label = label),
+                         par.iae = list(label = label))
 }
diff --git a/man/gs_score.Rd b/man/gs_score.Rd
index 6b569c0..e63b310 100644
--- a/man/gs_score.Rd
+++ b/man/gs_score.Rd
@@ -22,6 +22,9 @@ gs_score(data, features = NULL, slot = "score", suffix = "score")
 
 \item{slot}{a character, specify which slot to use when data is se object,
 optional, default 'score'}
+
+\item{suffix}{a character, specify the name suffix to save score when
+features is a named list}
 }
 \value{
 A vector of overall score for each sample
@@ -30,6 +33,6 @@ A vector of overall score for each sample
 compute overall score based on the given marker list
 }
 \examples{
-data <- matrix(rnorm(100), 10, dimnames = list(seq_len(3)))
+data <- matrix(rnorm(100), 10, dimnames = list(seq_len(10)))
 gs_score(data, features = seq_len(3))
 }
diff --git a/man/iae.Rd b/man/iae.Rd
index e1eb79a..78f4ef8 100644
--- a/man/iae.Rd
+++ b/man/iae.Rd
@@ -21,5 +21,5 @@ standard inverse average expression
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-iae(data)
+smartid:::iae(data)
 }
diff --git a/man/iae_hdb.Rd b/man/iae_hdb.Rd
index a239a7e..1200f79 100644
--- a/man/iae_hdb.Rd
+++ b/man/iae_hdb.Rd
@@ -27,6 +27,7 @@ a matrix of IAE score
 inverse average expression using hdbscan cluster as label
 }
 \examples{
+set.seed(123)
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-iae_hdb(data)
+smartid:::iae_hdb(data)
 }
diff --git a/man/iae_igm.Rd b/man/iae_igm.Rd
index e42b56a..b166999 100644
--- a/man/iae_igm.Rd
+++ b/man/iae_igm.Rd
@@ -25,5 +25,5 @@ labeled inverse average expression IGM
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-iae_igm(data, label = sample(c("A", "B"), 10, replace = TRUE))
+smartid:::iae_igm(data, label = sample(c("A", "B"), 10, replace = TRUE))
 }
diff --git a/man/iae_m.Rd b/man/iae_m.Rd
index e853a7a..d646156 100644
--- a/man/iae_m.Rd
+++ b/man/iae_m.Rd
@@ -21,5 +21,5 @@ inverse average expression: max
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-iae_m(data)
+smartid:::iae_m(data)
 }
diff --git a/man/iae_prob.Rd b/man/iae_prob.Rd
index 09daf97..aab7677 100644
--- a/man/iae_prob.Rd
+++ b/man/iae_prob.Rd
@@ -25,5 +25,5 @@ labeled inverse average expression: probability based
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-iae_prob(data, label = sample(c("A", "B"), 10, replace = TRUE))
+smartid:::iae_prob(data, label = sample(c("A", "B"), 10, replace = TRUE))
 }
diff --git a/man/iae_rf.Rd b/man/iae_rf.Rd
index 6534357..2f35563 100644
--- a/man/iae_rf.Rd
+++ b/man/iae_rf.Rd
@@ -25,5 +25,5 @@ labeled inverse average expression: relative frequency
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-iae_rf(data, label = sample(c("A", "B"), 10, replace = TRUE))
+smartid:::iae_rf(data, label = sample(c("A", "B"), 10, replace = TRUE))
 }
diff --git a/man/iae_sd.Rd b/man/iae_sd.Rd
index 5648f33..23adfa2 100644
--- a/man/iae_sd.Rd
+++ b/man/iae_sd.Rd
@@ -21,5 +21,5 @@ inverse average expression using standard deviation (SD)
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-iae_sd(data)
+smartid:::iae_sd(data)
 }
diff --git a/man/idf.Rd b/man/idf.Rd
index 8281ffa..a7a26d9 100644
--- a/man/idf.Rd
+++ b/man/idf.Rd
@@ -21,5 +21,5 @@ standard inverse cell frequency
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-idf(data)
+smartid:::idf(data)
 }
diff --git a/man/idf_hdb.Rd b/man/idf_hdb.Rd
index f98216e..d329845 100644
--- a/man/idf_hdb.Rd
+++ b/man/idf_hdb.Rd
@@ -27,6 +27,7 @@ a matrix of IDF score
 Title
 }
 \examples{
+set.seed(123)
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-idf_hdb(data)
+smartid:::idf_hdb(data)
 }
diff --git a/man/idf_igm.Rd b/man/idf_igm.Rd
index 3a56ee8..c6c6f99 100644
--- a/man/idf_igm.Rd
+++ b/man/idf_igm.Rd
@@ -25,5 +25,5 @@ labeled inverse document frequency IGM
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-idf_igm(data, label = sample(c("A", "B"), 10, replace = TRUE))
+smartid:::idf_igm(data, label = sample(c("A", "B"), 10, replace = TRUE))
 }
diff --git a/man/idf_m.Rd b/man/idf_m.Rd
index 65c36b4..c818958 100644
--- a/man/idf_m.Rd
+++ b/man/idf_m.Rd
@@ -21,5 +21,5 @@ inverse document frequency: max
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-idf_m(data)
+smartid:::idf_m(data)
 }
diff --git a/man/idf_prob.Rd b/man/idf_prob.Rd
index 0d4614c..f0d014d 100644
--- a/man/idf_prob.Rd
+++ b/man/idf_prob.Rd
@@ -25,5 +25,5 @@ labeled inverse cell frequency: probability based
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-idf_prob(data, label = sample(c("A", "B"), 10, replace = TRUE))
+smartid:::idf_prob(data, label = sample(c("A", "B"), 10, replace = TRUE))
 }
diff --git a/man/idf_rf.Rd b/man/idf_rf.Rd
index 55f3de2..e08ce2b 100644
--- a/man/idf_rf.Rd
+++ b/man/idf_rf.Rd
@@ -25,5 +25,5 @@ labeled inverse cell frequency: relative frequency
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-idf_rf(data, label = sample(c("A", "B"), 10, replace = TRUE))
+smartid:::idf_rf(data, label = sample(c("A", "B"), 10, replace = TRUE))
 }
diff --git a/man/idf_sd.Rd b/man/idf_sd.Rd
index f74dde9..08f0c65 100644
--- a/man/idf_sd.Rd
+++ b/man/idf_sd.Rd
@@ -21,5 +21,5 @@ inverse cell frequency using standard deviation (SD)
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-idf_sd(data)
+smartid:::idf_sd(data)
 }
diff --git a/man/tf.Rd b/man/tf.Rd
index 44a7d83..2677116 100644
--- a/man/tf.Rd
+++ b/man/tf.Rd
@@ -19,5 +19,5 @@ compute term/gene frequency within each cell
 }
 \examples{
 data <- matrix(rpois(100, 2), 10, dimnames = list(1:10))
-tf(data)
+smartid:::tf(data)
 }
diff --git a/man/top_markers.Rd b/man/top_markers.Rd
index c620c4c..f5afe94 100644
--- a/man/top_markers.Rd
+++ b/man/top_markers.Rd
@@ -49,6 +49,9 @@ top_markers(
 
 \item{n}{integer, number of returned top genes for each group}
 
+\item{use.glm}{logical, if to use \code{\link[stats:glm]{stats::glm()}} to compute group mean score,
+if TRUE, also compute mean score difference as output}
+
 \item{scale}{logical, if to scale data by row}
 
 \item{use.mgm}{logical, if to scale data using \code{\link[=scale_mgm]{scale_mgm()}}}
diff --git a/man/top_markers_init.Rd b/man/top_markers_init.Rd
index 819554b..32e1fd0 100644
--- a/man/top_markers_init.Rd
+++ b/man/top_markers_init.Rd
@@ -22,6 +22,9 @@ top_markers_init(
 
 \item{n}{integer, number of returned top genes for each group}
 
+\item{use.glm}{logical, if to use \code{\link[stats:glm]{stats::glm()}} to compute group mean score,
+if TRUE, also compute mean score difference as output}
+
 \item{scale}{logical, if to scale data by row}
 
 \item{use.mgm}{logical, if to scale data using \code{\link[=scale_mgm]{scale_mgm()}}}
@@ -29,9 +32,6 @@ top_markers_init(
 \item{softmax}{logical, if to apply softmax transformation on output}
 
 \item{...}{params for \code{\link[=top_markers_abs]{top_markers_abs()}} or \code{\link[=top_markers_glm]{top_markers_glm()}}}
-
-\item{use_glm}{logical, if to use \code{\link[stats:glm]{stats::glm()}} to compute group mean score,
-if TRUE, also compute mean score difference as output}
 }
 \value{
 a tibble with feature names, group labels and ordered processed scores