**2.3k**wrote:

I have RNA-seq data from two treatment groups (three replicates each) and I am trying to replicate the likelihood ratio test outlined on page 1516 of this paper by Marioni et al: http://www.ncbi.nlm.nih.gov/pubmed/18550803. So far I have been able to fit a Poisson regression model for each gene using the following function and command:

```
poissonGLM <- function(gene) {
data <- data.frame(treatment=pdata$treatment, count=gene)
model <- glm(count ~ treatment, data=data, family=poisson)
mean <- predict(model.pois, type="response"),
stdev <- predict(model.pois, type="response", se.fit=T)$se.fit
}
results <- apply(edata, 1, poissonGLM)
```

For each gene, this computes the expected mean and expected variance for each treatment group. However, I'm unsure how to use this model to compute the maximum likelihood estimates both under the null and alternative hypothesis (as described in the paper), and generate a P-value from the Chi-squared distribution.

**470**• written 22 months ago by James Ashmore •

**2.3k**