pmf_impr.Rdpmf_impr calculates the probability of observing population-level RNAseq data (reference- and variant allele counts per sample)
assuming an imprinted binomial mixture model with the mean reference allele fraction in heterozygotes fixed at 0.5 (no allelic bias).
More specifically, the formula used is (though under-the-hood not using R's built-in dbinom function but our own mathematical implementation):
pr * dbinom(x = ref_counts, size = ref_counts + var_counts, prob = 1-SE) +
0.5 * prv * dbinom(x = ref_counts, size = ref_counts + var_counts, prob = (0.5-(impr/2)) / (1-(impr/2)) * (1-SE) + (0.5/1-(impr/2)) * SE) +
0.5 * prv * dbinom(x = ref_counts, size = ref_counts + var_counts, prob = (0.5-(impr/2)) / (1-(impr/2)) * SE + (0.5/1-(impr/2)) * (1-SE)) +
pv * dbinom(x = ref_counts, size = ref_counts + var_counts, prob = SE)
pr, pv and prv are calculated assuming Hardy-Weinberg-Equilibrium with given input parameters (inbreeding coefficient inbr and reference allele frequency allelefreq)
pmf_impr(ref_counts, var_counts, allelefreq, impr, SE, inbr = 0)Number or Numeric vector Reference count(s).
Number or Numeric vector. Variant count(s).
Number. Allele frequency.
Number. Degree of imprinting.
Number. Sequencing error rate.
Number. Degree of inbreeding (default = 0).
Probability of observing ref_counts and var_counts
pmf_impr(10, 10, 0.5, 0, 0.002, 0.12)
#> [1] 0.0775267
pmf_impr(10, 10, 0.5, 1, 0.002)
#> [1] 1.854402e-22
pmf_impr(0, 10, 0.5, 1, 0.002, 0.12)
#> [1] 0.4900895
pmf_impr(c(5, 8, 10, 3, 5, 6, 23), c(8, 8, 6, 4, 4, 10, 0), 0.5, 1, 0.002, 0.12)
#> [1] 2.026483e-11 3.242372e-18 2.511768e-13 1.391617e-07 9.999602e-10
#> [6] 2.511768e-13 4.774990e-01
pmf_impr(c(5, 0, 0, 3, 5, 1, 23), c(1, 8, 6, 2, 0, 10, 0), 0.5, 1, 0.002)
#> [1] 5.940240e-03 4.920558e-01 4.940299e-01 1.992008e-05 4.950200e-01
#> [6] 1.078197e-02 4.774990e-01