[ Thoracic Oncology Original Research ] Independent Validation of Early-Stage Non-Small Cell Lung Cancer Prognostic Scores Incorporating Epigenetic and Transcriptional Biomarkers With Gene-Gene Interactions and Main Effects Ruyang Zhang, PhD; Chao Chen, BS; Xuesi Dong, MS; Sipeng Shen, PhD; Linjing Lai, MS; Jieyu He, MS; Dongfang You, BS; Lijuan Lin, MS; Ying Zhu, BS; Hui Huang, BS; Jiajin Chen, BS; Liangmin Wei, BS; Xin Chen, BS; Yi Li, PhD; Yichen Guo, PhD; Weiwei Duan, PhD; Liya Liu, PhD; Li Su, BS; Andrea Shafer, MPH; Thomas Fleischer, PhD; Maria Moksnes Bjaanæs, PhD; Anna Karlsson, PhD; Maria Planck, PhD; Rui Wang, PhD; Johan Staaf, PhD; Åslaug Helland, PhD; Manel Esteller, PhD; Yongyue Wei, PhD; Feng Chen, PhD; and David C. Christiani, MD, FCCP BACKGROUND: DNA methylation and gene expression are promising biomarkers of various cancers, including non-small cell lung cancer (NSCLC). Besides the main effects of biomarkers, the progression of complex diseases is also influenced by gene-gene (GG) interactions. RESEARCH QUESTION: Would screening the functional capacity of biomarkers on the basis of main effects or interactions, using multiomics data, improve the accuracy of cancer prognosis? STUDY DESIGN AND METHODS: Biomarker screening and model validation were used to construct and validate a prognostic prediction model. NSCLC prognosis-associated bio- markers were identified on the basis of either their main effects or interactions with two types of omics data. A prognostic score incorporating epigenetic and transcriptional biomarkers, as well as clinical information, was independently validated. RESULTS: Twenty-six pairs of biomarkers with GG interactions and two biomarkers with main effects were significantly associated with NSCLC survival. Compared with a model using clinical information only, the accuracy of the epigenetic and transcriptional biomarker-based prognostic model, measured by area under the receiver operating characteristic curve (AUC), increased by 35.38% (95% CI, 27.09%-42.17%; P ¼ 5.10  10–17) and 34.85% (95% CI, 26.33%-41.87%; P ¼ 2.52  10–18) for 3- and 5-year survival, respectively, which exhibited a superior predictive ability for NSCLC survival (AUC3 year, 0.88 [95% CI, 0.83-0.93]; and AUC5 year, 0.89 [95% CI, 0.83- 0.93]) in an independent Cancer Genome Atlas population. GG interactions contributed a 65.2% and 91.3% increase in prediction accuracy for 3- and 5-year survival, respectively. INTERPRETATION: The integration of epigenetic and transcriptional biomarkers with main effects and GG interactions significantly improves the accuracy of prognostic prediction of early-stage NSCLC survival. CHEST 2020; 158(2):808-819 KEY WORDS: early stage; interaction; multiomics; non-small cell lung cancer; prognostic score ABBREVIATIONS: AUC = area under the receiver operating charac- NSCLC = non-small cell lung cancer; QC = quality control; ROC = teristic curve; Cindex = concordance index; CPTAC = Clinical Proteo- receiver operating characteristic; SIS = sure independence screening; mic Tumor Analysis Consortium; GEO = Gene Expression Omnibus; TCGA = the Cancer Genome Atlas GG = gene-gene interactions; HR = hazard ratio; ISIS = iterative sure independence screening; KEGG = Kyoto Encyclopedia of AFFILIATIONS: From the Department of Biostatistics (Drs Zhang, Genes and Genomes; LUAD = lung adenocarcinoma; LUSC = lung Shen, Duan, Y. Wei, and F. Chen; Messrs C. Chen, You, Huang, J. squamous cell carcinoma; MGH = Massachusetts General Hospital; Chen, L. Wei, and X. Chen; and Mss Dong, Lai, He, Lin, and Zhu), 808 Original Research [ 158#2 CHEST AUGUST 2020 ] Lung cancer is a leading cause of cancer-related death disease, early-stage patients often have a considerably worldwide and was estimated to cause 1.76 million more favorable prognosis. However, significant deaths in 2018.1 The 5-year survival rate among patients heterogeneity in clinical prognosis is observed for with lung cancer remains relatively low, ranging from patients with early-stage non-small cell lung cancer 4% to 17% depending on clinical characteristics.2 (NSCLC) with similar clinical characteristics, which Compared with patients diagnosed with late-stage indicates the importance of understanding molecular mechanisms.3 Identifying molecular changes in oncogene and/or tumor suppressor genes that are Center for Global Health, School of Public Health, Nanjing Medical University, Nanjing, China; the Department of Environmental Health associated with NSCLC survival is helpful for (Drs Zhang, Shen, Guo, Y. Wei, and Christiani; Mss Dong and Lin; and developing targeted therapies to prolong patients’ Messrs You and Su), Harvard T. H. Chan School of Public Health, Boston, MA; the China International Cooperation Center for Envi- survival time. ronment and Human Health (Drs Zhang, Shen, Y. Wei, and F. Chen), Nanjing Medical University, Nanjing, China; the Department of DNA methylation is a heritable, reversible, epigenetic Medical Oncology (Drs Zhang and Wang), Jinling Hospital, School of Medicine, Nanjing University, Nanjing, China; the Department of modification that affects the spatial conformation of Epidemiology and Biostatistics (Ms Dong), School of Public Health, DNA and regulates gene expression.4,5 DNA Southeast University, Nanjing, China; the Department of Biostatistics methylation is a molecular biomarker and may be a (Dr Li), University of Michigan, Ann Arbor, MI; the Department of Biostatistics (Dr Guo), Harvard T. H. Chan School of Public Health, therapeutic target for the treatment of cancer.6,7 In Boston, MA; the Department of Bioinformatics (Dr Duan), School of addition, gene-gene (GG) interactions have long Biomedical Engineering and Informatics, Nanjing Medical University, Nanjing, China; the Department of Preventive Medicine (Dr Liu), been recognized to regulate the progression of Medical School of Ningbo University, Ningbo, China; the Pulmonary complex diseases, including NSCLC.8 The and Critical Care Division, Department of Medicine (Ms Shafer and Dr Christiani), Massachusetts General Hospital and Harvard Medical development of cancer may be related to interactions School, Boston, MA; the Department of Cancer Genetics (Drs between several key genes.9 Lung cancer prognosis- Fleischer, Moksnes Bjaanæs, and Helland), Institute for Cancer Research, Oslo University Hospital, Oslo, Norway; the Division of associated biomarkers have been proposed on the basis Oncology and Pathology, Department of Clinical Sciences (Drs of omics data, including DNA methylation,10 gene Karlsson, Planck, and Staaf), Lund and CREATE Health Strategic expression,11 microRNA,12 and long noncoding Center for Translational Cancer Research, Lund University, Lund, Sweden; the Institute of Clinical Medicine (Dr Helland), University of RNA.13 However, most studies are limited to a single Oslo, Oslo, Norway; the Josep Carreras Leukemia Research Institute type of omics data, which results in less accurate (Dr Esteller), Badalona, Barcelona, Spain; the Centro de Investigacion Biomedica en Red Cancer (Dr Esteller), Madrid, Spain; the Institucio prognostic models.14 For example, our previous Catalana de Recerca i Estudis Avançats (Dr Esteller), Barcelona, Spain; integrative omics study of the BTG2 gene showed that the Physiological Sciences Department (Dr Esteller), School of Medi- cine and Health Sciences, University of Barcelona, Barcelona, Spain; this gene could slightly improve the prediction the State Key Laboratory of Reproductive Medicine (Dr F. Chen), accuracy of early-stage NSCLC survival.6 However, a Nanjing Medical University, Nanjing, China; and the Jiangsu Key Laboratory of Cancer Biomarkers, Prevention and Treatment (Dr F. large-scale integrative analysis of multiomics data has Chen), Cancer Center, Collaborative Innovation Center for Cancer identified genes with either important main effects or Personalized Medicine, Nanjing Medical University, Nanjing, China. gene-gene (GG) interactions, based on which a more Dr Zhang, Mr C. Chen, and Ms Dong contributed equally to this work. accurate prognostic prediction model of NSCLC can Drs Zhang, Y. Wei, and F. Chen contributed equally to this work as corresponding authors. be constructed. Dr Christiani is the senior author who supervised the work. FUNDING/SUPPORT: This study was funded by the National Natural Specifically, we used a two-stage study design and Science Foundation of China [81530088 to F. C., 81973142 to Y. W.], the performed an integrative analysis of pan-cancer-related National Key Research and Development Program of China [2016YFE0204900 to F. C.], the Natural Science Foundation of genes to identify prognostic biomarkers with either a Jiangsu Province [BK20191354 to R. Z.], the Natural Science Foun- main effect or GG interactions using epigenome and dation of the Jiangsu Higher Education Institutions of China [18KJB310011 to R. Z.], the China Postdoctoral Science Foundation transcriptome data from multiple study centers. We [2018M633767 to R. Z.], the NIH [CA209414, CA092824, and then built a prognostic prediction model for early-stage ES000002 to D. C. C., AG056764 to Y. L.], and the Priority Academic NSCLC by incorporating both selected epigenetic and Program Development of Jiangsu Higher Education Institutions (PAPD). R. Z. was partially supported by the Outstanding Young transcriptional biomarkers. Teachers Training Program of Nanjing Medical University. CORRESPONDENCE TO: Feng Chen, PhD, SPH Bldg, Room 412, 101 Longmian Ave, Nanjing, Jiangsu 211166, China; e-mail: fengchen@njmu.edu.cn Methods Copyright Ó 2020 The Authors. Published by Elsevier Inc under li- cense from the American College of Chest Physicians. This is an open Only patients with early-stage (stage I or II) lung adenocarcinoma access article under the CC BY-NC-ND license (http:// (LUAD) and lung squamous cell carcinoma (LUSC) were included creativecommons.org/licenses/by-nc-nd/4.0/). in our study. DNA methylation data were harmonized from five DOI: https://doi.org/10.1016/j.chest.2020.01.048 international study centers, including Harvard, Spain, Norway, chestjournal.org 809 Sweden, and the Cancer Genome Atlas (TCGA). Gene expression data sequencing by expectation maximization. Level 3 gene quantification were composed of four datasets from the Gene Expression Omnibus data were downloaded from the TCGA data portal and were further (GEO) and TCGA. checked for quality. Gene probes were excluded if the missing rate > 80%, and the batch effect was corrected with ComBat. The Harvard: The Harvard cohort consisted of patients seen at Massachusetts expression value of each gene was transformed on a log2 scale and General Hospital (MGH), and histologically confirmed as having primary standardized before association analysis. NSCLC, recruited since 1992.15 We profiled 151 early-stage patients from this cohort. A lung pathologist at MGH evaluated each specimen for the DNA methylation and gene expression of 719 pan-cancer-related genes amount (tumor cellularity, > 70%) and quality of tumor cells. The were then used for subsequent association analysis. Gene symbols for specimens were classified histologically according to World Health the 719 pan-cancer-related genes were obtained from the Catalogue Organization criteria. The institutional review boards at the Harvard of Somatic Mutations in Cancer (COSMIC). After QC, there were T. H. Chan School of Public Health and MGH approved the study. All 12,806 CpG probes identified for association analysis. CpG probes patients provided written informed consent. from five genes (BTG2,6 KDM,7 EGLN2,8 LRRC3B,15 and SIPA1L320) Spain: The Spanish cohort included 226 patients with early-stage NSCLC reported in our previous study were also included. recruited from eight subcenters between 1991 and 2009.10 Patients provided written consent and tumors were surgically collected. This study was approved by the Bellvitge Biomedical Research Institute Statistical Analysis institutional review boards. The flow of analysis is depicted in Figure 1. Epigenetic and transcriptional analyses were performed simultaneously, and a Norway: The Norwegian cohort consisted of 133 patients with early- discovery phase and validation phase were used to identify NSCLC stage NSCLC from Oslo University Hospital, recruited between 2006 prognostic biomarkers. In each procedure, we conducted analysis of and 2011.16 The project was developed with the approval of the Oslo both the main effects and gene-gene interactions among biomarkers. University Institutional Review Board and regional ethics committee Patients having DNA methylation data from Harvard, Spain, (S-05307). All patients provided informed consent. Tumor tissues Norway, and Sweden, as well as patients having gene expression data were snap frozen in liquid nitrogen and stored at –80 C until DNA from GEO, were assigned to the discovery phase for epigenetic isolation. analysis and transcriptional analysis, respectively. Patients having Sweden: Tumor DNA was collected from 103 patients with early-stage two types of omics data from TCGA were assigned to the validation NSCLC, including 80 patients with LUAD and 23 patients with LUSC, phase. at the Skåne University Hospital in Lund, Sweden.17 The study was For the main effect analysis, we used sure independence screening (SIS) developed under the approval of the Regional Ethical Review Board and LASSO Cox penalized regression to screen biomarkers with main in Lund, Sweden (Registration nos. 2004/762 and 2008/702). effects that were relevant to survival, using the R package SIS. SIS TCGA: A total of 332 LUAD and 285 LUSC with full DNA methylation, LASSO is a two-stage procedure. At the first stage, SIS selects the survival time, and covariates data were included. Level 1 biomarkers with the strongest marginal associations with survival. At HumanMethylation450 DNA methylation data from patients with the second stage, LASSO was used to perform variable selection and early-stage NSCLC were downloaded on October 1, 2015. parameter estimation simultaneously among the biomarkers selected at the first stage. During the LASSO procedure, tuning parameter GEO: Transcriptome information from 425 patients with early-stage selection was based on Bayesian information criteria. To capture NSCLC was profiled using the Affymetrix Human Genome U133A biomarkers that might be missed at the first stage, we repeatedly Plus 2.0 Array (e-Table 1). Only data from patients with available applied the SIS LASSO algorithm to the remaining unselected survival time, clinical stage, and tumor tissue expression values were biomarkers until no new biomarkers can be recruited.21 This analyzed. iterative procedure is termed iterative SIS (ISIS) LASSO. To account for the biologic heterogeneity between LUAD and LUSC, we used a Quality Control for DNA Methylation Data histology-stratified multivariate Cox proportional hazards model to DNA methylation was assessed with Illumina Infinium test these biomarkers, using the R package survival. The stratified HumanMethylation450 BeadChips (Illumina Inc.). Raw image data model adjusted for the differences between LUAD and LUSC in were imported into GenomeStudio Methylation Module V1.8 baseline hazards. The other covariates adjusted in the model were (Illumina Inc.) to calculate methylation signals and to perform age, sex, study center, clinical stage, and smoking status. normalization, background subtraction, and quality control (QC). For the GG interaction analysis, a histology-stratified multivariate Unqualified probes were excluded if they fitted any one of the Cox proportional hazards model adjusted for the aforementioned following quality control criteria: (1) failed detection (P > .05) in covariates was applied to identify biomarkers with GG interactions. $ 5% samples; (2) coefficient of variance < 5%; (3) all samples were The P value thresholds for multiple testing were established by the methylated or all were unmethylated; (4) common single-nucleotide Bonferroni method, which set the significance level to .05 divided by polymorphisms located in probe sequence or in 10-bp flanking the number of tests. This way, the overall type I error would be regions; (5) cross-reactive probes18; or (6) data did not pass QC in controlled at the .05 level. In our study, the significance level of all centers. Samples with > 5% undetectable probes were excluded. GG interaction analysis of epigenetic and transcriptional Methylation signals were further processed for quantile biomarkers was defined as 6.10  10–10 ¼ 0.05/(12,806  12,805/2) normalization (betaqn function in R package minfi) as well as type I and 1.94  10–7 ¼ 0.05/(719  718/2), respectively. and II probe correction (BMIQ function in R package lumi). They were adjusted for batch effects (ComBat function in R package sva) Significant biomarkers observed in the discovery phase were further according to the best pipeline by a comparative study.19 Details of confirmed in the validation phase and were retained if the P value the QC process are described in e-Figure 1. was # .05 and there was consistent direction of the effect across two phases. We also performed a test of proportional hazards Quality Control for Gene Expression Data assumption for each significant biomarker. The hazard ratio (HR) The TCGA workgroup completed the mRNA sequencing data and 95% CI were described as per 1% level of DNA methylation or processing and QC. Raw counts were normalized using RNA- gene expression increment. Sensitivity analysis was performed to 810 Original Research [ 158#2 CHEST AUGUST 2020 ] Epigenetic analysis Transcriptional analysis Discovery phase Discovery phase • 613 early-stage NSCLC • 425 early-stage NSCLC • 12806 CpG probe • 719 gene probe Main effect analysis GxG interaction analysis Main effect analysis GxG interaction analysis Biomarker Screening 23 CpG probe 2495 interaction terms 13 gene probe 40 interaction terms • ISIS LASSO method • Histology-stratified Cox model • ISIS LASSO method • Histology-stratified Cox model • Bonferroni adjusted P < .05 • Bonferroni adjusted P < .05 Validation phase Validation phase • 617 early-stage NSCLC • 617 early-stage NSCLC • 1930 CpG probe • 79 gene probe 1 CpG probes 151 interation terms 1 gene probe 2 interaction terms • P < .05 • Interaction P < .05 • P < .05 • Interaction P < .05 • Consistent effect • Consistent effect • Consistent effect • Consistent effect 1 CpG probe 149 interaction terms 1 gene probe 2 interaction terms • Sensitivity analysis • Sensitivity analysis • Sensitivity analysis • Sensitivity analysis Discovery phase Candidate genes Discovery phase • Forward stepwise selection • 10 CpG probes • Forward stepwise selection • Sensitivity analysis • Sensitivity analysis Epigenetic score Integrative score Transcriptional score Model Validation • 25 pairs with interaction • Epigenetic score • 1 pairs with interaction • 1 probe with main effect • Transcriptic score • 1 probe with main effect Prognostic score • Integrative score • Clinical information Validation phase • Model validation in TCGA • Sensitivity analysis Figure 1 – Flow chart of study design and statistical analyses. In the epigenetic analysis, patients with lung adenocarcinoma and lung squamous cell carcinoma from the Harvard, Spain, Norway, and Sweden cohorts were used in the discovery phase for screening, whereas data from the Cancer Genome Atlas (TCGA) was used for validation. In transcriptional analysis, gene expression data from Gene Expression Omnibus and TCGA were used in the discovery phase and the validation phase, respectively. Both main effect and GG interaction analyses were performed. GG ¼ gene by gene; NSCLC ¼ non-small cell lung cancer. confirm these robustly significant biomarkers. Patients were excluded if Epigenetic and transcriptional scores were calculated on the their DNA methylation (logit2 transformed) or expression (log2 basis of a weighted linear combination of individual values of transformed) values were out of range, based on mean  3  SD. the DNA methylation and gene expression, with weights derived from the Cox model. Integrative scores were For those identified biomarkers, we applied a forward stepwise synthesized by epigenetic and transcriptional scores. Finally, regression strategy to build up a multibiomarker Cox proportional the prognostic score was defined as the linear combination of hazards model in the discovery phase, which was then validated clinical information and integrative score (see e-Appendix 1). in TCGA samples. In the forward stepwise regression, a likelihood ratio test was applied to test the main effect or GG interaction Kaplan-Meier survival curves adjusted for the covariates were drawn to of biomarkers if Pentry # .05 and Pelimination > .05. Sensitivity represent the survival difference among patients with different scores. analysis was also performed using two different thresholds: .10 We predicted 3- and 5-year overall survival of patients, using the and .15. nearest neighbor method for time-to-event data.22 The accuracy of chestjournal.org 811 the prediction is presented using a receiver operating characteristic expression, using the reverse-phase protein array from the TCGA (ROC) curve and was measured by area under the ROC curve database. Second, we performed differential expression analysis (AUC), computed by the R package survivalROC. The prediction between tumor and normal tissues, and further investigated the main accuracy was confirmed with an independent TCGA population in effects of genes and GG interactions between genes on LUAD the validation phase. The 95% CI and P value of the AUC survival, using the Clinical Proteomic Tumor Analysis Consortium improvement were calculated on the basis of 1,000-time bootstrap (CPTAC) database. Differential protein expression analysis was resampling. Stratification analysis of prognostic scores was carried performed with the R package limma, which generated a linear out within subgroups stratified by age, sex, smoking status, clinical model to estimate fold changes and SEs prior to empirical Bayes stage, and histology. The concordance index (Cindex), an average smoothing.24 Finally, Gene Ontology (GO) and Kyoto Encyclopedia accuracy of predictive survival across follow-up years, as well as the of Genes and Genomes (KEGG) pathway enrichment analysis was 95% CI, which ranges from 0.5 to 1.0, were calculated to estimate carried out with Metascape. Gene network analysis was conducted the predictive performance.23 A nomogram was generated with R with GeneMANIA,25 a plugin of the Cytoscape application. The package rms to facilitate application of our model. critical hubs, highly connected to nodes in a module, were defined as the highest connectivity degrees. We assessed the potential functions of the identified genes at the protein level by taking advantage of limited public resources. First, P values were two-sided. All statistical analyses were performed with R we evaluated the association between protein expression and gene version 3.5.1 (R Foundation), unless otherwise specified. Results transcriptional score (HRdiscovery ¼ 2.44 [95% CI, 1.78- After QC, 1,230 (Ndiscovery ¼ 613 and Nvalidation ¼ 617) 3.35]; P ¼ 2.79  10–8). The associations between survival patients with 12,806 CpG probes and 719 gene probes and each of these scores were independently confirmed in were included in the association analysis. The the validation phase when adjusted for covariates demographic and clinical information are described in (epigenetic score: HRvalidation ¼ 2.72 [95% CI, 2.31-3.20], e-Tables 2, 3. P ¼ 6.06  10–33; transcriptional score: HRvalidation ¼ 2.64 [95% CI, 1.73-4.04], P ¼ 7.51  10–6; integrative score: For the main effect analysis of DNA methylation and HRvalidation ¼ 2.72 [95% CI, 2.32-3.18], P ¼ 5.68  10–35; gene expression, 23 CpG probes (e-Tables 4-6) and 13 prognostic score: HRvalidation ¼ 2.72 [95% CI, 2.34-3.17], gene probes (e-Tables 7, 8) were selected by ISIS LASSO, P ¼ 5.04  10–38). respectively. However, only cg19286631TRIM27 was significantly associated with survival in both phases To evaluate the discriminative ability of these scores, (HRdiscovery ¼ 1.03 [95% CI, 1.01-1.05], P ¼ 1.43  10–2; samples in the validation phase were categorized into HRvalidation ¼ 1.03 [95% CI, 1.01-1.06], P ¼ 1.13  10–3) low-, medium-, and high-score groups based on the and remained significant in sensitivity analysis. Also, tertiles of epigenetic, transcriptional, integrative, and only one gene probe located in the NDRG1 gene prognostic scores, respectively. Compared with the remained significant in the validation phase epigenetic low-score group, the medium- and high- (HRdiscovery ¼ 1.41 [95% CI, 1.05-1.89], P ¼ 2.16  10–2; score groups had 4.39- and 21.24-fold mortality risk, HRvalidation ¼ 1.12 [95% CI, 1.01-1.42], P ¼ 4.33  10–2) respectively (HRMedium vs Low ¼ 4.39 [95% CI, 2.42- and sensitivity analysis. 7.99], P ¼ 1.22  10–6; HRHigh vs Low ¼ 21.24 [95% CI, 11.23-40.17], P ¼ 5.67  10–21) (Fig 2A). Patients with For the GG interaction analysis, we observed 2,495 a high transcriptional score had significantly worse and 40 GG interactions from epigenetic and survival (HRMedium vs Low ¼ 1.46 [95% CI, 0.92-2.33], transcriptional analysis, respectively, in the discovery P ¼ 1.04  10–1; HRHigh vs Low ¼ 2.26 [95% CI, 1.41- phase. Finally, 149 and 2 GG interactions were 3.60], P ¼ 6.52  10–4) (Fig 2B). The significant survival retained in the validation phase that were also significant difference was enhanced among patients with different in the sensitivity analysis (e-Tables 9-13). integrative scores (HRMedium vs Low ¼ 4.32 [95% CI, By forward stepwise regression analysis in the discovery 2.39-7.83], P ¼1.33  10–6; HRHigh vs Low ¼ 24.32 phase, we observed one CpG probe with a main effect [95% CI, 12.71-46.56], P ¼ 5.76  10–22) (Fig 2C). and 25 pairs of CpG probes with GG interactions in Moreover, when combined with clinical information, the multibiomarker model (e-Table 14), which was used including age, sex, study center, clinical stage, and to calculate the epigenetic score (e-Table 15) smoking status, the prognostic score significantly (HRdiscovery ¼ 2.71 [95% CI, 2.41-3.05]; P ¼ 1.15  10–61). discriminated NSCLC survival (HRMedium vs Low ¼ 7.32 One gene probe with a main effect and one pair of gene [95% CI, 3.50-15.33], P ¼ 1.29  10–7; HRHigh vs Low ¼ probes with a GG interaction were retained in the 28.85 [95% CI, 13.13-63.43], P ¼ 5.83  10–17) (Fig multibiomarker model and used to calculate the 2D). The discriminative ability of the prognostic score is 812 Original Research [ 158#2 CHEST AUGUST 2020 ] A B E Epigenetic score Transcriptional score Discriminative ability of prognostic score 100 HRM vs L = 04.39, P = 1.22e−06 100 HRM vs L = 01.46, P = 1.04e−01 97.6 97.6 HRH vs L = 21.24, P = 5.67e−21 HRH vs L = 02.26, P = 6.52e−04 91.7 80 80 84.1 Survival probability Survival probability 5-year survival probability 3-year survival probability 60 60 59.8 40 40 50.0 50.0 48.2 48.9 20 20 0 0 24.9 0 4 8 12 16 20 0 4 8 12 16 20 11.1 Time (years) Time (years) 0.0 0.0 Low Medium High Low Medium High 0 3 4 5 8 12 16 20 Survival years C D Level 1 Level 2 Level 3 Level 4 Level 5 Integrative score Prognostic score HRM vs L = 07.32, P = 1.29e−07 100 HRM vs L = 04.32, P = 1.33e−06 100 HRH vs L = 28.85, P = 5.83e−17 HRH vs L = 24.32, P = 5.76e−22 F 80 80 Survival probability Survival probability Median survival years 1.6 3.0 4.9 8.6 14.5 60 60 0 9.13 4.40 1.0 40 40 21.02 P = 2.53e-03 20 20 HR 50 P = 5.93e-06 66.09 0 0 P = 3.24e-10 0 4 8 12 16 20 0 4 8 12 16 20 100 Time (years) Time (years) P = 1.98e-17 Low Medium High Low Medium High Figure 2 – Estimated survival curves for patients grouped by various biomarker-based scores. A, Epigenetic score of DNA methylation. B, Tran- scriptional score of gene expression. C, Integrative score of DNA methylation and gene expression. D, Prognostic score of DNA methylation, gene expression, and clinical information. Patients were categorized into low-, medium-, and high-score groups by using the tertiles of each score as the cutoffs. E, Discriminative ability of the prognostic score. Results of 3- and 5-year survival rate, median survival time, and hazard ratio (HR) were compared across five groups, defined by using the quintiles of the prognostic score as the cutoffs. F, HR and P values were derived from the Cox proportional hazards model for patients with different quintile levels of the prognostic score. HRH vs L ¼ HRHigh vs Low; HRM vs L ¼ HRMedium vs Low. further illustrated by categorizing patients on the basis for 3- and 5-year survival, respectively, and exhibited a of the quintile level of the score. Figure 2E manifests an superior predictive ability for NSCLC survival (AUC3 ordering relation: patients in higher-quintile groups had year ¼ 0.88 [95% CI, 0.83-0.93]; AUC5 year ¼ 0.89 lower 3- and 5-year survival rates, as well as shorter [95% CI, 0.83-0.93]) (Fig 4). GG interactions median survival time. This indicates that patients with contributed an additional 65.2% for 3-year and higher mortality risks can be detected by using our 91.3% for the 5-year prediction accuracy increase. score system (HRLevel 5 vs 1 ¼ 66.09 [95% CI, 25.13- 173.80], P ¼ 1.98  10–17; HRLevel 4 vs 1 ¼ 21.02 In the sensitivity analysis, we reanalyzed the stepwise [95% CI, 8.13-54.31], P ¼ 3.24  10–10; HRLevel 3 vs 1 ¼ regression using two different thresholds (P ¼ .10 and 9.13 [95% CI, 3.51-23.78], P ¼ 5.93  10–6; HRLevel 2 .15) and found that the majority of the selected –3 biomarkers were the same as those in the original vs 1 ¼ 4.40 [95% CI, 1.68-11.53], P ¼ 2.53  10 ) (Fig 2F). The performance of the prognostic score was regression model (e-Table 16). We then recalculated further confirmed in the analysis stratified by covariates these scores, retested their associations with NSCLC (Fig 3). survival, and obtained similar results (e-Table 17). Meanwhile, the AUCs of our prognostic model using We then independently validated the predictive ability of different thresholds were comparable: 0.88P ¼ .05 these biomarkers. The model with only clinical vs 0.85P ¼ .10 vs 0.86P ¼ .15 for 3-year survival; 0.89P ¼ .05 information, as aforementioned, had very limited vs 0.83P ¼ .10 vs 0.86P ¼ .15 for 5-year survival (e-Figs 2 prediction ability (AUC3 year ¼ 0.65, AUC5 year ¼ 0.66). and 3). However, by adding biomarkers with either main effects or GG interactions, the AUCs significantly increased Moreover, we found that the effects of these four scores by 35.38% (95% CI, 27.09%-44.17%; P ¼ 5.10  10–17) did not differ significantly between patients with LUAD and 34.85% (95% CI, 26.33%-41.87%; P ¼ 2.52  10–18) and patients with LUSC (PEpigenetic score ¼ .6572; chestjournal.org 813 Subgroup N HR(95% CI) P Value Age ≤ 68 317 2.62(2.07−3.32) 1.80 × 10–15 Age > 68 291 3.59(2.69−4.48) 4.47 × 10–18 Male 252 3.25(2.33−4.53) 4.48 × 10–12 Female 356 2.92(2.32−3.67) 4.23 × 10–20 LUAD 328 2.86(2.22−3.69) 5.41 × 10–16 LUSC 280 3.20(2.44−4.20) 3.52 × 10–17 Stage I 388 3.50(2.70−4.57) 3.53 × 10–21 Stage II 220 2.41(1.84−3.16) 2.25 × 10–10 Never/former smoker 427 2.93(2.37−3.63) 7.54 × 10–23 Current smoker 163 3.16(2.14−4.66) 6.30 × 10–09 1.0 2.0 2.5 3.0 4.0 5.0 Figure 3 – Forest plots of results from stratification analysis of prognostic score. HR with 95% CI of the prognostic score on non-small cell lung cancer survival in various subgroups is stratified by clinical characteristics. LUAD ¼ lung adenocarcinoma; LUSC ¼ lung squamous cell carcinoma. See Figure 2 legend for expansion of other abbreviation. A B ROC for 3-year survival of NSCLC ROC for 5-year survival of NSCLC 1.0 1.0 Time–dependent True Positive Time–dependent True Positive 0.8 0.8 0.6 0.6 0.9 0.9 0.4 0.8 0.4 0.8 AUC AUC 0.7 0.7 0.6 0.6 0.2 0.2 0.5 0.5 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Time–dependent False Positive Time–dependent False Positive C : AUC1 = 0.65 C : AUC1 = 0.66 C + main effect of E : AUC2 = 0.67 C + main effect of E : AUC2 = 0.66 C + main effect of E + GxG interaction of E : AUC3 = 0.72 C + main effect of E + GxG interaction of E : AUC3 = 0.68 C + main effect of E + M + GxG interaction of E : AUC4 = 0.76 C + main effect of E + M + GxG interaction of E : AUC4 = 0.68 C + main effect of E + M + GxG interaction of M : AUC5 = 0.87 C + main effect of E + M + GxG interaction of M : AUC5 = 0.88 C + main effect of E + M + GxG interaction of E + M : AUC6 = 0.88 C + main effect of E + M + GxG interaction of E + M : AUC6 = 0.89 AUC6 vs AUC1: n35.38% (95% CI: 27.09%, 44.17%) P = 5.10e-17 AUC6 vs AUC1: n34.85% (95% CI: 26.33%, 41.87%) P = 2.52e-18 Figure 4 – Receiver operating characteristic curves for various predictive models using the clinical information (C), the main and interaction effects of DNA methylation (M), and gene expression (E). A, Three-year survival prediction. B, Five-year survival prediction. The AUC increase (%) was evaluated by comparing the model with that with only the clinical information. P values and 95% CIs were calculated by using 1,000 bootstrap samples. AUC ¼ area under the receiver operating characteristic curve; ROC ¼ receiver operating characteristic. See Figure 1 legend for expansion of other abbreviations. 814 Original Research [ 158#2 CHEST AUGUST 2020 ] PTranscriptional score ¼ .1823; PIntegrative score ¼ .5532; We conducted a two-stage integrative study of DNA PPrognostic score ¼ .9653) (e-Table 18). Our prognostic methylation and gene expression data from multiple model retained similar prediction ability in both the centers to propose a prognostic scoring method LUAD (AUC3 year ¼ 0.91, AUC5 year ¼ 0.89, and incorporating transomics biomarkers with main effects Cindex ¼ 0.82; 95% CI, 0.76-0.87) and LUSC (AUC3 and GG interactions. The prognostic score, which was year ¼ 0.85, AUC5 year ¼ 0.87, and Cindex ¼ 0.82; 95% CI, validated in an independent population, effectively 0.76-0.88) populations, indicating the usefulness of the discriminates survival outcomes for patients with early- selected biomarkers and their interactions in predicting stage NSCLC and significantly improves prediction the outcomes for patients with LUAD and patients with accuracy for their prognosis. LUSC (e-Fig 4). GG interactions are of interest because they provide To facilitate application of our prognostic prediction important clues regarding the biologic mechanisms of model, we combined clinical information and scores of complex diseases.27 It was suggested in previous studies biomarkers and developed a nomogram, which that identification of GG interactions would improve estimated well a patient’s 3- or 5-year survival (e-Fig 5). the predictive accuracy of statistical models.28,29 The Cindex of the prognostic score indicated acceptable However, interactions might not dramatically improve prediction accuracy (Cindex ¼ 0.82; 95% CI, 0.78-0.86) in prediction if their effects are weak or there are few an independent TCGA population. The calibration plots significant interactions, but might optimize statistical also showed good accordance between observed and modeling.30 Besides prediction, GG interactions could predicted survival time (e-Fig 6). increase the power to detect associations and then be leveraged for the identification of new biomarkers.27 Our In protein analysis, three of the four genes mapped in results showed that biomarkers with GG interactions TCGA had significant correlation between gene significantly and predominantly improved the expression and protein expression (e-Table 19). Most prognostic prediction accuracy of early-stage NSCLC, (77%) of the 47 genes mapped in CPTAC were which might be due to increased power. differentially expressed between tumor and normal tissue, with statistical significance (e-Fig 7). In addition, one To evaluate the prediction accuracy of our model, we gene with main effect and four pairs of genes with GG conducted a literature search to compare our studies interaction had a significant effect on LUAD survival (e- with others. The details of the literature screening Table 20). Among 49 genes identified in epigenetic process are summarized in e-Figure 8. The prediction analysis, five genes (FOXP1, AFF3, BCL6, MAPK1, and accuracy of our model is superior (e-Table 22), as the STAT3) were identified as hub genes with the highest one study with the best AUC (0.80) had a very small connectivity degrees, greater than 25 (Fig 5). These 49 sample size. Another study with the largest sample size, genes were enriched in cancer-related pathways including without independent validation, had unsatisfactory the non-small cell lung cancer pathway (e-Table 21). prediction capacity (Cindex ¼ 0.64). Our study has a relatively large sample size and provides a satisfactory prediction model that performed well in an independent Discussion population regardless of AUC (AUC3 year ¼ 0.88 and An accurate prognostic predictive model may aid AUC5 year ¼ 0.89) and Cindex (0.82). physicians in making clinical decisions or guiding adjuvant therapy, especially for the vulnerable patients Among the genes identified in transcriptional analysis, with high mortality risk. Although subject and tumor NDRG131 and RHOA32 have been reported to be characteristics have been commonly used as valid associated with lung cancer. In this study, among 49 predictors, increased evidence has indicated that genes identified in epigenetic analysis, five (AFF3, molecular biomarkers may provide early warning MAPK1, STAT3, FOXP1, and BCL6) were identified as signals. This is because tumor cells may metastasize even hub genes in the gene network. AFF3 is associated with when the tumor size is undetectable (< 0.01 cm3) and NSCLC prognosis.33 MAPK1 promotes NSCLC cell aberrations of biomarkers occur.26 Thus, there is added survival and is a therapeutic target for NSCLC value when a prognostic predictive model incorporates chemotherapeutic resistance.34 STAT3, one of the three both genetic and nongenetic factors, whose effects can major downstream pathways activated by EGFR be captured using approaches that are both biologically phosphorylation,35 is persistently activated in 22% to stable and technically reproducible. 65% of NSCLC.36-38 It is a strong predictor of poor chestjournal.org 815 A C Transcriptional misregulation in cancer Regulation of neuron differentiation Regulation of hemopoiesis Negtive regulation of cell differentiation Ossification Pathways in cancer Stem cell differentiation Positive regulation of DNA-binding transcription factor Non-small cell lung cancer Skeletal system development Tissue morphogenesis Neuron projection morphogenesis Positive regulation of protein kinase activity Stem cell population maintenance Cytokine production Endoderm development Myeloid cell differentiation Muscle structure development Peripheral nervous system development Head development B hsa05202: Transcriptional misregulation in cancer GO:0045664: Regulation of neuron differentiation GO:1903706: Regulation of hemopoiesis GO:0045596: Negative regulation of cell differentiation GO:0001503: Ossification hsa05200: Pathways in cancer GO:0048863: Stem cell differentiation GO:0051091: Positive regulation of DNA-binding transcription factor hsa05223: Non-small cell lung cancer GO:0001501: Skeletal system development GO:0048729: Tissue morphogenesis GO:0048812: Neuron projection morphogenesis GO:0045860: Positive regulation of protein kinase activity GO:0019827: Stem cell population maintenance GO:0001816: Cytokine production GO:0007492: Endoderm development GO:0030099: Myeloid cell differentiation GO:0061061: Muscle structure development GO:0007422: Peripheral nervous system development GO:0060322: Head development 0 2 4 6 8 10 12 –log10(P) Figure 5 – Gene network and gene enrichment analysis of 49 genes to which 25 pairs of CpG probes with interaction and one CpG probe with main effect are mapped. A, The gene network plot constructed by GeneMANIA. Central nodes with boldface outline represent hub genes, and the size represents the connectivity degree of each node. B, Barplot of gene pathways enriched with significant genes, and colored by P values. C, The pathway network plot of these pathways enriched with significant genes. Significant pathways with a similarity > 0.3 are connected by edges. Each node represents an enriched term and is colored by its cluster identification. The size of the node represents the number of genes in the pathway. The edge represents potential biologic relationships between two pathways. GO ¼ Gene Ontology. NSCLC prognosis and related to cisplatin resistance in enriched in the KEGG non-small cell lung cancer NSCLC cells.39-41 FOXP1 is an independent factor for pathway (hsa05223). The hub genes MAPK1 and STAT3 predicting poor NSCLC prognosis.42 BCL6 could inhibit in the network were also involved in this pathway. The cell apoptosis in lung cancer43 and plays a role in results indicated that, after functional confirmation, the sustaining NSCLC genomic instability.44 identified CpG probes are potential epigenetic targets for NSCLC chemotherapy. In enrichment pathway analysis, 49 genes were significantly enriched in pathways or processes that are Our study has some strengths, as follows: (1) Most cancer related. Notably, the identified genes were also studies focus only on main effects of biomarkers, 816 Original Research [ 158#2 CHEST AUGUST 2020 ] ignoring their GG interactions that account for or excellent computer hardware can finish GG missing heritability of complex diseases like NSCLC.15 interactions for time-to-even data on genome-wide Also, most studies focus on single omics data testing scales within weeks. We exhaustively tested all pairs of prognostic biomarkers.10-13 Taking advantage of pan-cancer-related genes; (2) limited clinical epigenomic and transcriptomic data and considering information was available for several cohorts that were both GG interactions and main effects,45 we built up initiated decades ago. However, in our study, a few transomics prognostic scores, which could improve easily accessible clinical predictors and dozens of prognostic value; (2) to identify reliable prognostic biomarkers exhibited considerable accuracy, which biomarkers for the prediction of early-stage NSCLC indicated potentiality for real-world application; (3) overall survival, we used stringent statistical criteria. In the event rate of survival time for TCGA population is the main effect analysis, candidate biomarkers, with relatively low (23%), which considerably reduced the effect sizes larger than a data-driven threshold in ISIS statistical power. However, through a conservative LASSO, must reach statistical significance to stay in two-stage strategy this study showed the robustness of the Cox regression model. For the GG interaction our findings; (4) our prognostic prediction model analysis, we applied the most conservative Bonferroni predicts survival outcome and distinguishes subgroups correction to control for false positives. In addition, of patients with high mortality risk accurately, which significant biomarkers observed in the discovery phase provides a potential opportunity for the delivery of must be further validated in an independent personalized medicine and interventions tailored to population. However, one consequence was that only a each individual’s level of risk. However, it requires few biomarkers were identified because of the limited information on 54 biomarkers, which might restrict its sample size of gene expression data, which therefore clinical translatability without testing of specimens. contributed a small proportion of improved accuracy Nevertheless, the history of cancer omics testing has of our model; (3) we used ISIS LASSO and stepwise taught us that, as technology improves and costs fall, regression to screen biomarkers with main effect and the trend is toward more convenient and interactions, respectively, and built multibiomarker comprehensive approaches to quickly capture models. These coefficients, used as weights to define biomarker information.49 In the coming years, scores, were derived from multibiomarker models advances in technology will facilitate our model’s rather than single-biomarker models. Single- usefulness through a customized biochip to enable biomarker models might result in biased estimates of widespread clinical application and maximize the effect sizes, whereas multibiomarker models are more benefit to patients; and (5) further studies with a large- beneficial to clarify the complex association and could scale population and extension of other ethics are improve prediction accuracy46,47; (4) the prediction warranted to confirm the results of our association accuracy of our prognostic model was robust toward study and verify the underline biologic mechanisms of different selection thresholds in stepwise regression as the genes and their interactions. Results of protein well as stratification by histology types; and (5) the analysis in public resources and in gene network and genes we identified as enriched in the non-small cell enrichment analyses might provide insight into the lung cancer pathway and most of the hub genes have functional mechanisms. been reported to be associated with NSCLC, indicating the reliability of our prognostic biomarkers. We also acknowledge some limitations of our study, as Conclusion follows: (1) We focused only on pan-cancer genes, The prognostic score incorporating transomics whereas most dysregulated genes represent the biomarkers with both main effects and GG consequences rather than the causes of neoplastic interactions significantly improves prognostic prediction process.48 Moreover, few powerful statistical methods accuracy for early-stage NSCLC survival. chestjournal.org 817 Acknowledgments prospects for epigenetic therapy. Nature. microarray. Epigenetics. 2013;8(2):203- 2004;429(6990):457-463. 209. Author contributions: R. Z. takes responsibility for the content of the 5. Feinberg AP, Tycko B. The history of 19. Marabita F, Almgren M, Lindholm ME, manuscript, including the data and analysis. cancer epigenetics. Nat Rev Cancer. et al. An evaluation of analysis pipelines 2004;4(2):143-153. for DNA methylation profiling using the R. Z., C. C., X. D., Y. W., F. C., and D. C. C. Illumina HumanMethylation450 contributed to the study design. R. Z., L. S., A. 6. Shen S, Zhang R, Guo Y, et al. A multi- BeadChip platform. Epigenetics. 2013;8(3): S., T. F., M. M. J., A. K., M. P., J. S., Å. H., M. omic study reveals BTG2 as a reliable 333-346. E., and D. C. C. contributed to data prognostic marker for early-stage non- collection. R. Z., C. C., X. D., S. S., L. Lai, Y. small cell lung cancer. Mol Oncol. 2018;12: 20. Zhang R, Lai L, Dong X, et al. SIPA1L3 L., Y. W., F. C., and D. C. C. performed 913-924. methylation modifies the benefit of statistical analysis and interpretation and smoking cessation on lung 7. Wei Y, Liang J, Zhang R, et al. Epigenetic adenocarcinoma survival: an epigenomic- drafted the manuscript. J. H., D. Y., L. Lin, Y. modifications in KDM lysine smoking interaction analysis. Mol Oncol. Z., H. H., J. C., L. W., X. C., Y. G., W. D., L. demethylases associate with survival of 2019;13(5):1235-1248. Liu, L. S., and R. W. revised the manuscript. early-stage NSCLC. Clin Epigenetics. All authors contributed to critical revision of 2018;10(1):41. 21. Saldana DF, Feng Y. SIS: an R package for the manuscript and approved its final sure independence screening in ultrahigh 8. Zhang R, Lai L, He J, et al. EGLN2 DNA dimensional statistical models. J Stat version. Financial support and study methylation and expression interact with supervision were provided by F. C. and Software. 2018;83(2):1-25. HIF1A to affect survival of early-stage D. C. C. NSCLC. Epigenetics. 2019;14(2):118-129. 22. Heagerty PJ, Lumley T, Pepe MS. Time- dependent ROC curves for censored Financial/nonfinancial disclosures: None 9. Lin Z, Hui L, Yufei H, et al. Cancer survival data and a diagnostic marker. declared. progression prediction using Gene Biometrics. 2000;56(2):337-344. Interaction Regularized Elastic Net. IEEE/ Role of sponsors: The sponsor had no role in 23. Brentnall AR, Cuzick J. Use of the ACM Trans Comput Biol Bioinform. the design of the study, the collection and 2017;14(1):145-154. concordance index for predictors of analysis of the data, or the preparation of the censored survival data. Stat Methods Med manuscript. 10. Sandoval J, Mendez-Gonzalez J, Nadal E, Res. 2018;27(8):2359-2373. et al. A prognostic DNA methylation Availability of data and materials: COSMIC signature for stage I non-small-cell lung 24. Ritchie ME, Phipson B, Wu D, et al. (https://cancer.sanger.ac.uk/cosmic/), which cancer. J Clin Oncol. 2013;31(32):4140- limma powers differential expression provided the symbol of 719 pan-cancer- 4147. analyses for RNA-sequencing and related genes; TCGA (https://tcga-data.nci. microarray studies. Nucleic Acids Res. 11. Shedden K, Taylor JM, Enkemann SA, 2015;43(7):e47. nih.gov), now hosted at GDC (https://portal. et al. Gene expression-based survival gdc.cancer.gov/); GEO (https://www.ncbi. prediction in lung adenocarcinoma: a 25. Warde-Farley D, Donaldson SL, Comes O, nlm.nih.gov/geo/), for which detailed et al. The GeneMANIA prediction server: multi-site, blinded validation study. Nat information is provided in e-Table 1; CPTAC biological network integration for gene Med. 2008;14(8):822-827. (https://cptac-data-portal.georgetown.edu/), prioritization and predicting gene which provided protein data; Metascape 12. Tan X, Qin W, Zhang L, et al. A 5- function. Nucleic Acids Res. 2010;38(Web (http://metascape.org), which was used to microRNA signature for lung squamous Server issue):W214-W220. perform pathway enrichment analysis. cell carcinoma diagnosis and hsa-miR-31 26. Hu Z, Ding J, Ma Z, et al. Quantitative for prognosis. Clin Cancer Res. evidence for early metastatic seeding in Publisher’s note: Springer Nature remains 2011;17(21):6802-6811. colorectal cancer. Nat Genet. 2019;51(7): neutral regarding jurisdictional claims in 13. Zhou M, Guo M, He D, et al. A potential 1113-1122. published maps and institutional affiliations. signature of eight long non-coding RNAs 27. Cordell HJ. Detecting gene-gene Other contributions: The authors thank predicts survival in patients with non- interactions that underlie human diseases. TCGA, GEO, and CPTAC for contributing small cell lung cancer. J Transl Med. Nat Rev Genet. 2009;10(6):392-404. clinical, DNA methylation, RNA sequencing, 2015;13(1):231. 28. Pharoah PD, Antoniou AC, Easton DF, and protein data; COSMIC for providing 14. Zhao Q, Shi X, Xie Y, et al. Combining et al. Polygenes, risk prediction, and pan-cancer-related genes; as well as all of the multidimensional genomic measurements targeted prevention of breast cancer. study subjects who participated in this study. for predicting cancer prognosis: N Engl J Med. 2008;358(26):2796-2803. observations from TCGA. Brief Bioinform. Additional information: The e-Appendix, 2015;16(2):291-303. 29. Khoury MJ, Yang Q, Gwinn M, Little J, e-Figures, and e-Tables can be found in the Dana Flanders W. An epidemiologic Supplemental Materials section of the online 15. Guo Y, Zhang R, Shen S, et al. DNA assessment of genomic profiling for article. Methylation of LRRC3B: a biomarker for measuring susceptibility to common survival of early-stage non-small cell lung diseases and targeting interventions. Genet cancer patients. Cancer Epidemiol Med. 2004;6(1):38-47. References Biomarkers Prev. 2018;27(12):1527-1535. 30. Aschard H, Chen J, Cornelis MC, 1. Bray F, Ferlay J, Soerjomataram I, et al. 16. Bjaanæs MM, Fleischer T, Halvorsen AR, Chibnik LB, Karlson EW, Kraft P. Global cancer statistics 2018: et al. Genome-wide DNA methylation Inclusion of gene-gene and gene- GLOBOCAN estimates of incidence and analyses in lung adenocarcinomas: environment interactions unlikely to mortality worldwide for 36 cancers in 185 association with EGFR, KRAS and TP53 dramatically improve risk prediction for countries. CA Cancer J Clin. 2018;68:394- mutation status, gene expression and complex diseases. Am J Hum Genet. 424. prognosis. Mol Oncol. 2016;10(2):330-343. 2012;90(6):962-972. 2. Hirsch FR, Scagliotti GV, Mulshine JL, 17. Karlsson A, Jonsson M, Lauss M, et al. 31. Azuma K, Kawahara A, Hattori S, et al. et al. Lung cancer: current therapies and Genome-wide DNA methylation analysis NDRG1/Cap43/Drg-1 may predict tumor new targeted treatments. Lancet. of lung carcinoma reveals one angiogenesis and poor outcome in 2017;389(10066):299-311. neuroendocrine and four adenocarcinoma patients with lung cancer. J Thorac Oncol. 3. Tang S, Pan Y, Wang Y, et al. Genome- epitypes associated with patient outcome. 2012;7(5):779-789. wide association study of survival in early- Clin Cancer Res. 2014;20(23):6127-6140. 32. Konstantinidou G, Ramadori G, Torti F, stage non-small cell lung cancer. Ann Surg 18. Chen Y-a, Lemire M, Choufani S, et al. et al. RHOA-FAK is a required signaling Oncol. 2015;22(2):630-635. Discovery of cross-reactive probes and axis for the maintenance of KRAS-driven 4. Egger G, Liang G, Aparicio A, et al. polymorphic CpGs in the Illumina lung adenocarcinomas. Cancer Discov. Epigenetics in human disease and Infinium HumanMethylation450 2013;3(4):444-457. 818 Original Research [ 158#2 CHEST AUGUST 2020 ] 33. Zhang DL, Qu LW, Ma L, et al. Genome- clinicopathologic features in lung 44. Marullo R, Ahn H, Cardenas M, wide identification of transcription factors adenocarcinoma and squamous cell Melnick A, Xue F, Cerchietti L. The that are critical to non-small cell lung carcinoma. Mol Diagn Ther. 2011;15(6): transcription factor BCL6 is a rational cancer. Cancer Lett. 2018;434:132-143. 347-352. target in non-small cell lung cancer (NSCLC) [abstract]. Cancer Res. 34. Vicent S, Lopez-Picazo JM, Toledo G, 39. Barre B, Vigneron A, Perkins N, et al. The 2016;76(14 suppl):Abstract 1271. et al. ERK1/2 is activated in non-small-cell STAT3 oncogene as a predictive marker of lung cancer and associated with advanced drug resistance. Trends Mol Med. 45. Ng SW, Mitchell A, Kennedy JA, et al. tumours. Br J Cancer. 2004;90(5):1047- 2007;13(1):4-11. A 17-gene stemness score for rapid 1052. determination of risk in acute leukaemia. 40. Ikuta K, Takemura K, Kihara M, et al. Nature. 2016;540(7633):433-437. 35. Mitsudomi T, Yatabe Y. Mutations of the Overexpression of constitutive signal epidermal growth factor receptor gene and transducer and activator of transcription 3 46. Kang J, Kugathasan S, Georges M, et al. related genes as determinants of epidermal mRNA in cisplatin-resistant human non- Improved risk prediction for Crohn’s growth factor receptor tyrosine kinase small cell lung cancer cells. Oncol Rep. disease with a multi-locus approach. Hum inhibitors sensitivity in lung cancer. 2005;13(2):217-222. Mol Genet. 2011;20(12):2435-2442. Cancer Sci. 2007;98(12):1817-1824. 41. Harada D, Takigawa N, Kiura K. The role 47. Segura V, Vilhjalmsson BJ, Platt A, et al. 36. Zimmer S, Kahl P, Buhl TM, et al. An efficient multi-locus mixed-model of STAT3 in non-small cell lung cancer. Epidermal growth factor receptor approach for genome-wide association Cancers (Basel). 2014;6(2):708-722. mutations in non-small cell lung cancer studies in structured populations. Nat influence downstream Akt, MAPK and 42. Feng J, Zhang X, Zhu H, et al. High Genet. 2012;44(7):825-830. Stat3 signaling. J Cancer Res Clin Oncol. expression of FoxP1 is associated with 48. Vogelstein B, Papadopoulos N, 2009;135(5):723-730. improved survival in patients with non- Velculescu VE, et al. Cancer genome 37. Looyenga BD, Hutchings D, Cherni I, small cell lung cancer. Am J Clin Pathol. landscapes. Science. 2013;339(6127):1546- et al. STAT3 is activated by JAK2 2012;138(2):230-235. 1558. independent of key oncogenic driver 43. Sun C, Li S, Yang C, et al. MicroRNA-187- 49. Kuo FC, Mar BG, Lindsley RC, et al. The mutations in non-small cell lung 3p mitigates non-small cell lung cancer relative utilities of genome-wide, gene carcinoma. PLoS One. 2012;7(2):e30820. (NSCLC) development through down- panel, and individual gene sequencing in 38. Jiang R, Jin Z, Liu Z, et al. Correlation of regulation of BCL6. Biochem Biophys Res clinical practice. Blood. 2017;130(4):433- activated STAT3 expression with Commun. 2016;471(1):82-88. 439. chestjournal.org 819