diff --git a/DESCRIPTION b/DESCRIPTION
index aae2c7d..7cd05bb 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: cytosignal
 Title: What the Package Does (One Line, Title Case)
-Version: 0.0.0.9000
+Version: 0.2.0
 Type: Package
 Authors@R: 
     person(given = "Jialin",
diff --git a/NAMESPACE b/NAMESPACE
index 61c4e51..54edeb5 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -10,10 +10,13 @@ S3method(graphNicheLR,CytoSignal)
 S3method(graphNicheLR,dgCMatrix)
 S3method(imputeNiche,CytoSignal)
 S3method(imputeNiche,dgCMatrix)
+S3method(imputeNicheVelo,CytoSignal)
 S3method(inferScoreLR,CytoSignal)
 S3method(inferScoreLR,dgCMatrix)
 S3method(inferSignif,CytoSignal)
 S3method(inferSignif,matrix_like)
+S3method(inferVeloLR,CytoSignal)
+S3method(inferVeloLR,matrix_like)
 S3method(normCounts,CytoSignal)
 S3method(normCounts,dgCMatrix)
 S3method(normCounts,list)
@@ -31,10 +34,13 @@ export(hex_bin)
 export(hex_coord)
 export(hex_pos)
 export(imputeNiche)
+export(imputeNicheVelo)
 export(inferEpsParams)
 export(inferScoreLR)
 export(inferSignif)
+export(inferVeloLR)
 export(normCounts)
+export(permuteImpEdge)
 export(permuteImpLR)
 export(permuteListImpLR)
 export(permuteListScoreLR)
@@ -56,6 +62,9 @@ exportClasses(ImpData)
 exportClasses(lrScores)
 exportClasses(lrVelo)
 exportMethods(show)
+importFrom(Matrix,colSums)
+importFrom(Matrix,rowSums)
+importFrom(Matrix,t)
 importFrom(RTriangle,pslg)
 importFrom(RTriangle,triangulate)
 importFrom(Rcpp,evalCpp)
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 109d6c9..e864b5f 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -57,6 +57,10 @@ graphMeanLR_cpp <- function(alpha_list, dge_lig, lig_index, lig_list, nb_index,
     .Call('_cytosignal_graphMeanLR_cpp', PACKAGE = 'cytosignal', alpha_list, dge_lig, lig_index, lig_list, nb_index, nb_list)
 }
 
+inferVeloLR_cpp <- function(dge_lig, dge_recep, dge_lig_velo, dge_recep_velo, lig_index, lig_list, recep_index, recep_list) {
+    .Call('_cytosignal_inferVeloLR_cpp', PACKAGE = 'cytosignal', dge_lig, dge_recep, dge_lig_velo, dge_recep_velo, lig_index, lig_list, recep_index, recep_list)
+}
+
 VelographNicheLR_cpp <- function(dge_lig, dge_recep, dge_velo, lig_index, lig_list, recep_index, recep_list, nb_index, nb_list) {
     .Call('_cytosignal_VelographNicheLR_cpp', PACKAGE = 'cytosignal', dge_lig, dge_recep, dge_velo, lig_index, lig_list, recep_index, recep_list, nb_index, nb_list)
 }
diff --git a/R/analysis.r b/R/analysis.r
index c060744..3c6cbda 100644
--- a/R/analysis.r
+++ b/R/analysis.r
@@ -159,9 +159,11 @@ findNNGauEB.CytoSignal <- function(
 		# imp.data.null = new("dgCMatrix"),
 		imp.data.null = list(),
 		# intr.data = new("dgCMatrix"),
+		imp.velo = new("matrix"),
 		nn.id = nn$id,
 		nn.dist = nn$dist,
 		scale.fac = new("numeric"),
+		scale.fac.velo = new("numeric"),
 		log = list(
 			"Parameters" = paste0("eps: ", eps, ", sigma: ", sigma),
 			"Num of neighbors" = paste0("Mean: ", mean(table(nn$id)), ", Median: ", median(table(nn$id)))
@@ -275,9 +277,11 @@ findNNDT.CytoSignal <- function(
 		# imp.data.null = new("dgCMatrix"),
 		imp.data.null = list(),
 		# intr.data = new("dgCMatrix"),
+		imp.velo = new("matrix"),
 		nn.id = nn$id,
 		nn.dist = nn$dist,
 		scale.fac = new("numeric"),
+		scale.fac.velo = new("numeric"),
 		log = list(
 			"Parameters" = "Delauany Triangulation",
 			"Num of neighbors" = paste0("Mean: ", mean(table(nn$id)), ", Median: ", median(table(nn$id)))
@@ -304,18 +308,28 @@ findNNRaw <- function(
 
 	message("Setting Imp obj using NO imputation...")
 
-	cells.loc <- object@cells.loc
-	nn <- findNNDT.matrix(cells.loc)
+	if ("DT" %in% names(object@imputation)){
+		message("DT has been done before, taking the same neighbors.")
+		nn <- new("list")
+		nn$id <- object@imputation[["DT"]]@nn.id
+		nn$dist <- object@imputation[["DT"]]@nn.dist
+	} else {
+		cells.loc <- object@cells.loc
+		nn <- findNNDT.matrix(cells.loc)
+	}
 
 	nn.obj <- methods::new(
 		"ImpData",
 		method = tag,
 		imp.data = object@counts,
-		# imp.norm.data = new("dgCMatrix"),
-		imp.norm.data = list(),
+		imp.data.null = new("dgCMatrix"),
+		# imp.norm.data = new("list"),
+		imp.velo = new("matrix"),
 		# intr.data = new("dgCMatrix"),
 		nn.id = nn$id,
 		nn.dist = nn$dist,
+		scale.fac = object@parameters[["lib.size"]],
+		scale.fac.velo = new("numeric"),
 		log = list(
 			"Parameters" = "Raw data without imputation",
 			"Num of neighbors" = paste0("Mean: ", mean(table(nn$id)), ", Median: ", median(table(nn$id)))
@@ -358,7 +372,8 @@ imputeNiche.dgCMatrix <- function(
 	dge.raw, 
 	nb.id.fac,
 	nb.dist.fac,
-	weights = c("mean", "counts", "dist", "none")
+	weights = c("mean", "counts", "dist", "none"),
+	return.graph = F
 ){
 	cat("Imputing...\n")
 
@@ -389,11 +404,17 @@ imputeNiche.dgCMatrix <- function(
 		weights.mtx@x = weights.mtx@x / rep.int(Matrix::colSums(weights.mtx), diff(weights.mtx@p))
 	}
 
+	if (return.graph){
+		return(weights.mtx)
+	}
+
 	dge.raw.imputed = dge.raw %*% weights.mtx
 
 	return(dge.raw.imputed)
 }
 
+
+
 #' Sub function for imputeNiche, input a Cytosignal object
 #' 
 #' @param object A Cytosignal object
@@ -441,7 +462,8 @@ imputeNiche.CytoSignal <- function(object,
 
 	# weighted sum of scale factors as well
 	scale.fac.imp <- imputeNiche.dgCMatrix(
-		Matrix(object@parameters$lib.size, sparse = T, nrow = 1), 
+		Matrix(object@parameters$lib.size, sparse = T, nrow = 1,
+					dimnames = list(NULL, colnames(dge.raw))), 
 		nb.id.fac,
 		nb.dist.fac,
 		weights = weights
@@ -460,6 +482,100 @@ imputeNiche.CytoSignal <- function(object,
 	return(object)
 }
 
+
+
+#' Impute the velocity mtx using the specified method
+#' 
+#' @param object A Cytosignal object
+#' @param method The method to use for imputation
+#' @param ... Other parameters
+#' 
+#' @return A Cytosignal object
+#' @export
+imputeNicheVelo <- function(
+  object,
+  ...
+) {
+  UseMethod(generic = 'imputeNicheVelo', object = object)
+}
+
+
+
+#' Sub function for imputeNicheVelo, input a Cytosignal object
+#' 
+#' @param object A Cytosignal object
+#' @param nn.type The type of neighbors
+#' @param weights The weight of the Delaunay triangulation
+#' 
+#' @return A Cytosignal object
+#' @export
+imputeNicheVelo.CytoSignal <- function(
+	object, 
+	nn.type = NULL, 
+	weights = c("mean", "counts", "dist", "none")
+){
+	# dge.raw <- object@counts
+	# Extract neighborhood factors
+
+	if (is.null(nn.type)){
+		nn.type <- object@imputation[["default"]]
+	}
+
+	if  (!weights %in% c("mean", "counts", "dist", "none")){
+		stop("Incorret weighting method!\n")
+	}
+
+	message(paste0("Imputing using ", nn.type, "..."))
+
+	dge.raw <- object@velo[["velo.s"]] + object@velo[["velo.u"]]
+
+	if (nn.type %in% names(object@imputation)){
+		nb.id.fac <- object@imputation[[nn.type]]@nn.id
+		nb.dist.fac <- object@imputation[[nn.type]]@nn.dist
+	} else {
+		stop("NN type not found.")
+	}
+
+	if (nn.type == "Raw") {
+		object@imputation[[nn.type]]@imp.velo <- dge.raw
+		return(object)
+	}
+
+	if (ncol(dge.raw) < length(levels(nb.id.fac))) {
+		stop("Number of index beads larger than the number of beads in DGE.")
+	}
+
+	dge.raw.imputed <- imputeNiche.dgCMatrix(
+		dge.raw, 
+		nb.id.fac,
+		nb.dist.fac,
+		weights = weights
+	)
+
+	# weighted sum of scale factors as well
+	scale.fac.imp <- imputeNiche.dgCMatrix(
+		Matrix::Matrix(object@parameters[["velo.lib.size"]], sparse = T, nrow = 1,
+						dimnames = list(NULL, names(object@parameters[["velo.lib.size"]]))), 
+		nb.id.fac,
+		nb.dist.fac,
+		weights = weights
+	)
+
+	scale.fac.imp <- as.numeric(scale.fac.imp)
+	names(scale.fac.imp) <- names(object@parameters[["velo.lib.size"]])
+	object@imputation[[nn.type]]@scale.fac.velo <- scale.fac.imp
+
+	res.density <- sum(dge.raw.imputed != 0)/length(dge.raw.imputed) # density 6.2%
+	cat(paste0("Density after imputation: ", res.density*100, "%\n"))
+
+	object@imputation[[nn.type]]@imp.velo <- dge.raw.imputed
+	object@imputation[[nn.type]]@log[["Density:"]] <- paste0(res.density*100, "%")
+
+	return(object)
+}
+
+
+
 #' Normalize the data using the specified method
 #' 
 #' @param object A Cytosignal object
@@ -535,13 +651,15 @@ normCounts.list <- function(
 #' @param object A Cytosignal object
 #' @param method The method to use for normalization
 #' @param slot.use The slot to use for normalization
+#' @param velo Whether to normalize the velocity data instead
 #' 
 #' @return A Cytosignal object
 #' @export
 normCounts.CytoSignal <- function(
 	object, 
 	method = c("default", "cpm"),
-	slot.use = NULL
+	slot.use = NULL,
+	velo = FALSE
 ){
 	if (is.null(slot.use)){
 		slot.use <- object@imputation[["default"]]
@@ -553,14 +671,23 @@ normCounts.CytoSignal <- function(
 
 	message(paste0("Normalizing on Imp slot: ", slot.use, "..."))
 
-	mat <- object@imputation[[slot.use]]@imp.data
-	scale.fac <- object@imputation[[slot.use]]@scale.fac
+	if (velo){
+		mat <- object@imputation[[slot.use]]@imp.velo
+		scale.fac <- object@imputation[[slot.use]]@scale.fac.velo
+	} else {
+		mat <- object@imputation[[slot.use]]@imp.data
+		scale.fac <- object@imputation[[slot.use]]@scale.fac
+	}
 
 	mat <- normCounts.list(
 		list(mat = mat, scale.fac = scale.fac), 
 		method = method)
 
-	object@imputation[[slot.use]]@imp.data <- mat
+	if (velo){
+		object@imputation[[slot.use]]@imp.velo <- mat
+	} else {
+		object@imputation[[slot.use]]@imp.data <- mat
+	}
 
 	# if (!all.equal(names(slot(object, slot.use)), names(object@imp.norm.data))){
 	# 	warning("Names of imp.data and imp.norm.data not equal!")
@@ -571,174 +698,6 @@ normCounts.CytoSignal <- function(
 }
 
 
-# normCounts_list <- function(object, 
-# 	slot.use = NULL, 
-# 	method = c("default", "cpm")
-# ){
-# 	if (is.null(slot.use)){
-# 		slot.use <- object@imputation[["default"]]
-# 	}
-
-# 	if (!slot.use %in% slotNames(object)) {
-# 		stop("Data not found.")
-# 	}
-
-# 	mat_list_norm <- lapply(object@imputation[[slot.use]], function(mat) {
-# 		if (method == "default") {
-# 			mat@x = mat@x / rep.int(Matrix::colSums(mat), diff(mat@p))
-# 			mat@x = lop1p(mat@x * 1e4)
-# 		} else if (method == "cpm") {
-# 			mat@x = mat@x / rep.int(Matrix::colSums(mat), diff(mat@p))
-# 			mat@x = lop1p(mat@x * 1e6)
-# 		} else {
-# 			stop("Method not found.")
-# 		}
-		
-# 		return(mat)
-# 	})
-	
-# 	object@imp.norm.data <- mat_list_norm
-	
-# 	if (!all.equal(names(slot(object, slot.use)), names(object@imp.norm.data))) {
-# 		warning("Names of imp.data and imp.norm.data not equal!")
-# 	}
-	
-# 	return(object)
-# }
-
-
-# set a new class union containing dgCMatrix and matrix
-setClassUnion(
-	"matrix_like",
-	c("dgCMatrix", "matrix")
-)
-
-#' Subset the imputed data (intr.data) by the specified intr index
-#' 
-#' @param object A Cytosignal object
-#' @param intr.index The intr index to use for subsetting
-#' @param slot.use The slot to use for subsetting
-#' 
-#' @return A Cytosignal object
-#' @export
-#' 
-changeUniprot <- function(
-  object,
-  ...
-) {
-  UseMethod(generic = 'changeUniprot', object = object)
-}
-
-#' Sub function for changeUniprot, input a sparse matrix
-#' 
-#' @param dge.raw A sparse matrix
-#' @param gene_to_uniprot A data frame with gene name and uniprot id
-#' @param verbose Whether to print out the progress
-#' 
-#' @return A sparse matrix
-#' @export
-#' 
-changeUniprot.matrix_like <- function(
-    dge.raw, 
-    gene_to_uniprot,
-    # mode = "unpt",
-    verbose = T
-){
-    # uppercase the gene to match database
-    rownames(dge.raw) = toupper(rownames(dge.raw))
-
-	unique.check = sapply(rownames(dge.raw), function(x){
-		length(unique(gene_to_uniprot$uniprot[which(gene_to_uniprot$gene_name == x)]))
-	})
-
-	if (max(unique.check) > 1){ # remove duplicated genes if nany
-		dup.num = sum(unique.check > 1)
-		stop(paste0("Must remove duplicated Uniprot items: ", dup.num, " genes."))
-	}
-    
-    uniprot.genes = rownames(dge.raw)[rownames(dge.raw) %in% gene_to_uniprot$gene_name] # 738/22683
-
-    if (verbose){
-        cat(paste0("Number of genes in the database: ", length(uniprot.genes), "\n"))
-    }
-
-    uniprot.names = sapply(uniprot.genes, function(x){
-        unique(gene_to_uniprot$uniprot[which(gene_to_uniprot$gene_name == x)])
-    })
-    
-    dge.raw = dge.raw[uniprot.genes, ]
-
-    # rownames(dge.raw) = unname(lowwords(rownames(dge.raw)))
-	rownames(dge.raw) = unname(uniprot.names)
-
-    # if (mode == "unpt"){
-    #     rownames(dge.raw) = unname(uniprot.names)
-    # }
-
-    if (anyDuplicated(rownames(dge.raw)) != 0){
-		del.num <- sum(duplicated(rownames(dge.raw)))
-		message(paste0("Removing ", del.num, " duplicated unpt genes."))
-
-        if (del.num > 10){
-            stop("More than 10 duplicated Uniprot items detected.")
-        }
-
-        dup.index = which(duplicated(rownames(dge.raw)) == T)
-        dge.raw = dge.raw[-dup.index, ]
-		uniprot.names = uniprot.names[-dup.index]
-		uniprot.genes = uniprot.genes[-dup.index]
-    }
-    
-	names(uniprot.names) <- lowwords(uniprot.genes)
-
-    return(list(
-		dge.raw,
-		uniprot.names
-	))
-}
-
-# #' Sub function for changeUniprot, input a cytosignal object
-# #' 
-# #' @param object A Cytosignal object
-# #' @param slot.use The slot to use for subsetting
-# #' @param verbose Whether to print out the progress
-# #' 
-# #' @return A Cytosignal object
-# #' @export
-# changeUniprot.CytoSignal <- function(
-# 	object,
-# 	slot.use = NULL, 
-# 	# mode = "unpt",
-# 	verbose = T
-# ){
-# 	if (is.null(slot.use)){
-# 		slot.use <- object@imputation[["default"]]
-# 	}
-
-# 	if (!slot.use %in% names(object@imputation)){
-# 		stop("Data not found.")
-# 	}
-
-# 	gene_to_uniprot <- object@intr.valid[["gene_to_uniprot"]]
-# 	dge.raw <- object@counts
-# 	# check duplicate items exists in database
-
-# 	unpt.list <- changeUniprot.matrix_like(dge.raw, gene_to_uniprot, verbose = verbose)
-	
-# 	object@imputation[[slot.use]]@intr.data <- unpt.list[[1]]
-# 	object@intr.valid[["symbols"]][[slot.use]] <- unpt.list[[2]]
-
-# 	scale.fac <- object@imputation[[slot.use]]@scale.fac$raw
-# 	scale.fac <- scale.fac[colnames(unpt.list[[1]])]
-# 	object@imputation[[slot.use]]@scale.fac$raw <- scale.fac
-
-# 	# if (!all.equal(names(slot(object, slot.use)), names(object@intr.data))) {
-# 	# 	warning("Names of imp.norm.data and intr.data not equal!")
-# 	# }
-
-# 	return(object)
-# }
-
 
 #' Compute the LR score for specific ligand-receptor imputation obj pairs
 #' 
@@ -878,14 +837,16 @@ inferScoreLR.CytoSignal <- function(
 #' 
 #' @param object A Cytosignal object
 #' @param nn.type The imputation method to use
-#' @param perm.size Size of the permutation test
+#' @param perm.size Size of the permutation test, by defualt NULL
+#' @param times Number of times to permute the whole dataset, by default 10
 #' 
 #' @return A Cytosignal object
 #' @export
 permuteImpLR <- function(
 	object,
 	nn.type = NULL,
-	perm.size = 100000
+	perm.size = NULL,
+	times = 10
 ){
 	if (is.null(nn.type)){
 		nn.type <- object@imputation[["default"]]
@@ -907,26 +868,43 @@ permuteImpLR <- function(
 
 	dimnames.list <- dimnames(object@imputation[[nn.type]]@imp.data)
 	cells.loc <- object@cells.loc
-	dge.raw <- object@raw.counts
+	dge.raw <- object@counts
 
-	# decide how many times to permute according to the size of the data
-	if (!is.numeric(perm.size)){
-		stop("perm.size must be a numeric value.")
-	}
+	if (!is.null(times)) {
+		if (!is.numeric(times)){
+			stop("times must be a numeric value.")
+		}
+
+		times <- ceiling(times)
 
-	if (perm.size < ncol(dge.raw)){
-		message("Permutation size too small, using ", ncol(dge.raw), " instead.")
-		perm.size <- ncol(dge.raw)
+		if (times < 1) {
+			stop("times must be larger than 1.")
+		}
+
+	} else {
+		# decide how many times to permute according to the size of the data
+		if (is.null(perm.size)){
+			message("Permutation size not specified, using 100000 instead.")
+			perm.size <- 100000
+		}
+
+		if (!is.numeric(perm.size)){
+			stop("perm.size must be a numeric value.")
+		}
+
+		if (perm.size < ncol(dge.raw)){
+			message("Permutation size too small, using ", ncol(dge.raw), " instead.")
+			perm.size <- ncol(dge.raw)
+		}
+		
+		times <- ceiling(perm.size / ncol(dge.raw))
 	}
-	
-	times <- ceiling(perm.size / ncol(dge.raw))
 
 	message("Permuting whole dataset ", times, " times...")
 
 	# permute cell locations
 	null.cells.loc.list = lapply(1:times, function(i){
-		perm.index = sample(seq_len(nrow(cells.loc)), nrow(cells.loc), replace = F)
-		null.cells.loc = cells.loc[perm.index,]
+		null.cells.loc = cells.loc[sample(nrow(cells.loc)), ]
 		rownames(null.cells.loc) = rownames(cells.loc)
 		return(null.cells.loc)
 	})
@@ -934,11 +912,15 @@ permuteImpLR <- function(
 	# Null imputations
 	null.dge.list <- lapply(1:times, function(i){
 		## MUST shuffule raw dge!!
-		null.dge.raw.all <- dge.raw[, sample(ncol(dge.raw))]
-		colnames(null.dge.raw.all) <- colnames(dge.raw)
+		perm.index <- sample(ncol(dge.raw))
+		null.dge.raw <- dge.raw[, perm.index]
+		scale.fac <- object@parameters[["lib.size"]][perm.index]
+		
+		colnames(null.dge.raw) <- colnames(dge.raw)
+		names(scale.fac) <- colnames(dge.raw)
 
-		scale.fac <- Matrix::Matrix(Matrix::colSums(null.dge.raw.all), nrow = 1, byrow = T, sparse = T) # computing scale factor
-		null.dge.raw <- changeUniprot.matrix_like(null.dge.raw.all, object@intr.valid[["gene_to_uniprot"]])[[1]]
+		scale.fac <- Matrix::Matrix(scale.fac, nrow = 1, byrow = T, sparse = T) # computing scale factor
+		# null.dge.raw <- changeUniprot.matrix_like(null.dge.raw.all, object@intr.valid[["gene_to_uniprot"]])[[1]]
 
 		if (use.gau) {
 			param.eps <- object@parameters$r.diffuse.scale
@@ -976,68 +958,6 @@ permuteImpLR <- function(
 
 }
 
-#' Permute LR score for specific ligand-receptor imputation obj pairs
-#' 
-#' This function is a follow-up function of inferScoreLR. It computes the NUL LRscores 
-#' using the NULL imputation results and stores the results in the LR score object.
-#' The null distribution of the LR scores can be used to test the significance of the LR scores.
-#' 
-#' @param object A Cytosignal object
-#' @param slot.use The imputation method to use
-#' 
-#' @return A Cytosignal object
-#' @export
-permuteNicheScoreLR <- function(
-	object,
-	slot.use = NULL
-){
-	if (is.null(slot.use)){
-		slot.use <- object@lrscore[["default"]]
-	}
-
-	message("Permuting scores on Score slot: ", slot.use, "...")
-
-	# score.obj <- object@lrscore[[slot.use]]
-	lig.slot <- object@lrscore[[slot.use]]@lig.slot
-	recep.slot <- object@lrscore[[slot.use]]@recep.slot
-	cells.loc <- object@cells.loc
-	intr.valid <- object@lrscore[[slot.use]]@intr.list
-
-	null.dge.eps.unpt <- object@imputation[[lig.slot]]@imp.data.null
-	null.dge.dt.unpt <- object@imputation[[recep.slot]]@imp.data.null
-
-    perm.index = sample(seq_len(nrow(cells.loc)), nrow(cells.loc), replace = F)
-    null.cells.loc = cells.loc[perm.index,]
-    rownames(null.cells.loc) = rownames(cells.loc)
-    null.nb.fac = findNNDT.matrix(null.cells.loc); gc() # Mean num of neighbors: 44, median: 36
-
-    null.lrscore.mtx = graphNicheLR.dgCMatrix(null.dge.eps.unpt, null.dge.dt.unpt, null.nb.fac[["id"]],
-						intr.valid[["ligands"]], intr.valid[["receptors"]]); gc()
-
-	# null.lrscore.mtx <- null.dge.eps.unpt * null.dge.dt.unpt
- 
-    if (sum(Matrix::colSums(null.lrscore.mtx) == 0) != 0){
-        message("A total of ", sum(Matrix::colSums(null.lrscore.mtx) == 0), " intr are empty in NULL scores.")
-        null.lrscore.mtx = null.lrscore.mtx[, !Matrix::colSums(null.lrscore.mtx) == 0]
-    }
-
-	intr.both <- intersect(colnames(null.lrscore.mtx), colnames(object@lrscore[[slot.use]]@score))
-
-	if (length(intr.both) != ncol(null.lrscore.mtx)) {
-		message("Removing ", ncol(null.lrscore.mtx) - length(intr.both), " more intr from NULL scores.")
-		null.lrscore.mtx = null.lrscore.mtx[, intr.both]
-	}
-
-	if (length(intr.both) != ncol(object@lrscore[[slot.use]]@score)) {
-		message("Removing ", ncol(object@lrscore[[slot.use]]@score) - length(intr.both), " corresponding intr from REAL scores.")
-		object@lrscore[[slot.use]]@score = object@lrscore[[slot.use]]@score[, intr.both]
-	}
-
-	object@lrscore[[slot.use]]@score.null <- null.lrscore.mtx
-
-    return(object)
-}
-
 
 #' Permute LR score for specific ligand-receptor imputation obj pairs
 #' 
@@ -1251,7 +1171,7 @@ inferSignif.CytoSignal <- function(
 	use.intr.slot.name <- object@lrscore[[slot.use]]@intr.slot
 	use.intr.db <- object@intr.valid[[use.intr.slot.name]]
 
-    res.list = inferSignif(object@raw.counts, object@lrscore[[slot.use]]@score, 
+    res.list = inferSignif.matrix_like(object@counts, object@lrscore[[slot.use]]@score, 
 				object@lrscore[[slot.use]]@score.null, nb.fac,
 				use.intr.db, object@intr.valid[["gene_to_uniprot"]], p.value,
 				reads.thresh, sig.thresh)
@@ -1337,14 +1257,14 @@ runPears.std <- function(
 #' @return A Cytosignal object
 #' @export
 #' 
-veloNicheLR <- function(
+inferVeloLR <- function(
   object,
   ...
 ) {
-  UseMethod(generic = 'veloNicheLR', object = object)
+  UseMethod(generic = 'inferVeloLR', object = object)
 }
 
-#' Sub function for veloNicheLR, input a sparse matrix
+#' Sub function for inferVeloLR, input a sparse matrix
 #' 
 #' @param dge.lig A sparse matrix for ligand
 #' @param dge.recep A sparse matrix for receptor
@@ -1355,19 +1275,19 @@ veloNicheLR <- function(
 #' @return A sparse matrix
 #' @export
 #' 
-veloNicheLR.matrix_like <- function(
+inferVeloLR.matrix_like <- function(
 	dge.lig,
     dge.recep,
-	velo.mtx,
-    nb.id.fac,
+	dge.lig.velo,
+	dge.recep.velo,
 	lig.fac,
 	recep.fac
 ){
 
 	### cavaet: remember to convert the Uniprot ids to indices!
 	# convert nb fac
-	nb.id.fac = sort(nb.id.fac)
-	nb.index = facToIndex(nb.id.fac)
+	# nb.id.fac = sort(nb.id.fac)
+	# nb.index = facToIndex(nb.id.fac)
 
 	if (max(as.integer(names(lig.fac))) > nrow(dge.lig)){
 		stop("Intr index out of dge bounds.")
@@ -1377,16 +1297,16 @@ veloNicheLR.matrix_like <- function(
 	recep.index = facToIndex(recep.fac)
 
 	# compute velos
-	res.mtx = VeloinferScoreLR_cpp(
+	res.mtx = inferVeloLR_cpp(
 		unname(as.matrix(dge.lig)),
 		unname(as.matrix(dge.recep)),
-		unname(as.matrix(velo.mtx)),
+		unname(as.matrix(dge.lig.velo)),
+		unname(as.matrix(dge.recep.velo)),
 		lig.index[[1]], lig.index[[2]], 
-		recep.index[[1]], recep.index[[2]],
-		nb.index[[1]], nb.index[[2]]
+		recep.index[[1]], recep.index[[2]]
 	)
 
-	dimnames(res.mtx) = list(colnames(dge.lig)[as.integer(levels(nb.id.fac))], levels(lig.fac))
+	dimnames(res.mtx) = list(colnames(dge.lig), levels(lig.fac))
 	# dimnames(res.mtx) = list(colnames(dge.lig), levels(lig.fac))
 	# res.mtx = Matrix(res.mtx, sparse = T)
 
@@ -1395,7 +1315,7 @@ veloNicheLR.matrix_like <- function(
 
 
 
-#' Sub function for veloNicheLR, input a CytoSignal object
+#' Sub function for inferVeloLR, input a CytoSignal object
 #' 
 #' @param object A Cytosignal object
 #' @param lig.slot The ligand slot to use
@@ -1407,7 +1327,7 @@ veloNicheLR.matrix_like <- function(
 #' 
 #' @return A Cytosignal object
 #' @export
-veloNicheLR.CytoSignal <- function(
+inferVeloLR.CytoSignal <- function(
 	object,
 	lig.slot,
 	recep.slot,
@@ -1439,78 +1359,51 @@ veloNicheLR.CytoSignal <- function(
 	
 	object@lrvelo[["default"]] <- tag
 
-	# if (is.null(slot.use)) {
-	# 	slot.use <- object@lrscore[["default"]]
-	# }
-
-	# if (!slot.use %in% names(object@lrscore)) {
-	# 	stop("lrvelo slot not found.")
-	# }
-
-	# if (is.null(nn.use)) {
-	# 	nn.use <- recep.slot
-	# }
-	
-	# if (is.character(nn.use)) {
-	# 	if (!nn.use %in% names(object@imputation)){
-	# 		stop("Imputation slot not found.")
-	# 	}
-	# 	nb.id.fac <- object@imputation[[nn.use]]@nn.id
-	# } else if (is.factor(nn.use)) {
-	# 	if (length(nn.use) != ncol(object@imputation[[lig.slot]]@intr.data))
-	# 		stop("nn.use must have the same length as the number of cells.")
-	# 	nb.id.fac <- nn.use
-	# } else {
-	# 	stop("nn.use must be either a factor or a character.")
-	# }
-
-	dge.lig <- object@imputation[[lig.slot]]@intr.data
-	dge.recep <- object@imputation[[recep.slot]]@intr.data
+	dge.lig <- object@imputation[[lig.slot]]@imp.data
+	dge.recep <- object@imputation[[recep.slot]]@imp.data
+	dge.lig.velo <- object@imputation[[lig.slot]]@imp.velo
+	dge.recep.velo <- object@imputation[[recep.slot]]@imp.velo
 
-	if (!all.equal(dim(dge.lig), dim(dge.recep))){
-		stop("dge.lig and dge.recep must have the same dimension.")
+	# compare the dimnames of all four matrices
+	if (!all.equal(dimnames(dge.lig), dimnames(dge.recep))){
+		stop("dge.lig and dge.recep must have the same dimension names.")
 	}
 
-	if (!all.equal(object@intr.valid[["symbols"]][[lig.slot]], 
-					object@intr.valid[["symbols"]][[recep.slot]])){
-		stop("Unpt symbols generated from imputations not equal!")
+	if (!all.equal(dimnames(dge.lig.velo), dimnames(dge.recep.velo))){
+		stop("dge.lig.velo and dge.recep.velo must have the same dimension names.")
 	}
-	
-	# velo and dge should have the same cells and genes
-	use.cells <- intersect(colnames(dge.lig), colnames(object@velo[["velo.s"]]))
-	use.genes <- intersect(rownames(dge.lig), rownames(object@velo[["velo.s"]]))
 
-	# infer new neighbor index since the cells are filtered
-	new.cells.loc <- object@cells.loc[use.cells, ]
-	velo.nn.list <- findNNDT.matrix(new.cells.loc)
+	use.genes <- intersect(rownames(dge.lig), rownames(dge.lig.velo))
 
-	# dge and velo should have the same cells and genes
-	dge.lig <- dge.lig[use.genes, use.cells]
-	dge.recep <- dge.recep[use.genes, use.cells]
-	velo.s <- object@velo[["velo.s"]][use.genes, use.cells]
-	velo.u <- object@velo[["velo.u"]][use.genes, use.cells]
-
-	message("Number of velo cells: ", ncol(dge.lig), " / ", ncol(object@imputation[[lig.slot]]@intr.data))
-	message("Number of velo genes: ", nrow(dge.lig), " / ", nrow(object@imputation[[lig.slot]]@intr.data))
+	dge.lig <- dge.lig[use.genes, ]
+	dge.recep <- dge.recep[use.genes, ]
+	dge.lig.velo <- dge.lig.velo[use.genes, ]
+	dge.recep.velo <- dge.recep.velo[use.genes, ]
 
+	message("Number of velo genes: ", length(use.genes), " / ", nrow(object@imputation[[lig.slot]]@imp.data))
+	
 	intr.db.list <- checkIntr(use.genes, object@intr.valid[[intr.db.name]])
 
-	res.mtx <- veloNicheLR.matrix_like(dge.lig, dge.recep, 
-				(velo.s + velo.u), velo.nn.list[["id"]], 
+	res.mtx <- inferVeloLR.matrix_like(dge.lig, dge.recep, 
+				dge.lig.velo, dge.recep.velo,
 				intr.db.list[["ligands"]], intr.db.list[["receptors"]])
 
+	intr.hq <- names(object@lrscore[[tag]]@res.list[["result.hq.pear"]])
+	intr.hq <- intr.hq[intr.hq %in% colnames(res.mtx)]
+
 	lrvelo.obj <- new(
 		"lrVelo",
 		lig.slot = lig.slot,
 		recep.slot = recep.slot,
 		intr.slot = intr.db.name,
 		intr.list = intr.db.list,
-		dim.valid = list(
-			"cells" = use.cells,
-			"genes" = use.genes),
 		intr.velo = res.mtx,
-		nn.id = velo.nn.list[["id"]],
-		nn.dist = velo.nn.list[["dist"]],
+		velo.gene = list(
+			"velo.genes" = use.genes,
+			"velo.intr" = colnames(res.mtx),
+			"velo.intr.hq" = intr.hq),
+		# nn.id = velo.nn.list[["id"]],
+		# nn.dist = velo.nn.list[["dist"]],
 		log = list(
 			"Used slot" = c(lig.slot, recep.slot)
 		)
diff --git a/R/archive.r b/R/archive.r
index a067cc1..4f48575 100644
--- a/R/archive.r
+++ b/R/archive.r
@@ -1,3 +1,386 @@
+#' Permute Imputation Results of specific imputation method
+#' 
+#' This function is a follow-up function of imputeNiche, which permutes the default or user-defined
+#' imputation method and stored the results in the ImpData object.
+#' 
+#' @param object A Cytosignal object
+#' @param nn.type The imputation method to use
+#' @param perm.size Size of the permutation test, by defualt NULL
+#' @param times Number of times to permute the whole dataset, by default 10
+#' 
+#' @return A Cytosignal object
+#' @export
+permuteImpEdge <- function(
+	object,
+	nn.type = NULL,
+	perm.size = NULL,
+	times = 10
+){
+	if (is.null(nn.type)){
+		nn.type <- object@imputation[["default"]]
+	}
+
+	if (!nn.type %in% names(object@imputation)){
+		stop("Cannot find corresponding imputation method.")
+	}
+
+	use.gau <- grepl("GauEps", nn.type)
+	use.dt <- grepl("DT", nn.type)
+	use.raw <- grepl("Raw", nn.type)
+
+	# if (!use.gau & !use.dt){
+	# 	stop("Cannot find corresponding imputation method.")
+	# }
+
+	message("Permuting imputation data on method ", nn.type, "...")
+
+	dimnames.list <- dimnames(object@imputation[[nn.type]]@imp.data)
+	cells.loc <- object@cells.loc
+	dge.raw <- object@raw.counts
+
+	if (!is.null(times)) {
+		if (!is.numeric(times)){
+			stop("times must be a numeric value.")
+		}
+
+		times <- ceiling(times)
+
+		if (times < 1) {
+			stop("times must be larger than 1.")
+		}
+
+	} else {
+		# decide how many times to permute according to the size of the data
+		if (is.null(perm.size)){
+			message("Permutation size not specified, using 100000 instead.")
+			perm.size <- 100000
+		}
+
+		if (!is.numeric(perm.size)){
+			stop("perm.size must be a numeric value.")
+		}
+
+		if (perm.size < ncol(dge.raw)){
+			message("Permutation size too small, using ", ncol(dge.raw), " instead.")
+			perm.size <- ncol(dge.raw)
+		}
+		
+		times <- ceiling(perm.size / ncol(dge.raw))
+	}
+
+	message("Permuting whole dataset ", times, " times...")
+
+	# # permute cell locations
+	# null.cells.loc.list = lapply(1:times, function(i){
+	# 	null.cells.loc = cells.loc[sample(nrow(cells.loc)), ]
+	# 	rownames(null.cells.loc) = rownames(cells.loc)
+	# 	return(null.cells.loc)
+	# })
+
+	null.cells.loc = cells.loc[sample(nrow(cells.loc)), ]
+	rownames(null.cells.loc) = rownames(cells.loc)
+
+	if (use.gau) {
+		param.eps <- object@parameters$r.diffuse.scale
+		param.sigma <- object@parameters$sigma.scale
+
+		null.eps.nb <- findNNGauEB.matrix(null.cells.loc, eps = param.eps, 
+						sigma = param.sigma, weight.sum = 2); gc()
+
+		null.graph <- imputeNiche.dgCMatrix(dge.raw, null.eps.nb$id, null.eps.nb$dist, 
+							weights = "none", return.graph = T); gc()
+	} else if (use.dt) {
+		null.nb.fac <- findNNDT.matrix(null.cells.loc)
+		null.graph <- imputeNiche.dgCMatrix(dge.raw, null.nb.fac$id, null.nb.fac$dist, 
+							weights = "none", return.graph = T); gc()
+	} else {
+		stop("Cannot find corresponding imputation method.")
+	}
+
+	null.dge.list <- lapply(1:times, function(i){
+		## MUST shuffule raw dge!!
+		null.dge.raw.all <- dge.raw[, sample(ncol(dge.raw))]
+		colnames(null.dge.raw.all) <- colnames(dge.raw)
+
+		scale.fac <- Matrix::Matrix(Matrix::colSums(null.dge.raw.all), nrow = 1, byrow = T, sparse = T) # computing scale factor
+		null.dge.raw <- changeUniprot.matrix_like(null.dge.raw.all, object@intr.valid[["gene_to_uniprot"]])[[1]]
+
+		if (use.gau) {
+			param.eps <- object@parameters$r.diffuse.scale
+			param.sigma <- object@parameters$sigma.scale
+
+			# permute graph and impute
+			
+
+			null.dge.eps <- null.dge.raw %*% null.graph; gc()
+
+			null.dge.eps <- imputeNiche.dgCMatrix(null.dge.raw, null.eps.nb$id, null.eps.nb$dist, weights = "none"); gc()
+			scale.fac.imp <- imputeNiche.dgCMatrix(scale.fac, null.eps.nb$id, null.eps.nb$dist, weights = "none")
+		} else if (use.dt) {
+			null.nb.fac <- findNNDT.matrix(null.cells.loc.list[[i]])
+        	null.dge.eps <- imputeNiche.dgCMatrix(null.dge.raw, null.nb.fac$id, null.nb.fac$dist, weights = "none"); gc() # den: 31%
+			scale.fac.imp <- imputeNiche.dgCMatrix(scale.fac, null.nb.fac$id, null.nb.fac$dist, weights = "none")
+		} else if (use.raw) {
+			null.dge.eps <- null.dge.raw
+			scale.fac.imp <- scale.fac
+		} else {
+			stop("Cannot find corresponding imputation method.")
+		}
+		
+		null.dge.eps = normCounts.list(list(mat=null.dge.eps, scale.fac=as.numeric(scale.fac.imp)), "default"); gc()
+		rm(null.dge.raw); gc()
+
+		return(null.dge.eps)
+	})
+
+	
+	# # Null imputations
+	# null.dge.list <- lapply(1:times, function(i){
+	# 	## MUST shuffule raw dge!!
+	# 	null.dge.raw.all <- dge.raw[, sample(ncol(dge.raw))]
+	# 	colnames(null.dge.raw.all) <- colnames(dge.raw)
+
+	# 	scale.fac <- Matrix::Matrix(Matrix::colSums(null.dge.raw.all), nrow = 1, byrow = T, sparse = T) # computing scale factor
+	# 	null.dge.raw <- changeUniprot.matrix_like(null.dge.raw.all, object@intr.valid[["gene_to_uniprot"]])[[1]]
+
+	# 	if (use.gau) {
+	# 		param.eps <- object@parameters$r.diffuse.scale
+	# 		param.sigma <- object@parameters$sigma.scale
+
+	# 		null.eps.nb <- findNNGauEB.matrix(null.cells.loc.list[[i]], eps = param.eps, 
+	# 						sigma = param.sigma, weight.sum = 2); gc()
+
+	# 		null.dge.eps <- imputeNiche.dgCMatrix(null.dge.raw, null.eps.nb$id, null.eps.nb$dist, weights = "none"); gc()
+	# 		scale.fac.imp <- imputeNiche.dgCMatrix(scale.fac, null.eps.nb$id, null.eps.nb$dist, weights = "none")
+	# 	} else if (use.dt) {
+	# 		null.nb.fac <- findNNDT.matrix(null.cells.loc.list[[i]])
+    #     	null.dge.eps <- imputeNiche.dgCMatrix(null.dge.raw, null.nb.fac$id, null.nb.fac$dist, weights = "none"); gc() # den: 31%
+	# 		scale.fac.imp <- imputeNiche.dgCMatrix(scale.fac, null.nb.fac$id, null.nb.fac$dist, weights = "none")
+	# 	} else if (use.raw) {
+	# 		null.dge.eps <- null.dge.raw
+	# 		scale.fac.imp <- scale.fac
+	# 	} else {
+	# 		stop("Cannot find corresponding imputation method.")
+	# 	}
+		
+	# 	null.dge.eps = normCounts.list(list(mat=null.dge.eps, scale.fac=as.numeric(scale.fac.imp)), "default"); gc()
+	# 	rm(null.dge.raw); gc()
+
+	# 	return(null.dge.eps)
+	# })
+
+	null.dge.eps.unpt = meanMat_cpp(null.dge.list, nrow(null.dge.list[[1]]), ncol(null.dge.list[[1]]))
+	dimnames(null.dge.eps.unpt) <- dimnames.list
+	rm(null.dge.list); gc()
+
+    object@imputation[[nn.type]]@imp.data.null <- null.dge.eps.unpt
+
+	return(object)
+
+}
+
+
+
+#' Compute the LR velo for specific ligand-receptor imputation obj pairs
+#' 
+#' @param object A Cytosignal object
+#' @param lig.slot The ligand slot to use
+#' @param recep.slot The receptor slot to use
+#' @param intr.db.name The intr database name to use
+#' @param nn.use The neighbor index as niche
+#' 
+#' @return A Cytosignal object
+#' @export
+#' 
+veloNicheLR <- function(
+  object,
+  ...
+) {
+  UseMethod(generic = 'veloNicheLR', object = object)
+}
+
+#' Sub function for veloNicheLR, input a sparse matrix
+#' 
+#' @param dge.lig A sparse matrix for ligand
+#' @param dge.recep A sparse matrix for receptor
+#' @param nb.id.fac A factor of neighbor indices
+#' @param lig.fac A factor of ligand indices
+#' @param recep.fac A factor of receptor indices
+#' 
+#' @return A sparse matrix
+#' @export
+#' 
+veloNicheLR.matrix_like <- function(
+	dge.lig,
+    dge.recep,
+	velo.mtx,
+    nb.id.fac,
+	lig.fac,
+	recep.fac
+){
+
+	### cavaet: remember to convert the Uniprot ids to indices!
+	# convert nb fac
+	nb.id.fac = sort(nb.id.fac)
+	nb.index = facToIndex(nb.id.fac)
+
+	if (max(as.integer(names(lig.fac))) > nrow(dge.lig)){
+		stop("Intr index out of dge bounds.")
+	}
+
+	lig.index = facToIndex(lig.fac)
+	recep.index = facToIndex(recep.fac)
+
+	# compute velos
+	res.mtx = VeloinferScoreLR_cpp(
+		unname(as.matrix(dge.lig)),
+		unname(as.matrix(dge.recep)),
+		unname(as.matrix(velo.mtx)),
+		lig.index[[1]], lig.index[[2]], 
+		recep.index[[1]], recep.index[[2]],
+		nb.index[[1]], nb.index[[2]]
+	)
+
+	dimnames(res.mtx) = list(colnames(dge.lig)[as.integer(levels(nb.id.fac))], levels(lig.fac))
+	# dimnames(res.mtx) = list(colnames(dge.lig), levels(lig.fac))
+	# res.mtx = Matrix(res.mtx, sparse = T)
+
+	return(res.mtx)
+}
+
+
+
+#' Sub function for veloNicheLR, input a CytoSignal object
+#' 
+#' @param object A Cytosignal object
+#' @param lig.slot The ligand slot to use
+#' @param recep.slot The receptor slot to use
+#' @param intr.db.name The intr database name to use
+#' @param nn.use slot that the neighbor index should be taken from, by default is the same as
+#' 			the recep.slot. For example, if velo.obj = GauEps-DT, then nn.use = "DT".
+#' 			nn.use could also be a user-defind factor.
+#' 
+#' @return A Cytosignal object
+#' @export
+veloNicheLR.CytoSignal <- function(
+	object,
+	lig.slot,
+	recep.slot,
+	intr.db.name,
+	tag = NULL
+){
+	if (!lig.slot %in% names(object@imputation)){
+		stop("Ligand slot not found.")
+	}
+
+	if (!recep.slot %in% names(object@imputation)){
+		stop("Receptor slot not found.")
+	}
+
+	message("Computing velos using ", intr.db.name, " database.")
+	message("Ligand: ", lig.slot, ", Receptor: ", recep.slot, ".")
+
+	if (!intr.db.name %in% c("diff_dep", "cont_dep")) {
+		stop("intr.db.name must be either 'diff_dep' or 'cont_dep'.")
+	}
+
+	if (is.null(tag)) {
+		tag <- paste0(lig.slot, "-", recep.slot)
+	}
+
+	if (tag %in% names(object@lrvelo)) {
+		stop("Tag already exists.")
+	}
+	
+	object@lrvelo[["default"]] <- tag
+
+	# if (is.null(slot.use)) {
+	# 	slot.use <- object@lrscore[["default"]]
+	# }
+
+	# if (!slot.use %in% names(object@lrscore)) {
+	# 	stop("lrvelo slot not found.")
+	# }
+
+	# if (is.null(nn.use)) {
+	# 	nn.use <- recep.slot
+	# }
+	
+	# if (is.character(nn.use)) {
+	# 	if (!nn.use %in% names(object@imputation)){
+	# 		stop("Imputation slot not found.")
+	# 	}
+	# 	nb.id.fac <- object@imputation[[nn.use]]@nn.id
+	# } else if (is.factor(nn.use)) {
+	# 	if (length(nn.use) != ncol(object@imputation[[lig.slot]]@intr.data))
+	# 		stop("nn.use must have the same length as the number of cells.")
+	# 	nb.id.fac <- nn.use
+	# } else {
+	# 	stop("nn.use must be either a factor or a character.")
+	# }
+
+	dge.lig <- object@imputation[[lig.slot]]@intr.data
+	dge.recep <- object@imputation[[recep.slot]]@intr.data
+
+	if (!all.equal(dim(dge.lig), dim(dge.recep))){
+		stop("dge.lig and dge.recep must have the same dimension.")
+	}
+
+	if (!all.equal(object@intr.valid[["symbols"]][[lig.slot]], 
+					object@intr.valid[["symbols"]][[recep.slot]])){
+		stop("Unpt symbols generated from imputations not equal!")
+	}
+	
+	# velo and dge should have the same cells and genes
+	use.cells <- intersect(colnames(dge.lig), colnames(object@velo[["velo.s"]]))
+	use.genes <- intersect(rownames(dge.lig), rownames(object@velo[["velo.s"]]))
+
+	# infer new neighbor index since the cells are filtered
+	new.cells.loc <- object@cells.loc[use.cells, ]
+	velo.nn.list <- findNNDT.matrix(new.cells.loc)
+
+	# dge and velo should have the same cells and genes
+	dge.lig <- dge.lig[use.genes, use.cells]
+	dge.recep <- dge.recep[use.genes, use.cells]
+	velo.s <- object@velo[["velo.s"]][use.genes, use.cells]
+	velo.u <- object@velo[["velo.u"]][use.genes, use.cells]
+
+	message("Number of velo cells: ", ncol(dge.lig), " / ", ncol(object@imputation[[lig.slot]]@intr.data))
+	message("Number of velo genes: ", nrow(dge.lig), " / ", nrow(object@imputation[[lig.slot]]@intr.data))
+
+	intr.db.list <- checkIntr(use.genes, object@intr.valid[[intr.db.name]])
+
+	res.mtx <- veloNicheLR.matrix_like(dge.lig, dge.recep, 
+				(velo.s + velo.u), velo.nn.list[["id"]], 
+				intr.db.list[["ligands"]], intr.db.list[["receptors"]])
+
+	lrvelo.obj <- new(
+		"lrVelo",
+		lig.slot = lig.slot,
+		recep.slot = recep.slot,
+		intr.slot = intr.db.name,
+		intr.list = intr.db.list,
+		dim.valid = list(
+			"cells" = use.cells,
+			"genes" = use.genes),
+		intr.velo = res.mtx,
+		nn.id = velo.nn.list[["id"]],
+		nn.dist = velo.nn.list[["dist"]],
+		log = list(
+			"Used slot" = c(lig.slot, recep.slot)
+		)
+	)
+
+	object@lrvelo[[tag]] <- lrvelo.obj
+
+	return(object)
+}
+
+
+
+
+
+
 # #' Permute Imputation Results of specific imputation method
 # #' 
 # #' This function is a follow-up function of imputeNiche, which permutes the default or user-defined
@@ -106,6 +489,69 @@
 # }
 
 
+#' Permute LR score for specific ligand-receptor imputation obj pairs
+#' 
+#' This function is a follow-up function of inferScoreLR. It computes the NUL LRscores 
+#' using the NULL imputation results and stores the results in the LR score object.
+#' The null distribution of the LR scores can be used to test the significance of the LR scores.
+#' 
+#' @param object A Cytosignal object
+#' @param slot.use The imputation method to use
+#' 
+#' @return A Cytosignal object
+#' @export
+permuteNicheScoreLR <- function(
+	object,
+	slot.use = NULL
+){
+	if (is.null(slot.use)){
+		slot.use <- object@lrscore[["default"]]
+	}
+
+	message("Permuting scores on Score slot: ", slot.use, "...")
+
+	# score.obj <- object@lrscore[[slot.use]]
+	lig.slot <- object@lrscore[[slot.use]]@lig.slot
+	recep.slot <- object@lrscore[[slot.use]]@recep.slot
+	cells.loc <- object@cells.loc
+	intr.valid <- object@lrscore[[slot.use]]@intr.list
+
+	null.dge.eps.unpt <- object@imputation[[lig.slot]]@imp.data.null
+	null.dge.dt.unpt <- object@imputation[[recep.slot]]@imp.data.null
+
+    perm.index = sample(seq_len(nrow(cells.loc)), nrow(cells.loc), replace = F)
+    null.cells.loc = cells.loc[perm.index,]
+    rownames(null.cells.loc) = rownames(cells.loc)
+    null.nb.fac = findNNDT.matrix(null.cells.loc); gc() # Mean num of neighbors: 44, median: 36
+
+    null.lrscore.mtx = graphNicheLR.dgCMatrix(null.dge.eps.unpt, null.dge.dt.unpt, null.nb.fac[["id"]],
+						intr.valid[["ligands"]], intr.valid[["receptors"]]); gc()
+
+	# null.lrscore.mtx <- null.dge.eps.unpt * null.dge.dt.unpt
+ 
+    if (sum(Matrix::colSums(null.lrscore.mtx) == 0) != 0){
+        message("A total of ", sum(Matrix::colSums(null.lrscore.mtx) == 0), " intr are empty in NULL scores.")
+        null.lrscore.mtx = null.lrscore.mtx[, !Matrix::colSums(null.lrscore.mtx) == 0]
+    }
+
+	intr.both <- intersect(colnames(null.lrscore.mtx), colnames(object@lrscore[[slot.use]]@score))
+
+	if (length(intr.both) != ncol(null.lrscore.mtx)) {
+		message("Removing ", ncol(null.lrscore.mtx) - length(intr.both), " more intr from NULL scores.")
+		null.lrscore.mtx = null.lrscore.mtx[, intr.both]
+	}
+
+	if (length(intr.both) != ncol(object@lrscore[[slot.use]]@score)) {
+		message("Removing ", ncol(object@lrscore[[slot.use]]@score) - length(intr.both), " corresponding intr from REAL scores.")
+		object@lrscore[[slot.use]]@score = object@lrscore[[slot.use]]@score[, intr.both]
+	}
+
+	object@lrscore[[slot.use]]@score.null <- null.lrscore.mtx
+
+    return(object)
+}
+
+
 
 #' Permute Imputation Results of specific imputation method
 #' 
@@ -705,4 +1151,48 @@ graphNicheLR.CytoSignal <- function(
 	object@lrscore[[tag]] <- score.obj
 
 	return(object)
-}
\ No newline at end of file
+}
+
+
+
+# #' Sub function for changeUniprot, input a cytosignal object
+# #' 
+# #' @param object A Cytosignal object
+# #' @param slot.use The slot to use for subsetting
+# #' @param verbose Whether to print out the progress
+# #' 
+# #' @return A Cytosignal object
+# #' @export
+# changeUniprot.CytoSignal <- function(
+# 	object,
+# 	slot.use = NULL, 
+# 	# mode = "unpt",
+# 	verbose = T
+# ){
+# 	if (is.null(slot.use)){
+# 		slot.use <- object@imputation[["default"]]
+# 	}
+
+# 	if (!slot.use %in% names(object@imputation)){
+# 		stop("Data not found.")
+# 	}
+
+# 	gene_to_uniprot <- object@intr.valid[["gene_to_uniprot"]]
+# 	dge.raw <- object@counts
+# 	# check duplicate items exists in database
+
+# 	unpt.list <- changeUniprot.matrix_like(dge.raw, gene_to_uniprot, verbose = verbose)
+	
+# 	object@imputation[[slot.use]]@intr.data <- unpt.list[[1]]
+# 	object@intr.valid[["symbols"]][[slot.use]] <- unpt.list[[2]]
+
+# 	scale.fac <- object@imputation[[slot.use]]@scale.fac$raw
+# 	scale.fac <- scale.fac[colnames(unpt.list[[1]])]
+# 	object@imputation[[slot.use]]@scale.fac$raw <- scale.fac
+
+# 	# if (!all.equal(names(slot(object, slot.use)), names(object@intr.data))) {
+# 	# 	warning("Names of imp.norm.data and intr.data not equal!")
+# 	# }
+
+# 	return(object)
+# }
\ No newline at end of file
diff --git a/R/objects.r b/R/objects.r
index f8ec7cf..efe4c90 100644
--- a/R/objects.r
+++ b/R/objects.r
@@ -1,5 +1,6 @@
-#sourceCpp("/home/alan/Documents/Welch_lab/pj19-cell_cell_int/cytosignal_v1_02142023/cytosignal/src/mat_exp.cpp")
-#sourceCpp("/home/alan/Documents/Welch_lab/pj19-cell_cell_int/cytosignal_v1_02142023/cytosignal/src/utils_velo.cpp")
+#' @importFrom Matrix colSums rowSums t
+NULL
+
 
 #' The CytoSignal Class
 #'
@@ -46,10 +47,17 @@ CytoSignal <- setClass(
   )
 )
 
+# set a new class union containing dgCMatrix and matrix
+setClassUnion(
+	"matrix_like",
+	c("dgCMatrix", "matrix")
+)
+
+
 # set a new class union containing dgCMatrix and matrix
 setClassUnion(
 	"imp_null_class",
-	c("dgCMatrix", "list")
+	c("dgCMatrix", "matrix", "list")
 )
 
 #' The ImpData Class
@@ -78,20 +86,21 @@ ImpData <- setClass(
 	Class = "ImpData",
 	slots = c(
 		method = "character",
-		imp.data = "dgCMatrix",
+		imp.data = "matrix_like",
 		# imp.norm.data = "dgCMatrix",
 		# intr.data = "dgCMatrix",
 		imp.data.null = "imp_null_class",
+		imp.velo = "matrix_like",
 		nn.id = "factor",
 		nn.dist = "factor",
 		scale.fac = "numeric",
+		scale.fac.velo = "numeric",
 		log = "list"
 	)
 )
 
 
 
-
 #' The lrScores Class
 #' 
 #' The lrScores object is created from one ST dataset. User could choose two imputation methods to calculate
@@ -130,7 +139,6 @@ lrScores <- setClass(
 )
 
 
-
 #' The lrvelo Class
 #' 
 #' The lrvelo object is created from one ST dataset. User could choose two imputation methods to calculate
@@ -162,10 +170,10 @@ lrVelo <- setClass(
 		recep.slot = "character",
 		intr.slot = "character",
 		intr.list = "list",
-		dim.valid = "list",
 		intr.velo = "matrix_like",
-		nn.id = "factor",
-		nn.dist = "factor",
+		velo.gene = "list",
+		# nn.id = "factor",
+		# nn.dist = "factor",
 		log = "list"
 	)
 )
@@ -385,7 +393,7 @@ createCytoSignal <- function(
 		log[["Dataset"]] <- name
 	}
 	if (is.null(version)) {
-		version <- "v1.0.0"
+		version <- packageVersion("cytosignal")
 	}
 
 	# check the class of raw.data
@@ -516,6 +524,91 @@ removeLowQuality <- function(object, counts.thresh = 300, gene.thresh = 50) {
 
 
 
+#' Subset the imputed data (intr.data) by the specified intr index
+#' 
+#' @param object A Cytosignal object
+#' @param intr.index The intr index to use for subsetting
+#' @param slot.use The slot to use for subsetting
+#' 
+#' @return A Cytosignal object
+#' @export
+#' 
+changeUniprot <- function(
+  object,
+  ...
+) {
+  UseMethod(generic = 'changeUniprot', object = object)
+}
+
+#' Sub function for changeUniprot, input a sparse matrix
+#' 
+#' @param dge.raw A sparse matrix
+#' @param gene_to_uniprot A data frame with gene name and uniprot id
+#' @param verbose Whether to print out the progress
+#' 
+#' @return A sparse matrix
+#' @export
+#' 
+changeUniprot.matrix_like <- function(
+    dge.raw, 
+    gene_to_uniprot,
+    # mode = "unpt",
+    verbose = T
+){
+    # uppercase the gene to match database
+    rownames(dge.raw) = toupper(rownames(dge.raw))
+
+	unique.check = sapply(rownames(dge.raw), function(x){
+		length(unique(gene_to_uniprot$uniprot[which(gene_to_uniprot$gene_name == x)]))
+	})
+
+	if (max(unique.check) > 1){ # remove duplicated genes if nany
+		dup.num = sum(unique.check > 1)
+		stop(paste0("Must remove duplicated Uniprot items: ", dup.num, " genes."))
+	}
+    
+    uniprot.genes = rownames(dge.raw)[rownames(dge.raw) %in% gene_to_uniprot$gene_name] # 738/22683
+
+    if (verbose){
+        cat(paste0("Number of genes in the database: ", length(uniprot.genes), "\n"))
+    }
+
+    uniprot.names = sapply(uniprot.genes, function(x){
+        unique(gene_to_uniprot$uniprot[which(gene_to_uniprot$gene_name == x)])
+    })
+    
+    dge.raw = dge.raw[uniprot.genes, ]
+
+    # rownames(dge.raw) = unname(lowwords(rownames(dge.raw)))
+	rownames(dge.raw) = unname(uniprot.names)
+
+    # if (mode == "unpt"){
+    #     rownames(dge.raw) = unname(uniprot.names)
+    # }
+
+    if (anyDuplicated(rownames(dge.raw)) != 0){
+		del.num <- sum(duplicated(rownames(dge.raw)))
+		message(paste0("Removing ", del.num, " duplicated unpt genes."))
+
+        if (del.num > 10){
+            stop("More than 10 duplicated Uniprot items detected.")
+        }
+
+        dup.index = which(duplicated(rownames(dge.raw)) == T)
+        dge.raw = dge.raw[-dup.index, ]
+		uniprot.names = uniprot.names[-dup.index]
+		uniprot.genes = uniprot.genes[-dup.index]
+    }
+    
+	names(uniprot.names) <- lowwords(uniprot.genes)
+
+    return(list(
+		dge.raw,
+		uniprot.names
+	))
+}
+
+
 #' Sub function for changeUniprot, input a cytosignal object
 #' 
 #' @param object A Cytosignal object
@@ -550,13 +643,18 @@ changeUniprot.CytoSignal <- function(
 #' @return a CytoSignal object
 #' 
 #' @export
-purgeBeforeSave <- function(object) {
+purgeBeforeSave <- function(object, purge.raw = T, purge.null = F) {
 	# get slot names
 	imp.names = setdiff(names(object@imputation), "default")
 
-	for (i in 1:length(imp.names)) {
-		object@imputation[[imp.names[i]]]@imp.data <- new("dgCMatrix")
-		object@imputation[[imp.names[i]]]@imp.norm.data <- new("dgCMatrix")
+	if (purge.null) {
+		for (i in 1:length(imp.names)) {
+			object@imputation[[imp.names[i]]]@imp.data.null <- new("dgCMatrix")
+		}
+	}
+
+	if (purge.raw) {
+		object@raw.counts <- new("dgCMatrix"); gc()
 	}
 
 	return(object)
@@ -615,12 +713,18 @@ addVelo <- function(
 	velo.s.intr = changeUniprot.matrix_like(velo.s, object@intr.valid[["gene_to_uniprot"]])
 	velo.u.intr = changeUniprot.matrix_like(velo.u, object@intr.valid[["gene_to_uniprot"]])
 
+	velo.s.intr[[1]] <- Matrix::Matrix(increase_columns(colnames(object@counts), velo.s.intr[[1]]), sparse = T)
+	velo.u.intr[[1]] <- Matrix::Matrix(increase_columns(colnames(object@counts), velo.u.intr[[1]]), sparse = T)
+
+	rownames(velo.s.intr[[1]]) <- rownames(velo.u.intr[[1]]) <- rownames(object@counts)
+
 	object@velo <- list(
 		"velo.s" = velo.s.intr[[1]],
 		"velo.u" = velo.u.intr[[1]]
 	)
 
 	object@intr.valid[["symbols"]][["velo"]] = velo.s.intr[[2]]
+	object@parameters[["velo.lib.size"]] <- Matrix::colSums(velo.s.intr[[1]]) + Matrix::colSums(velo.u.intr[[1]])
 
 	return(object)
 }
diff --git a/R/plotting.r b/R/plotting.r
index 1aa4077..339733a 100644
--- a/R/plotting.r
+++ b/R/plotting.r
@@ -338,12 +338,12 @@ plotVelo <- function(
     cat("Plotting the results...\n")
 
     plot.list <- lapply(num.plot, function(i){
-        cat("No.", i, ", ", sep = "")
         intr.rank <- velo.intr.index[i]
+        cat("No.", intr.rank, ", ", sep = "")
         use.intr = names(use.res.list)[intr.rank]
         
         # get value for the ranked #1 intr
-        plot.df = as.data.frame(object@cells.loc[velo.obj@dim.valid$cells, ])
+        plot.df = as.data.frame(object@cells.loc)
         plot.df$velo = velo.obj@intr.velo[rownames(plot.df), use.intr]
 
         pt.df = plot.df[, c(1,2)]
@@ -396,10 +396,12 @@ plotVelo <- function(
 
         intr.name = getIntrNames(use.intr, object@intr.valid$intr.index)
 
+        sig.len <- lengths(use.res.list)[intr.rank]
+
         if (plot.fmt == "png") {
-            png(paste0(plot_dir, "/", "Rank_", intr.rank, "_", intr.name, "_3d.png"), width = width, height = height, units = "in", res = 400)
+            png(paste0(plot_dir, "/", "Rank_", intr.rank, "_", intr.name, "-n_", sig.len, "-3D.png"), width = width, height = height, units = "in", res = 400)
         } else if (plot.fmt == "pdf") {
-            pdf(paste0(plot_dir, "/", "Rank_", intr.rank, "_", intr.name, "_3d.pdf"), width = width, height = height)
+            pdf(paste0(plot_dir, "/", "Rank_", intr.rank, "_", intr.name, "-n_", sig.len, "-3D.pdf"), width = width, height = height)
         } else {
             stop("Plotting format not supported.\n")
         }
@@ -452,4 +454,6 @@ plotVelo <- function(
         dev.off()
 
     })
+
+    cat("End.\nFinished!\n")
 }
diff --git a/R/utility.r b/R/utility.r
index 2747ea5..5309a13 100644
--- a/R/utility.r
+++ b/R/utility.r
@@ -243,10 +243,15 @@ graphSpatialFDR <- function(nb.fac, pval.mtx, spatial.coords=NULL, weighting='DT
 # in the cells, 2) low number of significant interactions.
 filterRes <- function(dge.raw, res.list, intr.db, gene_to_uniprot, reads.thresh = 100, sig.thresh = 100){
     low_genes = diff(Matrix::t(dge.raw)@p)
-    low_genes = toupper(rownames(dge.raw)[which(low_genes < reads.thresh)])
-    low_genes = intersect(low_genes, gene_to_uniprot$gene_name)
+    
+    # # if dge.raw is genes X cell
+    # low_genes = toupper(rownames(dge.raw)[which(low_genes < reads.thresh)])
+    # low_genes = intersect(low_genes, gene_to_uniprot$gene_name)
+    # low_unpt = gene_to_uniprot$uniprot[gene_to_uniprot$gene_name %in% low_genes]
+
+    # if dge is UNPT X cell
+    low_unpt = rownames(dge.raw)[which(low_genes < reads.thresh)]
 
-    low_unpt = gene_to_uniprot$uniprot[gene_to_uniprot$gene_name %in% low_genes]
     low_unpt = low_unpt[!duplicated(low_unpt)]
     all.intrs = intr.db[["combined"]]
     low.intrs = levels(all.intrs[names(all.intrs) %in% low_unpt, drop = T])
@@ -363,3 +368,54 @@ getIntrNames <- function(res.intr.list, inter.index){
 
     return(intr.names)
 }
+
+
+# x is the new column names of the sparse matrix
+# y is the sparse matrix
+increase_columns_sparse <- function(x, y) {
+	# Get the column names of the sparse matrix
+    colnames_y <- colnames(y)
+	
+	# Determine which columns are missing
+	missing_cols <- setdiff(x, colnames_y)
+	
+	# Add missing columns filled with 0s
+	if (length(missing_cols) > 0) {
+		n_missing_cols <- length(missing_cols)
+		n_rows <- nrow(y)
+		y_new <- cbind(y, sparseMatrix(i = integer(0), j = integer(0), x = double(0), 
+										dims = c(n_rows, n_missing_cols), 
+										dimnames = list(NULL, missing_cols)))
+	} else {
+		y_new <- y
+	}
+	
+	# Reorder columns to match x
+	y_new[, x]
+}
+
+
+
+# x is the new column names of the matrix
+# y is the matrix
+increase_columns <- function(x, y) {
+	# Get the column names of the matrix
+	colnames_y <- colnames(y)
+	
+	# Determine which columns are missing
+	missing_cols <- setdiff(x, colnames_y)
+	
+	# Add missing columns filled with 0s
+	if (length(missing_cols) > 0) {
+		n_missing_cols <- length(missing_cols)
+		n_rows <- nrow(y)
+		y_new <- cbind(y, matrix(0, nrow = n_rows, ncol = n_missing_cols,
+								 dimnames = list(NULL, missing_cols)))
+	} else {
+		y_new <- y
+	}
+	
+	# Reorder columns to match x
+	y_new[, x]
+}
+
diff --git a/man/changeUniprot.Rd b/man/changeUniprot.Rd
index 65e409b..c917436 100644
--- a/man/changeUniprot.Rd
+++ b/man/changeUniprot.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/analysis.r
+% Please edit documentation in R/objects.r
 \name{changeUniprot}
 \alias{changeUniprot}
 \title{Subset the imputed data (intr.data) by the specified intr index}
diff --git a/man/changeUniprot.matrix_like.Rd b/man/changeUniprot.matrix_like.Rd
index 8cc85f1..7b52e9d 100644
--- a/man/changeUniprot.matrix_like.Rd
+++ b/man/changeUniprot.matrix_like.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/analysis.r
+% Please edit documentation in R/objects.r
 \name{changeUniprot.matrix_like}
 \alias{changeUniprot.matrix_like}
 \title{Sub function for changeUniprot, input a sparse matrix}
diff --git a/man/imputeNiche.dgCMatrix.Rd b/man/imputeNiche.dgCMatrix.Rd
index ad6b435..57381b2 100644
--- a/man/imputeNiche.dgCMatrix.Rd
+++ b/man/imputeNiche.dgCMatrix.Rd
@@ -8,7 +8,8 @@
   dge.raw,
   nb.id.fac,
   nb.dist.fac,
-  weights = c("mean", "counts", "dist", "none")
+  weights = c("mean", "counts", "dist", "none"),
+  return.graph = F
 )
 }
 \arguments{
diff --git a/man/imputeNicheVelo.CytoSignal.Rd b/man/imputeNicheVelo.CytoSignal.Rd
new file mode 100644
index 0000000..8e75631
--- /dev/null
+++ b/man/imputeNicheVelo.CytoSignal.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/analysis.r
+\name{imputeNicheVelo.CytoSignal}
+\alias{imputeNicheVelo.CytoSignal}
+\title{Sub function for imputeNicheVelo, input a Cytosignal object}
+\usage{
+\method{imputeNicheVelo}{CytoSignal}(
+  object,
+  nn.type = NULL,
+  weights = c("mean", "counts", "dist", "none")
+)
+}
+\arguments{
+\item{object}{A Cytosignal object}
+
+\item{nn.type}{The type of neighbors}
+
+\item{weights}{The weight of the Delaunay triangulation}
+}
+\value{
+A Cytosignal object
+}
+\description{
+Sub function for imputeNicheVelo, input a Cytosignal object
+}
diff --git a/man/imputeNicheVelo.Rd b/man/imputeNicheVelo.Rd
new file mode 100644
index 0000000..2665959
--- /dev/null
+++ b/man/imputeNicheVelo.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/analysis.r
+\name{imputeNicheVelo}
+\alias{imputeNicheVelo}
+\title{Impute the velocity mtx using the specified method}
+\usage{
+imputeNicheVelo(object, ...)
+}
+\arguments{
+\item{object}{A Cytosignal object}
+
+\item{...}{Other parameters}
+
+\item{method}{The method to use for imputation}
+}
+\value{
+A Cytosignal object
+}
+\description{
+Impute the velocity mtx using the specified method
+}
diff --git a/man/inferVeloLR.CytoSignal.Rd b/man/inferVeloLR.CytoSignal.Rd
new file mode 100644
index 0000000..993bd38
--- /dev/null
+++ b/man/inferVeloLR.CytoSignal.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/analysis.r
+\name{inferVeloLR.CytoSignal}
+\alias{inferVeloLR.CytoSignal}
+\title{Sub function for inferVeloLR, input a CytoSignal object}
+\usage{
+\method{inferVeloLR}{CytoSignal}(object, lig.slot, recep.slot, intr.db.name, tag = NULL)
+}
+\arguments{
+\item{object}{A Cytosignal object}
+
+\item{lig.slot}{The ligand slot to use}
+
+\item{recep.slot}{The receptor slot to use}
+
+\item{intr.db.name}{The intr database name to use}
+
+\item{nn.use}{slot that the neighbor index should be taken from, by default is the same as
+the recep.slot. For example, if velo.obj = GauEps-DT, then nn.use = "DT".
+nn.use could also be a user-defind factor.}
+}
+\value{
+A Cytosignal object
+}
+\description{
+Sub function for inferVeloLR, input a CytoSignal object
+}
diff --git a/man/inferVeloLR.Rd b/man/inferVeloLR.Rd
new file mode 100644
index 0000000..903a6a3
--- /dev/null
+++ b/man/inferVeloLR.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/analysis.r
+\name{inferVeloLR}
+\alias{inferVeloLR}
+\title{Compute the LR velo for specific ligand-receptor imputation obj pairs}
+\usage{
+inferVeloLR(object, ...)
+}
+\arguments{
+\item{object}{A Cytosignal object}
+
+\item{lig.slot}{The ligand slot to use}
+
+\item{recep.slot}{The receptor slot to use}
+
+\item{intr.db.name}{The intr database name to use}
+
+\item{nn.use}{The neighbor index as niche}
+}
+\value{
+A Cytosignal object
+}
+\description{
+Compute the LR velo for specific ligand-receptor imputation obj pairs
+}
diff --git a/man/inferVeloLR.matrix_like.Rd b/man/inferVeloLR.matrix_like.Rd
new file mode 100644
index 0000000..a0add44
--- /dev/null
+++ b/man/inferVeloLR.matrix_like.Rd
@@ -0,0 +1,32 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/analysis.r
+\name{inferVeloLR.matrix_like}
+\alias{inferVeloLR.matrix_like}
+\title{Sub function for inferVeloLR, input a sparse matrix}
+\usage{
+\method{inferVeloLR}{matrix_like}(
+  dge.lig,
+  dge.recep,
+  dge.lig.velo,
+  dge.recep.velo,
+  lig.fac,
+  recep.fac
+)
+}
+\arguments{
+\item{dge.lig}{A sparse matrix for ligand}
+
+\item{dge.recep}{A sparse matrix for receptor}
+
+\item{lig.fac}{A factor of ligand indices}
+
+\item{recep.fac}{A factor of receptor indices}
+
+\item{nb.id.fac}{A factor of neighbor indices}
+}
+\value{
+A sparse matrix
+}
+\description{
+Sub function for inferVeloLR, input a sparse matrix
+}
diff --git a/man/normCounts.CytoSignal.Rd b/man/normCounts.CytoSignal.Rd
index 902d241..4c9f27d 100644
--- a/man/normCounts.CytoSignal.Rd
+++ b/man/normCounts.CytoSignal.Rd
@@ -4,7 +4,7 @@
 \alias{normCounts.CytoSignal}
 \title{Sub function for normCounts, input a Cytosignal object}
 \usage{
-\method{normCounts}{CytoSignal}(object, method = c("default", "cpm"), slot.use = NULL)
+\method{normCounts}{CytoSignal}(object, method = c("default", "cpm"), slot.use = NULL, velo = FALSE)
 }
 \arguments{
 \item{object}{A Cytosignal object}
@@ -12,6 +12,8 @@
 \item{method}{The method to use for normalization}
 
 \item{slot.use}{The slot to use for normalization}
+
+\item{velo}{Whether to normalize the velocity data instead}
 }
 \value{
 A Cytosignal object
diff --git a/man/permuteImpEdge.Rd b/man/permuteImpEdge.Rd
new file mode 100644
index 0000000..ca54840
--- /dev/null
+++ b/man/permuteImpEdge.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/archive.r
+\name{permuteImpEdge}
+\alias{permuteImpEdge}
+\title{Permute Imputation Results of specific imputation method}
+\usage{
+permuteImpEdge(object, nn.type = NULL, perm.size = NULL, times = 10)
+}
+\arguments{
+\item{object}{A Cytosignal object}
+
+\item{nn.type}{The imputation method to use}
+
+\item{perm.size}{Size of the permutation test, by defualt NULL}
+
+\item{times}{Number of times to permute the whole dataset, by default 10}
+}
+\value{
+A Cytosignal object
+}
+\description{
+This function is a follow-up function of imputeNiche, which permutes the default or user-defined
+imputation method and stored the results in the ImpData object.
+}
diff --git a/man/permuteImpLR.Rd b/man/permuteImpLR.Rd
index 8292ed1..ac92b91 100644
--- a/man/permuteImpLR.Rd
+++ b/man/permuteImpLR.Rd
@@ -4,14 +4,16 @@
 \alias{permuteImpLR}
 \title{Permute Imputation Results of specific imputation method}
 \usage{
-permuteImpLR(object, nn.type = NULL, perm.size = 1e+05)
+permuteImpLR(object, nn.type = NULL, perm.size = NULL, times = 10)
 }
 \arguments{
 \item{object}{A Cytosignal object}
 
 \item{nn.type}{The imputation method to use}
 
-\item{perm.size}{Size of the permutation test}
+\item{perm.size}{Size of the permutation test, by defualt NULL}
+
+\item{times}{Number of times to permute the whole dataset, by default 10}
 }
 \value{
 A Cytosignal object
diff --git a/man/permuteNicheScoreLR.Rd b/man/permuteNicheScoreLR.Rd
index f62b80b..9460fc0 100644
--- a/man/permuteNicheScoreLR.Rd
+++ b/man/permuteNicheScoreLR.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/analysis.r
+% Please edit documentation in R/archive.r
 \name{permuteNicheScoreLR}
 \alias{permuteNicheScoreLR}
 \title{Permute LR score for specific ligand-receptor imputation obj pairs}
diff --git a/man/purgeBeforeSave.Rd b/man/purgeBeforeSave.Rd
index d3ea565..1c62805 100644
--- a/man/purgeBeforeSave.Rd
+++ b/man/purgeBeforeSave.Rd
@@ -4,7 +4,7 @@
 \alias{purgeBeforeSave}
 \title{Remove imputed data and normalized imputed data from CytoSignal object to save disk space}
 \usage{
-purgeBeforeSave(object)
+purgeBeforeSave(object, purge.raw = T, purge.null = F)
 }
 \arguments{
 \item{object}{CytoSignal object}
diff --git a/man/veloNicheLR.CytoSignal.Rd b/man/veloNicheLR.CytoSignal.Rd
index a519cff..3a07df6 100644
--- a/man/veloNicheLR.CytoSignal.Rd
+++ b/man/veloNicheLR.CytoSignal.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/analysis.r
+% Please edit documentation in R/archive.r
 \name{veloNicheLR.CytoSignal}
 \alias{veloNicheLR.CytoSignal}
 \title{Sub function for veloNicheLR, input a CytoSignal object}
diff --git a/man/veloNicheLR.Rd b/man/veloNicheLR.Rd
index fe42519..4e40d73 100644
--- a/man/veloNicheLR.Rd
+++ b/man/veloNicheLR.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/analysis.r
+% Please edit documentation in R/archive.r
 \name{veloNicheLR}
 \alias{veloNicheLR}
 \title{Compute the LR velo for specific ligand-receptor imputation obj pairs}
diff --git a/man/veloNicheLR.matrix_like.Rd b/man/veloNicheLR.matrix_like.Rd
index 0f80603..8106d4e 100644
--- a/man/veloNicheLR.matrix_like.Rd
+++ b/man/veloNicheLR.matrix_like.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/analysis.r
+% Please edit documentation in R/archive.r
 \name{veloNicheLR.matrix_like}
 \alias{veloNicheLR.matrix_like}
 \title{Sub function for veloNicheLR, input a sparse matrix}
diff --git a/src/utils_velo.cpp b/src/utils_velo.cpp
index 42aab59..eeb1a32 100644
--- a/src/utils_velo.cpp
+++ b/src/utils_velo.cpp
@@ -129,6 +129,58 @@ arma::mat graphMeanLR_cpp(
 
 
 
+
+// [[Rcpp::export]]
+arma::mat inferVeloLR_cpp(
+    const arma::mat& dge_lig,
+    const arma::mat& dge_recep,
+    const arma::mat& dge_lig_velo,
+    const arma::mat& dge_recep_velo,
+    const arma::uvec& lig_index,
+    const arma::uvec& lig_list,
+    const arma::uvec& recep_index,
+    const arma::uvec& recep_list
+){
+    // lig_dge & recep_dge: genes X beads
+    // res_mtx: bead X interactions mtx
+    arma::mat res_mtx = arma::zeros<arma::mat>(dge_lig.n_cols, lig_index.size()-1);
+
+    for (uword i = 0; i < res_mtx.n_cols; i++){ // for each ligand
+
+        for (uword j = 0; j < res_mtx.n_rows; j++){ // for each bead
+            // All indices in Cpp starts with 0 --> all values - 1
+
+            arma::vec lig_vec = dge_lig.col(j);
+            lig_vec = lig_vec.elem(lig_list(span(lig_index(i), lig_index(i+1)-1))); // n_ligs X 1
+            double lig_val = arma::accu(lig_vec); // n_ligs X 1
+
+            // same for dge_velo
+            arma::vec lig_velo = dge_lig_velo.col(j);
+            lig_velo = lig_velo.elem(lig_list(span(lig_index(i), lig_index(i+1)-1))); // n_ligs X 1
+            double lig_velo_val = arma::accu(lig_velo); // n_ligs X 1
+
+            arma::vec recep_vec = dge_recep.col(j);
+            recep_vec = recep_vec.elem(recep_list(span(recep_index(i), recep_index(i+1)-1))); // n_ligs X 1
+            double recep_val = arma::accu(recep_vec); // 1 X n_receps
+
+            arma::vec recep_velo = dge_recep_velo.col(j);
+            recep_velo = recep_velo.elem(recep_list(span(recep_index(i), recep_index(i+1)-1))); // n_ligs X 1
+            double recep_velo_val = arma::accu(recep_velo); // 1 X n_receps
+
+            // res_mtx(j, i) = arma::accu(lig_vec * recep_vec.t());
+            res_mtx(j, i) = lig_val * recep_velo_val + recep_val * lig_velo_val;
+        }
+        
+        // cout << i << endl;
+    }
+
+    return res_mtx;
+}
+
+
+
+
+
 // [[Rcpp::export]]
 arma::mat VelographNicheLR_cpp(
     const arma::mat& dge_lig,