stageR example code

stager_example_code.R

ngene <- 1000
sig_ratio <- .1
niso <- rpois(ngene, lambda=5)
niso <- pmax(niso, 2)

set.seed(1)

gene <- factor(rep(1:ngene,niso))
pval <- runif(length(gene))
idx <- sapply(split(seq_along(gene), gene)[1:(ngene*sig_ratio)], `[`, 1)

# spike in small p-values for sig_ratio of genes
pval[idx] <- runif(length(idx), 1e-6, 1e-5)

library(dplyr)
library(tibble)

# apply some p-value aggregation method, here simple `min`
dat <- tibble(gene, pval)
screen <- dat %>%
  group_by(gene) %>%
  summarize(agg_pval = min(pval)) %>%
  deframe()

suppressPackageStartupMessages(library(stageR))

# note: tx2gene cannot have factors for columns
stageRObj <- stageRTx(
  pScreen=screen,
  pConfirmation=matrix(dat$pval,ncol=1,dimnames=list(seq_len(nrow(dat)))),
  tx2gene=data.frame(feature_id=seq_along(pval), gene_id=dat$gene)
)

stageRObj <- stageWiseAdjustment(stageRObj, method="dte", alpha=.05)

res <- getAdjustedPValues(stageRObj, order=FALSE, onlySignificantGenes=FALSE)

# some plots to see what happened
library(ggplot2)

res %>%
  as_tibble() %>%
  mutate(geneID = factor(geneID, levels=seq_len(ngene))) %>%
  group_by(geneID) %>%
  dplyr::slice(1) %>%
  ungroup() %>%
  mutate(geneID = as.numeric(geneID)) %>%
  ggplot(aes(geneID, -log10(gene))) +
  geom_point()

res %>%
  as_tibble() %>%
  mutate_at(c("geneID", "txID"), as.numeric) %>%
  filter(geneID < 20) %>%
  mutate(status = ifelse(txID %in% idx, "spike", "null")) %>%
  ggplot(aes(txID, -log10(transcript), color=status)) +
  geom_point()