Nonmonotonic Pathway Gene Expression Analysis Reveals Oncogenic Role of p27/Kip1 at Intermediate Dose.

The mechanistic basis by which the level of p27Kip1 expression influences tumor aggressiveness and patient mortality remains unclear. To elucidate the competing tumor-suppressing and oncogenic effects of p27Kip1 on gene expression in tumors, we analyzed the transcriptomes of squamous cell papilloma derived from Cdkn1b nullizygous, heterozygous, and wild-type mice. We developed a novel functional pathway analysis method capable of testing directional and nonmonotonic dose response. This analysis can reveal potential causal relationships that might have been missed by other nondirectional pathway analysis methods. Applying this method to capture dose-response curves in papilloma gene expression data, we show that several known cancer pathways are dominated by low-high-low gene expression responses to increasing p27 gene doses. The oncogene cyclin D1, whose expression is elevated at an intermediate p27 dose, is the most responsive gene shared by these cancer pathways. Therefore, intermediate levels of p27 may promote cellular processes favoring tumorigenesis-strikingly consistent with the dominance of heterozygous mutations in CDKN1B seen in human cancers. Our findings shed new light on regulatory mechanisms for both pro- and anti-tumorigenic roles of p27Kip1. Functional pathway dose-response analysis provides a unique opportunity to uncover nonmonotonic patterns in biological systems.

Introduction

In many human cancers, the expression level of the cyclin-dependent kinase (CDK) inhibitor p27Kip1 is strongly associated with prognosis and patient mortality.[1] Although p27 is a CDK inhibitor, it also has cytoplasmic functions and regulates cell migration.[2,3] As death from cancer is primarily due to metastatic spread of the disease, this suggests that p27 could contribute to patient mortality by affecting invasion and metastatic dissemination. However, the mechanistic basis by which p27 influences tumor aggressiveness and patient mortality remains unclear. p27 can affect cell proliferation,[4] migration,[2,3] apoptosis,[5-7] and differentiation,[8,9] through both CDK-dependent and CDK-independent mechanisms.[10,11] However, the relative contributions of these processes to cancer outcome are unknown. In addition, p27 was the first tumor suppressor gene shown to be haploinsufficient, a genetic concept that is now widely accepted.[12,13] Haploinsufficiency in this case means the loss of only 1 allele predisposes to tumor development. This property of p27 makes the study of its phenotypic effects even more complicated. Indeed, the low expression of p27 was widely observed to correlate with poor prognosis in various cancers, including prostate,[14] breast,[15] colorectal,[16,17] and ovarian cancers.[18] However, when its interaction with the cell cycle is disabled by deleting its CDK inhibitory domain, p27 has oncogenic properties, causing multiple tumor phenotypes by modulating prometastatic factors such as cell migration and invasion.[19,20] Increased motility and invasiveness are essential processes for cancer cell metastatic dissemination.[21] However, findings of p27’s roles on cell migration and invasion have been inconsistent. High expression levels of p27 were reported to promote migration in cerebral cortex[9] and more recently in bladder cancer cells.[22] Meanwhile, low expression of p27 is associated with invasion in neuroblastoma, bone marrow,[23] and early cervical carcinoma.[24] A nonmonotonic function of p27 was first proposed for its role in cell proliferation[25] but not yet for cell motility. Furthermore, the positive relation of p27 to cell migration and its negative relation to invasion suggest an intermediate dose of p27 at which both cell migration and invasion could be active. This hypothesis motivated us to analyze gene expression in chemically induced mouse squamous cell papillomas derived from p27 null, heterozygous, and wild-type mice to capture potential nonmonotonic effects not feasible with a 2-dose experiment. We characterized the transcriptomic responses using a novel informatics method for functional pathway analysis.

Pathway analysis, studying gene activities collectively, provides statistical robustness over single-gene approaches. Many analysis tools are effective at ranking pathways by their responses to experimental conditions. These approaches can be grouped into 3 categories: overrepresentation analysis (ORA), functional class scoring (FCS), and topology based.[26] Addressing several limitations of ORA—an attempt in early days, the FCS approaches, including gene set enrichment analysis[27] and its variations, compute a pathway statistic cumulative over all gene-level statistics. However, recent topology-based pathway analysis methods integrate gene interaction information to understand pathway rewiring[28-31] but require larger sample sizes than ORA or FCS. These pathway analysis approaches, however, were designed to determine the association between pathway responses and experimental conditions, not optimal for the complex functional dependency of pathway responses to dosage. Indeed, nonmonotonic dose-response curves are widely observed in biology such as in studies of hormones and endocrine-disrupting chemicals.[32]

