Greedy Pair (NPair=1)

GP1.R

# Input: 
#   x: gene expression (genes x samples)
#   bin: binary response (eg Resistant/Sensitive)
#   GENEID: gene id data frame (for displaying gene symbols only)

# Output:
#   two gene symbols
GP1 <- function(x, bin, GENEID) {
  require(sva)
  pheno <- data.frame(bin = bin)
  mod = model.matrix(~as.factor(bin), data = pheno)
  mod0 = model.matrix(~1, data = pheno)

  x <- as.matrix(x) # better not use a data frame as input
  pv <- f.pvalue(x, mod, mod0)
  gene1 <- which.min(pv)
  GENE1 <- GENEID[gene1, 'Symbol'] # [1] "ZBTB6"

  i1 <- bin == "Resistant"
  i2 <- bin == "Sensitive"
  n1 <- sum(i1)
  n2 <- sum(i2)
  var1 <- rowVars(x[, i1]) * (n1-1)
  var2 <- rowVars(x[, i2]) * (n2-1)
  var12 <- (var1 + var2)/(n1+n2-2)
  mu1 <- rowMeans(x[, i1])
  mu2 <- rowMeans(x[, i2])
  a <- (mu1 - mu2) / var12

  ax <- x * a
  xnew <- matrix(rep(a[gene1]*x[gene1, ], nrow(x)), byrow=T, nr=nrow(x)) + ax

  pv2 <- f.pvalue(xnew, mod, mod0)
  gene2 <- which.min(pv2)  # assume gene1 != gene2; o/w we need to exclude gene1 before searching gene2
  GENE2 <- GENEID[gene2, 'Symbol'] 

  out <- list(index = c(gene1, gene2), symbol = c(GENE1, GENE2))
}