Data-driven detection of age-related arbitrary monotonic changes in single-cell gene expression distributions

Proposed workflow

We developed a novel data analysis workflow to perform an association analysis of the distribution of single-cell gene expression profiles and age (Fig. 1). Our workflow introduces a method that can detect arbitrary monotonic changes in single-cell expression profiles. This is achieved by combining the quantification of distance within probability distributions using information theory, along with an association analysis based on distance matrices. The input data are single-cell RNA-seq data and donor age information from different donors. The single-cell expression data of a gene can be expressed as a probability distribution. By using the distance measure between probability distributions, which has been studied extensively in statistics and information theory, the dissimilarity of the single-cell expression profile between samples can be calculated (Okada et al., 2022).

In this workflow, The following procedure was used to quantify the distance between two probability distributions. First, the maximum and minimum values of the pooled data for each of the two samples were calculated. Between them, 1,000 equally spaced grids were set up, and the probability density of each grid was estimated with the density function of R. The probability density values of all grids were then normalized so that the sum was 1, and the single-cell expression profile was represented by a discrete probability mass function. Hellinger distances were calculated for these probability mass functions. Helligner Distance, a distance metric between two distributions P and Q could be expressed as: (1)${\displaystyle H\left(P,Q\right)=\frac{1}{\sqrt{2}}\times \sqrt{\sum _{m=1}^k{\left(\sqrt{P_m}-\sqrt{Q_m}\right)}^2}}$

The Hellinger distance is one of the typical distance metrics between two probability distributions and is also used in computational biology methods for cytometry or NGS data (Gingold et al., 2015; Cheng et al., 2023). The advantage of this method is that it can handle arbitrarily shaped probability distributions because it estimates probability distributions and calculates distances in a non-parametric manner.

In addition, an age distance matrix based on the absolute value of the age difference between the two samples was calculated. We used the association analysis between the two distance matrices proposed in a method previously proposed in the field of ecological data science (Somerfield, Clarke & Gorley, 2021). The Spearman coefficient is calculated between the elements of the upper triangular matrix of the distance matrix, which serves as an indicator of this association. Spearman’s correlation coefficient is calculated for each gene to comprehensively examine genes for which changes in the distribution of single-cell expression profiles are associated with aging. The statistical significance of the association was assessed by the Mantel test (Mantel, 1967) between the two distance matrices. In this study, the Mantel test was performed with the default settings for R packages “vegan”(v2.6.4), and the number of permutations is set to 1,000.

Simulation data generation

A simulation data analysis was performed to verify the accuracy of the proposed computational workflow. We established four ordered age categories (Age = 1, 2, 3, 4) and corresponding normal distribution, where the dispersion parameter σ follows three types of the age-related changing pattern (monotonically increasing with age, stabilizing with age, non-monotonic change with age). We considered a total of 20,000 genes, 5% of genes follow the monotonic increasing mode, 5% of genes follow the non-monotonic mode, whereas 90% of genes follow the stable mode. Under the true distribution defined for each age, ten distributions with slightly different means and standard deviations were generated for each age category. Under these true distributions, we generated ten distributions with slightly different mean and standard deviation for each age category.

For each gene, the normal distribution for the i-th sample (Normal(mean = μ_i, sd = σ_i)) was generated as follows. The mean (μ_i) and standard deviation (σ_i) of the normal distribution for the i-th sample are generated according to the following distribution. ${\displaystyle {\mu }_i\sim Normal\left(mean={\mu }_{true},sd=0.1\right)}$ where μ_true is the true mean value and set 10 in this simulation for all ages. ${\displaystyle {\sigma }_i\sim Normal\left(mean=\sigma \left[Age=ag{e}_i\right],sd=0.1\right)}$ where σ[Age = age_i] is the true sd assigned to the donor age of the ith sample and age_i is the i-th donor’s age. The σ[Age = age_i] is assigned different values depending on age. In the gene with monotonic increase pattern, σ[Age = 1], σ[Age = 2], σ[Age = 3], σ[Age = 4] are 0.5, 1, 1.5, 2, respectively. In the gene with a non-monotonic pattern, σ[Age = 1], σ[Age = 2], σ[Age = 3], σ[Age = 4] are 0.5, 1, 1, 0.5, respectively. In the genes with a stable pattern, σ[Age = 1], σ[Age = 2], σ[Age = 3], σ[Age = 4] are all 0.5. This procedure created a set of normal distributions with slightly different parameter values under each age-related parameter changing pattern. From the distribution of each gene of each sample, 1,000 cells were sampled to create statistical samples. We applied the proposed workflow to the simulated dataset to evaluate the performance of the methodology.

Real scRNA-seq data analysis

We then applied our method to the large-scale single-cell RNA-seq dataset for the aging study. We used a public dataset of aging mouse single-cell RNA-seq data (The-Tabula-Muris-Consortium, 2020). We downloaded the preprocessed data file (tabula-muris-senis-droplet-processed-official-annotations.h5ad) from the Tabula Muris Senis website and used it for downstream analysis. This dataset contains single-cell gene expression data for 20,138 genes in 245,389 cells derived from various organs of multiple mouse donors of different ages. The donor mice are 23 individuals with ages of 1 month, 3 months, 6 months, 18 months, 24 months, or 30 months. Of all cells, 164,027 are from male mice and 81,362 are from female mice. Since sex differences are known to occur in age-related changes in gene expression (Boheler et al., 2003; Hägg & Jylhävä, 2021), we only used a male dataset for the downstream analysis. The top 3,000 genes with the highest average expression in the cells for each age category were included in the analysis. We focused our analysis on four organs (kidney, limb muscle, lung, and marrow) that were measured in more than 10 mice in this dataset. We applied the proposed workflow to this preprocessed dataset and calculated the score of the age-related change of the single-cell expression profile for each gene.

Simulation data analysis

A simulation data analysis was performed to verify the accuracy of the proposed computational workflow. We established four ordered age categories (Age = 1, 2, 3, 4) and corresponding normal distribution, where the dispersion parameter σ follows three types of age-related changing pattern (monotonically increasing with age, stable to age, non-monotonic change with age). We considered a total of 20,000 genes, 5% of genes follow the monotonic increasing pattern, 5% of genes follow the non-monotonic pattern, whereas 90% of genes follow the stable pattern. The true distribution plots of each pattern are shown in Fig. 2. Under the true distribution defined for each age, 10 distributions with slightly different means and standard deviations were generated for each age category (example is see Fig. 3). From each distribution, 1,000 cells were sampled to create statistical samples.

In genes with monotonic parameter changes, consistently the highest correlation coefficients were observed (Fig. 4A). In the genes with non-monotonic changes, differences between groups were detected as moderate correlation coefficients (Fig. 4A). In the stable situation, the correlation coefficients were distributed around zero (Fig. 4A). While small Mantel test P-values are observed in the genes with monotonic increasing or non-monotonic parameter change, it is not in the gene with stable pattern (Fig. 4B). This result is consistent with the interpretation of the correlation coefficient results. Pseudo-bulk data analysis using the same dataset has shown that correlation coefficients near zero are observed in all cases (File S1). Pseudo-bulk expression values were defined as the average of the single-cell expression values in each sample. Conventional bulk gene expression analysis cannot directly capture such changes in variance because it is based on obtaining the average of the expression values of the cells in the sample. These results suggest that this method can detect the association between changes in the distribution of single-cell expression profiles and aging in a data-driven manner.

We performed the following computer experiments to verify the validity of this workflow. First, we performed an analysis using the Jensen–Shannon (JS) distance and the Kolmogorov–Smirnov (KS) distance instead of the Hellinger distance as a measure of the distance between probability distributions. JS distance is a symmetric distance measure based on KL divergence (Lin, 1991). In the calculation of KS distance, we applied the KS.diss functions in an R package “provenance” (Vermeesch, Resentini & Garzanti, 2016) to the two statistical samples. For both distance indices, we observed the same results as when using the Hellinger distance (File S2). The results have shown that the distance between probability distributions is interchangeable with other indices proposed in the field of information theory.

Real scRNA-seq data analysis

We applied the proposed workflow to the large-scale single-cell RNA-seq dataset for the aging study and presented the histograms of the correlation coefficients obtained for the four organs in Fig. 5. The shapes of the histograms were similar among the four organs, and genes with a correlation coefficient >0.85 cutoff value were considered as those showing changes in single-cell expression during aging (Table 1). In all tissues, the presence of genes with large positive correlation coefficients is observed. In limb muscle, 13 genes passed this threshold while no genes passed this criterion in other tissues.

For genes detected by correlation coefficients, visualization of the shape of the distribution in a histogram or QQ plot would be useful for the biological interpretation of the results. Histograms and QQ plot visualization for all the genes listed in Table 1 are shown in Files S3 and S4, respectively. Even among those genes that are strongly associated with age and single-cell expression profile, we observe that the specific patterns of change are diverse.

The gene with the largest correlation coefficient was Jhdm1d. Jhdm1d, also known as KDM7A gene, is primarily known for its involvement in epigenetic regulation through histone demethylation, impacting gene expression and cellular processes which indirectly related to the aging process (Yang et al., 2019). Figures 6 and 7 are histograms and QQ plots of the single-cell expression profile of the Jhdm1d gene as an example of other distributional aging changes. The distributions are also observed as a mixed distribution of cells negative and positive for this gene. At 30 months of age, the peak of the distribution of Jhdm1d-positive cells are observed to shift to the right.

As another example, the histograms and QQ plot of the single cell expression distribution of Aldh2 gene were present in Figs. 8 and 9, respectively. It is known that Aldh2 deficiency promotes age-related muscle loss, especially in oxidative fibers, which may be associated with an increased accumulation of oxidative stress via mitochondrial dysfunction (Kasai et al., 2022). The distributions are also observed as a mixed distribution of cells negative and positive for this gene. The number of cells in the bin with range between 0.5 to 1 decreases with age, suggesting that there may be a decrease in the number of intermediate cells between positive and negative cells. Screening using correlation coefficients, combined with visualization using histograms and QQ plots, allows us to examine changes in diverse single-cell expression profiles with aging.

Normally, changes in single-cell expression profiles are affected by both shifts in cell-level expression and changes in fractions of cellular subsets, but this method can detect them without distinguishing between them. If you want to focus only on shifts in cell-level expression, you can apply this method to a single cell type after cell type annotation with existing methods. We applied our workflow to each tissue and the major cell type. In this analysis, we used the following combinations of tissue and cell types as examples (Kidney: “B cell”, limb muscle: “mesenchymal stem cell”, Lung: “classical monocyte”, Marrow: “granulocyte”) as an example. Histograms of the correlation coefficients for each cell subset analysis are shown in File S5 . This analysis also suggested the presence of many genes associated with age. As with the whole tissue analysis, the age-related genes passing the threshold (0.85) have been observed only in limb muscle (21 genes, File S6). These genes were not consistent with the age-related genes in the analysis of whole tissue analysis. There are only two genes (Grb10, plac9) overlapped between two results.These results suggest that the detected genes will be different depending on the focus and range of cell populations.