To overcome the limitations of current pathway analysis methods for detecting nonmonotonic dose-response relationships, we present a novel method called functional pathway dose-response analysis to test whether pathway response is a mathematical function of gene dosage, without assuming a specific underlying functional form. This method is based on the recent FunChisq method which detects nonconstant and nonparametric functional dependency among random variables for causal inference.[33] FunChisq demonstrated outstanding empirical effectiveness in the recent HPN-DREAM Breast Cancer Network Inference Challenges.[34] A gene is functionally dependent on stimulant dosage if the gene’s expression value is uniquely determined by the dose. We show that our method is mathematically optimal in detecting functional dependency of pathway response to stimulant doses. We also found that the gamma distribution approximation of the null pathway response statistic is more accurate than the χ2 distribution. Next, we applied FunChisq and functional pathway dose-response analysis to characterize how genes and pathways respond to Cdkn1b gene dosage in mouse papillomas derived from 3 genotypes: wild type (p27+/+), heterozygous (p27+/−), and nullizygous (p27−/−). Our pathway analysis revealed that a disproportionally high number of genes in several KEGG cancer pathways[35] are more expressed at the intermediate gene dosage of p27. These include the oncogene Ccnd1 (cyclin D1). This result supports the notion that intermediate levels of p27 may have pro-oncogenic properties. It coincides with additional evidence showing the dominance of CDKN1B heterozygous somatic mutations in thousands of samples from diverse tumor types in The Cancer Genome Atlas (TCGA). Although p27 was previously suggested to contribute to nonmonotonic regulation of some cell cycle genes,[25] we present evidence for the genome-wide pathway effects of p27 dosage influencing many aspects of tumorigenesis. Therefore, complementing posttranscription modification and cellular localization explanations of p27’s oncogenic roles, our findings provide new insights toward p27 function and molecular mechanisms for gene regulation in tumor development.

Materials and Methods

Dose-response gene expression data and KEGG pathways

We collected samples of DMBA/TPA-induced skin papillomas (>5 mm) from C57BL/6J NIH F2 mice of 3 p27 genotypes: (1) nullizygous mutation (p27−/−), corresponding to a low dose; (2) heterozygous mutation for a null allele (p27+/−), corresponding to an intermediate dose; and (3) wild type (p27+/+), corresponding to a high dose.[36] Four papillomas from each p27 genotype were examined. Whole-genome microarray analysis was performed for the samples using Illumina MouseRef-8 v2.0 Expression BeadChips (Illumina, San Diego, CA, USA; 25 600 probe sets). Total RNA was isolated with Trizol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA) and cleaned with RNeasy Kit (Qiagen, Valencia, CA, USA). RNA integrity and purity were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Illumina BeadChips were hybridized with 1 µg of labeled complementary RNA, washed, stained, and scanned (Illumina BeadStation Scanner) according to Illumina’s protocols. Data values with detection scores were compiled using BeadStudio (Illumina) and raw intensity data were quantile normalized by Robust Multi-Array Analysis summarization for statistical analysis.[37]

We extracted 281 Mus musculus functional pathways and disease processes from KEGG. Among them are 2 interesting overlapped subsets: 14 pathways known to be affected by p27 and 17 pathways critical in cancer. The pathways involve a total of 18 443 genes.

The functional chi-square test for single-gene dose response

The functional χ2 (FunChisq) test determines whether discrete random variable Y is a nonconstant function of discrete random variable X. Let r and s be the numbers of discrete levels of X and Y, respectively. We construct an contingency table from observed data of X and Y, where X and Y are row and column variables, respectively. Let n be the sum of row i and n the sum of column j. Let n be the sample size. The functional χ2 statistic is defined as follows[33]:

and it follows a χ2 distribution with degrees of freedom under the null hypothesis that X and Y are statistically independent with the assumption that Y is uniformly distributed. Important properties of the test statistic, such as its null distribution, optimality for functional dependency, and asymmetry with respect to X and Y, were established previously.[33]

Functional pathway dose-response analysis

Let P be the set of genes on a given pathway and P also serves as an identifier of the pathway. To assess the functional dependency of responses of all genes on pathway P on dosage X, we design a novel pathway dose-response statistic by summing up the functional χ2 statistics of member genes in the pathway:

where F is the collection of dose-response functions for each gene in the pathway. We do not assume the same dose-response curve for different genes in the pathway. This is important for pathways that include genes which respond to the stimulus either positively or negatively. Based on the optimality of the functional χ2 statistic for functional dependency,[33] we derive the following theorem to establish the optimality of the pathway dose-response statistic.

