Four patient cohorts with published microarray transcriptional profiles were used in this study. The first cohort contains 442 lung adenocarcinoma patient samples from the Director’s Challenge Study (15). This study cohort is composed of 4 data sets (University of Michigan, H. Lee Moffitt Cancer Center, Memorial Sloan-Kettering Cancer Center and Dana-Farber Cancer Institute) contributed by 6 institutions. The clinical characteristics and smoking status of the patients are summarized in Table I.

The second patient cohort contains 130 squamous cell lung cancer samples from Raponi et al(16). The third cohort contains 111 NSCLC samples from Bild et al(17). The fourth cohort contains 164 airway epithelial cell lung tissue samples from current and former smokers published by Spira et al(2). This cohort has 60 lung cancer samples (48 NSCLC, 11 small cell lung cancer and 1 unknown histology) and 69 normal lung tissue samples. Patient gene expression profiles from Shedden et al(15), Raponi et al(16) and Spira et al(2) were quantified with Affymetrix HG-U133A. The dataset from Bild et al(17) was quantified with Affymetrix HG-U133 Plus 2. The raw microarray data were quantile-normalized and log2 transformed with dChip (18) for further analysis.

The implication induction algorithm (11) based on prediction logic (19) was used to derive coexpression between each pair of genes using software Genet (11,12,20). In the biological context, the six foremost implication rules relating two dichotomous variables are interpreted as follows: A ⇒ B: upregulation of gene A causes upregulation of gene B; A ⇒ ¬B: upregulation of gene A causes downregulation of gene B; ¬A ⇒ B: downregulation of gene A causes upregulation of gene B; ¬A ⇒ ¬B: downregulation of gene A causes downregulation of gene B; A ⇔ B: upregulation of gene A causes upregulation of gene B and upregulation of gene B causes upregulation of gene A; A ⇔ ¬B: upregulation of gene A causes downregulation of gene B and downregulation of gene B causes upregulation of gene A. Mean expression of each gene in the training set was used to define up- or downregulation. The minimum scope and the minimum precision of a derived implication relation were significantly greater than zero (P<0.05, one-sided z-tests).

