Combining Shapley value and statistics to the analysis of gene expression data in children exposed to air pollution

Model application

We applied CASh to the analysis of gene expression data published in [14] of 23 children from the polluted area of Teplice (TP) and 24 children from the rural control area of Prachatice (PR), in the Czech Republic. We addressed the problem of quantifying the relevance of genes in the TP area using the information provided by the microarray game defined when up-regulated genes are considered (${\overline{v}}^{TP+}$) and the microarray game defined when down-regulated genes are considered (${\overline{v}}^{TP-}$). The relevance of genes was computed as the Shapley values of games ${\overline{v}}^{TP+}$ and game ${\overline{v}}^{TP-}$ (see Methods for more computational details about CASh). The plot of the Shapley value distributions of genes in games ${\overline{v}}^{TP+}$ and ${\overline{v}}^{TP-}$ is shown in Figure 1 (for the formal definition of the notations used in this section see Methods).

Algorithm 1 was applied to compare the Shapley value ϕ of games ${\overline{v}}^{TP+}$ and ${\overline{v}}^{TP-}$ with the Shapley value computed on the microarray game defined when up-regulated genes in PR area are considered (${\overline{v}}^{PR+}$) and the Shapley value computed on the microarray game defined when down-regulated genes in PR area are considered (${\overline{v}}^{PR-}$), respectively. More precisely, Algorithm 1 was applied twice: first, to test the null hypothesis |ϕ _i (${\overline{v}}^{TP+}$) - ϕ _i (${\overline{v}}^{PR+}$)| = 0 against the alternative hypothesis |ϕ _i (${\overline{v}}^{TP+}$) - ϕ _i (${\overline{v}}^{PR+}$)|≠ 0; second, to test the null hypothesis |ϕ _i (${\overline{v}}^{TP-}$) - ϕ _i (${\overline{v}}^{PR-}$)| = 0 against the alternative hypothesis |ϕ _i (${\overline{v}}^{TP-}$) - ϕ(${\overline{v}}^{PR-}$)| ≠ 0.

Selection of significantly modified gene expressions in exposed versus non-exposed children

Note that for each predefined value of p in Table 1, the number of genes selected in the microarray game ${\overline{v}}^{TP+}$ is larger than that selected in the microarray game ${\overline{v}}^{TP-}$. Figure 2 shows the scatterplots of the p-values versus the absolute differences in Shapley values. These plots indicates that the 838 genes selected in ${\overline{v}}^{TP+}$ and the 889 genes selected in ${\overline{v}}^{TP-}$ are closed to the respective Pareto optimal frontiers.

The overlap between CASh and t-test in terms of the number of selected genes when ranked according to their p-values (p <0.05) in both CASh and t-test is represented in Figure 3 (red lines). Figure 3 also shows other three curves: the green line represents the overlap between a list of genes generated with the criterium of highest fold-change of average expression between TP and PR and a list of genes generated with the criterium of highest differential Shapley value between games on TP and PR; the blue line represents the overlap between a list of genes ranked according to their p-value (p < 0.05) in the t-test and a list of genes generated with the criterium of highest differential Shapley value between games on TP and PR; the orange line represents the overlap between a list of genes ranked according to their p-value (p < 0.05) in CASh and a list of genes generated with the criterium of highest differential Shapley value between games on TP and PR. The overlap between CASh p-values and t-test (red line) or between differential Shapley values and t-test (blue lines) is remarkable higher than the overlap between fold-changes and differential Shapley values (green lines). Using CASh instead of differential Shapley value does not determine any significant differences in terms of overlap with the t-test list.

The hierarchical clustering (Figure 4(a)) based on the set of 159 genes with highest Shapley value and unadjusted p-value < 0.01 returned a distinct separation of 22 (cluster B) exposed subjects and 15 (cluster A) non-exposed subjects (accuracy 78.7%). The hierarchical clustering (Figure 4(b)) based on the set of all of the 265 genes differentially expressed at p < 0.01 in the t-test returned a distinct separation of 22 (cluster A) exposed and 21 (cluster B) non-exposed subjects (accuracy 97.7%). The red/green bar on the top of the heat-maps is used to label the subjects of the two clusters provided by K-means clustering with K = 2. K-means clustering shows an accuracy of 88.4% for the list of genes selected by CASh and 97.7% for the list of genes selected by t-test.

