# Use of a negative binomial model for fitting alternative splicing event

I am working on RNA-Seq data (on alternative splicing). Let's say I am looking at a particular type of alternative splicing event - exon skipping. For each intron (or junction), I look if it is normally spliced or if there is an exon skipping event happening. I have three biological replicates ( for each junction, I have 3 set of values). To sum it up, my data set looks like this:

Junction 1:
Rep1    Rep2    Rep3
exonSkip    8      0       0
normal      12     6       8

Junction 2:
Rep1    Rep2    Rep3
exonSkip    5      9       8
normal      58     60      44
....
....


My objective is, for each junction, from these replicates, to find out if that particular exon-skipping event is statistically significant. Initially, I summed up all the values for exonSkip and normal separately (in the first case, 8 and 26) and then concluded there are at least 2 exonSkip events. However, I came to know its not the best and that there are better ways.

1) From the literature (on gene expression), I came to know that the biological replicates should be used to obtain an estimate. Since these are reads = count data, and they happen to have a high degree of dispersion among replicates, a negative binomial distribution is suggested.

So, I used a glm.nb model from the MASS library as follows:

# R-code
# Junction 1
require(MASS)
dat1 <- data.frame(y=as.numeric(c(8,0,0,12,6,8)),
exonSkip=as.factor(c("yes","yes","yes","no","no","no")))
out1 <- glm.nb( y ~ exonSkip, data=dat1)
summary(out1)
# Junction 2
dat2 <- data.frame(y=as.numeric(c(5,9,8,58,60,44)),
exonSkip=as.factor(c("yes","yes","yes","no","no","no")))
out2 <- glm.nb( y ~ exonSkip, data=dat2)
summary(out2)


For Junction 1: I got p=0.143 for exonSkipyes.

Question: Does this mean that the nb model fit can not be trusted?

For Junction 2: I got p<2e-16. However, the test glm.nb gave a

*warning: In theta.ml ... Iteration limit reached*.


Question: Is this okay?

Now, from the model, I then used the predict function to estimate the values for exonSkip and normal from these 3 replicates with the model fit.

exp(predict(out1, dat1))
# Result: 2.6667 and 8.6667 = 3 and 9
exp(predict(out2, dat2))
# Result: 7.333 and 54.000 = 7 and 54


Questions:

1. Is this method of estimating, assuming negative binomial glm model right? Particularly, the use of predict function.
2. Can I still use predict (as I have used above) in case of a non-significant fit? If not, then how else can I get an estimate?

2) Even if manage to get the estimate, I have again 2 numbers: 1 for exonSkip and other for Normal. From here, I would like to obtain a measure (or p-value) of how significant it is. How I can go about this?

I think the p-value from the glm.nb is how significant the model fits the data...
One way I thought of is: If there were totally X+Y (X = total exonSkip and Y = total Normal) events, then from a sample of b events, if I get a exonSkip, then it would follow a hypergeometric distribution and I could obtain the p-value as,

sum(dhyper(a:b, X, Y, b))


Is this correct?

I'd appreciate any feedback.

-
Before you can do significance testing, you need to specify the null and alternative hypotheses. If there is no alternative splicing, what do you expect to see? What parameter(s) should be equal to what value(s)? How this should change if there is alternative splicing? –  Aniko Nov 9 '11 at 14:46
@Aniko, I am looking at those places where there is exon-skipping event in either 1 or both time points. Therefore, I am interested in finding out, if there are exon-skipping events that are present in 1 time point but not in the other. I guess my null hypothesis (for each junction I test) would be: Given this occurrence of normal events for each time point, there is no difference in the occurrence of exon-skipping events between the two time points. –  Arun Nov 9 '11 at 15:11
A question about MASS package: In their documentation, after using glm.nb to fit the model, they also encourage using anova on the fitted model to obtain chisquare-test based p-values. What is the difference between the p-values obtained from glm.nb and this?? Anyone? –  Arun Nov 9 '11 at 15:14
I don't see anything relating to different time-points in your original question. Is it an additional data set like what you have above? Could you edit your question to clarify? –  Aniko Nov 9 '11 at 18:20
I'm sorry. I wrote it in a hurry while leaving my office! For the 2nd question, about p-values... Given that there are X exon skipping events in total and Y normal events in total, how probable is it to obtain "a" exon skipping events out of a sample of size "b"? This is my question. –  Arun Nov 9 '11 at 19:02

