Skip to content

Commit

Permalink
fix delayed array
Browse files Browse the repository at this point in the history
  • Loading branch information
YingxinLin committed Apr 12, 2024
1 parent 73e9694 commit 4706c3a
Show file tree
Hide file tree
Showing 8 changed files with 637 additions and 629 deletions.
149 changes: 76 additions & 73 deletions R/fastRUVIII.R
Original file line number Diff line number Diff line change
Expand Up @@ -39,88 +39,91 @@


fastRUVIII <- function(Y, M, ctl, k = NULL, eta = NULL,
svd_k = 50, include.intercept = TRUE, average = FALSE,
BPPARAM = SerialParam(), BSPARAM = ExactParam(),
fullalpha = NULL, return.info = FALSE, inputcheck = TRUE) {

m <- nrow(Y)
n <- ncol(Y)
M <- ruv::replicate.matrix(M)
ctl <- tological(ctl, n)

## Check the inputs
if (inputcheck) {
if (sum(is.na(Y)) > 0) {
stop("Y contains missing values. This is not supported.")
}
if (sum(Y == Inf, na.rm = TRUE) + sum(Y == -Inf, na.rm = TRUE) >
0) {
stop("Y contains infinities. This is not supported.")
}
svd_k = 50, include.intercept = TRUE, average = FALSE,
BPPARAM = SerialParam(), BSPARAM = ExactParam(),
fullalpha = NULL, return.info = FALSE, inputcheck = TRUE) {

m <- nrow(Y)
n <- ncol(Y)
M <- ruv::replicate.matrix(M)
ctl <- tological(ctl, n)

## Check the inputs
if (inputcheck) {
if (sum(is.na(Y)) > 0) {
stop("Y contains missing values. This is not supported.")
}


## RUV1 is a reprocessing step for RUVIII
Y <- ruv::RUV1(Y, eta, ctl, include.intercept = include.intercept)

if (inherits(BSPARAM, "ExactParam")) {
svd_k <- min(m - ncol(M), sum(ctl), svd_k, na.rm = TRUE)
} else {
svd_k <- min(m - ncol(M), sum(ctl), na.rm = TRUE)
if (sum(Y == Inf, na.rm = TRUE) + sum(Y == -Inf, na.rm = TRUE) >
0) {
stop("Y contains infinities. This is not supported.")
}
}


## RUV1 is a reprocessing step for RUVIII
Y <- ruv::RUV1(Y, eta, ctl, include.intercept = include.intercept)

if (inherits(BSPARAM, "ExactParam")) {
svd_k <- min(m - ncol(M), sum(ctl), svd_k, na.rm = TRUE)
} else {
svd_k <- min(m - ncol(M), sum(ctl), na.rm = TRUE)
}


## m represent the number of samples/observations ncol(M)
## represent the number of replicates If the replicate matrix
## is such that we have more replicates than samples, then
## RUV3 is not appropriate, thus, we return the Original input
## matrix
if (ncol(M) >= m | k == 0) {
newY <- Y
fullalpha <- NULL
} else {

if (is.null(fullalpha))
{
## The main RUVIII process Applies the residual operator of a
## matrix M to a matrix Y Y0 has the same dimensions as Y,
## i.e. m rows (observations) and n columns (genes).

Y0 <- my_residop(Y, M)

svdObj <- BiocSingular::runSVD(x = Y0, k = svd_k, BPPARAM = BPPARAM, BSPARAM = BSPARAM)

fullalpha <- DelayedArray::t(svdObj$u[, seq_len(svd_k), drop = FALSE]) %*% Y
} ## End is.null(fullalpha)
###############

## m represent the number of samples/observations ncol(M)
## represent the number of replicates If the replicate matrix
## is such that we have more replicates than samples, then
## RUV3 is not appropriate, thus, we return the Original input
## matrix
if (ncol(M) >= m | k == 0) {
newY <- Y
fullalpha <- NULL
} else {

if (is.null(fullalpha))
{
## The main RUVIII process Applies the residual operator of a
## matrix M to a matrix Y Y0 has the same dimensions as Y,
## i.e. m rows (observations) and n columns (genes).

Y0 <- my_residop(Y, M)

svdObj <- BiocSingular::runSVD(x = Y0, k = svd_k, BPPARAM = BPPARAM, BSPARAM = BSPARAM)

fullalpha <- DelayedArray::t(svdObj$u[, seq_len(svd_k), drop = FALSE]) %*% Y
} ## End is.null(fullalpha)
###############
alpha <- fullalpha[seq_len(min(k, nrow(fullalpha))), , drop = FALSE]
ac <- alpha[, ctl, drop = FALSE]
W <- Y[, ctl] %*% DelayedArray::t(ac) %*% solve(ac %*% DelayedArray::t(ac))
newY <- Y - W %*% alpha
} ## End else(ncol(M) >= m | k == 0)

## If the users want to get all the informations relating to
## the RUV, it can be done here.
if (!return.info) {
return(newY)
} else {
return(list(newY = newY, M = M, fullalpha = fullalpha))
}
alpha <- fullalpha[seq_len(min(k, nrow(fullalpha))), , drop = FALSE]
ac <- alpha[, ctl, drop = FALSE]
W <- Y[, ctl] %*% DelayedArray::t(ac) %*% solve(ac %*% DelayedArray::t(ac))
Y <- DelayedArray(Y)
newY <- Y - DelayedArray::DelayedArray(W %*% alpha)
} ## End else(ncol(M) >= m | k == 0)

## If the users want to get all the informations relating to
## the RUV, it can be done here.
if (!return.info) {
return(newY)
} else {
return(list(newY = newY, M = M, fullalpha = fullalpha))
}
}
############################
tological <- function(ctl, n) {
ctl2 <- rep(FALSE, n)
ctl2[ctl] <- TRUE
return(ctl2)
ctl2 <- rep(FALSE, n)
ctl2[ctl] <- TRUE
return(ctl2)
}


#' @importFrom DelayedArray t
my_residop <- function(A, B){
tBB = DelayedArray::t(B) %*% B
tBB_inv = Matrix::solve(tBB)
BtBB_inv = B %*% tBB_inv
tBA = DelayedArray::t(B) %*% A
BtBB_inv_tBA = BtBB_inv %*% tBA
return(A - BtBB_inv_tBA)
}
tBB = DelayedArray::t(B) %*% B
tBB_inv = Matrix::solve(tBB)
BtBB_inv = B %*% tBB_inv
tBA = DelayedArray::t(B) %*% A
BtBB_inv_tBA = DelayedArray::DelayedArray(BtBB_inv %*% tBA)
A <- DelayedArray(A)
return(A - BtBB_inv_tBA)
}
4 changes: 2 additions & 2 deletions R/scMerge2.R
Original file line number Diff line number Diff line change
Expand Up @@ -400,9 +400,9 @@ getAdjustedMat <- function(exprsMat, fullalpha,
W <- Y_stand[, ctl] %*% DelayedArray::t(alpha[, ctl]) %*% solve(alpha[, ctl] %*% DelayedArray::t(alpha[, ctl]))

if (!is.null(return_subset_genes)) {
newY <- exprsMat[, return_subset_genes] - W %*% alpha[, return_subset_genes]
newY <- exprsMat[, return_subset_genes] - DelayedArray::DelayedArray(W %*% alpha[, return_subset_genes])
} else {
newY <- exprsMat - W %*% alpha
newY <- exprsMat - DelayedArray::DelayedArray(W %*% alpha)
}

newY <- t(newY)
Expand Down
Loading

0 comments on commit 4706c3a

Please sign in to comment.