# Using the hypergeometric distribution for skipping events in transcriptome sequencing

My question is inspired by this post. However, its a bit more complicated than that to explain. I hope I succeed.

I work with RNA-Seq data on alternative splicing in plants. For this discussion, lets consider exon skipping events. Let $g_{1}, g_{2}, ..., g_{n}$ represent all genes that are expressed. A gene is expressed means that, when I map my RNA-Seq data, which are typically short reads (of length 80 bp in my case), back to the genome, a certain amount of reads would get mapped to that gene (location). The certain amount is based on an arbitrary threshold, if required, lets say, 20 short-reads.

Not used, but still, for sake of clarity, let $m >= n$ be the total number of genes in the genome.

1. The total amount of reads that are mapped to the genome are $R$, say (usually in the order of millions).
2. Let $r_{i}$ be the reads mapped to expressed gene $g_{i}$ (usually hundreds to thousands).
3. Let's say, in gene $g_{i}$, an exon $j$ = $e_{ij}$ has read $r1_{ij}$ times and it gets skipped $r2_{ij}$ times.

Now, I would like to find out the probability that these $r2_{ij}$ exon skipping events occur by chance. My question is, can I use a bivariate hypergeometric distribution which answers the question (modified from Karl's reply from the post linked above):

What is the probability of getting $r2_{ij}$ or more white balls in a sample of size $(r1_{ij}+r2_{ij})$ from an urn with $r_{i}$ white balls and $(R - r_{i})$ black balls?

I am confused because,

1. the reads that map to the gene $g_{i}$ depends on the how much the gene is expressed. Ideally (or biologically), I can not be certain that the gene expression of this gene (=reads mapped) is independent of another (gene) because, it might be that another gene has influenced its expression, but there's no way to know this from the data I have.

2. I have other genes where exon skipping events occur and genes where more than 1 exon skipping events occur. Can I obtain p-values from the hypergeometric distribution and then do a multiple testing correction? or does it become a multivariate problem?

Or this line of thinking doesn't apply to this problem maybe there is an entirely different approach.