I understand how the negative binomial distribution is derived from the poisson and binomial distributions to "model gene counts for differential expression analysis." However, I don't understand what ""model gene counts for differential expression analysis" means, exactly.
Are we trying to fit the negative binomial model for a singular gene across ALL samples? In this case, would the would the X-axis be a certain number of counts, and would the Y-axis be the number of samples in which that number of counts appears?
Or are we trying to fit the negative binomial model for a single gene per sample? In such a case, won't we get a different negative binomial model for every single sample for each gene? How can differential expression analysis be conducted if a different model exists for each gene?
Many thanks for the kind help.
Michael W. Sandel thanks for flagging!
Ak Krishnan, you can think of differential expression analysis for a single gene as a classical univariate regression problem where the response variable is gene expression (Y) and the group information is the covariate of interest (X). Typically, you are interested in estimating the beta coefficient that represents the strength of association between Y and X which is estimated based on the observed data on a per-gene basis. As a result, you get one beta estimate and one p-value per gene, which can be used to access differential expression (e.g., whether P
Gene counts are modelled as a random variable with a negBinom distribution.
This means that the likelihood of an observed gene count is taken from a negBinom distribution. This applies to each and every observed value, across all samples and conditions.
The statistical model considers that the expected value of the distibution depend on some covariables - the experimental factors (e.g. group membership, treatment, tissue, genotype, ...). The values of these covariables are connected via the predictor function with the expected value for the gene count. This predictor function is flexible in a way that it also contains parameters (often simply as coefficients) that determine how strongly the value of a covariable affects the expectation. Each observed gene count may be associated with its individual conmbination of values from these covariables and is thus modelled by a negBinom distribution with its own individual expectation.
The model fit is done by searching those values of the parameters (or coefficients) for which the likelihood of all the observed gene counts is maximal. Approximate standard errors of these parameter estimates are obtained from the curvature of the (log-)likelihood around its maximum, and one can test hypotheses about values of these covariates using approximate likelihood ratio tests (typically Chi² tests).
In a very simple experimental setup you have a couple of gene counts (Y) in a "control" group and another couple of gene counts in a "treated" group. So the only covariable you have is "treatment", which is conveniently coded numerically with X=0 for "control" and X=1 for "treated. The model will thus have a predictor function of the form E(Y|X) = µY(X) = b0 + b1X. In this model, the coefficient b0 is the "intercept" and gives the expected value of Y given that X=0 (expected gene count when measured in the control group), and b1 is the coefficient coding the difference between the expected gene count between X=1 vs. X=0 (the expected effect of the treatment). Practically, it is possible to use values for b0 and b1 that would lead to impossible (negative) expectations. To avoid this, one reformulates the equation as log(E(Y|X)) = b0 + b1X or, correspondingly, E(Y|X) = exp(b0 + b1X). This function "log" which links the predictor function to the expectation is called the link-function. The meaning of the coefficients are then a log expectation and a log ratio.
- Badges
- Science method