R/shadow.r

In viper: Virtual Inference of Protein-activity by Enriched Regulon analysis

#' Shadow analysis for msviper objects
#'
#' This function performs shadow analysis on msviper objects
#'
#' @param mobj msviper object generated by \code{msviper}
#' @param regulators This parameter represent different ways to select a subset of regulators for performing the shadow analysis, it can be either a p-value cutoff, the total number of regulons to be used for computing the shadow effect, or a vector of regulator ids to be considered
#' @param targets Integer indicating the minimum number of common targets to compute shadow analysis
#' @param shadow Number indicating the p-value threshold for the shadow effect
#' @param per Integer indicating the number of permutations
#' @param nullmodel Null model in marix format
#' @param minsize Integer indicating the minimum size allowed for the regulons
#' @param adaptive.size Logical, whether the target weight should be considered when computing the regulon size
#' @param iterative Logical, whether a two step analysis with adaptive redundancy estimation should be performed
#' @param seed Integer indicating the seed for the permutations, 0 for disable it
#' @param cores Integer indicating the number of cores to use (only 1 in Windows-based systems)
#' @param verbose Logical, whether progression messages should be printed in the terminal
#' @return An updated msviper object with an additional slot (shadow) containing the shadow pairs
#' @seealso \code{\link{msviper}}
#' @examples
#' data(bcellViper, package="bcellViper")
#' sig <- rowTtest(dset, "description", c("CB", "CC"), "N")$statistic
#' dnull <- ttestNull(dset, "description", c("CB", "CC"), "N", per=100) # Only 100 permutations to reduce computation time, but it is recommended to perform at least 1000 permutations
#' mra <- msviper(sig, regulon, dnull)
#' mra <- shadow(mra, regulators=10)
#' summary(mra)
#' @export

