18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Epigenome-wide three-way interaction study identifies a complex pattern between TRIM27, KIAA0226, and smoking associated with overall survival of early-stage NSCLC Xinyu Ji1, Lijuan Lin1, Juanjuan Fan1, Yi Li2, Yongyue Wei1,3,4, Sipeng Shen1, Li Su3, Andrea Shafer5, Maria Moksnes Bjaanæs6 , Anna Karlsson7, Maria Planck7, Johan Staaf7, Åslaug Helland6,8, Manel Esteller9,10,11,12 , Ruyang Zhang1,3,4 , Feng Chen1,4,13,14 and David C. Christiani3,5 1 Department of Biostatistics, Center for Global Health, School of Public Health, Nanjing Medical University, Nanjing, China 2 Department of Biostatistics, University of Michigan, Ann Arbor, MI, USA 3 Department of Environmental Health, Harvard T.H. Chan School of Public Health, Boston, MA, USA 4 China International Cooperation Center for Environment and Human Health, Nanjing Medical University, Nanjing, China 5 Pulmonary and Critical Care Division, Department of Medicine, Massachusetts General Hospital and Harvard Medical School, Boston, MA, USA 6 Department of Cancer Genetics, Institute for Cancer Research, Oslo University Hospital, Oslo, Norway 7 Division of Oncology, Department of Clinical Sciences Lund and CREATE Health Strategic Center for Translational Cancer Research, Lund University, Lund, Sweden 8 Institute of Clinical Medicine, University of Oslo, Oslo, Norway 9 Josep Carreras Leukaemia Research Institute, Barcelona, Spain 10 Centro de Investigacion Biomedica en Red Cancer, Madrid, Spain 11 Institucio Catalana de Recerca i Estudis Avançats, Barcelona, Spain 12 Physiological Sciences Department, School of Medicine and Health Sciences, University of Barcelona, Barcelona, Spain 13 State Key Laboratory of Reproductive Medicine, Nanjing Medical University, Nanjing, China 14 Jiangsu Key Lab of Cancer Biomarkers, Prevention and Treatment, Cancer Center, Collaborative Innovation Center for Cancer Personalized Medicine, Nanjing Medical University, Nanjing, China Keywords The interaction between DNA methylation of tripartite motif containing 27 DNA methylation; non-small-cell lung (cg05293407TRIM27) and smoking has previously been identified to reveal his- cancer; overall survival; prognostic tologically heterogeneous effects of TRIM27 DNA methylation on early- prediction; three-way interaction stage non-small-cell lung cancer (NSCLC) survival. However, to understand Correspondence the complex mechanisms underlying NSCLC progression, we searched three- R. Zhang, Department of Biostatistics, way interactions. A two-phase study was adopted to identify three-way inter- Center for Global Health, School of Public actions in the form of pack-year of smoking (number of cigarettes smoked Health, Nanjing Medical University, SPH per day × number of years smoked) × cg05293407TRIM27 × epigenome-wide Building Room 406, 101 Longmian Avenue, DNA methylation CpG probe. Two CpG probes were identified with FDR- Nanjing, Jiangsu 211166, China q ≤ 0.05 in the discovery phase and P ≤ 0.05 in the validation phase: Tel: +86-025-86868436 cg00060500KIAA0226 and cg17479956EXT2. Compared to a prediction model E-mail: zhangruyang@njmu.edu.cn and with only clinical information, the model added 42 significant three-way F. Chen, Department of Biostatistics, Center interactions using a looser criterion (discovery: FDR-q ≤ 0.10, validation: for Global Health, School of Public Health, P ≤ 0.05) had substantially improved the area under the receiver operating Nanjing Medical University, SPH Building characteristic curve (AUC) of the prognostic prediction model for both Room 412, 101 Longmian Avenue, Nanjing, Abbreviations AUC, area under the receiver operating characteristic curve; CI, confidence interval; C-index, concordance index; DEGs, differentially expressed genes; FDR, false discovery rate; GO, gene ontology; HR, hazard ratio; HSD, honest significant difference; iTWINS, integrated covariates and Three-Way INteraction Score; KEGG, Kyoto Encyclopedia of Genes and Genomes; LUAD, lung adenocarcinomas; LUSC, lung squamous cell carcinomas; NSCLC, non-small-cell lung cancer; QC, quality control; ROC, receiver operating characteristic; SD, standard deviation; SNP, single nucleotide polymorphisms; TCGA, The Cancer Genome Atlas; TIICs, tumor-infiltrating immune cells; TIME, tumor immune microenvironment. Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of 717 Federation of European Biochemical Societies. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Three-way interaction study of NSCLC X. Ji et al. Jiangsu 211166, China 3-year and 5-year survival. Our research identified the complex interaction Tel: +86-025-86868414 effects among multiple environment and epigenetic factors, and provided E-mail: fengchen@njmu.edu.cn therapeutic target for NSCLC patients. Xinyu Ji, Lijuan Lin and Juanjuan Fan contributed equally to this work. David C. Christiani Senior author who supervised the work. (Received 3 October 2021, revised 23 November 2021, accepted 20 December 2021, available online 7 January 2022) doi:10.1002/1878-0261.13167 diseases, such as NSCLC, often involves a multistep, 1. Introduction multigenic, and multicausal biological process [25]; Lung cancer, as one of the most commonly diagnosed interactions between two relevant factors may only pro- cancer, is the leading cause of cancer mortality world- vide limited information understanding the cancer pro- wide, with an estimated 1.8 million deaths every year cess driven by multiple factors [26]. Therefore, high- [1]. More than 85% of lung cancer cases are non- order interactions may be needed to characterize com- small-cell lung cancer (NSCLC), of which lung adeno- plex association patterns relating to cancer recurrence carcinoma (LUAD) and lung squamous cell carcinoma and survival [27,28], including NSCLC [29,30]. This is (LUSC) are the two major histological subtypes [2,3]. because some well-confirmed low-order interactions, Despite recent advances in early diagnosis and ther- such as cg05293407TRIM27 × smoking, may also be sub- apy, the 5-year survival rate is low, ranging from 4 to jected to the influence of a third epigenetic factor, which 17%, which depends on clinical characteristics [4]. necessitates the study of three-way interactions arising Even with the same histological type and at the same from confirmed two-way interactions [30]. To our clinical stage, there exists wide heterogeneity in overall knowledge, epigenetic analysis with three-way interac- survival because of varying therapeutic responses [5]. tion has seldom been conducted because of the complex However, molecular mechanisms underlying these ther- nature of such analysis. apeutic responses remain largely unclear [6]. Hereby, we performed a two-phase epigenome-wide DNA methylation is a reversible epigenetic modifi- study of three-way interactions on overall survival of cation [7,8], with aberrant methylation being an early early-stage NSCLC based on the previously confirmed precursor of cancer progression [9] and a possible ther- two-way interaction between cg05293407TRIM27 and apeutic target for NSCLC [10–17] and other cancers smoking, using patients from international consortium [18,19]. Furthermore, epigenetic features are modifiable and the Cancer Genome Atlas (TCGA) database. by various environmental exposures (e.g., cigarette smoking), and changes contribute to the development 2. Materials and methods and progression of cancer as well. Cigarette smoking is one of the leading risk factors to lung cancer morbid- ity and mortality [20,21]. Our previous study found 2.1. Study populations that a two-way interaction between cg05293407TRIM27 and smoking reveals the histologically heterogeneous Included in the study were early-stage (stages I and II) effect of TRIM27 DNA methylation on survival NSCLC patients, with harmonized DNA methylation among early-stage NSCLC patients [22]. data and gene expression data, from five international As is well known, gene–gene and gene–environment study centers, including Harvard [31], Spain [32], Nor- interactions provide important clues regarding biologic way [33], and Sweden [34], and TCGA. All of these mechanisms of complex diseases [23] and are associated studies were approved by each Institutional Review with the overall survival of NSCLC [10,11,14,22,24]. Board, and patients provided written informed con- However, tumorigenesis and progression of complex sent. The study methodologies conformed to the 718 Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License X. Ji et al. Three-way interaction study of NSCLC standards set by the Declaration of Helsinki and was methylation data were downloaded at October 1, 2015. approved by the local ethics committee. In addition, 221 LUAD and 235 LUSC patients had gene expression data. The TCGA workgroup com- pleted processing and quality control (QC) of mRNA 2.1.1. Harvard sequencing data. Raw counts were normalized by The Harvard Lung Cancer Study cohort was described RNA-seq expectation maximization (RSEM). Level-3 previously [31]. All patients were recruited at Mas- gene quantification data were downloaded to further sachusetts General Hospital since 1992. They were newly check quality. All of the gene expression value was diagnosed and histologically confirmed as primary transformed on a log2 scale and standardized before NSCLC at the time of recruitment. Snap-frozen tumor association analysis. samples were taken from patients during complete resec- tion of the therapeutic procedure. Total 149 early-stage 2.2. Quality control procedures for DNA patients, who had complete survival information, were methylation data enrolled in our study. Tumor DNA was extracted from 5-μm-thick histopathologic sections. Each specimen was DNA methylation was profiled using Infinium evaluated by a pathologist for amount (tumor cellularity HumanMethylation450 BeadChips (Illumina Inc., San- > 70%) and quality of tumor cells. Among them, 26 Diego, CA, USA). All studies followed the same QC patients had mRNA sequencing data. procedures prior to association study. GenomeStudio Methylation Module V1.8 (Illumina Inc.) was used to convert raw image data into beta values (continuous 2.1.2. Spain numbers ranging from 0 to 100%) for background As described previously [32], tumors were collected by subtraction and control normalization. Unqualified surgical resection from 207 patients. DNA extraction probes meeting any one of the following criteria were was performed on tumor specimens (tumor cellularity excluded: (a) failed detection (P > 0.05) in more than > 50%). However, Spanish study center did not profile 5% samples; (b) coefficient of variance (CV) < 5%; (c) mRNA sequencing data. methylated or unmethylated in all samples; (d) com- mon single nucleotide polymorphisms (SNP) located in the probe sequence or 10-bp flanking regions; (e) 2.1.3. Norway cross-reactive probes or cross-hybridizing probes [35]; The Norwegian study consisted of 132 LUAD patients or (f) did not pass quality control in all centers. Sam- with operable lung cancer tumors seen at Oslo Univer- ples with > 5% undetectable probes or missing of sity Hospital, Rikshospitalet, Norway, in 2006–2011 pack-year of smoking were excluded. Methylation sig- [33]. Tumor tissues collected during surgery were snap- nals were further processed for quantile normalization frozen in liquid nitrogen and stored at −80 °C until (betaqn function in R package minfi) as well as type I DNA isolation. None of the enrolled patients received and II probes correction (BMIQ function in R pack- chemotherapy or radiotherapy before surgery. In addi- age lumi), and were further adjusted for batch effects tion, 93 of these patients had mRNA sequencing data. (ComBat function in R package sva) according to the best pipeline by a comparative study [36]. The QC pro- cesses are detailed in Fig. S1. 2.1.4. Sweden Tumor tissue samples were collected from 36 patients 2.3. Statistical analysis with early-stage NSCLC, at Skane University Hospi- tal, Lund, Sweden [34]. The study was developed under the approval of the Regional Ethical Review 2.3.1. A two-phase study of three-way interaction Board in Lund, Sweden (Registration no. 2004/762 and 2008/702). Meanwhile, mRNA sequencing data The analysis flow was depicted in Fig. 1, showing a were collected from 34 of these patients. two-phase study to detect three-way interactions that affect overall survival for early-stage NSCLC on an epigenome-wide scale. In the discovery phase, a 2.1.5. TCGA histology-stratified Cox proportional hazards model TCGA was composed of 227 LUAD and 241 LUSC adjusted for age, sex, clinical stage, and study centers. cases with complete overall survival time and covari- The analysis flow was used to identify CpG probes that ates. Level-1 HumanMethylation450 DNA interact with pack-year of smoking ×cg05293407TRIM27. Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of 719 Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Three-way interaction study of NSCLC X. Ji et al. Fig. 1. Flow chart of study design and statistical analyses. Patients with lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC) from the Harvard, Spain, Norway, and Sweden cohorts were used in the discovery phase for screening, whereas data from the Cancer Genome Atlas (TCGA) were used for validation. Subgroup analyses were based on quartiles of iTWINS. That is, we tried to find three-way interactions in the aforementioned multivariable Cox proportional haz- form of pack-year of smoking × cg05293407TRIM27 × ards model. Specifically, CpG probe, by using samples from Harvard, Spain, Norway, and Sweden. Hazard ratios (HR) and 95% iTWINS ¼ ∑αi  Covariatei þ β1  smoking confidence intervals (CI) were computed for per 1% þ β2  cg05293407TRIM 27 þ β3  CpG level of methylation increment. Multiple testing correc- tions were performed by controlling the false discovery þ β12  smoking  cg05293407TRIM 27 rate (FDR) at a 5% level [37]. In the validation phase, the identified three-way interactions from the discovery þ β13  smoking  CpG þ β23 phase were further tested by using samples from TCGA.  cg05293407TRIM 27  CpG þ β123 Significant interactions were those meeting the following criteria: (a) FDR-q ≤ 0.05 in the discovery phase and  smoking  cg05293407TRIM 27  CpG P ≤ 0.05 in the validation phase; and (b) with consistent effect directions across two phases. Kaplan–Meier curves where αs and βs were the coefficients estimated in the were used to compare survival differences between Cox model. patients with high or low methylation level of each iden- We conducted the following analyses. (a) We tified CpG probe, under different levels of smoking divided patients into four subgroups by the quartiles intensity and cg05293407TRIM27 methylation. of iTWINS and compared their survival by Kaplan– Meier curves. (b) To explore the molecular characteris- tics of these subgroups, we preformed gene expression 2.3.2. iTWINS subgroup analysis differential analysis through ANOVA (aov function in An integrated covariates and Three-Way INteraction R package multcomp) and Tukey’s honest significant Score (iTWINS) was defined by a weighted linear com- difference (HSD) test (TukeyHSD function in R pack- bination of all components inside the significant three- age multcomp). Genes with FDR-q ≤ 0.05 and HSD- way interaction, with weights derived from the q ≤ 0.05 were regarded as differentially expressed 720 Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License X. Ji et al. Three-way interaction study of NSCLC genes (DEGs). (c) We performed functional annotation remained significant (P ≤ 0.05) in the validation phase and gene enrichment pathway analyses of Gene Ontol- and showed robust association in the combined data ogy (GO) and Kyoto Encyclopedia of Genes and Gen- (Table 2). Two CpG probes, together with pack-year omes (KEGG) on the DEGs by using the of smoking and cg05293407TRIM27, had significant three- WebGestaltR package in R [38–40]. FDR-q ≤ 0.05 was way interaction effects on NSCLC survival, including considered statistically significant. (d) To explore the cg00060500KIAA0226 (HRinteraction = 0.993, 95% CI: pattern of tumor immune microenvironment (TIME) 0.990-0.996, P = 7.79 × 10-6, FDR-q = 0.039 in the dis- among subgroups, we computed the immune score covery phase; HRinteraction = 0.992, 95% CI: 0.986-0.999, representing the level of infiltrating immune cells based P = 0.021 in the validation phase; HRinteraction = 0.993, on gene expression data by using the ESTIMATE 95% CI: 0.991-0.996, P = 8.91 × 10−7 in the combined algorithm [41]. (e) We used CIBERSORT, a deconvo- data) and cg17479956EXT2 (HRinteraction = 0.997, 95% lution algorithm based on linear support vector regres- CI: 0.996-0.998, P = 1.16 × 10−5, FDR-q = 0.046 in the sion, to quantify the compositions of 22 types of discovery phase; HRinteraction = 0.993, 95% CI: 0.987- tumor-infiltrating immune cells (TIICs) [42]. 0.998, P = 0.011 in the validation phase; HRinteraction = 0.997, 95% CI: 0.996-0.999, P = 5.72 × 10-5 in the com- bined data); however, these two probes did not have sig- 2.3.3. Development of a prognostic prediction model nificant marginal effects (Table S2) and were We further selected three-way interactions by using a looser significantly associated with the expression of their cor- criterion (discovery: FDR-q ≤ 0.10, validation: P ≤ 0.05). responding genes (KIAA0226 and EXT2), respectively Then, they were included in a prognostic prediction model (Fig. S2). of NSCLC together with cg05293407TRIM27 × pack-year of To visualize the three-way interactions, we generated smoking. The 3- and 5-year overall survival of patients a 3D figure illustrating, for example, the varying effects were estimated for each subgroup by using the Kaplan– of cg00060500KIAA0226 at different levels of pack-year of Meier method [43]. The prediction accuracy was illustrated smoking and cg05293407TRIM27 (Fig. 2A). As the smok- using a receiver operating characteristic (ROC) curve and ing intensity increased and methylation level of was measured by area under the ROC curve (AUC) by cg05293407TRIM27 decreased, we observed an enhanced the R package survivalROC. The 95% CI and P value of risk associated with cg00060500KIAA0226; the risk of the AUC increase were calculated by 1000-time bootstrap cg00060500KIAA0226 was elevated as smoking intensity resampling. The concordance index (C-index), an average decreased and methylation level of cg05293407TRIM27 accuracy of predictive survival across follow-up years, increased. which ranges from 0.5 to 1.0, were calculated to estimate To more clearly illustrate the effect modifiers, we the predictive performance. compared the effect of cg00060500KIAA0226 (high vs Continuous variables were expressed as mean  s- low defined by median) on survival in four subgroups, tandard deviation (SD), and categorical variables were defined by the median of pack-year of smoking and expressed in frequency (n) and proportion (%). Statis- cg05293407TRIM27 (Fig. 2B-E). For patients with low tical analysis was performed using R version 3.6.1 level of pack-year of smoking, the effect of (The R Foundation of Statistical Computing). cg00060500KIAA0226 on survival was harmful (HRHigh vs Low = 1.70, 95% CI: 1.02-2.83, P = 0.042) when cg05293407TRIM27 was high (Fig. 2B), and was protec- 3. Results tive (HRHigh vs Low = 0.34, 95% CI: 0.17-0.66, After quality control, our epigenome-wide DNA methyla- P = 1.37 × 10-3) when cg05293407TRIM27 was low tion data were composed of 311,891 CpG probes from (Fig. 2C); for patients with high level of pack-year of 992 early-stage NSCLC patients; 524 patients (NLUAD = smoking, cg00060500KIAA0226 had a protective effect 425 and NLUSC = 99) were in the discovery phase and (HRHigh vs Low = 0.50, 95% CI: 0.29-0.86, P = 0.012) 468 patients (NLUAD = 227 and NLUSC = 241) were in the when cg05293407TRIM27 was high (Fig. 2D), and a validation phase. Demographic and clinical information detrimental effect (HRHigh vs Low = 1.65, 95% CI: for these patients were detailed in Table 1 and Table S1. 1.02-2.68, P = 0.041) when cg05293407TRIM27 was low (Fig. 2E). Similar patterns were also observed for the three-way interaction of pack-year of smok- 3.1. Two significant three-way interactions ing × cg05293407TRIM27 × cg17479956EXT2 (Fig. S3). identified in the two-phase study As demonstrated in Fig. S4A, these two three-way In the discovery phase, 91 three-way interactions interaction modules were highly correlated (r = 0.82, (FDR-q ≤ 0.05) were identified, while only 2 of them P = 2.66 × 10-241), and KIAA0226 and EXT2 gene Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of 721 Federation of European Biochemical Societies. 722 Table 1. Demographic and clinical descriptions for early-stage NSCLC patients with DNA methylation data in five international study centers. Discovery Validation Combined Variable USA (N = 149) Spaina (N = 207) Norway (N = 132) Sweden (N = 36) All (N = 524) TCGA (N = 468) Overall (N = 992) Age (years) 67.78  9.92 65.77  10.66 65.4  9.28 72.25  7.22 66.7  10.04 66.52  9.3 66.61  9.69 Pack-year of smoking 51.49  40.62 43.24  31.72 26.93  17.62 22.64  25.19 40.06  33.00 47.05  28.40 43.36  31.10 Sex Female 66 (44.30%) 99 (47.83%) 70 (53.03%) 22 (61.11%) 257 (49.05%) 185 (39.53%) 442 (44.56%) Male 83 (55.70%) 108 (52.17%) 62 (46.97%) 14 (38.89%) 267 (50.95%) 283 (60.47%) 550 (55.44%) Three-way interaction study of NSCLC Smoking status Never 18 (12.08%) 28 (13.53%) 16 (12.12%) 18 (50.00%) 80 (15.27%) 0 (0.00%) 80 (8.06%) Former 81 (54.36%) 113 (54.59%) 74 (56.06%) 11 (30.56%) 279 (53.24%) 323 (69.02%) 602 (60.69%) Current 50 (33.56%) 66 (31.88%) 42 (31.82%) 7 (19.44%) 165 (31.49%) 145 (30.98%) 310 (31.25%) Clinical stage I 102 (68.46%) 167 (80.68%) 92 (69.70%) 34 (94.44%) 395 (75.38%) 301 (64.32%) 696 (70.16%) II 47 (31.54%) 40 (19.32%) 40 (30.30%) 2 (5.56%) 129 (24.62%) 167 (35.68%) 296 (29.84%) Histology LUAD 96 (64.43%) 169 (81.64%) 132 (100.00%) 28 (77.78%) 425 (81.11%) 227 (48.50%) 652 (65.73%) LUSC 53 (35.57%) 38 (18.36%) 0 (0.00%) 8 (22.22%) 99 (18.89%) 241 (51.50%) 340 (34.27%) Chemotherapy No 140 (93.96%) 166 (91.21%) 101 (76.52%) 25 (89.29%) 432 (87.98%) 151 (75.88%) 583 (84.49%) Yes 9 (6.04%) 16 (8.79%) 31 (23.48%) 3 (10.71%) 59 (12.02%) 48 (24.12%) 107 (15.51%) Unknown 0 25 0 8 33 269 302 Radiotherapy No 130 (87.25%) 172 (94.51%) 131 (99.24%) 28 (100.00%) 461 (93.89%) 190 (95.48%) 651 (94.35%) Yes 19 (12.75%) 10 (5.49%) 1 (0.76%) 0 (0.00%) 30 (6.11%) 9 (4.52%) 39 (5.65%) Unknown 0 25 0 8 33 269 302 Adjuvant therapyb No 125 (83.89%) 157 (86.26%) 100 (75.76%) 25 (89.29%) 407 (82.89%) 146 (73.37%) 553 (80.14%) Yes 24 (16.11%) 25 (13.74%) 32 (24.24%) 3 (10.71%) 84 (17.11%) 53 (26.63%) 137 (19.86%) Unknown 0 25 0 8 33 269 302 Survival yearc Median (95% CI) 6.61 (5.14–7.49) 3.83 (3.21–4.46) 5.40 (5.16–5.77) 3.25 (2.09–4.39) 5.01 (4.53–5.25) 0.58 (0.50–0.70) 2.31 (2.00–2.64) Censoring rate 18.79% 54.59% 68.94% 52.78% 47.90% 76.07% 61.19% TCGA, The Cancer Genome Atlas; 95% CI, 95% confidence interval; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma. a Spain center is a collaborative study center, containing samples from Spain, Italy, the UK, France, and the USA. b Adjuvant therapy included chemotherapy and/or radiotherapy. c Restricted mean survival time was given because median was not available; proportion of samples lost to follow-up or alive at the end of study. X. Ji et al. Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License X. Ji et al. Three-way interaction study of NSCLC expression were also correlated significantly (r = 0.61, 2.02E-03 6.09E-04 1.04E-07 1.30E-03 1.14E-07 6.62E-07 8.91E-07 5.23E-04 4.12E-05 1.74E-05 4.18E-05 4.69E-06 8.99E-05 5.72E-05 P = 7.71 × 10-52) (Fig. S4B). Therefore, the subsequent P analysis focused on the interaction involving cg00060500KIAA0226 that had a lower P value. 0.996 1.007 0.999 1..09 0.49 0.83 0.93 1.45 1.02 0.61 0.87 0.97 1.22 1.04 3.2. Significant heterogeneity among iTWINS 95% CI subgroups 0.991 1.003 0.996 0.04 0.51 0.86 1.09 1.04 1.01 0.17 0.68 0.93 1.07 1.02 Combined data We categorized patients into four subgroups based on the quartiles of iTWINS, with Groups 1 and 4 representing the lowest and highest risk patients, respectively. As shown 0.993 1.005 0.997 0.14 0.65 0.90 1.26 1.07 1.01 0.32 0.77 0.95 1.14 1.03 HR in Fig. 3A, patients with higher iTWINS had a shorter median survival year (0.97, 95% CI: 0.75–1.25) and a higher risk of mortality (HRGroup4 vs 1 = 3.98, 95% CI: 0.396 0.142 0.012 0.271 0.016 0.019 0.021 0.146 0.034 0.019 0.031 0.015 0.017 0.011 2.98–5.33, P = 1.74 × 10-20; HRGroup3 vs 1 = 2.23, 95% P CI: 1.67–2.98, P = 5.13 × 10-8; HRGroup2 vs 1 = 1.39, 95% CI: 1.03–1.87, P = 0.033) (Fig. 3B). For comparison, we 0.999 1.019 0.998 6.94 1.19 0.97 1.85 1.14 1.03 1.76 0.97 0.99 1.70 1.11 Table 2. Association results of two three-way interactions in the discovery phase, the validation phase, and the combined data. also divided patients into four subgroups based on the clinical score, a weighted linear combination of demo- 95% CI graphic and clinical variables, as shown in gray lines in 0.986 1.002 0.987 Validation phase 0.01 0.29 0.80 0.84 1.01 1.00 0.02 0.41 0.85 1.03 1.01 Fig. 3A. iTWINS obviously outperformed the clinical score in discriminability, and the association between iTWINS and overall survival remained significant in all 0.992 1.011 0.993 0.23 0.59 0.88 1.25 1.07 1.01 0.20 0.63 0.91 1.32 1.06 HR subgroup analyses stratified by various covariates (Fig. 4). Exploring the molecular characteristics of these iTWINS subgroups, we performed differential expres- FDR-q 0.039 0.046 sion analysis and identified a total of 3150 DEGs, which were significantly enriched in 106 KEGG path- ways, including NSCLC (Fig. 5A). In addition, GO 4.21E-03 4.24E-03 2.06E-06 3.91E-03 2.32E-06 4.23E-06 7.79E-06 4.97E-05 2.57E-06 1.37E-05 1.94E-06 3.15E-06 2.46E-05 1.16E-05 enrichment analysis identified 284 significant biological process pathways (Fig. 5B), 76 significant cellular com- P ponent pathways (Fig. 5C) and 90 significant molecu- lar function pathways (Fig. 5D). Moreover, iTWINS 0.996 1.008 0.998 0.53 0.88 0.94 1.48 1.10 1.02 0.46 0.83 0.97 1.27 1.04 was significantly but negatively correlated with the immune score (r = −0.096, P = 0.018) (Fig. S5A). For 95% CI example, compared to Group 4 of iTWINS that had Discovery phase 0.990 1.003 0.996 0.03 0.51 0.86 1.08 1.04 1.01 0.10 0.64 0.93 1.10 1.02 the worst survival, Group 1 with the best survival had the highest immune score (P = 0.0092) (Fig. S5B), indicating that abundant immune cells are beneficial 0.993 1.005 0.997 0.13 0.67 0.90 1.26 1.07 1.01 0.22 0.73 0.95 1.18 1.03 for NSCLC survival. Further, the compositions of 14 HR types of immune cells were differently distributed among iTWINS subgroups (Fig. 6A); the correlation cg00060500 × pack-year cg17479956 × pack-year of iTWINS with each immune cell composition varied, including negative (e.g., with T-cell subsets gamma delta, r = −0.29, P = 2.67 × 10-13) and positive corre- lations (e.g., with M2 macrophages, r = 0.32, cg00060500 cg17479956 pack-year pack-year pack-year pack-year P = 3.17 × 10-16) (Fig. 6B). 3.3. Three-way interaction empowered × × × × × × × × cg05293407 cg00060500 cg05293407 cg05293407 cg00060500 cg05293407 cg05293407 cg17479956 cg05293407 cg05293407 cg17479956 cg05293407 prognostic prediction model Pack-year pack-year We built a prognostic prediction model of NSCLC by incorporating significant three-way interactions. The Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of 723 Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Three-way interaction study of NSCLC X. Ji et al. Fig. 2. The three-way interaction pattern between pack-year of smoking, cg05293407TRIM27 and cg00060500KIAA0226. (A) The pattern illustrated by a 3D plot. (B-E) Kaplan–Meier survival analysis of cg00060500KIAA0226 stratified by pack-year of smoking and cg05293407TRIM27. The number of patients in each subgroup was 266, 282, 211 and 233. Hazard ratio (HR), 95% CI, and P-values were derived from a Cox proportional hazards regression model adjusted for age, sex, clinical stage and study centers. Fig. 3. Estimated survival curves for patients among iTWINS subgroups. (A) Kaplan–Meier survival curves for patients grouped by iTWINS subgroups. Patients were categorized into four subgroups by using the quantile of iTWINS as the cutoffs. The number of patients in groups 1-4 was 248. (B) Hazard ratio (HR) and P-values were derived from the Cox proportional hazards model for patients with different quartile levels of the iTWINS. model with only demographic and clinical variables had pack-year of smoking (AUC3-year = 0.67, AUC5- a limited prediction ability (AUC3-year = 0.64, AUC5- year = 0.71, C-index = 0.63) or the two significant three- year = 0.68, C-index = 0.62), and the prediction accuracy way interactions (AUC3-year = 0.72, AUC5-year = 0.74, C- slightly increased by adding the cg05293407TRIM27 × index = 0.66). However, by adding 42 three-way 724 Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License X. Ji et al. Three-way interaction study of NSCLC Fig. 4. Forest plots of results from stratification analysis of iTWINS. Hazard ratio (HR) with 95% CI of iTWINS on non-small-cell lung cancer (NSCLC) survival in various subgroups is stratified by clinical characteristics. LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma. Hazard ratio (HR), 95% CI, and P-values were derived from a Cox proportional hazards regression model for patients in different subgroups. interactions obtained by using a looser criterion (discov- by interacting with cg05293407TRIM27 × pack-year of ery: FDR-q ≤ 0.10, validation: P ≤ 0.05), the AUC smoking, displayed a significant three-way interaction increased by 32.8% (95% CI: 32.6%-33.1%, effect on the overall survival of NSCLC. The proposed P = 2.20 × 10-16) and 29.4% (95% CI: 29.2%-29.6%, iTWINS had some interesting properties. It enabled us to P = 2.20 × 10-16) for 3- and 5-year survival (AUC3- categorize patients with various risk levels of mortality year = 0.85, AUC5-year = 0.88, C-index = 0.76), respec- into four subgroups, which involved a series of significant tively (Fig. 7). DEGs belonging to several biological pathways. Also, iTWINS was negatively correlated with the immune score, representing the level of infiltrating immune cells. Finally, 4. Discussion prognostic models with three-way interactions substan- We were among the first few to perform an epigenome- tially improved the prediction accuracy. wide three-way interaction study on overall survival of It is known that gene–gene and gene–environment early-stage NSCLC. We identified two CpG probes interactions on an epigenetic level play a crucial role in (cg00060500KIAA0226 and cg17479956EXT2), each of which, revealing biological mechanisms of complex diseases Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of 725 Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Three-way interaction study of NSCLC X. Ji et al. Fig. 5. The functional enrichment analyses of DEGs. (A) Top 36 significant KEGG pathways. (B) Top 36 significant biological process pathways. (C) Top 36 significant cellular component pathways and (D) top 36 significant molecular function pathways. [23], including NSCLC [11,24]. Emerging evidence from Nevertheless, EXT2 is found to be the likely antigens in higher-order interaction studies with high-throughput the immune complexes and is possibly the target antigens data [44,45] indicates that accommodation of such com- in autoimmune membranous nephropathy [47]. Mean- plex association patterns can improve the predictive while, KIAA0226 is identified as a multifaceted immune power [46]. Building upon our previous study that iden- modulatory protein that regulates autophagy and phago- tified two-way interactions contributing to the improve- cytosis [48]. They were all related to immunity, and their ment of the prediction accuracy [14], the current study gene expression was significantly correlated. Thus, the found that inclusion of three-way interactions might lar- subsequent analysis only focused on the interaction gely enhance prediction performance for both 3- and 5- involving cg00060500KIAA0226. year survival. Therefore, complex association patterns Biologically, KIAA0226 was identified as a novel (e.g., higher-order interactions) among multiple factors Beclin-1-binding protein, and hence dubbed as Rubicon should be factored in for studies of complex diseases, (RUN domain and cysteine-rich domain containing, such as NSCLC. Beclin-1-interacting protein) [49,50]. Rubicon is a multi- To further elaborate on the three-way interactions faceted immune modulatory protein that regulates among pack-year of smoking, cg05293407TRIM27 and autophagy and phagocytosis [48]. Cigarette smoke cg00060500KIAA0226/cg17479956EXT2, we generated 3D induced extrinsic apoptosis partly dependent on autop- figures and observed varying effects of hagy protein Beclin-1 and caused a reduction of cg00060500KIAA0226/cg17479956EXT2 at different levels of Beclin-1 expression [51,52]. Furthermore, decreased pack-year of smoking and cg05293407TRIM27. expression of Beclin-1 may lead to changes in Rubicon 726 Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License X. Ji et al. Three-way interaction study of NSCLC Fig. 6. The association analysis between immune cells and iTWINS. (A) Comparisons of the abundances of 22 immune cells in iTWINS subgroups. ‘ns’ means P > 0.05, ∗ means P < 0.05, ∗∗∗ means P < 0.001 and ∗∗∗∗ means P < 0.0001. (B) Correlation analysis between immune cells and iTWINS. Yellow bars meant correlation coefficient > 0 and P ≤ 0.05, blue bars meant correlation coefficient< 0 and P ≤ 0.05, and gray bars meant P > 0.05. Correlation coefficients and hypothesis tests were based on Pearson correlation tests. Fig. 7. ROC curves for different predictive models using the covariates (clinical information), two-way interaction (cg05293407TRIM27 and pack-year of smoking) and three-way interactions (cg05293407TRIM27, pack-year of smoking and another CpG probe) with FDR-q ≤ 0.05 or FDR-q ≤ 0.10. (A) Three-year survival prediction. (B) Five-year survival prediction. The AUC increase (%) was evaluated by comparing the model with that with only the covariates. P-values and 95% CIs were calculated by using 1000 bootstrap samples and t-tests. AUC, area under the receiver operating characteristic curve; ROC, receiver operating characteristic. expression, thus affecting cross-regulation between Notably, besides the heterogeneous survival among apoptosis and autophagy [53]. That said, biological four epigenetic iTWINS subgroups, there were a series mechanisms underlying this three-way interaction need of significant transcriptional DEGs enriched in cancer to be more systematically studied. related pathways, including NSCLC pathway and Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of 727 Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Three-way interaction study of NSCLC X. Ji et al. EGFR tyrosine kinase inhibitor resistance pathway of our population was Caucasian (89%), generaliza- which directly affected the overall survival of NSCLC tion of the results to the other ethnicity groups should [54]. In addition, the subgroup with the lowest be cautioned. iTWINS had the best survival and the highest immune score, consistent with the previous study [55]. Further- 5. Conclusion more, among 13 immune cell types identified in NSCLC tumors [56], 7 of them were significantly cor- Our study identified two CpG probes (cg00060500KIAA0226 related with iTWINS in our study, including T cells and cg17479956EXT2), together with pack-year of smoking CD4, T cells CD8, B cells, macrophages, eosinophils, and cg05293407TRIM27, had genome-wide significant three- mast cells, and dendritic cells. For example, M2 way interaction effects on NSCLC survival. iTWINS, macrophages were positively correlated with iTWINS, which involved in TIME, can distinguish high risk NSCLC indicating the higher abundance of M2 macrophages, patients susceptible to mortality, by significantly improving the worse the survival of NSCLC, as confirmed by sev- prognostic prediction accuracy for NSCLC survival. Our eral studies [57,58]. Our iTWINS subgroup analysis findings suggested that lung cancer progression may be dri- further indicated that potential NSCLC-related path- ven by complex interacting effects that involve multiple ways and TIICs might be the driven factors of lung environmental and epigenetic factors, which can provide cancer progression. potential therapeutic targets of NSCLC. Our study has several strengths. First, to our knowl- edge, this is the first three-way interaction study Acknowledgements between DNA methylation of CpG probes and pack- year of smoking, shedding some light on the tumor The authors thank TCGA for contributing clinical, DNA recurrence and progression. Second, we controlled methylation, and RNA sequencing data, as well as all false positives in the discovery phase and also per- study subjects who participated in the five study cohorts. formed external validation using an independent popu- This study was funded by the National Key Research lation. In our study, the final significant signals were and Development Program of China (2016YFE0204900 defined as these with FDR-q ≤ 0.05 in the discovery to FC), Natural Science Foundation of Jiangsu Province phase and meanwhile with P ≤ 0.05 in the validation (BK20191354 to RZ), National Natural Science Founda- phase. Thus, the overall false positive rate is less than tion of China (81530088 to FC, 81973142 to YYW), Nat- 0.05, which is approximatively around 0.0025 ural Science Foundation of the Jiangsu Higher Education (0.05 × 0.05). Such two-phase design is a very popular Institutions of China (18KJB310011 to RZ), China Post- and quite conservative manner used in the omics anal- doctoral Science Foundation (2018M633767 to RZ), US ysis to control false positives. Third, 3D plots were National Institutes of Health (CA209414, CA092824, designed and used to illustrate the complex three-way and U01CA209414 to DCC and YL; AG056764 to YL), interactions, facilitating interpretation and dissemina- and Priority Academic Program Development of Jiangsu tion of our results. Fourth, by exploring the molecular Higher Education Institutions (PAPD). RZ was partially characteristics of different subgroups defined by supported by the Qing Lan Project of the Higher Educa- iTWINS, we gained a better understanding of the tion Institutions of Jiangsu Province and the Outstanding heterogeneous response to immunotherapy, and there- Young Teachers Training Program of Nanjing Medical fore, a potential cause of survival heterogeneity. University. Finally, our prognostic model by incorporating three- way interactions may aid physicians in making clinical Conflict of interest decisions or guiding immunotherapy. We acknowledge limitations. First, the study did not The authors declare no conflict of interest. elucidate the biological mechanism underlying the identified three-way interactions, although the results Peer review of immune-related and enrichment analyses might pro- vide insight into the functional mechanisms. Second, The peer review history for this article is available at the censoring rate (76.07%) of survival time for the https://publons.com/publon/10.1002/1878-0261.13167. TCGA cohort is relatively high, since early-stage NSCLC patients need longer follow-up time. Thus, the Data accessibility validation phase using TCGA population had low sta- tistical power. Nonetheless, the three-way interactions The DNA methylation image data of Harvard, Spain, remained significant in TCGA. Finally, as the majority Norway, and Sweden study cohort can be requested 728 Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License X. Ji et al. Three-way interaction study of NSCLC from DCC, ME, ÅH, and JS, respectively. Alterna- 7 Egger G, Liang G, Aparicio A, Jones PA. Epigenetics tively, it can be retrieved from gene expression omni- in human disease and prospects for epigenetic therapy. bus database (GSE39279, GSE66836 and GSE56044). Nature. 2004;429:457–63. TCGA: https://tcga-data.nci.nih.gov; now hosted at 8 Feinberg AP, Tycko B. The history of cancer GDC: https://portal.gdc.cancer.gov. epigenetics. Nat Rev Cancer. 2004;4:143–53. 9 Heyn H, Esteller M. DNA methylation profiling in the clinic: applications and challenges. Nat Rev Genet. Author contributions 2012;13:679–92. 10 Zhang R, Lai L, He J, Chen C, You D, Duan W, et al. XJ, LL, JF, RZ, FC, and DCC contributed to the EGLN2 DNA methylation and expression interact with study design. RZ, LS, AS, MMB, AK, MP, JS, ÅH, HIF1A to affect survival of early-stage NSCLC. ME, and DCC contributed to data collection. XJ, LL, Epigenetics. 2019;14:118–29. JF, YL, YW, SS, RZ, FC, and DCC performed statis- 11 Zhang R, Lai L, Dong X, He J, You D, Chen C, et al. tical analysis and interpretation and drafted the manu- SIPA1L3 methylation modifies the benefit of smoking script. XJ, LL, JF, YL, YW, and SS revised the cessation on lung adenocarcinoma survival: an manuscript. All authors contributed to critical revision epigenomic-smoking interaction analysis. Mol Oncol. of the manuscript and approved its final version. 2019;13:1235–48. Financial support and study supervision were provided 12 Dong X, Zhang R, He J, Lai L, Alolga RN, Shen S, by RZ, FC, and DCC. et al. Trans-omics biomarker model improves prognostic prediction accuracy for early-stage lung adenocarcinoma. Aging. 2019;11:6312–35. Consent for publication 13 Wei Y, Liang J, Zhang R, Guo Y, Shen S, Su L, et al. All participants or their surrogate care providers gave Epigenetic modifications in KDM lysine demethylases associate with survival of early-stage NSCLC. Clin written informed consent. All authors have reviewed Epigenetics. 2018;10:41. the manuscript and consented for publication. 14 Zhang R, Chen C, Dong X, Shen S, Lai L, He J, et al. Independent validation of early-stage NSCLC Web resources prognostic scores incorporating epigenetic and transcriptional biomarkers with gene-gene interactions TCGA: https://portal.gdc.cancer.gov and main effects. Chest. 2020;158:808. https://doi.org/ 10.1016/j.chest.2020.01.048. 15 Guo Y, Zhang R, Shen S, Wei Y, Salama SM, References Fleischer T, et al. DNA Methylation of LRRC3B: a 1 Sung H, Ferlay J, Siegel RL, Laversanne M, biomarker for survival of early-stage non-small cell lung Soerjomataram I, Jemal A, et al. Global Cancer cancer patients. Cancer Epidemiol Biomarkers Prev. Statistics 2020: GLOBOCAN estimates of incidence and 2018;27:1527–35. mortality worldwide for 36 cancers in 185 countries. CA 16 Shen S, Zhang R, Guo Y, Loehrer E, Wei Y, Zhu Y, Cancer J Clin. 2021;71:209–49. et al. A multi-omic study reveals BTG2 as a reliable 2 Chen Z, Fillmore CM, Hammerman PS, Kim CF, Wong prognostic marker for early-stage non-small cell lung KK. Non-small-cell lung cancers: a heterogeneous set of cancer. Mol Oncol. 2018;12:913–24. diseases. Nat Rev Cancer. 2014;14:535–46. 17 Shen S, Wei Y, Zhang R, Du M, Duan W, Yang S, 3 Molina JR, Yang P, Cassivi SD, Schild SE, Adjei AA. Non- et al. Mutant-allele fraction heterogeneity is associated small cell lung cancer: epidemiology, risk factors, treatment, with non-small cell lung cancer patient survival. Oncol and survivorship. Mayo Clin Proc. 2008;83:584–94. Lett. 2018;15:795–802. 4 Hirsch FR, Scagliotti GV, Mulshine JL, Kwon R, Curran 18 Joo JE, Dowty JG, Milne RL, Wong EM, Dugue PA, WJ Jr, Wu YL, et al. Lung cancer: current therapies and English D, et al. Heritable DNA methylation marks new targeted treatments. Lancet. 2017;389:299–311. associated with susceptibility to breast cancer. Nat 5 Zito Marino F, Bianco R, Accardo M, Ronchi A, Commun. 2018;9:867. Cozzolino I, Morgillo F, et al. Molecular heterogeneity 19 Widschwendter M, Zikan M, Wahl B, Lempiäinen H, in lung cancer: from mechanisms of origin to clinical Paprotka T, Evans I, et al. The potential of circulating implications. Int J Med Sci. 2019;16:981–9. tumor DNA methylation analysis for the early 6 Tang S, Pan Y, Wang Y, Hu L, Cao S, Chu M, et al. detection and management of ovarian cancer. Genome Genome-wide association study of survival in early- Med. 2017;9:116. stage non-small cell lung cancer. Ann Surg Oncol. 20 Flanders WD, Lally CA, Zhu BP, Henley SJ, Thun MJ. 2015;22:630–5. Lung cancer mortality in relation to age, duration of Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of 729 Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Three-way interaction study of NSCLC X. Ji et al. smoking, and daily cigarette consumption: results from 34 Karlsson A, Jonsson M, Lauss M, Brunnstrom H, Cancer Prevention Study II. Cancer Res. 2003;63:6556– Jonsson P, Borg A, et al. Genome-wide DNA 62. methylation analysis of lung carcinoma reveals one 21 Janjigian YY, McDonnell K, Kris MG, Shen R, Sima neuroendocrine and four adenocarcinoma epitypes CS, Bach PB, et al. Pack-years of cigarette smoking as associated with patient outcome. Clin Cancer Res. a prognostic factor in patients with stage IIIB/IV 2014;20:6127–40. nonsmall cell lung cancer. Cancer. 2010;116:670–5. 35 Chen YA, Lemire M, Choufani S, Butcher DT, 22 Ji X, Lin L, Shen S, Dong X, Chen C, Li Y, et al. Grafodatskaya D, Zanke BW, et al. Discovery of cross- Epigenetic-smoking interaction reveals histologically reactive probes and polymorphic CpGs in the Illumina heterogeneous effects of TRIM27 DNA methylation on Infinium HumanMethylation450 microarray. overall survival among early-stage NSCLC patients. Epigenetics. 2013;8:203–9. Mol Oncol. 2020;14:2759–74. 36 Marabita F, Almgren M, Lindholm ME, Ruhrmann 23 Trerotola M, Relli V, Simeone P, Alberti S. Epigenetic S, Fagerstrom-Billai F, Jagodic M, et al. An inheritance and the missing heritability. Hum Genomics. evaluation of analysis pipelines for DNA 2015;9:17. methylation profiling using the Illumina 24 Chen C, Wei Y, Wei L, Chen J, Chen X, Dong X, HumanMethylation450 BeadChip platform. et al. Epigenome-wide gene-age interaction analysis Epigenetics. 2013;8:333–46. reveals reversed effects of PRODH DNA methylation 37 Benjamini Y, Hochberg Y. Controlling the false on survival between young and elderly early-stage discovery rate: a practical and powerful approach NSCLC patients. Aging. 2020;12:10642–62. to multiple testing. J Royal Stat Soc. 1995;57:289–300. 25 Mucci LA, Wedren S, Tamimi RM, Trichopoulos D, 38 Ashburner M, Ball CA, Blake JA, Botstein D, Butler Adami HO. The role of gene-environment interaction in H, Cherry JM, et al. Gene ontology: tool for the the aetiology of human cancer: examples from cancers unification of biology. The Gene Ontology Consortium. of the large bowel, lung and breast. J Intern Med. Nat Genet. 2000;25:25–9. 2001;249:477–93. 39 Kanehisa M, Furumichi M, Tanabe M, Sato Y, 26 Kraft P, Hunter DJ. Genetic risk prediction–are we Morishima K. KEGG: new perspectives on genomes, there yet? N Engl J Med. 2009;360:1701–3. pathways, diseases and drugs. Nucleic Acids Res. 27 Tekin E, White C, Kang TM, Singh N, Cruz-Loya M, 2017;45:D353–61. Damoiseaux R, et al. Prevalence and patterns of higher- 40 Wang J, Vasaikar S, Shi Z, Greer M, Zhang B. order drug interactions in Escherichia coli. NPJ Syst WebGestalt 2017: a more comprehensive, powerful, Biol Appl. 2018;4:31. flexible and interactive gene set enrichment analysis 28 Basu S, Kumbier K, Brown JB, Yu B. Iterative random toolkit. Nucleic Acids Res. 2017;45:W130–7. forests to discover predictive and stable high-order 41 Yoshihara K, Shahmoradgoli M, Martı́n ez E, interactions. Proc Natl Acad Sci USA. 2018;115:1943–8. Vegesna R, Kim H, Torres-Garcia W, et al. Inferring 29 Woods NT, Monteiro AN, Thompson ZJ, Amankwah tumour purity and stromal and immune cell EK, Naas N, Haura EB, et al. Interleukin admixture from expression data. Nat Commun. polymorphisms associated with overall survival, disease- 2013;4:2612. free survival, and recurrence in non-small cell lung cancer 42 Newman AM, Liu CL, Green MR, Gentles AJ, Feng patients. Mol Carcinog. 2015;54(Suppl 1):E172–84. W, Xu Y, et al. Robust enumeration of cell subsets 30 Taylor MB, Ehrenreich IM. Higher-order genetic from tissue expression profiles. Nat Methods. interactions and their contribution to complex traits. 2015;12:453–7. Trends Genet. 2015;31:34–40. 43 Heagerty PJ, Lumley T, Pepe MS. Time-dependent 31 Asomaning K, Miller DP, Liu G, Wain JC, Lynch TJ, ROC curves for censored survival data and a diagnostic Su L, et al. Second hand smoke, age of exposure and marker. Biometrics. 2000;56:337–44. lung cancer risk. Lung Cancer. 2008;61:13–20. 44 Anastassiou D. Computational analysis of the synergy 32 Sandoval J, Mendez-Gonzalez J, Nadal E, Chen G, among multiple interacting genes. Mol Syst Biol. Carmona FJ, Sayols S, et al. A prognostic DNA 2007;3:83. methylation signature for stage I non-small-cell lung 45 Watkinson J, Liang KC, Wang X, Zheng T, Anastassiou cancer. J Clin Oncol. 2013;31:4140–7. D. Inference of regulatory gene interactions from 33 Bjaanaes MM, Fleischer T, Halvorsen AR, Daunay A, expression data using three-way mutual information. Busato F, Solberg S, et al. Genome-wide DNA Ann N Y Acad Sci. 2009;1158:302–13. methylation analyses in lung adenocarcinomas: 46 Ma S, Kosorok MR, Huang J, Dai Y. Incorporating Association with EGFR, KRAS and TP53 mutation higher-order representative features improves prediction status, gene expression and prognosis. Mol Oncol. in network-based cancer prognosis analysis. BMC Med 2016;10:330–43. Genomics. 2011;4:5. 730 Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of Federation of European Biochemical Societies. 18780261, 2022, 3, Downloaded from https://febs.onlinelibrary.wiley.com/doi/10.1002/1878-0261.13167 by Readcube (Labtiva Inc.), Wiley Online Library on [20/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License X. Ji et al. Three-way interaction study of NSCLC 47 Sethi S, Madden BJ, Debiec H, Charlesworth MC, prognostic signature in lung adenocarcinoma. Cancer Gross L, Ravindran A, et al. Exostosin 1/exostosin 2- Med. 2020;9:5960–75. associated membranous nephropathy. J Am Soc 56 Stankovic B, Bjørhovde HAK, Skarshaug R, Aamodt Nephrol. 2019;30:1123–36. H, Frafjord A, Müller E, et al. Immune cell 48 Yang CS, Lee JS, Rodgers M, Min CK, Lee JY, Kim composition in human non-small cell lung cancer. Front HJ, et al. Autophagy protein Rubicon mediates Immunol. 2018;9:3101. phagocytic NADPH oxidase activation in response to 57 Hwang I, Kim JW, Ylaya K, Chung EJ, Kitano H, microbial infection or TLR stimulation. Cell Host Perry C, et al. Tumor-associated macrophage, Microbe. 2012;11:264–76. angiogenesis and lymphangiogenesis markers predict 49 Zhong Y, Wang QJ, Li X, Yan Y, Backer JM, Chait BT, prognosis of non-small cell lung cancer patients. J et al. Distinct regulation of autophagic activity by Atg14L Transl Med. 2020;18:443. and Rubicon associated with Beclin 1-phosphatidylinositol- 58 Guo Z, Song J, Hao J, Zhao H, Du X, Li E, et al. M2 3-kinase complex. Nat Cell Biol. 2009;11:468–76. macrophages promote NSCLC metastasis by 50 Matsunaga K, Saitoh T, Tabata K, Omori H, Satoh upregulating CRYAB. Cell Death Dis. 2019;10:377. T, Kurotori N, et al. Two Beclin 1-binding proteins, Atg14L and Rubicon, reciprocally regulate autophagy at different stages. Nat Cell Biol. Supporting information 2009;11:385–96. Additional supporting information may be found 51 Kim HP, Wang X, Chen ZH, Lee SJ, Huang MH, online in the Supporting Information section at the end Wang Y, et al. Autophagic proteins regulate cigarette of the article. smoke-induced apoptosis: protective role of heme Table S1. Demographic and clinical descriptions of oxygenase-1. Autophagy. 2008;4:887–95. early-stage NSCLC patients with gene expression data 52 Pabón MA, Patino E, Bhatia D, Rojas-Quintero J, Ma in four international study centers. KC, Finkelsztein EJ, et al. Beclin-1 regulates cigarette Table S2. The association results of variables derived smoke-induced kidney injury in a murine model of from Cox proportional hazards model adjusted for chronic obstructive pulmonary disease. JCI Insight. covariates in NSCLC samples. 2018;3:18. Fig. S1. Quality control processes for DNA methyla- 53 Kang R, Zeh HJ, Lotze MT, Tang D. The Beclin 1 network regulates autophagy and apoptosis. Cell Death tion data. Differ. 2011;18:571–80. Fig. S2. Correlation between DNA methylation probe 54 Zhu L, Chen Z, Zang H, Fan S, Gu J, Zhang G, et al. and its corresponding gene expression. Targeting c-Myc to overcome acquired resistance of Fig. S3. The three-way interaction pattern between EGFR mutant NSCLC cells to the third generation pack-year of smoking, cg05293407TRIM27 and EGFR tyrosine kinase inhibitor, osimertinib. Cancer cg17479956EXT2. Res. 2021;81:4822–34. Fig. S4. Association between KIAA0226 and EXT2. 55 Sun S, Guo W, Wang Z, Wang X, Zhang G, Zhang H, Fig. S5. The association analysis between immune et al. Development and validation of an immune-related score and iTWINS. Molecular Oncology 16 (2022) 717–731 ª 2021 The Authors. Molecular Oncology published by John Wiley & Sons Ltd on behalf of 731 Federation of European Biochemical Societies.