Age-dependent co-dependency structure of biomarkers in the general population of the United States

National Health and Nutrition Examination Survey (NHANES) participants’ biomarkers are significantly correlated with age

We analyzed 50 clinical biomarkers on 27,508 participant samples from the National Health and Nutrition Examination Survey (NHANES), a dataset constituted of non-institutionalized individuals [13] aged between 20 years old and up to 80 years old (excluded), with a mean age of 46.3 years old (Tables 1 and 2). The distribution of the ages can be found in Figure S1 (quantiles: 0%:20, 25%:32, 50%:45, 75%:60, 100%:79). We found that that 92% of biomarkers are significantly associated (Bonferroni-corrected p-value < 0.05) with age (Figure 2 and Table 3).

Figure 2.

Changes in biomarkers with age. (A) Systolic blood pressure. Example of a biomarker that increases with age. (B) Upper leg length. Example of a biomarker that significantly decreases with age. (See discussion for explanation about the generational effect). (C) Volcano plot of the significance of the changes associated with age for the 50 biomarkers. The black horizontal line corresponds to the 0.05 significant threshold, after Bonferroni correction.

Age-dependent correlations between biomarkers reveal structure changes with age

Next, we estimated the overall correlation structure of the biomarkers ((Figure 3AB). The mean of the absolute values of the correlations is modest (0.028; SE: 0.11).

Figure 3.

Baseline correlation structure of the biomarkers. (A) Correlation matrix of the biomarkers.The biomarkers are ordered based on their clustering. Red is associated with a positive correlation, blue with a negative correlation. (B) Hierarchical clustering of the biomarkers. For each cluster, the number in green is the bootstrap probability (BP)--the percentage of bootstraps in which the cluster was present. The number in red is called the approximated-unbiased p-value (AU). AU is a better estimation of the unbiased p-value than BP, and the red boxes circle the significant clusters, based on this criterion, with alpha=0.95 (a cluster is marked as significant if its AU is greater than 95). The number in grey is the rank of the cluster, low numbers means the clustering happened early in the process. The height is the measure of the proximity between the two clusters being merged. The height is one minus the mean correlation between the two clusters, so two perfectly correlated biomarkers/clusters cluster at height zero, and two perfectly uncorrelated biomarkers/clusters cluster at height one. The first column compares males and females, the second column compares non-Hispanic Whites and Hispanics, the third column compares non-Hispanic Whites and non-Hispanic Blacks, and the fourth column compares the two control groups.

We found that 419 out of the 1,225 (33.2%) pairwise correlations show significant (Bonferroni-corrected p-value < 0.05) change with age after correcting for multiple testing (Figure 4). We classified these correlations with 6 different and general trends, as illustrated in Figures 4 and 5. 316 (25.8%) correlations significantly decreased with age (Figure 5A-C). They (1) start positive among young individuals and end negative among old individuals (189 findings [15.4%], such as standing height vs. systolic blood pressure, Figure 5A). Second, they (2) start positive and remain positive (123 findings [10.0%], such as hematocrit vs. total calcium, Figure 5B), Third, they (3) start negative and remain negative (four findings [0.3%], such as percentage of segmented neutrophils vs. number of lymphocytes, Figure 5C).

Figure 4.

Volcano plot of correlations changes with age. The black horizontal line corresponds to the Bonferroni-corrected significance threshold (0.05). The green dots are correlations that decrease with age but remain positive. The blue dots are correlations that decrease with age, starting positive and ending negative. The purple dots are correlations that decrease with age and remain negative. The yellow plots are correlations that increase with age but remain negative. The red dots are correlations that increase with age, starting negative and ending positive. The orange dots are correlations that increase with age and remain positive.

Figure 5.

Examples of different trajectories for changes in correlations with age. (A) Hematocrit vs. total calcium. Example of a correlation that decreases with age but remains positive. (B) Standing height vs. systolic blood pressure. Example of a correlation that decreases with age, starts positive and ends negative. In red, we replicated the analysis on the UK Biobank cohort. (C) Percentage of segmented neutrophils vs. number of lymphocytes. Example of a correlation that decreases with age and remains negative. (D) Hematocrit vs. chloride. Example of a correlation that increases with age but remains negative. (E) Percentage of segmented neutrophils vs. creatinine. Example of a correlation that increases with age, starts negative and ends positive. (F) Chloride vs. sodium. Example of a correlation that increases with age and remains positive.

Inversely, we found 103 correlations (8.4%) that increase with age (Figure 5D-F). They (1) start negative and end up positive (65 findings [5.3%], such as percentage of segmented neutrophils vs. creatinine, Figure 5D), (2) start negative and remain negative (22 findings [1.8%], such as hematocrit vs. chloride, Figure 5E), or (3) start positive and remain positive (16 findings [1.3%], such as chloride vs. sodium, Figure 5F). A summary table of the different age trajectories can be found in Table 4.

We successfully replicated our top finding (standing height vs. systolic blood pressure) using 470,895 samples from the UK Biobank cohort [14] (Figure 5B). The correlation decreases by 0.011 unit/year in the NHANES participants (0.39 at age 20, -0.28 at age 80), and by 0.009 unit/year on the UK Biobank dataset (0.23 at age 40, -0.06 at age 71).

We then tested for a global change in correlation amongst all 1,225 correlations. We tested for a change in the mean of the absolute values of the correlations with age (Figure 6). We confirmed a significant decrease in the means of the correlations with age, both when considering all the biomarkers (coefficient=-4.4x10^-4 unit/year, p-value=2.4x10^-20, Figure 6A), and those correlations that significantly changed with age only (coefficient=-1.5x10^-3 unit/year, p-value=2.6x10^-56, Figure 6B). Interestingly, we observed a non-linear trend. The mean absolute value of the correlations seems to decrease until the age of 50 and remain relatively stable after age 50. This suggests a non-linearity in the rate of aging.

Figure 6.

Correlations: Change of the means of the absolute values of the correlations with age. (A) Change of the mean of the absolute value of all the correlations with age. (B) Change of the mean of the absolute value of the correlations with age, limiting the analysis to correlations which significantly change with age.

We executed our pipeline on different sexes and ethnicities to compare the age-dependent correlation architecture between these groups (e.g., males versus females). We first established a baseline by comparing the correlation structure of the different demographics on the full age range (Figure 7).

Figure 7.

Correlations: Baseline differences between demographics. (A) Heatmap visualization of correlation structure. The upper left triangle shows the correlation structure for the first group of the comparison (e.g. males) while the lower right triangle of the matrix shows the difference between the second group and the first group of the comparison (e.g. difference between females and males). (B) Correlation between the 1,225 correlations of the first group (e.g. males) and the 1,225 correlations of the second group (e.g. females). The diagonal black line represents a perfect correlation. The further away from this line the points lie, the bigger the difference between the two groups, and the lower the correlation coefficient. (C) Volcano plots of the 1,225 differences in correlations between the groups. The horizontal black lines represent the threshold of significance of 0.05 for the Bonferroni corrected p-values. The vertical axis is not shared between the plots.

296 (24.2%) correlations are significantly different between males and females. 190 (64% of 296) are higher for females than males. For example, the correlation between alkaline phosphatase and waist circumference is 0.32 for females versus 0.04 for males (p-value= 1.6x10^-70), and the correlation between mean cell hemoglobin and red blood cell count is -0.40 for females and -0.54 for males (p-value=2.5x10^-13). In contrast, 106 correlations (36% of 296) are higher for males than females. For example, the correlation between mean cell hemoglobin and red cell distribution width is -0.26 for males and -0.55 for females (p-value=5.0x10^-51), and the correlation between standing height and weight is 0.44 for males and 0.28 females (p-value=4.0x10^-25).

We also detected differences between ethnicities (non-Hispanic Whites versus Hispanics: 58 findings, 25 larger for Hispanics, 33 larger for non-Hispanic Whites, Figure 7 columns 2. non-Hispanic Whites versus non-Hispanic Blacks: 113 findings, 26 larger for non-Hispanic Blacks, 87 larger for non-Hispanic Whites, Figure 7 column 3.). Those differences were stronger than the differences between the controls (2 spurious findings, Figure 7 column 4), which further confirms their significance.

We analyzed which correlations trend significantly differently with age in different demographics (Figure 8). We detected 36 significant sex differences in the age trajectories of correlations (Figure 8B column 1). 20 of them show an increase with age (1.6%). For example, cholesterol vs. albumin has a difference in correlation trajectory (p-value=4.0x10^-13). In females it increases from -0.22 at age 20 to 0.31 at age 80; in contrast, males exhibit no significant change in correlation (0.12 on average). 16 of the correlations’ differences showed a decrease with age (1.3%). For example, blood urea nitrogen vs. albumin has a difference in correlation trajectory (p-value=7.9x10^-9). In females it decreases from 0.38 at age 20 to -0.15 at age 80; in contrast, males exhibit no significant change in correlation (0.06 on average).

Figure 8.

Correlations: Differences in age trajectories between demographics. (A) Correlations between the rates at which the correlations are changing in the two compared groups, only taking into account the correlations for which significant changes were detected in both compared groups (e.g. males and females). (B) Volcano plots showing differences in the rates of change of correlations with age between different demographic groups. The horizontal black lines represent the threshold of significance of 0.05 for the Bonferroni corrected p-values.

We detected few differences in age trajectories of correlations between ethnicities. Specifically, we detected three differences in correlations’ trajectories between non-Hispanic Whites and Hispanics (Figure 8 columns 2) and eight differences between non-Hispanic Whites and non-Hispanic Blacks (Figure 8 columns 3).

The comparison between the two control groups, on the other hand, did not exhibit any differences. This suggests that the findings for the comparisons between the sexes and the ethnic groups are not spurious. A summary table of the trajectories of the correlations in different demographic groups can be found in Table 5.

Finally, we detected a difference between sexes in the rate at which the mean of the absolute values of the correlations that significantly change with age decreases. It decreases 37% faster for males (from 0.214 at age 20 to 0.165 at age 80) than for females (from 0.198 at age 20 to 0.162 at age 80).

Predictability of biomarkers is strongly age dependent

Next, we sought to estimate how predictions of biomarkers change with age. We hypothesized that first, biomarkers are predictive of other biomarkers, and second, that these predictions significantly change with age. We established baseline results using the full cohort (Figure 9A) by predicting each of the biomarkers using all the other biomarkers. The mean R² was 0.618, the standard deviation 0.294, the minimum 0.120 and the maximum 0.999.

Figure 9.

Predictabilities: Baseline predictabilities for the 50 biomarkers and significance of their changes with age. (A) Histogram of the predictabilities.Histogram of the predictabilities (R²) obtained on the 50 biomarkers using the full 27,508 samples. (B) Volcano plot of the changes in predictabilities with age.The horizontal black lines represent the threshold of significance of 0.05 for the Bonferroni corrected p-values.

Then, we tested for changes with age and found that the prediction accuracy of 22 (44%) of the biomarkers significantly change with age (Figure 9B). For example, we found that 55% of the variance of albumin can be explained by other biomarkers in people at age 20, whereas 0.00% of the variance can be explained when people are 80 (Figure 10A1). Inversely, we found that only 11% of the variance of serum glucose levels can be explained by other biomarkers at age 20, whereas 67% of it can be explained at age 80 (Figure 10A2). We found that the average predictability decreases by 8.7x10^-4 unit/year over all the biomarkers year over year (p-value=0.03) (Figure 10B1).

Figure 10.

Changes of the predictabilities (R²) with age. (A) Predictabilities: Examples of different age trajectories. (A1) Albumin. Example of a biomarker whose predictability decreases with age. (A2) Glucose, serum. Example of a biomarker whose predictability increases with age. (B) Change of the means of the predictabilities with age. (B1) Change of the mean of all the predictabilities with age. (B2) Change of the mean of the predictabilities with age, limiting the analysis to biomarkers whose predictability significantly change with age.

We investigated the difference in predictabilities between sexes and ethnicities by running our pipeline on different demographic groups. First, we established baseline differences, looking at the full age range (Figure 11A and 11B). We detected 15 differences between sexes. 8 biomarkers show higher predictability in males versus females, such as alanine aminotransferase ALT (R²=0.70 for males, R²=0.62 for females, p-value=1.5x10^-6). In contrast, 7 biomarkers are better predicted in females than males, such as red cell distribution width (R²=0.25 for males, R²=0.41 for females, p-value=7.2x10^-23).

Figure 11.

Predictabilities: Differences between demographics. (A) Correlation between the predictabilities of the 50 biomarkers in the first group and the predictabilities of the 50 biomarkers in the second group. The diagonal black line represents a perfect correlation. The further away from this line the points lie, the bigger the difference between the two groups, and the lower the correlation. (B) Volcano plots reporting the significance and the size of the baseline (calculated on the full age range 20-80) differences between the groups. The horizontal black lines represent the threshold of significance of 0.05 for the Bonferroni corrected p-values. (C) Volcano plots showing which predictabilities are changing at significantly different rates with age in different demographic groups. The horizontal black lines represent the threshold of significance of 0.05 for the Bonferroni corrected p-values.

We detected differences in predictabilities between ethnicities, as well. 3 show higher predictability in non-Hispanic Whites versus Hispanics, such as blood urea nitrogen (R²=0.76 for non-Hispanic Whites, R²=0.65 for Hispanics, p-value=3.2x10^-5), and 9 are better predicted for Hispanics, such as glycohemoglobin (R²=0.49 for non-Hispanic Whites, R²=0.65 for Hispanics, p-value 1.6x10^-10). 6 were better predicted for non-Hispanic Whites than for non-Hispanic Blacks, such as Blood urea nitrogen (R²=0.76 for non-Hispanic Whites, R²=0.68 for non-Hispanic Blacks, p-value=0.01), and 10 were better predicted for non-Hispanic Blacks, such as red cell distribution width (R²=0.25 for non-Hispanic Whites, R²=0.41 for non-Hispanic Blacks, p-value=2.0x10^-13). There were no significant findings detected when comparing the controls, confirming the signal of findings above.

Next, we tested for significant differences between demographic groups in the trajectories of predictabilities with age (Figure 11C). We detected different age trajectories of the predictabilities between males versus females. Four of them (10%) showed an increase of the signed difference between the R²s values with age, such as alanine aminotransferase. For males, the R² value decreases from 0.65 to 0.45 whereas for females, there was no significant change in R² value (0.47 on average). In contrast, in females, albumin has an R² value of 0.55 at age 20 and 0.05 at age 80; however, in males, the R² value does not significantly change (0.22 on average). We detected fewer differences between ethnicities (three differences between non-Hispanic Whites and Hispanics, and one difference between non-Hispanic Whites and non-Hispanic Blacks).

Features that drive prediction of biomarkers change with age

Finally, we investigated what exact biomarker variables influence the predictions for every age. We first estimated these coefficients on the whole cohort (Figure 12A), then examined what coefficients change with age (Figures 12B, 13). We found that only 40 out of the 2,450 (1.6%) coefficients significantly change with age. 20 (0.8%) increase and 20 (0.8%) decrease. For example, osmolality and sodium become weaker predictors of blood urea nitrogen with age (Figure 13A and 13D), while blood urea nitrogen becomes a stronger predictor of both osmolality and sodium with age (Figure 13C and 13F). We also observed that predictor selection is affected by age. For example, sodium and urine albumin are not selected as predictors of creatinine for young people but become good enough predictors to be selected at older age (Figure 13B and 13E).

Figure 12.

Regression coefficients: Baseline regression coefficients for the prediction of the 50 biomarkers and significance of their changes with age. (A) Heatmap of the regression coefficients on the full cohort. Each column corresponds to the 49 coefficients used to predict the values for one of the 50 biomarkers. (B) Volcano plot of the changes of the regression coefficients with age.The black horizontal line corresponds to the 0.05 significant threshold, after Bonferroni correction. The green dots are coefficients that decrease with age but remain positive. The blue dots are coefficients that decrease with age, starting positive and ending negative. The purple dots are coefficients that decrease with age and remain negative. The yellow plots are coefficients that increase with age but remain negative. The red dots are coefficients that increase with age, starting negative and ending positive. The orange dots are coefficients that increase with age and remain positive.

Figure 13.

Examples of different trajectories for changes in regression coefficients with age. (A) Target: Blood urea nitrogen. Predictor: Osmolality. Example of a coefficient that decreases with age but remains positive. (B) Target: Sodium. Predictor: Creatinine. Example of a coefficient that decreases with age, starts positive and ends negative. (C)Target: Sodium. Predictor: Blood urea nitrogen. Example of a coefficient that decreases with age and remains negative. (D) Target: Blood urea nitrogen. Predictor: Sodium. Example of a coefficient that increases with age but remains negative. (E) Target: Albumin, urine. Predictor: Creatinine. Example of a coefficient that increases with age, starts negative and ends positive. (F) Target: Osmolality. Predictor: Blood urea nitrogen. Example of a coefficient that increases with age and remains positive.

We found 70 (2.9%) significant baseline differences between males and females, 16 (0.65%) between non-Hispanic Whites and Hispanics, 29 (1.18%) between non-Hispanic Whites and non-Hispanic Blacks) and 1 (0.04%) false discovery between the two controls. We found even fewer differences for the age trajectories: 10 (0.4%) differences between sexes, 2 (0.08%) differences between non-Hispanic Whites and Hispanics, 3 (0.012%) differences between non-Hispanic Whites and non-Hispanic Blacks and 0 (0.00%) differences the controls (Figure 14).

Figure 14.

Regression coefficients: Differences between demographics. (A) Heatmaps displaying the baseline differences in the regression coefficients between demographic groups. (B) Correlation between the 2,450 regression coefficients in the first group and the 2,450 regression coefficients in the second group. The diagonal black line represents a perfect correlation. The further away from this line the points lie, the bigger the difference between the two groups is, and the lower the correlation coefficient is. (C) Volcano plots reporting the baseline (calculated on the full age range 20-80) differences in regression coefficients between the groups. The horizontal black lines represent the threshold of significance of 0.05 for the Bonferroni corrected p-values. (D) Volcano plots showing which regression coefficients are changing at significantly different rates with age in different demographic groups. The horizontal black lines represent the threshold of significance of 0.05 for the Bonferroni corrected p-values.

Summary

Table 6 summarizes which percentage of correlations, predictabilities and regression coefficients are significantly different from zero at baseline. It also reports which percentage of correlations, predictabilities and regression coefficients significantly change with age. Finally, it recapitulates the percentage of correlations, predictabilities and regression coefficients that are different between sexes and between ethnicities, both at baseline and in term of changes with age.

We also share our results on an interactive Shiny app: http://apps.chiragjpgroup.org/Aging_Biomarkers_Co-Dependencies/.

Dataset

We used the 1999-2014 cohorts of the National Health and Nutrition Examination Survey (NHANES) dataset (Figure 1A), a cross-sectional survey. We first selected the biomarkers under the categories “Laboratory” (blood biomarkers) or “Examination” (anthropometric biomarkers). We acknowledge the different definitions of biomarkers. For the scope of this paper, we refer to “biomarkers” when describing the variables selected, including phenotypes such as height. When a biomarker had its value recorded in both the American and the international metric system, we excluded the latter. For patients who had more than one measure for blood pressure, we took the average. We excluded from the analysis the individuals younger than 20 years old and the individuals aged 80 or above because of the limited sample size in those age ranges. We iteratively excluded the samples and biomarkers with the most missing values from the dataset until we were left with a complete matrix of 27,508 samples and 54 biomarkers, out of the initial 82,440 samples and 1,308 biomarkers (Figure 1B). Out of the initial 1,308 biomarkers, some correspond to repeated measures of a single biomarker, such as blood pressure. Some of those 1,308 biomarkers are measured in a large number of individuals, while some others were only measured on a subset of the cohort, hence the need to extract a complete matrix from this scarce starting matrix.

We calculated the pairwise correlations between those 54 biomarkers and filtered out four other biomarkers, which were more than 90% correlated with another (total cholesterol correlated with cholesterol, hemoglobin correlated with hematocrit, white blood cell count correlated with neutrophils number, mean cell volume correlated with mean cell hemoglobin). We performed logarithmic or exponential transformations on the biomarkers to make their distributions look more Gaussian then we centered and scaled the biomarkers using the survey weights. A description of the demographics can be found in Tables 1 and 2. We tested for a change of the biomarkers with age using a weighted linear regression (the significance threshold was 0.05, after Bonferroni correction.)

We replicated our top finding in term of changes in correlation with age using the UK Biobank dataset. After preprocessing the 502,628 samples for missing values, we obtained a sample size of 470,895. More details about the demographics of this cohort can be found in Tables S1 and S2.

Calculation of the correlation structure of the biomarkers

Calculation of the predictability for each biomarker

To evaluate the predictability for each biomarker, we first split the relevant cohort (in terms of demographics and age range) into a training and a testing set (50/50). We fit an elastic net with alpha=0.5 and we picked the optimal lambda among 100 values suggested by the glmnet R package using a 10 folds cross-validation on the training set. We then generated predictions for the testing set. Next, we switched the training and the testing set and used the same procedure to generate predictions on the other half of the dataset. We then merged the prediction of both halves and compared the predictions to the actual values over all the samples using a non-corrected coefficient of determination (R-squared) to report the predictability. We bootstrapped the calculation of the R-squared 1,000 times to obtain its standard error.

We observed that our models performed significantly better when trained on 600 samples than when trained on 200 samples. Because different age bins have different sample size, and that older age bins tend to have a smaller sample size, we would have obtained biased results that might have a decrease in the predictabilities of the biomarkers with age. We corrected for this bias by using the same number of samples for each age bin when comparing predictabilities for different ages. Similarly, we used the same number of samples to compare between sexes and between ethnicities. Taking the minimum of the sample sizes sometimes left us with too few samples to reach significance so, as a consequence, we used a larger age window for those analyses (Table 7). The windows did not overlap, and each age bin analysis is independent from the other others.

We investigated the age trajectories of predictabilities in two ways. First, as described above, we built a model for each age bin, and used it to estimate the predictability on this same age range. Secondly, we built a single model on the full age range, and then evaluated how this model performed to predict the values of the biomarkers on different age bins. The advantage of the first model using a sliding window is that each model is specific to the age bin and can be analyzed in the context of the correlations and regression coefficients of the same age bin. The downside is the small sample size available to train the models. The advantage of the second model is the larger training sample size. Its downside is that it is not age bin specific and can therefore not be analyzed in the context of the age bin specific correlations and regression coefficients. To be concise, we only present the results obtained using a sliding window in this paper, but the results obtained using a single model can be found on the website.

We examined the demographic differences and the age dependence of predictabilities following the same pipeline as we did when we examined the correlations between biomarkers. We first established a baseline on the full cohort and age range before performing the analysis on the difference age bins and testing for a significant change of the mean. Then we looked at demographic differences both for the baseline and for the age trajectories.

For some analyses, we obtained negative R² values. While it is unusual, this occurs when the performance of a model on the testing set is poorer than predicting the mean target of the testing set for each sample.

Calculation of the regression coefficients for the prediction of each biomarker

To calculate the regression coefficient, we built a model using the same protocol described above to evaluate the predictabilities, but without splitting the data into a training and a testing set. To obtain the variance of the coefficients, we used 100-fold parametric bootstrapping. The pipeline we followed was identical to the one used to explore age and demographics changes for correlations and predictabilities.

We performed the pipeline on “control groups”

In order to have a negative control for our findings with respect to sex and ethnicity differences, we built control groups. We randomly split the samples into two groups and performed the entire pipeline on those two groups like we did on each sex and ethnicity. By comparing the difference between those two groups, we estimated how much of the differences between sexes and ethnicities were due to random sampling, and how much were due to biological differences.

Report of corrected p-values

P-values are probabilities and therefore take values between zero and one, while the negative log p-values have an upper limit of one. However, we chose to report Bonferroni corrected p-values instead of p-values, so each p-value is multiplied by the number of tests we ran and can by bigger than one, and its corresponding negative log p-value is therefore not lower-bounded by zero.

Data availability

Our results can be found online in an interactive form: http://apps.chiragjpgroup.org/Aging_Biomarkers_Co-Dependencies/.