text:
Detecting gene signature activation in breast cancer in an absolute, single-patient manner

Background The ability to reliably identify the state (activated, repressed, or latent) of any molecular process in the tumor of a patient from an individual whole-genome gene expression profile obtained from microarray or RNA sequencing (RNA-seq) promises important clinical utility. Unfortunately, all previous bioinformatics tools are only applicable in large and diverse panels of patients, or are limited to a single specific pathway/process (e.g. proliferation). Methods Using a panel of 4510 whole-genome gene expression profiles from 10 different studies we built and selected models predicting the activation status of a compendium of 1733 different biological processes. Using a second independent validation dataset of 742 patients we validated the final list of 1773 models to be included in a de novo tool entitled absolute inference of patient signatures (AIPS). We also evaluated the prognostic significance of the 1773 individual models to predict outcome in all and in specific breast cancer subtypes. Results We described the development of the de novo tool entitled AIPS that can identify the activation status of a panel of 1733 different biological processes from an individual breast cancer microarray or RNA-seq profile without recourse to a broad cohort of patients. We demonstrated that AIPS is stable compared to previous tools, as the inferred pathway state is not affected by the composition of a dataset. We also showed that pathway states inferred by AIPS are in agreement with previous tools but use far fewer genes. We determined that several AIPS-defined pathways are prognostic across and within molecularly and clinically define subtypes (two-sided log-rank test false discovery rate (FDR) <5%). Interestingly, 74.5% (1291/1733) of the models are able to distinguish patients with luminal A cancer from those with luminal B cancer (Fisher’s exact test FDR <5%). Conclusion AIPS represents the first tool that would allow an individual breast cancer patient to obtain a thorough knowledge of the molecular processes active in their tumor from only one individual gene expression (N-of-1) profile. Electronic supplementary material The online version of this article (doi:10.1186/s13058-017-0824-7) contains supplementary material, which is available to authorized users.

Background
There are many pathway analysis approaches that seek to determine whether a specific molecular process or cellular response is activated, repressed, or latent in a given patient sample including pathway level analysis of gene expression (PLAGE) [1], zscore [2], single sample gene set enrichment analysis (ssGSEA) [3], functional analysis of individual microarray expression (FAIME) [4], gene set variation analysis (GSVA) [5], Pathifier [6], N-of-1 [7] and diversity arrays technology (DART) [8]. To make predictions, all of these methods require a suitable database of gene signatures, where each signature is composed of a set of genes with expression levels that correlate with different activation states of the specific molecular process or biological response. Typically, these approaches return a score that measures the level of activation of the process or response.
There are important limitations associated with these existing approaches. First, the vast majority (ssGSEA, GSVA, Pathifier, FAIME and DART) generate scores that are not easily interpreted in isolation. More specifically, the score for an individual pathway in a tumor must be interpreted relative to the distribution of scores for all patients in a large cohort. This is certainly problematic for situations where only the sample from the target patient is available without recourse to a large dataset of patients for interpretation, a scenario that often occurs in clinical investigations and in the clinic.
Second, the relativistic nature of some of the existing methods makes them highly sensitive to the specific composition of the patient's dataset [9,10]. For example, scores for a given pathway in a specific patient will vary depending on the clinicopathological and molecular characteristics of the other tumors/patients in the dataset. This is analogous to the instability that was previously described for methods that determine breast cancer subtype when the fraction of patients with estrogen receptor (ER)-positive or HER2-positive cancer varies between datasets [9,10]. Techniques such as N-of-1 [7] are potentially not as prone to these problems. However, they require both a normal and a tumor sample from the same patient, which is difficult to obtain and is not yet standard clinical practice.
Ideally, a tool to measure the activation of a pathway via a gene signature should have the following four essential properties: (1) it should be applicable to gene expression profiles of single patients without recourse to or need for a larger cohort for comparative analysis; (2) the calls should be "stable", i.e. predictions should not be influenced by the composition of a comparative patient cohort; (3) it should function across a large range of platforms (e.g. various microarrays and RNA-seq technologies); and (4) the prediction of the state of each gene signature in a patient's tumor should be justified statistically. Unfortunately, to our knowledge, there is no approach currently available that meets all these four criteria.
Here we present an approach entitled absolute inference of patient signatures (AIPS) that builds upon our previously described absolute intrinsic molecular subtyping (AIMS) method to predict the breast cancer subtype in a manner respecting the aforementioned properties 1-4 [9]. The absolute, stable and general properties of AIMS are generalized from markers of subtype classification to a large set of gene expression signatures representing an extremely diverse array of biological processes, pathways and functions, encompassing the hallmarks of breast and other cancers. The tool allows the molecular dissection of a single tumor for all relevant biological processes that are activated or repressed at the transcriptional level, and greatly extend the current technologies that measure just a single process (e.g. proliferation, immune-related information, etc.).

Methods
Supplemental methods are provided in Additional file 1.
Testing the stability of current approaches for assigning pathway activation scores This analysis took advantage of the implementations of ssGSEA, PLAGE, zscore and GSVA in the GSVA Bioconductor package [5]. We restricted our attention to the 4725 gene signature in C2 (curated gene sets) from the MSigDB [11]. Using the complete Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) gene expression dataset (n = 1992), we assigned gene signature activation scores using one of the four approaches for all gene sets in the C2 collections to all patients in METABRIC. This produces a matrix M ALL in which M ALL (P i ,GS j ) corresponds to the score assigned to patient P i using gene set GS j . We repeated the score assignment a second time by applying the approach to only the ER-positive patients in the METABRIC dataset and obtained a second matrix of scores M ER+ in which M ER+ (P i ,GS j ) corresponds to the score assigned to patient P i using gene set GS j . We repeated this analysis for the ER-negative patients to obtain a third matrix M ER-. The distributions presented in Fig. 1a are obtained from the absolute difference between the M ALL and M ER+ or M ERmatrix. For example, the distribution for the ER-positive patients are computed by |M ALL (P i ,GS j ) -M ER+ (P i ,GS j )| for all ER-positive patients and all the gene sets GS j .

Preparation of gene expression datasets
We utilized the same training set as described in Paquet et al. [9] with the exception of the dataset of Loi et al., because after careful examination, this dataset was mainly composed of ER-positive samples and had evidence of batch effects with ER-negative samples (Table 1). Briefly, we generated the training set by downloading all of the original normalized dataset from the Gene Expression Omnibus (GEO) and ArrayExpress ( Table 1). The normalized data are used directly for the region of independence (ROI) 95 assignments (see Additional file 1 for a description of the ROI 95 ), but as AIPS requires quantification of raw gene expression, we also downloaded the raw gene expression data and pre-processed them to remove background artifacts when applicable (mostly for the microarray platforms).
The McGill validation dataset was generated on the Human Affymetrix Gene ST platform as previously described in Tofigh et al. [12]. For the METABRIC dataset we kept the 12 replicate samples described by the authors of the original publication [13]. As these correspond to less than 0.2% of our training set, their inclusion does not significantly affect any of the presented results. The final analyses and models were restricted to the Entrez IDs present on all the platforms in the training and validation datasets. When multiple probes map to the same Entrez ID, the ROI 95 assignments used the most variable probe (using the interquartile range (IQR)). Unfortunately, the AIPS models cannot rely on such an approach because the assignments must be performed in the context of only a single sample. Therefore, the IQR cannot be used. The AIPS models used the probe with the highest raw gene expression in downstream analyses. We favored this solution over taking the mean of all probes because we believe that using the mean would bias the approach by favoring genes with more probes, as they would have less variance. On the other hand, it could also introduce noise because we cannot really expect all the isoforms of a gene to have similar levels of expression.

