|Home | About | Journals | Submit | Contact Us | Français|
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
The cell cycle machinery interprets oncogenic signals and reflects the biology of cancers. To date, various methods for cell cycle phase estimation such as mitotic index, S phase fraction, and immunohistochemistry have provided valuable information on cancers (e.g. proliferation rate). However, those methods rely on one or few measurements and the scope of the information is limited. There is a need for more systematic cell cycle analysis methods.
We developed a signature-based method for indexing cell cycle phase distribution from microarray profiles under consideration of cycling and non-cycling cells. A cell cycle signature masterset, composed of genes which express preferentially in cycling cells and in a cell cycle-regulated manner, was created to index the proportion of cycling cells in the sample. Cell cycle signature subsets, composed of genes whose expressions peak at specific stages of the cell cycle, were also created to index the proportion of cells in the corresponding stages. The method was validated using cell cycle datasets and quiescence-induced cell datasets. Analyses of a mouse tumor model dataset and human breast cancer datasets revealed variations in the proportion of cycling cells. When the influence of non-cycling cells was taken into account, "buried" cell cycle phase distributions were depicted that were oncogenic-event specific in the mouse tumor model dataset and were associated with patients' prognosis in the human breast cancer datasets.
The signature-based cell cycle analysis method presented in this report, would potentially be of value for cancer characterization and diagnostics.
A fundamental characteristic of all cancers is cell cycle deregulation . Although diverse factors such as point mutation, gene amplification, activation of oncogenes, inactivation of tumor suppressors, and hypermethylation are involved in cancer development, their influence ultimately is on the cell cycle machinery. Therefore, various methods of cell cycle phase estimation have been developed. The M phase indicator mitotic index, the number of mitotic bodies in a microscopic field, and the S-phase fraction, a DNA flow cytometry determination, are used to measure the tumor proliferation rate and are predictive for breast cancer prognosis [2-4]. Immunohistochemistry (IHC) against cell cycle markers is another tool. For example, the expression of G1-S transition marker cyclin E, S-G2 marker cyclin A, or S-G2-M marker geminin are predictive of poor prognosis of breast cancers [2-5]. However, these methods rely on one or few measurements and consequently provide a limited scope of information. There is a need for more systematic methods of cell cycle phase analysis, such as microarray-based techniques [3,4].
Gene expression signatures, which are capable of predicting the state of a sample from a given microarray dataset, are the emerging technology for developing cancer therapeutics. The "70-gene signature" from a breast cancer dataset has shown predictive power for the risk of recurrence . The "pathway deregulation signature" has shown the ability to predict pathway status and to characterize breast, lung and ovarian cancers . The "chemotherapy response signature" has accurately predicted clinical response to cytotoxic drugs for breast and ovarian cancers . Here, we report the development of the "cell cycle signature (CCS)" which indexes the cell cycle phase distribution from microarray profiles considering both cycling and non-cycling cells. The CCS method depicted "buried" cell cycle phase distributions that were oncogenic-event specific in a mouse tumor model dataset and were associated with patients' prognosis in human breast cancer datasets. The method has a potential to be of value in the characterization and diagnosis of cancers.
To analyze cell cycle phase distribution, a series of CCSs were created as described in Methods (Fig. (Fig.1A,1A, Additional file 1). The CCS masterset, 252 genes that express preferentially in cycling cells and in a cell cycle-regulated manner, represents the entire cell cycle and is henceforth denoted as CCScycling. Eighteen CCS subsets, each composed of genes whose expressions peak at a specific stage of the cell cycle, represent the phases of the cell cycle and are denoted using the subscript naming convention of CCSphase. For example, the CCS subsets for the G1 phase are expressed as CCSG1, for the G2-M phase as CCSG2-M, and so on.
Solid tumors are composed of various proportions of cycling and non-cycling cells , and cell cycle phase distributions can be assessed as per total cells or as per cycling cells. Since microarray measurements are the net expression of all cells in the sample, the data is generally per total cells. To obtain data per cycling cells from a given microarray dataset (Fig. (Fig.1B,1B, total gene dataset), a subdataset is created by extracting the expression values of CCScycling genes (Fig. (Fig.1B,1B, cycling gene dataset). Then, both the total and the cycling gene datasets undergo quantile normalization which gives the same expression value distribution for each sample . In the total gene dataset, normalization is done on all genes. On the other hand, in the cycling gene dataset, normalization is done only on the cycling genes. Because genes in the CCScycling preferentially express in cycling cells, the influence of non-cycling cells would be limited for the cycling gene dataset. Scores for each CCS are calculated for both datasets. CCScycling and CCSphase scores for the total gene dataset could index the proportion of cycling cells and of cells at the designated cell cycle phase per total cells, respectively. Similarly, CCSphase scores for the cycling gene dataset could index the proportion of cells at the cell cycle phase per cycling cells. CCScycling scores for the cycling gene dataset could index the proportion of cycling cells per cycling cells and thus would show constant values.
In the preliminary analysis of the Whitfiled et al. cell cycle dataset , CCS indexed cell cycle phase distribution as expected (Additional file 2). To confirm that the CCS method is valid for independent datasets, a cell cycle dataset of synchronized HCT116 cells was prepared and analyzed. As shown in Fig. Fig.2A,2A, similar heat map patterns were observed for the total and the cycling gene datasets. Differences in the CCScycling scores for both the total and the cycling gene datasets were slight in the situation where most cells were expected to be in the cell cycle. Peaks in the CCSphase scores shifted according to cell cycle progression (Fig. (Fig.2A,2A, DMSO 0–10 h), and peaks ceased around the M phase in cells treated with the mitosis inhibitor nocodazole (Fig. (Fig.2A,2A, Ncz 7–10 h), consistent with DNA flow cytometry measurements (Fig. (Fig.2B).2B). The CCS method was able to index cell cycle phase distribution even for an independent cell cycle dataset derived from a different cell line and a different platform.
Solid tumors are not solely composed of cycling cells but contain various numbers of non-cycling cells . Theoretically, changes in the proportion of cycling cells in the sample are expected to evenly change the proportion of cells in all cell cycle phases. To examine the influence of changes in the proportion of cycling cells on CCS scores, analysis was conducted on the Fournier et al. dataset  of profiles of human mammary epithelial cells (HMECs) cultured in leucine-rich extra cellular matrix. In this system, HMECs grow exponentially and then enter a quiescent state [12,13]. As shown in Fig. Fig.2C,2C, CCScycling and CCSphase scores for the total gene dataset uniformly decreased as the HMECs transitioned from cycling (day 3) to non-cycling state (day 7) (Fig. (Fig.2C,2C, upper panel). According to the DNA flow cytometry estimation in the original report, the S phase and G2+M phase fraction size decreased from 15% ± 5.1 (day 5) to 5.5% ± 0.5 (day 7), and from 12% ± 1.1 (day 5) to 7% ± 2.5 (day 7), respectively (day 3 data was not available) . On the other hand, the G0+G1 phase fraction size increased from 73% ± 6.3 (day 5) to 86% ± 4.6 (day 7). Due to the inability of DNA flow cytometry to distinguish cells in G0 from cells in G1, decisive conclusions cannot be made. However, from two situations in which 1) 3D cultured HMECs gradually underwent growth arrest and 2) CCSG1 scores decreased at day 7, this increase can be regarded as an increase in the number of cells at the G0 phase as well as a decrease in the number of cells at the G1 phase. To our surprise, the heat map for the cycling gene dataset showed increasing CCSG1 scores towards day 7 (Fig. (Fig.2C,2C, lower panel). This increase in CCSG1 scores could be due to the G1 phase prolongation which is known to occur under G0-inducing conditions, such as serum starvation and development [14,15]. For further confirmation, we analyzed the Cam et al. dataset  of profiles of growing and serum starved T98 breast cancer cells. Similar to the results for HMECs, a uniform decrease in CCScycling and CCSphase scores for the total gene dataset was observed in serum-starved cells (Fig. (Fig.2D,2D, upper panel). In addition, an increase in CCSG1 scores for the cycling gene dataset was observed (Fig. (Fig.2D,2D, lower panel), indicating prolongation of the G1 phase. Taken together, these results suggested that changes in the proportion of cycling cells in the sample can be presented as uniform changes in CCScycling and CCSphase scores for the total gene dataset.
The mammalian cell cycle is a highly regulated and conserved process . To investigate whether CCS derived from human datasets can be used to closely related species, the Yamamoto et al. dataset , cell cycle profiles (G0 to S) of NIH3T3 mouse fibroblasts, was analyzed. The heat map showed changes in the proportion of cycling cells (Additional file 3: upper panel) as well as cell cycle progression from G1 to S phase (Additional file 3: lower panel), as quiescent cells (FGF 0 h) re-enter the cell cycle, progress through G1 phase and enter S phase (FGF 12 h). These results showed that the human CCS created in this study can be applied for the analysis of mouse datasets.
The CCS method was applied to the Herschkowitz et al. dataset  which contains 122 profiles of 13 different mouse mammary carcinoma models and normal samples. The authors reported that some models developed similar tumors (homogeneous models) of gene expression and histological phenotype while other models showed heterogeneity (heterogeneous models) and gave "randomness of the molecular basis of tumor initiation" as the reason for the heterogeneity. As shown in Fig. Fig.3A,3A, CCScycling and CCSphase scores for the total gene dataset for the normal samples were consistently very low, while scores for tumors were varying degrees higher, indicating variation in the proportion of cycling cells. It is reasonable that heterogeneous models show variation in CCScycling and CCSphase scores. However, variation was also seen in each homogeneous model, although Tag models had a tendency towards higher scores and the Neu model had a tendency towards lower scores. In contrast, CCSphase scores for the cycling gene dataset were similar within the same homogeneous models, except in the Myc model (Fig. (Fig.3A,3A, lower panel). To illustrate this in detail, CCSphase scores of several models for both datasets were plotted as shown in Fig. Fig.3B.3B. It can be seen that each model has a specific cell cycle phase distribution. High CCSG1 and low CCSS-G2-M scores were seen in the Neu model. The opposite pattern was seen in one of the Tag models. The Myc model showed two different cell cycle phase distributions (Additional file 4) and the reason is not clear. However, because Myc has been reported to induce genomic instability and to contribute to tumorigenesis through a dominant mutator effect , additional oncogenic events may have been induced. In all cases, plots for the total gene dataset were vertically shifted in varying degrees which would be due to the influence of non-cycling cells, as presented in HMECs and T98 cells. On the other hand, plots for the cycling gene dataset showed minimal variation in alignment. These results indicated two findings: (i) the cell cycle phase distribution reflects the oncogenic events in tumors, and (ii) the cell cycle phase distribution can be better indexed when the influence of non-cycling cells is taken into account. The advantage of the CCS method can be underscored considering that the current cell cycle phase estimation methods relying on one or few measurements are not sufficient to depict cell cycle phase distribution or to distinguish non-cycling cells.
The CCS method was applied to the Ivshina et al. dataset  from a panel of 249 human breast cancers. The heat map for the total gene dataset showed various CCScycling scores, indicative of variations in the proportion of cycling cells in the sample (Fig. (Fig.4A,4A, upper panel). The CCSphase scores were not uniformly changed in some patients, suggesting that cell cycle phase distributions were also altered. The heat map for the cycling gene dataset displayed a rolling wave pattern (Fig. (Fig.4A,4A, lower panel). Patients with high CCScycling scores for the total gene dataset had high CCSS-G2-M and low CCSG1 scores for the cycling gene dataset, but several exceptions existed (Fig. (Fig.4A),4A), reminding the influence of non-cycling cells found in the analysis of mouse tumor models. Clinical annotations were available for this dataset and so the relevance between CCS scores and patient prognosis were tested. Patients were dichotomized by the median of each CCS score and then the risk differences between the two groups for disease free survival (DFS) were assessed using log-rank test and Cox univariate analysis (Fig. (Fig.4B).4B). The CCScycling score for the total gene dataset was significantly predictive of poor prognosis (Hazard ratio [HR] = 1.98, p = 0.00134) (Fig. (Fig.4B4B and Fig. Fig.4C,4C, CCScycling), consistent with the common view that a larger number of cycling cells correlates with worse clinical outcome. The CCSS-G2-M and several CCSG1 scores for the total gene dataset were also predictive of poor prognosis. On the other hand, CCSG1 scores for the cycling gene dataset had an adverse prognostic power and gave the highest prognostic value among the tests (HR = 0.41, p = 0.0000367) (Fig. (Fig.4B4B and Fig. Fig.4C,4C, CCSG1).
To exclude the possibility of dataset specificity, the CCS method was also applied to the Langerød et al. dataset  from a panel of 80 breast cancers. Similar results were obtained (Additional file 5). For the total gene dataset, variations in CCScycling scores and non-uniform changes in CCSphase scores in some patients were observed. Patients with high CCScycling scores for the total gene dataset had high CCSS-G2-M and low CCSG1 scores for the cycling gene dataset with some exceptions. CCSG1 scores for the cycling gene dataset were predictive for DFS as with the Ivshina et al. dataset and gave the highest prognostic value (HR = 0.41, p = 0.00553) (Additional file 5). Taken together, these results indicated that: (i) variations in the proportion of cycling cells exist among tumors, (ii) the proportion of cycling cells correlated to the cell cycle phase distribution per cycling cells with several exceptions, and (iii) the cell cycle phase distribution per cycling cells better associated with patients' prognosis.
In this study, we developed a signature-based method to index cell cycle phase distribution from microarray profiles under consideration of cycling and non-cycling cells, providing two sources of valuable information on cancers.
One source of information is the proportion of cycling cells in the sample. The rationale of most current cell cycle phase estimation methods, including mitotic index, S phase fraction and IHC against cell cycle markers, is that the high proliferative tumors leading to poor prognosis contain more cycling cells. In the analysis of the human breast cancer datasets, higher CCScycling scores for the total gene dataset, indicative of a larger number of cycling cells in the sample, did associate with poor prognosis. Naturally, it can be thought that an increase in the number of cycling cells leads to a uniform increase in the number of cells at all cell cycle phases. However, some patients showed non-uniform changes in CCSphase scores for the total gene dataset (Fig. (Fig.4A,4A, upper panel), suggesting that each cell cycle phase was not evenly changed. Similarly, Whitfield et al. observed that some cell cycle-regulated genes did not express in correlation with proliferation status in some breast cancers . Furthermore, although the G1 phase is a part of the cell cycle, G1 phase marker cyclin D1 often negatively correlates with poor prognosis of breast cancers [2-4,23]. Therefore, considering only the proportion of cycling cells seems insufficient.
The other source of information is cell cycle phase distribution. A number of oncogenic events are known to perturb the duration of cell cycle phases. For example, activation of oncogenes such as v-H-ras, v-Src, v-Raf, cyclin D1, cyclin E, and c-myc shortens the G1 phase [24-26]. Loss of tumor suppressor Pten shortens the G1 phase  and loss of Lzts1 and Lats2 shortens the M phase [28,29]. Viral infections such as SV40-Tag and HTLV-1 Tax also shorten the G1 phase [30,31]. Such perturbations in the cell cycle phase duration subsequently alter the cell cycle phase distribution. Thus, the cell cycle phase distribution per cycling cells would reflect the biology of cancers. Actually, in the analysis of mouse tumor models, oncogenic-event specific cell cycle phase distributions were observed. This suggests that the cell cycle phase distribution under consideration of both cycling and non-cycling cells has a potential for cancer characterization.
A model of tumors with different cell cycle phase distributions is proposed in Fig. Fig.5.5. Oncogenic events perturb the cell cycle each in a unique way which in turn alters the cell cycle phase distribution as well as the proliferation rate. High proliferative tumors grow rapidly and thereby produce a large number of cycling cells. The opposite is the true for low proliferative tumors. However, high proliferative tumors with a small number of cycling cells or low proliferative tumors with a large number of cycling cells would exist at a low probability. This model would account for non-uniform changes in CCSphase scores for the total gene dataset found in some breast cancer patients, the Whitfield et al.'s observation, and the adverse prognostic value of cyclin D1. Current cell cycle phase estimation methods are insufficient for detecting such cancers. Mitotic index and S-phase fraction do not recognize non-cycling cells. Combinatorial IHC  still needs improvement and validation. Shetty et al. reported a relationship between breast cancer grade and G1 phase length estimated from the ratio of geminin and Ki67 IHC measurements; however, it was not significant . The CCS method, on the other hand, indexed the cell cycle phase distribution under consideration of cycling and non-cycling cells, and showed a potential for characterizing cancers.
Previously, as an alternative microarray-based cell cycle analysis technique, Lu et al. introduced the "expression deconvolution" method . To predict the cell cycle phase distribution of yeast, they prepared about 700 equations with 5 variables representing 5 cell cycle phases and searched for the optimal solution. The method has comparable or even better potential to improve cancer characterization than the CCS method. However, it requires a tremendous amount of computational resources to find the optimal solution and avoid the local minimum, especially as the number of variables increases (18 + 1 phases were analyzed in our study). There are some hurdles that need to be overcome before high resolution cell cycle phase analysis is practical and we are currently tackling some of them.
The HCT116 colorectal cancer cell line (ATCC) was grown in McCoy's 5A medium modified (Sigma-Aldrich) with 10% FBS (JBS) and maintained at 37°C and 5% CO2. Synchronous culture was obtained by incubating cells for 19 h in 2 mM of thymidine, followed by a 9-h incubation in normal medium and a second 16-h incubation in thymidine (2 mM). Cells were washed with normal medium followed by treatment with DMSO for 0, 2, 4, 6, 7, 8, 9, and 10 h as a control or 0.1 mg/ml nocodazole (Sigma-Aldrich) for 7, 8, 9, and 10 h. Cells were stained with propidium iodide and analyzed with DNA flow cytometry.
Total RNA was reverse transcribed, labeled, and hybridized to Human Genome U133 Plus 2.0 arrays (Affymetrix) according to the manufacturer's instructions. The expression value for each probe was calculated using the GC-RMA algorithm. The microarray data were deposited in the GEO database (GEO number: GSE14103).
Two datasets were used to create the CCS. First, the Whitfield et al. dataset  of 47 profiles of synchronized Hela S3 cells for 0–46 h time points (1-h intervals) after release of double thymidine block was analyzed to identify genes which express in a cell cycle-regulated manner. Raw signal intensities from the Cy5 and Cy3 channels were quantile normalized for each sample. Cy5/Cy3 ratios were log-transformed and quantile normalized across the arrays. Resulting values were smoothened using a moving average with a window size of 3 and were standardized by Z-transformation. Then, Fourier transformations were applied to each probe for 1-40-h periods in 15-min increments to identify periodicity and phase offset. Fourier transformation magnitudes for the known 51 cell cycle-regulated genes (listed in Whitfield et al. ) demonstrated a peak at the 14.75-h periodicity (Additional file 6). Thus, probes were selected using the criterion of
where Pi is the Fourier transformation magnitude of the 14.75-h periodicity for probe i, i = 1,..., 44,160. The analysis yielded a list of 1,633 periodically expressed probes representing 976 genes. Second, the Bar-Joseph et al. dataset  of 17 profiles of synchronized primary human foreskin fibroblasts (FFs) for 0–32 h time points (2-h intervals) after release of double thymidine block and 2 profiles of serum starved FFs was investigated to identify genes which preferentially express in cycling cells. Serum starved cells are known to exit the cell cycle phase and to enter the non-cycling G0 phase , thus probes, whose expression is constantly higher throughout the cell cycle compared with non-cycling cells, were selected by the criterion
where eij is the expression value for probe i of serum-starved FFs sample j, j = 1, 2, and eik is the expression value for probe i of the synchronized FFs sample k, k = 1,..., 17. This yielded 2,304 out of 22,277 probes representing 1,779 genes. Then, from the intersection, a list of 335 probes representing 252 genes was obtained. These genes which preferentially express in cycling cells and in a cell cycle-regulated manner compose the CCS masterset (CCScycling). A number of well-known proliferation markers such as Ki67, geminin, TOP2A, aurora A, and PCNA [1-5,32] were included in this signature, while some cell cycle-regulated genes such as p21 and cyclin G1 whose expression can be up-regulated in non-cycling cells [36,37] were not. Lastly, according to their phase offsets, probes for CCScycling were assigned to 18 CCS subsets (CCSphase) which correspond to a 360° cell cycle evenly divided into 20° increments, so that each CCS subset contains at least 3 genes. Because some genes were represented by multiple probes, the same genes may appear in different CCS subsets. The CCS gene list is shown in Additional file 1.
The given microarray dataset was used as the total gene dataset. The cycling gene dataset was created by extracting the expression values for CCScycling constituents from the total gene dataset. Both total and cycling gene datasets then underwent the following steps independently to give CCS scores. Expression values were log-transformed, quantile normalized to achieve the same expression value distribution for each sample, and standardized with Z-transformation across the samples. The Z-scores of the probes for each CCS genes were averaged for each sample and used as the CCS scores. To obtain robust scores, each CCSphase score was adjusted by averaging with the neighboring CCS scores twice for a total of two cell cycle rounds. Heat maps were created by "Java Treeview" . In the analysis of the mouse tumor model dataset, gene ID mapping was done using human-mouse orthology information from HomoloGene . In the analysis of human breast cancer datasets, patients were ordered by peak in CCSphase scores for the cycling gene dataset.
Patients were dichotomized by the median of each CCS score. To assess the risk difference between two groups for DFS, Kaplan-Meier survival analysis, log-rank test and Cox univariate analysis were conducted using R "survival" package.
HM and KK designed the research. HM and YN performed the research. HM, NI, AS and KK participated in writing the manuscript. All authors read and approved the final manuscript.
The gene list for cell cycle signatures. The CCS genes and assigned CCS subset IDs are listed.
Validation of CCS method in the Whitfiled et al. cell cycle dataset. CCS scores were calculated for the total (upper panel) and the cycling (lower panel) gene dataset. The purple bars above the columns indicate Whitfield et al.'s estimations of the S phase.
Analysis of the Yamamoto et al. dataset. Serum starved NIH3T3 cells were stimulated with FGF to re-enter the cell cycle. Profiles of unstimulated cells (FGF 0 h) and FGF-stimulated cells (FGF 3–12 h) were analyzed.
CCS score plots for the WAP-Myc model. Same as for Fig. Fig.3B3B.
Analysis of the Langerød et al. breast cancer dataset. (A), (B) and (C) are the same as in Fig. Fig.44.
Power spectrum of the 51 cell cycle-regulated genes. The Hela S3 cell cycle dataset was processed as described in Methods. Fourier transformation magnitudes for the known 51 cell cycle-regulated genes for each periodicity were averaged and plotted.
We thank D. Schmitt, F. Ford, K. Takahashi, H. Ohmori, M. Haramura, M. Ashihara and M. Aoki of Chugai Pharmaceuticals for their helpful discussions and checking of the manuscript.