Immune infiltration-related N6-methyladenosine RNA methylation regulators influence the malignancy and prognosis of endometrial cancer

m6A regulator gene mutations and copy number variation in endometrial cancer patients

We obtained somatic mutation and CNV data of EC patients from the UCSC Xena database and analyzed changes in m6A regulatory genes. Among the EC samples from 433 cases, mutations in m6A regulatory genes were detected in 133 independent samples (Supplementary Table 4). Of 529 EC samples with CNV data, CNVs were frequently observed in 17 m6A regulator-encoding genes (Figure 1A). Among these, the frequency of CNV was highest in the m6A writer gene KIAA1429 (36.11%, 191/529, Table 1), followed by reader gene YTHDF3 (34.22%, 181/529, Table 1). We then evaluated the CNV pattern in EC samples and found that most CNV events were increases in copy number (1144/1086; Figure 1B). Finally, we determined the most common CNV type in m6A regulatory genes. Copy number gain of KIAA1429 was the most frequent alteration among all m6A regulatory gene CNVs (Figure 1C).

Figure 1. Copy number variation of m6A regulatory genes in endometrial cancer. (A) Analysis of the percentage of CNVs among m6A regulatory genes in TCGA data of 529 EC samples. (B) Incidence of m6A regulatory gene copy number increase or loss in EC samples. (C) The most common CNV among m6A regulatory genes in EC samples.

m6A gene mutations are associated with endometrial cancer clinicopathological and molecular characteristics

We assessed the relationship between alterations (mutation and/or CNV) in m6A regulatory genes and clinicopathological characteristics in 433 EC patients. The alterations of m6A regulatory genes were related to clinical stage, grade, infiltration degree, and lymph node metastasis in patients (Table 2). Since POLE, PTEN, PIK3CA, KRAS, PIK3R1, TP53, FBXW7, and PPP2R1A play an important role in EC pathogenesis, we further assessed whether mutations in m6A regulators were related to the changes in these eight genes. Changes in m6A regulatory genes were significantly associated with POLE, PTEN, KRAS, TP53, and PPP2R1A alterations Next, we analyzed the influence of different CNV types on the mRNA expression of m6A regulatory genes. Indeed, expression was significantly correlated with different CNV patterns in 529 EC samples. For all 17 genes, the gene copy number increase (amplification) was related to enhanced mRNA expression, while shallow or deep deletions reduced mRNA expression (Supplementary Figure 1).

Correlation between m6A gene copy number variants and survival of endometrial cancer patients

To explore the prognostic value of CNVs in m6A regulatory genes, the effects of CNVs on the overall survival (OS) and disease-free survival (DFS) of 412 patients with EC were analyzed. Different CNV types in m6A regulatory genes were found to be related to OS and DFS (Figure 2A-i, Figure 2B-i). A separate analysis of the 17 genes revealed that patients with deep or shallow deletions of IGF2BP1, ZC3H13, METTL14, ALKBH5, or YTHDC1 had lower OS and DFS than those who did not (Figure 2A-ii, Figure 2B-ii; and Supplementary Figures 2, 3). Univariate Cox analysis indicated that age, stage, grade, invasion, and lymph metastasis were independent risk factors for OS. The somatic mutation of POLE, PTEN, PIK3CA, KRAS, PIK3R1, TP53, FBXW7, and PPP2R1A genes were not significantly correlated with the survival rate of patients. Further, loss of eraser and writer genes was found to be an independent risk factor for DFS and OS (Table 3). Patients in the eraser loss (+) + and writer loss (+) groups had the shortest OS and DFS (Figure 2C). These results supported a relationship between decreased m6A regulator levels and poor survival.

Figure 2. Relationship between copy number variants of m6A regulatory genes and OS/DFS in 412 EC patients (A-i) OS of EC patients with or without mutations in m6A regulatory genes. (A-ii) OS of EC patients with different CNV types for IGF2BP1, ZC3H13, METTL14, ALKBH5, and YTHDC1. (B-i) DFS of EC patients with or without mutations in m6A regulatory genes. (B-ii) DFS of EC patients with different CNV types for IGF2BP1, ZC3H13, METTL14, ALKBH5, and YTHDC1. (C) OS and DFS of EC patients with simultaneous alterations in writer and eraser genes.

m6A regulates the expression of key gene mRNAs in endometrial cancer patients

