Inference of differentially expressed genes using generalized linear mixed models in a pairwise fashion

Preprocessing of the count matrices

This step is conducted before DEGs identification and it aims to reduce the technical variation in each experiment. It entails scale correction, matrix normalization, and filtering poorly expressed or non-expressed genes. The DEGs inference method proposed here is based upon differential expression analysis, that is, it models changes in gene expression levels across experimental conditions rather than simply measuring expression levels. The flowchart representation for this preprocessing step is displayed in Fig. 1.

We also used generalized linear models (GLMs) to validate the preprocessing step by comparing the DEGs detected to those produced by DESeq2 and edgeR. For that, the experimental design’s fixed effects were compared pairwisely, with the respective number of biological replicates for each experimental condition. The overdispersed genes, commonly modeled by GLM by means of negative binomial distribution with only fixed effects, were fitted using the MASS package (Venables & Ripley, 2002), version 7.3-53.1. The regression coefficients of the GLM derived for each gene were analyzed using the Wald test from the aod package (Lesnoff & Lancelot, 2012), version 1.3.1.

Scale correction and normalization of the count matrices

Scale correction and normalization are preprocessing steps that reduce experimental technical bias and are commonly used in transcriptomic studies (Li et al., 2017, 2014). Thus, the matrices must contain the raw read counts rather than any prior normalizing ones. If the user adds the converted counts using any normalizing approach, the DEG inference will not accurately estimate the mean-variance relationship in the data, significantly impacting the statistical analysis. Denoting [A] as the count matrix with rows representing genes and columns representing samples, scale correction for each sequenced library begins by summing its column’s elements into a vector $\overrightarrow{\nu }=({\nu }_1,{\nu }_2,\dots ,{\nu }_n)$, where ${\nu }_j={\sum }_{i=1}^n{a}_{ij},j\in \{1,...,n\}$ and $a_{ij}$ denotes entries in the matrix [A]. Next, its minimum element $\lambda$ is obtained, that is, $\lambda \stackrel{\mathrm{\Delta }}{=}\underset{j}{min}({\nu }_j).$

Next, we define the matrix $[A^{\ast }]$ with the same dimensions as [A], representing the corrected matrix. Its element in row $i$ and column $j$ is expressed by:

$$[A^{\ast }{]}_{ij}=a_{ij}^{\ast }=\frac{\lambda }{v_j}\cdot a_{ij},\mathrm{\forall }i\in \{1,2,...,m\}\mathrm{a}\mathrm{n}\mathrm{d}j\in \{1,2,...,n\}.$$

The relative log expression approach (Anders & Huber, 2010) is used to implement the normalizing step of matrix [A*]. Since this technique assumes that genes in the counting matrices are not differentially expressed, normalization aids in identifying those differentially expressed in a pairwise fashion.

Filtering of unexpressed and low-expressed genes

Genes with zero reads in all biological replicates are eliminated. To identify non-expressed genes, genes that also have null counts are compared across the replicates using Wilcoxon’s nonparametric hypothesis test with a 5% significance level. For this purpose, the stats package, version 3.6.2, from the R programming language (R Core Team, 2022), version 4.0.4 is used. Genes not significantly different from null counts are eliminated from the analysis, and all the remaining are grouped with genes expressed in all replicates.

Genes with the same number of reads across all biological replicates in all conditions are also deleted since they are not differentially expressed. The dataset is then filtered for genes with low counts by means of the gene’s average count per million (CPM). Those genes with values smaller or equal to one are removed because weakly expressed genes reduce the sensitivity techniques for DEGs detection (Sha, Phan & Wang, 2015). The final datasets from the preprocessing step are used for dispersion evaluation and to infer the biological coefficient of variation (BCV).

To classify between equi- or overdispersed, the read counts of each gene along the biological replicates are modeled using the generalized linear model with the Poisson distribution. This model is estimated by means of the parglm package (Christoffersen, 2021), version 0.1.7, used in parallel mode, included in R programming language. Subsequently, Pearson’s chi-square test ratio is used to validate gene dispersion. The gene is classified as equidispersed if its value is less than or equal to one and overdispersed otherwise (Payne et al., 2018; Roback & Legler, 2021). The BCV inference is estimated for each gene before and after the preprocessing step using the estimateCommonDisp function from the edgeR package (Robinson, McCarthy & Smyth, 2010), version 3.32.1.

DEGRE computational tool development

Primarily, DEGRE executes the preprocessing step previously described. Following this, a generalized linear mixed model (GLMM) with a negative binomial distribution is used to estimate the regression coefficients for inferring differentially expressed genes. For the GLMM fitting step in DEGRE, which employs both fixed and random effects, we used the glmmTMB package (Brooks et al., 2017), version 1.1.2.3. An example of the GLMM model (Stroup, 2012) is given below

$$g(\mu |b)=X\beta +Zb,$$where (i) $g(\cdot )$ is a link function, (ii) $\mu |b$ is the conditional expectation of the observations given $b$, (iii) X is the design matrix for fixed effects, (iv) $\beta$ is the regression coefficients vector for fixed effects, (v) Z is the design matrix for random effects, and (vi) $b$ is the random effects’ vector.

The regression coefficients of the models derived for each gene are analyzed using the Wald test (Wald, 1943), also from the glmmTMB package. The P-value attributed to each gene must be adjusted by Bonferroni (Dunn, 1961) or Benjamini-Hochberg (Benjamini & Hochberg, 1995) techniques. The final dataset identifying DEGs and non-DEGs from the previous analysis is used to calculate (i) log2 fold-change values, by means of which a gene is classified as up- or downregulated and (ii) average logCPM. DEGRE does not define a specific value for the log2 fold-change cutoff and the user can select the desired cutoff value.

Simulations of gene read count matrices in RNA-Seq

We created simulated count matrices of gene reads in which DEGs were known. The matrix rows represent the genes, while the matrix columns represent the biological replicates for two experimental conditions. These matrices were simulated with raw read counts for 38,293 genes in 2, 3, 4, 5, 6, 7, 8, 9, 10, and 15 biological replicates (using the “numberReps” parameter) in the mosim function from the MOSim package (Martínez-Mira, Conesa & Tarazona, 2018), version 1.4.0. The “times” parameter was set to zero corresponding to a non-longitudinal dataset and the “omics” parameter was defined as RNA-Seq. These numbers of replicates were chosen since they are commonly employed in differential expression analysis (Schurch et al., 2016). We investigated DEGs proportions ranging by 10%, 20%, and 30% for each matrix by setting the “diffGenes” parameter also in the mosim function. Moreover, for each matrix, the ratio of upregulated to downregulated genes is 1:1.

Estimating the random effects for the simulated count matrices

We modeled the random effects using a variant of the normal probability distribution whose mean was estimated using public data in the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/gds). The transcriptome studies chosen may have the same number of biological replicates as the previously described simulated datasets to prevent any bias in the mean parameter. The increased number of biological replicates might introduce additional biological variability into the datasets. Table 1 lists the transcriptome data identifiers (GEO IDs) utilized for random effects estimation.

To simulate random effects, the raw data of each GEO is submitted to the preprocessing step (Additional File 1: Fig. S1). In order to determine read differences between the experimental conditions, the gene count average along the biological replicates is calculated in each transcriptome data for the two experimental conditions. The average differences along the experimental conditions are calculated and its average is used to estimate the mean parameter under the assumption of normal distribution. Seven values for the random effect standard deviations were selected, namely 100, 300, 600, 900, 1,200, 2,000, and 3,000, to cover a wide range of variability situations. The normal distributions simulated were then added to matrices with only fixed effects. At this step, if the random effect reduces gene reads, decreasing it to less than zero, the count value is replaced by zero. All gene reads were rounded to the nearest integer. The raw matrices generated with fixed and random effects were then the input for DESeq2, edgeR, and DEGRE.

Evaluation of the identification of DEGs by the computational tools

The simulation of the count matrices is the basis for determining which genes are potential DEGs. This inference may be evaluated using recall, precision, and accuracy metrics.

Public RNA-Seq data from patients with bipolar disorder for evaluation of the DEGRE computational tool

For the biological validation of DEGRE, we selected the dataset coded GSE80336 (GEO ID) that comprises RNA expression of brain samples of healthy people (n = 18) and patients with BD (n = 18). The corresponding read count matrix was submitted to the BCV estimation before and after the preprocessing step. The results found by DEGRE were compared to those generated by DESeq2 in a previous study (Pacifico & Davis, 2017). We used 5% of the FDR and Benjamini-Hochberg P-value correction for this comparison. Here, we used disease vs. control as the fixed effect and gender as the random effect (Thoral et al., 2021; Donovan et al., 2020; Kleyheeg et al., 2018) since studies indicate that gender can impact gene expression in transcriptomic analyses (Oliva et al., 2020), even if they are samples from the brain (Raznahan & Disteche, 2021). For Gene Ontology (GO) analysis, we used Ensembl IDs as input for the function basicProfile set to level 4, from the goProfiles R package (Sanchez, Ocana & Salicru, 2022), version 1.58.0.

Filtering of the genes in the preprocessing step

The preprocessing step is assessed according to the filtering of (i) non-expressed genes, (ii) genes expressed in some biological replicates but with zero read counts in others, (iii) genes with equal read counts among biological replicates, and (iv) genes with low counts, in this order. We considered two scenarios: the first involves simulated datasets with only fixed effects, and the second considers simulated datasets with both fixed and random effects.

Scenario 1: Count matrices of RNA-Seq read counts only with fixed effects

The preprocessing selects genes that have reads in all biological replicates. The median number of genes expressed in all biological replicates are respectively 16,434, 14,163 and 13,092 for the cases of 10%, 20% and 30% of DEGs (Fig. 2A). All genes are kept in the datasets for further identification of low-expression genes.

The genes with zero read counts among the biological replicates are removed from the analysis (Fig. 2B). The median number of these genes are respectively 12,919, 14,141 and 15,083 for the cases of 10%, 20% and 30% of DEGs. The genes expressed in some biological replicates were also identified by the preprocessing (Fig. 2C). The median number of these genes were respectively 3,675, 4,080 and 3,912 for the cases of 10%, 20% and 30% of DEGs.

Genes with the same expression across biological replicates are removed from the analysis (Fig. 2D). The median of gene percentages that are not differentially expressed differed respectively 43, 59 and 81 for the cases of 10%, 20% and 30% of DEGs. Low-expressed genes were also identified (Fig. 2E), and their detection is known to be capable of improving the inference of DEGs. The median number of the low-expressed genes were respectively 8,442, 8,713 and 9,049 for the cases of 10%, 20% and 30% of DEGs. The loss of real DEGs during preprocessing can be considered low since its observed median rate has not exceeded 6.7% (Fig. 2F). The median proportion of these lost genes were 6.6%, 6.7%, and 5.7%, respectively, related to the total number of genes in matrices with 10%, 20% and 30% of DEGs.

Scenario 2: Count matrices of RNA-Seq read counts with the insertion of random effects

Prior datasets from scenario 1 were subjected to the insertion of random effects. The predicted impact of these effects on gene expression is related to the Standard Deviations of the Normal Distribution (SDND), which were selected as the vector t = 100, 300, 600, 900, 1,200, 2,000, 3,000. The mean parameter for the estimation of the random effects was 2.028 in matrices with two replicates, 9.937 for three replicates, 7.158 for four replicates, −3.429 for five to nine replicates, and −1.615 for ten and fifteen replicates. We present in Fig. 3 the expected impact of adding a simulated normal distribution varying standard deviation with the values on $t$ and with zero mean. The higher the random effect’s standard deviation, the more significant its deleterious impact on the datasets.

The results reveal that greater SDNDs reduce the number of genes expressed in all biological replicates (Additional File 1: Fig. S2A). Also, no gene is completely unexpressed, resulting in high gene expression variation across the transcriptome (Additional File 1: Fig. S2B). In the context of gene expression among biological replicates, the median number of expressed genes is reduced, meaning that gene expression was absent in some replicates (Additional File 1: Fig. S2C). In addition, the number of genes equally expressed in all replicates decreases (Additional File 1: Fig. S2D). Interestingly, for the case of SDND 100, we detected low-expressed genes, indicating greater heterogeneity among biological replicates (Additional File 1: Fig. S2E). For each SDND, the DEGs’ loss was insignificant, thus implying that the preprocessing step filters genes without affecting true DEGs removal (Additional File 1: Fig. S2F).

