Learn R Programming

powerEQTL (version 0.3.4)

powerEQTL.SLR: Power Calculation for EQTL Analysis Based on Simple Linear Regression

Description

Power calculation for eQTL analysis that tests if a SNP is associated to a gene probe by using simple linear regression. This function can be used to calculate one of the 4 parameters (power, sample size, minimum detectable slope, and minimum allowable MAF) by setting the corresponding parameter as NULL and providing values for the other 3 parameters.

Usage

powerEQTL.SLR(
  MAF,
  slope = 0.13,
  n = 200,
  power = NULL,
  sigma.y = 0.13,
  FWER = 0.05,
  nTests = 2e+05,
  n.lower = 2.01,
  n.upper = 1e+30)

Arguments

MAF

numeric. Minor allele frequency.

slope

numeric. Slope of the simple linear regression.

n

integer. Total number of subjects.

power

numeric. Power for testing if the slope is equal to zero.

sigma.y

numeric. Standard deviation of the outcome \(y_i\) in simple linear regression.

FWER

numeric. Family-wise Type I error rate.

nTests

integer. Number of tests (i.e., number of all (SNP, gene) pairs) in eQTL analysis.

n.lower

numeric. Lower bound of the total number of subjects. Only used when calculating total number of subjects.

n.upper

numeric. Upper bound of the total number of subjects. Only used when calculating total number of subjects.

Value

power if the input parameter power = NULL.

sample size (total number of subjects) if the input parameter n = NULL;

minimum detectable slope if the input parameter slope = NULL;

minimum allowable MAF if the input parameter MAF = NULL.

Details

To test if a SNP is associated with a gene probe, we use the simple linear regression $$y_i = \beta_0+\beta_1 x_i + \epsilon_i,$$ where \(y_i\) is the gene expression level of the \(i\)-th subject, \(x_i\) is the genotype of the \(i\)-th subject, and \(\epsilon_i\) is the random error term with mean zero and standard deviation \(\sigma\). Additive coding for genotype is used. That is, \(x_i=0\) indicates wildtype homozygotes; \(x_i=1\) indicates heterozygotes; and \(x_i=2\) indicates mutation heterozygotes.

To test if the SNP is associated with the gene probe, we test the null hypothesis \(H_0: \beta_1=0\) versus the alternative hypothesis \(H_1: \beta_1 = \delta\), where \(\delta\neq 0\).

Denote \(\theta\) as the minor allele frequency (MAF) of the SNP. Under Hardy-Weinberg equilibrium, we can calculate the variance of genotype of the SNP: \(\sigma^2_x=2 \theta (1-\theta)\), where \(\sigma^2_x\) is the variance of the predictor (i.e. the SNP) \(x_i\).

The exact power calculation formula can be derived as $$1-T_{n-2, \lambda}(t_{n-2}(\alpha/2)) + T_{n-2, \lambda}(-t_{n-2}(\alpha/2)),$$ where \(T_{n-2, \lambda}(a)\) is the value at \(a\) of cumulative distribution function of non-central t distribution with \(n-2\) degrees of freedom and non-centrality parameter \(\lambda=\delta/\sqrt{\sigma^2/[(n-1)\tilde{\sigma}^2_{x}]}\). And \(\tilde{\sigma}^2_{x}=\sum_{i=1}^n(x_i - \bar{x})^2/(n-1)\).

Dupont and Plummer (1998) mentioned the following relationship: $$\sigma^2 = \sigma^2_y - \beta_1^2 \sigma^2_x.$$ So we can plug in the above equation to the power calculation formula.

Under Hardy-Weinberg equilibrium, we have \(\sigma_x^2=2\theta(1-\theta)\), where \(\theta\) is the minor allele frequency (MAF).

Hence, the non-centrality parameter can be rewritten as $$\lambda=\frac{\delta}{\sqrt{ \left(\sigma_y^2 - \delta^2 2\left(1-\hat{\theta}\right)\hat{\theta}\right)/ \left[(n-1)2\left(1-\hat{\theta}\right)\hat{\theta}\right] }}$$

We adopted the parameters from the GTEx cohort (see the Power analysis" section of Nature Genetics, 2013; https://www.nature.com/articles/ng.2653), where they modeled the expression data as having a log-normal distribution with a log standard deviation of 0.13 within each genotype class (AA, AB, BB). This level of noise is based on estimates from initial GTEx data. In their power analysis, they assumed the across-genotype difference delta = 0.13 (i.e., equivalent to detecting a log expression change similar to the standard deviation within a single genotype class).

References

Dupont, W.D. and Plummer, W.D.. Power and Sample Size Calculations for Studies Involving Linear Regression. Controlled Clinical Trials. 1998;19:589-601.

Dong X, Li X, Chang T-W, Scherzer CR, Weiss ST, and Qiu W. powerEQTL: An R package and shiny application for sample size and power calculation of bulk tissue and single-cell eQTL analysis. Bioinformatics, 2021;, btab385

Examples

Run this code
# NOT RUN {
# calculate power
powerEQTL.SLR(
  MAF = 0.1,
  slope = 0.13,
  n = 179,
  power = NULL,
  sigma.y = 0.13,
  FWER = 0.05,
  nTests = 2e+05)
  
# calculate sample size (total number of subjects)
powerEQTL.SLR(
  MAF = 0.1,
  slope = 0.13,
  n = NULL,
  power = 0.8,
  sigma.y = 0.13,
  FWER = 0.05,
  nTests = 2e+05)
  
# calculate minimum detectable slope
powerEQTL.SLR(
  MAF = 0.1,
  slope = NULL,
  n = 179,
  power = 0.8,
  sigma.y = 0.13,
  FWER = 0.05,
  nTests = 2e+05)  
  
# calculate minimum allowable MAF
powerEQTL.SLR(
  MAF = NULL,
  slope = 0.13,
  n = 179,
  power = 0.8,
  sigma.y = 0.13,
  FWER = 0.05,
  nTests = 2e+05) 

# }

Run the code above in your browser using DataLab