An Information-Theoretic Framework for Identifying Age-Related Genes Using Human Dermal Fibroblast Transcriptome Data

Machine learning algorithms have surpassed classical signal processing methods and have gained dramatic attention of computational biologists due to multiple advantages such as accuracy, robustness and good generalization [12, 27, 30, 35, 37]. Specifically, machine learning has showed its capability for genomic data processing. One good example is the study of mechanisms associated with ageing through gene expression profiling in order to estimate the age [13]. The study of the mechanism of ageing, and in particular, the identification of age-related gene have attracted a lot of attention over the years [14], especially given its implications on the possibility for developing therapeutic approaches to slow down or control the human ageing process [15, 16]. It is hoped that such manipulations can be implemented at different scales, from whole body to organs to small tissues [1]. However, most attention has been directed toward costly experimental studies in laboratories. On the other hand, more recently, computational methods on selecting age-related genes have gained more attention, given their potential for facilitating more cost-effective understanding of mechanism of ageing, for instance, by providing a guide to experimental researchers for more focused investigations. For instance, Srivastava et al [17] employed machine learning tools to predict genes associated with age-related diseases. De Maclhaes et al [18] investigated age-related genes and provided a meta-analysis of work on gene expression data of mice, rats and humans. They considered the noisy nature of gene expression data, and studied how to associate transcriptional changes with the ageing. Their results identified a set of 56 genes as overexpressed over the life span, and another set of 17 as underexpressed. Uddin et. al. [19] performed an analysis age-related variations in the gene expression data of porcine tissues in order to identify some reference genes. They identified 9 genes as the reference genes, and used these for data normalization of the gene expression data.

Avelar et al [20] perform integrative computational analysis to select and validate a subset of age-related genes. Based on their comparative study they claim that human tissues’ gene over-expression with age can be used to characterize ageing and identifying age-related genes as well as age inhibiting genes.

For this work, we use the human dermal fibroblast transcriptome data reported in [13] which contains data on gene expression values. The database of gene expression values consists of a matrix of gene expression values from 27,142 genes across 143 individuals. Each row contains the expression values from the genes across 143 subjects in order of age (from less than 1 to 94 years old). Each column captures the expression values for each of 27,142 different genes for a given individual. We note that, the expression values for any given gene in the data set are measured from different individuals (143 subjects with age ranging from less than 1 to 94 years old), rather than across the lifetime of each individual. For a general study like the one proposed here, this could be important in smoothing out the effects on the analysis that may be due to the impact of one individual.

In the context of unsupervised clustering, we consider each row as a vector representing the expression values for a certain gene across the age, where it can potentially take a binary label with respect to the assumption on its association with ageing. However, for the semi-supervised approach, we have the label for a small number of genes known to be associated with ageing, and we seek to annotate the remaining genes.

We also used a database reported by Tacutu et al [14] which listed 307 genes that are known to be associated with ageing, and 243 genes that are most likely associated with aging.

We compute the entropy of each gene across all 143 subjects and also across the pair of age groups for each age threshold.

The entropy of a random variable was first defined by Claude Shannon in 1948 as the average amount of information or uncertainty in the variable’s possible outcomes [36]. In information theory literature, it is represented as follows:

where ${p_{i}}$ is the probability of the $i$-th outcome. However, we need to compute the entropy of histogram of expression values for each gene as the entropy of the probability distribution of the values within the age bins of the histogram, in which ${p_{i}}$ is defined as the number of values within a given age bin divided by the total number of expression values for that gene [29]. Here we empirically set the number of histogram bins to 7, after testing different values. We select the genes based on their corresponding contribution to aging, as measured by the entropy of expression values for each gene across all the subjects, or within the two separate groups of subjects for a given age threshold, for instance, Group 1 (1 to 25 years old) and Group 2 (26 to 94 years old) using an age threshold of 25. In fact, posing the gene selection problem within the context of age groups allows for a dynamic consideration in the selection process. That is, if gene functionality depends on the age of the subject, then, we expect the model to be able to recognize it. Also, by using multiple age thresholds, it becomes possible to see which genes play more significant roles at different stages in life, as captured by different age ranges.

Entropy Across All Subjects. As the first scenario, we compute the entropy of genes individually across all 143 subjects and then select the top $T$ genes with the highest entropy. The intuition behind this is that the genes conveying more information through their expression values across age range of 1 to 94, are more probable to be related to the ageing process.

