diff --git a/R/helper_celltype_expression.R b/R/helper_celltype_expression.R
index 3b5549b..7005c0a 100644
--- a/R/helper_celltype_expression.R
+++ b/R/helper_celltype_expression.R
@@ -79,7 +79,7 @@ remove_mito_ribo <- function(alldata) {
 #' function only calculates the HVG across all cells and returns a vector of HVGs.
 #' @noRd
 find_var_gene <- function(alldata, num_top_gene = 1500, ncores = 1, celltype = TRUE) {
-  BPparam <- generateBPParam(ncores)
+  BPparam <- scFeatures:::generateBPParam(ncores)
 
   if (celltype == TRUE) {
     # here calculates the HVG across all cells across all cell types
@@ -96,9 +96,9 @@ find_var_gene <- function(alldata, num_top_gene = 1500, ncores = 1, celltype = T
 
 
     # below calculates the HVG within each cell type
-    # thiscelltype <- unique( alldata$celltype)[1]
+    # thiscelltype <- unique( alldata$celltype)[28]
     gene <- BiocParallel::bplapply(unique(alldata$celltype), function(thiscelltype) {
-      this_data <- alldata$data[, alldata$celltype == thiscelltype]
+      this_data <- alldata$data[, alldata$celltype == thiscelltype, drop=FALSE ]
       this_data_sample <-  alldata$sample[ alldata$celltype == thiscelltype]
       thisgene <- c()
 
@@ -171,7 +171,7 @@ find_var_gene <- function(alldata, num_top_gene = 1500, ncores = 1, celltype = T
 #' The output is a returns a matrix of samples by features.
 #' @noRd
 helper_gene_mean_celltype <- function( alldata, genes = NULL, num_top_gene = NULL, ncores = 1) {
-  BPparam <- generateBPParam(ncores)
+  BPparam <- scFeatures:::generateBPParam(ncores)
 
 
   if (is.null(num_top_gene)) {
@@ -180,7 +180,7 @@ helper_gene_mean_celltype <- function( alldata, genes = NULL, num_top_gene = NUL
 
 
   if (is.null(genes)) {
-    all_marker <- find_var_gene(alldata ,
+    all_marker <-  find_var_gene(alldata ,
       num_top_gene = num_top_gene,
       ncores = ncores, celltype = TRUE
     )
@@ -337,8 +337,8 @@ helper_gene_cor_celltype <- function(alldata, genes = NULL, num_top_gene = NULL,
   #  thiscelltype <- unique( alldata$celltype) [1]
 
   cor_thiscelltype <- BiocParallel::bplapply(unique(all_marker$celltype), function(thiscelltype) {
-    thisdata <- alldata$data[, alldata$celltype == thiscelltype]
-    thisdata_sample <- alldata$sample[alldata$celltype == thiscelltype]
+    thisdata <- alldata$data[, alldata$celltype == thiscelltype , drop = FALSE ]
+    thisdata_sample <- alldata$sample[alldata$celltype == thiscelltype   ]
     gene <- all_marker[all_marker$celltype == thiscelltype, ]$marker
     thisdata <- thisdata[gene, ]
 
diff --git a/R/helper_proportion.R b/R/helper_proportion.R
index 513bbad..b1d2478 100644
--- a/R/helper_proportion.R
+++ b/R/helper_proportion.R
@@ -30,7 +30,7 @@ helper_proportion_raw <- function( alldata, logit = TRUE) {
     rownames(df) <- df$sample
     df <- df[, -1]
 
-    df <- df[unique(alldata$sample), ]
+    df <- df[unique(alldata$sample),]
 
     return(df)
 }
@@ -43,7 +43,7 @@ helper_proportion_raw <- function( alldata, logit = TRUE) {
 #' @noRd
 helper_proportion_ratio <- function( alldata, ncores = 1) {
   
-    BPparam <- generateBPParam(ncores)
+    BPparam <- scFeatures:::generateBPParam(ncores)
 
     allcelltype <- unique( alldata$celltype)
 
@@ -128,10 +128,10 @@ helper_proportion_ratio <- function( alldata, ncores = 1) {
 #' @noRd
 #' 
 helper_proportion_raw_st <- function(alldata, logit = TRUE, ncores = 1) {
-    BPparam <- generateBPParam(ncores)
+    BPparam <- scFeatures:::generateBPParam(ncores)
 
  
-    num_cell_spot <- get_num_cell_per_celltype(alldata)
+    num_cell_spot <- scFeatures:::get_num_cell_per_celltype(alldata)
 
     prop_table <- BiocParallel::bplapply(unique(alldata$sample), function(s) {
         index <- which(alldata$sample == s)
diff --git a/R/helper_spatial_metrics.R b/R/helper_spatial_metrics.R
index 6e38d18..1adf0e1 100644
--- a/R/helper_spatial_metrics.R
+++ b/R/helper_spatial_metrics.R
@@ -42,7 +42,7 @@ individual_celltype_interaction_sp <- function(this_sample) {
     })
 
     nn_list_cellTypes <- unlist(nn_list_cellTypes)
-    nn_list_cellTypes <- rearrange_string(nn_list_cellTypes)
+    nn_list_cellTypes <- scFeatures:::rearrange_string(nn_list_cellTypes)
     nn_list_cellTypes <- table(nn_list_cellTypes)
 
     return(nn_list_cellTypes)
@@ -57,7 +57,7 @@ individual_celltype_interaction_sp <- function(this_sample) {
 #' @noRd
 helper_celltype_interaction_sp <- function( alldata, ncores = 1) {
   
-    BPparam <- generateBPParam(ncores)
+    BPparam <- scFeatures:::generateBPParam(ncores)
 
     # s <- unique( alldata$sample)[1]
 
@@ -111,9 +111,10 @@ helper_celltype_interaction_sp <- function( alldata, ncores = 1) {
     }
 
 
+    
     rownames(temp) <- temp$nn_list_cellTypes
-    temp <- temp[, -1]
-
+    temp <- temp[, -1, drop=FALSE]
+ 
     colnames(temp) <- unique(alldata$sample)
 
     temp <- t(temp)
@@ -332,10 +333,13 @@ helper_L_stat_st <- function(alldata, ncores = 1) {
         colnames(temp) <- make.names(colnames(temp), unique = TRUE)
     }
 
-
-    rownames(temp) <- temp$rowname
-    temp <- temp[, -1]
-
+    if (ncol(temp == 2)){
+      temp <- temp[, -2, drop=FALSE]
+    }else{
+      rownames(temp) <- temp$rowname
+      temp <- temp[, -1]
+    } 
+     
     colnames(temp) <- unique(alldata$sample)
 
     temp <- t(temp)
@@ -373,7 +377,7 @@ individual_L_stat_sp <- function(this_sample) {
 
     L_patient <- list()
     for (i in seq_len(nrow(cellTypes_pair))) {
-        L_patient[[i]] <- L_stats(cell_points,
+        L_patient[[i]] <- scFeatures:::L_stats(cell_points,
             from = cellTypes_pair[i, 1],
             to = cellTypes_pair[i, 2],
             L_dist = 50
@@ -442,15 +446,21 @@ helper_L_stat_sp <- function( alldata, ncores = 1) {
     }
 
 
-    rownames(temp) <- temp$rowname
-    temp <- temp[, -1]
+    
+    if (ncol(temp == 2)){
+      temp <- temp[, -2, drop=FALSE]
+    }else{
+      rownames(temp) <- temp$rowname
+      temp <- temp[, -1]
+    } 
+      colnames(temp) <- unique(alldata$sample)
+      
+      temp <- t(temp)
+      
+      temp[is.na(temp)] <- 0
+      L_patient <- temp
 
-    colnames(temp) <- unique(alldata$sample)
 
-    temp <- t(temp)
-
-    temp[is.na(temp)] <- 0
-    L_patient <- temp
 
     return(L_patient)
 }
@@ -569,9 +579,14 @@ helper_nncorr_protein <- function(alldata, num_top_gene = NULL, ncores = 1) {
     }
 
 
-    rownames(temp) <- temp$rowname
-    temp <- temp[, -1]
-
+    if (ncol(temp == 2)){
+      temp <- temp[, -2, drop=FALSE]
+    }else{
+      rownames(temp) <- temp$rowname
+      temp <- temp[, -1]
+    } 
+    
+ 
     colnames(temp) <- unique(alldata$sample)
 
     temp <- t(temp)
@@ -698,10 +713,14 @@ helper_moran <- function( alldata, num_top_gene = NULL, ncores = 1) {
         colnames(temp) <- make.names(colnames(temp), unique = TRUE)
     }
 
-
-    rownames(temp) <- temp$rowname
-    temp <- temp[, -1]
-
+    if (ncol(temp == 2)){
+      temp <- temp[, -2, drop=FALSE]
+    }else{
+      rownames(temp) <- temp$rowname
+      temp <- temp[, -1]
+    } 
+    
+ 
     colnames(temp) <- unique(alldata$sample)
 
     temp <- t(temp)
diff --git a/R/run_scfeatures.R b/R/run_scfeatures.R
index a84ee4e..30f309e 100644
--- a/R/run_scfeatures.R
+++ b/R/run_scfeatures.R
@@ -48,9 +48,9 @@ run_proportion_raw <- function(data, type = "scrna", ncores = 1) {
                 "for feature `proportion_raw`"
             )
         }
-        X <- helper_proportion_raw(data, logit = FALSE)
+        X <- scFeatures:::helper_proportion_raw(data, logit = FALSE)
     } else if (type == "spatial_t") {
-        X <- helper_proportion_raw_st(data, logit = FALSE, ncores)
+        X <- scFeatures:::helper_proportion_raw_st(data, logit = FALSE, ncores)
     } else {
         cli::cli_abort(c(
             "Parameter {.var type} must be 'scrna', 'spatial_p' or 'spatial_t'",
@@ -112,7 +112,7 @@ run_proportion_logit <- function(data, type = "scrna", ncores = 1) {
                 "for feature proportion_logit."
             )
         }
-        X <- helper_proportion_raw(data, logit = TRUE)
+        X <- scFeatures:::helper_proportion_raw(data, logit = TRUE)
     } else if (type == "spatial_t") {
         X <- helper_proportion_raw_st(data, logit = TRUE, ncores)
     } else {