Genomics Score Based on Genome-Wide Network Analysis for Prediction of Survival in Gastric Cancer: A Novel Prognostic Signature

Purpose Gastric cancer (GC) is a product of multiple genetic abnormalities, including genetic and epigenetic modifications. This study aimed to integrate various biomolecules, such as miRNAs, mRNA, and DNA methylation, into a genome-wide network and develop a nomogram for predicting the overall survival (OS) of GC. Materials and Methods A total of 329 GC cases, as a training cohort with a random of 150 examples included as a validation cohort, were screened from The Cancer Genome Atlas database. A genome-wide network was constructed based on a combination of univariate Cox regression and least absolute shrinkage and selection operator analyses, and a nomogram was established to predict 1-, 3-, and 5-year OS in the training cohort. The nomogram was then assessed in terms of calibration, discrimination, and clinical usefulness in the validation cohort. Afterward, in order to confirm the superiority of the whole gene network model and further reduce the biomarkers for the improvement of clinical usefulness, we also constructed eight other models according to the different combinations of miRNAs, mRNA, and DNA methylation sites and made corresponding comparisons. Finally, Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were also performed to describe the function of this genome-wide network. Results A multivariate analysis revealed a novel prognostic factor, a genomics score (GS) comprising seven miRNAs, eight mRNA, and 19 DNA methylation sites. In the validation cohort, comparing to patients with low GS, high-GS patients (HR, 12.886; P < 0.001) were significantly associated with increased all-cause mortality. Furthermore, after stratification of the TNM stage (I, II, III, and IV), there were significant differences revealed in the survival rates between the high-GS and low-GS groups as well (P < 0.001). The 1-, 3-, and 5-year C-index of whole genomics-based nomogram were 0.868, 0.895, and 0.928, respectively. The other models have comparable or relatively poor comprehensive performance, while they had fewer biomarkers. Besides that, DAVID 6.8 further revealed multiple molecules and pathways related to the genome-wide network, such as cytomembranes, cell cycle, and adipocytokine signaling. Conclusion We successfully developed a GS based on genome-wide network, which may represent a novel prognostic factor for GC. A combination of GS and TNM staging provides additional precision in stratifying patients with different OS prognoses, constituting a more comprehensive sub-typing system. This could potentially play an important role in future clinical practice.