The increase in the observed DEGs proportion exhibits the same pattern described before. Additional File 1, Figs. S3 and S4, respectively, show gene filters for datasets containing 20% and 30% of DEGs. Filtering has produced comparable results for the three DEG proportions, although the number of expressed genes in all replicates has decreased for matrices with 20% and 30% of DEGs compared to 10%. The loss of DEGs was equivalent for the three cases, suggesting that preprocessing has removed genes that might interfere with the subsequent inference of DEGs.

Measuring the effect of the preprocessing step on gene dispersion

Under previously mentioned scenarios, the preprocessing step was used to prepare the data for DEGs’ inference.

Scenario 1: Count matrices of RNA-Seq read counts only with the fixed effects

We estimated the dispersion of each gene and then examined how it is influenced before and after the preprocessing step (Fig. 4).

The proportions corresponding to datasets with 10% (Fig. 4A), 20% (Fig. 4B), and 30% (Fig. 4C) of DEGs have a considerable amount (33–54%) of equidispersed genes before the preprocessing step. After the preprocessing, the number of equidispersed genes was practically null. Therefore, the preprocessing has preserved the number of overdispersed genes, which is useful in properly predicting DEGs.

Scenario 2: Count matrices of RNA-Seq read counts with the insertion of random effects

When SDND is 100, a substantial and a moderate decrease is respectively observed on equidispersed and overdispersed genes in all biological replicates (Fig. 5A). However, SDND 300 (Fig. 5B), 600 (Fig. 5C), 900 (Fig. 5D), 1,200 (Fig. 5E), 2,000 (Fig. 5F), and 3,000 (Fig. 5G) has resulted in a 100% reduction in fewer biological replicates for all equidispersed genes. A similar profile was found in datasets with 20% of DEGs (Additional File 1: Fig. S5) and 30% of DEGs (Additional File 1: Fig. S6).

The preprocessing step in modeling biological variability arising from random effects

The preprocessing’s influence on biological variability was first inferred using matrices with fixed effects (Fig. 6). The biological variation in the replicates before preprocessing on matrices containing 10% of DEGs ranges from 0.46 and 0.47 (Fig. 6A). When the preprocessing step was applied, the variation was between 0.30 and 0.31. Similarly, the decreasing proportion occurs in Fig. 6B, where the BCV before and after the preprocessing step ranges, respectively, 0.48–0.51 and 0.30–0.32. On the other hand, Fig. 6C shows that the decrease rate before and after the data preprocessing step is higher than the corresponding ones in Figs. 6A and 6B. The BCV before the preprocessing step ranged from 0.57 to 0.61, while after the preprocessing, it ranged from 0.32 to 0.34.

We also examined how the BCV changes in datasets with added random effects. For the case of SDND 100 (Fig. 7A), there is a decrease in BCV in matrices with 10%, 20%, and 30% of DEGs. There were no significant differences in BCV between the three DEG proportions. SDND 300 results in a BCV increase before the preprocessing step (Fig. 7B), and they are shown below for matrices with 10%, 20%, and 30% of DEGs. For the other SDNDs studied in this work the BCV increase presents the same pattern.

Evaluation of the preprocessing step in identifying differentially expressed genes using the generalized linear models

To evaluate how the preprocessing might affect the identification of DEGs, fixed effects count matrices were submitted to GLMs with Benjamini Hochberg (glm-BH) or Bonferroni (glm-BON) P-value adjustments. In datasets with 10% of DEGs, the increasing number of biological replicates submitted to all tools resulted in a precision decrease (Fig. 8A), highlighting the fact that all tools presented a high percentage of false positives. However, the recall has improved, indicating that the false negative rate has been reduced. The evaluation metrics in Fig. 8B show that all tools have improved their performance in datasets with 20% of DEGs. A similar result is observed in datasets containing 30% DEGs (Fig. 8C). Notice that the mean accuracies remain equivalent in all tools and DEG cases.

