R/deltaseq.R

In fhlsjs/Super-delta: Super-delta-sequence method

## This is the original delta-sequence function.  For simplicity, I only implemented t.test with two.sided hypothesis
deltaseq <- function(X, classlabel){
    ngenes <- nrow(X)
    ## Let's turn the classlabel into a binary vector
    classlabel <- as.integer(as.factor(classlabel)) - 1
    ## define Xa, Xb
    Xa <- X[, classlabel==0]; Xb <- X[, classlabel==1]
    n1 <- ncol(Xa); n2 <- ncol(Xb)
    ## gene-order is calculated via pooled sample variance/STD
    Sp <- (n1-1) *apply(Xa, 1, var) + (n2-1)* apply(Xb, 1, var)
    ## order and the reverse order
    o <- order(Sp); ro <- order(o)
    Xord <- X[o,]
    ndeltas <- ngenes %/%2
    if (ngenes %% 2==1){                #odd number of genes
        ## First kind of pairing.
        ## We take the row-differences and append the last one to it
        deltas1 <- rbind(Xord[seq(2, ndeltas*2, 2),] - Xord[seq(1, ndeltas*2-1, 2),], Xord[ngenes, ])
        ## compute the statistics at the deltas level
        delta.stats1 <- fastt(deltas1, classlabel)
        ## now translate them back to genes
        deltas2genes1 <- c(rep(1:ndeltas, each=2), ndeltas+1)
        ## multiply (-1,1, ..., -1, 1, 1) to ensure the sign of the
        ## t-stats are consistent.
        stats1 <- delta.stats1[deltas2genes1] * c(rep(c(-1,1), ndeltas), 1)
        ## the second kind of pairing
        deltas2 <- rbind(Xord[1,], Xord[seq(3, ndeltas*2+1, 2),] - Xord[seq(2, ndeltas*2, 2),])
        delta.stats2 <- fastt(deltas2, classlabel)
        deltas2genes2 <- c(1, rep(1:ndeltas, each=2))
        stats2 <- delta.stats2[deltas2genes2] * c(1, rep(c(-1,1), ndeltas))
    } else {                            #even number of genes
        ## First kind of pairing
        deltas1 <- Xord[seq(2, ndeltas*2, 2),] - Xord[seq(1, ndeltas*2-1, 2),]
        delta.stats1 <- fastt(deltas1, classlabel)
        deltas2genes1 <- rep(1:ndeltas, each=2)
        stats1 <- delta.stats1[deltas2genes1] * rep(c(-1,1), ndeltas)
        ## Second kind of pairing. Leave the top and bottom alone
        deltas2 <- rbind(Xord[1,], Xord[seq(3, ndeltas*2-1, 2),] - Xord[seq(2, (ndeltas-1)*2, 2),], Xord[ngenes,])
        delta.stats2 <- fastt(deltas2, classlabel)
        deltas2genes2 <- c(1, rep(2:ndeltas, each=2), ndeltas+1)
        stats2 <- delta.stats2[deltas2genes2] * c(1, rep(c(-1,1), ndeltas-1), 1)
    }
    ## The less extreme of the two tstats gets reported. Note that I
    ## have to reverse-order these stats so that they are consistent
    ## with the original gene list.
    ## 11/18/2015.  Multiply the test stats by sqrt(2) to boost power.    
    teststat <- sqrt(2)*apply(cbind(stats1, stats2), 1, function(tt) {at <- abs(tt); ifelse(at[1] <= at[2], tt[1], tt[2])})[ro]
    ## The raw pvalue. Only two-sided test is implemented at this moment.
    rawp=.t2p(teststat, df=ncol(X)-2, side="two.sided")
    rr <- data.frame(teststat=teststat, rawp=rawp)
    colnames(rr) <- c("teststat", "rawp")
    return(rr)
}