A Four-Gene Signature from NCI-60 Cell Line for Survival Prediction in Non–Small Cell Lung Cancer Free

Lung cancer is the most common cause of cancer deaths worldwide and non–small cell lung cancer (NSCLC) is the major dominant cell type (1, 2). For the early-stage lung cancer, surgery remains the treatment of choice. But the widespread metastasis of lung cancer often defeats this treatment and up to 40% of such patients may relapse within 5 years (3, 4). Therefore, identification of patients at the high risk of recurrence or metastasis after resection may help improve the management of NSCLC patients.

Migration and invasion ability of cancer cells are the most critical initial events of cancer metastasis (5). It is important to find gene signatures that may distinguish the invasion abilities of cancer cells to invade (6). Many studies on cancer metastasis have identified several molecules participating in tumor cell invasion and metastasis in different types of cancers (7). Microarray data have been widely used in discovering tumor subtypes, finding molecular characterization about tumor stage and metastatic status, and predicting patient survival (8–10). Various platforms of microarray have been applied to profile the gene expression patterns of the panel of NCI-60 cell lines (11).

We focused on the subset of the nine lung cancer cell lines in NCI-60 gene expression data set and used one set of public microarray data of western lung cancer patients to derive an invasion-associated gene signature for survival prediction. We validated the findings by three independent NSCLC cohorts. Our computational approach integrated both invasion profiles of lung cancer cell lines and gene expression profiles, and our results provided clues to interpret the mechanism of cancer metastasis.

A cohort that consisted of 69 patients from Taiwan was used for quantitative PCR (qPCR) validation conducted in our lab. Frozen specimens of lung cancer tissue came from 69 patients who underwent surgical resection of NSCLC at the Taichung Veterans General Hospital (VGHTC) between December 1999 and December 2003. The cohort included 27 adenocarcinomas, 36 squamous cell carcinomas, and 6 other types of cancer. The patients had not received adjuvant chemotherapy. The study was approved by the institutional review board of the hospital. Written informed consent was obtained from all patients.

The nine NCI-60 lung cancer cell lines (A549/ATCC, NCI-H460, NCI-H23, NCI-H322M, EKVX, NCI-H226, NCI-H522, HOP-62, and HOP-92) were purchased from the National Cancer Institute. The cells were grown in tissue culture flasks at 37°C in 5% CO₂ in RPMI 1640 with 2 mmol/L l-glutamine and 10% fetal bovine serum.

The invasion capacities for cell lines were examined by using membrane invasion culture system to confirm their invasion ability, respectively. A HTS FluoroBlok insert containing 8-μm pores (Falcon, Becton Dickinson) was coated with 30 μg Matrigel (BD). The cells were suspended in RPMI 1640 containing 10% fetal bovine serum and seeded into the upper wells of the chamber (2.5 × 10⁴ per well). After incubating for 24 h at 37°C, the membrane of the Transwell will be fixed for 10 min at room temperature with methanol and stained for 30 min with 50 μg/mL propidium iodide (Sigma) and the cell number in each blot was counted under a microscope with Analytical Imaging Station system (Imaging Research, Inc.).

We used two microarray gene expression profiling data sets for NCI-60 cell lines generated by Gene Logic using Affymetrix U95 and U133, respectively. The mapping of probes between U133 and U95 platforms was provided by Affymetrix.⁹

The NSCLC cohort we used to derive gene signature consisted of lung cancer patients from the Duke University. Gene expression of 111 lung tumor tissues was profiled by using HU133plus2 Affymetrix chip (12). We retrieved the data from the Gene Expression Omnibus (GSE 3141). Intensity values were log transformed to a base 2 scale.

The two western cohort data sets used for validation came from Shedden et al. (13). The patients were from four institutions, University of Michigan Cancer Center (UM), Moffitt Cancer Center (HLM), Memorial Sloan-Kettering Cancer Center (MSK), and the Dana-Farber Cancer Institute (CAN/DF), and divided into two cohorts, UM+HLM and MSK+CAN/DF, as in Shedden et al. We preprocessed the microarray data by normal score transformation.

The expression of four signature genes and two control genes (TBP for cell lines and GAPDH for tissue specimens) was measured by reverse transcription-qPCR (RT-qPCR) with specific Taqman probes and primer sets; the transcripts were amplified with Taqman One-Step RT-PCR Master Mix Reagent (Applied Biosystems) and a detection system (ABI Prism 7900HT, Applied Biosystems). Gene expression was quantified in relation to the expression of the respective control gene with the use of sequence detector software and the relative quantification method (Applied Biosystems). The choice of TBP and GAPDH as the internal control for real-time RT-PCR is based on their invariance of expression in cancer cell lines and clinical cancer specimens (14, 15).

The association between survival and each gene expression profile from microarray is evaluated by univariate Cox proportional hazard regression analysis. A linear combination of the gene expression levels weighted by the regression coefficients was used to calculate the risk score for each patient.

More specifically, a patient's risk score was calculated as the sum of the levels of expression of each gene measured by microarray multiplied by the corresponding Cox regression coefficients. Patients were classified as having a high-risk signature or a low-risk gene signature, with the median of risk score as the threshold value. Survival curves for both groups were obtained by the Kaplan-Meier method and compared using the log-rank test. Multivariate Cox proportional hazard regression analysis with stepwise selection was used to evaluate independent prognostic factors associated with patient survivals, taking the gene signature, age, sex, stage, and histology as the covariates.

Except for the two places noted below, all tests were two sided, and a P value of <0.05 was considered statistically significant. One-sided tests were used in showing the consistency between qPCR data and the microarray data (^‡P values; Table 1), and in one-sided log-rank test (P values of UM+HLM and MSK+CAN/DF cohorts) where the hypothesis of interest is the one-sided hypothesis that the survival time is longer for patients with the low-risk gene signature than that for patients with the high-risk signature, note the two-sided hypothesis that the survival time is different between the low-risk signature and the high-risk signature (16, 17).

The strategy of our approach is illustrated in Fig. 1. Based on invasion ability of the nine lung cancer cell lines, we divided them into two groups. We identified the invasion-associated genes from two Affymetrix NCI-60 microarray gene expression data sets by finding the genes with differential expression between the high and the low invasion groups. To reduce effect of microarray noises, we filtered out 50% of the probes with low average expression readings in data preprocessing.

To examine which genes from our cell line–based search may lead to useful clinical biomarkers, we obtained a large-scale microarray data set for tumor samples of western NSCLC patients together with patients' survival data from Gene Omnibus (GSE 3141), referred to as the Duke cohort (n = 111; ref. 12). We applied the univariate Cox proportional hazard method to each of the invasion-associated genes obtained from cell line assay for finding genes whose expressions across patients correlate well with patient survival. Because genes with higher (lower) expression levels for high invasion cell lines are expected to be risk (protective) genes, we further inspect the direction of correlation with survival to identify genes consistent with their invasion abilities. For the genes that bypassed the consistency checking, we confirmed their expression pattern on nine cell lines by qPCR assay on the nine cell lines. Finally, we validate this in silico–derived gene signature in three independent cohorts: two new U133plus2 microarray data of western NSCLC cohort (CAN/DF+MSK, n = 186; UM+HLM, n = 257; ref. 13) and a group of 69 Taiwanese NSCLC samples (VGHTC, n = 69) by qPCR assay.

We determined the invasion ability of nine human NSCLC cell lines in NCI-60 (Fig. 2A) and found three cell lines, A549, HOP-62, and HOP-92, with much higher invasion ability than the other six cell lines (Fig. 2A). To find invasion-associated genes from the U133plus2 NCI-60 microarray data set, gene expression profiles for the high and the low invasion groups were compared by two-tailed t tests and confirmed with the U95 NCI-60 microarray data set. A total of 110 genes were found differentially expressed and their expression profiles formed two clusters (see Fig. 2B, wherein the top panel showed the expression profiles of 72 genes that correlated negatively with the invasive ability of lung cancer cell lines, whereas the bottom panel gave the expression profiles of 38 genes that correlated positively).

The 110 genes associated with the invasion activity were identified by in vitro invasion assay on NSCLC cell lines and in silico microarray analysis. How do they affect the clinical outcome of lung cancer patients? To investigate this issue, we take advantage of the availability of the Duke lung cancer cohort (12) and evaluated the correlation between the expression of each gene and the patient survival by univariate Cox proportional hazard regression. This analysis was done to determine which genes may be associated with patient death. Protective genes have hazard ratio (HR) for death of <1, whereas risk genes have HR for death of >1. Of the 110 invasion-associated genes, 9 genes were found statistically significant. We further checked if the protective/risk gene assignment was consistent with the direction of the invasion ability of the gene and found five genes consistent. One gene probe was located at a pseudogene site, which was not considered further. For the remaining four genes, ANKRD49 and LPHN1 are protective and RABAC1 and EGLN2 are risk genes (Table 1).

Risk scores were calculated by using the above four genes. Risk prediction was calculated as the sum of the levels of expression of each gene measured by microarray, multiplied by the corresponding Cox regression coefficients. Patients in high- or low-risk group were classified by using the 50th percentile (median) of the risk score distribution (median, 3.84; range, 2.5-6.68) as the cut point.

The Kaplan-Meier method was used to estimate overall survival distribution for each group. Patients with high-risk gene signature as predicted by the four-gene signature had shortened median overall survival compared with patients with low-risk gene signature (24.33 versus 36.59 months; P = 0.001, log-rank test; Fig. 3A). The expression levels of these four candidate genes in the nine lung cancer cell lines were validated by doing the RT-qPCR assay and calculating the Pearson correlation between the qPCR data with the U133plus2 microarray data (Table 1).

Clinical characteristics of the NSCLC cohort are summarized according to their high or low risk of the four-gene signature (Supplementary Table S1). Multivariate Cox regression analysis with stepwise selection was used to search for independent prognostic factors associated with survival. The risk score by the four-gene signature was evaluated along with age, gender, tumor stage, and cell type. Patients with the high-risk signature had a higher risk for poor survival with adjusted HR of 2.35 [95% confidence interval (95% CI), 1.367-4.025; P = 0.002]. The four-gene signature was still a significant survival predictor after adjusting for tumor stage (Table 2A).

To investigate to what extent the four genes were also valuable predictors for the Asian patients, we investigated the VGHTC cohort. We took lung cancer specimens from each patient and quantified the gene expression of four candidate genes by RT-qPCR to compute the patient risk scores for each patient and then divided patients into high- and low-risk groups in the same way as was done in the Duke cohort. We found that Kaplan-Meier survival curves for the high- and low-risk groups in the VGHTC cohort were separable and the result was statistically significant (P = 0.0184, log-rank test; Fig. 3B). To show that inclusion of the four-gene signature in the clinical model had useful survival prediction performance independent of cancer stage, we conducted a multivariate Cox regression with backward deletion. From the initial set of five variables (age, gender, stage, cell type, and four-gene signature), the final model contains only age, stage, and four-gene signature (Table 2A).

Moreover, we validated the four-gene signature in two new western cohorts (UM+HLM and CAN/DF+MSK). Kaplan-Meier survival curves were significantly separable (UM+HLM: P = 0.002, log-rank test; CAN/DF+MSK: P = 0.017, log-rank test; Fig. 4A and B), and four-gene signature was significant by multivariate Cox regression in both cohorts (Table 2B).

This reassures the clinical usefulness of our finding. The four invasion-associated genes might be potential molecular targets for lung cancer therapy.

Many studies have developed useful gene signatures by examining only the gene profiling data of the patient specimens that they can access to (18–20). Differing from these works, our approach started with an in silico exploitation of the public domain microarray databases to derive the gene signature and concluded with the validation by our own clinical data as well as two public western clinical data. The public data we used in this report included two sets of NCI-60 cell line microarray gene expression data and three microarray gene expression data sets from a western patient cohort. We also conducted qPCR in our lab. As the genomic research gets deeper and broader, it has become more and more difficult for any research team to generate all sorts of data on their own. Consequently, it is worth investigating methods of how to use the public data such as Gene Expression Omnibus, which is growing at an astonishingly fast pace. More and more such in silico microarray and bioinformatic approaches are anticipated to foster the area of clinical genomics.

Previously, a five-gene signature (9), motivated by different cell line invasion assay and microarray platform, was also obtained for NSCLC. Although there is no overlap between the four genes identified in this article and the five-gene signature, both achieve the same goal. Both signatures had significant prediction power in multivariate Cox regression. Using both signatures, we stratified the VGHTC cohort into three groups and generated Kaplan-Meier plots to show their survival differences (Fig. 4C). The combined signature is significant under multivariate Cox regression with stage, gender, and age (Table 2B).

Among the four genes in our invasion-associated signature, the expression of RABAC1 and EGLN2 was found higher in high invasion lung cell lines, whereas that of LPHN1 and ANKRD49 was found higher in low invasion lung cell lines. RABAC1 (Rab acceptor 1) is a transmembrane protein, localized to the late Golgi, which also called PRA1 (21). It was involved with T-cell factor/β-catenin signaling and dephosphorylation of extracellular signal-regulated kinase 1/2 (22). RABAC1 gene was highly expressed in breast, colon, lung, and ovary cancers (23). Interestingly, RABAC1 activated NF-κB signaling through latent membrane protein 1, which is an EBV-encoded oncoprotein (24). Moreover, infection of EBV was associated with several cancers (25), consistent with risk effect of RABAC1 in our signature. The EGLN2 was known as prolyl-hydroxylase-1, which acted like an oxygen sensor to inhibit hypoxia-induced activation of hypoxia-inducible factor-1α (26, 27). EGLN2 was also characterized as estrogen-induced tag-6, which has been reported to promote cell growth in breast cancer (28). It is likely that the function of EGLN2 in tumor tissues may be dependent on the oxygen level. LPHN1 encoded a member of the latrophilin 1, which was a member of G protein–coupled receptors associated with activation of neurons (29). Latrophilin 1 regulated cell adhesion and neurotoxin secretion (30) and brain-specific angiogenesis inhibitor 1. Brain-specific angiogenesis inhibitor 1 was widely expressed in both glial cells and neurons but has low expression in glioma cell lines (31). The function of ANKRD49 (ankyrin repeat domain 49) is little known. However, ankyrin repeat, a motif of 30– to 34–amino acid residues (32), was first identified in the sequence of yeast Swi6p, Cdc10p, and Notch (33), and Notch signaling is one important pathway in tumorigenic processes. Families of ankyrin repeat proteins mediating protein-protein interaction have been associated with development of human cancer (34). The tumor suppressor gene p16 (INK4A) was an ankyrin repeat protein that inhibited cell cycle progression to reduce proliferation of cancer cells (35).

We concluded that an invasion-associated gene signature derived from NCI-60 lung cancer cell lines by our in silico microarray approach had good survival prediction power in NSCLC patients. The signature requires only four genes. For clinical applications, the smaller number of genes in the signature may have the greater potential of translating into practical methods of lung cancer patient risk stratification and developing personalized therapy.

No potential conflicts of interest were disclosed.

We thank Hao Ho for contribution in data assembly.