Theorem 1. (Optimality of pathway dose-response statistic.)

The pathway dose-response statistic is maximized when the responses of each gene g in a pathway P are nonconstant empirical functions of dosage X represented by function set F and is minimized to 0 when all gene responses are empirically independent of dosage.

A proof of the theorem is provided in Supplementary Note.

If the expression levels of each gene in a pathway are statistically independent, then follows a χ2 distribution under the null hypothesis of gene expression in the pathway being independent of dosage. The degrees of freedom of are the sum of those of over all g’s in P.

Because the member genes in a biological pathway are biochemically connected, their expression levels cannot be assumed independent in general. This implies that does not necessarily follow a χ2 distribution under the null hypothesis. To approximate the null distribution more accurately than the χ2 distribution, we developed a gamma approximation for the null distribution of using a bootstrap strategy. We permute samples instead of genes to maintain dependency among genes and randomly assign the permuted samples to stimulant doses, such that the experimental design is the same with the original except the swapped samples. For the 12 samples of transcriptome responding to 3 p27 doses, each permutation is a random assignment of 4 samples to each dose. The gamma distribution family of shape parameter k and scale parameter θ contains probability density functions in the form of

where the gamma function Γ(k) is defined as follows:

We first obtain the estimated mean and variance of by bootstrapping. Then, we calculate the shape and scale parameters by , . From the gamma distribution with the estimated parameters, we compute a P value to indicate the statistical significance of . In addition, even when is indeed χ2 distributed, gamma approximation is still applicable because χ2 distributions are a special case in the gamma distribution family.

Pathway pattern enrichment analysis

In each pathway, the member genes can differ greatly in responding to p27. We classify the responses into 4 patterns: increasing (INC), decreasing (DEC), high-low-high (HLH), and low-high-low (LHL). The enrichment analysis is to tell which pattern is overrepresented in a pathway, regarding the proportion of that pattern in the whole gene set. Let N be the total number of genes in the gene pool; K be the number of genes from the gene pool which respond in pattern P, n be the number of genes in the analyzed pathway, and k be the number of the genes in the pathway which respond in pattern P. Then, the probability of k in a pathway is defined by

where k follows a hypergeometric distribution if the pathway is equally enriched in pattern P as in the entire gene pool.

Results

Overview of the functional pathway dose-response analysis method

Figure 1 illustrates functional pathway dose-response analysis, using the autophagy pathway as an example. The input includes a list of genes on the pathway and gene expression data that were first discretized. The test statistic for a pathway dose-response is computed as the summation of gene-level functional χ2 statistics. Finally, the output of the analysis, the statistical significance of pathway response to stimulant dosage, is calculated based on an asymptotic gamma distribution whose parameters are estimated by bootstrapping through permutation.

Figure 1.

Overview of the functional pathway dose-response analysis. The input to the analysis includes transcriptome profiles responding to p27 dosage under 3 experimental conditions and gene sets from biological pathways. The output comprises functional pathway dose-response statistics and their statistical significance for each pathway. The input gene expression data were discretized and then used to form contingency tables with p27 dosage. The functional χ2 statistics for each table are computed and summed to produce the pathway dose-response statistic. The P value is calculated by a gamma distribution whose scale and shape parameters are estimated by bootstrapping.

Advantage of functional pathway dose-response analysis by simulation study

To validate the effectiveness of functional pathway dose-response analysis, we benchmarked it against an alternative method that replaces functional χ2 by the Pearson association χ2. The Pearson χ2 is a representative of many pathway statistics in the sense that it measures association but not mathematical functional dependency in a contingency table. We simulated pathway gene expression data modeled after 17 known cancer pathways from KEGG, at 3 levels of stimulant dosage in 4 biological replicates. The percentage of dose-responding genes in a pathway was controlled at 10%, 50%, and 90%. We also created null pathways whose genes do not functionally respond to stimulant dosage. Noise was generated using a house noise model[38] at the noise levels of 20%, 30%, and 40%. Figure 2 shows areas under the receiver operating characteristic (ROC) and precision-recall curves of the 2 pathway χ2 methods under all 9 configurations of response percentages and noise levels. The actual ROC and precision-recall curves for each configuration are shown in Supplementary Figures S1 and S2. Both methods improved when more genes on the pathways are responsive to the stimulus. It is evident that the functional pathway dose-response analysis performs consistently better at noise levels of 20% and 30%. Although the performance of the 2 methods is almost comparable at high noise levels of 40%, functional pathway dose-response analysis consistently outperformed nonfunctional Pearson χ2 test under all settings.