Assembling a large collection of harmonized gene signatures
We collected 6466 gene signatures from different sources including MSigDB (C1 to C7, n = 6183 [11]), GeneSigDB (n = 188, [14]), and breast-cancer-specific gene signatures obtained from the literature (n = 95). The number of signatures from each source is depicted in Fig. 2a. Entrez IDs were used in order to guarantee distinct, unambiguous gene names, identifications and symbols. As best possible, we determined the "directionality" of expression for each gene within each signature, and used this information to A B C D Fig. 1 Instability of current pathway activation inference tools and example application of the region of independence (ROI). a Absolute difference between the score obtained from a specific pathway activation tool (gene set variation analysis (GSVA), single sample gene set enrichment analysis (ssGSEA), zscore, or pathway level analysis of gene expression (PLAGE)) using all patients from the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) dataset, and the score obtained when this dataset is restricted to either estrogen receptor (ER)-negative (ER-) samples (left) or ER-positive (ER+) samples (right). b Heatmap depicts the ROI induced over the rank-sum-base ordering of patients in the METABRIC dataset using the estrogen activation gene signature from Doane et al. with the probability (1) and cumulative distribution function (2) for the random trials, in addition to the final assignments (3) into low, high and independent regions (defined as the 95% CI of the index of the random trials). c Distribution of the low, independent and high assignments defined by the ROI 95 in the function of the clinical subtypes (defined by ER and human epidermal growth factor receptor 2 (HER2) status). d Class assignments defined by the ROI 95 are prognostic with 5 years survival (log-rank test P < 0.00001) partition the genes into the overexpressed and underexpressed subsets.

Development of AIPS
The pipeline to generate AIPS models is summarized in Fig. 3a and contains two major steps: winnowing to informative gene signatures in breast cancer, followed by the selection of reliable absolute models.

Winnowing to informative gene signatures in breast cancer
The numbers refer to Fig. 3a.
(1)We applied the ROI 95 on the 6466 gene signatures in our collection using the largest of our gene expression datasets (METABRIC, n = 1992). The size of the dataset helps to ensure that the ROI 95 will have sufficient samples to precisely identify the region of independence (see Additional file 1 for a description of the ROI 95 ) (2)To estimate significance, we also performed the same analysis but first randomly permuted the genes in the METABRIC dataset prior to applying the ROI 95 .
In this manner, we are able to identify those gene signatures that behave essentially like random gene signatures. For every gene signature in the non-permuted versus permuted dataset we generated distributions from the percentage of patients assigned to the low, independent or high classes ( Fig. 3b) (3)We selected a value of 0.8 to apply to the percentage of patients in the independent class as a cutoff to filter out non-informative gene signatures. This conservative cutoff still induces a sufficient number of patients within each of the three classes (low, independent, and high) to enable training of our models. A 0.8 cutoff on the percentage of independent assigned patients for a given signature retains 3472 gene signatures, which we consider to be informative in breast cancer (Fig. 2c).
Building and selecting reliable absolute models from the informative gene signatures   Table 1) and the diverse range of platforms in our training set rules out an approach that builds a merged normalized meta-dataset  (5)Step (4) provides statuses (low, independent, and high) for 3472 signatures in each of 4510 patients. A model is learnt for each such signature as follows. First, we consider every possible pair of genes x and y in the expression dataset, and we ask how strongly the rule "if expression of x is greater than y, then classify signature as high" is associated with the ROI 95 assignments. This is repeated for a modification of the rule where x is less than y, and for the other activation statuses of independent and low. This process is performed on the entire training set and weights are assigned to every patient in a way to give equal weight to the individual datasets i.e. patients from smaller datasets are given more weight than patients from larger datasets. Then, for each signature, the K most highly associated rules are chosen from the ranked list of all possible gene pair rules. The optimal number of rules K is chosen using 10-fold cross-validation for each model. The final K rules are combined into a single probabilistic model via a naive Bayes' classifier (6)Following the generation of models corresponding to absolute versions that mimic the ROI 95 in step (5), we remove those models that do not achieve sufficiently high agreement with the gold standard ROI 95 (Bonferroni-adjusted P value <0.05 on kappa statistics, Fig. 3b). Given the fact that many or most of the gene signatures in our database were constructed in contexts outside of breast cancer or cancer in general, there is no reason to expect all gene signatures to be informative and to induce good absolute models. After this filtering, we have 1733 absolute models that well mimic the ROI 95 gold standard in the training set. This set constitutes the final version of AIPS. (7)The final list of 1733 models (AIPS) was applied to the validation set and their agreement compared to the agreement in the training set ( Fig. 3b).

Kallisto and RSEM on single RNA-seq profiles
We

Statistical analyses
All the statistical tests performed in this study are twosided. For the log-rank test and survival analyses we used the package "survival" in R version 2.15. For analysis of Cohen's kappa statistics and the significance we used the implementation available in the "fmsb" R package. All the visualizations were performed using custom R scripts. We ran the survival analyses for the 1733 AIPS models on different cohorts of patients defined using either the molecularly defined subtype or the clinically defined subtypes. We also ran association analyses to study the association between the clinical and molecular subtypes and the low, independent, and high partitioning obtained from AIPS. For the proliferation score we used the Kruskal-Wallis test to evaluate the association with AIPS partitioning. We combined the P values obtained from both the survival analyses and and analyses of association using the Benjamini-Hochberg method to obtain a final false discovery rate (FDR). We considered an event to be significant with FDR <5%.

Current single sample pathway assignment tools are unstable
The currently available tools including PLAGE [1], z-score [2], ssGSEA [3], and GSVA [5] infer activation scores for a given patient relative to a large cohort of patient profiles. We asked if, and to what degree, scores determined in this relativistic manner, were affected by perturbations to the composition of patients in the comparative dataset. Towards this end, the activation of gene signatures from MSigDB C2 (literature curated) [12] were measured using the METABRIC dataset [13] in three distinct ways: using (1) patients with only ER-positive samples, (2) patients with only ER-negative samples, and (3) all patients. Figure 1a depicts the distribution of the absolute differences between the activation scores obtained by each tool using only ER-positive samples versus the scores obtained for all patients from METABRIC. This is repeated to compare ER-negative-specific scores versus all patients from METABRIC.
We observed a significant difference in the inferred activation score for all of the tested approaches (Wilcoxon's P value < 0.0001 for all), with the largest differences observed for the ER-positive versus pan-METABRIC comparisons (Fig. 1a). Although the scores obtained from the four tools are incomparable, the results show that scores are influenced by perturbations in the patient composition of the dataset, establishing that none of these current approaches are absolute [9]. Analogous results were obtained when the dataset was stratified by HER2 status and grade (Additional file 2: Figure S1), suggesting that many clinicopathological factors may influence pathway activation tools in this manner. Concretely, if those approaches were to be used in a clinical context, conclusions about the activation status of any given pathway would be greatly influenced by the comparative "control" group used by the treating clinician.
A simple method to measure signature activation The instability described above is analogous to what we and others have previously reported in the context of bioinformatics tools such as prediction analysis of microarray 50 (PAM50) that infer breast cancer subtype from gene expression profiles [9,10]. In our previous work (AIMS), we developed a so-called absolute method that used only a gene expression profile from the target patient without recourse to or need of a larger compendium of patient profiles for comparative analysis [9]. Our  6) we selected the final list of models that constitute AIPS by requiring significant agreement with the ROI 95 assignments in all the individual datasets present in the training set; and (7) we validated the final list of 1733 AIPS models in the validation set. b Distributions of the kappa statistics for the selected 1733 models forming AIPS (green) and the 1739 models not forming AIPS (gray) in the entire training set (using the median of the individual training sets), the individual training sets, and the validation set. c Heatmaps depict the percentage of samples of a given class obtained from the ROI 95 (e.g. low, independent (ind.), or high) assigned to another class by AIPS in the training and validation sets. d Number of genes utilized in the ROI 95 versus the AIPS models. e The ROI 95 example for an epidermal growth factor receptor (EGFR) signature from MSigDB in the McGill validation dataset. AIPS assignments are presented at the top of the heatmap. f Heatmap ordered using the Euclidean distance and the Ward's linkage method presenting the different rules utilized in the AIPS-EGFR models (red means the rule is true and white means the rule is false). Underlined genes in rules marked by a star are enriched in genes upregulated by EGFR in MCF7 cell lines [22]. g Confusion matrix representing the agreement between the single sample AIPS-EGFR model and the whole-cohort ROI 95 assignments. h Confusion matrix representing the agreement between the AIPS-EGFR assignments performed on the same RNA extraction but different platforms (RNA sequencing (RNA-seq) versus microarray). i Boxplots depicting the distribution of the percentage of agreement for the AIPS partitions done on The Cancer Genome Atlas (TCGA) samples profiled on both microarray and RNA-seq. We also present a background distribution generated from shuffling the labels 100 times. ER estrogen receptor, HER2 human epidermal growth factor receptor 2, PAM50 prediction analysis of microarray 50 goal here is to generalize and improve this framework to allow such absolute assignments to be made for any given molecular pathway or process for which there is a suitable gene signature.
The construction of our tool requires several components: (1) a collection of gene signatures that cover relevant biological processes including the hallmark molecular pathways/processes of breast cancer; (2) a suitably large collection of gene expression profiles of clinical breast cancer samples (ideally generated via both microarray and RNAseq technologies) partitioned into training and validation datasets; and (3) a "gold standard" set of positive and negative patient profiles for each pathway with known activation status.
With respect to (1), we curated a collection of 6466 gene signatures from various databases of gene signatures including MSigDB [11] and GeneSigDB [14], which we complemented with our own in-house work to curate gene signatures especially applicable to breast cancer (n = 188, see "Methods" and Fig. 2a). With respect to (2), our training dataset comprised 4510 gene expression profiles obtained from 10 different cohorts on six different platforms (Table 1). We used our previously published dataset for validation (Affymetrix Gene ST, n = 429 patients [12]) and the most recent RNA-seq data from The Cancer Genome Atlas (TCGA) project (Table 1). We next describe how we estimate our gold standard.