The following pathway databases were used to evaluate the biological relevance of the derived coexpression networks, including NCBI Entrez Gene (21), Kyoto Encyclopedia of Genes and Genomes (KEGG) (22), NCI-Nature Pathway Interaction Database (http://pid.nci.nih.gov/), protein-protein interaction database STRING 8 (23), and Pathway Studio 7.0 (Ariadne Genomics, Rockville, MD, USA). In addition, five gene set collections [positional, curated, motif, computational and Gene Oncology (GO)] and canonical pathway databases from the MSigDB (24) were used in the network precision and FDR evaluation. Using these resources, a coexpression relation is considered a true positive (TP) if the pair of genes satisfy any of the following: i) on the same chromosome or cytogenetic band; ii) in the same curated or canonical pathway; iii) sharing a cis-regulator motif, binding motif, or transcription factor binding site; iv) annotated by the same GO term; v) having protein-protein interaction; or vi) within the same computational gene sets mined from cancer-oriented microarray data. The coexpression relation is considered a false positive (FP) if the gene pair do not satisfy all five conditions listed above (25). If at least one gene in the pair is not annotated, a coexpression relation is labeled as non-discriminatory (ND). Coexpression relations labeled as ND were excluded in this evaluation as they were not confirmed. Network precision is defined as: network _ precision = TP TP + FPThe portion of FP over all positive cases is defined as q-value: q – value = FP TP + FPThe FDR of the smoking-mediated coexpression networks was calculated by averaging the q-values obtained from the null distribution generated in 1,000 random permutations of the class labels in the test cohort.

The stability of the computationally derived smoking-mediated coexpression networks was evaluated using different subsets of patient samples from the training set in 100 iterations. The stability is defined as the portion of the smoking-mediated coexpression relations obtained from the original data that were retrieved by using only a random subset of the training data and the full test data.

NCI-H23 (ATCC no. CRL-5800) cells were cultured in RPMI-1640 medium (Mediatech, Manassas, VA, USA) supplemented with 10% FBS (Hyclone, Logan, UT, USA), 2 mM L-glutamine (Mediatech), 100 IU penicillin/ml (Sigma, St. Louis, MO, USA), and 100 μg streptomycin/ml (Sigma). BEAS-2B (ATCC no. CRL-9609) cells were cultured in Dulbecco’s modified Eagle’s medium (Mediatech) supplemented with 5% FBS (Hyclone), 2 mM L-glutamine (Mediatech), 100 IU penicillin/ml (Sigma) and 100 μg streptomycin/ml (Sigma).

H23 and BEAS-2B cells were treated with 100 nM NNK (Toronto Research Chemicals, North York, ON, Canada) for 15 min, 1 and 16 h. Four repeats (total of five samples) were performed on each cell line and for each time point. Following treatment, cells were harvested by trypsinization and protein was isolated. Cells were lysed in CLB lysis buffer (50 mM Tris-HCl, pH 7.4, 150 mM NaCl, 1% Triton X-100, 0.25% Na-deoxycholate, 5 mM EDTA and 1 mM NaF) supplemented with 1% (v/v) HALT Protease Inhibitor Cocktail, purchased from Thermo Scientific (Rockford, IL, USA), on ice for 15 min with occasional mixing by vortex. Lysates were centrifuged at 20,800 × g for 15 min to pellet insoluble debris and then supernatants were collected. Total protein concentration was determined by bicinchonic acid (BCA) protein assay purchased from Pierce Protein Research Products (Rockford, IL, USA).

Proteins (50 μg) were resolved on precast Mini-Protean TGX Gels (4–20%; Bio-Rad Laboratories, Hercules, CA, USA). After boiling for 5 min with reducing Laemmli buffer, proteins were separated and subsequently transferred to a PVDF membrane at 100 mV for 1 h at 4°C. After transfer, membranes were blocked in NET-gelatin solution (150 mM NaCl, 5 mM EDTA, 50 mM Tris-HCl, pH 7.5, 0.05% Triton X-100 and 0.25% gelatin) for 1 h at room temperature. Primary antibody was added to membranes in 15 ml NET-gelatin solution [1:500 dilution for anti-EGFR, 1:25,000 dilution for anti-GAPDH, 1:2,000 dilution for anti-p53, 1:1,000 dilution for phospho-p53 (phospho S15) and anti-Rb] and membranes were incubated for 2 h at room temperature with rocking. Membranes were then washed in NET-gelatin solution (3 x 20 min with shaking) with HRP-conjugated secondary monoclonal anti-mouse IgG antibody purchased from GE Healthcare UK Ltd. (Little Chalfont, UK). After 1 h of incubation, unbound secondary antibody was removed by washing in NET-gelatin solution (3 x 20 min with shaking). Signal was visualized using Immobilon chemiluminescent HRP substrate from Millipore (Billerica, MA, USA). Primary antibodies utilized included mouse monoclonal anti-EGFR, from Thermo Fisher Scientific (Fremont, CA, USA), and mouse monoclonal anti-GAPDH purchased from Fitzgerald Industries International Inc. (Acton, MA, USA). In addition, the following antibodies were used in western blotting: anti-p53 [Abcam, Mouse Monoclonal (ab26)], anti-phospho-p53/phospho S15 [Abcam, Rabbit Polyclonal (ab1431)] and anti-Rb [Abcam, Mouse Monoclonal (ab24)].

Relative EGFR, p53, phospho-p53, and Rb expression was determined by densitometric analysis using ImageJ software provided by NIH (http://rsb.info.nih.gov/ij/index.html). X-ray films were scanned at 300 and 600 DPI using a CanoScan (Canon, Lake Success, NY, USA) and images were imported into ImageJ for analysis. The raw signal intensity was determined by selecting the peak corresponding to each band and integrating the intensity within that peak. Local background intensity (calculated by averaging the background intensities at the upper and lower bounds of the peak) was integrated and subtracted from each raw intensity to give the background-corrected signal intensity. To account for loading differences, the corrected signal intensity for the assayed proteins was divided by the corrected GAPDH intensity. Samples treated with varying NNK exposure times were compared to untreated controls for the H23 and BEAS-2B cell lines.

Total-RNA was isolated from both cell lines using the mirVana™ miRNA Isolation kit and following the manufacturer’s protocol (Ambion, Austin, TX, USA). Total-RNA was eluted in 100 μl of nuclease-free water and stored at −80°C. RNA concentration was determined using the NanoDrop 1000 Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). RNA quality, 28S/18S ratio, and a visual image of the 28S and 18S bands were evaluated using the 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Total-RNA (1 μg) was converted into complementary DNA (cDNA) using the High Capacity cDNA Reverse Transcription Kit from Applied Biosystems (Life Technologies, Carlsbad, CA, USA). Thermal cycling conditions were as follows: 25°C for 10 min, 2 cycles of 37°C for 60 min and 85°C for 5 sec followed by a programmed hold at 4°C.

All qRT-PCR reactions were performed on a 7500 real-time PCR system from Applied Biosystems. The reports were generated using SDS2.3 software (Applied Biosystems). The Ct values obtained were normalized to the UBC housekeeping gene in each sample. Fold changes were computed using the 2^−ΔΔCt method of 5 biological replicates and 3 technical replicates (26). Statistical significance was computed using repeated ANOVA tests in R and is considered statistically significant at P≤0.05.

The coexpression relation of a gene pair derived with the implication induction algorithm was compared with the observed NNK-induced gene expression changes. The coexpression relation is confirmed when the observed NNK-induced gene expression changes are consistent with the predicted implication rule between the two genes. For example, if the rule between gene A and gene B is positive equivalence (A⇔B), it is confirmed if both gene A and gene B showed overexpression in the corresponding experiments.

Lung cancer survival genes were first selected from the whole genome on the training set (UM and HLM; n=256) from Shedden et al(15). A total of 2,310 genes were significantly associated with overall survival (P<0.05, univariate Cox model). Second, from this set of 2,310 genes, 217 genes exhibited significant differential expression (P<0.05, t-tests) in smokers versus non-smokers. These 217 survival and smoking-associated genes, as well as 6 major signaling pathway genes (EGF, EGFR, MET, KRAS, E2F3, and E2F5) were included in the network analysis. These signaling pathways are included in human NSCLC disease mechanisms delineated by the KEGG Pathway Database (http://www.genome.jp/kegg/pathway/hsa/hsa05223.html). They were selected based on their reported clinical relevance in NSCLC. These 6 signaling pathway genes were not significantly associated with survival nor were they differentially expressed in smokers.

Patient samples in the training set were separated into two groups: smokers (patients who smoked in the past or who are currently smoking) and non-smokers (patients who never smoked). For each smoking-defined patient group, a coexpression network among the 223 genes was constructed. Between each pair of the 223 genes, significant (P<0.05; z-tests) coexpression relations were retrieved in the smoker group and the non-smoker group, constituting smoking-mediated gene coexpression networks in NSCLC. By comparing the coexpression types between each pair of genes in the two networks, differential network components were identified and considered important for further evaluation. These differential components are interactions that were present in the smoker group but missing in the non-smoker group, or conversely, those present in the non-smoker group but absent in the smoker group. From the differential components associated with the smoker group and non-smoker group, genes having direct coexpression relations with all 6 lung cancer signaling pathway genes were identified as the signature genes (Fig. 1). As a result, 6 genes were identified from the smoker group and 1 gene was identified from the non-smoker group. This constituted the smoking-associated 7-gene signature for NSCLC (Table II).

We sought to investigate if the identified gene signature could provide accurate prognostic prediction of survival in lung adenocarcinoma patients. On the training cohort, the original microarray gene expression profiles of the identified 7 gene probes were fitted into a Cox model as covariates. A survival risk score was generated for each patient in the training set. A training model (Fig. 2A) was identified and applied to the test set (MSK and DFCI; n=186; Fig. 2B) without re-estimation of parameters. In both training and test sets, this scheme separated patients into two groups with different survival outcomes (P<0.007, Kaplan-Meier analyses). The hazard ratio of the 7-gene risk score [HR=1.89, 95% CI: (1.06, 3.38)] was higher than other lung cancer prognostic factors except cancer stage in the test set (Table III). There was no significant difference in prognostic value between the hazard ratio of the 7-gene risk score and cancer stage (II vs. I). The results demonstrate that the 7-gene risk score could provide a more accurate prognosis than some commonly used clinicopathological parameters.

The 7-gene signature gave accurate prognostic prediction in smokers in both training and test sets in Shedden’s cohorts (15) (P<0.01; Fig. 2C and D), but not in non-smokers in Kaplan-Meier analyses (P<0.12, results not shown). In addition, gene expression-defined high- and low-risk groups had significant association with smoking (P<0.02, χ² tests) and smoking cessation (P<0.00001, χ² tests; Table I). These results further confirmed the smoking association of the identified 7-gene signature.

The prognostic performance of the 7-gene signature was further evaluated on cohorts from Raponi et al(16) and Bild et al(17), which include another major subtype of NSCLC, squamous cell lung carcinoma. For robust validation, patient samples in these two studied cohorts were randomly partitioned into separate training and test sets. A prognostic classifier was constructed on the training set using the Cox model and validated on the test set without re-estimation of parameters.

In the Raponi cohort (16) of squamous cell carcinoma patients, the 7-gene signature stratified patients into two groups with distinct survival outcomes (log-rank P<0.005, Kaplan-Meier analysis) in the training set (Fig. 2E). This model generated borderline significant stratification (P=0.06, Kaplan-Meier analysis) in the test set (Fig. 2F). This might be owing to the fact that 10 patients (8%) of Raponi’s cohort were either non-smokers or their smoking status was not known, whereas the 7-gene signature provides refined prognosis specifically in smoking lung cancer patients.

In the Bild cohort (17) containing both lung adenocarcinoma and squamous cell carcinoma patients, the 7-gene signature stratified patients into two distinct survival groups in both training and test sets (P<0.04, Kaplan-Meier analyses) (Fig. 2G and H). Overall, these results demonstrate that the 7-gene signature could select high-risk NSCLC patients with a smoking history for chemotherapy.

We further investigated whether the 7-gene signature could be used for diagnostic screening of lung cancer in smokers. The smoking cohort from Spira et al(2) was separated into a training set (n=77) and two independent test sets (n=52 and n=35). Using a nearest neighbour algorithm implemented in WEKA (27), the 7-gene classifier could accurately identify lung cancer patients from normal patients with an overall accuracy of 73 and 74% in two test sets, respectively. The odds ratio of predicted lung cancer risk was highly significant in all three sets [OR=3.85, 95% CI: (1.45, 10.20), P<0.007 in training set; OR=7.35, 95% CI: (2.16, 25.04), P<0.001 in Test set 1; OR=8.45, 95% CI: (1.84, 38.75), P<0.006 in Test set 2; Table II]. These results indicate that the identified 7-gene signature has important implications in diagnostic screening of lung cancer risk in smokers.

BEAS-2B and H23 cells were treated with NNK for validation of smoking-associated gene expression. Each cell line was exposed to NNK for 15 min, 1 and 16 h. Ten signaling pathway genes and 7 signature genes were examined. Results showed that 9 genes (GPRC5C, LTF, SEMA3C, E2F1, E2F4, E2F5, EGF, EGFR and TP53) exhibited significant differential expression in the NNK-treated H23 cells (Fig. 3A). In BEASE-2B cells, all genes except CYP3A4 were expressed following NNK exposure, with 13 genes (GPRC5C, LTF, PIGN, SEMA3C, E2F1, E2F3, E2F5, EGF, EGFR, KRAS, MET, TP53 and RB1) exhibiting significant differential expression (Fig. 3B).

To further evaluate the NNK-induced protein expression, western blot assays were performed in NNK-treated BEAS-2B and H23 cells after 15 min, 1 and 16 h exposures. The results show that EGFR had consistent overexpression at both the mRNA and protein levels over the time course in BEAS-2B cells after NNK treatment. In H23 cells, EGFR exhibited over-expression at the mRNA level; however, protein expression was downregulated following NNK exposure (Fig. 3C and D). These results indicate that in normal lung epithelial cells, the EGFR gene is overexpressed upon NNK treatment, consistent with previous findings (28); whereas in lung adenocarcinoma cells, the NNK-induced transcriptional and translational regulation of EGFR are not concordant at the same time points.

p53 had consistent NNK-induced expression patterns at mRNA and protein levels, with short-term downregulation followed by upregulation in BEAS-2B cells. As H23 is p53 deficient, p53 protein was not expressed in these cells (Fig. 3E and F). Interestingly, downregulation of phospho-p53 was consistently observed in both NNK-treated BEAS-2B and H23 cells (Fig. 3G and H), concordant with its mRNA and total protein expression. These results are consistent with the report that NNK induces damage in the p53 gene (29).

In NNK-treated BEAS-2B cells, Rb gene expression was first significantly downregulated at 15 min, returned to its normal expression at 1 h, and then showed modest overexpression at 16 h; whereas the Rb protein showed a steady overexpression following the NNK treatment at all-time points. In NNK-treated H23 cells, both mRNA and protein expression of Rb were downregulated at all-time points (Fig. 3I and J). These results are consistent with the reported increased phosphorylation of the Rb Ser⁷⁹⁵ 6–15 h after NNK treatment of normal human bronchial epithelial cells (NHBE) and small airway epithelial cells (SAEC). This, in turn, promoted cells entering into the S phase (at 15–21 h) (30).

There were 17 gene coexpressions specifically associated with smokers, and one coexpression specifically associated with non-smokers significant (P<0.05) in both training and test sets from Shedden et al(15) (Fig. 4A). Among these 18 coexpression relations, 11 gene associations were confirmed with multiple biological databases (Fig. 4A; network precision=0.91; FDR<0.01), and the network structure was stable (Fig. 4B).

The smoking-associated gene coexpression network in lung adenocarcinoma patients (Fig. 4A) was further validated using the NNK-induced gene expression changes in lung adeno-carcinoma cell line H23 (Fig. 3A). Results show that most of the coexpressions in the smoking-associated network of lung adenocarcinoma tumors were confirmed by the coexpressions observed in NNK-treated H23 cells, at varying time points (Fig. 4C). All of the 17 smoking-associated gene coexpressions in NSCLC patients were observed in the cell experiments, except two coexpressions, one between EGF and LTF and one between CRTAC1 and GPRC5C. These two unobserved gene coexpressions could be related to other sources of carcinogens in tobacco, because NNK is only one of the carcinogens in tobacco, among about 54 others (14). Overall, the smoking-associated gene coexpression network in NSCLC patients was largely confirmed in NNK-treated cell experiments, elucidating a network of smoking-induced gene alterations in NSCLC.

This study presents novel implication network formalism for biomarker discovery. The ability to model cyclic relations in Genet overcomes the fundamental drawback of acyclic Bayesian networks in modelling molecular networks (31). In comparison with Bayesian belief networks, expression profiles of the identified 7 signature genes and 6 signalling pathway genes were used to build causal networks with TETRAD IV (http://www.phil.cmu.edu/projects/tetrad/current.html) for smoking and non-smoking lung adenocarcinoma patients in training and test sets, respectively (Fig. 2). There was only one interaction associated with smokers (between MET and SEMA3C; Fig. 2E) present in both training and test sets, which was considered a true positive when evaluated with MSigDB. In contrast, Genet generated significantly more biologically relevant gene coexpression relations that were validated by the biological databases and experimental results, confirming its topological advantage over the Bayesian belief networks.

Large-scale gene coexpression networks have been used in disease classification (32). These studies construct pair-wise gene coexpression networks by using correlation coefficients computed from gene expression profiles. Such networks indicate the distance or similarity between each pair of gene expression profiles but do not provide the direction or causal relations in the gene regulatory patterns. We have compared Genet with gene association networks based on Pearson’s correlation. In constructing smoking-mediated coexpression networks using 217 smoking and survival associated genes and 6 signalling hallmarks, both models had the same network precision and FDR. However, Genet generated significantly more biologically relevant gene association relations that were validated by the test set (20). These results indicate that prediction logic is more robust than Pearson’s correlation for inducting gene association networks.