Figure 2.

Advantage of functional pathway dose-response analysis over alternative nonfunctional methods. Data were simulated using 17 known cancer pathways in KEGG in response to 3 doses of a hypothetical stimulus at 4 biological replicates. The percentage of dose-responding genes in a pathway was controlled at 10%, 50%, and 90%, and the noise levels were 20%, 30%, and 40%. Areas under (left) ROC and (right) precision-recall curves were calculated for pathway analysis using functional χ2 versus nonfunctional Pearson χ2 on the simulated data. ROC indicates receiver operating characteristic.

Genes with expression patterns most responsive to p27 dosage

We first analyzed the papilloma gene expression data set for single-gene dose-response using FunChisq. The expression of each gene was quantized into 4 levels, using R package “Ckmeans.1d.dp”[39] relying on dynamic programming to guarantee optimal clustering for sample points of each gene. Then, we generated 3 × 4 contingency tables with each gene’s discrete expression level as the column variable and p27 gene dose as the row variable. Then, the FunChisq test was applied to identify genes which are most functionally responsive to p27 dosage. To correct for the multiple testing effect, we control the false discovery rate (FDR) using R package “qvalue”[40] while adjusting the FunChisq P values. When FDR is controlled at 10%, 21 genes are significantly responsive to p27 dosage. Figure 3 illustrates monotonic and nonmonotonic expression patterns of 16 genes with respect to p27 dosage.

Figure 3.

Representative p27 dose-response expression patterns of individual genes. The 3 p27 genotypes are null for nullizygous mutation (p27−/−), Het for heterozygous mutation (p27+/−), and WT for wild type (p27+/+). Null, Het, and WT represent increasing dosages of p27. Each color indicates a response pattern: red for increasing (INC), green for decreasing (DEC), and blue and yellow for nonmonotonic low-high-low (LHL) and high-low-high (HLH), respectively.

Table 1 shows 21 genes with the most significant FDR-adjusted P values, response patterns, and the pathways to which each gene belongs. We categorize the expression responses of each gene to 3 increasing p27 doses into 4 patterns: increasing (INC), decreasing (DEC), low-high-low (LHL), and high-low-high (HLH). We call the former 2 monotonic and the latter 2 nonmonotonic response patterns. The criterion to assign the expression responses of each gene to a pattern is based on its median expression level under each dose. If the medians are strictly increasing or decreasing as p27 dosage increases, the pattern is INC or DEC, respectively. Otherwise, the pattern is either LHL or HLH, if the median at the p27 intermediate dose is either the highest or the lowest among the 3 medians of responses to each p27 dose. It is worth pointing out that while either the correlation or Pearson χ2 test can detect INC and DEC patterns, they cannot capture LHL or HLH patterns.

Table 1.

Top individual genes that are most responsive to p27 dosage.

Symbol	Official gene name	FDR-adjusted P value	Response	Involved pathways
GPC3	Glypican 3	.0172	DEC	Proteoglycans in cancer
EEF2	Eukaryotic translation elongation factor 2	.0172	HLH	AMPK signaling pathway, oxytocin signaling pathway
LMNA	Lamin A/C	.0172	LHL	Dilated cardiomyopathy
ADCY4	Adenylate cyclase 4	.0172	DEC	Purine metabolism, dilated cardiomyopathy, oxytocin signaling pathway
TNFRSF1A	Tumor necrosis factor receptor superfamily, member 1A	.0172	HLH	Apoptosis, tuberculosis, Epstein-Barr virus infection
RYR1	Ryanodine receptor 1 (skeletal)	.0217	DEC	Calcium signaling pathway, oxytocin signaling pathway
PSMB10	Proteasome (prosome, macropain) subunit, beta type, 10	.0478	LHL	Proteasome
DHRS9	Dehydrogenase/reductase (SDR family) member 9	.0973	HLH	Retinol metabolism
CYB5R3	Cytochrome b5 reductase 3	.0973	LHL	Amino sugar and nucleotide sugar metabolism
AP1S1	Adaptor-related protein complex 1, sigma 1 subunit	.0973	LHL	Lysosome
RGL2	ral guanine nucleotide dissociation stimulator-like 2	.0973	INC	Ras signaling pathway
ELOVL1	ELOVL fatty acid elongase 1	.0973	HLH	Fatty acid elongation
SGPL1	Sphingosine-1-phosphate lyase 1	.0973	INC	Sphingolipid metabolism, sphingolipid signaling pathway
CHPT1	Choline phosphotransferase 1	.0973	HLH	Choline metabolism in cancer
MEF2C	Myocyte enhancer factor 2C	.0973	HLH	MAPK signaling pathway, oxytocin signaling pathway
GIT2	G protein–coupled receptor kinase-interacting ArfGAP 2	.0973	INC	Endocytosis
HSD11B1	Hydroxysteroid (11-β) dehydrogenase 1	.0973	DEC	Steroid hormone biosynthesis, chemical carcinogenesis
SERPINA1B	Serine (or cysteine) preptidase inhibitor, clade A, member 1B	.0973	DEC	Complement and coagulation cascades
SAR1B	Secretion associated, Ras-related GTPase 1B	.0973	HLH	Protein processing in endoplasmic reticulum
CYP4F18	Cytochrome P450, family 4	.0973	LHL	Arachidonic acid metabolism
PTGS1	Subfamily F, polypeptide 3 prostaglandin-endoperoxide synthase 1	.0973	HLH	Arachidonic acid metabolism, platelet activation