Estimation of the activation levels of each gene signature
Our goal is to train our AIPS algorithm to accurately infer the activation status of a given pathway within an expression profile of a patient. In order to train our algorithm, we require examples of patients that have activation statuses that we believe to be correct for each gene signature of interest. Furthermore, the examples must cover all possible states (e.g. high, low, and latent activation). The nature of human clinical samples, however, does not allow us to determine the activation status of a pathway in a direct, rigorous manner. Therefore our gold standard learning set must be comprised of estimations of statuses across the relevant expression datasets (item (2) as mentioned previously).
For each biological process of interest (item (1) as mentioned previously), we applied a de novo non-parametric rank-based method that partitions the patients in our dataset into three classes depending on the pattern of expression exhibited by the genes within the signature. The three classes correspond to those patients that appear to have high activation of the signature, low activation of the signature, and a set of patients where the expression of the genes within the signature lose their characteristic pattern of pairwise correlation (Fig. 1b provides an  example). The latter class is assigned to patients where the corresponding gene expression patterns are pairwise independent, thus supporting neither high nor low activation of the underlying pathway.
This de novo non-parametric test, referred to as the ROI at quantile q (ROI q ), proceeds as follows. In a univariate fashion, each gene within a given signature is used to rank all patients from the lowest to the highest expression. In some cases, the direction of expression of each gene within the signature is known a priori (e.g. the gene is overexpressed or underexpressed in samples with activation of the target pathway). Before ranking, we first negate any expression measurements for genes that are known to be underexpressed: such genes that are negatively correlated with activation of the signature, order the patients in the reverse order. Now for each patient, the ranks of all k genes from the signature are summed (see Additional file 1 for full details). The patients are then linearly ordered from the lowest to the highest rank. The approach of mapping expression data to a linear order, which has been used previously in breast cancer [15], makes intuitive sense as we can view the expression of each molecular process or pathway as having a state between "turned off" and "turned on" completely. Figure 1b depicts a proof-ofconcept linear ordering for an estrogen response signature from Doane et al. [16] using the METABRIC dataset (Table 1). Broadly speaking, such linear orders highlight patients at the left hand side that have low or negative expression of the signature, patients at the right hand side that have high or positive expression of the signature, and a region in the middle corresponding to patients with gene expression patterns that are independent. We refer to this as the observed linear order.
The second step in the ROI q procedure identifies the left and right boundaries of the low and high regions within the observed linear order. This is done via a permutation test where an "artificial" patient "n + 1" is created. Each of the k genes in the signature rank patient n + 1 with a uniformly randomly chosen number from (0… n + 1). Summing the randomized rank over all k genes in the signature, the position of patient n + 1 is computed within the observed linear order. This is repeated a suitably large number of times (e.g. n = 10,000). The ROI q is defined as the region that contains q% of the randomly generated samples (Fig. 1b bottom and see Additional file 1).
As expected, the patient ordering at ROI 95 for the estrogen response signature depicted in Fig. 1b is strongly associated with breast cancer subtype as defined by ER and HER2 status. In particular, the low activation region of the ordering (left) is enriched for ER-negative and/or HER2-positive tumors (Fisher's exact test, P < 0.000001, Fig. 1b, c), whereas the high activation region corresponds almost exclusively to ER-positive/HER2-negative tumors (Fisher's exact test, P < 0.000001, Fig. 1b, c). Given that ER-positive status is strongly associated with good outcome in breast cancer [12,17,18], the patient partition produced by the ROI 95 is strongly prognostic (log-rank P < 0.0001, Fig. 1d). Although only a single proof of concept example, the results suggest that the ROI q approach is capable of assigning pathway activation in breast cancer expression datasets. A more thorough investigation of the ROI 95 is presented in Additional file 1. The analyses suggest that the ROI 95 approach can faithfully recapitulate the low, independent, and high partitions of patients over a large range of biologically plausible parameters. For example, using simulated data, we tested the impact of several parameters including, for example, the gene signature size, fraction of patients in each category, and the strength of the signal on the capacity of the ROI 95 approach to correctly assign patients in the low, independent, and high categories. We confirmed using this simulated set of data that the ROI 95 is a robust approach within a wide range of parameters (see Additional file 1 and Additional file 2: Figure S2 and S3).

Identification of informative and non-informative gene signatures
To better ensure that the ROI 95 is accurately determining pathway status, we applied the method to all gene signatures in our collection (n = 6466, Fig. 2a) using the METABRIC dataset (n = 1992). The fraction of low, independent, and high samples across all signatures in our collection is presented in Fig. 2b. As a control, the gene labels of the METABRIC dataset were randomly permuted. This procedure should, with high probability, break the vast majority of gene-gene correlations within signatures, causing the fraction of uninformative genes to rise. We should then observe an increase in the independent partition of the ROI 95 with a concomitant decrease in the size of the low and high partitions. The results depicted in Fig. 2b confirm this, and suggest that for a large proportion of the signatures, the ROI 95 method is indeed assigning activation status in a very non-random fashion.
The results depicted in Fig. 2b also suggest that the ROI 95 method assigned almost every patient to the intermediate partition for some signatures. In other words, the ROI 95 method applied to these specific gene signatures was not distinguishable from random expression patterns. We removed all such gene signatures from further consideration, in particular, a gene signature was removed when the fraction of samples in the ROI 95 region exceed 0.8, as this is no better than partitions generated by random sampling. This led to a list of 3472 signatures that we considered informative in the context of breast cancer. A cutoff on the ROI 95 region will exclude gene signatures activated or repressed in less than 20% of samples. Although more liberal thresholds could be used when studying an individual gene signature, we chose this conservative threshold here to enable our high-throughput global analyses.
Given that our gene signatures were collected from various sources, we asked whether any particular source was enriched for uninformative signatures. Of the remaining informative signatures, pathway databases such as BioCarta, Kyoto Encyclopedia of genes and genomes (KEGG) and Gene Ontology (GO) have higher fractions of signatures that have near random behavior (bottom of Fig. 2c). We note that sources that contributing signatures from transcriptional profiling have a higher proportion of non-random signatures (top of Fig. 2c).