Conclusion:
We successfully developed a GS based on genome-wide network, which may represent a novel prognostic factor for GC. A combination of GS and TNM staging provides additional precision in stratifying patients with different OS prognoses, constituting a more comprehensive sub-typing system. This could potentially play an important role in future clinical practice.
Keywords: gastric cancer, genome-wide network, miRNA, mRNA, DNA methylation, nomogram BACKGROUND Gastric cancer (GC) is one of the most common malignant human tumors and the third leading cause of cancer-related mortalities worldwide. Reports estimate that nearly one million new cases and 800,000 deaths occur each year across the world (Torre et al., 2015). Despite the rapid research advancement, GCrelated impacts on human life remain high around the globe. According to the global cancer burden data, hundreds of billions of dollars in economic losses are incurred each year due to GC. At the same time, stomach cancer has been reported to cause 19.1 million disability-adjusted life years, with 98% of these resulting from years of life lost and 2% from years lived with disability (Global Burden of Disease Cancer Collaboration et al., 2019).
Despite major breakthroughs in GC prevention, diagnosis, and treatment therapies reported over the past decade, prognosis remains a challenge at different TNM stages (Jiang et al., 2018a;Sun et al., 2019a,b). Notably, patients with similar clinical features and at the same tumor stage who receive uniform treatment have exhibited varying clinical outcomes (Bang et al., 2012;Jiang et al., 2018b). Such evidence indicates the existing challenges to traditional TNM staging (Serra et al., 2019), possibly due to a lack of molecular tools for effectively predicting the prognosis and the therapeutic effect of GC patients. Therefore, more rigorous and reliable systems that accurately reflect the heterogeneity of different patients and guide the development of treatment approaches are urgently needed (Duarte et al., 2018;Serra et al., 2019).
Tumors are a product of multiple genetic mutations, including genetic (gene expression) and epigenetic (DNA methylation and histone modification) modifications, as well as deregulations of tumor-suppressor genes and proto-oncogenes (Anna et al., 2018;Choi et al., 2019). In addition, changes in a set of genetic materials have been closely associated with cancer outcomes (Anna et al., 2018;Choi et al., 2019). Therefore, to effectively predict the prognosis of tumors, such as GC, a single biomarker is insufficient, necessitating the need for a gene network.
A variety of mRNAs have been associated with GC prognosis (Camargo et al., 2019), with microRNAs (miRNAs) also implicated in tumor prediction in the recent years (Li et al., 2010;Ueda et al., 2010;Camargo et al., 2019). These small, non-coding RNAs, comprising 22 nucleotides, primarily function to regulate protein translation by inhibiting the expression of target messenger RNAs (mRNAs). Apart from genetics, epigenetics is currently receiving considerable research attention. DNA methylation is the most common epigenetic event associated with cancer development and progression (Anna et al., 2018). Consequently, numerous studies have implicated DNA methylation in the diagnosis and the prognosis of various tumors, including GC (Camargo et al., 2019;Choi et al., 2019). Although these studies have revealed several biomarkers that have proved to be prognostic predictors in GC, only a handful have been adopted in clinical therapies or are used to build predictive models for the disease (Anna et al., 2018;Duarte et al., 2018;Camargo et al., 2019;Choi et al., 2019;Serra et al., 2019).
Previous studies have identified and recommended numerous biomarkers for GC. However, since malignant tumors often involve multiple layers and different levels of genetic changes, including the genome, transcriptome, and proteome, or even epigenetic content, selecting reasonable candidate factors from tens of thousands of biomarkers and comprehensively analyzing them as an independent feature is imperative to effectively develop a suitable prognostic target. Therefore, genetic networks containing a panel of abnormal factors from different regulatory levels represent the best chance for achieving prognostic value.
The whole genome-wide network analysis is reported in several other cancers, such as colorectal cancer, breast cancer, and lung cancer (Hou et al., 2018;Zhang et al., 2018), and it shows great value in differentiating the prognosis of these patients. Therefore, it is feasible and advantageous to apply genome-wide network analysis to GC.
In the current study, we performed a series of sophisticated statistical analyses and identified 33 genetic molecules that were highly correlated with the prognosis of GC. Specifically, we screened The Cancer Genome Atlas (TCGA), a genome project with 33 types of cancer, including gene expression, and DNA methylation as well as other biological information. Furthermore, we extended these independent prognostic factors to the "omics" concept. Then, a genome-wide network was constructed. Interestingly, the genomics score (GS) obtained herein could supplement TNM staging and enhance the prognostic value of different patients. Moreover, we developed multiple prognostic models, then validated, and compared them to ascertain their strengths and weaknesses. Finally, we performed pathway enrichment and gene oncology annotation analyses to elucidate the function of this gene network.

Data Acquisition and Preprocessing
Level 3 data were downloaded from the TCGA database using TCGA-Assembler Module A, in January 2019, which was then pretreated with Module B. The dataset comprised of clinical variables from 443 patients, including age, sex, stage, primary site, grade, treatment, and survival, as well as associated genome-wide data. In addition, the expression levels of 1,871 miRNAs, 20,531 mRNA, and 485,577 DNA methylation sites (Illumina methylation 450) were obtained from 384, 377, and 394 patients, respectively. Afterward, an intersection with a total of 332 samples was eventually retained. Furthermore, patients with missing active follow-up data were excluded from the analysis, leaving 329 patients in the final cohort (Figure 1). Moreover, genome-wide level 3 data whose expression levels for miRNAs, mRNA, and DNA methylation sites were missing in more than 50% of all samples were excluded from the final analysis. Finally, 329 GC patients with 566 miRNAs, 17,963 mRNA, and 396,081 DNA methylation sites were chosen for further analysis.

Genome-Wide Network Analysis
Gene expression and DNA methylation data were normalized using R package before subsequent processing. Univariate and least absolute shrinkage and selection operator (LASSO) Cox regression models were combined and used to identify the most useful prognostic factors in miRNAs, mRNA, and DNA methylation sites associated with survival. Firstly, univariate Cox regression was performed on each candidate miRNA, mRNA, and DNA methylation site to elucidate its role in patient survival, then signatures with P value less than 0.05 were retained for subsequent analysis. Thereafter, the LASSO Cox regression model was applied to select and shrink the data (Supplementary Figure S3; Tibshirani, 1997). Finally, a GS, based on a genomewide network comprising seven miRNAs, eight mRNAs, and 19 DNA methylation sites, was constructed for predicting survival. A summary of the whole screening process is displayed in Supplementary Figure S1.