Abbreviations: DEC, decreasing; HLH, high-low-high; INC, increasing; LHL, low-high-low.

Of these 21 genes, 2 are involved in guanosine triphosphatase (GTPase) signaling, and the GTPase HRAS is the major oncogenic driver of DMBA/TPA-induced papillomas[41,42] and other tumor types. These genes also belong to diverse signaling pathways with oxytocin signaling and arachidonic acid metabolism being most prominent. This demonstrates a remarkable range of p27-driven effects on biological processes beyond cell cycle regulation.

Genes in cancer pathways nonmonotonically respond to p27 dosage

Although single-gene dose-response analysis has highlighted individually responsive genes, we extended the analysis to biologically functional pathways to provide a comprehensive map for the influence of p27 on cancer gene networks with improved statistical power. Therefore, we performed functional dose-response analysis on 281 KEGG pathways with the p27 dose-response transcriptome data set and identified 198 statistically significantly dose-responsive pathways and disease processes at an FDR of 10% (meaning about 20 pathways—10% of 198—might be false positives). The complete ranked list of all 198 pathways and disease processes with pathway dose-response significance P values is provided in Supplementary Table S1. The P values were estimated by gamma distributions. Figure 4 shows that gamma approximation becomes more accurate than the χ2 distribution after a few iterations of permutation and quickly converges to the true distribution afterward. Gamma approximation also converges faster to the true distribution than does the permutation test based on resampling.

Figure 4.

Convergence of gamma approximation to the null distribution in functional pathway dose-response analysis. Each box plot summarizes the deviation of the approximated P values from the true P values. The true P values were computed for all pathways by full permutation of the 12 samples from the papilloma gene expression data. The P values for the red box plot were computed for all pathways using a theoretical χ2 null distribution without permutation. Although both gamma approximation (green boxes) and random permutation test (yellow boxes) converge to the true P values after 5000 iterations, gamma approximation reduces errors faster than the random permutation test.

For each pathway, we also computed the number of member genes which responded in INC, DEC, HLH, and LHL patterns, respectively. We were interested in whether a certain pattern is enriched in that pathway. With the known numbers of genes in each pattern and in total, we computed the probability that a pattern is overrepresented in a pathway using the hypergeometric distribution for the null hypothesis that a pattern is equally represented in the pathway versus in the entire collection of pathways.

Table 2 shows the results for 14 pathways and disease processes known to be affected by p27 from the literature. We can detect all 14 at an FDR of 10%. This provides additional evidence to simulation studies to support the effectiveness of the functional pathway analysis. The overrepresentation statistics show that the LHL pattern is significantly enriched in 2 pathways of cell cycle and microRNAs in cancer, the DEC pattern is enriched in PI3K-Akt signaling pathway, and both HLH and LHL are significantly enriched (P < .05) in the Epstein-Barr virus infection pathway.

Table 2.

Response of 14 pathways and disease processes known to involve p27.

Pathway name	Response to p27	Pattern enrichment P value
	FDR-adjusted P value	DEC	INC	HLH	LHL
