OL Oncology Letters 1792-1074 1792-1082 D.A. Spandidos 10.3892/ol.2019.10366 OL-0-0-10366 Articles Predicting the survival of patients with lung adenocarcinoma using a four-gene prognosis risk model ZhangWei ShenYang FengGanzhu Department of Respiratory Medicine, The Second Affiliated Hospital of Nanjing Medical University, Nanjing, Jiangsu 210011, P.R. China Correspondence to: Dr Ganzhu Feng, Department of Respiratory Medicine, The Second Affiliated Hospital of Nanjing Medical University, 121 Jiangjiayuan Road, Nanjing, Jiangsu 210011, P.R. China, E-mail: zhangwei3125@sina.com 07 2019 17 05 2019 18 1 535 544 06052018 31012019 Copyright: © Zhang et al. 2019 This is an open access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

Lung adenocarcinoma (LAD) is difficult to diagnose as it tends to be small in size and metastasize early. The aim of the present study was to investigate prognostic factors for patients with LAD and establish a prognosis risk model. A training set consisting of clinical and RNA sequencing data from 503 patients with LAD, as well as expression data from a further 59 LAD and adjacent tissues, was obtained from The Cancer Genome Atlas. Additionally, a validation dataset was acquired from the Gene Expression Omnibus database (GSE26939), which included clinical and gene expression data from 115 patients. Using the DESeq2 package to compare expression between LAD and adjacent tissues, differentially expressed genes (DEGs) were identified. On the basis of survival and the random forests for survival, regression and classification package, genes for constructing the prognosis risk model were selected. The prognosis risk model was constructed and validated using the survival package. Subsequently, high- and low-risk groups were compared using the Limma package to identify DEGs, and enrichment analysis was performed using the web-based gene set analysis toolkit. A protein-protein interaction network was visualized using Cytoscape software. There were 18,567 DEGs between the LAD samples and the adjacent tissues, and 363 DEGs between the high- and low-risk groups. Of these, four genes were selected for constructing the prognosis risk model, myosin IE (MYO1E), endoplasmic reticulum oxidoreductase 1α (ERO1L), C1q and tumor necrosis factor-related protein 6 (C1QTNF6) and family with sequence similarity 83, member A (FAM83A). The survival time of high- and low-risk groups in the validation set were significantly different. Functional enrichment revealed that the genes that interacted with MYO1E, ERO1L, C1QTNF6 and FAM83A separately were enriched in ‘cell cycle regulation’, ‘synthesis and assembly of nucleic acids’, ‘histone modification and cell cycle progression’ and ‘cell secretion process’. The four-gene prognosis risk model could potentially be used for predicting the survival of patients with LAD.

 lung adenocarcinoma differentially expressed genes survival analysis prognosis risk model protein-protein interaction network Introduction

Lung adenocarcinoma (LAD), also termed pulmonary adenocarcinoma, is a common type of lung cancer that often occurs in the peripheral lung tissue (1,2). LAD is induced by smoking, but is also the most frequent type of lung cancer in non-smokers (3,4). The occurrence of LAD in the periphery of the lung may be due to the nature of cigarette filters, which keep larger particles out of the lungs (5). LAD is difficult to diagnose owing to its small mass, and it readily enters the regional lymph nodes and tends to metastasize at an early stage (6). LAD is a type of non-small cell lung carcinoma, and its response to radiation therapy is worse compared with that of small cell lung carcinoma (7). Further understanding of the mechanism of LAD development may result in improved prognosis for patients.

Sex-determining region box 2 can serve as a marker predicting poor outcome and may assist in identifying a high risk of recurrence following surgery in patients with stage I LAD (8,9). Zhang et al (10) reported an analysis of the clinicopathological features of LAD that identified the upregulation of programmed death 1 ligand 1 (PD-L1) and PD-L2, indicating that PD-L expression is an independent predictor of LAD (10). The high expression of metastasis-associated in colon cancer 1 may be a promising predictor of postoperative recurrence in patients with LAD following surgery (11,12). Upregulation of the cytochrome P450 family 24 subfamily A member 1 (CYP24A1) in patients with LAD is associated with poor survival, which may be due to the abrogation of anti-proliferative effects of the active form of vitamin D, 1α,25-dihydroxyvitamin D3 (13). Upregulation of histone deacetylase 1 in tumor cells can be used as a marker for poor prognosis in patients with LAD (14). Nevertheless, the previous studies have not comprehensively revealed prognostic factors for LAD. In the present study, RNA-sequencing data of LAD tissues was obtained from The Cancer Genome Atlas (TCGA), and the genes associated with the survival rate of patients with LAD were explored. On this basis, a prognosis risk model was constructed, and validated using microarray data from a separate cohort of patients with LAD. Furthermore, the biological significance of the prognostic differences in patients with LAD was investigated.

 Materials and methods <sec> <title>Expression profile data

