A Steiner tree-based method for biomarker discovery and classification in breast cancer metastasis

Datasets and PPI networks

In this study we used two microarray datasets herein referred as van de Vijver and Wang dataset [4, 5] respectively. The two datasets have 295 and 286 breast cancer patients where 78 and 106 patients have distant metastasis within five years of follow-up visit respectively. The microarray platform used for van de Vijver et al was Agilent Hu25K and for other dataset was Affymetrix HG-U133a. The first dataset was downloaded from the Netherland Cancer Institute website (http://bioinformatics.nki.nl/index.php) while the other dataset was obtained from GEO with the accession number GSE2034 [16]. SAM (Significant Analysis of Microarray) [17] was used to select genes that are significantly differentially expressed between metastatic and non-metastatic tumors (DE genes). We controlled the delta parameter in SAM to select a similar number of DE genes from each dataset. As a result a total of 333 and 319 DE genes were selected for van de Vijver and Wang datasets, corresponding to FDR 0.7% and 8.2% respectively. Varying the number of DE genes between 200 and 1000 only slightly changed the percent of overlap while the significance of the overlap is essentially not affected (data not shown).

Two human protein-protein interaction networks were used by this study. The first network was obtained from Protein Interaction Network Analysis (PINA) and contains 10,920 genes and 61,746 binary connections [12]. The second network was compiled by [10] from six different sources, and contains 57,235 interactions among 11,203 genes. In this study we only considered the largest connected component in each PPI network, which contains 10,794 genes and 56,864 connections for Chuang PPI network and 10,770 genes and 61,658 connections for PINA PPI network.

Randomized Steiner tree approach

The main idea of our approach (see Figure 1) is to identify a small number of genes that connect all DE genes in a PPI network. In graph theory this problem is known as Steiner tree problem. The Steiner tree for an edge-weighted graph G= (V, E, w) and a subset of vertices R ⊆ V is a minimum-weight connected tree T, with vertices U ⊆ V and edges S ⊆ E that spans all vertices in R. Here the vertices in R are known as terminal vertices and U\R as Steiner vertices. For an unweighted graph G, the problem then becomes to find the minimum number of vertices that can connect all the vertices in R through a tree in G. The Steiner tree problem is NP-hard [18]. We implemented a slightly improved version of a polynomial-time 2-approximation shortest path heuristic algorithm [19]. This algorithm starts with a forest $\left(T^{{}^{\prime }}\right)$ comprising the terminal vertices R. Then in each iteration it finds the two vertices in $T^{{}^{\prime }}$ that are closest in distance and adds the intermediate vertices to $T^{{}^{\prime }}$. This procedure is repeated until $T^{{}^{\prime }}$ becomes connected. Finally, a minimum spanning tree of those selected vertices are built and all leaf vertices that are not in R are removed. The running time of this algorithm is O(|V|³). The pseudocode is given below.

Steiner tree algorithm

Input: Weighted PPI network, G = (V,E,w); DE genes, R

In a Steiner tree the intermediate vertices (both Steiner vertices and non-leaf DE genes) play important roles in connecting DE genes together. We consider all the intermediate vertices as potential biomarkers (S teiner T ree-based M arkers or STM s) for breast cancer metastasis, as it is known that disease related genes play an important role in connecting DE genes in PPI network [10]. For example, for van de Vijver dataset with Chuang-PPI network, a single Steiner tree uses 136 vertices to connect 333 DE genes (Figure 1B-C). Among the DE genes, 264 are leaf nodes (degree = 1 in the tree) and the other 69 are internal nodes (degree > 1 in the tree). As these internal DE genes are important in connecting the remaining DE genes, we combine the Steiner vertices with these internal DE genes as potential biomarkers for breast cancer metastasis. Thus for this single tree we consider those 205 internal genes as potential biomarker (Figure 1B-C).

Next, we proposed a simple strategy to obtain multiple Steiner trees. The motivation is two-fold. First, as the heuristic algorithm does not guarantee optimality, by obtaining multiple solutions we increase the chance of finding the optimal Steiner vertices. Second, multiple solutions with similar qualities may represent alternative or redundant pathways that cannot be covered by a single Steiner tree. To obtain alternative Steiner trees without any modification to our Steiner tree algorithm, we assign to each edge in PPI network a random weight between 0.99 to 1 and run the standard Steiner tree algorithm. These random weights effectively break ties, so that if there are two paths with the same weight in the original network, one path will be chosen randomly. This procedure was repeated multiple times with different random weights from 0.99 to 1, until the total number of unique STMs converges approximately. Depending on the PPI network and microarray data, the rate of new coming STMs reduced significantly after 200-300 iterations (for example, see Figure 1D). After that, we take union of all internal nodes (or genes) of those trees and consider them as potential biomarkers. As previously mentioned, we called these genes as S teiner T ree-based M arkers (STM s). We obtained 1047 and 1100 STMs for Chuang PPI network and, 932 and 1135 STMs for PINA PPI network for van de Vijver and Wang dataset respectively (see Additional file 1 for complete gene list).

Statistical test of overlap significance

where C(n, k) is the binomial coefficient.

Classification

To evaluate the prediction ability of different features (STMs and DEs), we built logistic regression and support vector machine (SVM) classifiers to distinguish breast cancer patients who developed metastasis within five years after the date of the initial diagnosis from those who did not. We used the implementation in WEKA (version 3.6.3) and default parameter settings for this classification purpose [20]. To avoid overfitting and provide a realistic evaluation, we concentrated on cross-data classification where features obtained from one dataset were used to construct classifiers for the other dataset, because DEs were selected using the complete dataset and very specific to that particular dataset. Classification performance was estimated 100 times using 10-fold cross validation where iteratively one-tenth of the data were used for testing and nine-tenth were for training. Performance was measured by AUC (area under ROC curve).

Stability of STMs

Functional enrichment and pathway analysis of STMs

To reveal the biological functions of the candidate biomarkers we used online bioinformatics resource DAVID [21] to analyze the enriched gene ontology biological processes and KEGG pathways for STMs and DEs. We used the genes in the largest connected component of the PPI network as the background gene list in DAVID. This analysis showed that genes selected by our method are significantly enriched in functional terms that are well known to be involved in cancer carcinogenic process, including cell cycle, apoptosis, DNA repair as well as MAPK, ErbB and P53 signaling pathways. A list of the enriched terms and their enrichment scores are listed in Figure 2. These enriched functional terms and pathways have strong consistency with several previous studies [10, 22, 23]. DE genes clearly have much lower significance compared to STMs. For example, for van de Vijver dataset and Chuang PPI network DE genes are only weakly enriched in DNA repair (p-value<1.1E-1 vs 1.1E-10 in STMs) and apoptosis (p-value<4.2E-2 vs 3.0E-8 in STMs), and not significantly enriched in cell growth, ErbB, MAPK and P53 signaling pathways. A similar trend is found in Wang dataset (Figure 2). Using PINA PPI network for these two datasets we also have similar and significant results for STMs compared to DE genes (see Additional file 2). Thus starting with the DE genes with low significance our method is able to find a set of genes which are overrepresented in many known pathways involved in breast cancer metastasis.

STMs correspond to novel biomarkers of cancer

Next, to demonstrate that our method is able to identify cancer susceptibility genes, we compare our results with Chuang et al [10]. In that study they collected 60 known breast cancer marker genes from Online Mendelian Inheritance in Man (OMIM). As shown in Figure 3, for both datasets STMs have higher percentage of overlaps with those known markers than the genes selected by previous studies [4, 5, 10]. Also, DE genes cover a smaller portion of known cancer susceptibility genes than STMs. The enrichment of known breast cancer markers in STMs (p-value < 1.32E-09 and 5.28E-04, for van de Vijver and Wang datasets respectively, Fishers' exact test) are also more significant than in DE genes and genes selected by previous studies [4, 5, 10]. Some well-known biomarkers covered by STMs are BRCA1 (breast cancer 1, early onset), BRCA2 (breast cancer 2, early onset), ATM (Ataxia telangiectasia mutated), TP53 (tumor protein p53), ErbB2 (a.k.a. HER2, Human Epidermal Growth Factor Receptor 2), TFN (tumor necrosis factor) and Esr1 (Estrogen receptor 1) (for a complete list see Table 3). As the STMs are enriched in breast cancer related pathways and cover a significant number of known susceptibility genes, the rest of the STMs that are not known biomarkers can be considered as novel potential biomarkers for breast cancer metastasis.

We also collected 288 breast cancer susceptibility genes from Genetic Association Database of Disease from DAVID. STMs covered 26.39% and 26.04% of those known breast cancer susceptibility genes for Wang and van de Vijver datasets, respectively. For both datasets, genes selected by our method clearly outperformed DE genes (5.55% and 5.28% for Wang and van de Vijver datasets respectively). Evaluation results using PINA PPI show similar results.

Steiner tree-based markers improves the cross-dataset classification accuracy

We tested whether the STMs can be used to improve the prediction accuracy of breast cancer metastasis. To this end we use the STMs as features to train logistic regression and support vector machine (SVM) classifiers [20, 24] to separate metastatic from non-metastatic patients (see method section). For comparison, we also constructed classifiers using only DE genes and DE+STM genes. To evaluate the power of the candidate markers in predicting metastasis in an unbiased fashion, we focused on cross-dataset tests, where the genes selected from one dataset were used to construct classifiers for the other dataset. Figure 4 shows the classification accuracy, measured by AUC (area under ROC curve), achieved by different feature sets for logistic regression and support vector machine classifier. As can be seen, STMs resulted in better accuracy than DEs for both cases. On van de Vijver dataset, STM-based logistic regression classifier performed better than DE-based classifiers in 79 of the 100 runs (p-value = 4E-12, paired t-test), and STM-based SVM classifier outperformed DE-based SVM in 87/100 of runs (p-value = 6E-20). For Wang dataset, the corresponding numbers are 99/100 (p-value = 7.7E-44) and 86/100 (p-value = 5.5E-17). On the other hand, STMs performed better than STM+DE feature set in all cases except in van de Vijver dataset for SVM classifier. In that case STM+DE sets performed slightly better than the STMs. For all cases STM+DE set performed better than DEs. Thus STMs have better classification accuracy on predicting breast cancer metastatic potentials than DE genes.