An important preliminary to any experiment is assessing the power of the experimental design to optimize statistical sensitivity. In essence, a power analysis can tell us the number of replicates required to determine an effect size of a given magnitude for a given amount of experimental variability.

We'll use the powsimR package, which is not part of Bioconductor, to perform two types of power analysis. Both of these will be with a small real dataset, but first, we'll do it with two treatments—a test and control—then, we'll do it with just one. With each, we'll estimate the number of replicates we need to spot differences in gene expression of a particular magnitude—if they're present. powsimR takes a simulation-based approach, effectively generating many datasets and evaluating the detection power in each to create a distribution of detection power. The first step, then, is to estimate some parameters for these simulations—for this, we'll need some sample or preliminary data. After that, we can run simulations and assess power.

The dataset for this recipe will be a test or control RNAseq experiment from Arabidopsis with three replicates each. These are available as a prepared count matrix in datasets/ch1/arabidopsis.RDS in this book's data repository. In this section, we'll use a set of counts in a simple test or control experiment from Arabidopsis thaliana. The matrix has six columns (three mock treatments and three hrcc treatments) and 26,222 rows, each a gene feature. We'll need the dplyr, extRemes, and powsimR packages for this code.

Our package of interest, powsimR, isn't on CRAN; it's hosted as a source on GitHub at https://github.com/bvieth/powsimR. You'll need to use devtools to install it, which can be done using the following code:

install.packages("devtools")
devtools::install_github("bvieth/powsimR")

If you do this, there is a chance that this package will still fail to install. It has a lot of dependencies and you might need to install those manually; there is further information on the package GitHub repository and you should check that for the latest information. At the time of writing, you'll need to do the following two big steps. First, create the ipak function outlined here, then run the three different package installation steps with the ipak function:

ipak <- function(pkg, repository = c("CRAN", "Bioconductor", "github")) {
    new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
    # new.pkg <- pkg
    if (length(new.pkg)) {
        if (repository == "CRAN") {
            install.packages(new.pkg, dependencies = TRUE)
        }
        if (repository == "Bioconductor") {
            if (strsplit(version[["version.string"]], " ")[[1]][3] > "3.5.0") {
                if (!requireNamespace("BiocManager")) {
                  install.packages("BiocManager")
                }
                BiocManager::install(new.pkg, dependencies = TRUE, ask = FALSE)
            }
            if (strsplit(version[["version.string"]], " ")[[1]][3] < "3.5.0") {
                source("https://bioconductor.org/biocLite.R")
                biocLite(new.pkg, dependencies = TRUE, ask = FALSE)
            }
        }
        if (repository == "github") {
            devtools::install_github(new.pkg, build_vignettes = FALSE, force = FALSE, 
                dependencies = TRUE)
        }
    }
}

# CRAN PACKAGES
cranpackages <- c("broom", "cobs", "cowplot", "data.table", "devtools", "doParallel", 
    "dplyr", "drc", "DrImpute", "fastICA", "fitdistrplus", "foreach", "gamlss.dist", 
    "ggExtra", "ggplot2", "ggthemes", "grDevices", "glmnet", "grid", "gtools", 
    "Hmisc", "kernlab", "MASS", "MBESS", "matrixStats", "mclust", "methods", 
    "minpack.lm", "moments", "msir", "NBPSeq", "nonnest2", "parallel", "penalized", 
    "plyr", "pscl", "reshape2", "Rmagic", "rsvd", "Rtsne", "scales", "Seurat", 
    "snow", "stats", "tibble", "tidyr", "VGAM", "ZIM")

ipak(cranpackages, repository = "CRAN")

# BIOCONDUCTOR
biocpackages <- c("AnnotationDbi", "bayNorm", "baySeq", "Biobase", "BiocGenerics", 
    "BiocParallel", "DEDS", "DESeq2", "EBSeq", "edgeR", "IHW", "iCOBRA", "limma", 
    "Linnorm", "MAST", "monocle", "NOISeq", "qvalue", "ROTS", "RUVSeq", "S4Vectors", 
    "scater", "scDD", "scde", "scone", "scran", "SCnorm", "SingleCellExperiment", 
    "SummarizedExperiment", "zinbwave")
ipak(biocpackages, repository = "Bioconductor")

# GITHUB
githubpackages <- c("nghiavtr/BPSC", "cz-ye/DECENT", "mohuangx/SAVER", "statOmics/zingeR")
ipak(githubpackages, repository = "github")

When this is done, you should be able to install the package we're after with this code:

devtools::install_github("bvieth/powsimR", build_vignettes = TRUE, dependencies = FALSE)
library("powsimR")

We will do the power analysis using the following steps:

1. Estimate simulation parameter values:

arab_data <- readRDS(file.path(getwd(), "datasets", "ch1", "arabidopsis.RDS" ))
means_mock <- rowMeans(arab_data[, c("mock1", "mock2", "mock3")])
means_hrcc <- rowMeans(arab_data[, c("hrcc1", "hrcc2", "hrcc3")])
log2fc <- log2(means_hrcc / means_mock)
prop_de <- sum(abs(log2fc) > 2) / length(log2fc) 

2. Examine the distribution of the log2 fold change ratios:

finite_log2fc <-log2fc[is.finite(log2fc)]
plot(density(finite_log2fc))
extRemes::qqnorm(finite_log2fc )

3. Set up parameter values for the simulation run:

library(powsimR)
 library(dplyr)
 
 params <- estimateParam(
 countData = arab_data,
 Distribution = "NB",
 RNAseq = "bulk",
 normalization = "TMM" # edgeR method, can be others
 )
 
 de_opts <- DESetup(ngenes=1000,
 nsims=25,
 p.DE = prop_de,
 pLFC= finite_log2fc,
 sim.seed = 58673
 )
 
 sim_opts <- SimSetup(
 desetup = de_opts,
 params = params
 )
 
 num_replicates <- c(2, 3, 5, 8, 12,15)

4. Run the simulation:

 simDE <- simulateDE(n1 = num_replicates,
 n2 = num_replicates,
 sim.settings = sim_opts,
 DEmethod = "edgeR-LRT",
 normalization = "TMM",
 verbose = FALSE)

5. Run the evaluation of the simulation:

 evalDE <- evaluateDE(simRes = simDE,
 alpha.type = 'adjusted',
 MTC = 'BH',
 alpha.nominal = 0.1,
 stratify.by = 'mean',
 filter.by = 'none',
 strata.filtered = 1,
 target.by = 'lfc',
 delta = 0)

6. Plot the evaluation:

 plotEvalDE(evalRes = evalDE,
 rate='marginal',
 quick=FALSE, annot=TRUE)

Power analysis in powsimR requires us to do some pre-analysis so that we have estimates for some important parameters. To perform a simulation-based power analysis, we need to estimate the distribution of log fold changes between treatments and the proportion of features that are differentially expressed.

In step 1, we'll get the mean counts for each feature in the two treatments. After loading the expression data using the readRDS() function, we use the rowMeans() function on certain columns to get the mean expression counts of each gene in both the mock and hrcc1 treatments. We can then get the log2 ratio of those (by simply dividing the two vectors and, in the last line, use standard arithmetical operators to work out those that have a log2 fold change greater than 2). Inspecting the final prop_de variable gives the following output:

prop_de
## [1] 0.2001754

So, a proportion of about 0.2 of the features have counts changing by log2 twofold.

Step 2 looks at the distribution of the gene expression ratios. We first remove the non-finite ratios from the log2fc variable. We must do this because, when calculating ratios, we generate Inf values in R; this occurs when the denominator (the mock sample) has zero mean counts. We can remove them using indexing on the vector with the binary vector that comes from is.finite() function. With the Inf values removed, we can plot. First, we do a normal density plot using the density() function, which shows the distribution of ratios. Then, we use the qqnorm() function in the extRemes package, which plots the data against data sampled from an idealized normal distribution with the same mean. A strong, linear correlation indicates a normal distribution in the original data. We can see the output in the following screenshot:

They look pretty log-normally distributed, so we can assume a log-normal distribution.

The longest step here, step 3, is actually only four lines. We are basically setting up the parameters for the simulation, which requires us to specify a lot of values. The first set, params, which we create with the estimateParam() function needs the data source (countData), the distribution to use (we set Distribution = "NB", which selects the negative binomial); the type of RNAseq experiment—ours is a bulk RNAseq experiment (RNAseq = "bulk"), and normalization strategy—we use the edgeR style TMM (normalization = "TMM"). The second set, desetup, is created with the DESetup() function; in this, we choose the parameters relating to the number of genes for which to simulate differential expression. We set up 1,000 total gene simulations (ngenes) and 25 simulation runs (nsims). We set the proportion to be differentially expressed to that estimated in step 1 in prop_de. We use the vector of fold changes, finite_log2fc, as input for the pLFC parameter. Setting sim.seed is not necessary but will ensure reproducibility between runs. The third line uses the SimSetup() function to combine params and desetup into a single object, sim_opts. Finally, we create a num_replicates vector specifying the number of biological replicates (RNA samples) to simulate.

Step 4 is relatively straightforward: we run the differential expression simulation using the sim_opts parameters created in the previous steps, choosing "edgeR-LRT" as the differential expression method and "TMM" as the normalization. The simulation data is stored in the simDE variable.

In Step 5, we create an evaluation of the simulation—this analyzes and extracts various statistics. We pass the simDE simulation data to the evaluateDE() function along with values for things pertaining to grouping, filtering, and significance.

Finally, in Step 6, we can plot the evalDE object from Step 5 and see the results of the simulation. We get the following plot in which we can see the different powers at different replicate numbers. Note the x-axis indicates the number of replicate RNA samples used, and the metrics include FDR, False Negative/Positive Rate (FNR/FPR), and TNR/TPR (True Negative/Positive Rate):

When we have only one sample (or maybe even just one replicate), we have a hard time estimating the log2 fold change distribution and the number of differentially expressed genes. In place of estimates, we can use a callback function to generate numbers as needed. The body of the function just needs to return numbers from a specified distribution with parameters you decide. Here, we'll build a function that returns numbers with a normal distribution of mean 0 and standard deviation 2. This reflects that we think the log fold change distribution is normal with these parameters. When we've built the function, it gets used in the DESetup() function in place of the vector of log2 fold changes. For the proportion of genes differentially expressed, we just have to guess or take an estimate from something we already know about the experimental system:

log2fc_func <- function(x){ rnorm(x, 0, 2)} 
prop_de = 0.1
de_opts <- DESetup(ngenes=1000,
 nsims=25,
 p.DE = prop_de,
 pLFC= log2fc_func,
 sim.seed = 58673
 )