Entropy Across Subject Groups. As the second scenario, we divide the subjects into two groups for a given age threshold, based on their age. For instance, using an age threshold of 35 we will group the subjects into two age groups: Group 1: ages [1-35], and Group 2: ages [36-94]. Then the entropy of the genes across the subjects in each age group are computed individually and top $T$ genes are selected based on their entropy ranking. Thus, each age group will produce one gene ranking.

Entropy Difference. Given the binary grouping, for a given gene it becomes possible to consider its entropy difference across the subjects in the two groups. Thus, as part of the feature measurements, we also compute the entropy between the expression values for subjects in Group 1 and Group 2. We rank the genes based on this difference, and select the top $T$ genes for further consideration.

For our experiments, we have set $T$=1000.

Correlation is another feasible metric for exploring the underlying pattern of relationship between two variables – in this case, between subject age and expression values of each gene from the individual. In other words, correlation can be used to measure the possible association between expression values of the genes and the ageing process. The outcome of the correlation between two variables could be positive, neutral, or negative. With positive correlation both variables change in the same direction; with neutral correlation there is no relationship in the change patterns of the variables; and with negative correlation the variables change in opposite directions. Since we are seeking for any type of age-association, we use the absolute value of the correlation for further computation and inference.

There are two formulae for computing the correlation, namely Pearson’s correlation and Spearman’s correlation. The Pearson correlation coefficient can be used to measure the strength of the linear relationship between two variables when both variables follow a Gaussian or Gaussian-like distribution. While Spearman’s correlation can be used for two variables with a nonlinear relationship between them, and the variables may have a non-Gaussian distribution.

Since we are unsure of the distribution and possible relationships between these two variables, gene expression values and age, we need to use Spearman’s correlation as our metric to select the genes with most contribution. The formulation for Spearman’s correlation coefficient is presented as follows:

where ${x}$ and ${y}$ are our variables, ${R(x)}$ and ${R(x)}$ are rank of variable ${x}$ and ${y}$, ${S}$ is the standard deviation and ${C}$ represent the notation for covarience.

Similar to entropy, for each gene we also compute the Spearman’s correlation between the expression values and age, for all subjects, within each of the two age groups, and the difference in the correlation for subjects across the two age groups.

The Kullback-Leibler (KL) divergence (also called relative entropy) measures the difference between two probability distributions [36]. Suppose we have two discrete probability distributions ${P}$ and ${Q}$ which are defined on the same probability space, ${Z}$, the Kullback–Leibler divergence from ${Q}$ to ${P}$ is defined as follows:

For our purposes, given the age threshold, for each gene we compute the KL divergence for the two groups of expression values (as determined by the given age threshold). The genes are then ranked by their KL-divergence over the two groups, and we select the top $T$ genes with maximum divergence.

Given the information-theoretic measurements, we can combine the individual top $T$ rankings, to obtain one integrated set of genes for more refined analysis, for instance, using semi-supervised machine learning methods. For instance, we can combine the tip $T$ genes, by using simple intersection or union for the individual top-$T$ ranked gene sets. An alternative approach is to apply further unsupervised learning techniques on the data, for instance, using $k$-means clustering, but this time using the information-theoretic measurements as the features. The clustering results will then be passed to the semi-supervised learning algorithm for final refinement of the identified genes. In this work, we show results for the two approaches.

As it is presented in Table I, the measurements provide different performance in gene selection, where the best performance is associated with K-L divergence. This could be due to the fact that each gene relationship with aging would only be capture by only one a few of measurements, and not all of them. Considering this, one way to select a subset of genes from 27142 genes and pass them to next stage (semi-supervised approach) is two find the union between top genes selected by each measurement as the mentioned subset. The result of this union of top 1000 genes selected by each of the nine measurements in Table I provides us with 6799 genes. We pass it to semi-supervised gene selection stage where the results are presented in Table II.