With significantly enriched patterns
PI3K-Akt signaling pathway	.0772	.0119	.7687	.8553	.4929
Cell cycle	.0800	.9878	.2272	.9196	.0026
Chronic myeloid leukemia	.0800	.7180	.0504	.5568	.5549
Epstein-Barr virus infection	.0817	.9994	.8413	.0375	.0472
MicroRNAs in cancer	.0950	.3719	.5775	.9782	.0147
Significant in response to p27 at 10% FDR
Hepatitis B	.0635	.6704	.1319	.7249	.3233
Prostate cancer	.0635	.4537	.2326	.3973	.6476
HIF-1 signaling pathway	.0646	.6201	.0881	.6257	.5163
ErbB signaling pathway	.0646	.9792	.2778	.4075	.0912
Measles	.0772	.8567	.6647	.2469	.1889
FoxO signaling pathway	.0800	.5044	.5403	.1180	.7416
Transcriptional misregulation in cancer	.0800	.7843	.1229	.8571	.1357
Viral carcinogenesis	.0805	.9158	.1447	.3806	.3924
Small cell lung cancer	.0865	.8716	.1505	.5975	.2047

Abbreviations: DEC, decreasing; FDR, false discovery rate; HLH, high-low-high; INC, increasing; LHL, low-high-low.

Pattern enrichment P values in bold indicate that the respective patterns are significantly enriched for the pathway or disease process.

Table 3 lists the response patterns and statistical significance of 17 KEGG pathways known to play important roles in cancer. Of the 17 cancer-related pathways, 16 are responsive to p27 gene dosage at 10% FDR. Figure 5 presents a comprehensive view of cancer gene networks with these pathways highlighted for their connections to p27. The networks include important pathways and genes that contribute to cancer hallmarks such as sustained angiogenesis, evading apoptosis, tissue invasion, and metastasis. We note that there are 2 ovals of p27 representing nuclear and cytoplasmic p27, respectively. Cytoplasmic p27, likely to have CDK-independent functions, affects several signaling pathways. It has direct interactions with the PI3K-Akt signaling and mTOR signaling pathways (Figure 5), which mediate sustained angiogenesis and evading apoptosis. Via these 2 pathways, p27 has indirect interactions with the Jak-STAT signaling pathway, the ErbB signaling pathway, and, in turn, the cytokine-cytokine receptor interaction pathway. Another indirect interaction is with the focal adhesion and adherens junction pathways, essential for cell migration and metastasis.

Table 3.

Response of 17 cancer pathways to p27 dosage.

Pathway name	Response to p27	Pattern enrichment P value
	FDR-adjusted P value	DEC	INC	HLH	LHL
With significantly enriched patterns
Focal adhesion	.0635	.0266 [a]	.9250	.9804	.0373
Apoptosis	.0647	.8384	.1165	.9400	.0356 [a]
PPAR signaling pathway	.0706	.1473	.7086	.0072 [a]	.9915
ECM-receptor interaction	.0772	.0061 [a]	.8609	.9925	.0701
PI3K-Akt signaling pathway	.0772	.0119 [a]	.7687	.8553	.4929
Cell cycle	.0800	.9878	.2272	.9196	.0026 [a]
p53 signaling pathway	.0934	.9969	.0656	.7611	.0265 [a]
Significant in response to p27 at 10% FDR
Regulation of actin cytoskeleton	.0635	.5555	.7414	.7072	.0731[a]
ErbB signaling pathway	.0646	.9792	.2778	.4075	.0912[a]
mTOR signaling pathway	.0691	.8817	.1097	.8583	.0619[a]
MAPK signaling pathway	.0706	.5331	.9280	.2044[a]	.2559
Cytokine-cytokine receptor interaction	.0837	.1206[a]	.9058	.4630	.3972
Jak-STAT signaling pathway	.0886	.9541	.6199	.1599[a]	.1782
VEGF signaling pathway	.0886	.5046	.5441	.5445	.2010[a]
Adherens junction	.0934	.4420	.8315	.6284	.0698[a]
Wnt signaling pathway	.0960	.9465	.3935	.2090[a]	.2638
TGF-β signaling pathway	.1070	.2706	.8237	.8137	.0585[a]

Abbreviations: DEC, decreasing; FDR, false discovery rate; HLH, high-low-high; INC, increasing; LHL, low-high-low; TGF-β, transforming growth factor β; VEGF, vascular endothelial growth factor.

Pattern enrichment P values in bold indicate that the respective gene response patterns are significantly enriched for the pathway or disease process.

The most enriched gene response pattern for each pathway (smallest P value in pattern enrichment).

Figure 5.

p27-driven dose-responsive pathways in cancer gene networks. The solid outer box represents the cell membrane. The dashed inner box represents the nuclear membrane. Bold rectangles are known cancer hallmarks. Three blue rounded boxes are pathways known to link to p27. Eight red rounded boxes highlight pathways responsive to p27 according to our study. The figure was simplified from KEGG pathways in cancer. See Supplementary Figure S4 for details.