Development and Comparison of Individualized Prediction Models
The TCGA cohort with 329 cases was used as the training set, with a random 150 cases from the total cohort included as a validation group. The random number is 1,356. Firstly, we developed the original GS based on 34 biomarkers (seven miRNAs, eight mRNAs, and 19 DNA methylation sites). Then, considering the complexity of the original GS and difficult clinical application, in order to obtain a more concise and FIGURE 1 | A Venn diagram displays the patients' screening process.
Frontiers in Genetics | www.frontiersin.org effective GS, we also constructed eight other models according to the different combinations of miRNAs, mRNA, and DNA methylation and made corresponding comparisons. Finally, a total of nine GS models based on the genome-wide network from LASSO were adopted to screen for the most appropriate markers. These included the following models: genomics (seven miRNAs, eight mRNA, and 19 DNA methylation sites), miRNAs (seven miRNAs), mRNA (eight mRNA), methylation (19 DNA methylation sites), miRNAs + methylation (seven miRNAs and 19 DNA methylation sites), miRNAs + mRNA (seven miRNAs and eight mRNA), mRNA + methylation (eight mRNA and 19 DNA methylation sites), Cox-model 1 (two miRNAs, six mRNA, and nine DNA methylation sites), and Cox-model 2 (one miRNA, one mRNA, and seven DNA methylation sites). Among them, markers from Cox-model 1 were separately detected from miRNAs, mRNA, or DNA methylation sites using multivariate Cox regression analysis after LASSO (Supplementary Tables S2-S4). On the other hand, markers from Cox-model 2 depended on signatures from a multivariate Cox regression analysis combining the genomewide network and the clinical characteristics (Supplementary Table S5). Thereafter, we constructed several nomograms by incorporating significant (P < 0.05) GS variables and other clinical features following multivariate Cox regression (Iasonos et al., 2008), and a clinical nomogram was built as a blank control. The equations used for calculating the GS of these models are listed in Supplementary Table S6.
To calculate the discrimination and the stability of different Cox regression models, we applied C-statistics and calibration. Additionally, we performed an analysis of time-dependent receiver operator characteristics (ROC), based on the 1-, 3-, and 5-year survival endpoints, to assess the prognostic accuracy of the different nomograms. Furthermore, we evaluated the potential net benefit of different predictive models using decision curve analysis (DCA). DCA compares the clinical usefulness of different indicators by calculating the potential net benefit of each decision strategy at each threshold probability. Thus, DCA was a significant novel approach for comparing the old and the new models (Vickers and Elkin, 2006).

Screening for Potential miRNA Target Genes
We predicted the potential target genes of the seven miRNAs, from LASSO, by screening the miRTarBase, miRDB, and TargetScan databases. Common genes from each miRNA across the three databases were then used for subsequent studies. More than 90% of the miRNAs showed negative regulation to target genes. Consequently, the expression data from TCGA were used to perform a batch of correlation analysis of each miRNA, with corresponding target genes, and the three genes with the largest absolute negative correlation were retained as the most likely targets. Additionally, at least three potential target genes from miRTarBase, which is co-expressed with miRNAs, were considered as equally important and were subjected to Cytoscape (version 3.7.2) for identification of miRNA-target genes coexpression network analysis (Supplementary Figure S2).

Functional Enrichment Analysis
The potential target genes that were negatively correlated with miRNAs in TCGA, as well as the coding sequences for mRNA and DNA methylation sites, were used for functional enrichment analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) using DAVID 6.8 (Supplementary Figure S2). Functional enrichment analysis indicates why the gene network produces images on the survival of GC from a molecular mechanism. Visualization was then done using the "ggplot2" package implemented in R.

Statistical Analysis
The patients were divided into low-risk and high-risk groups by the median GS as the cutoff point. Survival estimates were obtained according to the Kaplan-Meier method and compared using the log-rank test. Variables that reached significance, with P < 0.05, were entered into the multivariable analyses using the Cox proportional hazards model, with an entry stepwise approach to identify covariates associated with increased allcause mortality, and then hazard ratio with 95% confidence intervals (CIs) of each variable was achieved. All the statistical significance values were set as two-sided (P < 0.05). LASSO Cox regression was performed through the "glmnet" package. Time-dependent ROC analysis at different follow-up times was implemented using the "timeROC" package of R project in order to further expound the performance of different GS models, and DCA was used to compare their clinical use by "rmda" package. Finally, nomogram based on the Cox regression model was constructed using the "rms" package. C-index and calibration to calculate the discrimination and the stability of these models were performed using c-statistics and Bootstrap sample. Harrell's concordance index (C-index) indicated a better prognostic model if its value was closer to 1, and the calibration diagram showed that the better the prediction if the closer the correction line was to the diagonal. All statistical methods are applied to both the training group and the validation group. Statistical analyses were performed using SPSS statistical software (version 18.0) and R software (version 3.5.3).

