BD-Func: a streamlined algorithm for predicting activation and inhibition of pathways

BD-Func algorithm

The basic principle behind BD-Func is to treat activated and inhibited genes as replicate observations in two populations. BD-Func is agnostic towards the type of signal used for analysis. For example, this paper uses BD-Func to study both fold-change values between two populations as well as raw signal intensities for a single column of signal values. In the paper, the t-test statistic is used to compare the expression patterns of activated and inhibited genes. However, BD-Func also allows users to compare the activated and inhibited distributions using a Mann–Whitney U test or Kolmogorov–Smirnov (K–S) test. BD-Func users have the option to calculate False Discovery Rates (FDR) using the method of Benjamini & Hochberg (1995) or the Storey q-value (Storey & Tibshirani, 2003).

BD-Func comes with four enrichment files: c2, c5, and c6 from MSigDB (Liberzon et al., 2011) and a list of functions defined directly from the human Gene Ontology (GO) annotations (Ashburner et al., 2000). All of these lists were created by searching for functions with both “up” and “down” (or “positive regulation” and “negative regulation”) gene lists. For the human GO file, a functional annotation needed to contain at least 10 positively regulated genes and 10 negatively regulated genes. We also encourage users to share their own custom models on the BD-Func discussion group: http://sourceforge.net/p/bdfunc/discussion/.

There are three different input files that can be analyzed using BD-Func (Fig. S1):

1-D Input File: If the user supplies expression values for a single column of data, a density plot is created for the signal for the activated and inhibited genes. In this case, the BD-Func algorithm works exactly as described above.

2-D Input File: If the input file contains multiple columns of data, BD-Func is first run separately for each sample (represented by a column in the data matrix), as described above. Next, box-plots are created for test-statistic scores for each group (labelled in the header of the input file; Fig. S2). Finally, an ANOVA p-value is provided to compare the test statistics between groups. Theoretically, test statistics could be used to make functional predictions for each sample in isolation. There are some examples in this paper where this strategy works OK. However, comparing test statistics across all samples within different groups is the only strategy that consistently works well for all the analysis presented in this paper.

2-D Input File for Classifier: If the input file contains multiple columns of data and two groups (called “positive” and “negative”), BD-Func will create a receiver operating characteristic (ROC) plot using the test-statistic as the classification score using the ROCR package (Sing et al., 2005). This is in addition to all the calculations and output files for a normal 2-D input file containing multiple columns for BD-Func to analyze (in this paper, this represents normalized signal intensity values).

Sample acquisition and processing

Microarray datasets were downloaded from GEO (Edgar, Domrachev & Lash, 2002) or Array-Express (Parkinson et al., 2009). When raw CEL files were available, samples were RMA normalized (Irizarry et al., 2003). Otherwise, processed intensity values were used for microarray analysis. Fold-Change values for all of the cell line datasets (TGFβ (Padua et al., 2008; Qin et al., 2009; Renzoni et al., 2004; Sartor et al., 2010; Scandura et al., 2004), mTOR (Wei et al., 2006), p53 (Elkon et al., 2005), and BRCA1 (Furuta et al., 2006)) and the MSigDB progesterone receptor dataset (Claus et al., 2008) were calculated using the method of the least-squares mean using Partek^® Genomics Suite™ (Partek Inc, 2012).

