# 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!

• Why not just add across libraries and do a $\chi^2$ test, of 5:251 vs 22:322? – Karl Oct 19 '11 at 2:44
• Yes, I could. And that's how I planned. However, analogous to my problem, from literature on gene expression, people used fisher test, with which statisticians aren't usually happy with. I may be confusing things here. But the general argument is that fisher test does not account for dispersion across biological replicates. That's the reason for model fitting and variance/dispersion estimation. I just assumed this is more or less relevant to my problem as well. – Arun Oct 19 '11 at 8:05
• I see. Well, if one were generally comparing the 251 and 322 with some GLM with a negative binomial model, you could use the same sort of model to compare the ratios. – Karl Oct 19 '11 at 11:09
• Hi Karl, Thank you again for your reply. Here, 251 and 322 are the total counts, which are similar to the models for gene expression. However, to fit the model using GLM, you know in this case that a negative binomial (or quasi-poisson) is deemed (or the right choice). When I take the ratios, do you think it is okay can estimate the model parameters without knowing what model (or distribution) this ratio would follow? I am not so informed about this part of statistics... If not, at the moment, as you say, chi-square seems to be the best option. – Arun Oct 19 '11 at 11:49
• I would take it that the 5 and 251 are draws from two independent negative binomials. I'll try to write a more complete answer later today, though perhaps someone else will pick up on this and answer before me. – Karl Oct 19 '11 at 13:12

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

2. 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.

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).

• that's great Karl. Thank you so much. Now I am spending sometime in understanding this summarized data. Also, I ran a glm model with poisson distribution using lmtest, and used lrtest to check between the two models, and this proves negative binomial is more appropriate in modeling dispersion. If I understand it right, the last row, last entry is the p-value I have to look for. I'll verify this anyways and should I have more questions, I hope you'll be able to answer them when you find time. Thank you very much! – Arun Oct 19 '11 at 19:26
• Okay, its harder to interpret this data than I had hoped. Suppose, y=(0,0,0,0,80,65,40,66,10,6,8,6,120,22,90,90), then, I get a value for interaction term as 22 (SE=7770) and Pr(>|Z|) = 0.998. I understand then, that the difference of log ratios is huge. But I see a very high p-value. I don't understand this part. Basically, from the estimation of variance, shouldn't I get a a p-value based on a test statistic? If the p-value was indeed 0.998, then, its strange to me as I would find this event(=y) very interesting. Am I anywhere close? – Arun Oct 19 '11 at 23:44
• @Arun: The four 0's are messing things up. If you change one of them to 1, things look great; with the 4 0's, you get estimated log ratios that are $\pm \infty$. – Karl Oct 19 '11 at 23:53
• once again thanks. It seems to solve the problem. I just read about this => zero-inflation(?) Adding 1 should be okay. I have couple more questions, please bear with me!! 1) I read on some other forums anova(out.nb) gives the likelihood ratio (and p-values, again). What's this for and what p-values are those? 2) how is the link function related to the independent variables here? (or how does y ~ speciesexonskip answer the difference between ratio of spA and spB?). I tried to grasp this... couldn't. 3) R question: How to extract *estimate, SE and p-value from summary? Thanks! – Arun Oct 20 '11 at 0:23