Next-generation sequencing technology has been widely used in genetics and genomics studies to elucidate the biology underlying complex human diseases and traits1. RNA sequencing (RNA-Seq) technology has been widely used to profile transcriptome-wide gene expression levels and has revolutionized transcriptome analyses2,3. Differential gene expression (DGE) analysis is one approach for studying RNA-Seq data to identify genes expressed differentially with respect to a trait of interest3,4,5. Due to the cost of RNA-Seq studies and the difficulty in obtaining large numbers of relevant tissue samples from individuals, most existing DGE methods have been developed to handle small sample sizes and dichotomous traits, e.g., DESeq26, edgeR7,8, Limma9, Voom10, and MACAU11. However, with the recently reduced cost of RNA-Seq technology, RNA-Seq data from hundreds of samples have been generated for studying both dichotomous and continuous quantitative traits.

Existing DGE methods generally need to dichotomize continuous phenotypes, thus failing to account for the continuous distribution of quantitative traits6,7,8,9,11. As a result, information could be lost, and power could be reduced by not characterizing the continuous characteristics of phenotypes. This loss of information by dichotomizing a continuous biologic process (i.e., a continuous trait) may homogenize individuals together who often, in fact, lie on a continuum of disease, especially for chronic conditions of aging. For example, the cognitive manifestation of AD related dementia unfolds over years to decades12. Further, in older persons, AD dementia is often due to a combination of mixed pathologies and resilience13,14, which is better characterized by continuous cognitive traits and neuropathologic traits. During the prolonged course of this chronic disease, individuals who initially show no cognitive impairment (NCI) may manifest mild cognitive impairment (MCI) for years before they finally develop Alzheimer’s dementia as a late final manifestation (Supplemental Table 1)15. Moreover, these categories themselves are not distinct but represent stages along a progressive continuum.

Table 1 Clinical and postmortem characteristics of the discovery analytic cohort.
Full size table

The decreasing cost of RNA-Seq technology and the recognition of the importance of RNA-seq data have led to the recent availability of much larger RNA-Seq sample sizes from individuals with both dichotomous and quantitative phenotypes16,17,18. For example, the Genotype-Tissues Expression (GTEx) project V8 profiled hundreds of samples per tissue for 53 human tissues (up to n = 803 for muscle tissue)19; the CommonMind Consortium sequenced RNA from dorsolateral prefrontal cortex (DLPFC) of people with schizophrenia (n = 258) and control subjects (n = 279) for studying schizophrenia and other psychological diseases20; the prospective cohort studies of Religious Orders Study (ROS) and the Rush Memory and Aging Project (MAP) sequenced RNA from DLPFC of ~ 1200 participants for study AD traits including the continuous cognitive decline and continuous markers of AD neuropathologic changes (AD-NC), i.e., neurofibrillary tangles (NFTs) and beta amyloid (Aβ))17. Enabling DGE of continuous quantitative traits (such as cognitive decline and AD-NC traits) with hundreds of sample sizes is crucial to advance our understanding and the development of targeted treatments for complex diseases.

Recently, methods based on the standard linear regression21 and robust regression19,21 were proposed for DGE analysis of quantitative traits, which takes the quantitative trait as the response variable and the log2 transformed RNA-Seq read counts per gene as the test covariate21 (see Methods). However, the methods based on standard linear regression and robust regression models often lead to inflated false positive rates by failing to account for unknown confounders11,22,23. To improve on existing DGE methods, we apply the linear mixed model (LMM) based method as implemented by the Genome-wide Efficient Mixed Model Association (GEMMA) tool22 to conduct DGE of quantitative traits. The LMM based method can account for shared confounding factors among test samples through the sample-specific mixed effect term, which has been widely used in large-scale genetic association testing to achieve calibrated false positive rates22,23,24. The Linear Mixed Model (LMM) implemented by GEMMA employs the full-rank sample-sample correlation matrix (based on all gene expressions) to model the sample-specific random effects (see Methods), which models unknown confounding factors and thus corrects the inflated false positives that occur in linear regression models without mixed effect terms22. To demonstrate the feasibility of this approach, we developed an analytic LMM pipeline to conduct DGE, and applied the pipeline to study four cognitive and pathologic AD traits—the rate of cognitive decline and three AD-NC traits (β-amyloid, tangle density, global AD pathology burden).

