Datasets
Models
Diff of
/R/internal-functions-samr-adapted.R
[000000]
..
[28e211]
Switch to side-by-side view
--- 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)) +}