In anthony-aylward/chenimbalance: Allelic imbalance mapping procedure from Chen et al. A uniform survey of allele-specific binding and expression over 1000-Genomes-Project individuals (2016)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

First, import the libraries we will use

library(VGAM)
library(chenimbalance)

Set some parameters

FDR_thresh <- 0.05
p <- 0.5
b <- 0.01855469

Prepare the data

cACGT <- accb[c("cA", "cC", "cG", "cT")]
lower <- apply(cACGT, 1, function(x) sort(x, partial = 3)[3])
higher <- apply(cACGT, 1, max)
total <- higher + lower
p_bin <- apply(
  data.frame(2 * mapply(pbinom, lower, total, p)),
  1,
  function(x) min(x, 1)
)
p_betabin <- apply(
  data.frame(
    2 * mapply(pbetabinom, lower, total, p, b)),
    1,
    function(x) min(x, 1)
)
accb[["p_betabin"]] <- p_betabin

Simulations

step <- 0.0001
p_thresh <- c(
  seq(0, 0.01, by = 0.001),
  seq(0.01, 0.1, by = 0.01)[-1],
  seq(0.1, 1, by = 0.1)[-1]
)

Table of empirical counts

w <- as.data.frame(table(total), stringsAsFactors = F)
w <- w[as.numeric(w[,1]) >= 6,]

FP

fp_binomial <- data.frame(
  pval = p_thresh,
  FP_bin = sapply(p_thresh, function(x) fp(w, p, x, "binomial"))
)
head(fp_binomial)

fp_betabinomial <- data.frame(
  pval = p_thresh,
  FP_betabin = sapply(p_thresh, function(x) fp(w, p, x, "betabinomial", b))
)
head(fp_betabinomial)

FDR.txt

tp_bin <- sapply(p_thresh, cutoff, y = p_bin) + 1
tp_betabin <- sapply(p_thresh, cutoff , y = p_betabin) + 1
fdr_bin <- fp_binomial[,2] / tp_bin
fdr_betabin <- fp_betabinomial[,2] / tp_betabin
p_choice_bin <- max(p_thresh[fdr_bin <= FDR_thresh])
p_choice_betabin <- max(p_thresh[fdr_betabin <= FDR_thresh])

fdr_choice_bin <- max(fdr_bin[fdr_bin <= FDR_thresh])
fdr_choice_betabin <- max(fdr_betabin[fdr_betabin <= FDR_thresh])

Bisection method to find p-value

p_choice_bin_1 <- as.data.frame(
  bisect(
    p_bin,
    fp_bin[,2],
    p_choice_bin,
    fdr_choice_bin,
    FDR_thresh,
    step,
    "binomial",
    b = 0,
    w,
    p
  )
)
p_choice_betabin_1 <- as.data.frame(
  bisect(
    p_betabin,
    fp_betabinomial[,2],
    p_choice_betabin,
    fdr_choice_betabin,
    FDR_thresh,
    step,
    "betabinomial",
    b = b,
    w,
    p
  )
)
head(p_choice_betabin_1)

p_choice_bin_1 <- p_choice_bin_1[p_choice_bin_1[,3] > 0,]
p_choice_bin_2 <- p_choice_bin_1[nrow(p_choice_bin_1), 1]
p_choice_betabin_1 <- p_choice_betabin_1[p_choice_betabin_1[,3] > 0,]
p_choice_betabin_2 <- p_choice_betabin_1[nrow(p_choice_betabin_1), 1]

Formatting FDR_txt

FDR_txt <- data.frame(
  pval = p_thresh,
  P_bin = tp_bin,
  FP_bin = fp_binomial[,2],
  FDR_bin = fdr_bin,
  P_betabin = tp_betabin,
  FP_betabin = fp_betabinomial[,2],
  FDR_betabin = fdr_betabin
)
FDR_txt = FDR_txt[-nrow(FDR_txt),]
head(FDR_txt)

Take in counts.txt and filter by p.betabin

interestingHets_betabinom = accb[accb[["p_betabin"]] <= p_choice_betabin,]
head(interestingHets_betabinom)

anthony-aylward/chenimbalance documentation built on May 17, 2019, 12:50 p.m.