Table 3 shows 4 pathways enriched for LHL, 3 for DEC, 1 for HLH, and none for INC patterns. We also compared the enrichment statistics among the patterns in the same pathway to reveal the most enriched (smallest P value) pattern. As a result, 10 cancer pathways have LHL as the most enriched, 4 pathways have HLH, 3 pathways have DEC, and none has INC as the most enriched. Evidently, the LHL response pattern is dominant in cancer pathways.

We took a further step to examine the member genes of the cancer pathways. Without prior expectation, Cyclin D1 (Ccnd1)—a cell cycle regulator—emerges as the most responsive gene in several p27-responsive cancer pathways. Human CCND1 is an oncogene which is overexpressed in many tumors, especially in breast cancer.[16,43,44] In terms of responding to p27 in this study, Ccnd1 ranks first in cell cycle, p53 signaling, Jak-STAT signaling, and Wnt signaling pathways. It also ranks third in focal adhesion pathway and fourth in the PI3K-Akt signaling pathway. More importantly, the response of Ccnd1 follows an LHL pattern, suggesting that Ccnd1 is active at the intermediate dosage of p27. Figure 6 shows the response of Ccnd1 to p27 dosage, which is consistent with the model proposed for prostate tumorigenesis by Gao et al.[25]

Figure 6.

Ccnd1 gene expression responses to increasing p27’s dose. The box plot from our experimental data shows that the response follows the pattern of Low-High-Low, suggesting that Ccnd1 is highly expressed at the intermediate level of p27 dosage. WT indicates wild type.

The dominance of nonmonotonic LHL response is evidence for the intermediate dose of p27 as cancer promoting: at this dose, a large number of cancer pathways become highly expressed in favor of tumor development. When p27 dose deviates from the intermediate range, the cellular context may be induced to transition from oncogenic to tumor suppressing. Taken together, these results reveal a deeper understanding of the seemingly contradictory roles of p27 in cancer reported so far in the literature.

Heterozygosity dominates somatic mutations in CDKN1B in human cancers

