R/standardScreeningBinaryTrait.R

In WGCNA: Weighted Correlation Network Analysis

# The function standardScreeningBinaryTrait computes widely used statistics for relating the columns of
# the input data frame (argument datExpr) to a binary sample trait (argument y). The statistics include
# Student t-test p-value and the corresponding local false discovery rate (known as q-value, Storey et al
# 2004), the fold change, the area under the ROC curve (also known as C-index), mean values etc. If the
# input option kruskalTest is set to TRUE, it also computes the kruskal Wallist test p-value and
# corresponding q-value. The kruskal Wallis test is a non-parametric, rank-based group comparison test.

standardScreeningBinaryTrait=function(datExpr, y, 
           corFnc = cor, corOptions = list(use = 'p'),
           kruskalTest=FALSE, qValues = FALSE, var.equal=FALSE, na.action="na.exclude",
           getAreaUnderROC = TRUE) 
{

   datExpr=data.frame(datExpr, check.names = FALSE)
   levelsy=levels(factor(y))
   if (length(levelsy)>2 ) 
     stop("The sample trait y contains more than 2 levels. Please input a binary variable y")
   if (length(levelsy)==1 ) 
     stop("The sample trait y is constant. Please input a binary sample trait with some variation.")
   yNumeric=as.numeric(factor(y));
   if (length(yNumeric) !=dim(datExpr)[[1]] ) 
     stop("the length of the sample trait y does not equal the number of rows of datExpr")
   pvalueStudent = t.Student = Z.Student = rep(NA, dim(datExpr)[[2]] ) 
   pvaluekruskal = stat.Kruskal = Z.Kruskal = sign.Kruskal = rep(NA, dim(datExpr)[[2]] ) 
   nPresent = rep(0, dim(datExpr)[[2]] ) 
   AreaUnderROC=rep(NA, dim(datExpr)[[2]] ) 

  if (var.equal) 
     printFlush(paste("Warning: T-test that assumes equal variances in each group is requested.\n",
                      "This is not the default option for t.test. We recommend to use var.equal=FALSE."));

  corFnc = match.fun(corFnc);
  corOptions$y = yNumeric;
  corOptions$x = datExpr;
  corPearson=as.numeric(do.call(corFnc, corOptions));
  nGenes = dim(datExpr)[[2]];
  nPresent1 = as.numeric( t(as.matrix(!is.na(yNumeric) & yNumeric==1)) %*% ! is.na(datExpr) );
  nPresent2 = as.numeric( t(as.matrix(!is.na(yNumeric) & yNumeric==2)) %*% ! is.na(datExpr) );
  nPresent = nPresent1 + nPresent2;

  for (i in 1:nGenes) {
        no.present1 = nPresent1[i];
        no.present2 = nPresent2[i];
        no.present = nPresent[i];
        if (no.present1<2 | no.present2<2 ) 
        {
           pvalueStudent[i]= t.Student[i] = NA 
        } else {
          tst = try(t.test( as.numeric(datExpr[,i])~yNumeric,var.equal=var.equal,na.action=na.action),
                    silent = TRUE)
          if (!inherits(tst, "try-error"))
          {
            pvalueStudent[i] = tst$p.value;
            t.Student[i] = -tst$statistic 
            # The - sign above is intentional to make the sign of t consistent with correlation
          } else {
            printFlush(paste("standardScreeningBinaryTrait: An error ocurred in t.test for variable", 
                             i, ":\n", tst));
            printFlush(paste("Will return missing value(s) for this variable.\n\n"));
          }
          
        }
        if (getAreaUnderROC) AreaUnderROC[i] = rcorr.cens(datExpr[, i], yNumeric, outx = TRUE)[[1]]
        if (kruskalTest) {
            if (no.present<5 ) 
            {
               pvaluekruskal[i] = stat.Kruskal[i] = NA
            } else {
               kt = try(kruskal.test(datExpr[, i] ~ factor(yNumeric),  na.action="na.exclude"), silent = TRUE)
               if (!inherits(kt, "try-error"))
               {
                 pvaluekruskal[i] = kt$p.value;
                 stat.Kruskal[i] = kt$statistic;
                 # Find which side is higher
                 r = rank(datExpr[, i]);
                 means = tapply(r, factor(yNumeric), mean, na.rm = TRUE);
                 sign.Kruskal[i] = 2 * ( (means[1] < means[2]) - 0.5);
                 # sign.Kruskal is 1 if the ranks in group 1 are smaller than in group 2
               } else {
                 printFlush(paste("standardScreeningBinaryTrait: An error ocurred in kruskal.test for variable",
                                  i, ":\n", kt));
                 printFlush(paste("Will return missing value(s) for this variable.\n\n"));
               }

           }
        } 
    }
    q.Student=rep(NA, length(pvalueStudent) )
    rest1= ! is.na(pvalueStudent) 
    if (qValues) 
    {
      x = try({ q.Student[rest1] = qvalue(pvalueStudent[rest1])$qvalues }, silent = TRUE);
      if (inherits(x, "try-error"))
        printFlush(paste("Warning in standardScreeningBinaryTrait: function qvalue returned an error.\n",
                         "calculated q-values will be invalid. qvalue error:\n\n", x, "\n")) 
      if (kruskalTest) {
         q.kruskal=rep(NA, length(pvaluekruskal) )
         rest1= ! is.na(pvaluekruskal) 
         xx = try( { q.kruskal[rest1] = qvalue(pvaluekruskal[rest1])$qvalues} , silent = TRUE);
         if (inherits(xx, "try-error"))
           printFlush(paste("Warning in standardScreeningBinaryTrait: function qvalue returned an error.\n",
                            "calculated q-values will be invalid. qvalue error:\n\n", xx, "\n")) 
      }
    }
    meanLevel1 = as.numeric(apply(datExpr[!is.na(y) & y == levelsy[[1]], ], 2, mean, na.rm = TRUE));
    meanLevel2 = as.numeric(apply(datExpr[!is.na(y) & y == levelsy[[2]], ], 2, mean, na.rm = TRUE));
  
    Z.Student = qnorm(pvalueStudent/2, lower.tail = FALSE) * sign(t.Student);
    if (kruskalTest)
      Z.Kruskal = qnorm(pvaluekruskal/2, lower.tail = FALSE) * sign(stat.Kruskal);

    
    stderr1=function(x) {no.present=sum(!is.na(x));
    if (no.present<2) out=NA else {out=sqrt(var(x,na.rm=TRUE)/no.present) }
    out } # end of function stderr1

    SE.Level1 = as.numeric(apply(datExpr[y == levelsy[[1]] & !is.na(y), ], 2, stderr1))
    SE.Level2 =as.numeric(apply(datExpr[y == levelsy[[2]] & !is.na(y), ], 2, stderr1))

    FoldChangeLevel1vsLevel2 = ifelse(meanLevel1/meanLevel2 > 1, meanLevel1/meanLevel2, -meanLevel2/meanLevel1)
    
    output = data.frame(ID = dimnames(datExpr)[[2]], corPearson = corPearson, 
        t.Student = t.Student,
        pvalueStudent = pvalueStudent, 
        FoldChange = FoldChangeLevel1vsLevel2, 
         meanFirstGroup = meanLevel1,
        meanSecondGroup = meanLevel2, 
         SE.FirstGroup = SE.Level1,
      SE.SecondGroup = SE.Level2);

    if (getAreaUnderROC)
       output$AreaUnderROC = AreaUnderROC;

    if (kruskalTest) {
        output = data.frame(output, stat.Kruskal = stat.Kruskal, 
                            stat.Kruskal.signed = sign.Kruskal * stat.Kruskal,
                            pvaluekruskal = pvaluekruskal);
    }

   if (qValues && !inherits(x, "try-error")) output=data.frame(output, q.Student)
   if (qValues &  kruskalTest ) {
      if ( !inherits(xx, "try-error")) output=data.frame(output, q.kruskal)
   }
   names(output)[3:5] = paste(names(output)[3:5], levelsy[[1]], "vs", levelsy[[2]], sep = ".")
   output = data.frame(output, nPresentSamples = nPresent);
   output
}