Multiple data sets

In this study, we consider a detection of gene set enrichment in discordant behaviors (or discordance gene set enrichment) for a series of two-sample genome-wide expression data sets. The term “enrichment in discordant behaviors" will be mathematically defined later. Let K be the number of data sets and let m be the number of common genes among this series of data sets. Each data set is collected for two given groups (same for all K data sets). In general, one group represents a normal status and the other group represents an abnormal status.

For a single two-sample genome-wide expression data set, differential expression analysis and gene set enrichment analysis are usually conducted. The purpose of analysis of differential expression is to identify genes showing significantly up-regulation or down-regulation when two sample groups are compared. The purpose of gene set enrichment analysis is to identify pathways (or gene sets) showing coordinate up-regulation or down-regulation, which may be considered as an extension of differential expression analysis.

Therefore, the following gene behaviors are usually of our research interest in two-sample expression data analysis: positive change (or up-regulation), negative change (or down-regulation) and null (or non-differentially expressed). However, these underlying behaviors are usually not observed and expression data are collected to make statistical inference about them.

Data pre-processing is important for both microarray and RNA-seq data and it has been well discussed in the literature [24–26]. In our study, the data can be downloaded from a well-known public database. We assume that the gene expression profiles have been appropriately pre-processed. In an analysis of multiple expression data sets, it is usually necessary to focus on common genes and gene identifiers can be useful for this purpose. In our study, we divide a relatively large data set into a series of non-overlapping subsets. Therefore, all the genes in the downloaded data are common.

z-Score

Many statistical tests have been proposed for analyzing a two-sample genome-wide expression data set [27, 28]. In this study, the traditional paired-two-sample t-test is chosen for its simplicity (although other statistics could be certainly considered, see below). For each gene in each data set (or subset), we perform the t-test to obtain a t-score. Its p-value is evaluated based on the permutation procedure (randomly switch the tumor/non-tumor labels for each pair of tissues) so that the normal distribution assumption is not assumed for the paired-difference data. All the permuted t-scores are pooled together so that tiny p-values can be calculated [29].

Discordance enrichment

