DMRInfer: P-value calculation given Wald statistics.

In haowulab/TRESS: Toolbox for mRNA epigenetics sequencing analysis

P-value calculation given Wald statistics.

Description

This function calculates p-values for candidate DMRs given Wald statistics, using respectively two-component mixtures of normal, truncated normal and standard normal distributions.

Usage

DMRInfer(stat, nullModel = "standN")

Arguments

stat	A vector containing Wald statistics of all candidate DMRs.
nullModel	A character to specify a method to calculate p-value based on the statistics. It can be "standN", "2mix" and "trunN" for standard normal, two-component mixed gaussian and truncated normal respectively. Defult is "standN".

Details

In addition to standard normal distribution, TRESS provides another two distributions to calculate p-values given Wald statistics, in case that statistics are inflated by potentially underestimated dispersion in data.

One is two-component mixtures of normal distribution, where TRESS assumes that,

T_w \sim p N(0, σ_0^2) + (1-p)N(0, σ_1^2)

where T_w is Wald statistics, σ_0 and σ_1 are standard deviations of distribution that T_w follows under null and alternative hypothesis in TRESS_DMRtest. P-values for T_w are calculated using N(0, \hat{σ}_0^2).

The other one is truncated normal distribution, where TRESS assumes that,

tT_w\sim p N(0, σ^2)

where tT_w = T_w \in [-b, b]. Here, T_w within range [-b, b] is assumed to sampled from null distribuion N(0, σ^2). For this truncated normal distribution, TRESS explores different values for boundary b ranging from 1.5 to 2 by step 0.1. TRESS estimates a \hat{σ} for each of 6 boundaries. If \hat{σ}_{max} - \hat{σ}_{min} > 0.5, TRESS calculates p-values for T_w using N(0, \hat{σ}_{min}^2). Otherwise, p-values are obtained using N(0, \hat{σ}^2), with \hat{σ} estimated under b=2.

Value

This function returns a dataframe containing p-values and Benjamini-Hochberg procedure adjusted p-values.

Examples

### use a randomly generated toy data as an
### illustrate of DMRInfer
set.seed(12345)
p = 0.8
nsites = 10000
flag.TP = rep(NA, nsites)
T_w = rep(NA, nsites)
for (i in seq_len(nsites)) {
  u = runif(1, min = 0, max = 1)
  if(u < p){
    flag.TP[i] = FALSE
    T_w[i] = rnorm(1, 0, sd = 1)
  }else{
    flag.TP[i] = TRUE
    T_w[i] = rnorm(1, 0, sd = 5)
  }
}
res = DMRInfer(stat = T_w, nullModel = "standN")
sum(res$padj < 0.05 & !flag.TP)/sum(res$padj < 0.05)

res = DMRInfer(stat = T_w, nullModel = "2mix")
sum(res$padj < 0.05 & !flag.TP)/sum(res$padj < 0.05)

res = DMRInfer(stat = T_w, nullModel = "trunN")
sum(res$padj < 0.05 & !flag.TP)/sum(res$padj < 0.05)

haowulab/TRESS documentation built on Aug. 27, 2022, 7:11 p.m.