Add core compute + plot func for VeloCytoSignal

NAMESPACE

@@ -1,7 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(changeUniprot,CytoSignal)
-S3method(changeUniprot,dgCMatrix)
+S3method(changeUniprot,matrix_like)
 S3method(findNNDT,CytoSignal)
 S3method(findNNDT,matrix)
 S3method(findNNGauEB,CytoSignal)
@@ -14,13 +14,19 @@ S3method(inferSignif,CytoSignal)
 S3method(inferSignif,dgCMatrix)
 S3method(normCounts,CytoSignal)
 S3method(normCounts,dgCMatrix)
+S3method(veloNicheLR,CytoSignal)
+S3method(veloNicheLR,matrix_like)
 export(addIntrDB)
+export(addVelo)
 export(changeUniprot)
 export(createCytoSignal)
 export(findNNDT)
 export(findNNGauEB)
 export(findNNRaw)
 export(graphNicheLR)
+export(hex_bin)
+export(hex_coord)
+export(hex_pos)
 export(imputeNiche)
 export(inferEpsParams)
 export(inferSignif)
@@ -35,8 +41,16 @@ export(showImp)
 export(showLog)
 export(showScore)
 export(showUnpt)
+export(showVelo)
 export(suggestInterval)
+export(veloNicheLR)
 exportClasses(CytoSignal)
 exportClasses(ImpData)
 exportClasses(lrScores)
+exportClasses(lrVelo)
 exportMethods(show)
+importFrom(RTriangle,pslg)
+importFrom(RTriangle,triangulate)
+importFrom(Rcpp,evalCpp)
+importFrom(plyr,count)
+useDynLib(cytosignal)

R/analysis.r

@@ -188,7 +188,7 @@ findNNDT <- function(
 #' 
 #' @param cells.loc A matrix of cells location
 #' @param weight.sum The sum of the weights
-#' 
+#' @importFrom RTriangle triangulate pslg
 #' @return A list of neighbors
 #' @export
 findNNDT.matrix <- function(
@@ -559,6 +559,13 @@ normCounts.CytoSignal <- function(
 # 	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
@@ -584,7 +591,7 @@ changeUniprot <- function(
 #' @return A sparse matrix
 #' @export
 #' 
-changeUniprot.dgCMatrix <- function(
+changeUniprot.matrix_like <- function(
     dge.raw, 
     gene_to_uniprot,
     # mode = "unpt",
@@ -1204,3 +1211,200 @@ runPears.std <- function(
 	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 = VelographNicheLR_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)
+}
+