The widespread increase in inter-individual variability of gene expression in the human brain with age

Analysis of the differentially expressed genes

First, we defined differentially expressed (DE) genes, based on the significance of the regression coefficients (FDR corrected p <= 0.05) for the linear model using the gene expression values as the dependent and age as the independent variable (see Methods). By applying two different pre-processing approaches on the data prior to DE estimation the number of DE genes was very different. SVA correction yielded 3499 DE genes, compared with 881 DE genes found by regression (Figure 1B, Supplemental Table 1). Nevertheless, 854 genes overlapped between the SVA and regression results, which constituted a quarter of genes found after SVA and 96% of the genes identified after regression correction. Quite a high overlap was consistent with the strong correlation between the expression level changes for the regression and SVA corrected data (Spearman ρ = 0.85, Supplemental Figure 3A).

To explore the biological processes affected by these changes in gene expression, gene set enrichment analysis was performed separately on the regression and SVA corrected data (Supplemental Table 2). It revealed 125 (out of 191 (regression) and 160 (SVA) categories) shared Gene Ontology Biological Process (GO BP) categories that were downregulated in the aging brain. Cognitive-function related GO terms, such as modulation of synaptic transmission, learning or memory, constituted a substantial fraction of these GO terms. In contrast, the number of the upregulated GO terms was much smaller and only 8 GO terms overlapped (out of 22 (regression) and 12 (SVA) categories) between the correction approaches, including detoxification, stress response to metal ions and cilium organization GO categories (see Supplemental Table 2).

Analysis of the differentially variable genes

Two different strategies were employed to measure change in the gene expression variability with age, namely continuous and grouped approaches. The continuous approach detects continuous monotonic change in variation from 20 to 80 years of age. The grouped approach compares the gene expression variation between two age groups: young (20 - 40 years old, N = 53) and old (60 - 80 years old, N = 22). Figure 2 illustrates the principles of these approaches and shows that the change in variability can be combined with any dynamics in the mean gene expression (upregulation, downregulation, no change). We checked if the changes in gene expression variability were confounded by the changes in gene expression level, but did not observe any relationship (Supplemental Figure 10, Fisher’s test p = 0.11, Odds ratio = 1.05).

Figure 2. Changes in gene expression and its variability with age for some individual genes, using the different approaches. Example genes are chosen that increase (A), or decrease (B) expression variability with age, when the mean gene expression either increases, does not change or decreases. The types of change, for expression and variability respectively, is shown in the parenthesis following the gene name, for each row. Genes were selected to have the biggest absolute values of Δvar(ρ) and Δvar(IQR) as well as demonstrate significant increase, decrease or no change in the expression level with age. The first column to the left illustrates mean expression level (regression-corrected) plotted against individual’s age on the x^0.25 – transformed scale. The regression line is colored in blue, with the β1 coefficient from the linear regression shown on the graph. The middle column illustrates the continuous approach to measure differential variability. Absolute values of the residuals (in red) from the regression line are plotted against age and the regression line between residuals and age (in blue) is drawn for illustrative purposes. The Spearman correlation estimates, Δvar(ρ) between the residuals and age are displayed on the graph and used in the subsequent analysis. The last column on the right illustrates the grouped approach to calculate differential variability. Gene expression levels (regression-corrected) of the individuals from the “young” (20 – 40 years old) and “old” (60 – 80 years old) groups are represented in the corresponding boxplots. A small random deviation (jitter) from the x-axis is applied for better visualization. Δvar(IQR), the fractional change in the variability in the “old” group, as compared to the “young”, is displayed on the graph.

In the continuous approach, we first fit a linear model to explain age-dependent change in expression (Figure 2, first column) and then used the residuals from this model to represent the variability. To measure change in the expression variability with age, we calculated the Spearman correlation coefficient (Δvar(ρ)) between the absolute value of residuals and age (Figure 2, middle column). The Δvar(ρ) measures ranged between -0.32 and 0.36 and were normally distributed (Shapiro-Wilk test, p > 0.05, see Methods) (Figure 3A). The distributions were significantly shifted towards positive values for both correction methods (median values range between 0.01 to 0.03, Wilcoxon test, p < 2.2e-16). Although the shift in the distribution was small, 57% to 63% percent of the genes showed increase in variability with age. However, we noted that the changes in variability calculated for each gene, using regression- and SVA-corrected data, were only weakly correlated, ρ(Δvar(ρ_regres), Δvar(ρ_SVA)) = 0.35 (Figure 3B).

Figure 3. Comparison of the change in the variability addressed using continuous and grouped approaches, regression and SVA correction. Distributions of the Δvar-measures for all the genes (n = 16675) obtained in the continuous (A) and grouped (C) approaches. Increase in the variability with age, Δvar > 0, is colored in orange, while decrease in variability, Δvar < 0, is marked in blue. The red dashed line depicts median of the distribution. The Δvar(ρ) distributions are normal with their mean and median values equal to 0.03 and 0.02 for regression and SVA, respectively; The Δvar(IQR) distributions are moderately skewed: skewness values are 0.66 and 0.68 for regression and SVA, respectively. The mean and median values of the Δvar(IQR) distribution are 0.05 and 0.02 for regression and 0.04 and 0.01 for SVA, respectively. Hexagonal heat maps illustrate relationship between regression and SVA-corrected measures of the variability for each gene, obtained in continuous - Δvar(ρ) (B) and grouped - Δvar(IQR) (D) approaches. The color gradient represents the density of the data. The linear regression line and the Spearman correlation estimate, ρ, for the corresponding variables are shown on each graph.

In the grouped approach, we first generated a distribution of expected variability in gene expression for the young individuals and treated it as a null distribution to compare with the variability from the old individuals. We used interquartile range (IQR) as a measure of variability because it is robust to outliers. In order to calculate a distribution of expected variability in the young group, we randomly selected a subsample of 22 individuals (the number of samples in the old group) from the 53 individuals in the young group 10 000 times and calculated IQR. The change in variability, Δvar (IQR), was measured as a fractional change in the IQR between old and young groups (see Methods). The p-value was determined by calculating how many times we observed value as extreme as IQR_old (see Methods). The distributions of change in variability, Δvar(IQR), were moderately skewed to the right and ranged from -0.70 up to 2.10 for the regression corrected data and from -0.78 up to 1.71 for the SVA corrected data (Figure 3C). The skew to the right was expected given that we calculate variability change as a fraction and, thus, it was more sensitive to increase in variability. In both cases, the distributions demonstrated a significant deviation from zero (Wilcox test, p-value < 2.2e-16 both for regression and SVA corrections). The data revealed that 6% and 2% more genes showed more variability in the old group, for regression and SVA approaches respectively. Similar to the continuous approach, the effect sizes calculated using regression and SVA corrected data correlated weakly ρ(Δvar(IQR_regres), Δvar(IQR_SVA)) = 0.24, Figure 3D).

Gene-level differential variability

We then asked if we could detect any genes with a significant change in variability. Using the continuous approach, we did not detect any significant change in variability with age after the multiple testing correction (Supplemental Table 3). The grouped approach leads to 741 and 746 differentially variable (DV) genes (FDR corrected p ≤ 0.05) using the regression and SVA correction, respectively (Figure 4A, Supplemental Table 4). However, the two sets of DV genes identified only have 83 genes in common (Figure 4A), one of which shows an opposite direction of change in the two sets. The correlation between Δvar (IQR) for regression and SVA corrected data is weak (ρ = 0.24), but correlation increases when we select only the common DV genes (ρ = 0.44) (Figure 4B). In agreement with our overview analysis above, we find twice as many DV genes with an increase in variability as those that decrease variability, using both correction methods: i) 533 genes increase and 208 decrease their variability in the regression correction, ii) 505 genes increase and 241 decrease their variability in the SVA correction (Figure 4A).

Figure 4. Differentially variable genes (grouped approach). (A) A bar plot of the differentially variable genes in aging identified using the grouped approach (regression, SVA correction and their overlap), direction of the change in the variability is shown in color: orange - increase in variability in aging, blue – decrease in variability (The single common gene, which shows inconsistency between approaches, is not displayed on the graph). (B) The relationship between the variability measures calculated using the grouped approach, Δvar (IQR), for regression and SVA-corrected data. Differentially variable (DV) genes identified in both corrections are highlighted in red (Spearman ρ = 0.44); DV genes identified in either regression, or SVA - in blue (Spearman ρ = 0.20); genes, that were not found to be differentially variable in any of the approaches – in grey (Spearman ρ = 0.20).

Differential variability of functional groups

Following the individual gene analysis, we explored whether genes that tend to increase or decrease variability with age are localized in particular functional groups. We performed multiple gene set enrichment analyses (GSEA) using the change in the variability with age (Δvar) measures obtained in the continuous and grouped approaches on the gene sets from KEGG and Biological Process GO categories (Supplemental Tables 5, 6). We observed no genome-level significant enrichment in particular functional groups on the data either from the continuous (SVA correction), or the grouped approach (Regression and SVA corrections). However, we found that 4 pathways, namely beta-Alanine metabolism, Ras signaling pathway, Phosphatidylinositol signaling system, Bacterial invasion of epithelial cells (FDR corrected p ≤ 0.05) were enriched among the genes showing more variability of expression in the continuous approach (Regression correction). These pathways had positive normalized enrichment scores (NES) i.e. enrichment for the genes that increase variability with age. Moreover, these pathways also had positive NES for other approaches, even though they were not significant (Supplemental Table 5).

Distribution of the DV genes in the pathways

The gene set enrichment analysis shows if there are particular gene sets that include the genes with the highest increase or decrease. Failing to detect such functional categories, we asked how the variability measures for the genes were distributed in the different functional groups of genes. For each of 310 KEGG pathways, encompassing 5922 unique genes, we analyzed the distributions of Δvar measures (Supplemental Table 7), focusing on the median value for the change in variability (Figure 5A, B). In line with the overall tendencies we observed (Figure 3A, C), the majority of pathways contained a larger number of genes that become more variable with age, irrespective of the approach or correction method used. Although the increase in variability is ubiquitous and is observed across the majority of the pathways (74-94%), the increase is small (the mean value for the shift in distributions range between 0.021 and 0.033) – in accordance with the small, but significant increase observed in the distribution for all genes. Since the pathways are not mutually exclusive, we checked if there are particular genes that are present in many different pathways and cause the shift. However, no significant correlation between the pathway membership of gene and its variability measure (Δvar) was detected (Supplemental Figure 11). We repeated the analysis using GO Biological Process categories and observed a similar trend (see Supplemental Information, Supplemental Table 7, 8).

Figure 5. Distributions of the variability measures (Δvar) obtained using a combination of continuous and grouped approaches with regression and SVA-correction for the individual pathways in KEGG database. The distribution of the variability measures (Δvar) for the genes within each pathway is represented as a box, encapsulating part of the distribution between 1^st and 3^rd quantile, median of the box is colored in black. (A) represents variability measure distributions for the continuous, while (B) for the grouped approaches. Pathways on Y-axis and corresponding them boxes are ordered by increasing median. Boxes are colored in orange if the corresponding pathways have median Δvar > 0, and in blue, if median Δvar < 0. Text label on the plot shows percentage of pathways with median Δvar > 0. Red dashed line marks Δvar = 0, while black straight-line marks median across all the pathways. The mean values for the median across all pathways are 0.033 for Δvar(ρ)_regression, 0.021 for Δvar(ρ)_SVA, 0.033 for Δvar(IQR)_regression, and 0.027 for Δvar(IQR)_SVA.

Data processing steps

Dataset selection: We utilized one of the largest age-series human brain expression datasets, featuring 269 prefrontal cortex samples from healthy individuals and spanning the whole lifespan from development (prenatal samples) through aging (80 years) [23]. These data were collected using microarray technology from people of both sexes and 4 races, namely African American (AA), Caucasian (CAUC), Hispanic (HISP) and Asian (AS). In the current analysis, we excluded fetal, childhood and early adulthood samples before the age of 20, thus limiting our sample size to 147. This was to exclude developmental processes taking place in the brain until the end of early adulthood, which exhibit discontinuous expression changes between early adulthood and aging [32]. Our main motivation was to study changes in gene expression variability during aging, considering 20 years old as a starting point.

Data characterization: The pre-processed data (loess normalization was applied on the background corrected log2 intensity ratios (sample/reference) [23]); sample and gene (probe set to Entrez gene mapping) annotations were obtained from the NCBI Gene Expression Omnibus (GEO) at accession number GSE30272. Samples were processed in 19 batches, had different quality measurements, namely pH and RNA integrity number (RIN), and differed in the time of collection after death (post-mortem interval (PMI)). Using a PCA, we found no sample outliers as judged by visual inspection of the first two principal components (Supplemental Figure 2). However, the relationship analysis between the above-mentioned factors (i.e. batch, RIN, PMI and others) and age yielded significant correlations for sex, post-mortem interval and RNA integrity, pointing to potential confounders in the data (Supplemental Figure 1). We further checked the overlap between significantly differentially variable genes in our analysis and previously reported genes that are affected by PMI and detected only a limited overlap (see Supplemental Information).

Probe set to Gene summarization: If one probe-set was mapped to several genes, it was deleted to avoid duplication. Conversely, when one gene had several probe-set expression values, they were averaged to obtain a unique gene expression value. In total 16675 genes were measured on the array.

Batch correction: To compensate for technical variation between samples, quantile normalization (QN) was performed using the ‘normalise.quantiles’ function from the ‘preprocessCore’ R library. To differentiate between the age effect and the effect of the unwanted technical and biological variability, we have applied different expression correction strategies: linear regression of the known covariates, unsupervised estimation of covariates using surrogate variable analysis (SVA) [25,33], and ComBat, a parametric empirical Bayesian framework for covariate adjustment [24]. As a result, we analyzed the same data four times, corrected using QN, QN+regression, QN+SVA, QN+ComBat. Different corrections work by adjusting for the different covariates in the linear model that explains the gene expression, namely: i) QN – no covariates were added; ii) QN+regression – 25 covariates considered: technical batches (N = 19), sex (N=2), race (N=4), post-mortem interval, RNA integrity number,pH; iii) QN + SVA – 20 surrogate variables (SV) were inferred from the expression data using the ‘sva’ function from “SVA” R library; iv) QN+ComBat: the 6 confounding factors: batch (N = 19), sex (N=2), race(N=4), PMI, RIN and pH were adjusted for, one at the time, by repeatedly applying the ComBat function from the “SVA” R library to the expression data.

Differential expression

A least squares linear regression model was used to model gene expression level change with age. Age^0.25 was used as an independent variable instead of age to account for the difference in the rate of gene expression changes between young (fast) and old (slow) as well as different density of the samples across ages. Nevertheless, the β₁- coefficients from the linear model, that uses age^0.25 correlate well with the one, that employs age (Supplemental Figure 4). Coefficients for the age covariate were used as a measure of the differential expression. P values for coefficients were adjusted using the FDR method with a threshold p ≤ 0.05 to account for multiple testing. Depending on the correction method applied, the linear model also accounted for different measured or unmeasured covariates (see Data processing steps) of the following general form:

$Y_i=\mathrm{ }{\beta }_{i0}+{\beta }_{i1}\mathrm{*}{age}^{0.25}+covariates+\mathrm{ }{\epsilon }_i$,

where $Y_i$ is the normalized log-expression level of a gene with i = 1,…,n, ${\beta }_{i0}$- intercept, ${\beta }_{i1}$- slope term and ${\epsilon }_i$ - residual (or error) term.

Differential variability

The continuous approach: First, a linear model to fit gene expression during aging, using age^0.25 and potential confounders, was constructed. Next, the Spearman correlation was calculated between the absolute values of the residuals, |ε_i| from the linear model and age. Consequently, Spearman correlation estimates were used as a measure of the change in variability, referred as Δvar_i(ρ). P values for the Spearman correlation estimates were corrected for multiple testing using FDR. FDR adjusted p ≤ 0.05 was used as a threshold to define significantly DV genes.

${\mathrm{Ä}var}_i\mathrm{ }\left(\rho \right)=\mathrm{ }Spearman\mathrm{ }\rho \mathrm{ }(\left|{\epsilon }_i\right|,\mathrm{ }age)$

The grouped approach: First, a corrected expression matrix was obtained by removing the effect of covariates (see data processing steps) from the data using the residuals from a linear regression model (Y_i = β_i0 + covariates + ε). The ‘grouped approach’ is a custom resampling-based test designed to compare gene expression variability between young (20 – 40 years old) and old (60-80 years old) groups using an interquartile range (IQR). IQR corresponds to the difference between the 75^th and 25^th percentiles of the distribution and is considered to be a robust measure of variability, meaning it is not susceptible to outliers and departure from normality in the data. In order to adjust for the unequal sample size of the young (N = 53) and old (N = 22) groups, we, first, calculated a null distribution of the IQR values for the young group by resampling it 10 000 times with the size of the old group. Next, we calculated significance as a percentage of samples where IQR_old was more extreme than IQR_young and corrected it for multiple testing using FDR correction, q ≤ 0.05. The ‘grouped’ measure of change in the variability, Δvar_i(IQR), for the gene $i,$corresponds to the difference between IQR value for the old, ${IQR}_{i old}$, and $\stackrel{-}{{IQR}_{i young}}$ (i.e. mean IQR value from the young distribution), which is then divided by the latter, see formula:

${\mathrm{Ä}var}_i\mathrm{ }\left(IQR\right)=\mathrm{ }\left({IQR}_{i\mathrm{ }old}-\mathrm{ }\stackrel{-}{{IQR}_{i\mathrm{ }young}}\right)/\stackrel{-}{{IQR}_{i\mathrm{ }young}}$

Gene Set Enrichment Analysis for KEGG pathways and GO categories

β₁ coefficients from the differential expression and Δvar measures from the differential variability analyses were used to perform gene set enrichment analysis, GSEA [34] using the “clusterProfiler” R library. KEGG pathways (N = 315) and BP GO terms (Biological processes Gene Ontology, N = 5822) with the size of between 10 and 500 genes were considered as gene sets for the GSEA.

Pathway distribution study

KEGG pathway to gene mapping was obtained from “KEGGREST” R library and pathways were pre-filtered to contain between 5 and 500 genes. As a result, 310 KEGG pathways that comprise 5922 unique genes were used for the subsequent analysis. The boxplots illustrated distributions of the Δvar measure for genes in each pathway. Pathways were sorted according to their median Δvar measure in ascending order. The percentage of pathways that have their median Δvar above zero was calculated. The analysis was replicated using BP GO terms (N = 5919) of a size between 10 and 500 genes, which in total contained 12538 unique genes. Mapping of GO terms to genes was obtained from “org.Hs.eg.db” R library.

Distribution tests

Distributions of the Δvar - measures for all the genes were tested for normality using the Shapiro-Wilk test in R (‘Shapiro.test’ function) on the multiple subsamples, consisting of 5000 measures. Skewness of the distributions was calculated using the ‘fBasics’ function from “BasicStatistics” R library.

Mean-variability relationship testing

To visualize and test if the change in gene expression variability is associated with the change in gene expression level, we plotted the difference in the means between the young and old groups against difference in the interquartile range (IQR) between the young and old groups. Mean and IQR for the old group were calculated once, while mean and IQR for the young group were calculated 10,000 times for the subsamples (see Grouped approach) and then means of the distributions of the corresponding values (mean and IQR) were used in the analysis. Fisher’s exact test was performed on the values used for the plotting.

Functional variation intolerance – variability relationship

To test if there is a relationship between age-dependent change in variability and how tolerant a gene is to functional variations, we used residual variation intolerance (RVI) percentages based on ExAC v2 data [29]. RVI percentage shows whether genes have more or less common functional genetic variation relative to the genome wide expectation. A gene with a higher percentage has more common functional variation whereas a lower percentage shows the genes that are intolerant to functional variation. We tested the association using the Spearman’s correlation coefficient between the Δvar measures and RVI percentages, for different approaches and correction strategies, separately.

Cell type-specific genes

In order to test if the changes in cell type composition can drive the changes in heterogeneity, we tested the association between cell-type specific and DV genes. First, we analyzed a cell-type specific transcriptome dataset (GSE9566) [35] which has expression profiles for FACS purified cell types in the mouse brain. Data is downloaded from GEO database [36], RMA corrected using ‘affy’ package in R [37], log2 transformed, and quantile normalized using ‘preprocessCore’ package in R [38]. Only the genes with 1 to 1 human orthologs in Ensembl Compara [39] are used. We used ‘biomaRt’ package in R [40] to retrieve the mapping between probesets and genes. When there is more than one gene that map to a single probeset, we discarded those from the analysis. If multiple probesets represent one gene, we used the mean expression value of these probesets to calculate gene expression level. To define cell-type specific genes, for each human ortholog, we first standardized the gene expression values, and calculated the effect size of each cell type (oligodendrocytes, myelinated oligodendrocytes, oligodendrocyte progenitors, astrocytes and neurons) and identified genes with an effect size higher than 2 to a particular cell type. Using this cutoff, cell type-specific gene lists did not overlap. We next tested for the overlaps between these cell type-specific and differentially variable gene lists, using Fisher’s exact test.

Software

R version 3.5.0 and “data.table” were used to perform the analyses, while “ggplot2” and “ggpubr” R libraries were used to create visualizations of the data.