Patient Characteristics
Among the 329 GC patients analyzed in this study, 212 (64.4%) were male, whereas 117 (35.6%) were female. The average age of the entire study population was 65.0 ± 10.6 years. In terms of pathological stage, 38 (11.6%) cases were identified as stage I, 117 (35.6%) were stage II, 155 (47.1%) were stage III, and 19 (5.8%) were at stage IV. With regards to treatment, 303 (92.1%) patients received surgery (280 cases of R0 surgery, 14 R1, and nine R2), whereas 146 (44.4%) were subjected to fluorouracil-based chemotherapy. The genomic nomogram classified 165 samples into low GS (GS ≤ -0.137) and 164 into high GS (GS > -0.137) groups based on the median cutoff (Figure 2). A detailed description of tumor location, pathology grade, and Lauren classification is outlined in Supplementary Table S1, while a heat map of the genomic scores layered by clinicopathological  Figure S5). Toward the last follow-up, a total of 129 deaths and 200 censoring were recorded. The estimated cumulative 1-, 3-, and 5-year OS in the total cohort were 78.9, 48.4, and 36.7%, respectively, although these rates increase to 95.1, 74.2, and 68.5%, respectively, in the low-GS group. Conversely, the 1-, 3-, and 5-year OS decreased to 63.3, 23.6, and 15.4%, respectively, in the high-GS group. The baseline information of the validation cohort is also listed in Supplementary Table S1 and Supplementary Figure S6.

Survival Analysis
We identified a basic genome-wide network comprising seven miRNAs, eight mRNAs, and 19 DNA methylation sites as the prognostic factor for OS, from hundreds of thousands of univariate Cox regression and LASSO analyses. This network was then classified as other models in the training and the validation groups. Among the 34 features identified, poor prognosis was significantly associated with a high expression of seven miRNAs (hsa-mir-100, hsa-mir-1304, hsa-mir-136, hsa-mir-193b, hsa-mir-22, hsa-mir-653, and hsa-mir-6808), six mRNAs (NRP1|8829, RNF144A|9781, ZNF22|7570, DUSP1|1843, CPNE8|144402, MAGED1|9500, and LOC91450|91450), and seven DNA methylation sites (cg07020967, cg08859156, cg12485556, cg15861578, cg15861578, cg25161386, and cg22740006). Conversely, poor prognosis was strongly associated with a low expression of SOX14|8403 and 12 DNA methylation sites, including cg02223323, cg00481239, cg14791193, cg15486740, cg20100408, cg20350671, cg22395807, cg24361571, cg25361506, cg22813794, cg26014401, and cg26856948 ( Table 2). Univariate analysis performed on clinical characteristics revealed a significant association between age, pathological stage, TNM, and surgery with OS (Table 1). On the other hand, results from multivariable Cox regression showed that age, pathological stage, and GS were significantly associated with all-cause mortality in GC (Table 1 and Figure 3). Furthermore, stratification of the pathological stage (I, II, III, and IV) revealed significant differences in survival rates between the high-GS and the low-GS groups (Figure 3). A similar result was found when the data were stratified by demographic variables (sex and age), clinical characteristics (primary site, grade, and Lauren classification) as well as treatments (surgery and chemotherapy; Supplementary Figures S7, S8). On the other hand, categorizing GS into high or low groups, using the median value across different models, indicated that the genomics nomogram had the highest HR value. Interestingly, HR was almost equal to miRNAs + methylation, mRNA + methylation, Cox-model 1, and Cox-model 2 nomograms, which contained fewer gene features. Moreover, the HR value showed a marked decrease in miRNAs, mRNA, methylation, and miRNAs + mRNA nomograms, which included the least characteristics (Table 3).