Figure 5(a) shows that most of the 160 genes with highest Shapley value difference have also a small p-value: 37 genes out of these 160 have a p-value from CASh bigger than 0.01, and only 2 have a p-value bigger then 0.05. Blue arrows show two genes, precisely A _23_P 166677 (MFSD1, major facilitator superfamily domain containing 1) and A _23_P 106002 (NFKBIA, nuclear factor of kappa light polypeptide gene enhancer in B-cells inhibitor, alpha), with the same Shapley value difference between ${\overline{v}}^{TP+}$ and ${\overline{v}}^{PR+}$ of 0.0063, and very different p-values of 0.029 and 0.004, respectively. As it is shown in the table of Figure 5(a), for gene A _23_P 166677 the difference between the medians of the sample statistics of the Shapley value in ${\overline{v}}^{TP+}$ and ${\overline{v}}^{PR+}$ is null, whereas the analogue difference of medians for gene A _23_P 106002 is 0.00098. The mean of the sample statistic of the Shapley value in a microarray game equals the Shapley value of the game [see Additional file 1]. Differently, the median of the sample statistic of the Shapley value of a microarray game has not an immediate game theoretical interpretation but it is more stable than the mean with respect to exceptional values.

Figure 5(b) shows that the list of 160 genes with the highest Shapley value difference has the same classification success (accuracy 78.7%) of the list of 159 genes selected by CASh at p < 0.01. The same classification accuracy of 78.7% is also achieved by the list of genes obtained from the 160 ones with the highest Shapley value difference after that 47 genes with p-value from CASh greater than 0.01 are removed (Figure 5(c)).

To compare the classification success of the lists of genes selected by CASh and t-test, we performed hierarchical clustering for equal numbers of genes (Figure 6). Table 3 shows a similar high level of classification accuracy for lists of genes by CASh and t-test, with a small increment in accuracy (4.2%) for t-test lists with 33 and 159 genes.

Upon t-test analysis DAVID returned more modified pathways (n = 84) than after CASh selection (n = 50). CASh-based pathway identification shares 11 annotation terms with t-test analysis-based pathway selection.

Simulation study

In microarray studies, the detection of differential gene expression under two different conditions is very important. On the other hand, also the detection of differential gene expression variance may allow to identify experimental variables that affect different biological processes and accuracy of DNA microarray measurements. So, in this simulation we compare the performance of CASh and t-test in selecting genes which differ between two conditions in terms of average expression or expression variance.

To assess the statistical power of CASh as a function of the sample size, we conducted a simulation study. We compared the performance of CASh against t-test on a simulated gene expression data-set of n = 1000 genes obtained by random samples from normal distributions under two simulated conditions: 90% of genes were sampled from the same distribution under the two conditions; the remaining 10% of genes was split in two groups of equal size: one group of genes 2-fold different in average expression between the two conditions and another group characterized by different variability across measures under the two conditions. Both CASh and t-test were applied on the simulated data-set, and genes with p-value smaller than predefined cutoffs, used to control the false positive rate at 0.1 in both methods, were selected. Figure 7 shows that the power of the t-test converges to 0.5 as expected, since half of the genes sampled from different distributions have a fold change not equal to 1. Differently, CASh converges to 1.

Game Theory

In this section we introduce some basic game theoretical notions and definitions. A coalitional game is a pair (N, v), where N denotes a finite set of players and v : 2 ^N → ℝ the characteristic function, with v(∅) = 0. Often we identify a coalitional game (N, v) with the corresponding characteristic function v. A group of players T ⊆ N is called a coalition and v(T) is the worth of the coalition T.