Absolute single-sample gene signature activation in breast cancer
Based on the aforementioned results, we used the ROI 95 method with the 3472 informative signatures to the training and validation datasets for calling signature activation levels using only the expression profile of a given single patient. The approach used here broadly follows our AIMS method that infers breast cancer subtype ( Fig. 3a; also [9]).
First, the ROI 95 is applied to each informative signature across 10 expression datasets generated from several microarrays (one-color and two-color) and RNA-seq platforms totaling 4510 samples (Table 1). This large and diverse training dataset provides us with more confidence that biases for specific clinicopathological or other patient variables are ablated, or at least reduced [9]. Our learning set consists of activation statuses (low, independent, and high) for 3472 signatures in each of 4510 patients from the training set. Now, for each signature a model is learnt as follows. First, we consider every possible pair of genes x and y in the expression dataset, and we ask how strongly the rule "if expression of x is greater than y, then classify signature as high" is associated with our gold standard learning set. This is repeated for a modification of the rule where x is less than y, and for the other activation statuses of independent and low. Then, for each signature, the K most highly associated rules are chosen from the ranked list of all possible gene pair rules. The K rules are combined into a single probabilistic model using a naive Bayes' classifier, and validated on an independent dataset (n = 742 samples, Table 1) [12].
The last step of our approach consists in selecting only those models with strong agreement with the ROI 95 approach using a cutoff of 0.05 on the Bonferroni adjusted Kappa's statistics P value. The two filtering steps that consist of first filtering out non-informative gene sets and then keeping models with significant agreement are essential to provide a set of reliable models (Fig. 3a).