We first used the LMM pipeline to conduct DGE analysis using discovery RNA-Seq data of DLPFC brain tissue (n = 632) with respect to cognitive decline and each of the AD-NC traits. Several previous studies have shown that Alzheimer Disease affects not only cognition but also non-cognitive traits such as motor functions that may be affected by the accumulation of AD-NC in tissues outside the brain25. Thus, to validate our findings with the discovery RNA-Seq data of DLPFC, we then conducted DGE analysis on the same continuous cognitive decline and AD-NC traits using additional RNA-Seq datasets from DLPFC brain tissue (n = 588) and three tissues relevant to motor functions a) supplementary motor area (SMA) within the brain (n = 234), (b) lumbar spinal cord (n = 232) neural tissues within the central nervous system (CNS) but outside the brain, and (c) muscle (n = 268), the final effector of all volitional movement, composted of non-neural tissue and located in the periphery outside the CNS. We showed that false positive rates were well calibrated by the LMM based method, compared to serious inflation of the DGE results by using the standard linear regression21, robust regression21, and Voom10. As a result, our DGE analyses by LMM identified 37 genes differentially expressed for either cognitive decline or at least one of the AD-NC traits, with 17 of those genes replicated in the additional RNA-Seq data from DLPFC, SMA, spinal cord, and muscle tissues.

Results

DGE in the discovery data of DLPFC tissue

We compared the results obtained from the DGE analyses of continuous cognitive decline and three AD-NC traits, by using four methods: LMM, standard linear regression model, robust regression, and Voom with quantitative traits dichotomized by their medians (see Methods). Our DGE analyses adjusted for various covariates, including sex, age, postmortem interval (PMI, time span between donor’s death and tissue harvest) and study group (ROS or MAP) in both models (Table 1). We constructed Quantile–Quantile plots (QQ-plots) to visualize the DGE p-values of all test genes per trait, and calculated the genomic control factor \(\lambda\).26 As depicted in Fig. 1, LMM-based test method well calibrated false positive rates (with \(\lambda \sim 1\)) for all traits (First Row in Fig. 1) with the discovery RNA-seq data of DLPFC tissue, while the standard linear regression method resulted in seriously inflated false positive rates associated (with \(\lambda >5\)) for all four traits (Second Row in Fig. 1). High inflated false positive rates were also observed in the results of all four traits with the discovery data by using the robust regression method (First Row in Supplemental Fig. 1) with \(\lambda >3\), and by using the Voom method (First Row in Supplemental Fig. 2) with \(\lambda >2\).

Figure 1
figure 1

QQ-plots and genomic control factors of DGE results by LMM (A–D) and standard linear regression model (E–H) with the discovery RNA-Seq data of DLPFC tissue of cognitive decline and three AD-NC traits.

Full size image

We identified a total of 37 potential statistically significant genes with differential expression (p-values < 0.0001 for at least one trait) by the LMM-based test method, including 2 for cognitive decline, 4 for \(\beta\)-amyloid associated genes, 2 for tangle density, and 4 for global AD pathology burden, with p-values less than the Bonferroni corrected significance threshold of 3.49 × 10–6 (Table 2; Figs. 2, 3, 4).

Table 2 LMM P-values of 37 potential DGEs (p-values < 0.0001) identified by LMM using the discovery ROS/MAP RNA-Seq data of DLPFC, for at least one trait of the cognitive decline and three AD pathologies.
Full size table
Figure 2
figure 2

Volcano plots of DGE results by LMM results by LMM with the discovery RNA-Seq data of DLPFC tissue of cognitive decline (A), tangle density (B), \(\beta\)-amyloid (C), and global AD pathology burden (D). Genes with effect size beta > 0.05 or <  − 0.05 (vertical red lines) and p-values < 0.05 (horizontal red line) were colored. Blue points were down regulated genes and red points were up regulated genes. Top five significant up and down regulated genes were labeled.

Full size image
Figure 3
figure 3

Manhattan plots of DGE p-values by LMM with the discovery RNA-Seq data of DLPFC tissue of cognitive decline (A) and tangle density (B). Top five significant DGEs were labeled. Red line indicates the significant threshold 3.49 × 10–6 with Bonferroni correction and blue line indicates p-value = 0.0001.

Full size image
Figure 4
figure 4

Manhattan plots of DGE p-values by LMM with the discovery RNA-Seq data of DLPFC tissue of \(\beta\)-amyloid (A) and global AD pathology burden (B). Top five significant DGEs were labeled. Red line indicates the significant threshold 3.49 × 10–6 with Bonferroni correction and blue line indicates p-value = 0.0001.

Full size image

