Metagenes Associated with Survival in Non-Small Cell Lung Cancer

Patients and Tumor Samples

The Ethical Committee of Tartu University approved the sample collection based on signed informed consent with the patients. Normal samples were collected from tissues adjacent to the patients’ tumors, and were confirmed to be non-cancerous by pathologists. The histological classification of the carcinomas was conducted according to the standards of the World Health Organization (WHO) classification method for carcinoma.^11–14 Eighty five patients underwent surgical resection and the tumors were pathologically confirmed as NSCLC pulmonary carcinoma in Tartu University Hospital between November 2002 and December 2006. Of the 85 tumors, 62 were SCCs and 23 were BCs and ACs. According to the guidelines of the American Joint Committee on Cancer ¹⁵ the patients were staged after the surgery (Supplemental Table 1).

Statistical analysis of gene expression data

Quantile-normalized and log-transformed expression data was obtained from our previous study, ¹⁶ which used Illumina Sentrix BeadChip (HumanWG-6_V2) to profile gene expression. Aforementioned array provides genome-wide transcriptional coverage of well-characterized genes, gene candidates and splice variants, with a significant portion targeting well-established sequences of 18,072 genes. Our initial dataset consisted of 85 lung tumor samples and 21 adjacent cancer-free lung samples (Supplemental Table 1).

Patients who had received preoperative chemotherapy were excluded from differential gene expression analyses leaving 78 tumor samples and 20 adjacent control samples.

Moderated two-sample t-test from R package LIMMA ¹⁷ was used to find differentially expressed genes between sample groups, using FDR-corrected P-value cut-off 0.05 and fold-change 2. The approach was used to analyze differential gene expression between all tumor samples (n = 78) and adjacent controls (n = 20), between SCC samples (n = 58) and AC/BC samples (n = 20), as well as between stage I tumor samples (n = 56) and all the other tumor samples (n = 22).

The differentially expressed genes and all samples were clustered hierarchically using Pearson correlation distance with average linkage and visualized using a heatmap. Pathway enrichment analyses were carried out using GeneCodis 2.0.^18,19

Bayesian analysis to assess risk genes in stage Ib patients

We restricted our analysis to stage Ib patients only, since preliminary analysis suggested that gene-expression patterns associated with survival differ in different stages and stage Ib group was the largest group (48 patients altogether, but clinical data available for only 46 of them). From the survival times we identified two distinct patient groups, a group with short survival (<1000 days, n = 20) and a group with long survival (>1000 days, n = 26) (Fig. 2A). To reduce the number of genes, 500 gene-clusters or “metagenes”, were formed from the 5000 genes with highest variance among the 46 patients. Clustering of the genes is performed using complete-linkage hierarchical clustering,^20,21 where the distance between two genes is defined as the correlation coefficient between gene-expression values. The constructed metagenes are summarized by the mean of the genes in the cluster. Thus genes demonstrating similar variation between the subjects are summarized as a single “metagene”.

The association between the metagenes and the two groups was analyzed with a sparse Bayesian probit model for binary response variables.^22–24 According to the model, the probability of being in the short survival group is given by Φ ( ∑ j = 1 500 B j x i , j ) , where Φ (·) is the cumulative distribution function of a standard normal distribution, B _j is the coefficient of the linear effect of metagene j, and x _i,j is the value of metagene j for patient i for i = 1, …, 46. The model parameters (B _j ) were estimated using a Markov chain Monte-Carlo (MCMC) approach together with sparsity inducing priors. For the details of the choice of prior distributions and sampling, see Hoti and Sillanpää 2006. ²³ For each of the 500 metagenes, we calculated the posterior probability that the metagene is included in the model, where being included in the model is defined as having a normalized effect greater than 0.1. The normalized effect of metagene j is defined as B _j divided by the standard deviation of the values of metagene j. Thus, for metagenes with a higher variance, a smaller effect is allowed. Using the same inclusion criteria, we calculated the posterior distribution of the number of genes included in the model. To take into account the possibility that the chain converges to a local solution, we combined 20 independent MCMC runs into a single analysis. Each model was run for 30,000 rounds, and the first 20,000 rounds were removed as burn-in rounds. From the remaining 10,000 rounds, every 10th round was saved and used in the analysis. The ability of the model to predict the survival outcome of new patients was evaluated by leave-one-out cross-validation, where one patient at a time is used to validate the model (test person), while the remaining 45 patients are used to build up the model. For each of the 46 models, the gene selection and metagene construction was performed based on the 45 patients only, and the test person was assigned to the survival group with the higher posterior probability (corresponds to cut off value 0.5).

The performance of the method was further evaluated by the ROC curve, constructed by plotting the false positive rate (sensitivity) against the true positive rate (1-specificity) for different cut off values in the range (0.1). In addition to drawing the empirical curve, we report the AUC (area under the ROC curve) value together with the P-value corresponding to the event AUC > 0.5. Kaplan-Meier survival curves were drawn to visually compare the survival rates in the model-predicted high-and low-risk groups. The significance of the difference was evaluated with the log-rank test.

Quantitative RT-PCR

To validate the gene expression levels detected with microarray analysis, qRT-PCR was performed for the four up-regulated genes: TTK (Dual specificity protein kinase), CCNB2 (Cyclin B2), BUB1B (budding uninhibited by benzimidazoles 1 homolog beta (yeast)), PTTG2 (pituitary tumor-transforming 2), and four down-regulated genes: CLIC5 (Chloride intracellular channel protein 5), GPR116 (G-protein coupled receptor 116), AGER (RAGE) (advanced glycosylation end product-specific receptor), TNNC1 (troponin C type 1). We chose ESD (esterase D/formylglutathione hydrolase) as the endogenous reference for qRT-PCR because it has been previously identified as invariant in clinical lung cancer specimens.^25,26 The transcripts were amplified using Maxima SYBR Green /ROX qPCR Master mix (Fermentas) and ABI Prism 7900HT (Applied Biosystems). Eight sample pairs (tumor sample with corresponding normal lung) that were present in array and three sample pairs which were not (independent validation set), were used in the qRT-PCR experiment. Relative gene expression levels were calculated using the relative quantification method (Applied Biosystems). For detailed information see supplemental method and supplemental Table 2.

Identification of two gene expression patterns and correlation with NSCLC subtypes

We identified 599 genes which were down-regulated and 402 genes which were up-regulated in NSCLC compared to the normal lung tissue (Supplemental Table 3 and supplemental Fig. 1). According to Genecodis 2.0 analyses the main up-regulated processes in cancer were not only related to mitosis, cell division, DNA replication, blood vessel development, keratinozyte differentiation and epidermis development (Supplemental Table 4). The number of various down-regulated processes is much larger including immune response, signal transduction, cell to cell adhesion, cell surface reception linked signaling pathway, cell differentiation and others (Supplemental Table 5). To find differentially expressed genes between lung cancer subtypes, we used only tumor samples of different subtypes as sample groups. Because of the small sample size and similar histology ²⁷ of AC and BC samples we used these in analyses as one group. We identified 112 genes which were up-regulated and 101 genes which were down-regulated in AC/BC compared to the SCC (Supplemental Fig. 2). Some of the genes showing largest fold changes in our dataset are belonging to keratin gene family, which overexpression is shown to be SCC specific. ²⁸

We also carried out analyses to identify genes which may distinguish NSCLC stages. Because most of our sample group consisted of early stage (Ia and Ib) tumor samples and there were limited number of later stage samples, we compared only I stage tumor samples with all the others. We did not find any genes which were significantly overexpressed in one of the sample groups after multiple testing correction. This may suggest that although TNM staging system reflects the clinical status of tumor, it may not be the best tool to assess the underlying biological properties of cancer caused by aberrant gene expression.

The expression of the four up-regulated and four down-regulated genes (see supplemental data Table 2) were confirmed using eight pairs of normal and tumor tissue samples from microarray analysis sample set and further validated using three sample pairs that were not analyzed using microarray previously (Fig. 1A). The results were quantified using 2^-ΔΔCt method and ESD as an endogenous reference. Significant correlation was observed between log2–transformed array signal intensities and qRT-PCR Ct values, showing the consistency between qRT-PCR and microarray data (Fig. 1B).

OPEN IN VIEWER

Figure 1.

A) Validation of microarray data with qRT-PCR of four up-regulated and down-regulated genes. ▪, average log fold-change for paired lung cancer samples on microarray (n = 21), □, qRT-PCR average log fold-change for the lung cancer sample pairs that were presented on microarray (n = 8),▫, qRT-PCR average log fold-change for the lung cancer sample pairs that were not presented on the microarray (n = 3). Error bars indicate the standard error of the mean (SEM). B) Correlation between array log2(signaltumor)-log2(signalnormal) and qRT-PCR ΔΔCt for validated genes using same sample pairs as previous graph (n = 8). Pearson correlation coefficients (R), correlation test p-values, and best-fitting (least squares) lines are shown.

Among the up-regulated genes we validated, CCNB2 and TTK are involved in M phase of the cell cycle,^29–31 which indicate cell proliferation is involved in carcinogenesis. AGER (RAGE), another gene we validated, showed a statistically significant down-regulation in the gene expression values between NSCLC and normal lung tissues, which is consistent with the earlier findings that down-regulation of AGER (RAGE) supports NSCLC. ³²

Identifying genes associated with survival in stage Ib patients

To identify prognostic markers for high-risk patients with early-stage disease, the Bayesian model was first applied to both Ia and Ib patients (results not shown). However, the model-based prediction into low and high survival risk groups was clearly less accurate for Ia patients than for Ib patients. This suggests that gene-expression profiles associated with survival are different between stages Ia and Ib. A possible explanation for why the model worked for stage Ib but not for stage Ia patients could be that the number of patients in stage Ib (n = 46) was larger than that in stage Ia (n = 14).

Since the largest group of samples in our dataset were from patients with tumor stage Ib, we next focused to this subset of patients. A sparse Bayesian regression model was constructed using all 46 stage Ib patients. The number of influential metagenes supported by the model was 4 (mode) with a 95% credible interval of 1–10. The top two metagenes (MAGEA9B (melanoma antigen family A, 9B), MAGEA10 (melanoma antigen family A, 10), and NLRP2 (NLR family, pyrin domain containing 2)) and (CYP3A5 (cytochrome P450, family 3, subfamily A, polypeptide 5), AKR1B1 (aldo-keto reductase family 1, member B1), FAM46B (family with sequence similarity 46, member B)) (Fig. 2A) were included into the model with 11% and 10% posterior probabilities, respectively. For the remaining metagenes, the posterior probabilities were less than 5%. According to our model, genes within these metagenes are associated with a survival prognosis of stage Ib patients, and are therefore potential targets for further research. In the SCC group, high gene expression values of the first top metagene are associated with long survival (≥1000 days), while SCC group with low values are associated with short survival (≤1000 days) (Fig. 2B).

OPEN IN VIEWER

Figure 2.

A) Leave-one-out cross-validation results demonstrating the ability of the model to classify new Ib patients into groups of low and high survival risk. The y-axis gives the posterior probabilities that a patient is classified into the long survival group (survival time > 1000 days) and the x-axis gives the observed survival time (including censored times). AC and BC patients are denoted by BA and SCC patients by S. Patients above the 0.5 posterior probability line (green dotted line) are classified into the long survival (low risk) group and patients below the line are classified into the short survival (high risk) group. B) A heatmap presentation of the two top metagenes most frequently found to be associated with the two risk groups in the Bayesian model. The genes included into the metagenes are listed on the y-axis and the patients sorted according to their survival times (days) are listed on the x-axis. High gene-expression values are shown with red color and low values with green ().

Metagenes MAGEA9B and MAGEA10 have been involved in congenital dyskeratosis, ³³ and NLRP2 may be involved in another type of familial imprinting disorder. ³⁴ CYP3A5, located in the second top metagene, contains one of the most commonly detected polymorphisms that influences the efficacy of vinorelbine-based therapies to treat NSCLC. ³⁵ However, the exact role of these metagenes in tumorigenesis requires further studies.

Identifying stage specific survival-associated genes would require fitting separate models or adopting a model that allows state-specific gene-effects. For example, the binary classification-tree approach ³⁶ allows one to combine very different gene-expression profiles (leaves of the tree) into the same class. We used the cut-off value of 1000 days, since at that point there exists a clear gap (757 days to 1225 days) in the survival times of Ib patients. As shown in Fig. 2A, there are 3 patients who died after the cut-off that may have a negative effect on the quality of the data, ie, if the censored times may turn out to be much longer. In Potti et al 2006, the analysis was based on a cohort consisting of NSCLC patients with a recurrence within 2.5 years and patients with no recurrence within 5 years. Also, their analysis differs from ours in the sense that their model considers both clinical data and gene data.

Predicting patient survival for stage Ib patients

The ability of the Bayesian model to correctly classify stage Ib patients into groups of short and long survival based on gene-expression profiles was evaluated using leave-one-out cross-validation and the ROC curve. The use of leave-one-out cross-validation provides an estimate of the ability of the model to predict survival of new patients (Fig. 2A). The cross-validation classification error was 33% (15 out of 46 were misclassified) for all patients, 29% (5/17) for the short survival group (<1000 days), and 34% (10/29) for the long survival group (>1000 days).

Kaplan-Meier analysis of the model-predicted risk groups for Ib patients was performed. The difference between the survival rates in the high and low risk groups was significant (P = 0.0007, using the log-rank test) (Fig. 4). Further, according to the ROC analysis the overall performance of the classification method is positive AUC = 0.728 (P-value 0.001) (Fig. 3).

OPEN IN VIEWER

Figure 3.

The empirical ROC curve (solid curve): The true positive rate plotted as a function of the false positive rate for different cut off values. Jumps in the curve correspond to changes in the classification outcome of the patients due to the use of different cut off values. For cut off value 0 all patients are classified into the long survival group, FTP = 0 and TPR = 0. In the other extreme, ie, cut off value 1 all patients are classified into the short survival group, FTP = 1 and FPR = 1. The area under the curve (AUC = 0.728, P-value 0.001) is an overall measure of the quality of the classification method. The dashed line (AUC = 0.5) is the expected ROC curve for a totally random classifier.

OPEN IN VIEWER

Figure 4.

Kaplan-Meier plot of the survival probability in the high and low risk groups predicted by the Bayesian model. The high risk group (short survival) consisting of 24 patients and the low risk group (long survival) consisting of 22 patients. Vertical drops indicate deaths and ticks on the solid lines are censored survival times. The survival rates of the two groups are significantly different (P = 0.0007).