Finding the distribution and/or fitting a model (to a biological problem) I am a bioinformatician and work on RNA-Seq data. The data contains a lot of reads (of length 80 bp in my case). These reads are fragments of those genes that were expressed. I map them back to my reference genome which helps answer many questions (gene expression, snp variants, alternative splicing etc...).
I work on alternative splicing, where I am interested in characterizing events where same genes express in more than 1 ways leading to different transcripts. Let's say I my data consists of 4 libraries (or biological replicates) for each of the 2 plant subspecies (closely related ones belonging to same species). And, I would like to look at exon skipping events (where certain exons from coding genes sometimes gets excluded). 
My question is about the statistical analysis towards characterizing these exon skipping events between the two subspecies over those replicates. The basis of my question (or line of thinking and therefore the complexity), however, comes from statistical analysis on RNA-Seq data for gene expression. So, I'll frame the problem from there. Please bear with me.
For every gene, $g_{i}$, its expression is = number of reads mapped to that gene (or count data). For one gene, $g_{1}$, say, the data would look like this.

     Lib1    Lib2    Lib3    Lib4
sA   400     420      600     250
sB   180     229       60     125

Since its count data, given the concentration $q_{1}$ of gene $g_{1}$, biostatisticians usually model the number of reads, $r_{1}$ that map to a $g_{1}$ as a poisson model. It can be shown that $r_{1} \propto q_{i}$, or $r_{1} = sq_{1}$ (poisson parameter). That is, $p(Y = r_{1} | sq_{1}) \sim poisson$. One question here: This means, that if I were to measure for the same gene, the read count and note down the concentration, then, I would end up in a poisson model, right? 
From here, since we can not know the concentration $q_{1}$, they make use of the biological replicates, to estimate the variation. If $Q_{1}$ represents the concentrations over (all) the biological replicates, then under the assumption that this follows a gamma distribution (for the sake of mathematical convenience), the read counts $R_{1}$ for gene $g_{1}$ would then follow a negative binomial distribution. This is an overdispersed poisson model, it seems. Then a statistical test is performed to estimate if the difference in expression is statistically significant under a given $\alpha$. There are R-packages available that implement this.
Now, coming back to my problem, as opposed to looking at every gene, for every exon, $E_{i}$, I find out the number of times this exon is spliced out (= not included in the transcript) and the number of times it is included, in all libraries, in both subspecies.
So, for an exon, $E_{1}$, the data would look like this, for example: The first number in each entry is the number of exon skipping events and the second is the normal events.

     Lib1    Lib2    Lib3    Lib4
sA   2, 80   1, 65   0, 40   2, 66
sB   10, 120 0, 22   8, 90   4, 90

The difference is that I have two values for each entry. I would like to find out, for each exon, if there is a difference between the two subspecies. My immediate thought was to calculate the ratio, for each entry, and then use existing model, however, they are on based on count data. So, my questions are these:
1) In general, is it known/ possible to estimate what distribution the ratio of two counts would follow?
2) Are there any other approaches ( for example, generalized linear models) which could from these data (including dispersion from replicates) help me calculate if the occurrence of exon skipping events between the two subspecies are statistically significant?
PS: In case somethings are not clear enough, I'd love to clarify them (I am by no means a statistician). I would appreciate any thoughts you guys have on this problem. Thanks again!
 A: *

*I think you can just say that the pair of numbers are independent draws from two negative binomial distributions.

*You can use a generalized linear model (with a negative binomial family and log link function) to compare the ratios in the same way that you would use such to compare the two groups.  In place of a model like y ~ species you would use y ~ species * exonskip; the interaction term would correspond to the difference between the log ratios.
Here's your data:
dat <- structure(list(y = c(2, 1, 0, 2, 80, 65, 40, 66, 10, 0, 8, 4, 
           120, 22, 90, 90), lib = structure(c(1L, 2L, 3L, 4L, 1L, 2L, 3L, 
           4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), .Label = c("Lib1", "Lib2", 
           "Lib3", "Lib4"), class = "factor"), species = structure(c(1L, 
           1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("sA", 
           "sB"), class = "factor"), exonskip = structure(c(2L, 2L, 2L, 
           2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L), .Label = c("no", 
           "yes"), class = "factor")), .Names = c("y", "lib", "species", 
           "exonskip"), row.names = c(NA, -16L), class = "data.frame")

That's easy for copy-and-paste but not so clear.  Here are the first few rows:
> dat[1:5,]
   y  lib species exonskip
1  2 Lib1      sA      yes
2  1 Lib2      sA      yes
3  0 Lib3      sA      yes
4  2 Lib4      sA      yes
5 80 Lib1      sA       no

The MASS package includes a function glm.nb for fitting a GLM with the negative binomial family and estimating the over-dispersion parameter.  You could use it as follows:
library(MASS)
out <- glm.nb(y ~ species * exonskip, data=dat)
summary(out)

Here are the key bits of the output:
Coefficients:
                      Estimate Std. Error z value         Pr(>|z|)
(Intercept)             4.1392     0.2334  17.731          < 2e-16
speciessB               0.2491     0.3288   0.758           0.4487
exonskipyes            -3.9160     0.5523  -7.091 0.00000000000133
speciessB:exonskipyes   1.2325     0.6742   1.828           0.0675

              Theta:  4.95 
          Std. Err.:  2.63 

The estimated coefficient for the interaction term, 1.2 (SE = 0.7), is the estimate of the difference (sB - sA) of the log ratio (exon skipped / not skipped).
