I have RNA-seq data that I'd like to test for allele-specific expression, and I have three biological replicates. For an individual locus, I have thousands of sequencing reads, and my null hypothesis is that the observed count of allele A, X, is given by X~Bin(n,p) where n is my total read count and p=0.5 (there are only two alleles). For a single (example) replicate, I have 8000 reads and they're all one allele, I understand how to calculate a p-value for the observed counts. E.g., in R
binom.test(x=8000,n=8000,p=0.5)
But how do I test for bias across all three samples? Say I have sample 1 (0 out of 8000 reads), sample 2 (0 out of 5000 reads), and sample 3 (0 out of 6000 reads). I would like to calculate a p-value for this result given my null that p=0.5 for every sample.
I've simulated this (3 samples of a few thousand reads with p=0.5), using
replicate(10000,rbinom(3,5000,0.5))
Of course I never observe a single 0 count let alone three (I never actually see anything less than 2000). But I don't know what the appropriate way to formally quantify this is or what the correct test would be.
