Approaches to identify significant pathways from high-throughput quantitative data have been developed in recent years. Still, the analysis of proteomic data stays difficult because of limited sample size. This limitation also leads to the practice of using a competitive null as common approach; which fundamentally implies genes or proteins as independent units. The independent assumption ignores the associations among biomolecules with similar functions or cellular localization, as well as the interactions among them manifested as changes in expression ratios. Consequently, these methods often underestimate the associations among biomolecules and cause false positives in practice. Some studies incorporate the sample covariance matrix into the calculation to address this issue. However, sample covariance may not be a precise estimation if the sample size is very limited, which is usually the case for the data produced by mass spectrometry. In this study, we introduce a multivariate test under a self-contained null to perform pathway analysis for quantitative proteomic data. The covariance matrix used in the test statistic is constructed by the confidence scores retrieved from the STRING database or the HitPredict database. We also design an integrating procedure to retain pathways of sufficient evidence as a pathway group. The performance of the proposed T2-statistic is demonstrated using five published experimental datasets: the T-cell activation, the cAMP/PKA signaling, the myoblast differentiation, and the effect of dasatinib on the BCR-ABL pathway are proteomic datasets produced by mass spectrometry; and the protective effect of myocilin via the MAPK signaling pathway is a gene expression dataset of limited sample size. Compared with other popular statistics, the proposed T2-statistic yields more accurate descriptions in agreement with the discussion of the original publication. We implemented the T2-statistic into an R package T2GA, which is available at https: //github. com/roqe/T2GA.
Great progress has been made toward the development of high-throughput technologies and their application to biological and clinical research. In a quantitative experiment, genes or proteins with significant changes in expression are potential to have important roles in a given phenotype or phenomenon. Therefore, the analysis of quantitative experimental data generally produces a list of differentially expressed genes or proteins in order. The list may share some insights if the aim of the experiments is restricted to few targets. As regards high-throughput data, the list hardly provides biological understanding of the mechanisms being studied, since the data involve complicated regulations among biomolecules and the number of biomolecules is too large to examine all candidates individually. Systematically investigating the underlying mechanisms from a high-throughput data therefore become a new challenge. To confront the challenge, one idea is to apply pathway analysis to identify the genes or proteins that are known to be involved in a biological process or interaction based on the existing knowledge. Many approaches of pathway analysis have been developed over years concerning different methodologies. Some general reviews of pathway analysis approaches can be found in [1–4]. These approaches can be broadly divided into two major factions: a competitive null with a gene sampling model and a self-contained null with a subject sampling model. The null hypothesis of a competitive test suggests that the target pathway (or a predefined gene set) is differentially expressed as well as the rest of all pathways. Practically a competitive null is closely tied to (although not necessarily) a gene sampling approach [5]. A gene sampling model principally implies the independence assumption; the assumption presupposes that genes are expressed independently of each other so that these genes can be manipulated as the sampling subject to produce the null distribution. The null hypothesis of a self-contained test, on the other hand, suggests that the target pathway is not differentially expressed between distinct phenotypes. Using the phenotypes of the experiments as the sampling subject is the setup of a subject sampling model. In other words, a competitive null with a gene sampling approach let pathways compete with each other in order to rank these pathways and use the number of genes as the sample size in the meanwhile; whereas a self-contained null with a subject sampling approach examine each single pathway to determine if the pathway is indeed differentially expressed between phenotypes and use the number of experiments as the sample size. Intuitively, a competitive null aims to find the pathway that is most significant among all pathways; a self-contained null aims to find the pathway that is the most significantly expressed between phenotypes. The classification of null hypotheses and sampling methods is firstly suggested by Geoman et al [5]. Even both categories have their own pitfalls and benefits, the authors suggested using a self-contained null with a subject sampling model in pathway analysis. Competitive null with gene sampling model usually implies the independence assumption which may produce inaccurate small p-values, cause serious false positives [6], and result misleading interpretations [5]. Further methodology issues of using a self-contained null can be found in [7–9], and of using a competitive null can be found in [10]. Another issue of pathway analysis comes from the construction of test statistics. These statistics can be divided into univariate tests and multivariate tests. Univariate statistics, such as a modification or a weighted summation of the t-scores [11–14], only focus on expression of genes or proteins and assume these biomolecules as independent units for statistical ease. The independence assumption ignores the associations among biomolecules with similar functions or cellular localization, as well as the interactions among them manifested as changes in expression ratios. In contrast, multivariate statistics take into consideration the associations among genes or proteins [15–20]. Some methodology studies [9,21–23] have evaluated univariate tests and multivariate tests with synthetic and experimental datasets. Compared with multivariate tests, univariate tests generally result a decrease in statistical power [21] and an increase in false positive rate [22] along with the rise of average correlation. Multivariate approaches usually incorporate the sample covariance matrix into calculation to address biological interaction. However, the sample covariance may not be precise enough to estimate the associations if the sample size is very limited. Compared with other gene expression data, proteomic data produced by mass spectrometry are more difficult to analyze systematically due to the limited number of experiments. This limitation causes current multivariate tests incompetent because the sample covariance will not be a robust statistic. The composition of pathway diagrams also become a challenge to pathway analysis. A pathway is a group of biomolecules that participate in a particular cellular process. The members of a pathway are usually defined by the tradition (i. e. , the history of pathway discovery) of molecular biology scientists. The structure of pathway diagrams is not standardized and therefore arises some issues to pathway analysis. First, the same pathway from different databases or other sources may have the same core members but different side members. Under different experimental conditions, the size of accessible biomolecules also changes. For example, the phosphorylation proteomic data may not provide information to the proteins not belong to phosphoproteome. Since the number of members within a pathway is not a constant, using this number as a parameter to determine if the pathway is significant or not may lead to inconsistent results. This issue arises with the assumption of a competitive null. Second, some molecules appear over pathways may play important roles as communication centers. For example, p53 appears in 38 pathways in the KEGG database. The shared members may interact or cooperate with each other and form a functional module. If this module is regulated (i. e. , the members within the module are differentially expressed as they cooperate together), the subsets of this module may present in abundant pathways and make these pathways seemingly significant. To identify the regulated pathway among the significant pathways of common modules, biologists usually utilize the distinctive molecules participating in that specific pathway; whereas pathways do not contain distinctive molecules may be irrelevant to the underlying mechanism. However, most of the current approaches do not take consideration of this issue. The suggestion of irrelevant pathways due to the redundancy over the significant pathways usually causes confusion to data description. A recent study [24] also focuses on the second issue. They demonstrate how the shared members affect p-values, and try to address this problem under a competitive null. In this study, we introduced a multivariate test, based on the Hotelling’s T2-statistic, to perform pathway analysis for quantitative proteomic data. The most serious problem of analyzing proteomic data produced by mass spectrometry is the limited number of experiments. We usually obtain only a few replicates (biological or technical) per experimental condition. To manage this issue, we had two special designs in our test. First, instead of using the sample covariance matrix (which is not robust when the sample size is limited), we use the covariance matrix that is constructed of the probabilistic confidence scores provided by the STRING and HitPredict databases. The proposed T2-statistic is then built of the protein expression profile and the covariance matrix to consider the expression level of individual proteins along with the associations among them. Second, we designed a self-contained model to produce a null distribution of altering protein expression while retaining the structure of protein associations. We are not capable of applying a subject sampling model because the number of experiments is too limited. In addition to the settlement of sample size issue, we designed an integrating procedure to categorize significant pathways as well as to avoid redundancy. The performance of the proposed T2-statistic is demonstrated using five public experimental datasets with different levels of biological complexity: the T-cell activation, the cAMP/PKA signaling, the effect of dasatinib on the BCR-ABL pathway, the differentiation process of myoblast, and the protective effect of myocilin via the MAPK signaling pathway. The first four datasets are proteomic data produced by mass spectrometry; the last dataset is a gene expression data of low sample size. We compared T2 with other popular statistics: DPA [25], GSEA [26], DAVID [27,28], and IPA [29]. For most of the situations, T2 yielded more accurate descriptions in agreement with the discussion of the original publication.
We took four proteomic datasets of different biological complexity and experimental properties to demonstrate our approach. To be comprehensive, the testing datasets include the case of pathway activation and inhibition; also the case of signaling phosphoproteome and cellular proteome. We only used the final ratios provided by the datasets since they may have different integration approaches (e. g. to combine the results of biological and/or technical repeats, to handle missing data or outliers) under different experimental designs. The summarized ratio is also the most available format for quantitative proteomic data. In this situation, the sample size of data becomes only one, calculating a covariance matrix is not even possible. Our approach provide a solution to undertake this difficulty. We also applied our approach on a gene expression dataset of three samples to demonstrate that the general idea is applicable to other quantitative data of low sample size. To provide more generalized results, we choose two pathway databases and two protein-protein interaction databases: KEGG and Reactome provide pathway categories served as predefined gene sets, STRING and HitPredict contributes the confidence scores to estimate the covariance between protein expressions. As we mentioned in Introduction, subsets of one common active module may cause a lot of pathways statistically significant. These pathways may only have slight relevance to the target mechanisms (Fig 1). To avoid misinterpretation due to irrelevant pathways being reported, we identify delegates to categorize pathways into pathway groups. A pathway group is a set of pathways, in which the pathway being the superset of other pathways is defined as the delegate. Since other pathways do not show any distinctive proteins to support themselves, the significance of other pathways may simply originate from the regulation of the delegate. In other words, we want to avoid the situation that a pathway is enriched only because of some common proteins that are shared with other pathways. Our pathway integration procedure operates as follows. The set of the pathways being evaluated is denoted by P and each pathway P ∈ P is associated with a p-value. We iteratively perform the following steps until all the pathways in P are assigned to a specific pathway group M. A demonstrative example can be found in S1 Fig. The pathway integration procedure aims to find the pathways with the most sufficient information to represent current data. Please notice that this idea is not similar to pathway hierarchy in which higher level pathways are defined as the supersets over lower level pathways, in that case only the pathways of highest level are able to be delegates. One can only apply the integration procedure to pathways without hierarchy relationship. In addition, the procedure is only applicable to the statistical test using a self-contained null since the measure of significance is independent for each pathway.
We discussed the flow of signal transduction in time order since the dataset is a time-series of 5 min, 15 min, and 60 min experiments. In the beginning, the treatment anti-CD3ϵ activated TCR signaling pathway. The TCR signaling pathway in KEGG and IPA depicts the signals from the TCR receptors all the way to the IL-2 expression. Usually, the response of signal transduction comes rapid, we expected that the TCR signaling pathway should be enriched in early time points; the downstream from the TCR signaling pathway to the IL-2 expression should be enriched in late time points. From Table 2 we found that T2×ST, T2×HP, DAVID, and IPA enriched the TCR signaling pathway in all experiments; GSEA enriched the TCR signaling pathway in 5 and 15 min experiments only. The ranks of the TCR signaling pathway are top in T2, DAVID, and IPA. The downstream of the TCR signaling pathway was illustrated in Fig 3, there are three possible routes directed from TCR signaling. The original publication focused on the TCR/Ras/MAPK route to the IL-2 expression and the cytoskeleton remodeling response. The TCR/PI3K-Akt/mTOR route, on the other hand, is also important and has been discussed as a cluster in literature [44–47]. From Fig 3 we found the enriched pathways by T2×ST were well-grounded for the following reasons: To sum up above, we found the results of T2×ST fit our expectation: T2×ST enriched the TCR signaling pathway in the early time data, the route TCR/Ras/MAPK and TCR/PI3K-Akt/mTOR in time order, and the response of actin regulation, in both pathway databases. T2×HP, DAVID, and IPA enriched most of the expected pathways in KEGG, although the ranks of these pathways are occasionally low in IPA; T2×ST, T2×HP and DAVID also enriched most of the expected pathways in Reactome (S1 Table); GSEA failed to enrich half of the expected pathways in both pathway databases. We also found that there is no distinguishable difference between the result using a self-contained null (T2×ST, T2×HP, DPA) or a competitive null (GSEA, DAVID, IPA) in the TCR dataset. The target of the original publication is PKA substrates. As the authors described in their results; PGE2 induced a rapid and maximal increase in phosphorylation level after 1 min, and the level remained high before the number of substrates gradually returned to near-basal conditions after 60 min (S2b Fig). From Table 2 we found T2×ST, T2×HP, DPA, and IPA enriched the cAMP signaling pathway for both 1 min and 60 min experiments. The downstream of the cAMP signaling pathway was illustrated in Fig 4, there are four possible routes directed from cAMP signaling. The proliferation route was suggested by T2×ST, T2×HP, and IPA; and the cytoskeleton remodeling route also, with an addition of DPA. The DNA repair route was enriched by T2×ST, T2×HP, and DPA, in both pathway databases. The glycogen synthesis route, on the other hand, was enriched only by IPA. The reason is that T2 takes expression ratios as an important feature, whereas the mapped proteins of the “Glycolysis / Gluconeogenesis” pathway are all of low ratios (min = −0. 23710, max = 0. 14950, mean = 0. 02742). Briefly, we found the results of T2 reasonable: T2×ST and T2×HP enriched the cAMP signaling pathway for both datasets; the PKA/Rap1/PI3K-Akt route in time order, and also the cytoskeleton remodeling route and the DNA repair route. DPA and IPA also enriched most of the expected pathways. GSEA and DAVID failed to enrich most of the expected pathways in both pathway databases. We also found that the result using a self-contained null (T2×ST, T2×HP, DPA) enriched more expected pathways than a competitive null (GSEA, DAVID, IPA). Since the pathways competes with each others under a competitive null, the success of some pathways will obstacle other pathways. For example, the top 1 enriched pathway provided by DAVID for the 1 min data is “Ribosome”. The pathway includes 87 proteins, and 51 of them were mapping by the data. The cAMP signaling pathway, on the other hand, includes 73 proteins, but only 8 of them were mapping by the data. The high mapping rate of “Ribosome” will makes it harder for DAVID to enrich the target cAMP signaling pathway. This study aims to characterize the changes in protein expression underlying the phenotype conversion from mononucleated muscle cells to multinucleated myotubes. According to the their analysis, five functional clusters were identified in this dataset: cell cycle withdrawal (72hr/0hr), cell adhesion and migration (24hr/0hr), RNA metabolism (both 24hr/0hr and 72hr/0hr), myofibril formation (72hr/0hr), and proteolysis, fusion, and ECM remodeling (both 24hr/0hr and 72hr/0hr). The corresponding KEGG pathways were illustrated in Fig 5. We chose the pathway “ECM-receptor interaction” as the target pathway in Table 2 because it is clearly stated to be differentially expressed in both experiments. From Table 2 we found the ECM-receptor interaction pathway was enriched by T2×ST, T2×HP and DAVID for both experiments. The original publication focused on the change of cellular phenotype accompanying myogenic differentiation and the development of myofibril. From Fig 5 we found the myofibril formation route was enriched by T2×ST, T2×HP, DAVID, and IPA. T2×ST and T2×HP further enriched muscle contraction related pathways, which were associated with myotube maturation as discussed in the original publication. The cell adhesion and migration play an essential role in the fusion of mononucleated myoblasts. Pathways related to adhesion and migration were enriched by T2×ST, T2×HP, and DAVID; those related to fusion were enriched by T2×ST, T2×HP, DAVID, and IPA. The elevation of lysosomal proteins contributed in remodeling intracellular components during the course of myotube formation. All tools enriched the lysosome pathway with an exception of IPA. The RNA metabolism and the cell cycle withdrawal routes represented the termination of proliferation since the growth factors and nutrition were removed from the medium. Related pathways were enriched by T2×ST, T2×HP and IPA. Briefly, T2×ST and T2×HP enriched all the pathways discussed in the original publication; the interpretation made by T2 fit the description of the dataset pretty well. In the meantime, DAVID and IPA also enriched most of the expected pathways, asides from muscle contraction pathways. DPA enriched most of the expected pathways in Reactome, but failed in KEGG; GSEA failed in both pathway databases. We also found that there is no distinguishable difference between the result using a self-contained null (T2×ST, T2×HP, DPA) or a competitive null (GSEA, DAVID, IPA) in the myogenesis dataset. This dataset, unlike the previous, is not a time series; it is a dose-comparison experiment. According to the original publication, two datasets shared nearly all identified proteins although 50nM dataset did down-regulate more phosphopeptides (S2d Fig). Consequently, the authors only paid attention to the proteins that are regulated by both 5nM and 50nM dasatinib. The main target of dasatinib is the BCR-ABL signaling pathway, which is described in the pathway “Chronic myeloid leukemia (CML) ” of KEGG. From Table 2 we found the CML pathway was enriched by T2×ST, T2×HP, DAVID and IPA for both experiments. The downstream of the CML pathway was illustrated in Fig 6, there are three possible routes directed from the inhibition of BCR-ABL. Since the two datasets shared nearly all proteins, we used the ranks of 5 nM dataset to represent the common pathways over two datasets. The original publication focused on the BCR-ABL/Ras/MAPK route and the connection between BCR-ABL signaling and apoptosis. From Fig 6 we found the enriched pathways by T2 were plausible for the following reasons: In short, we found the performance of T2 pleasant: T2×ST and T2×HP enriched the BCR-ABL/Ras/MAPK route, the BCR-ABL/PI3K-Akt/Apoptosis route, the JAK-STAT signaling pathway, and the pathways related to actin response. Both DAVID and IPA enriched some of the expected pathways. GSEA and DPA failed to enrich most of the expected pathways in both pathway databases. Generally there is no distinguishable difference between the result using a self-contained null or a competitive null, but T2 enriched more expected pathways than other methods. The dataset includes only one experiment, comparing the gene expression of myocilin expressed cells to control cells, under U0126 treatment. The authors concluded that myocilin has a protective effect to against apoptosis and further promotes cell survival and proliferation via the MAPK signaling pathways. They also experimentally confirmed that the Raf-MEK-ERK-MAPK cascade was activated by myocilin. From Table 2 we found T2×ST, T2×HP, DPA, and DAVID enriched the MAPK signaling pathway; all three methods under a self-contained null successfully enriched the target pathway, whereas only DAVID is under a competitive null. The upstream and downstream of the MAPK pathway was illustrated in Fig 7. There are two possible upstream receptors of the MAPK signaling pathway: the GPCR receptors and the Notch receptors. In KEGG, both upstreams were enriched by T2×ST, T2×HP and DAVID, DPA only enriched the GPCR/Ras/MAPK route and IPA only enriched the Notch/Ras/MAPK route. In Reactome, both upstreams were enriched by T2×ST and T2×HP, DAVID only enriched the Notch/Ras/MAPK route (S1 Table). The are two downstream routes, both are supported by the original publication. The results from both T2×ST and T2×HP suggested that the differentially expressed pathways are more upstream. This conclusion actually fit the discussion of the original publication, which suggests that myocilin may also regulate the upstream kinase of the MAPK signaling pathway. Briefly, both T2×ST and T2×HP successfully enriched the MAPK signaling pathway and its upstream; other tools enriched only some of the expected pathways. Generally speaking, we found that the proposed T2-statistic was able to enrich the pathways in agreement with the original publication, whereas the performances of DPA, GSEA, and DAVID were not stable. IPA, as a commercial software with high cost, also enriched most of the relevant pathways. Nevertheless some of the pathways are low-ranked and the numbers of enriched pathways are enormous. The results suggested that our multivariate design of the proposed T2-statistic does provide important information toward pathway analysis by considering the strength of interactions among proteins. In the meantime, our self-contained null hypothesis is capable of enrich relevant pathways by the significance of protein expression ratios, whereas the focus of the competition among pathways may neglect the clear distance between phenotypes. Both DAVID and IPA are based on the competitive null hypothesis, although DAVID performs a more stringent post hoc correction, they shared the failure of some relevant pathways. The tremendous numbers of enriched pathways also suggested IPA may report more false positive results. GSEA also applies a competitive null when the sample size is limited and its KS statistic is sensitive to small sample size. The unsatisfied results suggested that GSEA is not suitable for data of limited experiments. The design of DPA is similar to the proposed T2-statistic; both DPA and T2 are specifically designed for quantitative proteomic data, and they both use self-contained null hypotheses. The performance difference between DPA and T2 primarily comes from the aspect of statistic construction. The proposed T2-statistic outperformed DPA because it considers the strength of interactions among proteins. Briefly, in five testing datasets, the results using a self-contained null is generally more well founded than the results using a competitive null. The importance of applying the covariance matrix is to estimate accurate confidence interval. We illustrated an example in Fig 8 to demonstrate the situation that may cause inaccurate estimation. Both M1 and M2 are accurate null distributions since the data are normalized using proper covariance matrix and the distribution hence follows χ2. M4 indicates that situation that an independent data are misinterpreted as a correlated data. This happens when we have false positive protein-protein interactions in the databases. In order to minimize the risk, we only use the confidence scores derived from directly experimental evidence. M3 represents the case that a correlated data are misinterpreted as an independent data. This may actually happen due to our incomplete knowledge of the biology system. In this case, the null distribution will not follow χ2 and the estimation of p-value will be inaccurate. Even so, for pathways of vary high expression ratios, applying the covariance matrix or not does not change its p-value dramatically. Here we demonstrated the impact toward p-values using permuted and purged confidence scores. We performed 100 experiments with 30% and 60% permuted confidence scores (i. e. we reassign the score using the same score distribution) and another 100 experiments with 30% and 60% purged confidence scores (i. e. we randomly remove the scores from the PPI databases), and we checked if the expected pathways are still significant under current significance level (α = 0. 05). From Table 3 we observed that: The construction of the proposed T2-statistic showed T2 is heavily dependent on expression ratios. After all, the null hypothesis for T2 is to test if the mean vector equals to zero. The contribution of applying the covariance matrix is to estimate p-values in a more accurate manner: to rescue some pathways with moderate expression ratios but their regulation directions are consistent with current knowledge of protein interaction, and to discard some pathways with inconsistency.
In this study, we presented a knowledge-based T2 approach to perform pathway analysis for quantitative proteomic data of a limited number of experiments. The proposed T2 is constructed as a multivariate statistic and the test of significance is under a self-contained null. We use the probabilistic confidence score provided by the STRING or HitPredict databases to approximate the covariance matrix of the protein profiles. The proposed T2-statistic is therefore able to reveal the influence of protein-protein interactions while performing the analysis. In addition, our pathway integration procedure is able to categorize pathways into pathway groups as well as to avoid redundancy. We performed the T2-statistic on five published quantitative proteomic dataset. In all cases, T2 was able to eliminate irrelevant pathways, as well as correctly identify relevant pathways that had been discussed in the original publication. The idea of incorporating biological evidence into conventional statistic can be widely applied to the analysis of quantitative proteomic data.
Pathway analysis is a common approach to quickly access the pathways being regulated in the experiments. There are numerous statistics to perform pathway analysis; most of them assume that the genes or proteins are independent of each other for statistical ease. This assumption, however, is unrealistic to the real biological system and may cause false positives in practice. A standard way to address this issue is to measure the associations among genes or proteins. Unfortunately, the estimation of associations requires sufficient sample size, which is usually not available for proteomic data produced by mass spectrometry. In this study, we propose a T2-statistic, which estimates the associations among gene products, to perform pathway analysis for quantitative proteomic data. Instead of calculating the associations directly from data, we use the confidence scores retrieved from protein-protein interaction databases. We also design an integrating procedure to reserve pathways of sufficient evidence as a regulated pathway group. We compare the proposed T2-statistic to other popular statistics using five published experimental datasets, and the T2-statistic yields more accurate descriptions in agreement with the discussion of the original papers.
phosphorylation
random variables
covariance
protein expression
mathematics
tcr signaling cascade
molecular biology techniques
cellular structures and organelles
mapk signaling cascades
cytoskeleton
research and analysis methods
proteins
gene expression
biological databases
proteomics
molecular biology
probability theory
molecular biology assays and analysis techniques
gene expression and vector techniques
biochemistry
signal transduction
proteomic databases
cell biology
post-translational modification
database and informatics methods
genetics
biology and life sciences
physical sciences
cell signaling
signaling cascades
A knowledge-based T2-statistic to perform pathway analysis for quantitative proteomic data