Nomograms Based on Genome-Wide Network
A genomics nomogram was first constructed based on the genome-wide network, comprising 34 gene features (Figure 4). To obtain a more concise and effective nomogram, we also built a Cox-model 1 (17 gene features) and Cox-model 2 (nine gene features) nomograms (Supplementary Figures S9, S10). Next, a clinical nomogram, based on stage and age, was built as a control (Supplementary Figure S11). Thereafter, we performed internal and external validation to evaluate the feasibility of all nomograms using a three-grouped random bootstrap sampling ( Figure 5 and Supplementary Figures S9-S11). We observed good predictive performance in the first three nomograms, but not in the simple clinical model.
Additionally, DCA showed that the genomics, Cox-model 1, mRNA + methylation, and methylation nomograms had a significant net benefit compared to other GS models and the clinical nomogram (Figures 6C,D).

Functional Analysis of Genome-Wide Network
We imported the 301 potential target genes, mRNA, and DNA methylation site-coding sequences, identified above, into DAVID for KEGG and GO analyses and identified biological processes, molecular functions as well as cellular components (Figures 7A-C). Their corresponding KEGG pathways were also plotted ( Figure 7D).

DISCUSSION
GC can be divided into two types or four main categories, according to the Lauren and World Health Organization (WHO) classifications (Lauren, 1965;Nagtegaal et al., 2019), although neither of these classifications is based on molecular markers. In the last decade, however, three novel molecularbased classification systems have been suggested for GC. The Singapore-Duke Group was the first to describe a classification with two intrinsic genomic subtypes, G-INT, and G-DIF, which had different gene expression (Tan et al., 2011;Serra et al., 2019). The subtypes have different levels of resistance to various chemotherapy drugs and show limited prognostic value. Later, TCGA used molecular evaluation to propose a new classification with four subtypes: EBV, MSI, GS, and CIN. The identification of these subtypes has provided a roadmap for patient stratification  as well as targeted therapeutic trials (Cancer Genome Atlas Research, 2014). However, initial data on disease outcomes from this cohort did not show differences in survival among the four groups. A series of positive studies on prognosis based on TCGA classification was also reported (Sohn et al., 2017). In addition, the Asian Cancer Research group divided GC into four subtypes, MSI, EMT, MSS/TP53+, and MSS/ TP53-, based on gene expression data and found significantly different survival outcomes across them (Cristescu et al., 2015;Serra et al., 2019). Despite the significant milestones of these studies, they are all mainly based on the analysis of gene expression (mRNA). Besides that, a 2019 study proposed a five-miRNA model, while it had a C-index of 0.72 only ). In the current study, we included methylation data and performed functional enrichment analysis, making our work stronger. The aforementioned classifications are also complicated and need further optimization to increase clinical applicability.
Currently, focus has been directed on identifying prognostic miRNAs for GC. Particularly, one miRNA can regulate multiple targets, while multiple miRNAs can regulate a single mRNA. Therefore, a single miRNA may play an opposite role in cancer progression by regulating different target genes. For example, Mir-22 and Mir-100 were found to be tumor suppressors in various cancers, including GC (Chen et al., 2014;Zuo et al., 2015). Similarly, a high expression of Mir-136 was found to promote proliferation and invasion in GC cell lines by inhibiting PTEN expression , while a contrasting result was reported when HOXC10 was targeted (Zheng et al., 2017). Similarly, Mir-193b reportedly induced GC proliferation or apoptosis by mediating different mRNA expressions (Mu et al., 2014;Song et al., 2018), whereas a high Mir-1304 expression in GC was reported as a negative predictor for prognosis of lung and thyroid cancers (Kurata and Lin, 2018;Liu et al., 2019;Pan et al., 2019). However, the function of Mir-653 and Mir-6808 has not been previously reported. In the current study, we found an association between a high expression of all miRNAs and poor survival. Different outcomes may be observed in our study, relative to previous reports, owing to the huge number of corresponding miRNA target genes herein and lack of evidence on their role in GC development. Messenger RNAs have been reported to play an essential role in GC cancer. For example, high NRP1 expression and hypermethylation were associated with poor GC prognosis (Wang et al., 2019), whereas another study indicated that it could be an anti-tumor target (Grandclement and Borg, 2011). In addition, high DUSP1 expression levels were found to promote progression, drug resistance, and poor prognosis of GC (Teng et al., 2018). On the other hand, SOX14 showed opposite prognostic values in cervical cancer and leukemia, with anti-tumor and carcinogenic effects, respectively (Stanisavljevic et al., 2017;Tosic et al., 2018). Studies have also implicated CPNE8, MAGED1, and RNF144A in ovarian and breast cancers (Zeng et al., 2012;Nagasawa et al., 2019;Yang et al., 2019). However, LOC91450 and ZNF22 have not been reported in cancer.
Accumulating evidence indicates that DNA methylation plays a significant role in cancer progression. However, only a handful of studies have described the relationship between levels of FIGURE 4 | Genomics nomogram to predict the probability of 1-, 3-, and 5-year overall survival (OS) in the training cohort (A) and the validation cohort (B): to determine how many points for each variable to the probability of OS, locate the variable on its axis, draw a line straight upward to the point axis, repeat this process for each variable, sum up the points achieved for each of the risk factors, locate the final sum on the total point axis, and draw a line straight down to find the patient's probability of OS.
single-site methylation and GC prognosis. Particularly, high expressions of MAP7D2, ACOT13, EID1, RPS4X, and TTRAP have been associated with poor prognosis in gastric, lung, and pancreatic cancers as well as hepatic carcinoma, respectively, while a high TTRAP expression reportedly inhibits the growth of osteosarcoma (Kamio et al., 2010;Zhou et al., 2013;Kuang et al., 2017;Liu K.T. et al., 2018;Liu X. et al., 2018). Notably, the relationship between methylation levels and corresponding gene expression profiles is unknown, necessitating further research. Furthermore, the remaining DNA methylation sites and their corresponding genes have not been reported. Lastly, no study has described the prognostic significance using a genomewide network.
Last but not least, in general, no study concerning their prognostic significance as a genome-wide network has been reported yet.
Tumorigenesis involves multiple interacting biological processes. In addition, an integrated genetic network is better at reflecting intra-tumor heterogeneity compared to a single biomarker. In the current study, we identified a novel, prognostic, signature genome-wide network, consisting of seven miRNAs, eight mRNA, and 19 DNA methylation sites after screening the entire TCGA cohort using training and random cohorts. This network was further divided into several other models.
Our results revealed that the integrative signature was an independent prognostic factor for survival in GC patients and performed better than any single biomarker or clinical characteristic. Moreover, stratification by other clinicopathological features, such as stage, age, sex, primary site, pathology grade, Lauren classification, and treatments, revealed significantly different prognosis values based on different GSs. In addition, staging was still an effective prognostic factor after dividing into low-and high-genomics-score groups, suggesting that GS and traditional staging can complement each other, and the genetic network could add prognostic value to traditional staging. Exclusion of patients with I staging showed that chemotherapy is a significant prognostic factor because I staging does not always need additional chemotherapy for effective prognosis.
We also developed and validated nomograms based on the GS. Particularly, results from ROC and DCA indicated that all of them had significantly better predictive performances than the traditional clinical nomogram. Comprehensive property (similar C-index) was not significantly different in genomics nomogram and Cox-model 1 nomogram, and compared to the genomics nomogram, Cox-model 1 nomogram had fewer biomarkers. In addition, Cox-model 1 nomogram performed well, with a higher positive net reclassification improvement (NRI). Therefore, Cox-model 1 nomogram might be more suitable for clinical application, which deserved further study. Besides that, the Cox-model 2 nomogram had the least feature (nine biomarkers) including miRNAs, mRNA, and DNA methylation sites for constructing a genome-wide network, while it had a lower C-index and a negative NRI. The other six   models showed a relatively poor performance in ROC or DCA, with limited application value. What is more, it is possible that DNA methylation was the highest contributor to the survival prediction of this gene network. We suspect that this may be related to the larger number of DNA methylation sites compared to miRNA and mRNA. We adopted GO and KEGG analyses to assess the influence of genome-wide network in the prognosis of GC. Generally, biological processes mainly involve various biological functions, such as methylation, phosphorylation, and endocrine regulation. Methylation pathway was related to the occurrence and the development of GC, which was consistent with our results. Besides that, functional enrichment analysis revealed that phosphorylation pathway was significantly enriched as well, which got more and more attention these years. On the other hand, the main components of participation included organelles, cytomembranes, extrinsic to membranes, nuclear and synapses, whereas molecular functions comprise nucleoside, ATP, RNA, and transcription factor binding as well as activity of various enzymes. Abnormal cell composition is closely related to the development of tumor. The abnormal protein may act on the nucleus, membrane, or cell matrix, thereby leading to the progression of cancer, such as NRP1 protein (Wang et al., 2019). In the current study, KEGG analysis indicated that the gene network function was a relevant pathway in cancer, cell cycle, and adipocytokine signaling, while the other pathways had been reported in small cell lung and bladder cancers. Further experiments to reveal the biological function of this gene network are needed.
We also employed a series of complex statistical analyses to construct and validate a genome-wide network based on different biomarkers and then divided it into different models. We recommend the resulting GS despite it not being an absolute representative of tumor heterogeneity. This network could complement the deficiency of traditional staging and generate a more accurate prediction of survival rates in GC patients. Additionally, it effectively distinguishes patients who could benefit from chemotherapy, thereby reducing unnecessary treatments. It is also possible that the network could be used to identify novel therapeutic targets for GC, although this requires further investigation.