TCGA data was used to analyze the expression of IGF2BP1, ZC3H13, METTL14, ALKBH5, and YTHDC1 mRNA in 35 normal and 546 EC tissues. The UALCAN online tool (http://ualcan.path.uab.edu/index.html) was employed to assess corresponding protein expression in 31 normal and 100 EC tissues. ZC3H13, METTL14, ALKBH5, YTHDC1, and IGF2BP1 expression were downregulated. We performed univariate survival analysis of the EC TCGA data to determine the relationship between the expression of these five genes and EC prognosis. High expression of IGF2BP1 was associated with unfavorable prognosis, while high expression of ZC3H13, METTL14, and ALKBH5 was associated with improved prognosis (Figure 3A–3E).

Figure 3. Candidate m6A regulator protein and mRNA expression. (A–E) TCGA and UALCAN website analysis of IGF2BP1, ZC3H13, METTL14, ALKBH5, and YTHDC1 mRNA expression in 35 normal and 546 EC tissues, protein expression in 31 normal and 100 EC tissues, and Kaplan-Meier survival curve of the prognostic signature based on five m6A-related genes. *P < 0.05, ** P < 0.01, *** P < 0.001. (F) Expression of the five genes in 34 EC tissues and 20 normal tissues, as determined by qRT-PCR. (G) Protein expression of ZC3H13, METTL14, and YTHDC1 in 35 EC tissues and 20 normal tissues was determined by immunohistochemistry. DFS curves for ZC3H13, METTL14, and YTHDC1 protein expression in 35 endometrial carcinoma cases. *P < 0.05, ** P < 0.01, *** P < 0.001.

We then detected the expression of IGF2BP1, ZC3H13, METTL14, ALKBH5, and YTHDC1 in 20 normal endometrial tissues and 34 EC tissues via reverse transcription polymerase chain reaction (RT-PCR). ZC3H13, METTL14, and YTHDC1 expression was lower in EC tissues (Figure 3F). The protein expression of ZC3H13, METTL14, and YTHDC1 was quantified in 20 normal endometrial tissues and 35 EC tissues using immunohistochemistry. All three proteins were primarily localized within the nucleus, and positive expression was indicated by brownish yellow or uniform brown granules. Moreover, the expression rates of ZC3H13, METTL14, and YTHDC1 in the EC group [17.1% (6/35); 48.6% (17/35); and 40.0% (14/35), respectively] were significantly lower than those in the normal control group [55.0% (11/20); 70.0% (14/20); and 60.0% (12/20), respectively; P < 0.05]. We then analyzed the association of ZC3H13, METTL14, and YTHDC1 protein expression with DFS of EC patients. High expression of ZC3H13, METTL14, and YTHDC1 correlated with a higher survival rate. However, these associations were not statistically significant (Figure 3G).

Prognostic markers related to the immune microenvironment of endometrial cancer

The TIMER database was used to analyze the correlation between EC immune infiltration level and the expression, mutation, and somatic copy number alterations (SCNAs) of METTL14, ZC3H13, and YTHDC1 in 545 EC tissues. METTL14, ZC3H13, and YTHDC1 expression were significantly related to immune cell infiltration in various tumors. In EC, the expression of METTL14 was positively correlated with CD4+ T cell infiltration and negatively correlated with macrophage infiltration (Figure 4A). ZC3H13 and YTHDC1 expression were positively correlated with CD4+ T cell infiltration, while being negatively correlated with macrophage and NK cell infiltration (Figures 5A, 6A). METTL14, ZC3H13, and YTHDC1 mutations were prevalent in EC (Supplementary Figures 4A, 5A, 6A). Violin chart analysis indicated that METTL14 mutations were associated with CD8+ T cell, CD4+, T cell, and B cell infiltration (Figure 4B). ZC3H13 mutations were associated with CD8+ T cell, CD4+ T cell, B cell, neutrophil, macrophage, and NK cell infiltration (Figure 5B). YTHDC1 mutations were associated with CD8+ T cell, CD4+ T cell, Treg, B cell, macrophage, and NK cell infiltration (Figure 6B). The SCNA status of METTL14 was associated with CD8+ T cell, CD4+ T cell, Treg, macrophage, and NK cell infiltration (Figure 4C). ZC3H13 SCNA was associated with CD8+ T cell infiltration (Figure 5C). The SCNAs of YTHDC1 had no statistically significant relationship with immune cell infiltration (Figure 6C). Supplementary Figures 4B, 5B, 6B depict the SCNA status of METTL14, ZC3H13, and YTHDC1 in diverse cancer types. The outcome module indicated the relevance of infiltrating immune subsets to tumor clinical stage (Figures 4D, 5D, 6D). Spearman analysis revealed that, in various cancers, the expression of METTL14, ZC3H13, and YTHDC1 was positively correlated with programmed death-ligand 1 (PD-L1) expression (Figure 7A). The same was observed for EC (Figure 7B).

Figure 4. Correlation between immune infiltration and the expression, mutation, SCNA status, and outcome module of METTL14. (A) Heatmap depicting the correlation of METTL14 expression with six tumor-infiltrating immune cell types (CD4+ T cells, CD8+ T cells, Treg, B cells, neutrophils, macrophages, and NK cells) and the immune infiltration level in diverse cancer types. Spearman’s correlation was used for this analysis. (B) Violin plots to visualize the effect of METTL14 gene mutations on immune cell infiltration and different infiltrating immune cell types in endometrial cancer. (C) Violin plots visualize the effect of METTL14 SCNA, including ‘deep deletion’, ‘arm-level deletion’, ‘diploid/normal’, ‘arm-level gain’, and ‘high amplification’, on immune cell infiltration and different immune cell types in endometrial cancer, as determined by GISTIC2.0. (D) Outcome module showing the clinical stage relevance of tumor immune subsets as well as the hazard ratio and P-value for the Cox model. The log-rank P-value for the KM curve is shown on the KM curve plot.

Figure 5. Correlation between immune infiltration level and the expression, mutation, SCNA status, and outcome module of ZC3H13. (A) Heatmap depicting the correlation of ZC3H13 expression with immune infiltration level in diverse cancer types. (B) Violin plots examining the effect of ZC3H13 gene mutations on immune cell infiltration and immune cell types in endometrial cancer. (C) Violin plots depicting the effect of ZC3H13 gene SCNA status on immune cell infiltration and immune cell types in endometrial cancer. (D) Outcome module showing the clinical relevance of tumor immune subsets.

Figure 6. Correlation between the immune infiltration level and the expression, mutation, SCNA status, and outcome module of YTHDC1. (A) Heatmap depicting the correlation of YTHDC1 expression with the immune infiltration level in diverse cancer types. (B) Violin plots depicting the effect of YTHDC1 gene mutations on immune cell infiltration and immune cell types in endometrial cancer. (C) Violin plots depicting the effect of YTHDC1 SCNA status on immune cell infiltration and immune cell types in endometrial cancer. (D) Outcome module showing the clinical relevance of tumor immune subsets.

Figure 7. The correlation of METTL14, ZC3H13, and YTHDC1 expression with PD-L1. (A) Heatmap depicting the correlation of METTL14, ZC3H13, and YTHDC1 expression with PD-L1 in diverse cancer types. (B) Spearman’s correlation analysis was employed.

ZC3H13 and YTHDC1 knockdown promotes the proliferation and invasion of Ishikawa and HEC-1A cells

EC is divided into two types, of which type I is estrogen-dependent, while type II is non-estrogen-dependent. The Ishikawa cell line is derived from type I tumors, and HEC-1A cells derived from type II tumors [21]. To investigate the carcinogenic effect of ZC3H13 and YTHDC1 in endometrial carcinoma, we knocked them down in Ishikawa and HEC-1A cell lines. The transfection efficiency of ZC3H13- and YTHDC1-targeting small interfering RNA (siRNA) was determined via quantitative RT-PCR (qRT-PCR; Figure 8A). ZC3H13 and YTHDC1 knockdown both promoted the proliferation of EC cell lines, as determined through EdU assays (Figure 8B). Furthermore, both si-ZC3H13- and si-YTHDC1-transfected cells exhibited stronger invasive ability compared to control cells (Figure 8C).

Figure 8. Knockdown of ZC3H13 or YTHDC1 promotes the malignant behavior of Ishikawa and HEC-1A cells. (A) qRT-PCR was used to determine transfection efficiency. (B) The effect of ZC3H13 or YTHDC1 expression on the proliferation of Ishikawa and HEC-1A cells was determined via EdU assays. (C) Transwell assays were used to determine the number of invading cells. Data are presented as the mean ± SEM (n = 3 per group), *P < 0.05, ** P < 0.01, *** P < 0.001 vs. the NC group.

Potential regulatory mechanisms of ZC3H13 and YTHDC1

In view of ZC3H13 and YTHDC1 suppression promoting proliferation and invasion, we sought to further explore the underlying molecular mechanisms of both in EC. To this end, we downloaded lncRNA and miRNA data of 529 EC samples and 35 normal samples from the TCGA and cBioPortal for Cancer Genomics databases. A total of 846 differentially expressed genes were identified for ZC3H13, including 146 upregulated and 700 downregulated genes, in addition to 55 differentially expressed non-coding transcripts, 8 of which were upregulated and 47 downregulated. A total of 1107 differentially expressed genes were identified for YTHDC1, comprising 177 upregulated genes and 930 downregulated genes. In addition, 66 differentially expressed non-coding transcripts were identified, including 5 upregulated and 61 downregulated molecules. The bar plot depicts significantly enriched GO terms in patients with low ZC3H13 and YTHDC1 expression. Functional enrichment analyses revealed significantly enriched terms, including the humoral immune response and B cell-mediated immunity (Figure 9A–9D). We constructed a competing endogenous RNA (ceRNA) regulatory network to explain the potential regulatory mechanism of low ZC3H13 and YTHDC1 expression in EC. The ceRNA network related to ZC3H13 expression included 3 lncRNAs, 124 miRNAs, and 35 mRNAs, while that related to YTHDC1 expression included 1 lncRNA, 121 miRNAs, and 125 mRNAs (Figure 9E, 9F). There were 1301 differentially expressed genes associated with METTL14, including 186 upregulated genes and 1115 downregulated genes, in addition to 143 differentially expressed non-coding transcripts, of which 5 were upregulated and 138 were downregulated (Supplementary Figure 7A). Functional enrichment results of low METTL14 expression in EC revealed several immune-related GO terms (Supplementary Figure 7B). The ceRNA regulatory network related to METTL14 expression included 3 lncRNAs, 12 miRNAs, and 8 mRNAs (Supplementary Figure 7C).

Figure 9. Potential regulatory mechanism of ZC3H13 and YTHDC1 dysregulation in endometrial cancer. (A) Volcano map showing the differential mRNA, microRNA, and lncRNA expression in patients with low ZC3H13 expression. An absolute log2-fold change (FC) > 1.5 and FDR-adjusted P < 0.01 were used as the cutoff for significantly differentially expressed mRNAs and lncRNAs. Red represents significantly upregulated genes. Blue represents genes that are significantly downregulated. Grey represents genes that are not differentially expressed. (B) The bar graph shows that ZC3H13-regulated expression was enriched for various GO terms. (C) Volcano plot showing differentially expressed transcripts in patients with low YTHDC1 expression. (D) The bar graph shows that YTHDC1-regulated expression was enriched for various GO terms. (E, F) The ceRNA regulatory network in patients with low expression of ZC3H13 and YTHDC1. Red indicates lncRNA; yellow indicates mRNA; blue indicates microRNA.

Ligand preparation and molecular docking

Finally, we explored drugs targeting METTL14 and YTHDC1 with potential therapeutic effects via virtual screening. In particular, we analyzed the molecular docking binding force of small-molecule compounds. The RNA-binding domain of METTL14 lies within both METTL3 and METTL14 [22]. We targeted this domain of the human METTL14 protein (key amino acids: R245/R249/R254/R255/K297/R298). Studies have reported the co-crystal structure of human YTHDC1 (PDB ID: 4R3I) and m6A, detailing the conserved m6A-binding region of YTHDC1 and the key amino acid residues (N363/N367/W377/S378) [23], which were used in the present study. From the resulting 2D and 3D model diagrams, we identified the five most relevant compounds (Figure 10).

Figure 10. The 2D and 3D structures of five top-scoring small-molecule drugs. (A) Top five compounds acting on the RNA-binding domain of METTL14. (B–F) 2D and 3D mapping of the binding interaction between the top five compounds (Z1791109186, Z1445387105, Z89321655, HY-N0055, and HY-102053) and METTL14. (G) The five preferred compounds targeting the m6A-binding domain of the YTHDC1 protein. (H–L) 2D and 3D mapping of the binding between the top five compounds (Z1246434517, Z2283664368, Z1455079275, HY-N0343, and HY-B0069) and YTHDC1.

Data sources and processing

We obtained relevant clinical, somatic mutation, CNV, and RNA-seq data from the UCSC Xena database (https://xenabrowser.net/datapages/). Protein expression in EC tissues was analyzed using the UALCAN (http://ualcan.path.uab.edu/index.html) online tool. The GISTIC algorithm was used to examine CNV and assigned estimated values to -2, -1, 0, 1, and 2, representing homozygous deletion, single copy deletion, diploid normal copy, low-level copy number amplification, or high-level copy number amplification. To investigate the clinicopathological significance of CNV and mutations, we divided TCGA EC cohort data into two groups: with mutation and/or CNV and without CNV or mutation. The mRNA expression data were calculated from RNA-Seq V2 RSEM, applying log scale before analyzing the relationship between mRNA expression and CNV. Scores for six types of immune cells (B cell, CD4 T cell, CD8 T cell, Treg, neutrophil, macrophage, and natural killer (NK) cell) were obtained based on mRNA expression data using the TIMER package (https://cistrome.shinyapps.io/timer/).

Differential expression analyses and gene set enrichment analysis

We downloaded EC RNA-seq, miRNA-seq, and clinical data from TCGA database. According to the Ensembl (Homo sapiens; https://asia.ensembl.org/index.html) database, gene IDs were converted into gene names. miRNA-seq data were analyzed via the same method. Differentially expressed lncRNAs, mRNAs, and miRNAs were screened using the edgeR program package in R software. Differentially expressed transcripts had an absolute fold change (FC) >1.5 and P < 0.01. We used the R program to construct a volcano map of differentially expressed genes (DEGs). An absolute log2-FC > 1.5 and false discovery rate (FDR)-adjusted P value < 0.01 were used as the cutoff for significantly differentially expressed mRNAs (DEMs) and lncRNAs (DELs).

Construction of the ceRNA regulatory network

GO and KEGG enrichment analyses were carried out using the clusterProfiler R package. A Benjamini and Hochberg-adjusted P < 0.05 was considered statistically significant. DEGs were compared in public expression databases miRcode, miRTar Base, TargetScan, and miRDB. The lncRNA-miRNA-mRNA interaction was divided into respective lncRNA-miRNA and miRNA-mRNA interactions. By comparing lncRNA and miRNA with the miRcode database, predictions were made regarding lncRNA-miRNA interactions. We then cross-referenced the prediction results with the miRTarBase, TargetScan, and miRDB databases in order to obtain differentially expressed miRNA regulatory genes. Cytoscape v3.5.117 was used to visualize and map the network.

Patients and samples

We selected patients who had undergone surgical removal of the uterus at the Department of Gynecology, Shengjing Hospital, affiliated to China Medical University, from 2011 to 2017. We collected 34 samples of EC tissue and 20 samples of normal endometrial fresh tissue. In addition, 40 paraffin-embedded specimens, including 35 EC tissues and 20 normal endometrial tissues, were obtained. No patient had received anti-cancer treatments, including radiotherapy, chemotherapy, or immunotherapy, before surgery. The age range was 23–69 years, and histopathological analysis was performed by two pathologists. This study was approved by the Ethics Committee of Shengjing Hospital, affiliated to China Medical University (ethical approval number: 2018PS251K).

Cell culture

Ishikawa cells were cultured in RPMI 1640 medium (Gibco, Carlsbad, CA, USA), while HEC-1A cells were cultured in McCoy’s 5A medium (Gibco). Both cell lines were obtained from the Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences (Shanghai, China).

Transfection of Ishikawa and HEC-1A cells

Cells were transfected with siRNAs against ZC3H13 (si-ZC3H13) and YTHDC1 (si-YTHDC1) (GenePharma, Shanghai, China) using Lipofectamine™ 3000 (ThermoFischer, Carlsbad, CA, USA) as per the manufacturer instructions. si-NC (GenePharma, Shanghai, China) was used to transfect control cells. The si-RNA sequences are listed in Supplementary Table 1.

Cell proliferation assay

Ishikawa and HEC-1A cells at the logarithmic growth phase were uniformly suspended and transferred to 96-well plates. A blank control with three parallel wells was included in each group. An EdU cell proliferation detection kit (RiboBio, Guangzhou, China) was used to assess the cell proliferation capacity as per the manufacturer’s instructions. Briefly, 50 μM of EdU mixed reagent was prepared, and the cells were incubated with it for 2 h. After the cells were fixed, DNA staining was performed, and cells were washed with PBS. Image acquisition and analysis were performed using a fluorescence microscope (Nikon, Japan) at 20× the original magnification.

Cell invasion assay

After the Transwell chamber was sterilized, Matrigel (pore size 8 μm; Corning, NY, USA) was added and incubated at 37° C overnight for the gel to solidify. Cells in the logarithmic growth phase were harvested, centrifuged, and the supernatant was discarded. The cells were resuspended in serum-free medium. Subsequently, 800 μL of culture medium (containing 10% fetal bovine serum) was added to another 24-well plate, and the Transwell chamber was transferred into the 24-well plate. The cell suspension (5 × 104 cells) was added to the upper chamber and incubated for 24 h at 5% CO2 and 37° C. The Transwell chamber was then removed. The non-adherent cells in each well were removed, and the upper chamber liquid was discarded. After the well was air-dried, it was fixed in 4% paraformaldehyde. Crystal violet dye solution was added, and cells were rinsed with PBS. Images were acquired under a microscope, and cells were counted.

RNA extraction and quantitative RT-PCR

TRIzol reagent (Takara, Shiga, Japan) was used to extract total RNA from cultured cells and tissues. Reverse transcription of total RNA into cDNA was performed with PrimeScript™ RT-PCR Kit (Takara) and SYBR® TB Green™ Premix Ex Taq II (Takara) for real-time PCR analysis. Specific PCR primers were designed by Sangon Biotech Co., Ltd. (Shanghai, China). The fold change in expression was calculated using the 2^-∆∆Ct method, with GAPDH as an internal control. The primer sequences are listed in Supplementary Table 2.

Immunohistochemistry

Immunohistochemistry was carried out to detect the protein expression of ZC3H13, YTHDC1, and METTL14. The primary antibodies used are listed in Supplementary Table 3. Staining was primarily localized within the nucleus, seen as brown foci. The scoring method for positive expression was based on whether the cytoplasm had a brownish yellow or brown color. Samples were scored as follows: no staining, 0; light yellow, 1; yellow, 2; brown or sepia, 3. The scores were assigned according to the percentage of positive cells, with a negative count denoted as 0. The percentage of positive cells was scored as follows: < 10 % = 1, ≥ 10–50 % = 2, > 50–75 % = 3, and ≥ 75 % = 4. The product of the two scores was considered as the total score, and the results were interpreted as follows: ≤ 2 = negative, 3–4 = weak positive (+), 5–8 = medium positive (++), and 9–12 = strong positive (+++). Low expression was indicated by -/+, and high expression was indicated by ++/+++. The results were evaluated by two senior pathologists who were blinded to the patients’ data. Each sample was independently observed to determine the positive cell count and evaluate the background. In cases of disagreement, a third pathologist made the judgement.

Virtual screening process

The crystal structure of the target protein was obtained from the protein database (PDB database, https://www.rcsb.org/). Next, we used the Protein Preparation Wizard module to hydrogenate the protein, remove water molecules, SAM, ethylene glycol, and SEP, repair missing residues, and add side chains. The LigPrep Module was used to perform hydrogenation, energy optimization, and to construct the 3D structure. The Virtual Screening Workflow module was used for virtual screening. After importing the prepared compound, the Glide module was used for molecular docking, that is, the receptor and ligand molecules were aligned with each other through geometric and energy matching. The higher the absolute value of the docking score, the stronger the binding force.

Statistical analysis

We used GraphPad Prism 8 software (GraphPad, Inc., La Jolla, CA, USA) for statistical analysis. The difference in count data was assessed via the χ² test, and data are expressed as mean ± SEM. Wilcoxon’s signed-rank test was used to compare the immune cell infiltration between mutant and wild-type expression data. The difference between different CNV types were assessed using the Kruskal-Wallis test. Kaplan-Meier analysis and the log-rank test were used to compare the OS and DFS between patients with different CNV types. DFS was defined as the time from the date of diagnosis to the time of progression/death or last follow-up. All statistical tests were two-sided and implemented in the R programing language. P < 0.05 was considered statistically significant.

Ethics approval and consent to participate

The study protocol was reviewed and approved by the Scientific Research and New Technology Ethical Committee of the Shengjing Hospital of China Medical University. Ethical number: 2018PS251K.

Consent for publication

We have obtained consent to publish this paper from all the participants of this study.

Availability of data and materials

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