To determine the mutational profile of p27 in human cancers, we examined somatic mutations in CDKN1B from data in TCGA (http://cancergenome.nih.gov/), where thousands of tumors and matched normal samples were collected for 30 major cancer types. We call a point mutation heterozygous when exactly 1 allele is different from the reference nucleotide, and a homozygous point mutation refers to both alleles differing from the reference nucleotide. The frequency of observed somatic mutations of CDKN1B in human cancers is shown in Figure 7A. A total of 114 tumor samples carry 15 homozygous and 99 heterozygous mutations in CDKN1B. Notably, all 5 observed CDKN1B mutations are heterozygous in 53 colon adenocarcinoma tumor samples, the highest ratio among all cancers (Figure 7) and in consonance with known correlation between p27 expression and poor prognosis of colorectal cancer.[16,17] We further compared the occurrence of mutations in CDKN1B with 5 other known cancer-related genes TP53, BRCA1, BRCA2, RB1, and APC. Because these genes have large variations in the total number of mutations observed from the TCGA data sets, we show the percentages of heterozygous mutations over all mutations in tumor samples for each gene in Figure 7B. The bar plots indicate that CDKN1B heterozygosity is the highest among the 6 cancer genes. The increased percentage of CDKN1B heterozygous mutations in tumor samples relative to these genes are all statistically significant, as indicated by the heterogeneity χ2 test P values in Figure 7C. This result reveals the dominance of CDKN1B heterozygosity in human cancers and is consistent with the intermediate dose of p27 being oncogenic.

Figure 7.

Heterozygosity of CDKN1B mutations in samples of major human tumors. All somatic mutations are obtained from the publicly accessible portion of The Cancer Genome Atlas Web site. (A) A total of 114 CDKN1B mutations were found among tumor samples from 14 of 30 major human cancers. The proportion of heterozygous mutations (orange) over the number of tumor samples is much higher than that of homozygous mutations (blue) in several cancers including colon adenocarcinoma. (B) The percentages of heterozygous mutations over all mutations are shown for each gene. Heterozygous mutations are relatively more prevalent in tumor samples for CDKN1B than 5 other known cancer-related genes TP53, BRCA1, BRCA2, RB1, and APC when examined cumulatively in all major tumor types. (C) Heterogeneity χ2 test P values confirm the strongest CDKN1B heterozygosity in tumors, relative to TP53, BRCA1, BRCA2, RB1, and APC.

We also compared copy number variation (CNV) of CDKN1B in cancer with 2 tumor suppressor genes, PTEN and TP53, as well as 2 oncogenes, MYC and EGFR. The CNV data were extracted from the Catalogue Of Somatic Mutations In Cancer (COSMIC) at cancer.sanger.ac.uk.[45] Pair-wise differences between the copy numbers of the 5 genes are statistically significant (P < .001 after Bonferroni adjustment). Figure 8 shows that the copy number ranges of tumor suppressor genes PTEN and TP53 are lower, which is expected because the low expression of these genes correlates with the loss of cancer suppression. On the other side, the copy number ranges of MYC and EGFR are higher, which agrees with their oncogenic role. CDKN1B, once again, has a medium range of copy numbers, consistent with our observation that CDKN1B has oncogenic properties at the intermediate dosage.

Figure 8.

Copy number variation of PTEN, TP53, CDKN1B, MYC, and EGFR in cancer. Tumor suppressor genes PTEN and TP53 (green) have lower copy numbers, oncogenic genes MYC and EGFR (red) have higher copy numbers, whereas the copy number of CDKN1B (blue) is in between. COSMIC indicates Catalogue Of Somatic Mutations In Cancer.

Discussion

We have presented a novel pathway dose-response analysis method to uncover important cellular processes that might be regulated by p27. This method is unique for pathway analysis because it is functional and nonparametric. Functionality has long been used for causal inference, and nonparametric function dependency enables detection of previously unknown forms of nonmonotonic dose-responses. The pathway-level dose-response statistic, improved on statistical power of single-gene functional χ2, can identify most responsive pathways functionally and nonparametrically. Once a dose-response relationship has been validated, a parametric function[46] can be estimated for a more precise characterization of a specific dose-response curve.

Some human tumors show prominent cytoplasmic staining of p27.[47] Mouse models show exceptions to strict dose dependence of tumor suppression by p27,[25,48] and tissue culture studies show that p27 has numerous cell cycle–independent functions.[2,49] Together, these previous studies suggest that in addition to its well-known CDK inhibitory/tumor suppressor function, p27 may also have oncogenic activity in the cytoplasm, perhaps related to cell motility and metastasis. However, in the absence of in vivo data, the relevance of this remains unclear. Most critically, we do not know if p27 is prometastatic or antimetastatic, and more generally, we do not know the mechanistic basis for the association between p27 expression patterns in human tumors and patient mortality. This study has revealed, in an in vivo tumor model, that the largest proportion of genes on all known cancer pathways responded nonmonotonically to p27 dosage. Our pathway dose-response analysis not only detects direct links between p27 and a pathway but also suggests potential indirect interactions via cascading pathway responses. Among key cancer pathways (Table 3), our analysis established new links from p27 to 8 of them. Here, a link is considered new if it has not been specified in KEGG database. These 8 pathways are adherens junction, apoptosis, cytokine-cytokine receptor interaction, focal adhesion, Jak-STAT signaling pathway, MAPK signaling pathway, and vascular endothelial growth factor signaling pathway. These connections support the CDK-independent functions of p27. Moreover, focal adhesion and adherens junction could be the link via which p27 influences cell migration and metastasis. Indeed, p27 is known to affect cell migration and metastasis via cyclin D1 and RhoA signaling pathway.[49,50] These links to key cancer pathways suggest a more complex role of p27 in cancer. p27 could be much more active in its CDK-independent function than cell cycle inhibition.

We have yet to propose a mathematical model for the biochemical mechanism that is consistent with the observed nonmonotonic dose-response. In control theory, the functional dependency of other genes on p27 with nonmonotonic patterns is decomposable to monotonic subsystems with negative feedback.[51] This is a promising direction for future modeling work. A limitation of this study is that the gene expression data are not specific to protein modifications or subcellular localization. When additional location-specific protein data become available, the presented functional pathway dose-response analysis is readily applicable to deepen the biological findings. Thus, one direction is to explore the connection between the 2 cancer-induced/repressed promoter groups and the abundance and activity of cytoplasmic/nuclear p27 protein.

The presented functional pathway dose-response analysis is a novel method to capture functional dependency of the response of biological systems to stimuli. The main strength of this method is its capability to test for directional influence from dose to response, hence, suggesting potential causal relationships that might have been missed by undirected association analyses. Our study on p27Kip1 shows that it has nonmonotonic effects on expression of genes in almost all cancer pathways beyond its inhibitory role in cell cycle regulation. This complies with findings in the literature that p27 has both tumor suppressing and oncogenic activities, and other data support the generalization of this conclusion to human cancers.