A popular clustering method is the $k-means$ algorithm. Though this algorithm is well known, given the significance of clustering in our approach, we briefly describe the algorithm below, for completeness. The $k$-means clustering algorithm, takes a training set ${x^{(1)},...,x^{(m)}}$, and the goal is to group the data points into a few predefined number of coherent clusters. To describe a customised version of the algorithm for this work, say we are given feature vectors for each gene (data point) $x^{(i)}\in\mathbb{R}^{n}$ where $n$ is the dimension of the feature vector. The algorithm initially selects $k$ (here k=2) random data points as the centroid of the $k$ clusters, and then according to the distance of the remaining data points to each of these centroids, assigns a cluster to the data point. Next, the mean of each cluster is computed as the new centroid of the cluster, and the update process of assigning clusters to the data points repeats until all data points remain at the assigned centroids and consequently the centroids could not be updated anymore.

More specifically, here rather than 143 measurements, we represent genes individually by nine information-theoretic measurements as described above, and apply ${k}$-means algorithm on the gene expression dataset to assign each of the 27,142 genes to one of the binary clusters. As explained eariler, we here we set $k=2$ since we wish to classify the genes as being associated with age, or otherwise. Thus, the goal is to predict $k=2$ centroids and a label $c^{(i)}$ for each gene as our data points. The k-means clustering algorithm is as follows:

1.Randomly initialize the cluster centroids ${\mu_{1}}$ and ${\mu_{2}}$.