The unanimity game (N, u _R ) based on R ⊆ N is the game described by u _R (T) = 1 if R ⊆ T and u _R (T) = 0, otherwise. Every coalitional game (N, v) can be written as a linear combination of unanimity games in a unique way, i.e. v = ∑_{S⊆N, S≠∅}λ _S (v) u _S (see for instance [33]. The coefficients ${({\lambda }_S(v))}_{S\in 2^N\backslash \{\varnothing \}}$ are called unanimity coefficients or dividends of the game (N, v).

Given a coalitional game (N, v), an allocation or payoff vector is a vector (x _i )_i∈N∈ ℝ ^N assigning to player i ∈ N the amount x _i .

A solution for a class of coalitional games is a function ψ that assigns a payoff vector ψ(v) to every coalitional game in the class. A well-known solution for coalitional games is the Shapley value, introduced by [7].

where π is a permutation of the players, P(π; i) is the set of players that precede player i in the permutation π and n is the cardinality of N.

In [6], the definition of microarray game was introduced as a coalitional game (N, $\overline{v}$) with the objective to stress the relevance ('sufficiency') of groups of genes in relation to a specific condition. Let N = {1,...,n} be a set of genes. On a single microarray experiment on N, a sufficient requirement to realize in a coalition S ⊆ N the association between a condition and an expression property is that all the genes showing that expression property belong to coalition S (sufficiency principle). Different expression properties for genes might be considered like, e.g., under- or over-expression, strong variation, abnormal expression etc. A group of genes S ⊆ N which realizes the association between the expression property and the condition on a single array is called a winning coalition for that array. For example, consider a single microarray experiment on a set of genes N = {1, 2,...,10} under a given condition (e.g., exposure to air pollution) and suppose that only genes 1, 3 and 7 show the expression property (e.g., over-expression). Then, each set of genes S ⊆ N with 1, 3, 7 ∈ S is a winning coalition in such an experiment.

Moving to k ≥ 1 microarray experiments on N, we refer to a Boolean matrix B ∈ {0,1}^{n × k}, where the Boolean values 0 – 1 represent two complementary expression properties, for example the property of normal expression (coded by 0) and the property of abnormal expression (coded by 1). Let B. _j be the j-th column of B. We define the support of B. _j , denoted by sp(B. _j ), as the set sp(B. _j ) = {i ∈ N | B _ij = 1}.

where Θ(T) = {j ∈ K | sp(B. _j ) ⊆ T, sp(B. _j ) ≠ ∅}, with K = {1,...,k} and where card(Θ(T)) is the cardinality of Θ(T). Finally, we define $\overline{v}$(∅) = 0.

Note that the expression of a gene is a continuous variable which hypothetically may assume whatever value across different samples, then it is not at all easy to identify good criteria to discriminate between different expression properties. The binarization method used in this work is described in section Data processing for CASh. For alternative binarization methods in gene expression analysis, see for instance [34, 35].

Comparative Analysis of Shapley Value

Consider two groups of microarray experiments on the same set of genes N = {1,...,n}, respectively collected under two different conditions 1 and 2. Let B¹ ∈ {0, 1}^{n × k}and B¹ ∈ {0, 1}^{n × h}be two Boolean matrices, where B¹ is obtained via a discretization procedure from an expression data set with k biological samples under condition 1, and B² is obtained via a discretization procedure from an expression data set with h biological samples under condition 2.

Let ${\overline{v}}^1$, ${\overline{v}}^2$ be the microarray games corresponding to the Boolean matrices B¹ and B², respectively. Let ϕ(${\overline{v}}^1$) be the Shapley value on the game ${\overline{v}}^1$ and let ϕ(${\overline{v}}^2$) be the Shapley value on the game ${\overline{v}}^2$.

for each i ∈ N, where ϕ _i (${\overline{v}}^1$) is the Shapley value of gene i in the microarray game corresponding to the Boolean matrix B¹ and ϕ _i (${\overline{v}}^2$) is the Shapley value of gene i in the microarray game corresponding to the Boolean matrix B².

We formally present a procedure (Algorithm 1) to test the null hypothesis that each gene in N has the same Shapley value between the two conditions 1 and 2. In fact we want to test the null hypothesis δ _i (ϕ(${\overline{v}}^1$), ϕ(${\overline{v}}^2$)) = 0 against the alternative hypothesis δ _i (ϕ(${\overline{v}}^1$), ϕ(${\overline{v}}^2$)) ≠ 0. More precisely, we introduce a test procedure based on a non-parametric Bootstrap method of re-sampling with replacement (see [12, 13] as general introduction to Bootstrap methods; see [36] as a Bootstrap application to microarray analysis), which is able to test the null hypothesis of no difference between the means of two random samples without assuming under the null hypothesis that the probability distributions in the populations are the same.

Algorithm 1

END.

In Additional files, a more detailed version of the pseudo-code of Algorithm 1 [see Additional file 1] and its implementation [see Additional file 4] are given. A discussion on the generation of raw p-values using bootstrap method and the related procedures to adjust p-values for multiple hypothesis testing is provided [see Additional file 2]. Calculations of CASh on a numerical instance are also illustrated [see Additional file 3].

Description of data processing

We analyzed the microarray gene expression data published in [14]. Study subjects were children from the Teplice (TP) area in the north and from the rural Prachatice (PR) in the south of the Czech Republic, for a total of 47 children; 23 from the TP area and 24 from the PR area. For details on study population, collection and processing of blood, RNA isolation and microarray analysis of gene expression see [14].

Data processing for CASh

The final matrix X of 5873 genes and 47 samples that was distilled from the data filtering and preparation as described above, was split in two distinct expression matrices, X ^TP and X ^PR , whose columns were selected from X accordingly to the 23 subjects from TP area and the 24 subjects from PR area.

First, a procedure aimed to discriminate over-regulated levels of gene expression with respect to expressions measured in the PR area was applied. Each continuous value in the vector X_i.= (X_{i 1},...,X_{i 47}) which was equal to or greater than Mean [$X_{i.}^{PR}$] + Stdev[$X_{i.}^{PR}$] was coded as 1, or as 0 if otherwise. Consequently, a Boolean matrix B⁺ with n rows and 47 columns and with values {0, 1} was generated from X. Separately, a procedure aimed to discriminate under-regulated levels of gene expression with respect to expressions measured in the PR area was applied. Each continuous value in the vector X _i . = (X_{i 1},...,X_{i 47}) which was equal to or smaller than Mean[$X_{i.}^{PR}$] - Stdev[$X_{i.}^{PR}$] was coded as 1, or as 0 if otherwise. Consequently, a Boolean matrix B^- with n rows and 47 columns with values {0, 1} was also generated from X. According to the distinction between PR and TP biological samples, the Boolean matrix B⁺ was split in two different Boolean matrices B^TP+and B^PR+, and the Boolean matrix B^- was split in two other Boolean matrices B^TP-and B^PR-. By relation (2), from the Boolean matrix B^TP+the microarray game ${\overline{v}}^{TP+}$ is defined and, in a similar way, the microarray game ${\overline{v}}^{TP-}$ from the Boolean matrix B^TP-is also defined.

In order to remove those genes whose high level of Shapley value may be attributed to chance, we applied the Bootstrap-based Algorithm 1. In practice, we applied Algorithm 1 (b = 1000) with B^TP+in the role of B¹ and B^PR+in the role of B² and the un-adjusted p-values were computed. In a similar manner, we applied Algorithm 1 (b = 1000) with B^TP-in the role of B¹ and B^PR-in the role of B² and the corresponding un-adjusted p-values were computed.

As further criterium for filtering, genes with Shapley value smaller than the mean plus the standard deviation in both microarray games ${\overline{v}}^{TP+}$ and ${\overline{v}}^{TP-}$ were filtered out. Following this criterium, 838 genes were selected in game ${\overline{v}}^{TP+}$ and 889 genes were selected in game ${\overline{v}}^{TP-}$) (see the highlighted intervals of Figure 1).

for n ≥ 1.

Simulation study

A simulated gene expression data-set of of n = 1000 genes obtained by random samples from normal distributions under two simulated conditions which are denoted by a class variable y ∈ {1, 2}. Nine hundred genes were randomly sampled from a normal distribution with mean= 1 and stdev = 1 under both conditions 1 and 2. The remaining 100 genes were slpit in two sets of target genes (i.e., to be discovered genes), since the parameters of the normal distributions from which the expression values are sampled change with the conditions: 50 genes were randomly sampled from a normal distribution with mean = 2 and stdev = 1 under condition 1, and from a normal distribution with mean = 1 and stdev = 1 under condition 2; remaining 50 genes were randomly sampled from a normal distribution with mean = 1 and stdev = 1 under condition 1, and from a normal distribution with mean = 1 and stdev = 2 under condition 2. To be processed by CASh, each randomly sampled continuous value of gene i, i = 1,..., 1000, under condition 1 (condition 2) which was equal to or greater than the average expression of gene i plus its standard deviation under condition 2 (condition 1) was coded as 1, or as 0 if otherwise.