Stratified pathway analysis to identify gene sets associated with oral contraceptive use and breast cancer.

Cancer biomarker discovery can facilitate drug development, improve staging of patients, and predict patient prognosis. Because cancer is the result of many interacting genes, analysis based on a set of genes with related biological functions or pathways may be more informative than single gene-based analysis for cancer biomarker discovery. The relevant pathways thus identified may help characterize different aspects of molecular phenotypes related to the tumor. Although it is well known that cancer patients may respond to the same treatment differently because of clinical variables and variation of molecular phenotypes, this patient heterogeneity has not been explicitly considered in pathway analysis in the literature. We hypothesize that combining pathway and patient clinical information can more effectively identify relevant pathways pertinent to specific patient subgroups, leading to better diagnosis and treatment. In this article, we propose to perform stratified pathway analysis based on clinical information from patients. In contrast to analysis using all the patients, this more focused analysis has the potential to reveal subgroup-specific pathways that may lead to more biological insights into disease etiology and treatment response. As an illustration, the power of our approach is demonstrated through its application to a breast cancer dataset in which the patients are stratified according to their oral contraceptive use.

Introduction

Combining high-throughput microarray data with pathway information has proved to be a fruitful approach to uncovering biological networks from genomics data. Compared to single-gene based analysis, pathway-based methods can identify more subtle changes in expression.1,35 Furthermore, pathway-based methods can help generate biological hypotheses, which can be readily tested using complementary approaches, such as proteomics and metabolomics technologies. Each pathway used in such analyses generally serves a particular cellular or physiologic function, and these annotated pathways usually come from various external databases, such as KEGG2 and BioCarta (http://www.biocarta.com/). Wang et al.3 recently reviewed pathway-based methods, including enrichment analysis, nonparametric regression, discriminant analysis, partial least squares, and random forests. Despite these recent progresses in pathway-based analysis, not much attention has been paid to sample heterogeneity that may be partly captured by observable clinical variables. The incorporation of such clinical information may help identify relevant pathways unique to a subpopulation, leading to novel insights into disease etiology and more specific treatment schemes.

In this paper, we propose to perform stratified pathway analysis for subpopulations, with attention to the analysis results that overlap across strata as well as those that are unique to one of the stratum. If the pathways are found to be significant only to one of the subpopulations, it means that potentially those pathways are switched on for that subgroup. If the pathways are found to be significant for both subpopulations, it suggests that the genes are consistently switched on for both subpopulations, ie, regardless of the stratifier. We primarily study a breast cancer dataset to illustrate the usefulness of this approach. Estrogen receptor (ER) and progesterone receptor (PR) statuses are commonly used to estimate the risk of breast cancer, design therapy, and predict survival rate.4–10 For example, breast cancer patients are usually treated either with hormone therapy or chemotherapy depending on the hormone status of ER and PR. However, not all breast carcinomas are responsive to the treatment, and pathway analysis may help identify novel therapeutic targets and develop new agents. Pathways or genes that can predict the PR status could potentially tell us more about the biological mechanism of the disease. PRs have been found to provide prognostic information as well.11,12 Oral contraceptives are known to increase the risk of pre-menopausal breast cancer.13 Additional evidence supporting hormonal use as a confounding factor toward the risk of breast cancer was found based on the data from the Women’s Health Initiative.14 Another study found that women who have taken contraceptive pills are less likely to die of cancer or heart disease.15 This is strong evidence that oral contraceptive use can be a confounder in classifying breast cancer samples using hormonal status. Therefore, we perform stratified pathway analysis based on oral contraceptive use to explore potentially different pathways involved in breast cancer. More specifically, we aim to identify pathways involved in distinguishing PR status for users and non-users of oral contraceptives based on gene expression data. Top pathways may contain genes with expression that are good at distinguishing receptor status. For example, progesterone expression may predict progesterone status as it is related to genetic loss of heterozygosity.16 Other genes may serve as surrogates for that process.

For pathway analysis, we use the Random Forests approach, which has been found to perform well among a number of machine learning methods in pathway-based analysis.17,18

The rest of the paper is organized as follows. The detailed methodology is discussed in the Materials and methods section. In the Results section, we demonstrate the usefulness of this approach through the application of our method to a breast cancer microarray dataset. We conclude the paper in the Discussion section.

Materials and Methods

Our approach is built on the previous proposal of adopting the Random Forests approach for pathway analysis.17 We describe below how we perform stratified pathway analysis, build pathway connections, compare overlapped pathway similarities, and discover how genes are shared among them.

Stratified pathway analysis

The stratified pathway analysis considers important covariates in data analysis. In the case of the breast cancer dataset to be analyzed, we use information about the oral contraceptive use of each patient. In the simplest case, we partition all the samples into subgroups based on oral contraceptive use status, and analyze each subpopulation separately. For each pathway, we build a Random Forest to predict an individual PR+/PR− status based on the gene expression levels within this pathway. To understand Random Forests, we first need to understand how to build a classification tree.

A classification tree is built as follows:

Step (I): For each pathway, draw a bootstrap sample from the original data.

Step (II): A classification tree is grown for each bootstrap sample.

Step (IIa): At each node of the tree, select predictors (√p for classification) at random for splitting.

Step (IIb): Using the gini impurity criterion described later, a node is split using the single predictor from step (IIa). Gini impurity criterion for a binary classification problem is 1 − p12 − p22, with p1 and p2 being the proportions of individuals in class 1 and 2, respectively.

Step (IIc): Repeat steps (IIa) and (IIb) until each terminal node contains samples in the same class or only one sample.

Random Forests construct many classification trees and thus the name ‘forest’. Every tree is built using a deterministic algorithm that differs from ordinary tree algorithms in two regards. First, at each node, a best split is chosen from a random subset of predictors rather than all of them. Second, every tree is built using a bootstrap sample of the original observations. For more details, see Breiman.19 The default parameters in R’s Random Forest implementation are used, except for running 20,000 trees.

For calculating the classification error, we employ five-fold cross-validation to make it more stable across the stratified subgroup. For each cross-validation, four-fifths of the samples are used to build Random Forests, and the other samples are used to estimate the classification error. In doing so, each subject from the left out set of the cross-validation iteration is put down every tree in the forest for classification using the input vector of gene expression for genes within a particular pathway. Each tree gives a classification for this subject, and the forest chooses the class that gives the majority votes for this subject. Small classification error based on genes in a given pathway would indicate the pathway as potentially interesting.17 Multiple-testing adjusted permutation P-values are calculated. The permutation P-value is the proportion of observed cross-validation errors smaller than the cross-validation errors obtained from 500 Random Forest runs of randomly permuted labels of patients. The top pathways with a permutation P-value of less than 0.05 from both users and non-users of oral contraceptives will be presented for overlapped and non-overlapping pathways. Important genes, ie, those having strong discrimination power and selected based on the importance measure in Random Forests, are also investigated.

We will provide some biological interpretations of our results with the help of PubMatrix,20 software for comparing a list of terms against any other list of terms in PubMed.

Datasets

Pathways

A total of 285 pathways from BioCarta were used for the analysis. Most of these are related to signal transduction in humans with a smaller group of metabolic pathways.

Microarray data

A breast cancer microarray dataset accompanied by oral contraceptive use data was analyzed. This dataset used Affymetrix GeneChip® hgu-133 plus 2.0 with 54,613 probe sets. The INTEGEN (http://www.integen.org/) dataset (accession number GSE2109 from GEO) consists of 199 breast tissue samples with clinical status of PR. A total of 123 of the 199 samples were taken from patients with oral contraceptive use and the remaining 76 from patients who did not take oral contraceptives.

We chose the breast cancer dataset and PR positive/negative status (PR+/PR−) to study because breast cancer has been extensively studied in the literature and tumor samples are normally classified on the basis of hormone receptor status. A recent publication described PR as a stronger predictor of treatment response of adjuvant tamoxifen than the ER.21 More recently, additional studies suggested that measurement of PR status in conjunction with ER status may help identify patients that benefit from therapy.22,23 Daniel et al.24 provide strong rationale for targeting PR and ER in combination. The PR status has also been used to guide breast cancer therapy, breast cancer survival rate, and estimate breast cancer risk.4,5,9

Software

R was used to perform stratified pathway analysis. The R code is based on pathway analysis using Random Forests and it is available here: http://people.duke.edu/~hp44/r_code.htm. For pathway visualizations, Cytoscape25 was used.

Results

Top pathways

In this section, we show the results of the analyses. Top pathways for both users and non-users of oral contraceptives that are good at distinguishing between PR+ and PR− with permutation P-values of less than 0.05 are presented. Table 1 lists the top four pathways for non-users of oral contraceptives while the top five pathways for users of oral contraceptives are shown in Table 2. The top pathways in Table 1 are: (i) Eph Kinases and ephrins support platelet aggregation and (ii) IL 10 Anti-inflammatory Signaling Pathway. It has been found that Eph–ephrin signaling and IL 10 signaling are related to cancer26,27 and to breast cancer in particular.28 The following genes are important for the classification between PR+ and PR− in the pathways in Table 1: ADCY1, MAP3K7, and PAK2. And for Table 2, these are the important genes that are shared in some of the pathways listed: JUN, MAPK3, PIK3C, PIK3R1, and SOS1. Moreover, oral ethinylestradiol, an active estrogen compound found in oral contraceptives, decreased expression of chemokine receptors such as CCR3.29 This potentially explains why CCR3 signaling in eosinophils pathway was found as one of the top pathways among oral contraceptive users. About 27% of the important genes among the top pathways for non-users of oral contraceptives have literature citations compared with 39% of the important genes among the top pathways for users of oral contraceptives. While some pathway names and relationships to cancer are less apparent, the top genes in them may be biologically meaningful. The most important genes in Table 3 will be explored in more detail in the next section. Given that the sample size for oral contraceptives users is larger than that for the non-users, we investigated whether the sample size had an impact on the number of top pathways by sampling only 76 from the 123 oral contraceptive patients. All the pathways listed in Table 2 remain significant and have permutation P-values less than or equal to 0.05. If we did not stratify by the users and non-users of oral contraceptives, then four of the five top pathways are different from what have been found, including MAPKinase Signaling Pathway; Telomeres, Telomerase, Cellular Aging, and Immortality; CARM1 and Regulation of the Estrogen Receptor pathway; and Cell to Cell Adhesion Signaling pathway. Keratinocyte Differentiation pathway is the only one found in Table 3.

Table 1

Top pathways for non-users of oral contraceptives.

PATHWAY NAME	P-VALUE*
Eph Kinases and ephrins support platelet aggregation	0.030
IL 10 Anti inflammatory signaling pathway	0.010
Regulation of spermatogenesis by CREM	0.012
TNFR1 signaling pathway	0.018

Note:

From permutation.

Table 2

Top pathways for users of oral contraceptives.

PATHWAY NAME	P-VALUE*
CCR3 signaling in eosinophils pathway	0.012
Neuropeptides VIP and PACAP inhibit the apoptosis of activated T cells	0.004
PDGF signaling pathway	0.008
The IGF 1 receptor and longevity	0.024
Transcription factor CREB and its extracellular signals	0.010

Note:

From permutation.

Table 3

Overlapping top pathways for both users and non-users of oral contraceptives.

PATHWAY NAME	P-VALUE 1	P-VALUE 2
IL 2 receptor beta chain in T cell activation	0.006	0.002
IL 6 signaling pathway	0.022	0.032
Keratinocyte differentiation	0.002	0.002
Pelp1 modulation of estrogen receptor activity	0.004	0.002
Rho cell motility signaling pathway	0.014	0.006
Estrogen-responsive Efp controls cell cycle breast tumors growth	0.008	0.026
Role of ERBB2 in signal transduction and oncology	0.002	0.002

Note:

From permutation P-value 1 and P-value 2 for non-users and users, respectively.

Important genes and biological implications

Figure 1 contains two plots, each showing the set of pathways listed in Table 3 as well as the corresponding important genes. The top and bottom halves of Figure 1 correspond to non-users and users of oral contraceptives, respectively. The genes are hexagon shaped and are shaded according to their discriminative power in distinguishing PR+/PR− samples, with darker red indicating more discriminative power. Clearly as an overview of the plot, we can see that users of oral contraceptives have darker red genes than non-users of oral contraceptives. Moreover, it is important to note that there are genes that are important for distinguishing PR status in non-users but not in users of oral contraceptives and vice versa.

Figure 1

Top overlapped pathways for non-users (top) and users (bottom) of oral contraceptives.

Notes: Hexagon shaped are genes. Dark red as most important, white as least important.

The genes that are colored red and common among them include: ESR1, SFN, TP53, MAPK3, IRS1, IL6ST, ACTR2, ACTR3, SRC1, and PELP1. Of these genes, ESR1, IL6ST, and MAPK3 are shared among pathways, and may therefore help facilitate pathway crosstalk. Now, let us look at the unique top important genes separately. We make use of PubMatrix, web-based software, to compare a list of terms against any other list of terms in PubMed. The top important genes SOCS1, CREBBP, IKBKB, CDC2, CDK2, FAS, JUN, and FOS of the bottom pathway cluster in Figure 1 have 102 literature citations relating them to oral contraceptives from PubMatrix. Whereas for non-users of oral contraceptives, the top important genes, PRKCD, CFLAR, EGFR, BCL2L1, and ERBB2, only have 30 literature citations relating them to oral contraceptives from PubMatrix. Interestingly, these genes are unique among pathways and are not shared. A Mann–Whitney U test gives a two-sided P-value of 0.075 when comparing the number of literature in user and non-user groups for the top five genes. Furthermore, using GeneGO, we identified therapeutic drug targets for eight of the top genes identified, see Table 4.

Table 4

Therapeutic drug target for top genes.

KNOWN DRUGS (EXCLUDING EXPERIMENTAL)	GENE SYMBOL
Afatinib, Neratinib, Lapatinib, Canertinib, Gefitinib	ERBB2 & EGFR
Roniciclib, Alvocidib	CDK1 & CDK2
Erlotinib, Falnidamol, Vandetanib, Genistein, Cediranib, Varlitinib, Pelitinib, Suramin	EGFR
Ingenol	PRKCD
Raloxifene, Afimoxifene, Megestrol, Diethylstilbestrol, Clomifene	ESR1
Riviciclib	CDK1
Bardoxolone, Xanthohumol	IKBKB
Seliciclib	CDK2
Navitoclax	BCL2L1
Masoprocol	ERBB2

Discussion

Stratified pathway analysis refers to performing pathway analysis and measuring prediction accuracy within subgroups or strata of the experimental population. Clinical covariates that may confound the pathway analysis or gene-set enrichment analysis should be taken into account. If necessary, stratified pathway analysis should be performed. In this paper, we highlight the importance of incorporating covariates in pathway analysis. We have described a Random Forests-based approach to perform stratified pathway analysis in distinguishing PR positive and negative breast cancer patients. The novel method presented in the article was able to identify unique pathways specific to oral contraceptives users and non-users as well as shared pathways among those groups. In addition, we were able to tease out the important genes that relate to outcome of interests that are biologically meaningful. Although we used Random Forests for classification, other methods, such as support vector machines, can also be employed here for stratified analysis.

We have demonstrated the biological relevance of our approach using PubMatrix. The number of important genes identified with literature agreed well with those shortlisted by our analyses. Furthermore, with the aid of network visualization tools, we can allow biologists to investigate how the important genes are related to each other within a set of pathways. One limitation of our approach is that the computational intensity of our approach increases linearly with the number of levels of the confounding variable.

Bioinformaticians and biologists can make use of this method to analyze specific subgroup of patients, focus on a few sets of genes, identify pathway targets, and find out how important genes are shared among the top pathways. This allows researchers to obtain results that are more closely tied to the biological mechanism of diseases. This analysis can also be applied to ER status. In this case, a weighted random forests algorithm should be used to deal with the unbalanced proportion of ER positive and negative groups.30 Other machine learning gene selection strategies may be incorporated.31,32 Moreover, it may be applied to genotyping data as well.33,34