R/SXTuni.R

In jaspershen/SXTuni: Give the PLS score plot and other values.

SXTuni <- function(subject = subject,
                   test.method = "t",
                   adjust.method = "fdr",
                   log.scale = FALSE,
                   info = info,
                   path = NULL) {
  # packages <- library()[[2]][,1]
  # if (!any(packages == "pROC")) {
  #   install.packages("pROC")}
  # library(pROC)
  # browser()

  if (is.null(path)) {
    path <- getwd()
  }else{
    dir.create(path)
  }

  idx <- NULL
  for (i in 1:length(info)) {
    idx1 <- match(info[[i]], colnames(subject))
    idx <- c(idx, idx1)
  }
  idx <- sort(idx[!is.na(idx)])
  subject <- subject[,idx]

  if (test.method == "t" | test.method == "wilcox") {
    group1.index <- match(info[[1]], colnames(subject))
    group1.index <- group1.index[!is.na(group1.index)]
    group2.index <- match(info[[2]], colnames(subject))
    group2.index <- group2.index[!is.na(group2.index)]

    if (length(info) == 2){
      fc <- apply(subject, 1, function(x) {median(x[group2.index])/median(x[group1.index])})
      save(fc, file = file.path(path, "fc"))
    }


    Y <- NULL
    Y[group1.index] <- 0
    Y[group2.index] <- 1

    ##log transformation
    if (log.scale == FALSE) {
      subject <- subject
    }

    if (log.scale == "e") {
      subject <- log(subject + 1)
    }

    if (log.scale != FALSE & log.scale != "e"){
      subject <- log(subject + 1, as.numeric(log.scale))
    }

    if (test.method == "t"){
      subject.test <- apply(subject, 1, function(x) {t.test(x[group1.index], x[group2.index])})
    }
    if (test.method == "wilcox"){
      subject.test <- apply(subject, 1, function(x) {wilcox.test(x[group1.index], x[group2.index])})
    }

    p <- unlist(lapply(subject.test, function(x) x$p.value))
    p.cor <- p.adjust(p = p, method = adjust.method)
    save(p, file = file.path(path, "p"))
    save(p.cor, file = file.path(path, "p.cor"))
  }

  if (test.method == "anova") {

    subject.test <- list()
    aov.p <- NULL
    subject.tuk <- list()

    tuk.p <- matrix(nrow = nrow(subject),ncol = choose(length(info),2))

    Y <- NULL
    name <- colnames(subject)

    for (i in 1:length(info)) {
      Y[match(info[[i]], name)] <- i-1
    }
    cat("ANOVA %: \n")
    for (i in 1:nrow(subject)) {
      count <- floor(nrow(subject)*c(seq(0,1,0.1)))
      if (any(i == count)) {cat(ceiling(i*100/nrow(subject)))
        cat(" ")}
      subject.test[[i]] <- aov(as.numeric(subject[i,]) ~ as.factor(Y))
      aov.p[i] <- summary(subject.test[[i]])[[1]][,5][1]
      subject.tuk[[i]] <- TukeyHSD(subject.test[[i]])
      tuk.p[i,] <- subject.tuk[[i]][[1]][,4]
      colnames(tuk.p) <- names(subject.tuk[[i]][[1]][,4])
    }
    save(aov.p, file = file.path(path, "aov.p"))
    save(tuk.p, file = file.path(path, "tuk.p"))
  }

  feature.auc <- apply(subject, 1, function(x) auc(Y,x))
  save(feature.auc, file = file.path(path, "feature.auc"))
}