R/vca.R

In ziyili20/TOAST: Tools for the analysis of heterogeneous tissues

vca <- function(R, p, SNR=NULL, verbose=FALSE) {
     L <- nrow(R)
     N <- ncol(R)
     if (p < 0 || p > L) {
          stop("p is out of range (negative or too big)")
     }
     SNRth <- 15 + 10 * log10(p)

     if (!is.null(SNR) && (SNR < SNRth)) {
          if (verbose) message("Select the projective projection")
          d <- p - 1

          if (exists("xp")) {
               Ud <- Ud[, 1:d]
          } else {
               rm <- apply(R, 1, mean)
               R0 <- R - rm
               svdObj <- svd(R0 %*% t(R0), nu = p, nv = p)
               Ud <- svdObj$u
               xp <- t(Ud) %*% R0
          }

          Rp <- Ud %*% xp[1:d, ] + rm
          x <- xp[1:d, ]
          c <- sqrt(max(colSums(x^2)))
          y <- rbind(x, rep(c, N))
     } else {
          if (verbose) message("Select projection to p-1")
          d <- p
          Ud <- svd(R %*% t(R) / N, nu = d, nv = d)$u

          xp <- t(Ud) %*% R
          Rp <- Ud %*% xp[1:d, ]

          x <- xp
          u <- rowMeans(x)
          y <- x / matrix(kronecker(colSums(x * u), rep(1, d)), nrow=d)
     }

     ## VCA itself

     indice <- rep(0, p)
     A <- matrix(0, nrow=p, ncol=p)
     A[p, 1] <- 1

     for (i in 1:p) {
          w <- matrix(runif(p), ncol=1)
          f <- w - A %*% pseudoinverse(A) %*% w;
          f <- f / sqrt(sum(f^2))

          v <- t(f) %*% y
          indice[i] <- which.max(abs(v))
          A[, i] <- y[, indice[i]]
     }
     Ae = Rp[, indice]
     return(Ae)
}