All other clinical samples (Anders et al., 2008; Bild et al., 2006; Chin et al., 2006; Finak et al., 2008; Huang et al., 2003; Ivshina et al., 2006; Sotiriou et al., 2003; The Cancer Genome Atlas Network, 2012) were downloaded and analyzed for differential expression using BRAVO (http://bravo.coh.org/) (X Deng, C Warden, Z Liu, I Zhang, Y-C Yuan, unpublished data). The novel progesterone receptor gene signature presented in this paper was produced by identifying genes in the expO dataset (GEO Series GSE2109) with a |fold-change| > 2 and an False Discovery Rate (FDR) < 0.05 (where the FDR is calculated using the method of Benjamini & Hochberg (1995) to analyze the distribution of t-test p-values). This is how the samples in this particular paper were processed, but users are not required to use this particular set of tools for preparing BD-Func input files and/or creating gene lists for custom signatures.

Reads Per Kilobase per Million mapped reads (RPKM (Mortazavi et al., 2008)) values for RNA-Seq data was downloaded from the TCGA web portal (The Cancer Genome Atlas Network, 2012). RPKM values were transformed by addition of 0.1 (to avoid large fold-change values for low coverage reads) followed by a log₂ transformation (to normalize the signal distribution).

GSEA comparison

With the exception of the progesterone receptor signature (which utilized the CLAUS_PGR_POSITIVE_MENINGIOMA (Claus et al., 2008) signature from MSigDB-c2, version 3.1), oncogenic signatures were defined using the following gene lists from MSigDB-c6 (Liberzon et al., 2011, version 3.1) for GSEA analysis: TGFB_UP.V1 (Padua et al., 2008), MTOR_UP.N4.V1 (Wei et al., 2006), P53_DN.V2 (Elkon et al., 2005), and BRCA_DN.V1 (Furuta et al., 2006). GSEA (Subramanian et al., 2005, version 2.0) calculated p-values by permutation over phenotypes whenever possible (Anders et al., 2008; Chin et al., 2006; Claus et al., 2008; Elkon et al., 2005; Finak et al., 2008; Huang et al., 2003; Padua et al., 2008; Sartor et al., 2010; Scandura et al., 2004; The Cancer Genome Atlas Network, 2012; Wei et al., 2006), although there were a few datasets with less than 3 replicates for which gene sets had to be permuted instead of phenotypes (Furuta et al., 2006; Qin et al., 2009; Renzoni et al., 2004). For recovery of known perturbations of oncogenic regulators, GSEA results must either show a FWER p-value < 0.25 or a NOM p-value < 0.05, which are the default cut-offs.

IPA comparison

Ingenuity Pathway Analysis (IPA; Ingenuity^® Systems, www.ingenuity.com) contains an “Upstream Regulator” module that compares the enrichment of activated and inhibited genes among up- and down-regulated genes. So, the underlying principle is similar to BD-Func except it utilizes Ingenuity’s propriety database of regulatory interactions and uses a z-score to calculate statistical significance between activated and inhibited genes. Gene lists in IPA were filtered for those genes showing |fold-change| > 1.5 while the entire gene list is used to define background enrichment. For recovery of known perturbations of oncogenic regulators, the upstream regulator must be identified as “activated” or “inhibited” (|z-score| > 2), which are the default cut-offs.

LBH589 signature

Activated and inhibited genes were defined using overlapping gene lists from 3 cell line treatments that have been previously published (Kubo et al., 2013). That same study showed that LBH589 treatment significantly decreased tumor volume in exemestane (EXE) resistant MCF-7aro xenografts in mice. This study analyzes novel microarray data from EXE-resistant tumors treated with (EXE + LBH589) or without (EXE only) treatment of LBH589. All animal research procedures were approved by the City of Hope Institutional Animal Care and Use Committee. This novel microarray data is available in GEO series GSE47346.

In order to be included in the novel BD-Func signature, genes must show differential expression in all 3 cell lines. Genes were defined as differentially expressed if they showed a |fold-change| > 1.5 and p-value < 0.05 , and the BD-Func signature genes had to meet these conditions for each of the 3 LBH589 cell line treatments (with consistent direction of fold-change). P-values were calculated via 1-way ANOVA with appropriate linear contrast was used to compare data sets using Partek^® Genomics Suite™ (Partek, Inc., St. Louis, MO). Fold-change values were calculated based upon the least-squares mean, and data was normalized using robust multichip average (RMA) normalization (Irizarry et al., 2003).

BD-Func shows equal or greater performance to GSEA and IPA for functional enrichment

Given the relative ease by which samples can be classified as having positive or negative activity for an individual biomarker, the accuracy of BD-Func was first tested by applying several MSigDB oncogenic signatures (Liberzon et al., 2011) to the datasets from which the signatures were defined (Claus et al., 2008; Elkon et al., 2005; Furuta et al., 2006; Padua et al., 2008; Wei et al., 2006). BD-Func was able to detect the activation or inhibition of all of these oncogenic signatures (Table 1, Fig. 2). GSEA could detect all of the signatures except the Claus et al. (2008) progesterone receptor signature. Among these 5 test datasets, IPA could only detect the activity of 2 of these genes; however, this is not a completely fair comparison because we would expect some overfitting of the MSigDB signatures for the GSEA and BD-Func analysis. Nevertheless, the significance of this analysis is that BD-Func can accurately detect perturbation of all of these biomarkers on datasets where we know that these specific genes will be altered.

In order to test the performance of BD-Func, GSEA, and IPA on novel datasets (which were not used by MSigDB to define gene lists), we applied the 3 algorithms to four datasets with TGFβ treatments (Table 2). All 3 algorithms showed roughly equal performance for predicting TGFβ treatment in the appropriate samples. Overall, analyses of these nine datasets indicate that BD-Func can provide similar quality results as GSEA and IPA.

BD-Func can determine the accuracy of novel predictive models

Six breast cancer datasets were also used to test the robustness of the progesterone receptor (PGR) signature (Anders et al., 2008; Chin et al., 2006; Finak et al., 2008; Huang et al., 2003; The Cancer Genome Atlas Network, 2012). Unfortunately, neither BD-Func, GSEA, nor IPA could predict progesterone status in all seven patient populations (Table 3). To be fair, the original Claus et al. (2008) dataset was used to define progesterone receptor status in meningioma samples whereas the novel datasets tested were all breast cancer samples (where testing for over-expression of progesterone receptor is common (Bardou et al., 2003)). Nevertheless, BD-Func is designed to be able to utilize custom gene signatures with activated and inhibited, so we defined a novel progesterone receptor signature using the expO dataset (GEO Series GSE2109). This signature can identify progesterone receptor positive and negative patients for all 7 cohorts (1 meningioma and 6 breast cancer datasets), so it is robust enough for application to multiple cancer types.

Another unique feature of BD-Func is the ability to use the activation versus inhibition test statistic as a classifier to define a predictive model. If a t-test statistic of 2 is used for the cut-off of distinguishing the positive and negative populations (roughly corresponding to a p-value < 0.05), it is clear that the MSigDB signature is extremely accurate at predicting PGR status in the original dataset but not in the breast cancer datasets (Table S1). Likewise, the TGFβ signature could differentiate between the treated and untreated groups if the test statistic of 2 was used as the threshold to distinguish the groups (Fig. S2). However, this threshold does not work well in all circumstances: unlike the analysis of fold-change values, the p-value (for any statistical method) is not always the ideal statistic for assessment of functional enrichment on intensity values. For example, the mTOR and BRCA1 signatures (Fig. 2B) show appropriate changes in test statistics that clearly distinguish treated and untreated groups, but activation and inhibition can’t be defined based upon a pre-defined cut-off for the test statistic value (e.g., 2 or −2). For this reason, we provide an ANOVA p-value to quantify the difference in test-statistic between groups, where the test statistic serves as a score for a second calculation of statistical significance.

Additionally, we believe that predictive statistics are a useful method for accessing BD-Func scores for individual samples within large patient populations. In order to quantify the accuracy of the model without depending on a pre-defined cut-off, BD-Func produces receiver operating characteristic (ROC) plots for each cohort and the area under the curve (AUC) is calculated for each ROC plot (where a 100% accurate model would have an AUC = 1.00) (Fig. S3, Table S1). The superior performance of the novel PGR signature on the breast cancer cohorts becomes even clearer when these predictive statistics are compared for the two models (Fig. 3, Fig. S4, and Table S2). Importantly, the novel signature showed the same level of accuracy on the TCGA breast cancer dataset as the original expO dataset. This is significant because the TCGA dataset is over twice as large as the expO dataset, and the TCGA dataset utilizes RNA-Seq while the expO dataset utilizes microarrays to quantify gene expression. In other words, this shows that BD-Func is capable of producing very robust predictions that translate across different genomic technologies.

BD-Func accurately quantifies LBH589 activity using a novel signature

LBH589 (panobinostat) is a histone deacetylase inhibitor that has been previously shown to suppress the proliferation of aromatase inhibitor resistant breast cancer cells, which was a conclusion supported in part by functional enrichment analysis of commonly affected genes from 3 cell line experiments (Kubo et al., 2013). Gene lists derived from these cell line experiments can be easily used to define a BD-Func signature, so we hypothesized that the results from this previous in vitro cell culture study could be used to predict LBH589 activity in vivo in an animal study. Specifically, we asked if the signature defined based upon LBH589 treatment in 3 cell lines (H295R, MCF-7her2, HeLa) could detect LBH589 activity in a mouse xenograft from a different cell line (MCF-7aro xenograft treated with EXE). Indeed, BD-Func correctly used the cell line LBH589 signature to identify common gene expression changes in the tumours treated with LBH589 and EXE compared to the mice that were only exposed to EXE treatment (Fig. 4).