#cybert

Cybert vs. Doug.

1  Introduction

CyberT works by calculating a bayesian regularized t-statistic for each gene. This is a lot of language to basically say that we calculate a t-statistic for each gene, but we tweak the variance of each group by using genes with similar expression values.

The Mann-Whitney (or Wilcoxon Rank-Sum, or Mann-Whitney-Wilcoxon (MWW)) test is the non-parametric version of the t-test. Non-parametric in this case means that the statistic works on the ranks of the data, rather than assuming some parameterized underlying distribution. In contrast, proper analysis of the t-statistic is dependent on assuming that each group is from a normal distribution, i.e., is from a distribution parameterized by a mean and variance. The MWW statistic is very closely related to the area under the receiver operating characteristic curve, i.e., the area under the ROC curve, or AUC. In fact, if x is a vector of predictions on a set of positive examples, and y is a vector of predictions on a set of negative examples, then:

AUC(x,y) = MWW(x, y)/(|x| * |y|)

In this post, I'm going to demonstrate how to calculate the AUC for each gene in a microarray dataset. This basically says, if we used a single gene's expression value as a predictor for two different phenotypes, what would the AUC of this predictor be? Or another way to think of it is as a MWW statistic that is normalized to be between 0 and 1, or really between 0.5 and 1. One can also calculate p-values for these AUC's based on the p-values of the MWW stats.

2  The Function

Here is a function to do this calculation (in R of course). This actually calculates a one-vs-rest for all groups.

> ###################################################################### > #   calcAUCArray > # > #   Function to calculate the AUC by row for a matrix of predictions > #    > #   Input: > #       x           -   the data.frame or array  > #       repPerCond  -   a vector of the number per condition (assumes  > #                       the data.frame is in order). > ###################################################################### > calcAUCArray <- function(x, repPerCond, calcPVals = T){ +     classFactor <- factor(rep(1:length(repPerCond), repPerCond)) +      +     #The Auc is equiv to the Mann-Whitney or Wilcox test.   +     #For generalizations beyond 2 class, do a one-vs-rest scheme. +     classIndices <- tapply(1:ncol(x), classFactor, list) +      +     theRes <- sapply(classIndices, function(inds){ +         sapply(as.data.frame(t(x)), function(y){ +             maxW <- length(inds)*length(y[-inds]); +             theTest <-  wilcox.test(y[inds], y[-inds]); +             res <- theTest$statistic/maxW; +              +         }) +     }) +      +     ## This line flips the prediction (AUC < 0.5 is worse than random) +     theRes <- apply(theRes, 2, function(x){ifelse(x > 0.5, x, 1-x)}) +     colnames(theRes) <- paste("AUC_Cond_", 1:length(repPerCond), sep=''); +     if (calcPVals){ +         thePVals <- sapply(classIndices, function(inds){ +             sapply(as.data.frame(t(x)), function(y){ +                 theTest <-  wilcox.test(y[inds], y[-inds]); +                 theTest$p.value; +             }) +         }) +         colnames(thePVals) <- paste("AUC_PVal_Cond_", 1:length(repPerCond), sep=''); +         theRes <- cbind(theRes, thePVals) +     }   +     theRes; + }   

3  An Example

Here, let's fake up some data and show how the MWW could be useful in finding some genes of interest. First, make some fake completely normal data:

> data <- data.frame(matrix(rnorm(100), nrow=10)) > colnames(data) <- c(sprintf('Ctrl%d', 1:5),  +                     sprintf('Expt%d', 1:5)) > head(data)  

        Ctrl1       Ctrl2       Ctrl3      Ctrl4      Ctrl5      Expt1 1  1.05075025 -1.30088148  2.04442575  0.2396374 -0.1640639 0.46255951 2 -0.04274524  0.14674350  0.19292496 -3.0934266  1.3382232 1.44391446 3  0.32498054 -0.38788724 -1.18156093  0.4588947  0.4427220 1.74483569 4  1.61588825 -0.78319953 -0.07895495 -0.3512329  1.7633533 0.08020991 5 -0.15815873 -0.71218677 -0.54358175  1.1996652 -0.9436978 0.64061365 6 -1.68811312 -0.01021523  0.81717540  1.0148944 -1.6237105 1.39887543        Expt2      Expt3       Expt4        Expt5 1  0.1977868  1.2951139 -0.94278682 -0.458105214 2  1.1318220 -0.3654534  0.45353342 -0.514351083 3 -0.9689614 -0.7129198  0.07177276 -1.175522844 4 -1.7997450 -0.4597961  0.56408353  1.044657360 5 -0.4004646  0.5248907  0.49722154  0.009011202 6  0.5671171  0.9390662 -0.12485235  0.362181267  

Then, let's randomly pick a gene and just right shift the expression of just the Expt cols here by 2.

> fakeGeneId <- sample(1:10, 1) > fakeGeneId  

[1] 10  

> data[fakeGeneId, 6:10] <- data[fakeGeneId, 6:10] + 2;  

Now, let's calculate the AUCs and p-values here.

> repPerCond <- c(5,5) # 5 in each group > aucs <- calcAUCArray(data, repPerCond) > aucs  

      AUC_Cond_1 AUC_Cond_2 AUC_PVal_Cond_1 AUC_PVal_Cond_2 V1.W        0.56       0.56     0.841269841     0.841269841 V2.W        0.60       0.60     0.690476190     0.690476190 V3.W        0.60       0.60     0.690476190     0.690476190 V4.W        0.60       0.60     0.690476190     0.690476190 V5.W        0.76       0.76     0.222222222     0.222222222 V6.W        0.68       0.68     0.420634921     0.420634921 V7.W        0.52       0.52     1.000000000     1.000000000 V8.W        0.80       0.80     0.150793651     0.150793651 V9.W        0.56       0.56     0.841269841     0.841269841 V10.W       1.00       1.00     0.007936508     0.007936508  

> aucs[fakeGeneId,]  

     AUC_Cond_1      AUC_Cond_2 AUC_PVal_Cond_1 AUC_PVal_Cond_2      1.000000000     1.000000000     0.007936508     0.007936508   

There you go! AUCs and you should see a much lower p-value for the fakeGeneId!

File translated from TEX by TTH, version 3.85. On 7 Apr 2011, 17:54.