Approximately 50% of informative signatures are amenable to absolute assignment
To ensure that our models are applicable across different technologies, we only retained models that significantly agreed with the gold standard in all 10 of the training datasets (kappa statistics Bonferroni-adjusted P value <0.05). This resulted in the retention of 1733 models (1733/3472, approximately 50%). We observed that the retained models had better agreement with the gold standard in the validation dataset in comparison with the models that were removed (Fig. 3b). This observation suggests that our training procedure did not introduce any significant over-fitting as selected models behave similarly in the training and validation sets and also all the models obtain a significant kappa statistic in the validation set (kappa statistics P values <0.01 for all).
We also stress here that the validation dataset contains data generated on a microarray platform (Affymetrix Gene ST) not present in the training dataset, suggesting that AIPS assignments are applicable on technologies not utilized in the training procedure. AIPS correctly assigns activation status for samples either assigned low or high activation (mostly red (approximately 80%) for the low-low and high-high lines in Fig. 3c). About 60% of samples were assigned to the independent class over all the 1773 models, For the samples in the gold standard that were assigned independent status, AIPS correctly assigned this status in 60% of the cases, suggesting that predictions made for samples in the independent class are generally less reliable than predictions made for the low or high classes (mostly yellow for the "ind.-ind." tagged line in Fig. 3c). Importantly, we rarely observed cases were AIPS predicted high activation when the gold standard was low, and vice versa (mostly blue for the low-high and high lines in Fig. 3c).
Last, we note that the AIPS models used fewer genes to infer pathway activation status than the original ROI 95 method to generate the gold standard (median 50 versus 200 genes for AIPS versus ROI 95 respectively; Wilcoxon's test P < 0.0001, Fig. 3d).
Overall, these analyses confirmed that AIPS could accurately recapitulate the assignments of the gold standard. The 1733 AIPS models are listed in Additional file 3: Table S1 and pathway activation assignments can be computed for new individual samples using our AIPS R package [19].
Absolute assignment of EGFR pathway activation using AIPS Epidermal growth factor receptor (EGFR) is well-studied in breast cancer with high activation of this pathway associated with poor patient outcome [20,21]. We examined the behavior of our AIPS-EGFR model in the McGill validation dataset (Fig. 3e-g). We observed that the activation of samples at the far left and right (low and high respectively) are nearly perfectly inferred by AIPS (kappa = 0.54, P < 0.0001) with the majority of disagreements related to samples in the independent region ( Fig. 3e-g). Figure 3f depicts the simple binary rules used by the AIPS model for the EGFR signature across the patient samples. There is a large cluster of EGFR-high patients associated with the PAM50 basallike (BasalL) subtype, and a second large cluster of EGFR-low patients associated with luminal A and luminal B subtypes, a finding consistent with previous studies [20,21].
Interestingly, gene set enrichment analysis of the genes selected to participate in the binary rules revealed an enrichment for genes upregulated by EGFR in MCF7 cells (FDR q value = 1.45e-17) [22]. Furthermore, all of these genes are on the right side (or left side) of binary rules associated (or not-associated) with high EGFR activation (rules marked by an asterisk in the heatmap of Fig. 3f ). Although AIPS selects gene pair rules for each model from the large space of all possible gene pairs, it still surprisingly often selects genes that were present in the original signature, and therefore are likely good markers of the underlying biological processes. The enrichment of genes from the original signature was also reported for other "absolute" models such as AIMS [9]. Almost all of our AIPS models had such enrichment (1335 out of 1733 models (77%), Additional file 2: Figure S5). It is important to note that although most models statistically significantly overlapped the original signature, the number of genes from the original signature was still below 10%, suggesting that AIPS models do require many other genes to mimic the ROI95 assignments.
We asked if the absolute nature of the AIPS method would result in a more consistent EGFR model across gene expression platforms. In particular, we asked if our AIPS model inferred the same activation status for the EGFR pathway in both the microarray and RNAseq platform for the same patient. Using TCGA data for 398 patients [23], AIPS assignments agreed on 87% of patients between both platforms (Fig. 3h, kappa = 0.81, P < 0.0001). Systematic analysis over the entire partitions (n = 1733 models) revealed that this agreement value is representative of almost all the partitions induced by AIPS and is significantly different from a random distribution (Fig. 3i, Wilcoxon's test P < 0.0001, all kappa statistics P < 0.0001) supporting the argument that absolute assignments are robust across multiple platforms [9]. Together these results suggest that AIPS is capable of inferring signature activation levels with comparable performance to relativistic tools but with the added benefits of an absolute single sample approach.

AIPS assignments agree with whole-cohort inferred pathway scores
Our goal was to compare AIPS assignments with a second approach from the literature that takes full advantage of an entire dataset to assign signature activation scores. In particular, we used 21 non-redundant scores from the publication of Gatza et al. generated from breast cancer RNA-seq expression data from the TCGA project (n = 456) [24]. Concomitantly, we estimated activation status using our AIPS models for these signatures on the same patients. Overall good agreement between AIPS and the pathway scores from Gatza et al. was observed (Fig. 4a) although the two approaches are quite dissimilar. Figure 4a suggests that well-known breast cancer biological processes are recapitulated by AIPS assignments. For example, patients with the luminal A or B subtype (LumA or LumB) are mostly assigned to the AIPS-high class for the ER gene signature, consistent with the fact that subtypes are enriched form ER-positive patients [9,25]. Also, the AIPS assignments are in good agreement with the proliferation, ESC human and MYC DUKE pathway scores as these processes are known to be associated with the highly proliferative basal-like (BasalL), Her2 and LumB subtypes [25,26]. We also observed a significant proportion of Her2-positive patients assigned to the AIPS-high Her2 gene signature. The interferon alpha and gamma, STAT1 and TNF alpha pathway scores are in good agreement with the AIPS assignments; these processes are associated with the BasalL subtype [27]. The P53 WT signature from AIPS is in good agreement with the pathway scores and is enriched for the LumA subtype that has been shown to be depleted of P53 mutations [23].
Generally, if AIPS modules and the pathway scores of Gatza et al. are in good agreement, then the patients within the high class of AIPS should also have the highest pathway scores according to Gatza et al. We tested the agreement between these two approaches and observed a strong relationship (Wilcoxon's P < 1.4e-14 for all, Additional file 2: Figure S4A). Overall this analysis suggests that the "single-sample" AIPS approach is in good agreement with an approach that uses an entire cohort of samples to judge activation.
Sample partitions induced by AIPS are prognostic and associated with breast cancer subtypes As there are many pathways and process that are known to vary in their expression across breast cancer subtypes, we investigated the relationship between patient subtype (called using the AIMS tool [9] or using clinical information) and the entire patient partitions generated by AIPS on the McGill dataset (Table 1).
We first studied the relationship between the partitions induced by AIPS and survival for the different molecular and clinical subtypes (Fig. 5a, b). We noticed that almost half (42.4%) of the partitions are significantly associated with survival if the analysis is performed on the entire cohort (Fig. 5a, b; Additional file 3: Table S1). This number drops drastically if we restrict the survival analysis to patients of given subtype. For example, only around 5% of the partitions are significantly associated with survival for the BasalL, Her2 and Luminal A and to close to nothing for the LumB and normal-like (NormalL) subtypes. Similarly, for the clinical subtypes, we found between 30 and 50% of partitions associated with the ER-positve, Her2-negative and ER-positive/Her2-negative subtypes. Those numbers drops between 3 and 6% for the ER-negative, Her2-positive and ER-negative/ Her2-negative subtypes. We found almost no partitions A Fig. 4 Assignment by absolute inference of patient signatures (AIPS) agrees with whole-cohort pathway scoring tools. a Comparison between the pathway scores retrieved from Gatza et al. and AIPS. The ordering of the rows and columns is preserved between the two heatmaps. BasalL, basal-like, HER2 human epidermal growth factor receptor 2, LumA luminal A, LumB luminal B, ER estrogen receptor, TGFB transforming growth factor beta, EGFR epidermal growth factor receptor, VEGF vascular endothelial growth factor, IFN interferon  Table S1).
We also studied the association between AIPS partitions and the molecular, ER, Her2 and proliferation score [25] using Fisher's exact test (Fig. 5a, b). We found that between 85 and 92% of the partitions are associated with the different grouping with the exception of the Her2 subtype and the clinical features (Fig. 5a, b). Almost all (92.4%) partitions are associated with proliferation.
We also examined the frequency of patients classified as low over the entire set of 1733 AIPS models (Fig. 5c). Patients with the LumA and NormalL subtyptes obtained a significantly higher number of assignments compared to the those with the remaining subtypes, with an increase of between 200 and 250 in low assignments (Fig. 5c).
AIPS applied on single RNA-seq samples can be performed in a timely fashion AIPS should enable a significant amount of information to be extracted from single RNA-seq samples and the growing number of single-cell sequencing datasets. We measured the time necessary to obtain AIPS partitions from a single RNA-seq profile using Kallisto, a fast pseudoalignment program used to obtain transcript quantification from sequencing data [28].
Using 100-bp paired-end sequencing data (median of 155 million paired-ends per patient [28], n = 96, Fig. 5d) we monitored the time taken to obtain transcript quantification and AIPS partitions using a single central processing unit (CPU). Overall, it required a median 21.6 minutes to obtain AIPS partitions for an individual patient directly from raw paired-end sequencing data (Fig. 5e).
Given that different transcript quantification pipelines return slightly different results, we evaluated the agreement between AIPS partitions made using Kallisto quantification to partitions made using an alternative approach (Bowtie2 [29] + RSEM [30]). Overall, we found significant agreement between AIPS partitions made from the two quantification approaches (Wilcoxon's test P < 0.0001 and all kappa statistics P < 0.0001) with a median agreement of approximately 85% (Fig. 5f). Together this establishes that AIPS could be applied in a time-effective manner on single-sample RNA-seq data with the aid of a sufficiently fast pseudoalignment program e.g. Kallisto.

Discussion
The work presented here is predicated upon the observation that existing pathway analysis tools are relativistic in nature. In essence, the tools make use of a large panel of samples to "judge" whether there is evidence that the given pathway is underexpressed or overexpressed relative to the panel. We showed here that the scores returned by these tools are sensitive to the composition of patients within the dataset, using a large breast cancer gene expression dataset. More precisely, we showed that the scores returned by these tools vary when the frequency of specific clinicopathological variables is perturbed. Although we have shown this is true for grade, ER and HER2 status, it is likely that many other variables that were assessed (e.g. age, claudinlow status, tumor size) and or not assessed (e.g. tumoral heterogeneity, clonal complexity, lifestyle and information on exposure of the patient) can also affect the estimations of pathway activity using these tools [31][32][33]. This is nonunintuitive: if the gene expression profile of a patient is included in two different datasets with marked differences in the overall composition of patients (e.g. they differ on the fraction of ER-positive patients), in one dataset a target pathway could be assessed as having high activation but in the second dataset as low. As these, or similar variables, are involved in most, if not all cancers, it is highly likely that this degenerate behavior is not restricted to studies of breast carcinoma. Clearly an "absolute" tool that ablates this "relativistic" behavior would be a step in the right direction.
The main focus of this study is the development of a de novo framework to estimate the activity of a given pathway using only a single sample. Here we have trained and validated predictive models for 1733 gene signatures for these pathways using a large compendium of breast cancer gene expression profiles. The profiles originate from several distinct microarray and RNA-seq platforms. In order to develop a suitable training and validation dataset, we developed the notion of the q% region of independence (ROI q ), which assigns simple discrete levels of activation for a given gene signature and a sufficiently large dataset. Using a large number of synthetic datasets we have shown that the method is robust and can faithfully retrieve low and high activation for many gene signatures within realistic configurations (see Additional file 1). By extending our previous AIMS methodology [9] with this gold standard, we were able to generate 3472 absolute single-sample gene signature activation models of which approximately 50% (1733/3472) performed sufficiently well as to be included in AIPS. We have shown that the AIPS models are more compact, their assignments are reproducible when the same patient is profiled using two distinct platforms, and the models are highly prognostic. Moreover we showed that our ability to estimate the activation of most pathways is not reduced when switching from a relativistic to an absolute method.
AIPS provides 1733 models that are immediately applicable to new breast cancer samples even when they are profiled in isolation, outside of a large cohort to make comparative assessments. Furthermore, we have shown that AIPS models are prognostic and compare favorably to other whole-cohort approaches, and that AIPS could be applied effectively to RNA-seq data. The term "absolute" expresses the idea that pathway assessment made on such a new "isolated" sample is a function only of the learning phase for each of the 1733 models, and is not done relative to a comparative cohort. This is, to the best of our knowledge, a marked difference from all other pathway tools currently available. The power of this approach is that it allows us to completely define and control the learning set, removing biases and potential confounding variables in downstream analysis. This is not possible with other current approaches where pathway analyses are affected by the other patient profiles in the cohort. Of course, here "absolute" does not imply that pathway assessments are perfect nor can the method judge in all cases the state of a pathway in a tumor relative to a healthy normal control. However, the presence of normal-like samples in our training sets allows us to assess such an "absolute state" of the pathway (the AIPS R package includes this analysis). Of course, the quality and definition of the AIPS models is still a function of the training set. In other words, modifications to the training set might impact model definitions, but here we trained our models over multiple large cohorts in order to minimize the risk that small specialized datasets would skew model parameters. Further refinement and curation of the learning dataset might potentially lead to absolute models.
Although there is a steadily increasing number of breast cancer gene signatures derived from microarray and RNA-seq based studies, at best a dozen of these signatures are currently directly available to clinicians and patients [34,35] and almost all of them are suspected to be essentially sophisticated multigene predictors of proliferation [12,26,36,37]. Given an expression profile of a patient sample, AIPS is able to estimate the activation of 1733 of the pathways, molecular processes and functions simultaneously in a timely fashion. This represents a step towards a clinically feasible tool that would provide healthcare providers and clinicians with important information on many aspects of the tumor beyond proliferation.


extfieldsofstudy: Biology Medicine

text:
UCHL1-HIF-1 axis-mediated antioxidant property of cancer cells as a therapeutic target for radiosensitization

Hypoxia-inducible factor 1 (HIF-1) has been recognized as an important mediator of the reprogramming of carbohydrate metabolic pathways from oxidative phosphorylation to accelerated glycolysis. Although this reprogramming has been associated with the antioxidant and radioresistant properties of cancer cells, gene networks triggering the HIF-1-mediated reprogramming and molecular mechanisms linking the reprogramming with radioresistance remain to be determined. Here, we show that Ubiquitin C-terminal hydrolase-L1 (UCHL1), which we previously identified as a novel HIF-1 activator, increased the radioresistance of cancer cells by producing an antioxidant, reduced glutathione (GSH), through HIF-1-mediated metabolic reprogramming. A luciferase assay to monitor HIF-1 activity demonstrated that the overexpression of UCHL1, but not its deubiquitination activity-deficient mutant (UCHL1 C90S), upregulated HIF-1 activity by stabilizing the regulatory subunit of HIF-1 (HIF-1α) in a murine breast cancer cell line, EMT6. UCHL1 overexpression induced the reprogramming of carbohydrate metabolism and increased NADPH levels in a pentose phosphate pathway (PPP)-dependent manner. The UCHL1-mediated reprogramming elevated intracellular GSH levels, and consequently induced a radioresistant phenotype in a HIF-1-dependent manner. The pharmacological inhibition of PPP canceled the UCHL1-mediated radioresistance. These results collectively suggest that cancer cells acquire antioxidant and radioresistant phenotypes through UCHL1-HIF-1-mediated metabolic reprogramming including the activation of PPP and provide a rational basis for targeting this gene network for radiosensitization.

rather than mitochondrial oxidative phosphorylation, not only under hypoxic but also normoxic conditions 6,14,15 . Glucose-6-phosphate, an intermediate metabolite of glycolysis, is the initial substrate of the pentose phosphate pathway (also known as the phosphogluconate pathway and hexose monophosphate shunt), which generates NADPH and pentoses (5-carbon sugars) as well as ribose-5-phosphate [16][17][18] . A recent study demonstrated that the pentose phosphate pathway is associated with the radioresistance of cells 19 because its byproduct, NADPH, is essential for the production of an antioxidant, reduced glutathione (GSH), from glutathione-S-S-glutathione (GSSG), and because ribose-5-phosphate is used in the de-novo synthesis of nucleotides, which are essential for repairing DNA damage. However, a gene network triggering the reprogramming of carbohydrate metabolism and the subsequent pentose phosphate pathway has yet to be fully elucidated.
In the present study, we focused on the UCHL1-HIF-1 axis and investigated whether it induced the metabolic reprogramming, antioxidant property, and radioresistant phenotype of cancer cells using murine breast cancer-derived EMT6 cells.

Results
UCHL1 deubiquitinates HIF-1α protein and upregulates HIF-1 activity in murine breast cancer-derived EMT6 cells. Because HIF-1 activity is known to be regulated at multiple steps, we first aimed to identify the key regulatory mechanism in the UCHL1-mediated activation of HIF-1 in EMT6 cells. First, we performed a luciferase assay using the 5HRE-luc reporter gene 43 , which expresses luciferase bioluminescence in a HIF-1-dependent manner, in order to test whether UCHL1 enhanced HIF-1 activity in EMT6 cells. We transfected the cells with the reporter gene and UCHL1 expression vector, cultured them under normoxic or hypoxic conditions, and performed the luciferase assay. The forced expression of UCHL1 significantly enhanced HIF-1 activity regardless of the oxygen conditions (Fig. 1a). We next examined the impact of UCHL1 overexpression on the translational initiation of HIF-1α protein. We employed the HIF-1α 5′UTR-luc reporter gene 39 , in which a pyrimidine tract of HIF-1α 5′UTR was inserted between the constitutively active CMV promoter and luciferase coding sequence 27,39 , and confirmed that UCHL1 exhibited no influence on the translational initiation of the HIF-1α gene (Fig. 1b).
Because HIF-1 activity is markedly influenced by the balance between the degradation and stabilization of HIF-1α protein, we next analyzed the possibility that UCHL1 affected this balance. Western blotting using an antibody to detect the hydroxylated proline residue 564 of HIF-1α, demonstrated that the hydroxylation, which could be detected in the presence of the proteasome inhibitor MG132 in an oxygen-and PHD-dependent manner, was not affected by the forced expression of UCHL1 (Fig. 1c). On the other hand, by utilizing the SV40p-ODD-Luc reporter gene 44 , which expressed a fusion protein composed of a HIF-1α ODD domain (HIF-1α 548-604) and luciferase from the constitutively active SV40 promoter, we confirmed that UCHL1 was able to actively stabilize the ODD-luciferase fusion protein (Fig. 1d). This finding was supported by Western blotting for HIF-1α protein; namely, expression levels of HIF-1α protein were significantly increased under hypoxic conditions when the UCHL1 expression vector was introduced into the cells ( Fig. 1e and f). However, whether this also occurred under normoxic conditions as well remained unclear, as the basal HIF-1α expression levels were below the detection limits of Western blotting ( Fig. 1e and f). Since UCHL1 has been recognized as a deubiquitinating enzyme 45 , we tested whether its ubiquitinating activity is critical for the stabilization of HIF-1α protein. The luciferase assay using the 5HRE-luc reporter gene showed that a catalytically inactive mutant of UCHL1, UCHL1 C90S 46 , failed to upregulate HIF-1 activity (Fig. 1g). The luciferase assay using the SV40p-ODD-Luc reporter gene also demonstrated that UCHL1 C90S overexpression did not stabilize the ODD-fusion protein (Fig. 1h). In addition, when we tested whether UCHL1 influenced the ubiquitination status of HIF-1α by performing an immunoprecipitation experiment, the forced expression of UCHL1 markedly decreased the amount of ubiquitinated HIF-1α (Fig. 1i). Finally, although we examined the function of UCHL1 in the upregulation of the transactivation activity of HIF-1 by utilizing another luciferase assay system 40 , the forced expression of UCHL1 had no effect on it (Fig. 1j). These results collectively indicate that UCHL1 stabilizes HIF-1α protein, and elicits HIF-1 activity in a deubiquitinating activity-dependent manner in EMT6 cells. The UCHL1-HIF-1 axis induces reprograming of the glucose metabolic pathway and subsequent production of an antioxidant, GSH. Although HIF-1 has been reported to induce reprogramming of the glucose metabolic pathway from mitochondrial oxidative phosphorylation to glycolysis 23,24 , it remains unclear whether UCHL1 acts as a molecular trigger for this switch. To examine this possibility, we analyzed the influence of UCHL1 overexpression on the choice of metabolic pathway by quantifying levels of an end-metabolite of glycolysis, lactate, and primary metabolites in the TCA cycle, citrate and isocitrate (Fig. 2). LC/ MS-based metabolite analyses demonstrated that the overexpression of UCHL1 significantly increases the flux of metabolism from glucose to lactate (Fig. 2a) and, on the other hand, it decreased glucose metabolism to both citrate and isocitrate ( Fig. 2b and c). The observed changes in flux were partially but significantly suppressed by silencing the expression of HIF-1α. Importantly, the suppressive impacts of HIF-1α silencing were detected in the presence of UCHL1 expression, but not in its absence (Fig. 2a-c). All of the data suggest an important role of UCHL1 as a trigger for the HIF-1-dependent metabolic reprogramming from mitochondrial oxidative phosphorylation to aerobic glycolysis ( Fig. 2a-c).
We next tested whether reprogramming was accompanied by activation of the pentose phosphate pathway and whether this led to the antioxidant properties of cancer cells. Quantitative analyses of the ratio of NADPH to NADP + (NADPH/NADP + ) and that of GSH to GSSG (GSH/GSSG) revealed that the aberrant overexpression of UCHL1 increased intracellular levels of both NADPH and GSH ( Fig. 3a and b). The increases in both NADPH and GSH were partially but significantly suppressed by silencing the expression of a key element of the Figure 1. Importance of the deubiquitinating activity of UCHL1 in the upregulation of HIF-1 activity in murine breast cancer EMT6 cells. (a-h) EMT6 cells were transfected with either pcDNA4/UCHL1 (UCHL1) or its empty vector (EV), cultured under the indicated oxygen conditions, and subjected to both the luciferase assay using a co-transfected reporter plasmid, p5HRE-luc (a,g), pGL3/HIF-1α-5′UTR-Luc (b), or pGL3/ ODD-luc (d,h), and Western blotting using the indicated antibodies (a-h). The relative intensity of the bands in Fig. 1e (HIF-1α vs. β-actin) was quantified using ImageJ (f). (i) Cells transiently transfected with both pMT132 and either pcDNA4/myc-His A (EV) or pcDNA4A/HIF-1α-myc were additionally transfected with pcDNA4/ myc-His A (EV) or pcDNA4/UCHL1 (UCHL1), and cultured with MG132 (10 μM) for 6 h. Ubiquitinated HIF-1α was immunoprecipitated using an anti-myc antibody and detected with an anti-HA antibody (left). One-twentieth of the whole cell lysate was subjected to Western blotting with the indicated antibodies as input controls (right). (j) EMT6 cells were cotransfected with four kinds of plasmids: pG5H1bLuc containing five Gal4-binding sites upstream of an adenovirus E1b promoter and firefly luciferase CDS 33, 50 , pRL-SV40 as an internal control for normalization, pcDNA6/Gal4 DBD-HIF-1αTAD P564A, and either pcDNA4/UCHL1 (UCHL1) or its EV. Cells were cultured under normoxic or hypoxic conditions and subjected to the luciferase assay to evaluate the transactivation activity of HIF-1α TADs. All results (except for c and e) are presented as means ± s.d. n = 3. *P < 0.05. **P < 0.01. NS, not significant (Student's t-test).
pentose phosphate pathway, the glucose-6-phosphate dehydrogenase X-linked (G6pdx) gene ( Fig. 3a and b). In agreement with these results, luciferase assay-based quantification experiments also confirmed that overexpression of UCHL1 significantly increased the intracellular levels of both NADPH and GSH in a pentose phosphate pathway-dependent manner ( Fig. 3c and d). On the other hand, forced expression of the UCHL1 C90S mutant neither induced the carbohydrate metabolic reprogramming nor increased the levels of antioxidant GSH (Figs 2d-f and 3e and f). Moreover, the UCHL-dependent increases in the levels of both NADPH and GSH were almost completely suppressed by silencing the HIF-1α gene ( Fig. 3g and h). Taken together, these results strongly suggest the possibility that activation of the UCHL1-HIF-1 axis causes the production of the antioxidant GSH by reprograming the glucose metabolic pathway and stimulating the pentose phosphate pathway.
The UCHL1-HIF-1 axis functions in the induction of the radioresistant phenotype of cancer cells. We then performed conventional in vitro clonogenic cell survival assays to investigate whether UCHL1 causes the radioresistance of cancer cells in a HIF-1-dependent manner. EMT6 cells were transfected with the UCHL1 expression vector or its empty vector as a negative control and subjected to various doses of X-irradiation. The surviving fraction, calculated as described previously 47 , demonstrated that the UCHL1 overexpression significantly increased clonogenic survival after X-irradiation (Fig. 4a). The dose of radiation needed to decrease the number of colonies by 90% (D 10 value) was significantly increased after UCHL1 overexpression from 3.25 ± 0.08 to 5.54 ± 0.78 Gy in EMT6 cells (P < 0.05; Table 1). The UCHL1-mediated radioresistance was confirmed using another cell line, HeLa (Supplementary Figure 1). We conducted similar clonogenic survival assays after silencing HIF-1α expression in parallel with UCHL1 overexpression in order to determine whether the UCHL1-mediated radioresistance was dependent on HIF-1. We found that HIF-1 knockdown completely abrogated the radioresistance of EMT6 cells mediated by UCHL1, as expected ( Fig. 4b and Table 1). D 10 values were equivalent regardless of the expression levels of UCHL1 when HIF-1α expression was silenced; D 10 values were 3.49 ± 0.12 and 3.64 ± 0.17 Gy in the EV and UCHL1 groups, respectively (P = 0.13; Table 1). In order to investigate whether the UCHL1-mediated radioresistance was attributable to the increased antioxidant property of cells, we performed an additional clonogenic cell survival assay with X-irradiation in the presence or absence of one of the most representative antioxidants, N-acetyl cysteine (NAC). We confirmed that, although the pH of the culture medium, which was originally 7.62 ± 0.03, decreased to 7.39 ± 0.03 at 6 min after NAC addition, it immediately returned to its original value, 7.60 ± 0.01, and remainded at that level until cells were treated with radiation; therefore, the change in pH did not directly influence the radiosensitivity of cells (Supplementary Figure 2). The surviving fraction was significantly increased by 5 mM NAC treatment in cases without UCHL1 overexpression ( Fig. 4c and d). However, in cases with UCHL1 overexpression, the NAC treatment exhibited no further influence on cell survival ( Fig. 4c and d), collectively suggesting that UCHL1 has the potential to induce the maximum cellular antioxidant effect on radiosensitivity. The UCHL1-dependent increase in cellular radioresistance was markedly decreased when intracellular levels of the antioxidant GSH were decreased by a G6pdx inhibitor, 6AN (Fig. 4e). All of these results strongly suggest that the aberrant overexpression of UCHL1 induces the antioxidant and radioresistant properties of cancer cells in a HIF-1-and G6pdx-mediated PPP-dependent manner.

Discussion
In the present study, we found that the UCHL1-mediated activation of HIF-1 through the deubiquitination of HIF-1α protein induced the antioxidant and radioresistant properties of cancer cells by producing an antioxidant, GSH, through the so-called carbohydrate metabolic reprogramming and subsequent activation of the pentose phosphate pathway.
The luciferase assay using a deubiquitinating activity-deficient mutant of UCHL1 (C90S mutant) demonstrated that the ubiquitination activity of UCHL1 was essential to stabilize the ODD-fusion protein and upregulated HIF-1 activity in breast cancer-derived EMT6 cells. This result is consistent with a previous report that UCHL1 stabilized HIF-1α protein when VHL functioned as a key component of E3 ubiquitin ligase in the ubiquitination of HIF-1α protein 39 . In addition to such a molecular mechanism, we recently revealed the possibility that UCHL1 increases the expression levels of HIF-1α by upregulating the efficiency of the transcriptional initiation of the HIF-1α gene (data not shown). In order to fully elucidate the molecular mechanisms underlying the UCHL1-mediated upregulation of HIF-1 activity, further investigation is needed.
In the present study, although the constitutively active cytomegalovirus (CMV) promoter was exploited in the UCHL1 expression vector, UCHL1 protein levels were significantly increased under hypoxic conditions. These observations suggest that UCHL1 expression was upregulated at a post-transcriptional level, such as at mRNA stability levels, translational initiation levels, and/or protein stability levels. The increase in the UCHL1 protein levels under hypoxia might contribute to the rapid accumulation of HIF-1α protein in response to acute hypoxic stimuli. Alternatively, it may suggest the existence of a positive feedforward loop that boosts the accumulation of HIF-1α protein in the case that the hypoxia-dependent increase in the UCHL1 levels is HIF-1-dependent.
The radioresistance of cancer cells is influenced by various intrinsic and extrinsic factors, such as DNA damage repair activity, the cell cycle status, oxygen availability, and pH. Especially, gene networks, which induce the antioxidant property of cancer cells, have drawn considerable attention in recent years. Production of the most representative antioxidant, reduced glutathione (GSH), is mediated by multiple regulatory steps: cysteine uptake by the cystine/glutamate antiporter (system xc-) 48 , glutathione synthesis by the glutathione synthetase  ; a-h), UCHL1 (a-d,g,h), or UCHL1 C90S (e,f), or for shRNA against HIF-1α (shHIF-1α; g,h) or scramble shRNA (shNC; g,h), or with the indicated siRNA for G6pdx (siG6pdx; a-d) or scramble siRNA (siNC; a-d) were subjected to quantitative analyses of the ratio of NADPH to NADH + (a,c,e,g) and GSH to GSSG (b,d,f,h) using the NADP/NADPH Quantification Colorimetric Kit (a,e,g), the NADP/NADPH-Glo Assay Kit (c), the Glutathione Colorimetric Assay Kit (b,f,h), and the GSH/ GSSG-Glo Assay Kit (d). Means ± s.d. n = 3. *P < 0.05, **P < 0.01. NS, not significant (Student's t-test).
(GSS) 49 , and the reduction of glutathione-S-S-glutathione (GSSG) to GSH by glutathione-disulfide reductase (GSR), which uses NADPH as an electron donor 49 . Because NADPH is known to be provided as a byproduct of the pentose-phosphate pathway (PPP) 17,18 , our result that UCHL1 overexpression increased radioresistance by producing GSH through the accelerated glycolysis and PPP is reasonable.
The UCHL1-dependent increase in the intracellular GSH levels was significantly suppressed by silencing the expression of a key molecule of the pentose phosphate pathway, the glucose-6-phosphate dehydrogenase X-linked (G6pdx) gene 17,18 . Moreover, silencing the HIF-1α gene completely abrogated the UCHL1-mediated radioresistance of cancer cells. Based on these findings, our study provides an insight into a novel strategy targeting G6pdx and HIF-1α to overcome the UCHL1-dependent radioresistance of cancer cells.
Our clonogenic cell survival assays and the quantitative analysis of metabolite levels collectively demonstrated that the UCHL1-mediated radioresistance was at least in part dependent on the antioxidant property of cancer  cells elicited by HIF-1. However, whether this mechanism is fully responsible for the radioresistance is questionable because HIF-1 is known to have numerous functions that potentially influence the radiosensitivity/radioresistance of cells, such as cell cycle regulation 6,10,11 . Further investigation is needed to fully understand the downstream effectors of the UCHL1-HIF-1 axis, which play critical roles in increasing the radioresistance of cancer cells.

Luciferase assay and
were transfected with the indicated plasmids, they were harvested in 250 μL Cell Lytic Buffer (Sigma-Aldrich). The HIF-1α-myc protein was immunoprecipitated using the Immunoprecipitation Kit Dynabeads Protein G (Life Technologies) with anti-myc antibody according to the manufacturer's instructions. Western blotting was performed using anti-HA antibody.

Quantifications of metabolites and reduced glutathione. LC/MS-based metabolome analysis
to quantify the levels of 13 C 6 -labeled lactate, citrate, and isocitrate was performed as described previously 40 . Intracellular NADPH levels were quantified using the NADP/NADPH Quantification Colorimetric Kit (BioVision Inc.) and NADP/NADPH-Glo Assay Kit (Promega) according to the manufacturers' instructions. Intracellular GSH levels were quantified using the Glutathione Colorimetric Assay Kit (BioVision Inc.) and GSH/ GSSG-Glo Assay Kit (Promega) according to the manufacturers' instructions.
Clonogenic survival assay. The indicated cells (100 and 1,000 cells/60-mm dish for 0 and 2/4/8 Gy, respectively) were precultured for 24 hours with or without the indicated concentrations of NAC (4c and 4d) or 6AN (4e), treated with the indicated dose of X-radiation (Acrobio Co., Tokyo, Japan), and cultured for 2 additional weeks. NAC and 6AN were removed 24 and 1 hour after the radiation, respectively. Surviving colonies were fixed with 70% ethanol and stained with Giemsa solution. Colonies consisting of more than 50 cells were counted as surviving colonies. The plating efficiency and surviving fraction were calculated as described previously 47 .
Statistical analyses. The significance of differences was determined using Student's t-test. A P-value < 0.05 was considered to be significant.


s2fieldsofstudy:
Biology
Computer Science

text:
Tissue-Specific Mitochondrial Decoding of Cytoplasmic Ca2+ Signals Is Controlled by the Stoichiometry of MICU1/2 and MCU

SUMMARY Mitochondrial Ca2+ uptake through the Ca2+ uniporter supports cell functions, including oxidative metabolism, while meeting tissue-specific calcium signaling patterns and energy needs. The molecular mechanisms underlying tissue-specific control of the uniporter are unknown. Here, we investigated a possible role for tissue-specific stoichiometry between the Ca2+-sensing regulators (MICUs) and pore unit (MCU) of the uniporter. Low MICU1:MCU protein ratio lowered the [Ca2+] threshold for Ca2+ uptake and activation of oxidative metabolism but decreased the cooperativity of uniporter activation in heart and skeletal muscle compared to liver. In MICU1-overexpressing cells, MICU1 was pulled down by MCU proportionally to MICU1 overexpression, suggesting that MICU1:MCU protein ratio directly reflected their association. Overexpressing MICU1 in the heart increased MICU1:MCU ratio, leading to liver-like mitochondrial Ca2+ uptake phenotype and cardiac contractile dysfunction. Thus, the proportion of MICU1-free and MICU1-associated MCU controls these tissue-specific uniporter phenotypes and downstream Ca2+ tuning of oxidative metabolism.

SUPPLEMENTAL EXPERIMENTAL PROCEDURES Chemicals
Standard chemicals were purchased from Fisher Scientific or Sigma-Aldrich. Thapsigargin and CGP-37157 were from Enzo Life Sciences (Plymouth Meeting, PA). Fluorescent Ca 2+ -indicator dyes, and Pluronic F-127 were from Molecular Probes or Teflabs. Chelex 100 Sodium form was purchased from BioRad.
AAV and DNA constructs AAV9-Luciferase was from Penn Vector Core (Philadelphia, PA) and AAV9-MICU1 was custom prepared by Vector Biolabs (Philadelphia, PA). FLAG tagged-MCU was described in (Patron et al., 2014).

Cell culture and transfection
HeLa cells were grown in DMEM, 10% FBS and 1% penicillin/streptomycin and were transiently transfected with siRNA or plasmid DNA as described earlier (Csordas et al., 2013). MICU1-KO HEKs were kindly provided by Dr. Vamsi Mootha and grown as previously described . Transfection in HEKs cells was performed using Lipofectamine 3000 (Thermo Fisher).

Heart, skeletal muscle and liver mitochondria isolation
Mice were euthanized by cervical dislocation. Heart was rapidly excised and washed in cold isolation buffer containing 225mM Mannitol, 75mM Sucrose, 20mM Hepes and 0.1mM EGTA, at pH 7.4. Heart tissue was finely minced and incubated for 1min in the isolation buffer with 0.5% BSA and 3U/g of tissue of proteases (Sigma). Then, the cardiac tissue was homogenized with a glass Teflon potter homogenizer and the homogenate centrifuged at 500g/5min/4°C to remove nuclei and cellular debris. The supernatant was next centrifuged at 9,000g/10min/4°C and the mitochondrial pellet was resuspended in 225mM Mannitol, 75mM Sucrose, and 10mM Hepes at pH 7.4 and kept on ice. After liver harvesting, liver mitochondria were isolated by differential centrifugations. Briefly, liver was minced in a cold isolation buffer containing 250mM Sucrose, 10mM Tris/HCl, 0.1mM EGTA and 0.5% BSA, at pH 7.4 and homogenized in a glass Teflon potter. Cellular debris were removed by two consecutive centrifuges at 370g/10min/4°C. Supernatant was then centrifuged at 7,500g/10min/4°C and mitochondrial pellet washed once in isolation buffer without BSA before a new centrifugation at 7,500g/10min/4°C. Finally, liver mitochondria were resuspended in 250mM Sucrose and 10mM Tris/HCl at pH 7.4 and kept on ice. Skeletal muscles were carefully dissected, minced and homogenized in a glass Teflon potter in the homogenization buffer containing 140mM KCl, 20mM Hepes, 5mM MgCl 2 , 2mM EGTA, 1mM ATP(Di-K) at pH 7.0 with 1% BSA and 2U/g of tissue of proteases (Sigma). After centrifugation at 500g/10min/4°C, the collected supernatant was further centrifuged at 9,000g/8min/4°C. Pellet was resuspended in the isolation medium (140mM KCl, 20mM Hepes, 5mM MgCl 2 , 1mM EGTA at pH 7.0) and incubated on ice for 5min for myofibrillar re-polymerization. The two previous centrifuge steps were repeated and the final mitochondrial pellet resuspended in the isolation medium. Protein content was determined by the Lowry method (BioRad).

Live cell calcium imaging
RCaMP and ratiometric pericam were described in (Akerboom et al., 2013;Nagai et al., 2001). For imaging experiments, the cells were pre-incubated in a serum-free extracellular medium (ECM, 120 mM NaCl, 5 mM NaHCO 3 , 10 mM Na-HEPES, 4.7 mM KCl, 1. Fura2/furaFF and TMRM fluorescence were recorded simultaneously using 340-380 nm excitation and 500 nm emission, and 545 nm excitation and 580 nm emission, respectively. Complete dissipation of ΔΨ m was elicited by addition of 2 µM FCCP. Calibration of the fura signal was carried out at the end of each measurement, adding 1 mM CaCl 2 , followed by 10 mM EGTA/Tris, pH 8.5.

RNA isolation and qPCR
Total RNA were extracted from cells using TRIzol reagent (Ambion) and then treated with RQ1 DNase (Promega). Mice liver and heart RNA were isolated using the animal tissue RNA purification kit (Norgen) according to the manufacturer's instructions. Quantification was done on the nanophotometer Pearl (Implen).Then cDNA was synthesized using SuperScript III (Invitrogen) and used for qPCR reactions using SYBR Green (Invitrogen) on an ABI Prism 7000 sequence detection system (Life Technologies). Data were analyzed using the comparative 2-ΔΔCt method. Ct of the gene of interest was normalized to that of β-actin for heart and liver comparison. Primers used are listed below: