Introduction

EpipwR is a quasi-simulation based approach to power analysis for epigenome-wide association studies (EWAS) with a continuous or binary outcome variable. It is fast, efficient and allows the exact level of precision to be set by the user. EpipwR is the first open source tool able to calculate power for EWAS with a continuous outcome. For binary outcomes, it re-implements the pwrEWAS methodology (Graw et al., 2019 - currently available on GitHub) with a few changes and significantly improves the computation time. EpipwR is simple in that it has three main functions: get_power_cont() and get_power_cc calculate power for various settings with continuous and binary outcomes (respectively) and EpipwR_plot() plots the output. The following vignette describes in detail the EpipwR workflow for an analysis of a continuous outcome.

Methodology (Continuous Outcome)

Suppose an upcoming EWAS plans to analyze the methylation levels at \(K\) CpG sites and some continuous outcome \(Y\). Of the \(K\) CpG sites, it is assumed that only a subest of these (\(K_m\) ) have non-zero correlation to \(Y\). The goal of EpipwR is to determine the level of power achieved by this study at a given sample size \(n\), where power is defined as the proportion of the \(K_m\) CpG sites are found to have a statistically significant association with \(Y\).

To begin, EpipwR generates data for the \(K_m\) CpG methylation levels and for \(Y\). Users specify \(\mu_\rho\) and \(\sigma_\rho\), the average and standard deviation of the correlation between each of the \(K_m\) methylation levels and \(Y\). If users prefer that all \(K_m\) methylation levels have equal correlation with \(Y\), then \(\sigma_\rho\) can be set to 0. Like many power analysis tools, EpipwR assumes no correlation across CpG sites (i.e., they are conditionally independent given \(Y\)). Users also select which empirical EWAS data to use to generate initial methylation levels. Full details of the data generation process can be found in the corresponding EpipwR article (Barth and Reynolds, 2024).

Once a dataset has been generated, p-values are calculated by performing the specified hypothesis test (Pearson, Spearman or Kendall correlation) for each \(K_m\) site. Please note that in absence of any “vertical” correlation or additional covariates, performing a standard analysis in limma is equivalent to using the Pearson correlation test. Users can specify either a false discovery rate or family-wise type I error rate correction for determining statistical significance. If the false discovery rate is selected, p-values for the remaining null tests are generated via a \(Uniform(0,1)\) distribution.

Since the process described above is based heavily on random sampling, it must be repeated on several data sets to produce a reliable power estimate. Whereas traditional tools require the user to specify the number of repetitions, EpipwR sets 20 to the minimum number of data sets, allows users to set a maximum number of data sets and a level of precision in the form of a 95% confidence interval margin of error. Anytime after the first 20 datasets, EpipwR terminates once one of these two conditions are met. The default margin of error is 3% (it is recommended that this be between 1 and 5%). For example, if the target margin of error is 3% and the maximum number of data sets is 1000, then the algorithm stops when \(N=1000\) or once the following conditions are met:

\[1.96\cdot\frac{s_{(N)}}{\sqrt{N}}<.03,\qquad N\ge20 \]

where \(s_{(N)}\) is the standard deviation of average power over \(N\) data sets.

Function Inputs

The below table describes the details of the get_power_cont inputs.

Parameter	Description
`dm`	The total number of non-null tests (i.e., the number of CpG sites with methylation having a non-zero correlation with the phenotype).
`Total`	The total number of tests (i.e., total number of CpG sites in the study).
`n`	The sample size of the study - accepts a vector of values.
`fdr_fwer`	Either the false discovery rate (FDR) or the family-wise type I error rate (FWER), depending on `use_fdr.`
`rho_mu`	The average correlation of the non-null tests - accepts a vector of values.
`rho_sd=0`	The standard deviation of the correlation for non-null tests.
`Tissue="Saliva"`	Tissue corresponding to empirical EWAS data (see `?get_power_cont` for valid options).
`Nmax=1000`	The maximum number of data sets that EpipwR will use to determine power. Cannot be less than 20.
`MOE=.03`	The target margin of error of a 95% confidence interval for average power.
`test="Pearson"`	Type of correlation/test utilized. Acceptable inputs are `"Pearson"`, `"Kendall"` and `"Spearman"` .
`use_fdr=T`	Indicates that the false discovery rate should be used as the correction. Otherwise, uses FWER.
`Suppress_Updates=F`	If `T` , removes intermediate status updates.

Methodology (Binary outcome)

For the case-control setting, EpipwR re-implements the pwrEWAS methodology (Graw et al., 2019) with a few key changes to improve speed and interpretation, namely (1) keeping the precision equal between groups rather than variance, (2) requiring users to specify mean and variance for an effect size distribution (rather than a target maximal difference) and (3) generating p-values directly for null tests. Like in get_power_cont, get_power_cc calculates the optimal amount of data sets needed to determine power at a user-specified level of precision.

EpipwR Workflow

Suppose the user is planning an EWAS of 100,000 CpG sites with an estimated 100 of those having a non-zero correlation with continuous phenotype. They would like to test sample sizes of 70, 80, 90, 100 and fixed correlations of 0.3, 0.4 and 0.5. They plan to use an FDR of 5% on Pearson correlation tests, the Saliva empirical EWAS and 1,000 as the maximum number of data sets. Since they hope to have very accurate results, the margin of error is set to be 1%. They would then execute the following code:

output <- get_power_cont(
  dm=100
  ,Total=100000
  ,n=c(70,80,90,100)
  ,fdr_fwer=.05
  ,rho_mu = c(0.3, 0.4, 0.5)
  ,rho_sd = 0
  ,Tissue = "Saliva"
  ,Nmax=1000
  ,MOE=.01
  ,test="pearson"
  ,use_fdr=T
  ,Suppress_updates=F
)

output

which produces the following data frame.

Note that executing this code block will produce very similar but not identical results, since EpipwR is based on random sampling. The first two columns track the setting for which power was calculated (each row is a unique combination of n and rho_mu ). The next four columns show the average power for this setting, the standard deviation of power, the number of data sets used for this setting, and finally the standard error of the average power (column 4 divided by the square root of column 5). Note that when power is very close to 0, often the minimum number of data sets is sufficient to achieve the desired level of precision - the same is true when power is close to 1.

Finally, users can also summarize the output with error bar plots with EpipwR_plot(output). For our previous example, this produces the following plot:

Please note that EpipwR_plot can only run with the data frame output of get_power_cont or get_power_cc - it will not work with data in any other format.

# Load sessioninfo package
library(sessioninfo)

# Display session information
session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R Under development (unstable) (2024-10-21 r87258)
#>  os       Ubuntu 24.04.1 LTS
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language (EN)
#>  collate  C
#>  ctype    en_US.UTF-8
#>  tz       America/New_York
#>  date     2024-10-29
#>  pandoc   3.1.3 @ /usr/bin/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package          * version date (UTC) lib source
#>  AnnotationDbi      1.69.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  AnnotationHub      3.15.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  Biobase            2.67.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  BiocFileCache      2.15.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  BiocGenerics       0.53.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  BiocManager        1.30.25 2024-08-28 [2] CRAN (R 4.5.0)
#>  BiocVersion        3.21.1  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  Biostrings         2.75.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  bit                4.5.0   2024-09-20 [2] CRAN (R 4.5.0)
#>  bit64              4.5.2   2024-09-22 [2] CRAN (R 4.5.0)
#>  blob               1.2.4   2023-03-17 [2] CRAN (R 4.5.0)
#>  bslib              0.8.0   2024-07-29 [2] CRAN (R 4.5.0)
#>  cachem             1.1.0   2024-05-16 [2] CRAN (R 4.5.0)
#>  cli                3.6.3   2024-06-21 [2] CRAN (R 4.5.0)
#>  crayon             1.5.3   2024-06-20 [2] CRAN (R 4.5.0)
#>  curl               5.2.3   2024-09-20 [2] CRAN (R 4.5.0)
#>  DBI                1.2.3   2024-06-02 [2] CRAN (R 4.5.0)
#>  dbplyr             2.5.0   2024-03-19 [2] CRAN (R 4.5.0)
#>  digest             0.6.37  2024-08-19 [2] CRAN (R 4.5.0)
#>  dplyr              1.1.4   2023-11-17 [2] CRAN (R 4.5.0)
#>  EpipwR           * 1.1.0   2024-10-29 [1] Bioconductor 3.21 (R 4.5.0)
#>  EpipwR.data        0.99.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  evaluate           1.0.1   2024-10-10 [2] CRAN (R 4.5.0)
#>  ExperimentHub      2.15.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  fansi              1.0.6   2023-12-08 [2] CRAN (R 4.5.0)
#>  fastmap            1.2.0   2024-05-15 [2] CRAN (R 4.5.0)
#>  filelock           1.0.3   2023-12-11 [2] CRAN (R 4.5.0)
#>  generics           0.1.3   2022-07-05 [2] CRAN (R 4.5.0)
#>  GenomeInfoDb       1.43.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  GenomeInfoDbData   1.2.13  2024-10-23 [2] Bioconductor
#>  glue               1.8.0   2024-09-30 [2] CRAN (R 4.5.0)
#>  htmltools          0.5.8.1 2024-04-04 [2] CRAN (R 4.5.0)
#>  httr               1.4.7   2023-08-15 [2] CRAN (R 4.5.0)
#>  IRanges            2.41.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  jquerylib          0.1.4   2021-04-26 [2] CRAN (R 4.5.0)
#>  jsonlite           1.8.9   2024-09-20 [2] CRAN (R 4.5.0)
#>  KEGGREST           1.47.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  knitr              1.48    2024-07-07 [2] CRAN (R 4.5.0)
#>  lifecycle          1.0.4   2023-11-07 [2] CRAN (R 4.5.0)
#>  magrittr           2.0.3   2022-03-30 [2] CRAN (R 4.5.0)
#>  memoise            2.0.1   2021-11-26 [2] CRAN (R 4.5.0)
#>  pillar             1.9.0   2023-03-22 [2] CRAN (R 4.5.0)
#>  pkgconfig          2.0.3   2019-09-22 [2] CRAN (R 4.5.0)
#>  png                0.1-8   2022-11-29 [2] CRAN (R 4.5.0)
#>  R6                 2.5.1   2021-08-19 [2] CRAN (R 4.5.0)
#>  rappdirs           0.3.3   2021-01-31 [2] CRAN (R 4.5.0)
#>  rlang              1.1.4   2024-06-04 [2] CRAN (R 4.5.0)
#>  rmarkdown          2.28    2024-08-17 [2] CRAN (R 4.5.0)
#>  RSQLite            2.3.7   2024-05-27 [2] CRAN (R 4.5.0)
#>  S4Vectors          0.45.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  sass               0.4.9   2024-03-15 [2] CRAN (R 4.5.0)
#>  sessioninfo      * 1.2.2   2021-12-06 [2] CRAN (R 4.5.0)
#>  tibble             3.2.1   2023-03-20 [2] CRAN (R 4.5.0)
#>  tidyselect         1.2.1   2024-03-11 [2] CRAN (R 4.5.0)
#>  UCSC.utils         1.3.0   2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  utf8               1.2.4   2023-10-22 [2] CRAN (R 4.5.0)
#>  vctrs              0.6.5   2023-12-01 [2] CRAN (R 4.5.0)
#>  xfun               0.48    2024-10-03 [2] CRAN (R 4.5.0)
#>  XVector            0.47.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>  yaml               2.3.10  2024-07-26 [2] CRAN (R 4.5.0)
#>  zlibbioc           1.53.0  2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> 
#>  [1] /tmp/Rtmp4iXNNy/Rinstba1c95f5ee82b
#>  [2] /home/biocbuild/bbs-3.21-bioc/R/site-library
#>  [3] /home/biocbuild/bbs-3.21-bioc/R/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────

References

Barth, J., Reynolds, A.W. (2024). EpipwR: Efficient Power Analysis for EWAS with Continuous Outcomes [preprint]. bioRxiv, 2024.09.06.611713; doi: https://doi.org/10.1101/2024.09.06.611713

Graw, S., Henn, R., Thompson, J. A., and Koestler, D. C. (2019). pwrEWAS: A user-friendly tool for comprehensive power estimation for epigenome wide association studies (EWAS). BMC Bioinformatics, 20(1):218.