Identification of the differentially expressed genes using the DEGRE tool

When evaluating the DEGRE performance of DEGs’ detection in matrices with the added random effects, those containing 10% of DEGs have shown precision and recall rates above those values observed in DESeq2 up to six biological replicates (Fig. 9A). The edgeR displayed the highest precision among all replicates. Still, its recall rate remained below the DEGRE with Benjamini Hochberg (DEGRE-BH) P-value correction, also up to six biological replicates. For matrices with 20% of DEGs, the DEGRE tool exhibits a high precision but below those observed in the edgeR case (Fig. 9B). However, we notice that DEGRE-BH kept its recall above the registered by the edgeR, indicating that the false negatives are more controlled by DEGRE-BH than by edgeR. For matrices with 30% of DEGs, we notice that precision and recall were above those presented by DESeq2 and edgeR, indicating that in datasets heavily impacted by random effects, DEGRE may help the detection of real DEGs (Fig. 9C). We highlight that the mean accuracy for 10% and 20% of DEGs cases are the same for all tools. However, for 30% of DEGs, the accuracy of the DEGRE-BH and DEGRE with Bonferroni (DEGRE-BON) P-value correction was above those measured on DESeq2 and edgeR tools.

Identification of the differentially expressed genes in a case study using the DEGRE tool

Here, we examined a public dataset of RNA-Seq involving healthy individuals and patients with bipolar disorder (BD), showing a BCV of 0.37. After the preprocessing step, the BCV measured was 0.32, indicating that the preprocessing has reduced the technical variation in this experiment. Of the total 47,886 genes, 29,734 were deleted, resulting in a dataset of 18,152 genes for further DEGs’ discovery. DEGRE identified 133 genes as differentially expressed candidates (Additional File 2).

We performed a GO analysis to identify the biological processes in which these 133 potential DEGs might be involved (Additional File 3). DEGRE has detected genes related to metabolic processes, with a higher frequency number in the macromolecule metabolic process (GO:0043170, n = 58) and also in the immune response, with a higher frequency number in the regulation of the immune system process (GO:0002682, n = 11). When the DEGRE tool was compared with the DESeq2 used in the previous study (Pacifico & Davis, 2017), we identified 128 different gene candidates that were differentially expressed but not noticed by the above reference (Fig. 10).

The nine genes not identified by the DEGRE tool but identified by Pacifico & Davis (2017) have Q-values above the 5% FDR. The 128 genes only identified by DEGRE as possible DEGs have an ontology that relates them to biological processes involved in patients with BD (Beech et al., 2010). An association with BD in 20.31% (n = 26) of the genes has been identified in our study (Additional File 4). Of these, the genes EIF4EBP1, SLC5A3, CYP7B1, and AHCTF1 are here discussed as examples of our method’s capability in discovering potential candidate genes linked to BD.

The EIF4EBP1 gene encodes a translation repressor protein that is the major substrate of the mTOR signaling pathway (Akbarian et al., 2020; Cha et al., 2015). This gene was found to be downregulated in patients with bipolarity (log2 fold-change = −1.15). The SLC5A3 gene is also downregulated (log2 fold-change = −1.05) and encodes the protein sodium-coupled myo-inositol transporter 1, which is related to the transport of myo-inositols for the cell (Schneider, 2015; Dai et al., 2016). The CYP7B1 gene encodes a member protein of the cytochrome P450 superfamily of enzymes (Chiang, 2004). These enzymes have functions related to the metabolism of oxysterols and steroid hormones, including neurosteroids (Chiang, 2004). DEGRE classifies the CYP7B1 gene expression profile as upregulated in patients with bipolar disorder compared to healthy individuals (log2 fold-change = 0.24). The AHCTF1 gene encodes the nucleoporin ELYS that regulates the nucleus size in mammalian cells by controlling the number of nuclear pore complexes and the ability to import molecules into the nucleus (Jevtić et al., 2019). The AHCTF1 gene was downregulated in patients with bipolarity compared to healthy individuals (log2 fold-change = −0.22).