Importantly, many of the potential significant genes associated with cognitive decline and AD-NC traits identified by the LMM method have been reported in prior studies. This confluence of findings bolsters the credibility of our outcomes. Notably, for example, gene MEIS3 (P-value = 1.16 × 10–8 for cognitive decline) and gene NPNT (p-value = 1.49 × 10–6 for tangle density) have previously identified as differentially expressed genes in the context of AD in earlier studies27,28. Likewise, Gene DDAH2 (p-value = 4.18 × 10–7 for global AD pathology burden) has been linked to increased levels of oxidative stress in AD brains29. Furthermore, mutations in the ALDH family genes such as ALDH6A1 (p-value = 8.36 × 10–7 for global AD pathology burden) were identified as significant risk variants for AD30. Additionally, PLCD3 (p-value = 7.64 × 10–7 for β-amyloid) is a notable protein that is cross-correlated with β-amyloid and Tau proteins in AD brains31. These results indicated the effectiveness of LMM based test for conducting DGE analyses of quantitative traits with well calibrated false positive rates. These intriguing results further underscore the efficacy of the LMM-based approach in facilitating DGE analyses involving quantitative traits, while concurrently maintaining well-calibrated false positive rates.

In summary, 8 of these 37 genes had significant LMM P-values with Bonferroni correction, and 5 of these (NPNT, ALDH6A1, DDAH2, PLCD3, MEIS3) were also reported in previous studies27,28,30. Interestingly, both NPNT and MEIS3 showed significant differential expression in cognitive decline and tangle density. We also found that 3 of these 37 genes had significant differential expression in cognitive decline and at least one of the AD pathology traits, suggesting a shared gene regulatory mechanism between cognition and AD pathology.

DGE in the replication RNA-Seq datasets

We applied the LMM, standard linear regression, robust regression, and Voom methods to conduct DGE analyses of the same cognitive decline and AD-NC traits, using additional replication RNA-Seq data of DLPFC (n = 588), SMA (n = 234), spinal cord (n = 232), and muscle (n = 268) tissues (Supplemental Tables 1–4). Confounding covariates such as sex, age, and postmortem interval were adjusted for in all analyses. Study group (ROS or MAP) was only adjusted in DLPFC datasets as participants of SMA, spinal cord, and muscle tissues are all from MAP. QQ-plots (Supplemental Figs. 1–8) still showed that LMM-based test results with these validation datasets were better calibrated than those obtained by standard linear regression, robust regression, and Voom, especially for studying the validation data of DLPFC (Supplemental Figs. 1–3). Thus, we only present the validation results obtained by LMM method here.

With validation RNA-Seq data of DLPFC, we replicated 10 of these 37 potential significant genes identified in the discovery analyses with validation p-values<\(1.35\times {10}^{-3}\) for either cognitive decline or at least one of the AD-NC traits (Table 3). For example, PLCD3 differentially expressed for \(\beta\)-Amyloid and global AD pathology was replicated with P-value = 5.42 × 10–5 for global AD pathology; TRIP6 differentially expressed for global AD pathology was replicated with p-value = 6.34 × 10–4 for global AD pathology; PLCE1 differentially expressed for \(\beta\)-Amyloid, tangle density, and global AD pathology was replicated with p-value = 4.51 × 10–4 for global AD pathology.

Table 3 LMM P-values of 10 replicated DGEs (p-value < 0.5/37 = 0.00135) using the validation ROS/MAP RNA-Seq data of DLPFC.
Full size table

Since the validation RNA-Seq datasets of the motor function related tissues (SMA, spinal cord, muscle) have sample sizes of only ~ 100 (Supplemental Table 3), we used a more liberal p-value threshold (nominal p-value < 0.05 for either cognitive decline or at least one of the AD-NC traits) to identify replicated genes. As a result, for the replicated differentially expressed genes in the CNS tissues, we found 8 in SMA with 5 overlapped in DLPFC, 2 in spinal cord, 3 in muscle with one overlapped in spinal cord, and one overlapped in SMA (Table 4). For example, ALDH6A1 differentially expressed for global AD pathology in the discovery data was replicated with P-value = 0.008 for cognitive decline in SMA and muscle. Additionally, HRSP12, a differentially expressed gene related to cognitive decline in the discovery analyses was replicated with significant p-values < 0.0001 in SMA. Interestingly, several differentially expressed genes including ADAMTS2 for cognitive decline, NPNT for all four traits, RERG for \(\beta\)-Amyloid and global AD pathology, as well as MEIS3 for cognitive decline and tangle density that were identified in the discovery data were replicated in the validation datasets of DLPFC and CNS tissues. It is noteworthy that replicated gene ADAMTS2 in both DLPFC and SMA tissues was suggested to be a therapeutic target for AD32, while replicated gene HRSP12 in SMA is also known by its alias RIDA, which was found as a GWAS risk loci for blood protein levels by previous studies33.

Table 4 LMM P-values of 11 replicated DGEs (p-value < 0.05) using the validation ROS/MAP RNA-Seq data of SMA, spinal cord and muscle tissues.
Full size table

To further illustrate the reason why differentially expressed genes in the DLPFC tissue could be replicated in the SMA, spinal cord, and muscle tissues, we created correlation heatmaps of the gene expression levels of these validated genes. For each trait, we sorted samples based on their trait values and divided sorted samples into ten equal parts (i.e., decile). For each replicated gene, we calculated the average gene expression level of the discovery DLPFC and replication tissues, for samples in each decile of the discovery and replicated traits. Then we calculated the correlations between these two vectors of average gene expression. Heatmaps of these correlations (Supplemental Figs. 9,10) show that most correlations are > 0.15, demonstrating that these replicated differentially expressed genes are not tissue specific, but likely to be shared across these motor function related tissues. For example, differentially expressed genes ADAMTS2 and HRSP12 that were replicated with all four traits in SMA all have gene expression correlations > 0.15.

In conclusion, our analysis revealed that all 5 differentially expressed genes (NPNT, ALDH6A1, RERG, PLCD3, MEIS3) with Bonferroni corrected significant p-values in the discovery data were replicated in the validation data of DLPFC tissue. Furthermore, we validated a total of 17 unique genes across the validation RNA-Seq data of DLPFC and three motor function related tissues, including 10 in DLPFC, 8 in SMA, 2 in spinal cord, and 3 in muscle. The identification of shared differentially expressed genes in two different tissue types, such as DLPFC and SMA, DLPFC and non-neural muscle, as well as spinal cord and non-neural muscle outside the brain, suggests a possible shared molecular mechanism between motor and cognition functions.

Pathway enrichment analysis

To illustrate the underlying pathways and biological functions of our identified differentially expressed genes of all 4 AD traits. We selected top 100 significantly differentially expressed genes identified by using the discovery RNA-Seq data of DLPFC tissue for each of the 4 traits to conduct pathway enrichment analyses by pathDIP34. Databases of NetPath35, Panther Pathway36, and Spike37 were used in the enrichment analyses. Significant enrichment in several biological pathways were identified with the top 100 significantly differentially expressed genes of cognitive decline and tangle density (Fig. 5). These significant pathways were reported by previous studies to be relevant with AD.

Figure 5
figure 5

Significant pathways with FDR < 0.05 that are enriched with top 100 differentially expressed genes of cognitive decline (A) and tangle density (B) with the discovery RNA-Seq data of DLPFC tissue.

Full size image

For example, for the pathways significantly enriched with top 100 differentially expressed genes of cognitive decline (Fig. 5A), the Thyrotropin releasing hormone (TRH) receptor signaling pathway (FDR = 3.54 × 10–2) has been associated with aging and neurodegenerative diseases, such as Alzheimer's disease and Parkinson's disease38. Similarly, the 5HT2 type receptor mediated signaling pathway (FDR = 4.37 × 10–2) could influence the behavioral and psychological symptoms of dementia (BPSD) in Alzheimer's disease (AD)39. The Oxytocin receptor medicated signaling pathway (FDR = 4.37 × 10–2) was suggested to be a novel protective target for vascular dementia and mixed dementias40. The toll-like receptor (TLR) signaling pathway (FDR = 1.42 × 10–2) may be involved in clearance of amyloid β-protein (Aβ) in the brain making it a potential therapeutic target for AD41. The Renin-Anigotensin System (RAS) and tumorigenesis pathway (FDR = 1.52 × 10–2) is known to play a key role in interacting with pathophysiological mechanisms of AD42. Several evidences suggest that enhancing Wnt pathway (FDR = 1.88 × 10–2) can boost synaptic function during aging, and ameliorate synaptic pathology in AD which could be novel therapeutic for restoration in the brain43. The Epidermal growth factor receptor (EGFR1) pathway (FDR = 2.01 × 10–2), a preferred target for treating memory loss induced by amyloid-beta (Aβ)44, is also enriched in β-amyloid.

Also, for the pathways significantly enriched with top 100 differentially expressed genes of tangle density (Fig. 5B), Death-Associated Protein Kinase 1 in DAPk family (FDR = 5.98 × 10–3) that plays a critical role in deregulation in AD thus manipulating DAPK1 activity and/or expression could be a promising drug target in AD45. The additional protective mechanism of AndrogenReceptor (FDR = 3.92 × 10–2) might enhance neural health and deter the progression of AD46.

Discussion

Most existing DGE analysis methods6,7,8,9,11 are developed for dichotomous traits with small sample sizes. However, with increased sample sizes in RNA-Seq datasets, there is a huge demand for methods for studying quantitative traits in population-based RNA-Seq studies. As shown by the MACAU11 method paper, incorporating a mixed term into DGE analysis can help to control for false positive rates in RNA-Seq studies. In this study, we develop an analytic pipeline for implementing the GEMMA tool22, enabling the DGE analysis of quantitative traits by LMM, and apply to real ROS/MAP RNA-Seq datasets of DLPFC, SMA, spinal cord, and muscle tissues for studying continuous cognitive decline and AD-NC traits. The pipeline is freely available from https://github.com/tangjiji199645/LMM_DGE_Pipeline.

Our application studies found that DGE analyses results obtained by LMM-based tests were all well calibrated for false positive rate, especially in our discovery RNA-Seq data of DLPFC, while the DGE results obtained by the alternative standard linear regression, robust regression, and Voom methods all have inflated false positive rates. A list of 37 potential differentially expressed genes were identified by LMM in the discovery data, and 17 of these were replicated in the additional RNA-seq data of DLPFC, SMA, spinal cord, and muscle tissues. AD relevant biological pathways were also found to be enriched with top differentially expressed genes of cognitive decline and tangle density.

However, the LMM-based test still has several limitations for studying RNA-Seq data. First, the LMM assumes normal distributions for both the response variable and covariates, while the raw RNA-Seq data are read counts per gene. Even after log2 transformation for the raw RNA-Seq read counts, the transformed quantitative gene expression level per gene may not be normally distributed47, which could lead to a biased estimation of effect sizes. One might apply both the MACAU (mixed Poisson model-based test that is suitable for modeling RNA-Seq read counts) and our LMM analytic pipeline for DGE analysis and examine the results by QQ-plots. Second, the effect of the gene expression on the quantitative trait of interest may be heterogeneous, with effect sizes varying across the quantiles of the quantitative trait. The LMM-based method only tests the association between gene expression and quantitative trait of interest in expectation, ignoring the possible heterogeneous effects on different quantiles of the quantitative trait48. Therefore, further studies are needed to develop a DGE method based on quantile regression with a mixed effect term to account for the possible heterogeneous effect of gene expression across all the quantiles of the quantitative trait of interest, while controlling for false positive rates.

Overall, we provide a useful LMM pipeline for conducting DGE analysis with quantitative traits and large sample sizes, which is shown well calibrating for false positive rates in real studies. Our real application studies not only demonstrated the effectiveness of the LMM approach for DGE analysis, but also identified a list of differentially expressed genes for cognitive decline and AD-NC traits in DLPFC that were validated in DLPFC, SMA, spinal cord, and muscle tissues. Our findings have important implications for understanding the underlying biological mechanisms of the continuous AD traits of cognitive decline and neuropathologic changes, and may provide insights into the development of new therapeutic approaches for AD.

Methods

RNA-Seq data normalization

Preprocessing and normalization of raw read counts is a critical step in DGE analysis. Generally, samples with total mapped reads < 10 million are suggested to be excluded, and genes with expression levels < 0.1 transcript per million (TPM) in > 20% samples are also suggested to be excluded. We use DESeq26 to normalize raw RNA-Seq data.

Let \({{\varvec{K}}}_{{\varvec{i}}{\varvec{j}}}\) denote the read count for gene j and sample i, following a negative binomial distribution with mean \({{\varvec{\mu}}}_{{\varvec{i}}{\varvec{j}}}\) and dispersion \({\boldsymbol{\alpha }}_{{\varvec{j}}}\) given by the following formulas:

$${\varvec{K}}_{{{\varvec{ij}}}} \sim {\varvec{NB}}\left( {{\varvec{\mu}}_{{{\varvec{ij}}}} ,{\varvec{\alpha}}_{{\varvec{j}}} } \right),$$
(1)
$${\varvec{\mu}}_{{{\varvec{ij}}}} = {\varvec{s}}_{{{\varvec{ij}}}} {\varvec{x}}_{{{\varvec{ij}}}} \user2{ }\user2{.}$$
(2)

The normalization factor \({{\varvec{s}}}_{{\varvec{i}}{\varvec{j}}}\) is assumed to be shared per sample, \({{\varvec{s}}}_{{\varvec{i}}{\varvec{j}}}={{\varvec{s}}}_{{\varvec{i}}}\), and \({{\varvec{s}}}_{{\varvec{i}}}\) is estimated by the median (across all genes) of the ratios of raw read count per gene and its corresponding geometric mean \({{\varvec{K}}}_{{\varvec{i}}}^{{\varvec{R}}}\) (across all samples) as in the following formula49,50:

$${\varvec{s}}_{{\varvec{i}}} = {\varvec{median}}_{{{\varvec{j}}:{\varvec{K}}_{{\varvec{j}}}^{{\varvec{R}}} \ne 0}} \frac{{{\varvec{K}}_{{{\varvec{ij}}}} }}{{{\varvec{K}}_{{\varvec{j}}}^{{\varvec{R}}} }}\;with\;{\varvec{K}}_{{\varvec{i}}}^{{\varvec{R}}} = \left( {\mathop \prod \limits_{{{\varvec{i}} = 1}}^{{\varvec{m}}} {\varvec{K}}_{{{\varvec{ij}}}} } \right)^{{1/{\varvec{m}}}} .$$
(3)

Normalized read counts \({{\varvec{x}}}_{{\varvec{i}}{\varvec{j}}}\) are given by \(\frac{{K}_{ij}}{{s}_{ij}}\) and then log2 transformed with an offset of 1 and taken as the test variable in the LMM, standard linear regression, robust regression, and Voom methods.

Standard linear regression model for DGE analysis of quantitative traits

As proposed by previous study21, testing if gene \(j\) with expression levels \({{\varvec{X}}}_{{\varvec{j}}}\) (normalized and log2 transformed) is differentially expressed with respect to a quantitative trait \({{\varvec{Y}}}_{{\varvec{n}}\times 1}\) can be done based on the following standard linear regression :

$${\varvec{Y}} = \user2{ W\alpha } + \user2{ X}_{{\varvec{j}}} \beta_{j} \user2{ } + \user2{ \varepsilon },\user2{ }$$
(4)
$$\user2{\varepsilon }\sim \user2{ MVN}\left( {0,{\varvec{I}}_{{\varvec{n}}} \sigma^{2} } \right),$$
(5)

where n is the number of test samples; \({\varvec{W}}\) is a \({\varvec{n}}\times {\varvec{c}}\) covariate matrix; \(\boldsymbol{\alpha }\) is a \({\varvec{c}}\times 1\) vector of covariate effects including the intercept; \({{\varvec{\beta}}}_{{\varvec{j}}}\) is the effect size of gene j; and \({\varvec{\epsilon}}\) denotes the error term following a Multivariate Normal distribution (MVN). The DGE analysis is to test the null hypothesis of \({H}_{0}: {\beta }_{j}=0\boldsymbol{ }vs. {H}_{a}:{\beta }_{j}\ne 0\), which can be conducted by using the Wald test statistic,

$$\frac{\widehat{{\beta }_{j}}}{se(\widehat{{\beta }_{j}})}\sim N(0, 1) \,{\text u \text n\text d\text e \text r}\, {H}_{0},$$

with the maximizing likelihood estimator \(\widehat{{\beta }_{j}}\) and its standard error \(se(\widehat{{\beta }_{j}})\).

Robust regression for DGE analysis of quantitative traits

Robust regression21 assumes the same model as the standard linear regression model (4), yet it furnishes robust coefficient estimates when the test samples contain influential outliers that could heavily impact standard linear regression estimates. Different from the standard linear regression method where each sample contributes equally to the ordinary least squared estimation of regression coefficients, robust regression incorporates Huber’s M-estimation51 that is obtained by minimizing the following objection function through a numerical method called iteratively reweighted least squares (IRLS):

$$\mathop \sum \limits_{i = 1}^{n} \rho \left( {Y_{i} - W_{i} \alpha - X_{i,j} \beta_{j} } \right),\; with\; \rho \left( e \right) = \left\{ {\begin{array}{*{20}l} {\frac{1}{2}e^{2} } \hfill & { if \left| e \right| < k} \hfill \\ {k\left| e \right| - \frac{1}{2}k^{2} } \hfill & {if \left| e \right| \ge k} \hfill \\ \end{array} } \right., \;w\left( e \right) = \left\{ {\begin{array}{*{20}c} 1 & {if \left| e \right| \le k} \\ {\frac{k}{\left| e \right|}} & { if \left| e \right| \ge k} \\ \end{array} } \right. ,$$

where \(k\) is the tuning constant (often taken as 1.345 \(\sigma\)); \(\rho (e)\) is Huber’s objective function; and \(w(e)\) is Huber’s weighted function. The weighted least squares estimate with sample weights given by \(w({Y}_{i}-{W}_{i}\widehat{\alpha }-{X}_{i,j}\widehat{{\beta }_{j}})\) will be iteratively calculated given coefficient estimates from last iteration, until a stopping criterion is met.

As described in the previous study that proposed using the robust regression for DGE analysis21, the R packages of “rlm” and “sfsmisc” can be used to conduct the statistical test of \({H}_{0}: {\beta }_{j}=0\boldsymbol{ }vs. {H}_{a}:{\beta }_{j}\ne 0\). The robust regression method is expected to provide more reliable estimates of \({\beta }_{j}\) when outliers are present in the test samples.

Voom for DGE analysis

Since the Voom method10 is developed for detecting differentially expressed genes between two or more conditions, the quantitative trait \({{\varvec{Y}}}_{{\varvec{n}}\times 1}\) needs to be dichotomized for testing if gene \(j\) with expression levels \({{\varvec{X}}}_{{\varvec{j}}}\) (normalized and log2 transformed read counts) is differentially expressed. Generally, the quantitative trait is dichotomized to \({{\varvec{Y}}}_{{\varvec{n}}\times 1}^{{\varvec{d}}}\) by taking the median as a cut-off. That is, if \({{\varvec{Y}}}_{{\varvec{i}}}^{{\varvec{d}}}\) is greater than the median, it is assigned a value of 1, otherwise \({{\varvec{Y}}}_{{\varvec{i}}}^{{\varvec{d}}}\) is assigned a value of 0. The Voom method assumes the following model:

$${\varvec{E}}\left( {{\varvec{X}}_{{{\varvec{ij}}}} } \right) = \user2{ W}_{{\varvec{i}}} \user2{\alpha } + \user2{ Y}_{{\varvec{i}}}^{{\varvec{d}}} {\varvec{\beta}}_{{\varvec{j}}} \user2{ }$$

where \({\varvec{W}}\) is a \({\varvec{n}}\times {\varvec{c}}\) matrix of confounding covariates; \(\boldsymbol{\alpha }\) is a \({\varvec{c}}\times 1\) vector of covariate effects including an intercept term; \({\beta }_{j}\) is coefficient of gene j representing log2-fold-changes between two conditions of the dichotomous trait. The same hypothesis with \({H}_{0}: {\beta }_{j}=0\boldsymbol{ }\; vs. {H}_{a}:{\beta }_{j}\ne 0\) is tested by Voom. Different from the ordinary least squared estimates based on (10), the Voom method robustly estimates the mean–variance relationship of the log2 transformed read counts, generates a precision weight for each sample, and enters these into the limma empirical Bayes analysis pipeline10.

LMM by GEMMA

To test if gene \(j\) with expression levels \({{\varvec{X}}}_{{\varvec{j}}}\) (normalized and log2 transformed read counts) is differentially expressed with respect to a quantitative trait \({{\varvec{Y}}}_{{\varvec{n}}\times 1}\) with \(n\) samples in DGE analysis, the following LMM is assumed:

$${\varvec{Y}} = \user2{ W\alpha } + \user2{ X}_{{\varvec{j}}} \beta_{j} \user2{ } + \user2{ Zu } + \user2{ \varepsilon },\user2{ }j = 1, ..., p$$
(6)
$${\varvec{u}}\sim {\varvec{MVN}}\left( {0,\user2{\gamma \tau }^{ - 1} {\varvec{M}}} \right),$$
(7)
$${\varvec{\varepsilon}}\sim {\varvec{MVN}}\left( {0,{\varvec{\tau}}^{ - 1} {\varvec{I}}_{{\varvec{n}}} } \right),$$
(8)

where \({\varvec{W}}\) is a \({\varvec{n}}\times {\varvec{c}}\) matrix of confounding covariates; \(\boldsymbol{\alpha }\) is a \({\varvec{c}}\times 1\) vector of covariate effects including an intercept term; \({\beta }_{j}\) is the effect size of gene j; Z is a \({\varvec{n}}\times {\varvec{n}}\) loading matrix which is taken as an identify matrix for DGE analysis; \({\varvec{\varepsilon}}\) is a \({\varvec{n}}\times 1\) vector of independent errors following a normal distribution with mean 0 and variance \({{\varvec{\tau}}}^{-1}\); u is a \({\varvec{n}}\times 1\) vector denoting random effects of all samples following a Multivariate Normal distribution with mean 0 and variance–covariance matrix \({\varvec{\gamma}}{{\varvec{\tau}}}^{-1}{\varvec{M}}\); \({\varvec{\gamma}}\) is the ratio of variance components between random effects and errors; \({\varvec{M}}\) is a \({\varvec{n}}\times {\varvec{n}}\) sample-sample correlation matrix with all gene expressions; and \({{\varvec{I}}}_{{\varvec{n}}}\) is an \({\varvec{n}}\times {\varvec{n}}\) identity matrix.

GEMMA tool22 is a C++ programed software that can be used to conduct tens of thousands of Wald test for \({H}_{0}:{\beta }_{j}=0 (j=1, \dots , p)\). Under the above LMM, GEMMA efficiently calculates the REstricted Maximum Likelihood (REML) estimates of \({\varvec{\gamma}},\boldsymbol{ }\boldsymbol{ }\boldsymbol{ }{\beta }_{j},\) and the standard error of \({\beta }_{j}\). Although GEMMA is originally developed for genome-wide association study to test the association between a genetic variant and a quantitative trait, it can be applied to DGE if both genetic variant and log2 transformed gene expression follow normal distributions.

We demonstrate the feasibility and effectiveness of conducting DGE analyses using the LMM method through application studies with the ROS/MAP data17. Our analyses were conducted in two steps. First, we normalized raw RNA-Seq data using DESeq26. Second, we tested DGE using the LMM method as implemented in the GEMMA tool22.

ROS/MAP data

The Religious Order Study (ROS) and the Rush Memory and Aging Project (MAP) are two prospective community-based harmonized cohort studies of aging, which recruit senior individuals without known dementia at study entry52. All ROSMAP participants agree to structured annual clinical testing and autopsy and brain donation upon their death. RNA-Seq data (DLPFC, SMA, spinal cord, and muscle tissues) and AD pathologies were profiled from decedents. Both studies were approved by the Institutional Review Board of Rush University Medical Center, and all participants signed informed and repository consents and an Anatomic Gift Act.

Alzheimer’s disease clinical and neuropathologic traits

Our study focused on cognitive decline and three AD-NC traits, including β-amyloid, tangle density, and global AD pathology burden. The cognitive decline (annual rate of cognitive decline) is the estimated person-specific rate of change in the global cognition variable over all follow-ups generated by a mixed effects model53,54. The AD-NC trait of tangle density was quantified using molecularly specific immunohistochemistry. It was profiled as the average PHFtau tangle density within two or more 20 µm sections from eight brain regions—hippocampus, entorhinal cortex, midfrontal cortex, inferior temporal, angular gyrus, calcarine cortex, anterior cingulate cortex, and superior frontal cortex. These two are identified by molecularly specific immunohistochemistry. Trait β-Amyloid quantifies the average percent area of cortex occupied by β-Amyloid protein in adjacent sections from the same eight brain regions. The global AD pathology burden is a quantitative summary of AD pathology derived from counts of three AD pathologies: neuritic plaques, diffuse plaques, and neurofibrillary tangles with 5 brain regions midfrontal cortex, midtemporal cortex, inferior parietal cortex, entorhinal cortex, and hippocampus (total 15 regional counts). To improve normality, these three quantitative AD pathology traits were transformed by taking the square root. Futher details have been previously reported55,56.

RNA-Seq data of DLPFC, SMA, spinal cord, and muscle tissues

RNA-Seq data were profiled from deceased ROS/MAP participants for DLPFC tissue within the brain (n = 1220, n = 632 as discovery data Table 1, n = 588 as validation data; Supplemental Table 1) and three validation tissues (Supplemental Tables 2–4)––SMA in the brain (n = 234), contralateral ventral horn in the lumbar spinal cord (outside the brain, n = 232), and non-neural quadriceps muscle ipsilateral to the ventral horn (outside the brain, n = 268). The raw RNA-Seq fastq data were first aligned to the reference human genome and then quantified by the number of reads mapped to gene regions. Raw read counts were first normalized and log2 transformed by DESeq2, and then used as the test gene expression covariates in the LMM model. The ROS/MAP RNA-Seq data of DLPFC (n = 632) were analyzed as discovery data, while RNA-Seq datasets another 588 DLPFC samples, and samples of SMA, spinal cord, and muscle tissues were analyzed as validation data. Participants are not overlapped between the DLPFC samples and samples of SMA, spinal cord, and muscle tissues, while participants of RNA-Seq data of SMA, spinal cord, and muscle tissues are largely overlapped. Technical details of RNA-Seq data profiling can be found in the Supplemental Text.

Ethics declarations

The Religious Order Study (ROS) and the Rush Memory and Aging Project (MAP) were approved by the Institutional Review Board of Rush University Medical Center, and all participants signed informed and repository consents and an Anatomic Gift Act. All data analyzed in this study were de-identified which are not considered as human data according to NIH protocols. We confirm that all analytical methods were performed in accordance with the relevant guidelines and regulations.