Bulk RNA-seq is usually used to study gene expression patterns across different conditions or groups. The technology involves sequencing the entire transcriptome of a biological sample, which allows researchers to identify differentially expressed genes and analyze their functional pathways. By comparing gene expression profiles between different conditions, we can gain insights into the molecular mechanisms underlying these conditions and identify potential therapeutic targets, biomarkers for diseases or drug responses, which can be used for diagnosis, prognosis, or personalized medicine. In this topic, we are discussing the algorithm for identifying differentially expressed genes and tools performing this task.
Differentially expressed genes and experimental design
In the previous topic, we briefly discussed what are differentially expressed genes (DEGs). DEGs are genes that significantly change their expression levels across different conditions. The process of identifying DEGs is called differential gene expression analysis.
But firstly let's discuss experimental design. When conducting any experiment, we should always ask a question about what we want to know. Usually, in bulk rna-seq, we are interested in the mechanism underlying a particular phenotype. To answer this question we should take two groups of samples. The first group is control. Samples from this group represent normal phenotype — healthy or untreated tissue. The second group contain samples of tissue with the phenotype of interest. If you are interested in more than one phenotype, add corresponding groups, but the control group should always persist. Moreover, since identifying DEGs basically represents a statistical test, each group should contain more than one sample, so the test has statistical significance. After choosing sample groups of interest we conduct bulk RNA-seq experiment and obtain data with read counts as discussed in the previous topic.
In DGE analysis we statistically compare gene expression levels between two groups of samples (control vs treated). But the problem here is that we can't use common statistical tests. The parametrical tests, like the t-test or z-test, require data to have a normal distribution, while non-parametric tests are too weak since we usually have few samples in each group. That’s why tools for DGE analysis use generative linear models (GLM) to find the statistical confidence of the difference between the mean expression levels of two groups. Moreover, these models can be applied to data distributed not normally. It may sound complicated, so let’s look at this in more detail.
Linear models
Before we proceed, let’s revise some material about linear regression models. Linear regression is a statistical approach that models the relationship between a dependent variable (y) and one or more independent variables (x). It is used to predict the value of the dependent variable based on the values of the independent variables. The model assumes a linear relationship between the variables, meaning that the change in the dependent variable is proportional to the change in the independent variable(s). Therefore, the relationship between x and y looks as follows: y=ax+b. The goal of linear regression is to find the best-fit line (a and b coefficients) that represents the relationship between the variables, which can then be used for prediction or inference. The coefficients found can be then tested for statistical significance.
Now let’s explore linear regression application on a simplified example.
Imagine we have data about gene expression levels at different ages. We can use a linear regression model to find the relationship between the gene expression level and the age (age=a*expression+b). Knowing the coefficients a and b it becomes possible to estimate the change between expression levels in two ages. However, in real life, we have samples from discrete conditions like young and old, healthy and diseased. But we still can apply linear regression. In this case, we would have samples from two groups young and old, and the expression level of a particular gene would look as follows expression=b1*young+b2*old. The previous equation can also be present like this expression=b0+b1*old. картинки In this case, b1 represents the difference between gene expression levels in young and old samples. And that’s it.
In DEG analysis we wanted to find genes that significantly change their expression and calculate the magnitude of the change. Therefore, we take bulk RNA-seq data about particular gene expression levels in different samples, fit a linear regression model (in other words, find coefficients b0 and b1), calculate the statistical significance of the coefficients and determine the magnitude of the change. Then repeat this procedure for all genes present in the data.
However, the process is not so simple. Firstly, in the real world gene expression levels depend not only on age, they also can be affected by a lot of other variables, such as sex, health state, environmental conditions etc. Ideally, all these variables should be included in the model.
Secondly, linear regression can be applied if residues are normally distributed and have constant variance (a.k.a. homoscedasticity). However, in the case of bulk RNA-seq data, this criterion is not met.
Let’s have a closer look at this problem. We are considering linear regression for categorical variables, which means our data is obtained from one of those groups. And each gene in each group has its distribution of expression. When talking about expression we imply read counts of a particular gene. In the first approximation, the more read counts the sample has the higher the gene expression level is. Thus, distributions of read counts between samples would differ by mean and variance. But we still haven’t said what is this distribution. It was shown that the mean-variance relationship of read counts in bulk RNA-seq data resembles negative binomial distribution. Thus, bulk RNA-seq data does not follow Normal distribution and therefore does not meet the conditions of linear regression.
And now let’s look into how we could use generalized linear models. In simple terms, GLM first transforms data so it follows a normal distribution and then fits a linear regression to the transformed data just as usual linear regression. The complicated step is data transformation. This procedure varies between different tools.