1. You need a testable hypothesis before getting p-values. That means that you need to describe what the data would look like without and with the effect of interest. You can do parameter estimation without hypothesis testing, but then you need to specify what parameter is of interest, and there are no p-values.
2. While your values are counts, they represent the number of reads falling into each category with a fixed total per experimental run. Or more precisely, the total is random, but you want to condition on it (you keep saying things like "out of a sample of b events"). If so, your counts do not have a negative binomial distribution, but rather a binomial distribution.

The following assumes such a conditioning on the total number of reads per experimental run. The right analysis approach still depends on whether your goal is to estimate hte proportion of exon-skipping reads for each junction (with a confidence interval) or do you want to be able to estimate the entire distribution of the number of exon-skipping reads given the total number of reads. You seem to be asking for the latter, but I am not sure whether that's what you really need, as it is not clear what are you planning to do with the results of the analysis.

You have very little information to inform the model, so you have to make major assumptions. Based on a preliminary analysis using quasi-binomial regression it appears that there might be overdispersion: the probability of an exon-skipping read varies somewhat between the replicates. If that's true, binomial regression cannot be used, but you could consider beta-binomial regression. Here I will show the binomial regression.

First, you have to set up the data so that counts from the same replicate are in the same row.

d1 <- data.frame(Rep=1:3, Skip=c(8,0,0), Normal=c(12,6,8))


And then we can use a binomial glm:

m1 <- glm(cbind(Skip,Normal) ~ 1, data=d1, family=binomial)
summary(m1)

Call:
glm(formula = cbind(Skip, Normal) ~ 1, family = binomial, data = d1)

Deviance Residuals:
1       2       3
1.634  -1.794  -2.072

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.1787     0.4043  -2.915  0.00355 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 10.18  on 2  degrees of freedom
Residual deviance: 10.18  on 2  degrees of freedom
AIC: 15.613

Number of Fisher Scoring iterations: 4


Note that the p-value shown is probably meaningless to you. It tests whether the parameter is 0, and that corresponds to testing whether the probability of exon-skipping is 0.5.

The predict function can be used to get the predicted probability of exon-skipping for each replicate, and the dbinom function can can be used to get the individual response probabilities:

p1 <- predict(m1, type="response")
> p1
1         2         3
0.2352941 0.2352941 0.2352941
> newtotal <- 10
> dbinom(0:newtotal, size=newtotal, p=p1[1])
[1] 6.838240e-02 2.104074e-01 2.913333e-01 2.390427e-01 1.287153e-01 4.752565e-02
[7] 1.218606e-02 2.142605e-03 2.472236e-04 1.690418e-05 5.201286e-07


So, for example, the probability of having 0 exon-skipping reads out of 10 reads at junciton 1 is estimated to be 0.06828.

-
Aniko, Thanks for your reply. I have a couple of things to say and ask. So, this is going to take more than 1 comment. 1) Yes, you're right. There is no p-value here. From my RNA-Seq data, I have some million reads. Suppose you have a gene region E1,I1,E2,I2,E3 (I=intron, E=exon, this is not the whole gene, just a portion). Every E and I have a start and end coordinate. Some spliced reads map to I1 (meaning part maps to E1 and other part to E2), and some spliced reads map to I2 (meaning a part maps to E2 and other part to E3). But sometimes a read maps across I1 and I2 (reads across E1 and E3) –  Arun Nov 10 '11 at 7:39