Our proposed method is a type of gene set enrichment analysis. As it has been discussed by Lai et al. [15], we defined “enrichment” as “the number of events of interest is larger than expected” and our “event of interest” in this study is “a list of clearly discordant behaviors” from a gene. If we know whether the expression profile of a gene is up-regulated (simplified as “up”), down-regulated (simplified as “down”) or non-differentially expressed (simplified as “null”) in a data set, then a list of concordant behaviors among K data sets for this gene could be (up, up, …, up), (down, down, …, down) or (null, null, …, null). In this study, we focused on a list with at least one “up” and at least one “down” among K data sets. For example, a list like (down, up, up, …, up) is an event of interest but a list like (null, up, up, …, up) is not. The reason is “down” and “up” can be visually distinguished by the negative (“-") and positive (“+”) signs in z-scores, respectively. However, zero z-scores are rarely observed. Therefore, it is less clear to distinguish “null” from “up” (or “null” from “down”).

In our mixture model, we used normal distributions to model the z-scores. A novel contribution is that the parameter space of our model increases linearly when the number of data sets is increased. This is due to the two-level structure of our model. (The parameter space of a general model for this analysis increases exponentially when the number of data sets is increased). For each gene in each data set, we considered three normal distribution components that represent up-regulation (positive distribution mean), down-regulation (negative distribution mean) and null (zero mean). (Theoretically, p-values under the null hypothesis are uniformly distributed. Therefore, z-scores under the null hypothesis are normally distributed with mean zero and variance one). The mathematical details are described below.

A two-level mixture model

In the above model, \(\phantom {\dot {i}\!}\phi _{\mu, \sigma ^{2}}(\cdot)\) is the probability density function (p.d.f.) of a normal distribution with mean μ and variance σ ². Three components represent up-regulation with μ _1,k >0, down-regulation with μ _2,k <0 and null with μ _0,k =0 (also recall that \(\sigma ^{2}_{0,k}=1\)). For this model, an assumption is that the p.d.f. of z _i,k is simply \(\phi _{\mu _{j_{k},k}, \sigma ^{2}_{j_{k},k}}(z_{i,k})\) if we know the underlying component information j _k for the i-th gene in the k-th data set. However, the component information is usually not observed in practice. Then, we have this one-dimensional mixture model after the introduction of component proportion parameters \(\left \{ \rho _{j_{k},k}, j_{k}=0,1,2 \right \}\) for the k-th data set.

When we extend the above mixture model to a higher dimension (i.e. K data sets), without a structure consideration, the parameter space increases exponentially due to the 3 ^K different component combinations (3 components in each of K data sets). Therefore, when K is not a small number (i.e. K>4), we need a more efficient model [15]. Biologically, when different data sets are collected for the same or similar research purpose, some genes are likely to show consistent behaviors across different data sets and some genes are likely to show different behaviors. For genes likely showing consistent behaviors across K data sets, we consider a complete concordance (CC) multivariate model to approximate the distribution of {z _i,k ,k=1,2,…,K}. For genes likely showing different behaviors across K data sets, we consider a complete independence (CI) multivariate model to approximate the distribution of {z _i,k ,k=1,2,…,K}. (Notice that there is no overlap among multiple data sets. If the component information among these data sets is known, then z-scores are independent.) We first describe the CI model and CC model as below.

This model is simply a product of K one-dimensional three-component mixture-models.

This model has three components and each component is a product of K normal probability density functions.

Notice that this two-level model is still a mixture model. We further assume that \(\{ \mu _{j_{k},k}, \sigma _{j_{k},k}^{2}, j_{k}=0,1,2, k=1,2,\ldots,K \}\) are shared by both CI and CC models. It is evident that the model parameter space increases linearly when the number of data sets (K) increases.

We can use the well-established Expectation-Maximization (E-M) algorithm [19] for parameter estimation. First, it is necessary to introduce some indicator variables (for component information) for the z-scores {z _i,k ,k=1,2,…,K} of the i-th gene. Then, we describe the E-step and M-step.

E-step and M-step are iterated until a numerical convergence is achieved. In this study, the numerical convergence is defined as that the difference between the current log-likelihood and the previous one is within a given tolerance value (e.g. 10⁻⁴).

Enrichment score

As we have discussed in Discordance enrichment, in this study, we focus on genes’ behaviors with at least one up-regulation and at least one down-regulation among K data sets (our event of interest: a gene with clearly discordant behaviors). However, we do not need to enumerate all these combinations (among 3 ^K in total). The related computing can be simplified if we enumerate the compliment events instead. There are three combinations for complete concordance: (up, up,..., up), (down, down,..., down) and (null, null,..., null). They will be excluded. There are \(\sum _{l=1}^{K-1}{K \choose l}\) combinations with both “null" and “up" (without “down") and there are \(\sum _{l=1}^{K-1}{K \choose l}\) combinations with both “null" and “down" (without “up"). They will also be excluded. Then, the remaining combinations are our events of interest (at least one “up" and at least one “down").

False discovery rate

Computational approximation

As discussed in Lai et al. [15], the exact calculation of DES can be difficult due to the complexity of heterogeneous Bernoulli process. A Monte Carlo approximation has been suggested as follows. First, set an integer variable X=0. For the i-th gene in S, simulate a Bernoulli random variable with probability of event ζ _S,i . Then, count the number of events from all genes in S, and increase X by one if this number is larger than m _S θ. Repeat the simulation and counting B times and report X/B as the approximated DES. B=2000 was suggested by Lai et al. [15].

Genome-wide expression data and KEGG pathway collection

Zhang et al. [23] recently conducted a genome-wide expression study for forty-five matched pairs of pancreatic tumor and adjacent non-tumor tissues. The data were collected by the microarray technology (Affymetrix GeneChip Human Gene 1.0 ST arrays) and were made publicly available in the NCBI GEO database [23]. The collections of gene sets or pathways can be downloaded from the Molecular Signature Database [7, 8]. At the time of study, the collections have been updated to version 4.0. In this study, we focus on 186 KEGG pathways for our data analysis. There are 28677 genes available for our discordance enrichment analysis. As we have explained in the Methods, we expect to identify pathways with enrichment in clearly discordant gene behaviors among a series of pre-defined genome-wide expression data sets. (Notice that a pathway with DES∼1 is significantly enriched in clearly discordant behaviors; and a pathway with DES∼0 is evidently not enriched in clearly discordant behaviors).

Data division based on gene PNLIP

The hierarchical clustering tree (with Euclidean distance and the “median” agglomeration method) for the log2-transformed ratio values of gene PNLIP is included in Fig. 1 a. Several major clusters of subjects can be generated if we cut the tree at 0.15. After including these isolated subjects into their nearby clusters, we can obtain seven clusters (subgroups of tumor/non-tumor pairs). Therefore, seven subsets of genome-wide expression data were defined accordingly with sample size 7+7, 7+7, 6+6, 4+4, 6+6, 9+9, or 6+6 (see Fig. 2 a). Figure 3 b shows the paired expression ratio values of gene PNLIP [log2-transformation applied here for the convenience of visualization of up-regulation (positive sign) or down-regulation (negative sign)]. Figure 3 a shows the individual expression values for gene PNLIP in different subsets. Notice that, from Fig. 3 b, subsets 1 represents a clear down-regulation of gene PNLIP, and subsets 6 and 7 represents null and up-regulation of gene PNLIP, respectively.

z-Scores based on gene PNLIP

Figure 4 shows pair-wise scatterplot for comparing z-scores from the seven subsets defined by the paired-ratio of gene PNLIP. Most scatterplots for adjacent or close-to-adjacent subsets are showing a relatively regular positive correlation pattern (implying overall consistent gene behaviors). The scatterplots for far-from-adjacent subsets are mostly showing an irregular weak correlation pattern (implying a considerable amount of inconsistent gene behaviors). As mentioned above, subsets 1, 6 and 7 are representative for down-regulation, null and up-regulation of gene PNLIP, respectively. It is clear that the scatterplot for subsets 7 vs. 1 is showing the most irregular pattern, which implies that many genes have clearly discordant behaviors when gene PNLIP changes its behavior from down-regulation to up-regulation.

Significant pathways based on gene PNLIP

Table 2 lists the significant KEGG pathways identified by the discordance enrichment analysis (with DES>0.80, also the related maximum FDR<0.05). Among these eleven pathways, there are neuroactive ligand receptor interaction, olfactory transduction, alpha-linolenic-acid metabolism and linoleic-acid metabolism pathways. The literature support for the association between pancreatic cancer and each of these pathways will be discussed later. For the olfactory transduction and neuroactive ligand receptor interaction pathways, Fig. 5 shows their z-score pattern changes when all the adjacent subsets are pair-wisely compared and three representative subsets (1, 6, 7, see above for their details) are also pair-wisely compared. For the pairs of subsets 2 vs. 1, 3 vs. 2, concordant behaviors can be overall observed for the genes in these two pathways. Discordant behaviors can be overall observed for the pairs 6 vs. 5, 7 vs. 6, 6 vs. 1 and 7 vs. 1. Particularly for the pair 7 vs. 1 (up-regulation vs. down-regulation for gene PNLIP), the genes in olfactory transduction pathway are mostly down-regulated in subset 1 but evenly up-regulated or down-regulated in subset 7, and the genes in neuroactive ligand receptor interaction pathway are almost evenly up-regulated or down-regulated in both subsets.

KEGG pathway	DES	FDR
Gene PNLIP based analysis
Neuroactive ligand receptor interaction	>0.99	<0.01
Olfactory transduction	>0.99	<0.01
Ribosome	>0.99	<0.01
Maturity onset diabetes of the young	>0.99	<0.01
alpha-Linolenic acid metabolism	>0.99	<0.01
Glycine serine and threonine metabolism	0.98	<0.01
Steroid hormone biosynthesis	0.97	<0.01
Pentose and glucuronate interconversions	0.96	0.01
Ascorbate and aldarate metabolism	0.95	0.02
Linoleic acid metabolism	0.84	0.03
Proximal tubule bicarbonate reclamation	0.81	0.04
Gene TP53 based analysis
Neuroactive ligand receptor interaction	>0.99	<0.01
Olfactory transduction	>0.99	<0.01
Cardiac muscle contraction	0.89	0.04
alpha-Linolenic acid metabolism	0.85	0.07
Linoleic acid metabolism	0.80	0.09

KEGG pathways with DES >0.80 are listed by decreasing order of DES with related FDR. KEGG pathways in bold font are identified by both analysis (one based on gene PNLIP and the other based on gene TP53)

Data division based on gene TP53

The hierarchical clustering tree (with Euclidean distance and the “median" agglomeration method) for the log2-transformed ratios values of gene TP53 is included in Fig. 1 b. Several major clusters of subjects can be generated if we cut the tree at 0.03. After including these isolated subjects into their nearest clusters, we can obtain six clusters (subgroups of tumor/non-tumor pairs). Therefore, six subsets of genome-wide expression data were defined accordingly with sample size 4+4, 7+7, 6+6, 13+13, 10+10, or 5+5 (see Fig. 2 a). Figure 6 b shows the paired expression ratio values of gene TP53 [log2-transformation applied here for the convenience of visualization of up-regulation (positive sign) or down-regulation (negative sign)]. Figure 6 a shows the individual expression values for gene TP53 in different subsets. Notice that, from Fig. 6 b, subsets 1 represents a clear down-regulation of gene TP53, and subsets 3 and 6 represents null and up-regulation of gene TP53, respectively.

z-Scores based on gene TP53

Figure 7 shows pair-wise scatterplot for comparing z-scores from the six subsets defined by the paired-ratio of gene TP53. Many scatterplots for adjacent or close-to-adjacent subsets are still showing a relatively regular positive correlation pattern (implying overall consistent gene behaviors). Almost all the scatterplots for far-from-adjacent subsets are showing an irregular weak correlation pattern (implying a considerable amount of inconsistent gene behaviors). As mentioned above, subsets 1, 3 and 6 are representative for down-regulation, null and up-regulation of gene TP53, respectively. All the pair-wise scatterplots for these three subsets are showing irregular patterns (with the scatterplot for subsets 6 vs. 1 the most irregular), which implies that many genes have clearly discordant behaviors when gene TP53 change its behavior from down-regulation to null, and then to up-regulation.

Significant pathways based on gene TP53

Table 2 list the significant KEGG pathways identified by the discordance enrichment analysis (with DES>0.80, also the related maximum FDR<0.10). Among these five pathways, there are neuroactive ligand receptor interaction, olfactory transduction, alpha-linolenic-acid metabolism and linoleic-acid metabolism pathways (which have been identified above by the analysis based on gene PNLIP). For the olfactory transduction and neuractive ligand receptor interaction pathways, Fig. 8 shows their z-score pattern changes when all the adjacent subsets are pair-wisely compared and three representative subsets (1, 3, 6, see above for their details) are also pair-wisely compared. For the pairs of subsets 6 vs. 5, 5 vs. 4 and 4 vs. 3, concordant behaviors can be overall observed for the genes in these two pathways. Discordant behaviors can be overall observed for the pairs 2 vs. 1, 3 vs. 2, 3 vs. 1, 6 vs. 1 and 6 vs. 3. Particularly for the pair 6 vs. 1 (up-regulation vs. down-regulation for gene TP53), the genes in olfactory transduction pathway are mostly down-regulated in subset 6 but evenly up-regulated or down-regulated in subset 1, and the genes in neuractive ligand receptor interaction pathways are somewhat evenly up-regulated or down-regulated in both subsets.

Literature support

We have conducted a discordance enrichment analysis based on gene PNLIP and a discordance enrichment analysis based on gene TP53. Among two lists of identified pathways, there are four in common: neuroactive ligand receptor interaction, olfactory transduction, alpha-linolenic-acid metabolism and linoleic-acid metabolism pathways (see Table 2). To further understand these pathways, we have checked the related biomedical literature.

The genome-wide expression data analyzed in this study were collected based on the microarray technology for RNA profiling. Genome-wide association study (GWAS) data have also been collected for pancreatic cancer research based on the microarray technology for DNA profiling (single nucleotide polymorphism, or SNP). Wei et al. [32] recently conducted a pathway analysis for a large GWAS data on pancreatic cancer research. They reported only two pathways. Interestingly, these two pathways are neuroactive ligand receptor interaction and olfactory transduction pathways (top two identified from both of our analysis results, see above for details). Notice that their findings were based on a different type of molecular data. This is a strong support for the discordance enrichment analysis results.

We also found at least one support for both alpha-linolenic-acid metabolism and linoleic-acid metabolism pathways. Wenger et al. [33] conducted a study on the roles of alpha-linolenic acid (ALA) and linoleic acid (LA) on pancreatic cancer and they observed an association between the disease and these two fatty acids.

Insignificant pathways

Figure 9 shows 186 DES based on PNLIP vs. DES based on TP53. These two lists of DES’s are highly correlated (Spearman’s rand correlation 0.642), although some pathways identified in the analysis results based on PNLIP are not significant in the analysis results based on TP53. (Notice that a pathway with DES∼1 is significantly enriched in clearly discordant behaviors; and a pathway with DES∼0 is evidently not enriched in clearly discordant behaviors.) Only a small number of pathways were identified by the discordance enrichment analysis. The histograms in the figure show that most pathways are showing insignificant DES’s. For each of two analysis results, there are more than 140 pathways (among 186) with DES<0.05. The number of pathways with both DES<0.01 or both DES<0.05 is 111 (60%) or 138 (74%), respectively. For both DES<0.25, 0.5 or 0.75, there are 154 (83%), 164 (88%) or 173 (93%) pathways, respectively. Therefore, most pathways are evidently not enriched in clearly discordant behaviors among the series of subsets defined by the paired expression ratio of gene PNLIP; neither are they among the series of subsets defined by the paired expression ratio of gene TP53. Many disease related pathways have been listed by KEGG (http://www.genome.jp/kegg/pathway.html). The collection of pancreatic cancer related pathways (or KEGG pancreatic cancer) and the collection of cancer related pathways (or KEGG pathways in cancer) are not enriched from both analysis results (DES<0.001). Among the pathway components of these two collections (e.g. cell cycle pathway, apoptosis pathway, etc.), the highest DES value is <0.01 for the PPAR signaling pathway from the analysis results based on PNLIP, and the highest DES value is <0.05 for the cytokine-cytokine receptor interaction pathway from the analysis results based on TP53. Pathways like hedgehog signaling, proteasome, and primary immunodeficiency are also showing low DES values (all <0.05).

Expression profiles of PNLIP vs. TP53

PNLIP is a gene shown recently its association with pancreatic cancer [22]. TP53 is a well-known tumor suppressor gene. From the above comparison, it is interesting that the discordance enrichment analysis results based on PNLIP are highly correlated with the discordance enrichment analysis results based on TP53. To further understand this correlation, we compared the expression profile of PNLIP with the expression profile of TP53. Figure 2 a shows a relatively weak negative correlation (Spearman’s rank correlation -0.250) between two lists of paired-ratios but the correlation is not statistically significant (p-value=0.098). In the non-tumor group (Fig. 2 b), the negative correlation (Spearman’s rank correlation -0.318) achieves a p-vlaue 0.033. In the tumor group (Fig. 2 c), the negative correlation (Spearman’s rank correlation -0.276) is again not statistically significant (p-value=0.066). Furthermore, the ratio cutoff values for defining subsets were added to Fig. 2 a. A contingency table can be generated according to these grids (for example, the cell number is one for row one and column one in the table). The chi-square test for this sparse contingency table is not statistically significant (simulation based p-value >0.3). Therefore, in summary, gene PNLIP may be negatively associated with gene TP53 but no clear statistical significance has been observed in this study.

Comparison to gene set analysis

Efron and Tibshirani [34] have proposed a gene set analysis (GSA) method for analyzing enrichment in pathways (or gene sets). It was suggested by Maciejewski [35] that this method is preferred in a gene set enrichment analysis. In some situations of integrative data analysis, different data sets cannot be simply pooled together. For each data set, the p-value of enrichment in up-regulation can be obtained for each gene set. To integrate the p-values from multiple data sets (for the same gene set), we can consider Fisher’s method (Fisher’s combined probability test). log-Transformed p-values are summed up and then multiplied by -2, which is well-known to follow a chi-squared distribution under the null hypotheses. In this way, we can perform an integrative gene set enrichment analysis of multiple data sets (when different data sets cannot be pooled together). Gene sets (or pathways) can be ranked by their chi-squared p-values. (Similarly, the p-value of enrichment in down-regulation can also be obtained by GSA for each gene set and each data set. Then, the related chi-squared p-values can be calculated by Fisher’s method.) Notice that, our analysis purpose is to detect discordance enrichment among multiple data sets. However, the discordance feature is usually not considered in a traditional integrative analysis.

In this study, our analysis results were based on several subsets divided from a genome-wide expression data set with a relatively large sample size. These subsets could be pooled back (to be the original large data set). Therefore, we applied GSA to the original data (so that we could take the advantage of its relatively large sample size). However, after considering the adjustment for multiple hypothesis testing, no pathways (or gene sets) could be identified even at the false discovery rate 0.3 (or FDR<30%). (Therefore, the detail of GSA results is not reported).