Small GSEA implementation

enrichment.R

ES <- function(p, w, pws, both=FALSE) {
    n <- length(pws)
    nr <- sum(abs(w[pws == p]))
    nh <- sum(pws == p)

    scores <- vector(length=n)
    scores[pws == p] <- abs(w[pws == p])/nr
    scores[pws != p] <- -1/(n - nh)
    r <- range(cumsum(scores))
    i <- which.max(abs(r))
    if (both) i <- 1:2

    r[i]
}

NES <- function(p, w, pws, n=200) {
    es <- ES(p, w, pws)
    norm <- replicate(n, ES(p, w, sample(pws), both=TRUE))

    if (es < 0) {
        no <- norm[norm < 0]
        nes <- -es/mean(no)
        pval <- sum(no < es)/length(no)
    } else {
        no <- norm[norm >= 0]
        nes <- es/mean(no)
        pval <- sum(no > es)/length(no)
    }

    return(c(nes, pval))
}