RNA-sequencing data, including raw count and fragments per kilobase of transcript per million fragments mapped (FPKM) and the relevant clinical data were downloaded from the Genomic Data Commons database (GDC; http://portal.gdc.cancer.gov). There were expression data and clinical data for 503 tissues from patients with LAD in the TCGA dataset (TCGA-LAD), which were taken as the training set for constructing the prognosis risk model for LAD. Additionally, the expression data of a further 59 pairs of LAD tissues and adjacent tissues were obtained.

Gene expression data and clinical data were acquired from the Gene Expression Omnibus database (GEO; http://www.ncbi.nlm.nih.gov/geo) (GSE26939). A total of 115 samples were selected as the validation set for the prognosis risk model. The clinical data of the training set and the validation set are presented in Table I.

 Screening for differentially expressed genes (DEGs) between LAD and adjacent tissues

Using raw data from 59 pairs of LAD and adjacent tissues, DEGs between the two groups were identified using the R 3.1.0 software package DESeq2 1.4.5 (http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html) (15). The thresholds were log2 fold change (|log2 FC|) >1 and false discovery rate (FDR) <0.05.

 Identification of genes for constructing the prognosis risk model

In combination with the clinical data and the gene expression data (FPKM values) for each patient, the Cox univariate regression analysis in the survival analysis package for R (https://cran.r-project.org/web/packages/survival/index.html) (16) was used to identify the survival time-associated genes from the DEGs. The log-rank test (17), sc-test (Rao's score test) (18) and Wald test (19) were utilized to screen genes with P<0.01. Subsequently, the random forests for survival, regression and classification package in R (https://cran.r-project.org/web/packages/randomForestSRC/index.html) (20) was applied for ranking the survival time-associated DEGs by survival contribution degrees. Finally, the genes with survival contribution degree >0.002 were selected as the genes for constructing the prognosis risk model.

 Construction and evaluation of prognosis risk model

According to the FPKM values of the selected genes and the survival status of the patients, the risk scoring formula for prognosis of LAD was constructed using the Cox multivariate regression analysis in R package survival (16). Subsequently, the risk scores of the 503 patients with LAD were calculated. To assess the predictive independence of the risk scores in the training set, the risk scores and eight clinical factors, including age, sex, smoking status, stage, epidermal growth factor receptor (EGFR), KRAS proto-oncogene, serine/threonine kinase 11 and subtype were subjected to Cox univariate and multivariate regression analyses. In addition, stratified analysis was performed for other survival-associated factors. Thereafter, the patients were divided into high- and low-risk groups according to the median of risk scores. The Kaplan-Meier (KM) survival curve (21) was used to evaluate the survival difference between the two groups of patients. Significant differences in the tumor stages of patients in the high- and low-risk groups were identified using a χ2 test with a threshold value of P<0.05.

 Validation of the prognosis risk model

Using the expression and clinical data from the GSE26939 dataset, the risk scores of 115 patients with LAD were calculated using the risk scoring formula. Similarly, the 115 patients were classified into high- and low-risk groups. The KM survival curve (21) was used to assess the survival difference between the two groups of patients.

 Biological analysis of DEGs between the high- and low-risk groups

On the basis of the prognosis risk model, the 503 LAD samples were divided into two groups. Using the R package Limma (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (22), DEGs between the two groups were identified, using the thresholds of |log2 FC| >1 and FDR <0.05. Next, the functional and pathway enrichment analyses were conducted for the DEGs using the web-based gene set analysis toolkit (WebGestalt) (http://www.webgestalt.org/option.php) (23).

 Protein-protein interaction (PPI) network analysis

Using the PPI information in the InWeb_InBioMap database (https://www.intomics.com/inbio/map/#home) (24), the PPI network for the genes involved in the prognosis risk model was constructed using Cytoscape software (http://www.cytoscape.org) (25). The WebGestalt tool (23) was used to perform functional enrichment analysis for the genes identified to interact with the genes of the prognosis risk model.

 Results <sec> <title>Screening for DEGs between LAD and adjacent tissues

There were a total of 18,567 DEGs between the LAD and adjacent tissues, including 15,483 genes that were significantly upregulated in the LAD tissues and 3,084 genes that were significantly downregulated.

 Identification of genes for constructing the prognosis risk model

Cox univariate regression analysis showed that a total of 279 DEGs were correlated with survival time of patients. Subsequently, the 279 DEGs were ranked by survival contribution degrees, and the top four genes were selected for constructing the prognosis risk model (Table II). These four genes were family with sequence similarity 83, member A (FAM83A; contribution degree=0.00309), myosin IE (MYO1E; contribution degree=0.00231), C1q and tumor necrosis factor-related protein 6 (C1QTNF6; contribution degree=0.00226) and endoplasmic reticulum oxidoreductase 1α (ERO1L; contribution degree=0.00211).

 Construction and evaluation of prognosis risk model

According to the FPKM values of the four selected genes and the survival status of the patients, the following risk scoring formula for prognosis of LAD was constructed: Risk score = (0.00635 × FAM83A) + (0.02114 × MYO1E) + (0.04234 × C1QTNF6) + (0.00634 × ERO1L).

The risk scores of the 503 patients with LAD were calculated, and then the risk scores and the eight clinical factors were subjected to correlation and regression analyses. The results indicated that the risk score and the stage of LAD were significantly correlated with patient survival time (Table III). Based on the median risk score, the patients were divided into high- and low-risk groups. The risk scores and survival status of the 2 groups were statistically analyzed, showing that the high-risk group had more mortality cases (high-risk group, 114/252 patients; low-risk group, 68/251 patients) and a significantly lower overall survival time (high-risk group, 25.9±24.8 months; low-risk group, 34.7+33.6 months; P=4.44×10−4) compared with the low-risk group (Fig. 1). Using KM survival analysis, the survival difference between the 2 groups of patients was evaluated (Fig. 2A).

Tumor stage can also result in a difference in the patient survival time (Fig. 2B); therefore, the number of patients at each stage was evaluated. There were 109, 74, 53 and 15 patients in stage I, II, III and IV in the high-risk group, respectively. In the low-risk group, there were 161, 45, 28 and 10 patients at stage I, II, III and IV, respectively. The tumor stages of patients in the high- and low-risk groups were identified to be significantly different using a χ2 test (P=1.10×10−5).

In addition, stratified analysis was performed for the tumor stage (stage I, II, III and IV). The risk scores of the patients at each stage were calculated, and the patients were divided into high- and low-risk groups. The survival difference between the two groups at each stage was assessed using KM survival analysis (Fig. 2C-F). The results suggested that the prognosis risk model could significantly differentiate between patients at early stages (P<0.05), but could not distinguish patients at later stages (P>0.05). Thus, the prognosis risk model was more efficient for patients in the early stages of LAD.

 Validation of the prognosis risk model

Using the prognosis risk model, the risk scores of the 115 patients with LAD from GSE26939 were calculated. Similarly, regression analyses were performed for the risk scores and the eight clinical factors, demonstrating that the risk score was significantly correlated with the patient survival time (Table IV).

On the basis of the risk scores, the patients were divided into high- and low-risk groups. Statistical analysis revealed a significantly lower overall survival time in the high-risk group compared with that in the low-risk group (high-risk group, 34.2±31.3 months; low-risk group, 46.9±37.9 months; P=0.026; Fig. 3). There was a significant difference in the survival time between the high- and low-risk groups in the validation set, indicating that the prognosis risk model was efficient and reliable (Fig. 4).

 Biological analysis of DEGs between the high- and low-risk groups

On the basis of the prognosis risk model, the 503 LAD samples were divided into high- and low-risk groups. A total of 363 DEGs between high- and low-risk groups were identified, including 188 genes that were upregulated in the high-risk group and 175 genes that were downregulated.

Gene Ontology enrichment analysis demonstrated that the DEGs were involved in extracellular structure (extracellular matrix organization, extracellular structure organization, collagen fibril organization), the circulatory system (circulatory system development, blood vessel development, vasculature development, angiogenesis), cell adhesion and cell migration (Fig. 5A). Pathway enrichment analysis indicated that the DEGs were mainly enriched in focal adhesion, ECM-receptor interaction and protein digestion and absorption pathways (Fig. 5B).

No significant difference was identified by χ2 test between the high- and low-risk groups in terms of the occurrence of relapse (high-risk group, 82/252 patients; low-risk group, 70/251 patients; P=0.2989). There was also no significant difference between the high- and low-risk groups in terms of the occurrence of metastasis (high-risk group, 17/188 patients with relevant information; low-risk group, 7/170 patients with relevant information; P=0.09914). This may be due to limitations in the clinical factors.

 PPI network analysis

In the PPI network for the four genes involved in the prognosis risk model, there were 266, 154, 10 and 8 nodes separately interacting with MYO1E, ERO1L, C1QTNF6 and FAM83A (Fig. 6). Functional enrichment analysis indicated that the interacting genes of MYO1E, ERO1L, C1QTNF6 and FAM83A were mainly involved in M/G1 transition of mitotic cell cycle, ribosome biogenesis, regulation of centrosome duplication and mitotic cell cycle, respectively (Fig. 7).

 Discussion

In the present study, a total of 18,567 DEGs (including 15,483 upregulated and 3,084 downregulated genes) between LAD and adjacent tissues were screened. According to the survival contribution degrees, the top four genes (MYO1E, ERO1L, C1QTNF6 and FAM83A) were selected to construct the prognosis risk model. The survival time of the high- and low-risk groups in the validation set were significantly different, indicating that the prognosis risk model was efficient and reliable. The 503 LAD samples were divided into high- and low-risk groups on the basis of the prognosis risk model, and 363 DEGs (including 188 upregulated and 175 down--regulated genes in the high-risk group) between the two groups were identified. There was no significant difference between the high- and low-risk groups in terms of the occurrence of relapse or metastasis, suggesting that the prognosis risk model is limited.

The oxidizing enzyme encoded by the ERO1L gene is expressed in the endoplasmic reticulum, and its expression in a variety of cancer cells is higher compared with that in normal tissue (26). ERO1L contributes to the formation of disulfide bonds on the cell surface and secreted proteins, in combination with the protein disulfide-isomerase, and may be a potential target for cancer immunotherapy (27). ERO1L overexpression is detected in various tumor types, and it suppresses the T cell response in the tumor by promoting the production of myeloid-derived suppressor cells through the regulation of cytokines and chemokines (28,29). ERO1L serves an important role in mediating tumor growth factors, including vascular endothelial growth factor, and its expression has been correlated with poor prognosis in patients with cancer (30). Therefore, ERO1L may be involved in the prognosis of patients with LAD.

FAM83 family members (FAM83A, FAM83B and FAM83D) exhibit oncogenic potential and are involved in tumor development, and may serve as promising targets for cancer therapy (31). FAM83 proteins are transforming oncoproteins that regulate EGFR/RAS signaling, and they may prove to be targets for novel tumor-targeting therapeutic strategies (32). Overexpression of FAM83 proteins leads to higher tumor grade and decreased overall survival, and thus FAM83 proteins are oncoproteins that can be targeted for tumor treatment aimed at inhibiting mitogen-activated protein kinase signaling (33,34). The mRNA expression of FAM83B in lung squamous cell carcinoma (LSCC) is significantly higher compared with that in LAD or normal lung, indicating that FAM83B may serve as a diagnostic and prognostic biomarker for LSCC (35). This indicates that FAM83A may serve an important role in the prognosis of patients with LAD.

MYO1E regulates the migration, adhesion, endocytosis and invadosome dynamics in podocyte cells, which may be used as a prognostic biomarker in clinical applications (36). Cyclin-dependent kinases mediate the onset of mitosis, cell cycle commitment and DNA synthesis, which are involved in tumor development and are potential targets for cancer treatment (37,38). Tumor-secreted factors are precursors for tumor diffusion and migration, which serve critical roles in tumor metastasis and deterioration (39). Functional enrichment analysis revealed that the genes interacting with MYO1E, ERO1L, C1QTNF6 and FAM83A were involved in cell cycle regulation, synthesis and assembly of nucleic acids, histone modification and cell cycle progression, and cell secretion process, suggesting that MYO1E, ERO1L, C1QTNF6 and FAM83A may be prognostic biomarkers for LAD.

There are a number of limitations to the present study. The results obtained in the present study were not validated through in vitro or in vivo laboratory experiments or clinical study. Therefore, further studies should be designed and performed to confirm these results. Future studies are required to analyze the correlation between MYO1E, ERO1L, C1QTNF6 and FAM83A expression and prognosis of lung adenocarcinoma in clinical samples.

In conclusion, a total of 18,567 DEGs between LAD and adjacent tissues were screened. On the basis of the prognosis risk model, 363 DEGs between the high- and low-risk groups were identified. Additionally, the four genes (including MYO1E, ERO1L, C1QTNF6 and FAM83A) utilized in the construction of the prognosis risk model may be used for predicting the prognosis of LAD. Therefore, in future clinical practice, it may be possible to predict the prognosis of patients by calculating the risk score using data from surgical specimens. In the future, this risk prognosis model should be further validated using peripheral blood from patients with LAD.

 Acknowledgements

Not applicable.

 Funding

No funding was received.

 Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

 Authors' contributions

WZ performed data analysis and wrote the manuscript. YS contributed significantly to the analysis of the data and manuscript revision. GF conceived and designed the study. All authors read and approved the final manuscript.

 Ethics approval and consent to participate

Not applicable.

 Patient consent for publication

Not applicable.

 Competing interests

The authors declare that they have no competing interests.

 References TravisWDPathology of lung cancerClin Chest Med32669692201110.1016/j.ccm.2011.08.00522054879 McGuireSWorld cancer report 2014. Geneva, Switzerland: World health organization, international agency for research on cancer, WHO press, 2015Adv Nut7418419201610.3945/an.116.012211 LongoDLFauciASKasperDLHauserSLJamesonJLLoscalzoJHarrison's principles of internal medicine18thMcGraw Hill ProfessionalNew York, NYC, USA2872011 SubramanianJGovindanRLung cancer in never smokers: A reviewJ Clin Oncol25561570200710.1200/JCO.2006.06.801517290066 MarugameTSobueTNakayamaTSuzukiTKuniyoshiHSunagawaKGenkaKNishizawaNNatsukawaSKuwaharaOTsuburaEFilter cigarette smoking and lung cancer risk; a hospital-based case-Control study in JapanBr J Cancer90646651200410.1038/sj.bjc.660156514760379 SpiroSGTannerNTSilvestriGAJanesSMLimEVansteenkisteJFPirkerRLung cancer: Progress in diagnosis, staging and therapyRespirology154450201010.1111/j.1440-1843.2009.01674.x20199634 FukuiTMitsudomiTSmall peripheral lung adenocarcinoma: Clinicopathological features and surgical treatmentSurgery Today40191198201010.1007/s00595-008-4100-420180071 ShollLMBarlettaJAYeapBYChirieacLRHornickJLSox2 protein expression is an independent poor prognostic indicator in stage I lung adenocarcinomaAm J Surg Pathol3411931198201010.1097/PAS.0b013e3181e5e02420631605 NakatsugawaMTakahashiAHirohashiYTorigoeTInodaSMuraseMAsanumaHTamuraYMoritaRMichifuriYSOX2 is overexpressed in stem-like cells of human lung adenocarcinoma and augments the tumorigenicityLab Invest9117961804201110.1038/labinvest.2011.14021931300 ZhangYWangLLiYPanYWangRHuHLiHLuoXYeTSunYChenHProtein expression of programmed death 1 ligand 1 and ligand 2 independently predict poor prognosis in surgically resected lung adenocarcinomaOnco Targets Ther7567573201410.2147/OTT.S5995924748806 ShimokawaHUramotoHOnitsukaTChundongGHanagiriTOyamaTYasumotoKOverexpression of MACC1 mRNA in lung adenocarcinoma is associated with postoperative recurrenceJ Thorac Cardiovasc Surg141895898201110.1016/j.jtcvs.2010.09.04421093878 WangZLiZWuCWangYXiaYChenLZhuQChenYMACC1 overexpression predicts a poor prognosis for non-small cell lung cancerMed Oncol31790201410.1007/s12032-013-0790-624310811 ChenGKimSHKingANZhaoLSimpsonRUChristensenPJWangZThomasDGGiordanoTJLinLCYP24A1 is an independent prognostic marker of survival in patients with lung adenocarcinomaClin Cancer Res17817826201110.1158/1078-0432.CCR-10-178921169243 MinamiyaYOnoTSaitoHTakahashiNItoMMitsuiMMotoyamaSOgawaJExpression of histone deacetylase 1 correlates with a poor prognosis in patients with adenocarcinoma of the lungLung Cancer74300304201110.1016/j.lungcan.2011.02.01921466904 LoveMIHuberWAndersSModerated estimation of fold change and dispersion for RNA-seq data with DESeq2Genome Biol15550201410.1186/s13059-014-0550-825516281 TherneauTMGrambschPMModeling survival data: Extending the cox modelSpringerNew York200010.1007/978-1-4757-3294-8 BlandJMAltmanDGThe logrank testBMJ3281073200410.1136/bmj.328.7447.107315117797 BeraAKBiliasYRao's score, Neyman's C(α) and Silvey's LM tests: An essay on historical developments and some new resultsJ Stat Plan Inference97944200110.1016/S0378-3758(00)00343-8 HarrellFERegression modeling strategies: With applications to linear models, logistic regression, and survival analysisSpringerNew York2001 IshwaranHKogalurUBBlackstoneEHLauerMSRandom survival forestsAnn Appl Stat2841860200810.1214/08-AOAS169 StelVSDekkerFWTripepiGZoccaliCJagerKJSurvival analysis I: The Kaplan-Meier methodNephron Clin Pract119c83c88201110.1159/00032475821677442 RitchieMEPhipsonBWuDHuYLawCWShiWSmythGKlimma powers differential expression analyses for RNA-sequencing and microarray studiesNucleic Acids Res43e47201510.1093/nar/gkv00725605792 WangJDuncanDShiZZhangBWEB-based GEne SeT anaLysis Toolkit (WebGestalt): Update 2013Nucleic Acids Res41W77W83201310.1093/nar/gkt43923703215 LiTWernerssonRHansenRBHornHMercerJSlodkowiczGorkmanCTRiginaORapackiKStaerfeldtHHA scored human protein-protein interaction network to catalyze genomic interpretationNat Methods146164201710.1038/nmeth.408327892958 ShannonPMarkielAOzierOBaligaNSWangJTRamageDAminNSchwikowskiBIdekerTCytoscape: A software environment for integrated models of biomolecular interaction networksGenome Res1324982504200310.1101/gr.123930314597658 KutomiGTamuraYTanakaTKajiwaraTKukitaKOhmuraTShimaHTakamaruTSatomiFSuzukiYHuman endoplasmic reticulum oxidoreductin 1-α is a novel predictor for poor prognosis of breast cancerCancer Sci10410911096201310.1111/cas.1217723578220 TanakaTKutomiGKajiwaraTKukitaKKochinVKanasekiTTsukaharaTHirohashiYTorigoeTOkamotoYCancer-associated oxidoreductase ERO1-α promotes immune escape through up-regulation of PD-L1 in human breast cancerOncotarget82470624718201710.18632/oncotarget.1496028160557 TanakaTKajiwaraTTorigoeTOkamotoYSatoNTamuraYCancer-associated oxidoreductase ERO1-α drives the production of tumor-promoting myeloid-derived suppressor cells via oxidative protein foldingJ Immunol19420042010201510.4049/jimmunol.140253825595776 KajiwaraTTanakaTKukitaKKutomiGSaitoKOkuyaKTakayaAKochinVKanasekiTTsukaharaTHypoxia augments MHC class I antigen presentation via facilitation of ERO1-α-mediated oxidative folding in murine tumor cellsEur J Immunol4628422851201610.1002/eji.20164652527667124 TanakaTKutomiGKajiwaraTKukitaKKochinVKanasekiTTsukaharaTHirohashiYTorigoeTOkamotoYCancer-associated oxidoreductase ERO1-α drives the production of VEGF via oxidative protein folding and regulating the mRNA levelBr J Cancer11412271234201610.1038/bjc.2016.10527100727 SnijdersAMLeeSYHangBHaoWBissellMJMaoJHFAM83 family oncogenes are broadly involved in human cancers: An integrative multi-omics approachMol Oncol11167179201710.1002/1878-0261.1201628078827 BartelCAParameswaranNCiprianoRJacksonMWFAM83 proteins: Fostering new interactions to drive oncogenic signaling and therapeutic resistanceOncotarget75259752612201610.18632/oncotarget.954427221039 CiprianoRMiskimenKLBrysonBLFoyCRBartelCAJacksonMWConserved oncogenic behavior of the FAM83 family regulates MAPK signaling in human cancerMol Cancer Res1211561165201410.1158/1541-7786.MCR-13-028924736947 LinBChenTZhangQLuXZhengZDingJLiuJYangZGengLWuLFAM83D associates with high tumor recurrence after liver transplantation involving expansion of CD44+ carcinoma stem cellsOncotarget77749577507201610.18632/oncotarget.1271527769048 OkabeNEzakiJYamauraTMutoSOsugiJTamuraHImaiJItoEYanagisawaYHonmaRFAM83B is a novel biomarker for diagnosis and prognosis of lung squamous cell carcinomaInt J Oncol469991006201510.3892/ijo.2015.281725586059 ScholzTRumpAUtaPHinrichsMThielCTsiavaliarisGMyosin-1 function extends to microtubule-dependent processesBiophys J98(Suppl)561a201010.1016/j.bpj.2009.12.3040 MalumbresMBarbacidMCell cycle kinases in cancerCurr Opin Genet Dev176065200710.1016/j.gde.2006.12.00817208431 AsgharUWitkiewiczAKTurnerNCKnudsenESThe history and future of targeting cyclin-dependent kinases in cancer therapyNat Rev Drug Discov14130146201510.1038/nrd450425633797 PeinadoHLavotshkinSLydenDThe secreted factors responsible for pre-metastatic niche formation: Old sayings and new thoughtsSemin Cancer Biol21139146201110.1016/j.semcancer.2011.01.00221251983

Risk scores and survival status of the patients in the high- and low-risk groups of the training set. (A) Risk scores of the patients in the high (red) and low (blue) -risk groups. (B) Survival status of the patients in the two groups. Blue and red represent surviving and deceased patients, respectively. (C) Expression of ERO1L, FAM83A, MYO1E and C1QTNF6 in each patient. ERO1L, endoplasmic reticulum oxidoreductase 1α; FAM83A, family with sequence similarity 83, member A; MYO1E, myosin IE; C1QTNF6, C1q and tumor necrosis factor-related protein 6.

Kaplan-Meier survival curves comparing overall survival rates on the basis of tumor stage and risk score. (A) Patients in the high- and low-risk groups on the basis of median risk score of the prognosis risk model. (B) Patients with tumors at all stages. Patients with tumors at (C) stage I, (D) stage II, (E) stage III and (F) stage IV in the high- and low-risk groups according to risk score.

Risk scores and survival status of patients in the high- and low-risk groups of the validation set. (A) Risk scores of the patients in the high (red) and low (blue) -risk groups. (B) Survival status of the patients in the two groups. Blue and red represent surviving and deceased patients, respectively. (C) Expression of ERO1L, FAM83A, MYO1E and C1QTNF6 in each patient. ERO1L, endoplasmic reticulum oxidoreductase 1α; FAM83A, family with sequence similarity 83, member A; MYO1E, myosin IE; C1QTNF6, C1q and tumor necrosis factor-related protein 6.

Kaplan-Meier survival curve for the high- and low-risk groups of the validation set.

(A) GO terms and (B) pathways enriched for differentially expressed genes between the high- and low-risk groups. GO, Gene Ontology; ECM, extracellular matrix.

Protein-protein interaction network constructed for the four genes, ERO1L, FAM83A, MYO1E and C1QTNF6. Red diamonds and green circles represent the four genes and the nodes interacting with them, respectively. The lines represent interaction between the nodes. ERO1L, endoplasmic reticulum oxidoreductase 1α; FAM83A, family with sequence similarity 83, member A; MYO1E, myosin IE; C1QTNF6, C1q and tumor necrosis factor-related protein 6.

Functional GO terms enriched for the interacting genes of (A) MYO1E, (B) ERO1L, (C) FAM83A and (D) C1QTNF6. GO, Gene Ontology; MYO1E, myosin IE; ERO1L, endoplasmic reticulum oxidoreductase 1α; FAM83A, family with sequence similarity 83, member A; C1QTNF6, C1q and tumor necrosis factor-related protein 6; ncRNA, non-coding RNA; ER, endoplasmic reticulum; COPII, coatomer protein II.

Clinical characteristics of patients in the training set and validation set.

Clinical characteristic Training set (n=503) Validation set (n=115)
Age, years ± SD 63.0±16.4 63.9±10.9
Sex, n
  Female 233 50
  Male 270 49
  Unknown 0 16
Smoking status, n
  Smoker 173 100
  Non-smoker 28 12
  Unknown 302 3
Stage, n
  I 270 55
  II 119 16
  III 81 14
  IV 25 2
  Unknown 8 28
Epidermal growth factor receptor, n
  WT 189 70
  Mutated 31 11
  Unknown 283 34
KRAS proto-oncogene, n
  WT 154 88
  Mutated 66 20
  Unknown 283 7
Serine/threonine kinase 11, n
  WT 182 75
  Mutated 38 6
  Unknown 283 34
Subtype, n
  Bronchioid 85 46
  Magnoid 61 40
  Squamoid 74 29
  Unknown 283 0
Survival time, months ± SD 30.3±29.9 40.5±35.1
Survival status, n
  Alive 320 49
  Deceased 183 66

WT, wild-type; SD, standard deviation.

Univariate and multivariate Cox regression analyses for FAM83A, MYO1E, C1QTNF6 and ERO1L.

Univariate Cox regression Multivariate Cox regression
Gene Coefficient HR P-value Coefficient HR P-value
FAM83A 9.94×10−3 1.01 1.29×10−13 6.35×10−3 1.01 9.05×10−5
MYO1E 4.18×10−2 1.04 2.86×10−7 2.11×10−2 1.02 5.87×10−2
C1QTNF6 8.59×10−2 1.09 1.28×10−8 4.23×10−2 1.04 2.22×10−2
ERO1L 1.27×10−2 1.01 2.00×10−8 6.34×10−3 1.01 2.29×10−2

FAM83A, family with sequence similarity 83, member A; MYO1E, myosin IE; C1QTNF6, C1q and tumor necrosis factor-related protein 6; ERO1L, endoplasmic reticulum oxidoreductase 1α; HR, hazard ratio.

Univariate and multivariate Cox regression analyses for risk score and eight clinical factors in the training set.

Univariate Cox regression Multivariate Cox regression
Clinical factor P-value HR Lower 95% CI Upper 95% CI P-value HR Lower 95% CI Upper 95% CI
Risk score <0.001 1.66 1.48 1.86 <0.001 1.73 1.36 2.19
Age 0.05 1.01 1.00 1.02 0.02 1.01 1.00 1.03
Sex (female/male) 0.78 1.04 0.78 1.39 0.49 1.20 0.71 2.04
Smoking status (smoker/non-smoker) 0.30 0.73 0.42 1.29 0.15 0.63 0.33 1.19
Stage (I/II/III/IV) <0.001 3.33 2.14 5.20 0.02 2.68 1.29 5.72
Epidermal growth factor receptor (WT/mutated) 0.42 0.79 0.44 1.39 0.07 0.55 0.28 1.06
KRAS proto-oncogene (WT/mutated) 0.39 1.23 0.76 2.02 0.73 1.10 0.64 1.89
Serine/threonine kinase 11 (WT/mutated) 0.21 0.72 0.44 1.18 0.14 0.63 0.34 1.16
Subtype (bronchioid/magnoid/squamoid) 0.07 1.69 1.01 2.83 0.55 0.81 0.42 1.57

HR, hazard ratio; CI, confidence interval; WT, wild-type.

Univariate and multivariate Cox regression analyses for risk score and eight clinical factors in the validation set.

Univariate Cox regression Multivariate Cox regression
Clinical factor P-value HR Lower 95% CI Upper 95% CI P-value HR Lower 95% CI Upper 95% CI
Risk score <0.001 248.55 25.21 2450.60 0.02 73.67 2.22 2445.24
Age 0.02 1.03 1.00 1.05 0.15 1.03 0.99 1.07
Sex (female/male) 0.07 1.64 0.96 2.78 0.18 1.67 0.79 3.51
Smoking status (smoker/nonsmoker) 0.96 0.98 0.44 2.15 0.87 1.21 0.12 11.89
Stage (I/II/III/IV) 0.03 3.16 1.06 10.76 0.16 4.50 0.99 23.49
Epidermal growth factor receptor (WT/mutated) 0.22 0.58 0.23 1.47 1.00 0.00 0.00 Infinity
KRAS proto-oncogene (WT/mutated) 0.41 1.30 0.70 2.41 0.93 0.95 0.28 3.24
Serine/threonine kinase 11 (WT/mutated) 0.40 0.63 0.19 2.03 0.81 0.84 0.21 3.36
Subtype (bronchioid/magnoid/squamoid) 0.84 1.13 0.63 2.06 0.46 0.67 0.23 1.94

HR, hazard ratio; CI, confidence interval; WT, wild-type.