--- a
+++ b/R/internal-functions-samr-adapted.R
@@ -0,0 +1,2238 @@
+# This script contains customized versions of functions found in the samr
+# package. This is necessary because samr seems to have been abandoned, so an
+# upstream collaboration doesn't seem possible at the time of writing.
+# ATTENTION: The source code in this file is licensed under LGPL-3.
+
+# ==============================================================================
+# Constants
+# ==============================================================================
+samr.const.twoclass.unpaired.response <- "Two class unpaired"
+samr.const.twoclass.paired.response <- "Two class paired"
+samr.const.oneclass.response <- "One class"
+samr.const.quantitative.response <- "Quantitative"
+samr.const.multiclass.response <- "Multiclass"
+samr.const.twoclass.unpaired.timecourse.response <-
+  "Two class unpaired timecourse"
+samr.const.twoclass.paired.timecourse.response <- "Two class paired timecourse"
+samr.const.oneclass.timecourse.response <- "One class timecourse"
+samr.const.survival.response <- "Survival"
+samr.const.patterndiscovery.response <- "Pattern discovery"
+
+# ==============================================================================
+# Functions
+# ==============================================================================
+
+#' @title Significance analysis of microarrays
+#' @description This function is an adaptation of `samr::samr`
+#' @param data Data object with components x- p by n matrix of features, one
+#' observation per column (missing values allowed); y- n-vector of outcome
+#' measurements; censoring.status- n-vector of censoring censoring.status
+#' (1= died or event occurred, 0=survived, or event was censored), needed for a
+#' censored survival outcome
+#' @param resp.type Problem type: "Quantitative" for a continuous parameter
+#' (Available for both array and sequencing data); "Two class unpaired" (for
+#' both array and sequencing data); "Survival" for censored survival outcome
+#' (for both array and sequencingdata); "Multiclass": more than 2 groups (for
+#' both array and sequencing data); "One class" for a single group (only for
+#' array data); "Two class paired" for two classes with paired observations
+#' (for both array and sequencing data); "Two class unpaired timecourse" (only
+#' for array data), "One class time course" (only for array data),
+#' "Two class.paired timecourse" (only for array data), or "Pattern discovery"
+#' (only for array data)
+#' @param assay.type Assay type: "array" for microarray data, "seq" for counts
+#' from sequencing
+#' @param s0 Exchangeability factor for denominator of test statistic; Default
+#' is automatic choice. Only used for array data.
+#' @param s0.perc Percentile of standard deviation values to use for s0; default
+#' is automatic choice; -1 means s0=0 (different from s0.perc=0, meaning
+#' s0=zeroeth percentile of standard deviation values= min of sd values.
+#' Only used for array data.
+#' @param nperms Number of permutations used to estimate false discovery rates
+#' @param center.arrays Should the data for each sample (array) be median
+#' centered at the outset? Default =FALSE. Only used for array data.
+#' @param testStatistic Test statistic to use in two class unpaired case.Either
+#' "standard" (t-statistic) or ,"wilcoxon" (Two-sample wilcoxon or Mann-Whitney
+#' test). Only used for array data.
+#' @param time.summary.type Summary measure for each time course: "slope", or
+#' "signed.area"). Only used for array data.
+#' @param regression.method Regression method for quantitative case: "standard",
+#' (linear least squares) or "ranks" (linear least squares on ranked data).
+#' Only used for array data.
+#' @param return.x Should the matrix of feature values be returned? Only useful
+#' for time course data, where x contains summaries of the features over time.
+#' Otherwise x is the same as the input data
+#' @param knn.neighbors Number of nearest neighbors to use for imputation of
+#' missing features values. Only used for array data.
+#' @param random.seed Optional initial seed for random number generator
+#' (integer)
+#' @param nresamp For assay.type="seq", number of resamples used to construct
+#' test statistic. Default 20. Only used for sequencing data.
+#' @param nresamp.perm For assay.type="seq", number of resamples used to
+#' construct test statistic for permutations. Default is equal to nresamp and it
+#' must be at most nresamp. Only used for sequencing data.
+#' @param xl.mode Used by Excel interface
+#' @param xl.time Used by Excel interface
+#' @param xl.prevfit Used by Excel interface
+#' @importFrom impute impute.knn
+sammy <- function(
+    data, resp.type = c(
+      "Quantitative", "Two class unpaired",
+      "Survival", "Multiclass", "One class", "Two class paired",
+      "Two class unpaired timecourse", "One class timecourse",
+      "Two class paired timecourse", "Pattern discovery"
+    ), assay.type = c(
+      "array",
+      "seq"
+    ), s0 = NULL, s0.perc = NULL, nperms = 100, center.arrays = FALSE,
+    testStatistic = c("standard", "wilcoxon"), time.summary.type = c(
+      "slope",
+      "signed.area"
+    ), regression.method = c("standard", "ranks"),
+    return.x = FALSE, knn.neighbors = 10, random.seed = NULL,
+    nresamp = 20, nresamp.perm = NULL, xl.mode = c(
+      "regular",
+      "firsttime", "next20", "lasttime"
+    ), xl.time = NULL, xl.prevfit = NULL) {
+  this.call <- match.call()
+  resp.type.arg <- match.arg(resp.type)
+  assay.type <- match.arg(assay.type)
+  xl.mode <- match.arg(xl.mode)
+  set.seed(random.seed)
+  if (is.null(nresamp.perm)) nresamp.perm <- nresamp
+  nresamp.perm <- min(nresamp, nresamp.perm)
+  if (xl.mode == "regular" | xl.mode == "firsttime") {
+    x <- NULL
+    xresamp <- NULL
+    ttstar0 <- NULL
+    evo <- NULL
+    ystar <- NULL
+    sdstar.keep <- NULL
+    censoring.status <- NULL
+    pi0 <- NULL
+    stand.contrasts <- NULL
+    stand.contrasts.star <- NULL
+    stand.contrasts.95 <- NULL
+    foldchange <- NULL
+    foldchange.star <- NULL
+    perms <- NULL
+    permsy <- NULL
+    eigengene <- NULL
+    eigengene.number <- NULL
+    testStatistic <- match.arg(testStatistic)
+    time.summary.type <- match.arg(time.summary.type)
+    regression.method <- match.arg(regression.method)
+    x <- data$x
+    y <- data$y
+    argy <- y
+    if (!is.null(data$eigengene.number)) {
+      eigengene.number <- data$eigengene.number
+    }
+    if (sum(is.na(x)) > 0) {
+      x <- impute.knn(x, k = knn.neighbors)
+      if (!is.matrix(x)) {
+        x <- x$data
+      }
+    }
+    are.blocks.specified <- FALSE
+    cond <- resp.type %in% c(
+      "One class", "Two class unpaired timecourse",
+      "One class unpaired timecourse", "Two class paired timecourse",
+      "Pattern discovery"
+    )
+    if (assay.type == "seq" & cond) {
+      stop(paste("Resp.type=", resp.type, " not allowed when assay.type='seq'"))
+    }
+    if (assay.type == "seq" & min(x) < 0) {
+      stop(paste("Negative values not allowed when assay.type='seq'"))
+    }
+    if (assay.type == "seq" & (sum(x %% 1 != 0) != 0)) {
+      stop("Non-integer values not alled when assay.type='seq'")
+    }
+    if (assay.type == "seq" & center.arrays) {
+      stop(paste("Centering  not allowed when assay.type='seq'"))
+    }
+    if (assay.type == "seq" & regression.method == "ranks") {
+      stop(paste("regression.method==ranks not allowed when assay.type='seq'"))
+    }
+    if (center.arrays) {
+      x <- scale(x, center = apply(x, 2, median), scale = FALSE)
+    }
+    depth <- rep(NA, ncol(x))
+    if (assay.type == "seq") {
+      message("Estimating sequencing depths...")
+      depth <- samr.estimate.depth(x)
+      message("Resampling to get new data matrices...")
+      xresamp <- resa(x, depth, nresamp = nresamp)
+    }
+    if (resp.type == samr.const.twoclass.unpaired.response) {
+      if (substring(y[1], 2, 6) == "Block" | substring(
+        y[1],
+        2, 6
+      ) == "block") {
+        junk <- parse.block.labels.for.2classes(y)
+        y <- junk$y
+        blocky <- junk$blocky
+        are.blocks.specified <- TRUE
+      }
+    }
+    if (resp.type == samr.const.twoclass.unpaired.response |
+      resp.type == samr.const.twoclass.paired.response |
+      resp.type == samr.const.oneclass.response |
+      resp.type == samr.const.quantitative.response |
+      resp.type == samr.const.multiclass.response
+    ) {
+      y <- as.numeric(y)
+    }
+    sd.internal <- NULL
+    if (resp.type == samr.const.twoclass.unpaired.timecourse.response |
+      resp.type == samr.const.twoclass.paired.timecourse.response |
+      resp.type == samr.const.oneclass.timecourse.response) {
+      junk <- parse.time.labels.and.summarize.data(
+        x, y,
+        resp.type, time.summary.type
+      )
+      y <- junk$y
+      x <- junk$x
+      sd.internal <- sqrt(rowMeans(junk$sd^2))
+      if (min(table(y)) == 1) {
+        warning(
+          "Only one timecourse in one or more classes;\n",
+          "SAM plot and FDRs will be unreliable;",
+          "only gene scores are informative"
+        )
+      }
+    }
+    if (resp.type == samr.const.twoclass.unpaired.timecourse.response) {
+      resp.type <- samr.const.twoclass.unpaired.response
+    }
+    if (resp.type == samr.const.twoclass.paired.timecourse.response) {
+      resp.type <- samr.const.twoclass.paired.response
+    }
+    if (resp.type == samr.const.oneclass.timecourse.response) {
+      resp.type <- samr.const.oneclass.response
+    }
+    stand.contrasts <- NULL
+    stand.contrasts.95 <- NULL
+    if (resp.type == samr.const.survival.response) {
+      censoring.status <- data$censoring.status
+    }
+    check.format(
+      y,
+      resp.type = resp.type, censoring.status = censoring.status
+    )
+    if (resp.type == samr.const.quantitative.response & regression.method ==
+      "ranks") {
+      y <- rank(y)
+      x <- t(apply(x, 1, rank))
+    }
+    n <- nrow(x)
+    sd <- NULL
+    numer <- NULL
+    if (resp.type == samr.const.twoclass.unpaired.response &
+      testStatistic == "standard" & assay.type == "array") {
+      init.fit <- ttest.func(x, y, sd = sd.internal)
+      numer <- init.fit$numer
+      sd <- init.fit$sd
+    }
+    if (resp.type == samr.const.twoclass.unpaired.response &
+      testStatistic == "wilcoxon" & assay.type == "array") {
+      init.fit <- wilcoxon.func(x, y)
+      numer <- init.fit$numer
+      sd <- init.fit$sd
+    }
+    if (resp.type == samr.const.oneclass.response & assay.type ==
+      "array") {
+      init.fit <- onesample.ttest.func(x, y, sd = sd.internal)
+      numer <- init.fit$numer
+      sd <- init.fit$sd
+    }
+    if (resp.type == samr.const.twoclass.paired.response &
+      assay.type == "array") {
+      init.fit <- paired.ttest.func(x, y, sd = sd.internal)
+      numer <- init.fit$numer
+      sd <- init.fit$sd
+    }
+    if (resp.type == samr.const.survival.response & assay.type ==
+      "array") {
+      init.fit <- cox.func(x, y, censoring.status)
+      numer <- init.fit$numer
+      sd <- init.fit$sd
+    }
+    if (resp.type == samr.const.multiclass.response & assay.type ==
+      "array") {
+      init.fit <- multiclass.func(x, y)
+      numer <- init.fit$numer
+      sd <- init.fit$sd
+    }
+    if (resp.type == samr.const.quantitative.response & assay.type ==
+      "array") {
+      init.fit <- quantitative.func(x, y)
+      numer <- init.fit$numer
+      sd <- init.fit$sd
+    }
+    if (resp.type == samr.const.patterndiscovery.response &
+      assay.type == "array") {
+      init.fit <- patterndiscovery.func(x)
+      numer <- init.fit$numer
+      sd <- init.fit$sd
+    }
+    if ((resp.type == samr.const.quantitative.response &
+      (testStatistic == "wilcoxon" | regression.method ==
+        "ranks" & assay.type == "array") | resp.type ==
+      samr.const.patterndiscovery.response) | resp.type ==
+      samr.const.twoclass.unpaired.response & assay.type ==
+      "array" & testStatistic == "wilcoxon" | (nrow(x) <
+      500) & is.null(s0) & is.null(s0.perc)) {
+      s0 <- quantile(sd, 0.05)
+      s0.perc <- 0.05
+    }
+    if (is.null(s0) & assay.type == "array") {
+      if (!is.null(s0.perc)) {
+        if ((s0.perc != -1 & s0.perc < 0) | s0.perc >
+          100) {
+          stop(
+            "Illegal value for s0.perc: must be between 0 and 100, ",
+            "or equal\nto (-1) (meaning that s0 should be set to zero)"
+          )
+        }
+        if (s0.perc == -1) {
+          s0 <- 0
+        }
+        if (s0.perc >= 0) {
+          s0 <- quantile(init.fit$sd, s0.perc / 100)
+        }
+      }
+      if (is.null(s0.perc)) {
+        s0 <- est.s0(init.fit$tt, init.fit$sd)$s0.hat
+        s0.perc <- 100 * sum(init.fit$sd < s0) / length(init.fit$sd)
+      }
+    }
+    if (assay.type == "seq") {
+      s0 <- 0
+      s0.perc <- 0
+    }
+    if (resp.type == samr.const.twoclass.unpaired.response &
+      testStatistic == "standard" & assay.type == "array") {
+      tt <- ttest.func(x, y, s0 = s0, sd = sd.internal)$tt
+    }
+    if (resp.type == samr.const.twoclass.unpaired.response &
+      testStatistic == "wilcoxon" & assay.type == "array") {
+      tt <- wilcoxon.func(x, y, s0 = s0)$tt
+    }
+    if (resp.type == samr.const.oneclass.response & assay.type ==
+      "array") {
+      tt <- onesample.ttest.func(x, y, s0 = s0, sd = sd.internal)$tt
+    }
+    if (resp.type == samr.const.twoclass.paired.response &
+      assay.type == "array") {
+      tt <- paired.ttest.func(x, y, s0 = s0, sd = sd.internal)$tt
+    }
+    if (resp.type == samr.const.survival.response & assay.type ==
+      "array") {
+      tt <- cox.func(x, y, censoring.status, s0 = s0)$tt
+    }
+    if (resp.type == samr.const.multiclass.response & assay.type ==
+      "array") {
+      junk2 <- multiclass.func(x, y, s0 = s0)
+      tt <- junk2$tt
+      stand.contrasts <- junk2$stand.contrasts
+    }
+    if (resp.type == samr.const.quantitative.response & assay.type ==
+      "array") {
+      tt <- quantitative.func(x, y, s0 = s0)$tt
+    }
+    if (resp.type == samr.const.patterndiscovery.response &
+      assay.type == "array") {
+      junk <- patterndiscovery.func(
+        x, s0 = s0, eigengene.number = eigengene.number
+      )
+      tt <- junk$tt
+      eigengene <- junk$eigengene
+    }
+    if (resp.type == samr.const.twoclass.unpaired.response &
+      assay.type == "seq") {
+      junk <- wilcoxon.unpaired.seq.func(xresamp, y)
+      tt <- junk$tt
+      numer <- junk$numer
+      sd <- junk$sd
+    }
+    if (resp.type == samr.const.twoclass.paired.response &
+      assay.type == "seq") {
+      junk <- wilcoxon.paired.seq.func(xresamp, y)
+      tt <- junk$tt
+      numer <- junk$numer
+      sd <- junk$sd
+    }
+    if (resp.type == samr.const.quantitative.response & assay.type ==
+      "seq") {
+      junk <- quantitative.seq.func(xresamp, y)
+      tt <- junk$tt
+      numer <- junk$numer
+      sd <- junk$sd
+    }
+    if (resp.type == samr.const.survival.response & assay.type ==
+      "seq") {
+      junk <- cox.seq.func(xresamp, y, censoring.status)
+      tt <- junk$tt
+      numer <- junk$numer
+      sd <- junk$sd
+    }
+    if (resp.type == samr.const.multiclass.response & assay.type ==
+      "seq") {
+      junk2 <- multiclass.seq.func(xresamp, y)
+      tt <- junk2$tt
+      numer <- junk2$numer
+      sd <- junk2$sd
+      stand.contrasts <- junk2$stand.contrasts
+    }
+    if (
+      resp.type == samr.const.quantitative.response |
+        resp.type == samr.const.multiclass.response |
+        resp.type == samr.const.survival.response
+    ) {
+      junk <- getperms(y, nperms)
+      perms <- junk$perms
+      all.perms.flag <- junk$all.perms.flag
+      nperms.act <- junk$nperms.act
+    }
+    if (resp.type == samr.const.twoclass.unpaired.response) {
+      if (are.blocks.specified) {
+        junk <- compute.block.perms(y, blocky, nperms)
+        permsy <- matrix(junk$permsy, ncol = length(y))
+        all.perms.flag <- junk$all.perms.flag
+        nperms.act <- junk$nperms.act
+      } else {
+        junk <- getperms(y, nperms)
+        permsy <- matrix(y[junk$perms], ncol = length(y))
+        all.perms.flag <- junk$all.perms.flag
+        nperms.act <- junk$nperms.act
+      }
+    }
+    if (resp.type == samr.const.oneclass.response) {
+      if ((length(y) * log(2)) < log(nperms)) {
+        allii <- 0:((2^length(y)) - 1)
+        nperms.act <- 2^length(y)
+        all.perms.flag <- 1
+      } else {
+        nperms.act <- nperms
+        all.perms.flag <- 0
+      }
+      permsy <- matrix(NA, nrow = nperms.act, ncol = length(y))
+      if (all.perms.flag == 1) {
+        k <- 0
+        for (i in allii) {
+          junk <- integer.base.b(i, b = 2)
+          if (length(junk) < length(y)) {
+            junk <- c(
+              rep(0, length(y) - length(junk)),
+              junk
+            )
+          }
+          k <- k + 1
+          permsy[k, ] <- y * (2 * junk - 1)
+        }
+      } else {
+        for (i in 1:nperms.act) {
+          permsy[i, ] <- sample(c(-1, 1),
+            size = length(y),
+            replace = TRUE
+          )
+        }
+      }
+    }
+    if (resp.type == samr.const.twoclass.paired.response) {
+      junk <- compute.block.perms(y, abs(y), nperms)
+      permsy <- junk$permsy
+      all.perms.flag <- junk$all.perms.flag
+      nperms.act <- junk$nperms.act
+    }
+    if (resp.type == samr.const.patterndiscovery.response) {
+      nperms.act <- nperms
+      perms <- NULL
+      permsy <- NULL
+      all.perms.flag <- FALSE
+    }
+    sdstar.keep <- NULL
+    if (assay.type != "seq") {
+      sdstar.keep <- matrix(0, ncol = nperms.act, nrow = nrow(x))
+    }
+    ttstar <- matrix(0, nrow = nrow(x), ncol = nperms.act)
+    foldchange.star <- NULL
+    if (resp.type == samr.const.twoclass.unpaired.response |
+      resp.type == samr.const.twoclass.paired.response) {
+      foldchange.star <- matrix(0, nrow = nrow(x), ncol = nperms.act)
+    }
+    if (resp.type == samr.const.multiclass.response) {
+      stand.contrasts.star <- array(NA, c(
+        nrow(x), length(table(y)),
+        nperms.act
+      ))
+    }
+  }
+  if (xl.mode == "next20" | xl.mode == "lasttime") {
+    evo <- xl.prevfit$evo
+    tt <- xl.prevfit$tt
+    numer <- xl.prevfit$numer
+    eigengene <- xl.prevfit$eigengene
+    eigengene.number <- xl.prevfit$eigengene.number
+    sd <- xl.prevfit$sd - xl.prevfit$s0
+    sd.internal <- xl.prevfit$sd.internal
+    ttstar <- xl.prevfit$ttstar
+    ttstar0 <- xl.prevfit$ttstar0
+    n <- xl.prevfit$n
+    pi0 <- xl.prevfit$pi0
+    foldchange <- xl.prevfit$foldchange
+    y <- xl.prevfit$y
+    x <- xl.prevfit$x
+    xresamp <- xl.prevfit$xresamp
+    censoring.status <- xl.prevfit$censoring.status
+    argy <- xl.prevfit$argy
+    testStatistic <- xl.prevfit$testStatistic
+    foldchange.star <- xl.prevfit$foldchange.star
+    s0 <- xl.prevfit$s0
+    s0.perc <- xl.prevfit$s0.perc
+    resp.type <- xl.prevfit$resp.type
+    resp.type.arg <- xl.prevfit$resp.type.arg
+    assay.type <- xl.prevfit$assay.type
+    sdstar.keep <- xl.prevfit$sdstar.keep
+    resp.type <- xl.prevfit$resp.type
+    stand.contrasts <- xl.prevfit$stand.contrasts
+    stand.contrasts.star <- xl.prevfit$stand.contrasts.star
+    stand.contrasts.95 <- xl.prevfit$stand.contrasts.95
+    perms <- xl.prevfit$perms
+    permsy <- xl.prevfit$permsy
+    nperms <- xl.prevfit$nperms
+    nperms.act <- xl.prevfit$nperms.act
+    all.perms.flag <- xl.prevfit$all.perms.flag
+    depth <- xl.prevfit$depth
+    nresamp <- xl.prevfit$nresamp
+    nresamp.perm <- xl.prevfit$nresamp.perm
+  }
+  if (xl.mode == "regular") {
+    first <- 1
+    last <- nperms.act
+  }
+  if (xl.mode == "firsttime") {
+    first <- 1
+    last <- 1
+  }
+  if (xl.mode == "next20") {
+    first <- xl.time
+    last <- min(xl.time + 19, nperms.act - 1)
+  }
+  if (xl.mode == "lasttime") {
+    first <- nperms.act
+    last <- nperms.act
+  }
+  for (b in first:last) {
+    message(c("perm = ", b))
+    if (assay.type == "array") {
+      xstar <- x
+    }
+    if (assay.type == "seq") {
+      xstar <- xresamp[, , 1:nresamp.perm]
+    }
+    if (resp.type == samr.const.oneclass.response) {
+      ystar <- permsy[b, ]
+      if (testStatistic == "standard") {
+        ttstar[, b] <- onesample.ttest.func(xstar, ystar,
+          s0 = s0, sd = sd.internal
+        )$tt
+      }
+    }
+    if (resp.type == samr.const.twoclass.paired.response) {
+      ystar <- permsy[b, ]
+      if (assay.type == "array") {
+        ttstar[, b] <- paired.ttest.func(xstar, ystar,
+          s0 = s0, sd = sd.internal
+        )$tt
+        foldchange.star[, b] <- foldchange.paired(
+          xstar,
+          ystar, data$logged2
+        )
+      }
+      if (assay.type == "seq") {
+        ttstar[, b] <- wilcoxon.paired.seq.func(
+          xstar,
+          ystar
+        )$tt
+        foldchange.star[, b] <- foldchange.seq.twoclass.paired(
+          x,
+          ystar, depth
+        )
+      }
+    }
+    if (resp.type == samr.const.twoclass.unpaired.response) {
+      ystar <- permsy[b, ]
+      if (assay.type == "array") {
+        if (testStatistic == "standard") {
+          junk <- ttest.func(xstar, ystar, s0 = s0, sd = sd.internal)
+        }
+        if (testStatistic == "wilcoxon") {
+          junk <- wilcoxon.func(xstar, ystar, s0 = s0)
+        }
+        ttstar[, b] <- junk$tt
+        sdstar.keep[, b] <- junk$sd
+        foldchange.star[, b] <- foldchange.twoclass(
+          xstar,
+          ystar, data$logged2
+        )
+      }
+      if (assay.type == "seq") {
+        ttstar[, b] <- wilcoxon.unpaired.seq.func(
+          xstar,
+          ystar
+        )$tt
+        foldchange.star[, b] <- foldchange.seq.twoclass.unpaired(
+          x,
+          ystar, depth
+        )
+      }
+    }
+    if (resp.type == samr.const.survival.response) {
+      o <- perms[b, ]
+      if (assay.type == "array") {
+        ttstar[, b] <- cox.func(
+          xstar, y[o],
+          censoring.status = censoring.status[o], s0 = s0
+        )$tt
+      }
+      if (assay.type == "seq") {
+        ttstar[, b] <- cox.seq.func(
+          xstar, y[o],
+          censoring.status = censoring.status[o]
+        )$tt
+      }
+    }
+    if (resp.type == samr.const.multiclass.response) {
+      ystar <- y[perms[b, ]]
+      if (assay.type == "array") {
+        junk <- multiclass.func(xstar, ystar, s0 = s0)
+        ttstar[, b] <- junk$tt
+        sdstar.keep[, b] <- junk$sd
+        stand.contrasts.star[, , b] <- junk$stand.contrasts
+      }
+      if (assay.type == "seq") {
+        junk <- multiclass.seq.func(xstar, ystar)
+        ttstar[, b] <- junk$tt
+        stand.contrasts.star[, , b] <- junk$stand.contrasts
+      }
+    }
+    if (resp.type == samr.const.quantitative.response) {
+      ystar <- y[perms[b, ]]
+      if (assay.type == "array") {
+        junk <- quantitative.func(xstar, ystar, s0 = s0)
+        ttstar[, b] <- junk$tt
+        sdstar.keep[, b] <- junk$sd
+      }
+      if (assay.type == "seq") {
+        junk <- quantitative.seq.func(xstar, ystar)
+        ttstar[, b] <- junk$tt
+      }
+    }
+    if (resp.type == samr.const.patterndiscovery.response) {
+      xstar <- permute.rows(x)
+      junk <- patterndiscovery.func(
+        xstar,
+        s0 = s0, eigengene.number = eigengene.number
+      )
+      ttstar[, b] <- junk$tt
+      sdstar.keep[, b] <- junk$sd
+    }
+  }
+  if (xl.mode == "regular" | xl.mode == "lasttime") {
+    ttstar0 <- ttstar
+    for (j in seq_len(ncol(ttstar))) {
+      ttstar[, j] <- -1 * sort(-1 * ttstar[, j])
+    }
+    for (i in seq_len(nrow(ttstar))) {
+      ttstar[i, ] <- sort(ttstar[i, ])
+    }
+    evo <- apply(ttstar, 1, mean)
+    evo <- evo[seq(length(evo), 1)]
+    pi0 <- 1
+    if (resp.type != samr.const.multiclass.response) {
+      qq <- quantile(ttstar, c(0.25, 0.75))
+    }
+    if (resp.type == samr.const.multiclass.response) {
+      qq <- quantile(ttstar, c(0, 0.5))
+    }
+    pi0 <- sum(tt > qq[1] & tt < qq[2]) / (0.5 * length(tt))
+    foldchange <- NULL
+    if (resp.type == samr.const.twoclass.unpaired.response &
+      assay.type == "array") {
+      foldchange <- foldchange.twoclass(x, y, data$logged2)
+    }
+    if (resp.type == samr.const.twoclass.paired.response &
+      assay.type == "array") {
+      foldchange <- foldchange.paired(x, y, data$logged2)
+    }
+    if (resp.type == samr.const.oneclass.response & assay.type ==
+      "array") {
+    }
+    stand.contrasts.95 <- NULL
+    if (resp.type == samr.const.multiclass.response) {
+      stand.contrasts.95 <- quantile(
+        stand.contrasts.star,
+        c(0.025, 0.975)
+      )
+    }
+    if (resp.type == samr.const.twoclass.unpaired.response &
+      assay.type == "seq") {
+      foldchange <- foldchange.seq.twoclass.unpaired(
+        x,
+        y, depth
+      )
+    }
+    if (resp.type == samr.const.twoclass.paired.response &
+      assay.type == "seq") {
+      foldchange <- foldchange.seq.twoclass.paired(
+        x, y,
+        depth
+      )
+    }
+    if (return.x == FALSE) {
+      x <- NULL
+    }
+  }
+  return(
+    list(
+      n = n, x = x, xresamp = xresamp, y = y, argy = argy,
+      censoring.status = censoring.status, testStatistic = testStatistic,
+      nperms = nperms, nperms.act = nperms.act, tt = tt, numer = numer,
+      sd = sd + s0, sd.internal = sd.internal, s0 = s0, s0.perc = s0.perc,
+      evo = evo, perms = perms, permsy = permsy, nresamp = nresamp,
+      nresamp.perm = nresamp.perm, all.perms.flag = all.perms.flag,
+      ttstar = ttstar, ttstar0 = ttstar0, eigengene = eigengene,
+      eigengene.number = eigengene.number, pi0 = pi0, foldchange = foldchange,
+      foldchange.star = foldchange.star, sdstar.keep = sdstar.keep,
+      resp.type = resp.type, resp.type.arg = resp.type.arg,
+      assay.type = assay.type, stand.contrasts = stand.contrasts,
+      stand.contrasts.star = stand.contrasts.star,
+      stand.contrasts.95 = stand.contrasts.95,
+      depth = depth, call = this.call
+    )
+  )
+}
+
+#' @title Estimate sequencing depths
+#' @param x data matrix. nrow=#gene, ncol=#sample
+#' @return depth: estimated sequencing depth. a vector with len sample.
+samr.estimate.depth <- function(x) {
+  iter <- 5
+  cmeans <- colSums(x) / sum(x)
+  for (i in 1:iter) {
+    n0 <- rowSums(x) %*% t(cmeans)
+    prop <- rowSums((x - n0)^2 / (n0 + 1e-08))
+    qs <- quantile(prop, c(0.25, 0.75))
+    keep <- (prop >= qs[1]) & (prop <= qs[2])
+    cmeans <- colMeans(x[keep, ])
+    cmeans <- cmeans / sum(cmeans)
+  }
+  depth <- cmeans / mean(cmeans)
+  return(depth)
+}
+
+#' @title Resampling
+#' @param x data matrix. nrow=#gene, ncol=#sample
+#' @param d estimated sequencing depth
+#' @param nresamp number of resamplings
+#' @return xresamp: an rank array with dim #gene*#sample*nresamp
+#' @description Corresponds to `samr::resample`
+#' @importFrom stats rpois runif
+resa <- function(x, d, nresamp = 20) {
+  ng <- nrow(x)
+  ns <- ncol(x)
+  dbar <- exp(mean(log(d)))
+  xresamp <- array(0, dim = c(ng, ns, nresamp))
+  for (k in 1:nresamp) {
+    for (j in 1:ns) {
+      xresamp[, j, k] <- rpois(n = ng, lambda = (dbar / d[j]) *
+        x[, j]) + runif(ng) * 0.1
+    }
+  }
+  for (k in 1:nresamp) {
+    xresamp[, , k] <- t(rankcols(t(xresamp[, , k])))
+  }
+  return(xresamp)
+}
+
+#' @title Rank columns
+#' @description Ranks the elements within each col of the matrix x and returns
+#' these ranks in a matrix
+#' @note this function is equivalent to `samr::rankcol`, but uses `apply` to
+#' rank the colums instead of a compiled Fortran function which was causing our
+#' DEGanalysis functions to freeze in large datasets.
+#' @param x x
+rankcols <- function(x) {
+  # ranks the elements within each col of the matrix x
+  # and returns these ranks in a matrix
+  n <- nrow(x)
+  p <- ncol(x)
+  mode(n) <- "integer"
+  mode(p) <- "integer"
+  mode(x) <- "double"
+  xr <- apply(x, 2, rank)
+  return(xr)
+}
+
+#' @title Check format
+#' @param y y
+#' @param resp.type resp type
+#' @param censoring.status censoring status
+check.format <- function(y, resp.type, censoring.status = NULL) {
+  # here i do some format checks for the input data$y
+  # note that checks for time course data are done in the
+  #   parse function for time course;
+  #  we then check the output from the parser in this function
+  if (resp.type == samr.const.twoclass.unpaired.response |
+    resp.type == samr.const.twoclass.unpaired.timecourse.response) {
+    if (sum(y == 1) + sum(y == 2) != length(y)) {
+      stop(paste(
+        "Error in input response data: response type ",
+        resp.type, " specified; values must be 1 or 2"
+      ))
+    }
+  }
+  if (resp.type == samr.const.twoclass.paired.response | resp.type ==
+    samr.const.twoclass.paired.timecourse.response) {
+    if (sum(y) != 0) {
+      stop(paste(
+        "Error in input response data: response type ",
+        resp.type, " specified; values must be -1, 1, -2, 2, etc"
+      ))
+    }
+    if (sum(table(y[y > 0]) != abs(table(y[y < 0])))) {
+      stop(paste(
+        "Error in input response data:  response type ",
+        resp.type, " specified; values must be -1, 1, -2, 2, etc"
+      ))
+    }
+  }
+  if (resp.type == samr.const.oneclass.response | resp.type ==
+    samr.const.oneclass.timecourse.response) {
+    if (sum(y == 1) != length(y)) {
+      stop(paste(
+        "Error in input response data: response type ",
+        resp.type, " specified;  values must all be 1"
+      ))
+    }
+  }
+  if (resp.type == samr.const.multiclass.response) {
+    tt <- table(y)
+    nc <- length(tt)
+    if (sum(y <= nc & y > 0) < length(y)) {
+      stop(
+        "Error in input response data: response type ", resp.type,
+        " specified; values must be 1,2, ... number of classes"
+      )
+    }
+    for (k in 1:nc) {
+      if (sum(y == k) < 2) {
+        stop(paste(
+          "Error in input response data: response type ",
+          resp.type, " specified; there must be >1 sample per class"
+        ))
+      }
+    }
+  }
+  if (resp.type == samr.const.quantitative.response) {
+    if (!is.numeric(y)) {
+      stop(paste(
+        "Error in input response data: response type",
+        resp.type, " specified; values must be numeric"
+      ))
+    }
+  }
+  if (resp.type == samr.const.survival.response) {
+    if (is.null(censoring.status)) {
+      stop(paste(
+        "Error in input response data: response type ",
+        resp.type, " specified; error in censoring indicator"
+      ))
+    }
+    if (!is.numeric(y) | sum(y < 0) > 0) {
+      stop(
+        "Error in input response data:  response type ", resp.type,
+        " specified; survival times  must be numeric and nonnegative"
+      )
+      if (sum(censoring.status == 0) + sum(censoring.status ==
+        1) != length(censoring.status)) {
+        stop(
+          "Error in input response data: response type ", resp.type,
+          " specified; censoring indicators must be ",
+          "0 (censored) or 1 (failed)"
+        )
+      }
+    }
+    if (sum(censoring.status == 1) < 1) {
+      stop(paste(
+        "Error in input response data:   response type ",
+        resp.type, " specified; there are no uncensored observations"
+      ))
+    }
+  }
+  return()
+}
+
+#' @title Twoclass Wilcoxon statistics
+#' @param xresamp an rank array with dim #gene*#sample*nresamp
+#' @param y outcome vector of values 1 and 2
+#' @return the statistic.
+wilcoxon.unpaired.seq.func <- function(xresamp, y) {
+  tt <- rep(0, dim(xresamp)[1])
+  for (i in seq_len(dim(xresamp)[3])) {
+    tt <- tt + rowSums(xresamp[, y == 2, i]) - sum(y == 2) *
+      (length(y) + 1) / 2
+  }
+  tt <- tt / dim(xresamp)[3]
+  return(list(tt = tt, numer = tt, sd = rep(1, length(tt))))
+}
+wilcoxon.paired.seq.func <- function(xresamp, y) {
+  tt <- rep(0, dim(xresamp)[1])
+  for (i in seq_len(dim(xresamp)[3])) {
+    tt <- tt + rowSums(xresamp[, y > 0, i]) - sum(y > 0) *
+      (length(y) + 1) / 2
+  }
+  tt <- tt / dim(xresamp)[3]
+  return(list(tt = tt, numer = tt, sd = rep(1, length(tt))))
+}
+getperms <- function(y, nperms) {
+  total.perms <- factorial(length(y))
+  if (total.perms <= nperms) {
+    perms <- permute(seq_len(length(y)))
+    all.perms.flag <- 1
+    nperms.act <- total.perms
+  }
+  if (total.perms > nperms) {
+    perms <- matrix(NA, nrow = nperms, ncol = length(y))
+    for (i in 1:nperms) {
+      perms[i, ] <- sample(seq_len(length(y)), size = length(y))
+    }
+    all.perms.flag <- 0
+    nperms.act <- nperms
+  }
+  return(list(
+    perms = perms, all.perms.flag = all.perms.flag,
+    nperms.act = nperms.act
+  ))
+}
+foldchange.twoclass <- function(x, y, logged2) {
+  m1 <- rowMeans(x[, y == 1, drop = FALSE])
+  m2 <- rowMeans(x[, y == 2, drop = FALSE])
+  if (!logged2) {
+    fc <- m2 / m1
+  }
+  if (logged2) {
+    fc <- 2^{
+      m2 - m1
+    }
+  }
+  return(fc)
+}
+#' @title Foldchange of twoclass unpaired sequencing data
+#' @param x x
+#' @param y y
+#' @param depth depth
+foldchange.seq.twoclass.unpaired <- function(x, y, depth) {
+  x.norm <- scale(x, center = FALSE, scale = depth) + 1e-08
+  fc <- apply(x.norm[, y == 2], 1, median) /
+    apply(x.norm[, y ==
+      1], 1, median)
+  return(fc)
+}
+foldchange.seq.twoclass.paired <- function(x, y, depth) {
+  nc <- ncol(x) / 2
+  o1 <- o2 <- rep(0, nc)
+  for (j in 1:nc) {
+    o1[j] <- which(y == -j)
+    o2[j] <- which(y == j)
+  }
+  x.norm <- scale(x, center = FALSE, scale = depth) + 1e-08
+  d <- x.norm[, o2, drop = FALSE] / x.norm[, o1, drop = FALSE]
+  fc <- lapply(d, 1, function(x) median(x, na.rm = TRUE))
+  return(fc)
+}
+permute <- function(elem) {
+  # generates all perms of the vector elem
+  if (!missing(elem)) {
+    if (length(elem) == 2) {
+      return(matrix(c(elem, elem[2], elem[1]), nrow = 2))
+    }
+    last.matrix <- permute(elem[-1])
+    dim.last <- dim(last.matrix)
+    new.matrix <- matrix(0, nrow = dim.last[1] * (dim.last[2] +
+      1), ncol = dim.last[2] + 1)
+    for (row in 1:(dim.last[1])) {
+      for (col in 1:(dim.last[2] + 1)) {
+        new.matrix[row + (col - 1) * dim.last[1], ] <-
+          insert.value(last.matrix[row, ], elem[1], col)
+      }
+    }
+    return(new.matrix)
+  } else {
+    message("Usage: permute(elem)\n\twhere elem is a vector")
+  }
+}
+insert.value <- function(vec, newval, pos) {
+  if (pos == 1) {
+    return(c(newval, vec))
+  }
+  lvec <- length(vec)
+  if (pos > lvec) {
+    return(c(vec, newval))
+  }
+  return(c(vec[1:pos - 1], newval, vec[pos:lvec]))
+}
+parse.block.labels.for.2classes <- function(y) {
+  # this only works for 2 class case- having form jBlockn,
+  #   where j=1 or 2
+  n <- length(y)
+  y.act <- rep(NA, n)
+  blocky <- rep(NA, n)
+  for (i in 1:n) {
+    blocky[i] <- as.numeric(substring(y[i], 7, nchar(y[i])))
+    y.act[i] <- as.numeric(substring(y[i], 1, 1))
+  }
+  return(list(y.act = y.act, blocky = blocky))
+}
+parse.time.labels.and.summarize.data <- function(
+    x,
+    y, resp.type, time.summary.type) {
+  # parse time labels, and summarize time data for each
+  #   person, via a slope or area
+  # does some error checking too
+  n <- length(y)
+  last5char <- rep(NA, n)
+  last3char <- rep(NA, n)
+  for (i in 1:n) {
+    last3char[i] <- substring(y[i], nchar(y[i]) - 2, nchar(y[i]))
+    last5char[i] <- substring(y[i], nchar(y[i]) - 4, nchar(y[i]))
+  }
+  if (sum(last3char == "End") != sum(last5char == "Start")) {
+    stop("Error in format of time course data: a Start or End tag is missing")
+  }
+  y.act <- rep(NA, n)
+  timey <- rep(NA, n)
+  person.id <- rep(NA, n)
+  k <- 1
+  end.flag <- FALSE
+  person.id[1] <- 1
+  if (substring(y[1], nchar(y[1]) - 4, nchar(y[1])) != "Start") {
+    stop(
+      "Error in format of time course data: ",
+      "first cell should have a Start tag"
+    )
+  }
+  for (i in 1:n) {
+    message(i)
+    j <- 1
+    while (substring(y[i], j, j) != "T") {
+      j <- j + 1
+    }
+    end.of.y <- j - 1
+    y.act[i] <- as.numeric(substring(y[i], 1, end.of.y))
+    timey[i] <- substring(y[i], end.of.y + 5, nchar(y[i]))
+    if (nchar(timey[i]) > 3 & substring(timey[i], nchar(timey[i]) -
+      2, nchar(timey[i])) == "End") {
+      end.flag <- TRUE
+      timey[i] <- substring(timey[i], 1, nchar(timey[i]) -
+        3)
+    }
+    if (nchar(timey[i]) > 3 & substring(timey[i], nchar(timey[i]) -
+      4, nchar(timey[i])) == "Start") {
+      timey[i] <- substring(timey[i], 1, nchar(timey[i]) -
+        5)
+    }
+    if (i < n & !end.flag) {
+      person.id[i + 1] <- k
+    }
+    if (i < n & end.flag) {
+      k <- k + 1
+      person.id[i + 1] <- k
+    }
+    end.flag <- FALSE
+  }
+  timey <- as.numeric(timey)
+  # do a check that the format was correct
+  tt <- table(person.id, y.act)
+  junk <- function(x) {
+    sum(x != 0)
+  }
+  if (sum(apply(tt, 1, junk) != 1) > 0) {
+    num <- seq_len(nrow(tt))[apply(tt, 1, junk) > 1]
+    stop(paste(
+      "Error in format of  time course data, timecourse #",
+      as.character(num)
+    ))
+  }
+  npeople <- length(unique(person.id))
+  newx <- matrix(NA, nrow = nrow(x), ncol = npeople)
+  sd <- matrix(NA, nrow = nrow(x), ncol = npeople)
+  for (j in 1:npeople) {
+    jj <- person.id == j
+    tim <- timey[jj]
+    xc <- t(scale(t(x[, jj, drop = FALSE]), center = TRUE, scale = FALSE))
+    if (time.summary.type == "slope") {
+      junk <- quantitative.func(xc, tim - mean(tim))
+      newx[, j] <- junk$numer
+      sd[, j] <- junk$sd
+    }
+    if (time.summary.type == "signed.area") {
+      junk <- timearea.func(x[, jj, drop = FALSE], tim)
+      newx[, j] <- junk$numer
+      sd[, j] <- junk$sd
+    }
+  }
+  y.unique <- y.act[!duplicated(person.id)]
+  return(list(y = y.unique, x = newx, sd = sd))
+}
+ttest.func <- function(x, y, s0 = 0, sd = NULL) {
+  n1 <- sum(y == 1)
+  n2 <- sum(y == 2)
+  m1 <- rowMeans(x[, y == 1, drop = FALSE])
+  m2 <- rowMeans(x[, y == 2, drop = FALSE])
+  if (is.null(sd)) {
+    sd <- sqrt(((n2 - 1) * varr(x[, y == 2], meanx = m2) +
+      (n1 - 1) * varr(x[, y == 1], meanx = m1)) * (1 / n1 +
+      1 / n2) / (n1 + n2 - 2))
+  }
+  numer <- m2 - m1
+  dif.obs <- (numer) / (sd + s0)
+  return(list(tt = dif.obs, numer = numer, sd = sd))
+}
+
+wilcoxon.func <- function(x, y, s0 = 0) {
+  n1 <- sum(y == 1)
+  n2 <- sum(y == 2)
+  p <- nrow(x)
+  r2 <- rowSums(t(apply(x, 1, rank))[, y == 2, drop = FALSE])
+  numer <- r2 - (n2 / 2) * (n2 + 1) - (n1 * n2) / 2
+  sd <- sqrt(n1 * n2 * (n1 + n2 + 1) / 12)
+  tt <- (numer) / (sd + s0)
+  return(list(tt = tt, numer = numer, sd = rep(sd, p)))
+}
+
+onesample.ttest.func <- function(x, y, s0 = 0, sd = NULL) {
+  n <- length(y)
+  x <- x * matrix(y, nrow = nrow(x), ncol = ncol(x), byrow = TRUE)
+  m <- rowMeans(x)
+  if (is.null(sd)) {
+    sd <- sqrt(varr(x, meanx = m) / n)
+  }
+  dif.obs <- m / (sd + s0)
+  return(list(tt = dif.obs, numer = m, sd = sd))
+}
+
+patterndiscovery.func <- function(x, s0 = 0, eigengene.number = 1) {
+  a <- mysvd(x, n.components = eigengene.number)
+  v <- a$v[, eigengene.number]
+  # here we try to guess the most interpretable orientation
+  #   for the eigengene
+  om <- abs(a$u[, eigengene.number]) > quantile(
+    abs(a$u[, eigengene.number]),
+    0.95
+  )
+  if (median(a$u[om, eigengene.number]) < 0) {
+    v <- -1 * v
+  }
+  aa <- quantitative.func(x, v, s0 = s0)
+  eigengene <- cbind(seq_len(nrow(a$v)), v)
+  dimnames(eigengene) <- list(NULL, c("sample number", "value"))
+  return(
+    list(tt = aa$tt, numer = aa$numer, sd = aa$sd, eigengene = eigengene)
+  )
+}
+
+paired.ttest.func <- function(x, y, s0 = 0, sd = NULL) {
+  nc <- ncol(x) / 2
+  o <- 1:nc
+  o1 <- rep(0, ncol(x) / 2)
+  o2 <- o1
+  for (j in 1:nc) {
+    o1[j] <- seq_len(ncol(x))[y == -o[j]]
+  }
+  for (j in 1:nc) {
+    o2[j] <- seq_len(ncol(x))[y == o[j]]
+  }
+  d <- x[, o2, drop = FALSE] - x[, o1, drop = FALSE]
+  if (is.matrix(d)) {
+    m <- rowMeans(d)
+  }
+  if (!is.matrix(d)) {
+    m <- mean(d)
+  }
+  if (is.null(sd)) {
+    if (is.matrix(d)) {
+      sd <- sqrt(varr(d, meanx = m) / nc)
+    }
+    if (!is.matrix(d)) {
+      sd <- sqrt(var(d) / nc)
+    }
+  }
+  dif.obs <- m / (sd + s0)
+  return(list(tt = dif.obs, numer = m, sd = sd))
+}
+
+cox.func <- function(x, y, censoring.status, s0 = 0) {
+  # find the index matrix
+  Dn <- sum(censoring.status == 1)
+  Dset <- seq_len(ncol(x))[censoring.status == 1] # the set of observed
+  ind <- matrix(0, ncol(x), Dn)
+  # get the matrix
+  for (i in 1:Dn) {
+    ind[y > y[Dset[i]] - 1e-08, i] <- 1 / sum(y > y[Dset[i]] -
+      1e-08)
+  }
+  ind.sums <- rowSums(ind)
+  x.ind <- x %*% ind
+  # get the derivatives
+  numer <- x %*% (censoring.status - ind.sums)
+  sd <- sqrt((x * x) %*% ind.sums - rowSums(x.ind * x.ind))
+  tt <- numer / (sd + s0)
+  return(list(tt = tt, numer = numer, sd = sd))
+}
+
+multiclass.func <- function(x, y, s0 = 0) {
+  ## assumes y is coded 1,2...
+  nn <- table(y)
+  m <- matrix(0, nrow = nrow(x), ncol = length(nn))
+  v <- m
+  for (j in seq_len(length(nn))) {
+    m[, j] <- rowMeans(x[, y == j])
+    v[, j] <- (nn[j] - 1) * varr(x[, y == j], meanx = m[
+      ,
+      j
+    ])
+  }
+  mbar <- rowMeans(x)
+  mm <- m - matrix(mbar, nrow = length(mbar), ncol = length(nn))
+  fac <- (sum(nn) / prod(nn))
+  scor <- sqrt(fac * (apply(matrix(nn,
+    nrow = nrow(m), ncol = ncol(m),
+    byrow = TRUE
+  ) * mm * mm, 1, sum)))
+  sd <- sqrt(rowSums(v) * (1 / sum(nn - 1)) * sum(1 / nn))
+  tt <- scor / (sd + s0)
+  mm.stand <- t(scale(t(mm), center = FALSE, scale = sd))
+  return(list(tt = tt, numer = scor, sd = sd, stand.contrasts = mm.stand))
+}
+
+est.s0 <- function(tt, sd, s0.perc = seq(0, 1, by = 0.05)) {
+  ## estimate s0 (exchangeability) factor for denominator.
+  ## returns the actual estimate s0 (not a percentile)
+  br <- unique(quantile(sd, seq(0, 1, len = 101)))
+  nbr <- length(br)
+  a <- cut(sd, br, labels = FALSE)
+  a[is.na(a)] <- 1
+  cv.sd <- rep(0, length(s0.perc))
+  for (j in seq_len(length(s0.perc))) {
+    w <- quantile(sd, s0.perc[j])
+    w[j == 1] <- 0
+    tt2 <- tt * sd / (sd + w)
+    tt2[tt2 == Inf] <- NA
+    sds <- rep(0, nbr - 1)
+    for (i in 1:(nbr - 1)) {
+      sds[i] <- stats::mad(tt2[a == i], na.rm = TRUE)
+    }
+    cv.sd[j] <- sqrt(var(sds)) / mean(sds)
+  }
+  o <- seq_len(length(s0.perc))[cv.sd == min(cv.sd)]
+  # we don;t allow taking s0.hat to be 0th percentile when
+  #   min sd is 0
+  s0.hat <- quantile(sd[sd != 0], s0.perc[o])
+  return(list(s0.perc = s0.perc, cv.sd = cv.sd, s0.hat = s0.hat))
+}
+
+permute.rows <- function(x) {
+  dd <- dim(x)
+  n <- dd[1]
+  p <- dd[2]
+  mm <- runif(length(x)) + rep(seq(n) * 10, rep(p, n))
+  matrix(t(x)[order(mm)], n, p, byrow = TRUE)
+}
+
+foldchange.paired <- function(x, y, logged2) {
+  nc <- ncol(x) / 2
+  o <- 1:nc
+  o1 <- rep(0, ncol(x) / 2)
+  o2 <- o1
+  for (j in 1:nc) {
+    o1[j] <- seq_len(ncol(x))[y == -o[j]]
+  }
+  for (j in 1:nc) {
+    o2[j] <- seq_len(ncol(x))[y == o[j]]
+  }
+  if (!logged2) {
+    d <- x[, o2, drop = FALSE] / x[, o1, drop = FALSE]
+  }
+  if (logged2) {
+    d <- x[, o2, drop = FALSE] - x[, o1, drop = FALSE]
+  }
+  if (!logged2) {
+    fc <- rowMeans(d)
+  }
+  if (logged2) {
+    fc <- 2^rowMeans(d)
+  }
+  return(fc)
+}
+foldchange.seq.twoclass.unpaired <- function(x, y, depth) {
+  x.norm <- scale(x, center = FALSE, scale = depth) + 1e-08
+  fc <- apply(x.norm[, y == 2], 1, median) /
+    apply(x.norm[, y == 1], 1, median)
+  return(fc)
+}
+integer.base.b <- function(x, b = 2) {
+  xi <- as.integer(x)
+  if (xi == 0) {
+    return(0)
+  }
+  if (any(is.na(xi) | ((x - xi) != 0))) {
+    print(list(ERROR = "x not integer", x = x))
+  }
+  N <- length(x)
+  xMax <- max(x)
+  ndigits <- (floor(logb(xMax, base = 2)) + 1)
+  Base.b <- array(NA, dim = c(N, ndigits))
+  for (i in 1:ndigits) {
+    Base.b[, ndigits - i + 1] <- (x %% b)
+    x <- (x %/% b)
+  }
+  if (N == 1) {
+    Base.b[1, ]
+  } else {
+    Base.b
+  }
+}
+varr <- function(x, meanx = NULL) {
+  n <- ncol(x)
+  p <- nrow(x)
+  Y <- matrix(1, nrow = n, ncol = 1)
+  if (is.null(meanx)) {
+    meanx <- rowMeans(x)
+  }
+  ans <- rep(1, p)
+  xdif <- x - meanx %*% t(Y)
+  ans <- (xdif^2) %*% rep(1 / (n - 1), n)
+  ans <- drop(ans)
+  return(ans)
+}
+quantitative.func <- function(x, y, s0 = 0) {
+  # regression of x on y
+  my <- mean(y)
+  yy <- y - my
+  temp <- x %*% yy
+  mx <- rowMeans(x)
+  syy <- sum(yy^2)
+  scor <- temp / syy
+  b0hat <- mx - scor * my
+  ym <- matrix(y, nrow = nrow(x), ncol = ncol(x), byrow = TRUE)
+  xhat <- matrix(b0hat, nrow = nrow(x), ncol = ncol(x)) + ym *
+    matrix(scor, nrow = nrow(x), ncol = ncol(x))
+  sigma <- sqrt(rowSums((x - xhat)^2) / (ncol(xhat) - 2))
+  sd <- sigma / sqrt(syy)
+  tt <- scor / (sd + s0)
+  return(list(tt = tt, numer = scor, sd = sd))
+}
+timearea.func <- function(x, y, s0 = 0) {
+  n <- ncol(x)
+  xx <- 0.5 * (x[, 2:n] + x[, 1:(n - 1)]) * matrix(diff(y),
+    nrow = nrow(x), ncol = n - 1, byrow = TRUE
+  )
+  numer <- rowMeans(xx)
+  sd <- sqrt(varr(xx, meanx = numer) / n)
+  tt <- numer / sqrt(sd + s0)
+  return(list(tt = tt, numer = numer, sd = sd))
+}
+cox.seq.func <- function(xresamp, y, censoring.status) {
+  # get the dimensions
+  ns <- ncol(xresamp)
+  # prepare for the calculation
+  # find the index matrix
+  Dn <- sum(censoring.status == 1)
+  Dset <- seq_len(ns)[censoring.status == 1] # the set of died
+  ind <- matrix(0, ns, Dn)
+  # get the matrix
+  for (i in 1:Dn) {
+    ind[y >= y[Dset[i]] - 1e-08, i] <- 1 / sum(y >= y[Dset[i]] -
+      1e-08)
+  }
+  ind.sums <- rowSums(ind)
+  # calculate the score statistic
+  tt <- apply(xresamp, 3, function(x, cen.ind, ind.para, ind.sums.para) {
+    dev1 <- x %*% cen.ind
+    x.ind <- x %*% ind.para
+    dev2 <- (x * x) %*% ind.sums.para - rowSums(x.ind * x.ind)
+    dev1 / (sqrt(dev2) + 1e-08)
+  }, (censoring.status - ind.sums), ind, ind.sums)
+  tt <- rowMeans(tt)
+  return(list(tt = tt, numer = tt, sd = rep(1, length(tt))))
+}
+compute.block.perms <- function(y, blocky, nperms) {
+  # y are the data (eg class label 1 vs 2; or -1,1, -2,2 for
+  #   paired data)
+  # blocky are the block labels (abs(y) for paired daatr)
+  ny <- length(y)
+  nblocks <- length(unique(blocky))
+  tab <- table(blocky)
+  total.nperms <- prod(factorial(tab))
+  # block.perms is a list of all possible permutations
+  block.perms <- vector("list", nblocks)
+  # first enumerate all perms, when possible
+  if (total.nperms <= nperms) {
+    all.perms.flag <- 1
+    nperms.act <- total.nperms
+    for (i in 1:nblocks) {
+      block.perms[[i]] <- permute(y[blocky == i])
+    }
+    kk <- 0:(factorial(max(tab))^nblocks - 1)
+    # the rows of the matrix outerm runs through the 'outer
+    #   product'
+    # first we assume that all blocks have max(tab) members;
+    #   then we remove rows of outerm that
+    #  are illegal (ie when a block has fewer members)
+    outerm <- matrix(0, nrow = length(kk), ncol = nblocks)
+    for (i in seq_len(length(kk))) {
+      kkkk <- integer.base.b(kk[i], b = factorial(max(tab)))
+      if (length(kkkk) > nblocks) {
+        kkkk <- kkkk[(length(kkkk) - nblocks + 1):length(kkkk)]
+      }
+      outerm[i, (nblocks - length(kkkk) + 1):nblocks] <- kkkk
+    }
+    outerm <- outerm + 1
+    # now remove rows that are illegal perms
+    ind <- rep(TRUE, nrow(outerm))
+    for (j in seq_len(ncol(outerm))) {
+      ind <- ind & outerm[, j] <= factorial(tab[j])
+    }
+    outerm <- outerm[ind, , drop = FALSE]
+    # finally, construct permutation matrix from outer product
+    permsy <- matrix(NA, nrow = total.nperms, ncol = ny)
+    for (i in 1:total.nperms) {
+      junk <- NULL
+      for (j in 1:nblocks) {
+        junk <- c(junk, block.perms[[j]][outerm[i, j], ])
+      }
+      permsy[i, ] <- junk
+    }
+  }
+  # next handle case when there are too many perms to enumerate
+  if (total.nperms > nperms) {
+    all.perms.flag <- 0
+    nperms.act <- nperms
+    permsy <- NULL
+    block.perms <- vector("list", nblocks)
+    for (j in 1:nblocks) {
+      block.perms[[j]] <- sample.perms(y[blocky == j], nperms = nperms)
+    }
+    for (j in 1:nblocks) {
+      permsy <- cbind(permsy, block.perms[[j]])
+    }
+  }
+  return(list(
+    permsy = permsy, all.perms.flag = all.perms.flag,
+    nperms.act = nperms.act
+  ))
+}
+sample.perms <- function(elem, nperms) {
+  # randomly generates  nperms of the vector elem
+  res <- permute.rows(matrix(elem,
+    nrow = nperms, ncol = length(elem),
+    byrow = TRUE
+  ))
+  return(res)
+}
+mysvd <- function(x, n.components = NULL) {
+  # finds PCs of matrix x
+  p <- nrow(x)
+  n <- ncol(x)
+  # center the observations (rows)
+  feature.means <- rowMeans(x)
+  x <- t(scale(t(x), center = feature.means, scale = FALSE))
+  if (is.null(n.components)) {
+    n.components <- min(n, p)
+  }
+  if (p > n) {
+    a <- eigen(t(x) %*% x)
+    v <- a$vec[, 1:n.components, drop = FALSE]
+    d <- sqrt(a$val[1:n.components, drop = FALSE])
+    u <- scale(x %*% v, center = FALSE, scale = d)
+    return(list(u = u, d = d, v = v))
+  } else {
+    junk <- svd(x, LINPACK = TRUE)
+    nc <- min(ncol(junk$u), n.components)
+    return(list(u = junk$u[, 1:nc], d = junk$d[1:nc], v = junk$v[
+      ,
+      1:nc
+    ]))
+  }
+}
+quantitative.seq.func <- function(xresamp, y) {
+  tt <- rep(0, dim(xresamp)[1])
+  for (i in seq_len(dim(xresamp)[3])) {
+    y.ranked <- rank(y, ties.method = "random") - (dim(xresamp)[2] +
+      1) / 2
+    tt <- tt + (xresamp[, , i] - (dim(xresamp)[2] + 1) / 2) %*%
+      y.ranked
+  }
+  ns <- dim(xresamp)[2]
+  tt <- tt / (dim(xresamp)[3] * (ns^3 - ns) / 12)
+  return(list(tt = as.vector(tt), numer = as.vector(tt), sd = rep(
+    1,
+    length(tt)
+  )))
+}
+multiclass.seq.func <- function(xresamp, y) {
+  # number of classes and number of samples in each class
+  K <- max(y)
+  n.each <- rep(0, K)
+  for (k in 1:K) {
+    n.each[k] <- sum(y == k)
+  }
+  # the statistic
+  tt <- temp <- rep(0, dim(xresamp)[1])
+  stand.contrasts <- matrix(0, dim(xresamp)[1], K)
+
+  for (i in seq_len(dim(xresamp)[3])) {
+    for (k in 1:K) {
+      temp <- rowSums(xresamp[, y == k, i])
+      tt <- tt + temp^2 / n.each[k]
+      stand.contrasts[, k] <- stand.contrasts[, k] + temp
+    }
+  }
+  # finalize
+  nresamp <- dim(xresamp)[3]
+  ns <- dim(xresamp)[2]
+  tt <- tt / nresamp * 12 / ns / (ns + 1) - 3 * (ns + 1)
+  stand.contrasts <- stand.contrasts / nresamp
+  stand.contrasts <- scale(stand.contrasts,
+    center = n.each * (ns + 1) / 2,
+    scale = sqrt(n.each * (ns - n.each) * (ns + 1) / 12)
+  )
+  return(list(
+    tt = tt, numer = tt, sd = rep(1, length(tt)),
+    stand.contrasts = stand.contrasts
+  ))
+}
+
+# ==============================================================================
+# samr.compute.delta.table
+# ==============================================================================
+## Jun added starts
+samr.compute.delta.table <- function(
+    samr.obj, min.foldchange = 0,
+    dels = NULL, nvals = 50) {
+  res <- NULL
+  if (samr.obj$assay.type == "array") {
+    res <- samr.compute.delta.table.array(
+      samr.obj, min.foldchange,
+      dels, nvals
+    )
+  } else if (samr.obj$assay.type == "seq") {
+    res <- samr.compute.delta.table.seq(
+      samr.obj, min.foldchange,
+      dels
+    )
+  }
+  return(res)
+}
+## Jun added ends
+
+## Jun added the first row below, and commented the row
+#   after it
+samr.compute.delta.table.array <- function(
+    samr.obj,
+    min.foldchange = 0, dels = NULL, nvals = 50) {
+  # samr.compute.delta.table <- function(samr.obj,
+  #   min.foldchange=0, dels=NULL, nvals=50) {
+  # computes delta table, starting with samr object 'a', for
+  #   nvals values of delta
+  lmax <- sqrt(max(abs(sort(samr.obj$tt) - samr.obj$evo)))
+  if (is.null(dels)) {
+    dels <- (seq(0, lmax, length = nvals)^2)
+  }
+  res1 <- NULL
+  foldchange.cond.up <- matrix(TRUE,
+    nrow = nrow(samr.obj$ttstar),
+    ncol = ncol(samr.obj$ttstar)
+  )
+  foldchange.cond.lo <- matrix(TRUE,
+    nrow = nrow(samr.obj$ttstar),
+    ncol = ncol(samr.obj$ttstar)
+  )
+  if (!is.null(samr.obj$foldchange[1]) & (min.foldchange >
+    0)) {
+    foldchange.cond.up <- samr.obj$foldchange.star >= min.foldchange
+    foldchange.cond.lo <- samr.obj$foldchange.star <= 1 / min.foldchange
+  }
+  cutup <- rep(NA, length(dels))
+  cutlow <- rep(NA, length(dels))
+  g2 <- rep(NA, length(dels))
+  errup <- matrix(NA, ncol = length(dels), nrow = ncol(samr.obj$ttstar0))
+  errlow <- matrix(NA, ncol = length(dels), nrow = ncol(samr.obj$ttstar0))
+  cat("", fill = TRUE)
+  cat("Computing delta table", fill = TRUE)
+  for (ii in seq_len(length(dels))) {
+    cat(ii, fill = TRUE)
+    ttt <- detec.slab(samr.obj, dels[ii], min.foldchange)
+    cutup[ii] <- 1e+10
+    if (length(ttt$pup > 0)) {
+      cutup[ii] <- min(samr.obj$tt[ttt$pup])
+    }
+    cutlow[ii] <- -1e+10
+    if (length(ttt$plow) > 0) {
+      cutlow[ii] <- max(samr.obj$tt[ttt$plow])
+    }
+    g2[ii] <- sumlengths(ttt)
+    errup[, ii] <- colSums(samr.obj$ttstar0 > cutup[ii] &
+      foldchange.cond.up)
+    errlow[, ii] <- colSums(samr.obj$ttstar0 < cutlow[ii] &
+      foldchange.cond.lo)
+  }
+  gmed <- apply(errup + errlow, 2, median)
+  g90 <- apply(errup + errlow, 2, quantile, 0.9)
+  res1 <- cbind(
+    samr.obj$pi0 * gmed, samr.obj$pi0 * g90, g2,
+    samr.obj$pi0 * gmed / g2, samr.obj$pi0 * g90 / g2, cutlow,
+    cutup
+  )
+  res1 <- cbind(dels, res1)
+  dimnames(res1) <- list(NULL, c(
+    "delta", "# med false pos",
+    "90th perc false pos", "# called", "median FDR", "90th perc FDR",
+    "cutlo", "cuthi"
+  ))
+  return(res1)
+}
+
+#######################################################################
+# \tcompute the delta table for sequencing data
+#######################################################################
+samr.compute.delta.table.seq <- function(
+    samr.obj,
+    min.foldchange = 0, dels = NULL) {
+  res1 <- NULL
+  flag <- TRUE
+  ## check whether any gene satisfies the foldchange
+  #   restrictions
+  if ((samr.obj$resp.type == samr.const.twoclass.unpaired.response |
+    samr.obj$resp.type == samr.const.twoclass.paired.response) &
+    (min.foldchange > 0)) {
+    sat.up <- (samr.obj$foldchange >= min.foldchange) & (samr.obj$evo >
+      0)
+    sat.dn <- (samr.obj$foldchange <= 1 / min.foldchange) &
+      (samr.obj$evo < 0)
+    if (sum(sat.up) + sum(sat.dn) == 0) {
+      flag <- FALSE
+    }
+  }
+  if (flag) {
+    if (is.null(dels)) {
+      dels <- generate.dels(samr.obj, min.foldchange = min.foldchange)
+    }
+    cat("Number of thresholds chosen (all possible thresholds) =",
+      length(dels),
+      fill = TRUE
+    )
+    if (length(dels) > 0) {
+      ## sort delta to make the fast calculation right
+      dels <- sort(dels)
+      ## get the upper and lower cutoffs
+      cat("Getting all the cutoffs for the thresholds...\n")
+      slabs <- samr.seq.detec.slabs(samr.obj, dels, min.foldchange)
+      cutup <- slabs$cutup
+      cutlow <- slabs$cutlow
+      g2 <- slabs$g2
+      ## get the number of errors under the null hypothesis
+      cat("Getting number of false positives in the permutation...\n")
+      errnum <- samr.seq.null.err(
+        samr.obj, min.foldchange,
+        cutup, cutlow
+      )
+      res1 <- NULL
+      gmed <- apply(errnum, 2, median)
+      g90 <- apply(errnum, 2, quantile, 0.9)
+      res1 <- cbind(samr.obj$pi0 * gmed, samr.obj$pi0 *
+        g90, g2, samr.obj$pi0 * gmed / g2, samr.obj$pi0 *
+        g90 / g2, cutlow, cutup)
+      res1 <- cbind(dels, res1)
+      dimnames(res1) <- list(NULL, c(
+        "delta", "# med false pos",
+        "90th perc false pos", "# called", "median FDR",
+        "90th perc FDR", "cutlo", "cuthi"
+      ))
+    }
+  }
+  return(res1)
+}
+
+# ==============================================================================
+# samr.plot
+# ==============================================================================
+samr.plot <- function(samr.obj, del = NULL, min.foldchange = 0) {
+  ## make observed-expected plot
+  ## takes foldchange into account too
+  if (is.null(del)) {
+    del <- sqrt(max(abs(sort(samr.obj$tt) - samr.obj$evo)))
+  }
+  b <- detec.slab(samr.obj, del, min.foldchange)
+  foldchange.cond.up <- rep(TRUE, length(samr.obj$evo))
+  foldchange.cond.lo <- rep(TRUE, length(samr.obj$evo))
+  if (!is.null(samr.obj$foldchange[1]) & (min.foldchange >
+    0)) {
+    foldchange.cond.up <- samr.obj$foldchange >= min.foldchange
+    foldchange.cond.lo <- samr.obj$foldchange <= 1 / min.foldchange
+  }
+  col <- rep(1, length(samr.obj$evo))
+  col[b$plow] <- 3
+  col[b$pup] <- 2
+  if (!is.null(samr.obj$foldchange[1]) & (min.foldchange >
+    0)) {
+    col[!foldchange.cond.lo & !foldchange.cond.up] <- 1
+  }
+  col.ordered <- col[order(samr.obj$tt)]
+  ylims <- range(samr.obj$tt)
+  xlims <- range(samr.obj$evo)
+  plot(samr.obj$evo, sort(samr.obj$tt),
+    xlab = "expected score",
+    ylab = "observed score", ylim = ylims, xlim = xlims,
+    type = "n"
+  )
+  points(samr.obj$evo, sort(samr.obj$tt), col = col.ordered)
+  abline(0, 1)
+  abline(del, 1, lty = 2)
+  abline(-del, 1, lty = 2)
+}
+
+# ==============================================================================
+# samr.compute.siggenes.table
+# ==============================================================================
+samr.compute.siggenes.table <- function(
+    samr.obj, del,
+    data, delta.table, min.foldchange = 0, all.genes = FALSE,
+    compute.localfdr = FALSE) {
+  ## computes significant genes table, starting with samr
+  #   object 'a' and 'delta.table'
+  ##  for a  **single** value del
+  ## if all.genes is true, all genes are printed (and value
+  #   of del is ignored)
+  if (is.null(data$geneid)) {
+    data$geneid <- paste("g", seq_len(nrow(data$x)), sep = "")
+  }
+  if (is.null(data$genenames)) {
+    data$genenames <- paste("g", seq_len(nrow(data$x)), sep = "")
+  }
+  if (!all.genes) {
+    sig <- detec.slab(samr.obj, del, min.foldchange)
+  }
+  if (all.genes) {
+    p <- length(samr.obj$tt)
+    pup <- (1:p)[samr.obj$tt >= 0]
+    plo <- (1:p)[samr.obj$tt < 0]
+    sig <- list(pup = pup, plo = plo)
+  }
+  if (compute.localfdr) {
+    aa <- localfdr(samr.obj, min.foldchange)
+    if (length(sig$pup) > 0) {
+      fdr.up <- predictlocalfdr(aa$smooth.object, samr.obj$tt[sig$pup])
+    }
+    if (length(sig$plo) > 0) {
+      fdr.lo <- predictlocalfdr(aa$smooth.object, samr.obj$tt[sig$plo])
+    }
+  }
+  qvalues <- NULL
+  if (length(sig$pup) > 0 | length(sig$plo) > 0) {
+    qvalues <- qvalue.func(samr.obj, sig, delta.table)
+  }
+  res.up <- NULL
+  res.lo <- NULL
+  done <- FALSE
+
+  # two class unpaired or paired  (foldchange is reported)
+  if ((samr.obj$resp.type == samr.const.twoclass.unpaired.response |
+    samr.obj$resp.type == samr.const.twoclass.paired.response)) {
+    if (!is.null(sig$pup)) {
+      res.up <- cbind(
+        sig$pup + 1, data$genenames[sig$pup],
+        data$geneid[sig$pup], samr.obj$tt[sig$pup], samr.obj$numer[sig$pup],
+        samr.obj$sd[sig$pup], samr.obj$foldchange[sig$pup],
+        qvalues$qvalue.up
+      )
+      if (compute.localfdr) {
+        res.up <- cbind(res.up, fdr.up)
+      }
+      temp.names <- list(NULL, c(
+        "Row", "Gene ID", "Gene Name",
+        "Score(d)", "Numerator(r)", "Denominator(s+s0)",
+        "Fold Change", "q-value(%)"
+      ))
+      if (compute.localfdr) {
+        temp.names[[2]] <- c(temp.names[[2]], "localfdr(%)")
+      }
+      dimnames(res.up) <- temp.names
+    }
+    if (!is.null(sig$plo)) {
+      res.lo <- cbind(
+        sig$plo + 1, data$genenames[sig$plo],
+        data$geneid[sig$plo], samr.obj$tt[sig$plo], samr.obj$numer[sig$plo],
+        samr.obj$sd[sig$plo], samr.obj$foldchange[sig$plo],
+        qvalues$qvalue.lo
+      )
+      if (compute.localfdr) {
+        res.lo <- cbind(res.lo, fdr.lo)
+      }
+      temp.names <- list(NULL, c(
+        "Row", "Gene ID", "Gene Name",
+        "Score(d)", "Numerator(r)", "Denominator(s+s0)",
+        "Fold Change", "q-value(%)"
+      ))
+      if (compute.localfdr) {
+        temp.names[[2]] <- c(temp.names[[2]], "localfdr(%)")
+      }
+      dimnames(res.lo) <- temp.names
+    }
+    done <- TRUE
+  }
+
+  # multiclass
+  if (samr.obj$resp.type == samr.const.multiclass.response) {
+    if (!is.null(sig$pup)) {
+      res.up <- cbind(
+        sig$pup + 1, data$genenames[sig$pup],
+        data$geneid[sig$pup], samr.obj$tt[sig$pup], samr.obj$numer[sig$pup],
+        samr.obj$sd[sig$pup], samr.obj$stand.contrasts[sig$pup, ],
+        qvalues$qvalue.up
+      )
+
+      if (compute.localfdr) {
+        res.up <- cbind(res.up, fdr.up)
+      }
+
+      collabs.contrast <- paste(
+        "contrast-", as.character(seq_len(ncol(samr.obj$stand.contrasts))),
+        sep = ""
+      )
+      temp.names <- list(NULL, c(
+        "Row", "Gene ID", "Gene Name",
+        "Score(d)", "Numerator(r)", "Denominator(s+s0)",
+        collabs.contrast, "q-value(%)"
+      ))
+
+      if (compute.localfdr) {
+        temp.names[[2]] <- c(temp.names[[2]], "localfdr(%)")
+      }
+      dimnames(res.up) <- temp.names
+    }
+    res.lo <- NULL
+    done <- TRUE
+  }
+
+  # all other cases
+  if (!done) {
+    if (!is.null(sig$pup)) {
+      res.up <- cbind(
+        sig$pup + 1, data$genenames[sig$pup],
+        data$geneid[sig$pup], samr.obj$tt[sig$pup], samr.obj$numer[sig$pup],
+        samr.obj$sd[sig$pup], samr.obj$foldchange[sig$pup],
+        qvalues$qvalue.up
+      )
+      if (compute.localfdr) {
+        res.up <- cbind(res.up, fdr.up)
+      }
+      temp.names <- list(NULL, c(
+        "Row", "Gene ID", "Gene Name",
+        "Score(d)", "Numerator(r)", "Denominator(s+s0)",
+        "q-value(%)"
+      ))
+      if (compute.localfdr) {
+        temp.names[[2]] <- c(temp.names[[2]], "localfdr(%)")
+      }
+      dimnames(res.up) <- temp.names
+    }
+    if (!is.null(sig$plo)) {
+      res.lo <- cbind(
+        sig$plo + 1, data$genenames[sig$plo],
+        data$geneid[sig$plo], samr.obj$tt[sig$plo], samr.obj$numer[sig$plo],
+        samr.obj$sd[sig$plo], samr.obj$foldchange[sig$plo],
+        qvalues$qvalue.lo
+      )
+      if (compute.localfdr) {
+        res.lo <- cbind(res.lo, fdr.lo)
+      }
+      temp.names <- list(NULL, c(
+        "Row", "Gene ID", "Gene Name",
+        "Score(d)", "Numerator(r)", "Denominator(s+s0)",
+        "q-value(%)"
+      ))
+      if (compute.localfdr) {
+        temp.names[[2]] <- c(temp.names[[2]], "localfdr(%)")
+      }
+      dimnames(res.lo) <- temp.names
+    }
+    done <- TRUE
+  }
+  if (!is.null(res.up)) {
+    o1 <- order(-samr.obj$tt[sig$pup])
+    res.up <- res.up[o1, , drop = FALSE]
+  }
+  if (!is.null(res.lo)) {
+    o2 <- order(samr.obj$tt[sig$plo])
+    res.lo <- res.lo[o2, , drop = FALSE]
+  }
+  color.ind.for.multi <- NULL
+  if (
+    samr.obj$resp.type == samr.const.multiclass.response & !is.null(sig$pup)
+  ) {
+    condition_1 <-
+      samr.obj$stand.contrasts[sig$pup, ] > samr.obj$stand.contrasts.95[2]
+    condition_2 <-
+      samr.obj$stand.contrasts[sig$pup, ] < samr.obj$stand.contrasts.95[1]
+    color.ind.for.multi <- 1 * condition_1 + (-1) * condition_2
+  }
+  ngenes.up <- nrow(res.up)
+  if (is.null(ngenes.up)) {
+    ngenes.up <- 0
+  }
+  ngenes.lo <- nrow(res.lo)
+  if (is.null(ngenes.lo)) {
+    ngenes.lo <- 0
+  }
+  return(
+    list(
+      genes.up = res.up, genes.lo = res.lo,
+      color.ind.for.multi = color.ind.for.multi, ngenes.up = ngenes.up,
+      ngenes.lo = ngenes.lo
+    )
+  )
+}
+generate.dels <- function(samr.obj, min.foldchange = 0) {
+  dels <- NULL
+  ## initialize calculation
+  tag <- order(samr.obj$tt)
+  if ((samr.obj$resp.type == samr.const.twoclass.unpaired.response |
+    samr.obj$resp.type == samr.const.twoclass.paired.response) &
+    (min.foldchange > 0)) {
+    res.mat <- data.frame(
+      tt = samr.obj$tt[tag], fc = samr.obj$foldchange[tag],
+      evo = samr.obj$evo, dif = samr.obj$tt[tag] - samr.obj$evo
+    )
+    res.up <- res.mat[res.mat$evo > 0, ]
+    res.lo <- res.mat[res.mat$evo < 0, ]
+    res.up <- res.up[res.up$fc >= min.foldchange, ]
+    res.lo <- res.lo[res.lo$fc <= 1 / min.foldchange, ]
+  } else {
+    res.mat <- data.frame(
+      tt = samr.obj$tt[tag], evo = samr.obj$evo,
+      dif = samr.obj$tt[tag] - samr.obj$evo
+    )
+    res.up <- res.mat[res.mat$evo > 0, ]
+    res.lo <- res.mat[res.mat$evo < 0, ]
+  }
+  ## for the upper part
+  up.vec <- rep(NA, nrow(res.up))
+  if (nrow(res.up) > 0) {
+    st <- 1e-08
+    i.cur <- 1
+    for (i in seq_len(nrow(res.up))) {
+      if (res.up$dif[i] > st) {
+        st <- res.up$dif[i]
+        up.vec[i.cur] <- st
+        i.cur <- i.cur + 1
+      }
+    }
+  }
+  ## for the lower part
+  lo.vec <- rep(NA, nrow(res.lo))
+  if (nrow(res.lo) > 0) {
+    st <- -1e-08
+    i.cur <- 1
+    for (i in seq(nrow(res.lo), 1)) {
+      if (res.lo$dif[i] < st) {
+        st <- res.lo$dif[i]
+        lo.vec[i.cur] <- st
+        i.cur <- i.cur + 1
+      }
+    }
+  }
+  ## combine them
+  vec <- c(up.vec, -lo.vec)
+  vec <- vec[!is.na(vec)]
+  vec <- vec - 1e-08
+  dels <- sort(unique(vec))
+  return(dels)
+}
+samr.seq.detec.slabs <- function(samr.obj, dels, min.foldchange) {
+  ## initialize calculation
+  tag <- order(samr.obj$tt)
+  if ((samr.obj$resp.type == samr.const.twoclass.unpaired.response |
+    samr.obj$resp.type == samr.const.twoclass.paired.response) &
+    (min.foldchange > 0)) {
+    res.mat <- data.frame(
+      tt = samr.obj$tt[tag], fc = samr.obj$foldchange[tag],
+      evo = samr.obj$evo, dif = samr.obj$tt[tag] - samr.obj$evo
+    )
+    res.up <- res.mat[res.mat$evo > 0, ]
+    res.lo <- res.mat[res.mat$evo < 0, ]
+    res.up <- res.up[res.up$fc >= min.foldchange, ]
+    res.lo <- res.lo[res.lo$fc <= 1 / min.foldchange, ]
+  } else {
+    res.mat <- data.frame(
+      tt = samr.obj$tt[tag], evo = samr.obj$evo,
+      dif = samr.obj$tt[tag] - samr.obj$evo
+    )
+    res.up <- res.mat[res.mat$evo > 0, ]
+    res.lo <- res.mat[res.mat$evo < 0, ]
+  }
+  ## begin calculating
+  cutup <- rep(1e+10, length(dels))
+  cutlow <- rep(-1e+10, length(dels))
+  g2.up <- g2.lo <- rep(0, length(dels))
+  if (nrow(res.up) > 0) {
+    res.up <- data.frame(
+      dif = res.up$dif, tt = res.up$tt,
+      num = seq(nrow(res.up), 1)
+    )
+    ## get the upper part
+    j <- 1
+    ii <- 1
+    while (j <= nrow(res.up) & ii <= length(dels)) {
+      if (res.up$dif[j] > dels[ii]) {
+        cutup[ii] <- res.up$tt[j]
+        g2.up[ii] <- res.up$num[j]
+        ii <- ii + 1
+      } else {
+        j <- j + 1
+      }
+    }
+  }
+  if (nrow(res.lo) > 0) {
+    res.lo <- data.frame(
+      dif = res.lo$dif, tt = res.lo$tt,
+      num = seq_len(nrow(res.lo))
+    )
+    ## get the lower part
+    j <- nrow(res.lo)
+    ii <- 1
+    while (j >= 1 & ii <= length(dels)) {
+      if (res.lo$dif[j] < -dels[ii]) {
+        cutlow[ii] <- res.lo$tt[j]
+        g2.lo[ii] <- res.lo$num[j]
+        ii <- ii + 1
+      } else {
+        j <- j - 1
+      }
+    }
+  }
+  g2 <- g2.up + g2.lo
+  return(list(cutup = cutup, cutlow = cutlow, g2 = g2))
+}
+sumlengths <- function(aa) {
+  length(aa$pl) + length(aa$pu)
+}
+
+samr.seq.null.err <- function(
+    samr.obj, min.foldchange,
+    cutup, cutlow) {
+  errup <- matrix(NA, ncol = length(cutup), nrow = ncol(samr.obj$ttstar0))
+  errlow <- matrix(NA, ncol = length(cutlow), nrow = ncol(samr.obj$ttstar0))
+  cutup.rank <- rank(cutup, ties.method = "min")
+  cutlow.rank <- rank(-cutlow, ties.method = "min")
+  for (jj in seq_len(ncol(samr.obj$ttstar0))) {
+    keep.up <- keep.dn <- samr.obj$ttstar0[, jj]
+    if ((samr.obj$resp.type == samr.const.twoclass.unpaired.response |
+      samr.obj$resp.type == samr.const.twoclass.paired.response) &
+      (min.foldchange > 0)) {
+      keep.up <- keep.up[samr.obj$foldchange.star[, jj] >=
+        min.foldchange]
+      keep.dn <- keep.dn[samr.obj$foldchange.star[, jj] <=
+        1 / min.foldchange]
+    }
+    errup[jj, ] <- length(keep.up) - (rank(c(cutup, keep.up),
+      ties.method = "min"
+    )[seq_len(length(cutup))] - cutup.rank)
+    errlow[jj, ] <- length(keep.dn) - (rank(c(-cutlow, -keep.dn),
+      ties.method = "min"
+    )[seq_len(length(cutlow))] - cutlow.rank)
+  }
+  errnum <- errup + errlow
+  return(errnum)
+}
+detec.slab <- function(samr.obj, del, min.foldchange) {
+  ## find genes above and below the slab of half-width del
+  # this calculation is tricky- for consistency, the slab
+  #   condition picks
+  # all genes that are beyond the first departure from the
+  #   slab
+  # then the fold change condition is applied (if applicable)
+  n <- length(samr.obj$tt)
+  tt <- samr.obj$tt
+  evo <- samr.obj$evo
+  tag <- order(tt)
+  pup <- NULL
+  foldchange.cond.up <- rep(TRUE, length(evo))
+  foldchange.cond.lo <- rep(TRUE, length(evo))
+  if (!is.null(samr.obj$foldchange[1]) & (min.foldchange >
+    0)) {
+    foldchange.cond.up <- samr.obj$foldchange >= min.foldchange
+    foldchange.cond.lo <- samr.obj$foldchange <= 1 / min.foldchange
+  }
+  o1 <- (1:n)[(tt[tag] - evo > del) & evo > 0]
+  if (length(o1) > 0) {
+    o1 <- o1[1]
+    o11 <- o1:n
+    o111 <- rep(FALSE, n)
+    o111[tag][o11] <- TRUE
+    pup <- (1:n)[o111 & foldchange.cond.up]
+  }
+  plow <- NULL
+  o2 <- (1:n)[(evo - tt[tag] > del) & evo < 0]
+  if (length(o2) > 0) {
+    o2 <- o2[length(o2)]
+    o22 <- 1:o2
+    o222 <- rep(FALSE, n)
+    o222[tag][o22] <- TRUE
+    plow <- (1:n)[o222 & foldchange.cond.lo]
+  }
+  return(list(plow = plow, pup = pup))
+}
+
+#' @importFrom stats smooth.spline
+localfdr <- function(
+    samr.obj, min.foldchange, perc = 0.01,
+    df = 10) {
+  ## estimates compute.localfdr at score 'd', using SAM
+  #   object 'samr.obj'
+  ## 'd' can be a vector of d scores
+  ## returns estimate of symmetric fdr  as a percentage
+  # this version uses a 1% symmetric window, and does not
+  #   estimate fdr in
+  # windows  having fewer than 100 genes
+  ## to use: first run samr and then pass the resulting fit
+  #   object to
+  ## localfdr
+  ## NOTE: at most 20 of the perms are used to estimate the
+  #   fdr (for speed sake)
+  # I try two window shapes: symmetric and an assymetric one
+  # currently I use the symmetric window to estimate the
+  #   compute.localfdr
+  ngenes <- length(samr.obj$tt)
+  mingenes <- 50
+  # perc is increased, in order to get at least mingenes in a
+  #   window
+  perc <- max(perc, mingenes / length(samr.obj$tt))
+  nperms.to.use <- min(20, ncol(samr.obj$ttstar))
+  nperms <- ncol(samr.obj$ttstar)
+  d <- seq(sort(samr.obj$tt)[1], sort(samr.obj$tt)[ngenes],
+    length = 100
+  )
+  ndscore <- length(d)
+  dvector <- rep(NA, ndscore)
+  ind.foldchange <- rep(TRUE, length(samr.obj$tt))
+  if (!is.null(samr.obj$foldchange[1]) & min.foldchange > 0) {
+    ind.foldchange <- (samr.obj$foldchange >= min.foldchange) |
+      (samr.obj$foldchange <= min.foldchange)
+  }
+  fdr.temp <- function(temp, dlow, dup, pi0, ind.foldchange) {
+    return(sum(pi0 * (temp >= dlow & temp <= dup & ind.foldchange)))
+  }
+  for (i in 1:ndscore) {
+    pi0 <- samr.obj$pi0
+    r <- sum(samr.obj$tt < d[i])
+    r22 <- round(max(r - length(samr.obj$tt) * perc / 2, 1))
+    dlow.sym <- sort(samr.obj$tt)[r22]
+    r22 <- min(r + length(samr.obj$tt) * perc / 2, length(samr.obj$tt))
+    dup.sym <- sort(samr.obj$tt)[r22]
+    oo <- samr.obj$tt >= dlow.sym & samr.obj$tt <= dup.sym &
+      ind.foldchange
+    if (!is.null(samr.obj$foldchange[1]) & min.foldchange >
+      0) {
+      temp <- as.vector(samr.obj$foldchange.star[, 1:nperms.to.use])
+      ind.foldchange <- (temp >= min.foldchange) | (temp <=
+        min.foldchange)
+    }
+    temp <- samr.obj$ttstar0[, sample(1:nperms, size = nperms.to.use)]
+    fdr.sym <- median(apply(
+      temp, 2, fdr.temp, dlow.sym,
+      dup.sym, pi0, ind.foldchange
+    ))
+    fdr.sym <- 100 * fdr.sym / sum(oo)
+    dlow.sym <- dlow.sym
+    dup.sym <- dup.sym
+    dvector[i] <- fdr.sym
+  }
+  om <- !is.na(dvector) & (dvector != Inf)
+  aa <- smooth.spline(d[om], dvector[om], df = df)
+  return(list(smooth.object = aa, perc = perc, df = df))
+}
+
+predictlocalfdr <- function(smooth.object, d) {
+  yhat <- predict(smooth.object, d)$y
+  yhat <- pmin(yhat, 100)
+  yhat <- pmax(yhat, 0)
+  return(yhat)
+}
+
+qvalue.func <- function(samr.obj, sig, delta.table) {
+  # returns q-value as a percentage (out of 100)
+  LARGE <- 1e+10
+  qvalue.up <- rep(NA, length(sig$pup))
+  o1 <- sig$pup
+  cutup <- delta.table[, 8]
+  FDR <- delta.table[, 5]
+  ii <- 0
+  for (i in o1) {
+    o <- abs(cutup - samr.obj$tt[i])
+    o[is.na(o)] <- LARGE
+    oo <- seq_len(length(o))[o == min(o)]
+    oo <- oo[length(oo)]
+    ii <- ii + 1
+    qvalue.up[ii] <- FDR[oo]
+  }
+  qvalue.lo <- rep(NA, length(sig$plo))
+  o2 <- sig$plo
+  cutlo <- delta.table[, 7]
+  ii <- 0
+  for (i in o2) {
+    o <- abs(cutlo - samr.obj$tt[i])
+    o[is.na(o)] <- LARGE
+    oo <- seq_len(length(o))[o == min(o)]
+    oo <- oo[length(oo)]
+    ii <- ii + 1
+    qvalue.lo[ii] <- FDR[oo]
+  }
+  # any qvalues that are missing, are set to 1 (the highest
+  #   value)
+  qvalue.lo[is.na(qvalue.lo)] <- 1
+  qvalue.up[is.na(qvalue.up)] <- 1
+  # ensure that each qvalue vector is monotone non-increasing
+  o1 <- order(samr.obj$tt[sig$plo])
+  qv1 <- qvalue.lo[o1]
+  qv11 <- qv1
+  if (length(qv1) > 1) {
+    for (i in 2:length(qv1)) {
+      if (qv11[i] < qv11[i - 1]) {
+        qv11[i] <- qv11[i - 1]
+      }
+    }
+    qv111 <- qv11
+    qv111[o1] <- qv11
+  } else {
+    qv111 <- qv1
+  }
+  o2 <- order(samr.obj$tt[sig$pup])
+  qv2 <- qvalue.up[o2]
+  qv22 <- qv2
+  if (length(qv2) > 1) {
+    for (i in 2:length(qv2)) {
+      if (qv22[i] > qv22[i - 1]) {
+        qv22[i] <- qv22[i - 1]
+      }
+    }
+    qv222 <- qv22
+    qv222[o2] <- qv22
+  } else {
+    qv222 <- qv2
+  }
+  return(list(qvalue.lo = 100 * qv111, qvalue.up = 100 * qv222))
+}