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
- Is this method of estimating, assuming negative binomial glm model right? Particularly, the use of
- 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.