shadow <- function(mobj, regulators=.01, targets=10, shadow=.01, per=1000, nullmodel=NULL, minsize=NULL, adaptive.size=NULL, iterative=NULL, seed=1, cores=1, verbose=TRUE) {
    if (seed>0) set.seed(round(seed))
    if (is.null(minsize)) minsize <- mobj$param$minsize
    if (is.null(adaptive.size)) adaptive.size <- mobj$param$adaptive.size
    if (is.null(iterative)) iterative <- mobj$param$iterative
    if (is.null(nullmodel)) nullmodel <- mobj$nullmodel
    epval <- mobj$es$p.value
    pos <- grep("--", names(mobj$es$p.value))
    if (length(pos)>0) epval <- mobj$es$p.value[-pos]
    if (length(regulators)>1) tfs <- names(mobj$regulon)[names(mobj$regulon) %in% regulators]
    else {
        if (regulators>1) tfs <- names(epval)[order(epval)[1:round(regulators)]]
        else tfs <- names(epval)[epval<regulators]
    }
    if (length(tfs)<2) return(mobj)
    regind <- t(combn(tfs, 2))
    tmp <- apply(regind, 1, function(x, regul1) {
        length(which(names(regul1[[x[1]]]$tfmode) %in% names(regul1[[x[2]]]$tfmode)))
    }, regul1=mobj$regulon)
    regind <- filterRowMatrix(regind, tmp>targets)
    if (length(regind)==0) return(mobj)
    regind <- rbind(regind, regind[, 2:1])
    ss1 <- apply(abs(mobj$signature), 2, rank)/(nrow(mobj$signature)+1)
    ss2 <- apply(mobj$signature, 2, rank)/(nrow(mobj$signature)+1)
    ss2 <- 2*(ss2-.5)
    pb <- NULL
    if (verbose) {
        message("\nPerforming shadow analysis on ", nrow(regind), " regulon pairs by permutation testing.")
        message("Process started at ", date())
    }
    if (cores>1) {
        res <- mclapply(1:nrow(regind), function(i, regind, ss1, ss2, regul1, per) {
            x <- regind[i, ]
            tfm <- regul1[[x[2]]]$tfmode
            if (all(names(tfm) %in% names(regul1[[x[1]]]$tfmode))) return(1)
            ll <- regul1[[x[2]]]$likelihood
            pos <- match(names(tfm), rownames(ss1))
            ss1 <- filterRowMatrix(ss1, pos)
            ss2 <- filterRowMatrix(ss2, pos)
            dnull <- sapply(1:per, function(i, ll, samp) {
                ll[sample(length(ll), samp)] <- 0
                return(ll)
            }, ll=ll, samp=length(which(names(tfm) %in% names(regul1[[x[1]]]$tfmode))))
            dnull <- cbind(ll, dnull)
            dnull[names(tfm) %in% names(regul1[[x[1]]]$tfmode), 1] <- 0
            dnull <- t(t(dnull)/colSums(dnull))
            sum1 <- t(ss2) %*% (dnull*tfm)
            es <- colMeans(abs(sum1)+(t(ss1) %*% (dnull * (1-abs(tfm)))))
            x1 <- es[1]
            es <- es[-1]
            iqr <- quantile(es, c(.5, 5/length(es)))
            pd <- ecdf(es)
            a <- list(x=knots(pd), y=pd(knots(pd)))
            filtro <- a$x<iqr[1] & a$x>=iqr[2] & a$y<1
            spl <- smooth.spline(a$x[filtro], -log(a$y[filtro]), spar=.75)
            p <- exp(-predict(spl, x1)$y)
            pos <- which(x1>iqr[1])
            if (x1>iqr[1]) p <- pd(x1)
            return(p)
        }, regind=regind, ss1=ss1, ss2=ss2, regul1=mobj$regulon[which(names(mobj$regulon) %in% unique(unlist(regind)))], per=per, mc.cores=cores)
        res <- sapply(res, function(x) x)    
    }
    else {
        if (verbose) pb <- txtProgressBar(max=nrow(regind), style=3)
        res <- sapply(1:nrow(regind), function(i, regind, ss1, ss2, regul1, per, pb, verbose) {
            x <- regind[i, ]
            if (verbose) setTxtProgressBar(pb, i)
            tfm <- regul1[[x[2]]]$tfmode
            if (all(names(tfm) %in% names(regul1[[x[1]]]$tfmode))) return(1)
            ll <- regul1[[x[2]]]$likelihood
            pos <- match(names(tfm), rownames(ss1))
            ss1 <- filterRowMatrix(ss1, pos)
            ss2 <- filterRowMatrix(ss2, pos)
            dnull <- sapply(1:per, function(i, ll, samp) {
                ll[sample(length(ll), samp)] <- 0
                return(ll)
            }, ll=ll, samp=length(which(names(tfm) %in% names(regul1[[x[1]]]$tfmode))))
            dnull <- cbind(ll, dnull)
            dnull[names(tfm) %in% names(regul1[[x[1]]]$tfmode), 1] <- 0
            dnull <- t(t(dnull)/colSums(dnull))
            sum1 <- t(ss2) %*% (dnull*tfm)
            es <- colMeans(abs(sum1)+(t(ss1) %*% (dnull * (1-abs(tfm)))))
            x1 <- es[1]
            es <- es[-1]
            iqr <- quantile(es, c(.5, 5/length(es)))
            pd <- ecdf(es)
            a <- list(x=knots(pd), y=pd(knots(pd)))
            filtro <- a$x<iqr[1] & a$x>=iqr[2] & a$y<1
            spl <- smooth.spline(a$x[filtro], -log(a$y[filtro]), spar=.75)
            p <- exp(-predict(spl, x1)$y)
            pos <- which(x1>iqr[1])
            if (x1>iqr[1]) p <- pd(x1)
            return(p)
        }, regind=regind, ss1=ss1, ss2=ss2, regul1=mobj$regulon[which(names(mobj$regulon) %in% unique(unlist(regind)))], per=per, pb=pb, verbose=verbose)
    }
    regind <- filterRowMatrix(regind, res<shadow)
    if (nrow(regind)==0) return(mobj)
    regind <- filterRowMatrix(regind, !(paste(regind[, 1], regind[, 2], sep="-") %in% paste(regind[, 2], regind[, 1], sep="-")))
    if (nrow(regind)==0) return(mobj)
    regul1 <- mobj$regulon
    for (i in 1:nrow(regind)) {
        reg1 <- regul1[[regind[i, 1]]]
        reg2 <- regul1[[regind[i, 2]]]
        reg2$likelihood[names(reg2$tfmode) %in% names(reg1$tfmode)] <- 0
        regul1[[regind[i, 2]]] <- reg2
    }
    regul1 <- lapply(regul1, function(x) {
        filtro <- x$likelihood>0
        x$tfmode <- x$tfmode[filtro]
        x$likelihood <- x$likelihood[filtro]
        return(x)
    })
    regul1 <- regul1[sapply(regul1, function(x) length(x$tfmode))>0]
    res <- msviper(ges=mobj$signature, regulon=regul1, nullmodel=nullmodel, pleiotropy=FALSE, minsize=minsize, adaptive.size=adaptive.size, ges.filter=FALSE, synergy=0, cores=cores, verbose=verbose)
    mobj$regulon <- c(res$regulon, mobj$regulon[grep("--", names(mobj$regulon))])
    pos <- grep("--", names(mobj$es$nes))
    mobj$es$es <- c(res$es$es, mobj$es$es[pos])
    mobj$es$nes <- c(res$es$nes, mobj$es$nes[pos])
    mobj$es$nes.se <- c(res$es$nes.se, mobj$es$nes.se[pos])
    mobj$es$size <- c(res$es$size, mobj$es$size[pos])
    mobj$es$p.value <- c(res$es$p.value, mobj$es$p.value[pos])
    mobj$shadow <- regind
    return(mobj)
}