R/pauc_out.R In rmRNAseq: RNA-Seq Analysis for Repeated-Measures Data

Documented in pauc_out

```#' Evaluation of  Differential Expression Analysis Methods
#'
#'This function calculates FDR (false discovery rate), pauc (partial area
#'under ROC curve),  when we know the vector of
#'p-values obtained from a particular differential expression analysis method and the
#'true status of each gene. The function requires that the first EE genes are
#'true EE, and the last DE genes are true DE. This requirement can be fulfilled
#'by reorder the rows of gene expression data set.
#'
#'@param p a vector of p-values.
#'@param EE number of EE genes (the first EE genes in p-value vector p).
#'@param DE number of DE genes (the last DE genes in p-value vector p).
#'@return a vector including  V (the
#'  number of false positives), R (the number of declared positives), FDR (the
#'  false discovery rate), pauc (the partial area under ROC curve with respect
#'  to false positive rate fpr less than or equal to a specified level), and
#'  auc.
#'  Here we consider 3 fpr: 0.05, 0.10, 0.20. So the output includes these 15 elements
#'  and the total auc.
#'
#' @examples
#' set.seed(1)
#' EE <- 1000
#' DE <- 500
#' p1 <- runif(EE)
#' p2 <- rbeta(DE, shape1 = .5, shape2 = 1)
#' p <- c(p1, p2)
#' rmRNAseq:::pauc_out(p, EE, DE)

pauc_out <- function(p, EE, DE) {
lab <- as.factor(c(rep(0, EE), rep(1, DE)))
roc.out <- AUC::roc(1 - p, lab)

fpr <- 0.05
roc.ind <- sum(roc.out\$fpr <= fpr)
roc.min <- roc.out\$cutoffs[roc.ind]
pauc <- AUC::auc(roc.out, min = roc.min)
if (length(pauc) == 0)
pauc <- 0
qv <- jabes.q(p)
R <- sum(qv <= fpr)
V <- sum(which(qv <= fpr) <= EE)
FDR <- V/max(1, R)
R.05 <- R; V.05 <- V; FDR.05 <- FDR; pauc.05 <- pauc

fpr <- 0.10
roc.ind <- sum(roc.out\$fpr <= fpr)
roc.min <- roc.out\$cutoffs[roc.ind]
pauc <- AUC::auc(roc.out, min = roc.min)
if (length(pauc) == 0)
pauc <- 0
qv <- jabes.q(p)
R <- sum(qv <= fpr)
V <- sum(which(qv <= fpr) <= EE)
FDR <- V/max(1, R)
R.10 <- R; V.10 <- V; FDR.10 <- FDR; pauc.10 <- pauc

fpr <- 0.20
roc.ind <- sum(roc.out\$fpr <= fpr)
roc.min <- roc.out\$cutoffs[roc.ind]
pauc <- AUC::auc(roc.out, min = roc.min)
if (length(pauc) == 0)
pauc <- 0
qv <- jabes.q(p)
R <- sum(qv <= fpr)
V <- sum(which(qv <= fpr) <= EE)
FDR <- V/max(1, R)
R.20 <- R; V.20 <- V; FDR.20 <- FDR; pauc.20 <- pauc

auc <- AUC::auc(roc.out)
out <- c(R.05 = R.05, V.05 = V.05, FDR.05 = FDR.05, pauc.05 = pauc.05,
R.10 = R.10, V.10 = V.10, FDR.10 = FDR.10, pauc.10 = pauc.10,
R.20 = R.20, V.20 = V.20, FDR.20 = FDR.20, pauc.20 = pauc.20,
auc =auc)
return(out)
}
```

Try the rmRNAseq package in your browser

Any scripts or data that you put into this service are public.

rmRNAseq documentation built on Nov. 8, 2019, 5:06 p.m.