Limitation
This study had several limitations. Firstly, information relating to patient co-morbidities and performance status was not available in the TCGA database. Secondly, the systemic chemotherapy regimens were not uniform, and most of them were based on fluoropyrimidines. Thirdly, the gene network contains too many biomarkers, increasing the difficulty of clinical use. Lastly, this was a retrospective study, without any independent external patient datasets as test. Despite some limitations, it was the first, to the best of our knowledge, to integrate miRNAs, mRNA, and DNA methylation sites as a genome-wide network to predict the OS of patients with GC, and we would try to design a validation in our hospital.

CONCLUSION
In summary, we used a TCGA cohort to develop and validate a novel genome-wide network comprising seven miRNAs, eight mRNAs, and 19 DNA methylation sites for the prognosis of GC. A combination of GS and TNM staging enhances its prognostic value, proposing a more comprehensive sub-typing system. The developed network is expected to aid in predicting GC patients who may benefit from chemotherapy to some degree.

DATA AVAILABILITY STATEMENT
The datasets generated and analyzed during the current study are available in the TCGA database [https://portal.gdc.cancer.gov/].

ETHICS STATEMENT
Since TCGA was a public-use database, no additional permission was required from the Ethics Committee. In addition, this study was deemed exempt from institutional review board approval by Nanfang Hospital of Southern Medical University (Guangzhou, China).

AUTHOR CONTRIBUTIONS
All the authors listed had made a substantial contribution to this work. YJ and GL put forward the conception and designed the study. WH, TL, and MZ collected and collated the data. ZS, HC, and ZH analyzed the data and wrote the manuscript together. HL, JY, and YH made contribution to proofread the article. Finally, all the authors took responsibility of the final manuscript and approved it for publication.

ACKNOWLEDGMENTS
We would like to thank the staff members of TCGA program.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020. 00835/full#supplementary-material FIGURE S1 | Screening process for the genome-wide network.     FIGURE S7 | Kaplan-Meier curve of overall survival for the low-and the high-genomics-score groups divided by age (young and old), sex (male and female), primary site (cardia, fundus/body, and antrum), pathology grade (I-II and III-IV), Lauren classification (intestinal type and diffused type), and treatment (R0 surgery, R1 surgery, chemotherapy, and no chemotherapy) in the training group.
FIGURE S8 | Kaplan-Meier curve of overall survival for the low-and the high-genomics-score groups divided by age (young and old), sex (male and female), primary site (cardia, fundus/body, and antrum), pathology grade (I-II and III-IV), Lauren classification (intestinal type and diffused type), and treatment (R0 surgery, R1 surgery, chemotherapy, and no chemotherapy) in the validation group.
FIGURE S9 | Cox-model 1 nomogram and its calibration plot in the training group and the validation group.
FIGURE S10 | Cox-model 2 nomogram and its calibration plot in the training group and the validation group.
FIGURE S11 | Clinical nomogram and its calibration plot in the training group and the validation group.
FIGURE S12 | Time-dependent receiver operating characteristic curves on 1, 3, and 5 years for each nomogram in the training group.
FIGURE S13 | Time-dependent receiver operating characteristic curves on 1, 3, and 5 years for each nomogram in the validation group.
FIGURE S14 | Venn diagrams for the target genes of each miRNA.
FIGURE S17 | Kaplan-Meier curve of overall survival in the low-genomics-score group stratified by clinical features (e.g., stage).
FIGURE S18 | Kaplan-Meier curve of overall survival in the high-genomics-score group stratified by clinical features (e.g., stage).
FIGURE S19 | Kaplan-Meier curve of overall survival in stages II and III patients.