2. Repeat until convergence ${\{}$

For each ${i\in\{1,2,...,27142\}}$, set:

For every ${j\in\{1,2\}}$, set:

The goal is to perform sub-space search and clustering of the set of 307 known genes by employing $k$-means for ${K=2,3}$, and then comparing the hypothesized age-related genes resulting from the unsupervised learning stage against the members of the respective clusters of the known genes for further similarity analysis. As it is formerly presented, the $k$-means algorithm provides us with $k$ clusters of data points, here clusters of genes. We set number of clusters $k=2$ and $k=3$ and refine the gene selection using the Jensen-Shannon divergence (JSD) and cosine similarity. For $k=2$ as binary clustering, the $k$-means algorithm provided us with two clusters consisting of 176 and 131 known genes, respectively. For $k=3$, we had clusters consisting of 51, 117 and 139 known genes, respectively. Latter the results corresponding to each number of clusters will be assessed.

Age-related genes are not expected to all function in the same way. Some may have a causal relationship with ageing, that is, influencing the aging process, while the function of some others will be influenced by ageing. Even so, some genes with causal relationship with ageing could speed up ageing or slow it down. Moreover, some studies have suggested the possibility of age-related dynamic changes in form of monotonic over-expression, monotonic under-expression, or a mixture of these over the life span. Accordingly, our rationale behind the sub-space search is that these types of different associations, functionalities or monotonic/non-monotonic patterns of expression should manifest in terms of ”different” clusters of age-related genes. Later in this work, we will empirically examine the role of this type of subspace analysis in terms of clustering of the known genes [18].

Jensen-Shannon divergence (JSD) or information radius is an information-theoretic approach to measure the similarity between two probability distributions [36]. This similarity measurement which has been widely employed in bioinformatics and genome comparison, is represented as follows:

where ${A_{i}}$ and ${B_{i}}$ represent the probability distribution for two given vectors of gene expression data, and D (.——.) is K-L divergence as was introduced earlier in Section III-C.

Cosine similarity is a well-known approach for measuring similarity between two vectors in inner product space. With this measure of similarity, we do not need to estimate any probability distribution, we only need to measure the cosine of the angle between the two vectors. Let $\mathbf{A}$ and $\mathbf{B}$ be two non-zero vectors of gene expression values; the cosine similarity between $A$ and $B$ is computed using the dot product and the magnitude as follows:

where ${A_{i}}$ and ${B_{i}}$ represent the elements of the corresponding vectors and ${n}$ is the number of elements of each vector. It is noteworthy to mention that we perform the cosine similarity on both raw vectors of expression data for genes (as presented in Table II and III) as well as a vectors of nine measurements for genes individually. However, we observed that the results for the raw vectors are better. This could be due to the fact that in order to extract the phase information from vectors, we should solely consider the original raw values of the vectors rather than handcrafted features of that vectors.

The unsupervised learning stage will produce two initial sets of selected genes, namely, one with 6492 (6799-307) genes using union of top $T=1000$ ranked genes, and the other one with 5362 (5669-307) genes from $k$-means. We will then refine these initial selections using semi-supervised methods to select those that are most to the set of known genes, using cosine similarity and JSD.

As noted earlier, age-related genes are not necessarily monolithic – they could still cluster into different groups, depending on the basis of the clustering. Thus, we perform the similarity based refinement for our gene selection for three scenarios based on the clustering on the known genes, namely using one cluster, two clusters, and three clusters, respectively,

First Scenario: One Cluster. To perform gene selection for this scenario, we compute the average JSD as well as average cosine similarity of each gene with all 307 genes and then select the top $Q$ genes with highest similarity values. We do this for both JSD and cosine. Accordingly we will have:

where ${|S_{j}|}$ is the absolute value of similarity between gene number $j$ (out of 53622) and the whole set of 307 genes known to be age-related. Similarly, ${JSD_{j}}$ is the average of the JSD between the $j$-th gene and each of the 307 known genes. We use absolute value of cosine similarity as were are interested in any type of association, whether positive or negative similarity. To select the most similar genes in the sense of expression values, we rank the genes based on cosine similarity as well as JSD separately, and select first $Q$ genes for each method. For our experiments, we set $Q$ to 500, which is close to the ${2\%}$ of genes in the human skin presumably with gene expression changes associated with ageing [1, 10].

Second Scenario: Two Clusters.

where ${S_{jc1}}$ and ${S_{jc2}}$ are the similarity between genes number $j$ (out of 5452) and cluster 1 with 176 genes and cluster 2 with 131 genes, respectively. Next, we rank the genes separately, based on ${S_{j}}$, and ${JSD_{j}}$, respectively. Then, for each ranking, we select the top $Q$ genes.

Third Scenario: Three Clusters. Similarly, we perform clustering to divide the 307 genes into 3 clusters, and compute the similarity values:

where ${S_{jc1}}$, ${S_{jc2}}$ and ${S_{jc3}}$are the similarity between genes number $j$ (out of 5452) and clusters 1, 2 and 3, respectively. We select top $Q$ ranked genes separately, based on ${S_{j}}$ and ${JSD_{j}}$ respectively.

In order to evaluate the unsupervised clustering approach, we contrast the results with the list of the 307 genes [14] known to be associated with aging. The results using each of the measurements are presented in Table I. As the table shows, clearly, KL divergence provided the best performance – recognizing 254 genes out of 307 known age-related genes, within its top $T=1000$ ranked genes. The best age threshold appears to be at age 40, followed by age 35. Incidentally, these are quite close to the age range 38 to 40, which Belsky et al [32] found to be significant in considering the pace of ageing in young adults. Rahman and Adjeroh [33] also observed that human physical activity and various biomarkers of biological aging show a turning point around this age range. These show that the age range 38 to 40 could represent a major dividing point for most individuals in terms of age progression.

Besides, the $k$-means algorithm provides us with two clusters of 5669 and 21,473 genes respectively. Although the second cluster contains even genes with all measurements equal to zero, the first cluster (5669 genes) contains all the genes that are known to be age-related. Therefore, we consider this smaller cluster as the one that most probably would contain other age-related genes that are not yet known. Hence, we pass this cluster to the semi-supervised algorithm to leverage knowledge from known genes to further refine the selection of new genes that are most likely to be age-related.

To this end, besides the list of 307 known genes used to train the semi-supervised learning algorithm, we use an additional 243 genes also known to be associated with aging to test our proposed methods for identifying novel genes with a high likelihood of association with aging.

Table II shows the results using the combined top $T$ ranked genes from the individual information-theoretic measurements, using unsupervised learning as input to the semi-supervised method. With $T$ set at 1000, taking a union of all the top $T$ ranked genes (at $\tau_{A}=40$, see Table I) resulted in a total of 6799 unique genes, which were then fed to the semi-supervised learning algorithm. The combined cosine distance and JSD results were obtained by simply computing the average rank from the two methods.

For the semi-supervised scenarios using different number of clusters, the results show that clustering using $k=2$ clusters tend to lead to better results for similarity measurements. Tables II and III show that final selection refinement using inputs as the gene clusters based on clustering in the high-dimensional feature space performed better that using the union of the top $T$ ranked genes from individual information theoretic measurements. The results indicate that, at $k=2$ clusters, the top $Q=500$ identified genes contained most, if not all of the 243 known genes. Thus, each method identified over 250 other genes that are potentially associated with aging. Further analysis showed that, aside from the 243 known genes, taking the intersection of the results from the two approaches (combined results using $k=2$ clusters) resulted in a total of 82 unique genes. Thus, these represent novel age-related genes identified by our proposed approach.