Identification of a Six-Gene Signature for Predicting the Overall Survival of Cervical Cancer Patients

Introduction

Cervical cancer (CC) is one of the main causes of cancer-related death in women and the third leading malignancy among women worldwide, after breast and colorectal cancers, and accounts for 569,000 new cancer cases each year.^1,2 CC is still a major problem for healthy women worldwide despite significant efforts. Although the development of CC screening has led to a drastic reduction in the incidence of CC, there are still many challenges in reducing the incidence of advanced lesions. Therefore, searching for potential biomarkers for CSCC treatment and prognosis prediction with multiomics data is necessary. However, few studies have been conducted to explore the relationship between genomic factors and CSCC prognosis.

At present, the biomarkers used for the prognosis of CSCC are mainly divided into two categories. One category consists of clinical biomarkers for CSCC survival prediction, such as squamous cell antigen (SCC). The other category consists of prognosis-related signatures constructed from several prognostic genes by high-throughput sequencing data analysis. For instance, Shen et al discovered that the CD28 and PTEN genes play important roles in the occurrence and development of CSCC through methylated microarray data analysis.³ Mao et al identified long noncoding RNAs (lncRNAs) related to the prognosis of CSCC from the TCGA database and developed a 15-lncRNA signature risk score to comprehensively assess the prognostic role of lncRNAs.⁴ Ma et al found a potential microRNA (miRNA) expression signature capable of predicting survival time for CSCC patients.⁵ Therefore, the identification of genes of prognostic value by bioinformatics analysis will likely be crucial and helpful in CSCC treatment and prognosis prediction. High-throughput multiomics sequencing data have laid a solid foundation for identifying genes related to cancer prognosis. Multiomics data analysis can reveal the mechanism of cancer development from multiple perspectives.

In this study, we analyzed transcriptome data to obtain genes related to the prognosis of CSCC. Then, we analyzed the copy number variation (CNV) data and mutation data of CSCC patients to obtain genes related to the occurrence and development of CSCC. Subsequently, we integrated the results and proposed a prognostic gene signature by using the random survival forest algorithm. We proved that the model performed effectively in the test dataset and the external independent dataset and that the signature shows strong clinical independence. Overall, the prognostic gene signature can effectively predict the prognosis of CSCC and provide a basis for a better understanding of the pathogenesis of CSCC.

Materials and Methods

Collection and Preprocessing of Publicly Available Data

The CSCC RNA sequencing (RNA-seq) fragments per kilobase of transcript per million mapped reads (FPKM) dataset and corresponding clinical follow-up information were downloaded from the publicly available TCGA database by the University of California, Santa Cruz (UCSC) Genome Browser (https://xenabrowser.net/datapages/); these datasets contain information on 309 samples and 306 samples, respectively. The CNV dataset of the SNP 6.0 array contains information for a total of 297 samples. The mutation annotation file (MAF) downloaded from the GDC client contains information for a total of 289 samples. The GSE44001 expression profile data and clinical follow-up information were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) and contains information for a total of 300 samples. To process the GSE44001 data, the software package GenomeStudio was applied to extract raw data. The intensity of the probes was log2-transformed and then normalized using the quantile normalization method.

The TCGA transcriptome data were standardized, and the samples were randomly divided into two groups: the training dataset (n=146) and the test dataset (n=147). GSE44001 was used as an independent validation dataset. The age, survival status, grade and tumor stage of patients who provided the samples are expressed in Table S1.

Univariate Cox Proportional Hazards Regression Analysis

Univariate Cox proportional hazards regression analysis was performed for each gene with transcriptome data in TCGA training dataset, and genes with a significance level of p < 0.05 were selected as prognostic genes.

CNV Analysis

GISTIC 2.0⁶ software was used to identify the genes with significant amplification or deletion. The threshold criteria were set as an amplification or a deletion length greater than 0.1 and z significance level of p < 0.05.

Gene Mutation Analysis

MUTSIG 2.0 software was used to identify the genes in MAF with significant mutations with p < 0.05.

Gene Signature Construction

First, the candidate genes significantly related to survival were screened for significant amplification, deletion, and mutation. Then, we used the random survival forest algorithm⁷ to rank the importance of the genes. The number of Monte Carlo iterations was 100, and the number of advances was 5 (nrep = 100, nstep = 5). Genes with a relative importance greater than 0.5 were identified as characteristic genes. Finally, we performed a multivariate Cox regression analysis to build the risk score model:

Risk score =

where N is the number of prognostic genes, is the expression value of prognostic genes, and is the estimated regression coefficient of genes in the multivariate Cox regression analysis.

Functional Enrichment Analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of genes were performed by using the R package clusterProfiler⁸ to identify overrepresented GO biological process terms and KEGG pathways. For this analysis, a false discovery rate (FDR) <0.05 was considered to denote statistical significance.

Immunohistochemical Staining (IHC)

We collected a total of 160 human cervical tissue samples, 130 of which had accompanying follow-up information, and 30 cancer-adjacent cervical tissue samples from paraffin-embedded tissue archives between January 2010 and January 2014 at the Department of Pathology of Peking Union Medical College Hospital. The follow-up was performed until December 2018. The pathological diagnoses were reconfirmed by a pathologist. The project was approved by the Ethical Committee (Peking Union Medical College Hospital), and informed consent was acquired from patients or family members (As some patients have passed away when signing the reconfirmed consent). IHC was performed as previously described.⁹ Antibodies against the following targets were used: CCL17 (1:200); DES (desmin; 1:100); DPAGT1 (1:200); FURIN (1:200); SLC19A3 (1:500); and SLC22A3 (1:500). The scoring details have been described previously.¹⁰ Positive staining was defined as follows: five random microscope fields were selected for semiquantification of the staining. The intensity of immunostaining was graded as follows: 1+, weak; 2+, moderate; 3+, strong or 4+, very strong. The area of positive cancer cells in each microscopic field was categorized as follows: 1+, 0 to 25%; 2+, 25 to 50%; 3+, 50 to 75% or 4+, 75 to 100%. A total score between 5 and 80 was obtained by multiplying the two scores by 5. A total score from 0 to 42 was assigned as “low expression” and that from 43 to 80 was assigned as “high expression”. Two pathologists independently reviewed the slides using their individual criteria, and consensus interpretations were used as the final interpretations. Discrepant interpretations were adjudicated by a third pathologist.

Results

Selection and Identification of the Prognostic Gene Set and the Mutated Gene Set

From the TCGA training dataset, a total of 1412 genes were selected as prognostic candidate genes based on univariate Cox regression analysis (p < 0.05). The top 20 most significant prognostic genes with the lowest p values are listed in Table 1 (more details in Table S2)

Table 1 Top 20 Most Significant Prognostic Genes

Moreover, genes with significant amplification or deletion were identified by GISTIC 2.0. Figure 1A shows the significantly amplified segments in the CC genome. The significantly amplified genes on each segment are listed in Table S3. For instance, ERBB2 was significantly amplified on the 17q12 segment (q value = 5.62E-12), MYC was significantly amplified on the 8q24.21 segment (q value = 7.70E-09), and EGFR was significantly amplified on the 7p11.2 segment (q value = 0.01188). In total, 624 genes were identified with significant amplification. Figure 1B shows the significantly deleted segments in the CC genome. The significantly deleted genes in each segment are recorded in Table S4. For instance, CD3D was significantly deleted in the 11q24.2 segment (q value = 9.11E-26), SMAD4 was significantly deleted in the 18q21.2 segment (q value = 0.0010189), and PTEN was significantly deleted in 10q23.31 (q value = 0.0077083). A total of 1038 genes were identified with significant deletions.

Figure 1 (A) Significantly amplified segments in the CC genome. (B) Significantly deleted segments in the CC genome. (C) The top 50 genes with the most significant p values. The upper histogram represents the total number of synonymous and nonsynonymous mutations in the 50 genes from each patient, and the right histogram represents the mutation proportion of 50 genes in all samples. (D) A total of 839 genes with significant copy number variation and mutations were included in the KEGG pathway analysis. (E) A total of 839 genes with significant CNV and mutations were included in the GO biological process analysis. (F) The relationship between prediction error rates and the number of separate trees. (G) The out-of-bag importance rank of the top 6 genes.

Subsequently, we used MUTSIG 2.0 to identify genes with significant mutations from the TCGA MAF data. P < 0.05 was set as the threshold. A total of 163 genes were identified with significant mutation frequencies. Figure 1C shows the top 50 genes with the highest mutation frequencies and the distribution of the different mutation forms, including synonymous mutations, missense mutations, frame insertions or deletions, frame movements, nonsense mutations, splice sites and other nonsynonymous mutations, in the training dataset. The upper histogram represents the total number of synonymous and nonsynonymous mutations among the 50 genes in each patient, and the right histogram represents the mutation proportion of the 50 genes in all samples. Some of these genes, such as FBXW7, PIK3CA, PTEN, KRAS, and TP53, were reported to be closely related to the occurrence and development of cancer in previous studies.

Furthermore, we analyzed the pathways and biological processes involved in genes with CNVs and mutations to investigate the potential mechanism. After the integration of genes with CNVs and genes with mutations, a total of 1643 genes were obtained for functional enrichment analysis. As shown in Figure 1D, the 1643 genes were significantly enriched in human T-cell leukemia virus 1 infection, apoptosis, the FoxO signaling pathway, the ErbB signaling pathway, central carbon metabolism in cancer, endometrial cancer, etc. As shown in Figure 1E, the 1643 genes were significantly associated with developmental processes, protein metabolic processes, and epidermal cell differentiation.

Prognostic Signature Construction for CSCC Patients

In this study, we aimed to obtain a prognostic signature for CSCC prognosis prediction. Therefore, we selected 75 genes with significant CNVs and mutations among the prognostic candidate gene set for stepwise importance ranking by the random survival forest algorithm. Finally, we identified 6 genes with relative importance greater than 0.5 as final candidates for the prognostic signature. The hazard risk (HR), z-score, p value, importance and relative importance of the 6 genes are shown in Table 2. Figure 1F shows the relationship between the prediction error rates and the number of separate trees, and Figure 1G expresses the out-of-bag importance rank of the top 6 genes. Multivariant Cox regression analysis was conducted to develop a 6-gene signature, and the risk model was constructed as follows:

Table 2 6-Genes Significantly Associated with the Overall Survival in the Training-Set Patients

Then, the 6-gene signature risk score for each sample was calculated. According to the median value of all the samples, patients were divided into a high-risk group and a low-risk group (cutoff=0.02513276). The scores from this 6-gene signature model in the TCGA training dataset are shown in Figure 2. Figure 2A shows that 73 patients each were divided into the high-risk group and the low-risk group. Figure 2B demonstrates the receiver operating characteristic (ROC) curve, and the area under the curve (AUC) at 60 months was 0.82. Figure 2C shows that as the risk score increased, the survival period decreased significantly, and the majority of patients in the high-risk group died within 60 months. We also analyzed the expression level of 6 genes in patients with high versus low-risk scores. Patients in the high-risk group tended to have higher expression levels of SLC19A3, FURIN, SLC22A3 and DPAGT1 and lower expression levels of CCL17 and DES. We also compared the calibration of the model in 1, 3, and 5 years. Figure 2D shows that the difference between the predicted value and the true value is small, which indicates the robust prediction performance of the model.

Figure 2 (A) Kaplan–Meier survival curve of patients stratified according to the 6-gene signature in the TCGA training dataset. (B) ROC curve and AUC of the 6-gene signature. (C) Risk score, survival period and gene expression level. (D) The nomogram prediction of OS.

We also analysis of the expression of 6 genes in normal tissues and tumor tissues as shown in Figure S1. None of the 6 genes showed significant differences. In the TCGA data set, there are only 3 normal samples and 304 primary tumor samples. Therefore, statistical tests may be biased.

Detection of the Robustness of the Gene Signature

To verify the accuracy of the 6-gene signature, the test dataset was assessed with the same model and cutoff value as training dataset. The Kaplan–Meier survival curve, ROC curve and AUC of the 6-gene signature and the relationships between the risk score, survival period and gene expression levels are shown in Figure 3. Figure 3A demonstrates the significant prognostic difference between the high-risk group and the low-risk group (log-rank p = 0.009019433). Figure 3B and C shows that SLC19A3, FURIN, SLC22A3 and DPAGT1 might be risk factors, while CCL17 and DES might be protective factors, which is consistent with the findings in the training dataset.

Figure 3 (A) Kaplan–Meier survival curve of patients stratified according to the 6-gene signature in the TCGA test dataset. (B) ROC curve and AUC of the 6-gene signature. (C) Risk score, survival period and gene expression level.

Furthermore, we validated the 6-gene signature in an external independent dataset downloaded from the GEO to confirm the stability with different platforms. The same formula and cutoff value were applied. Consistent with our previous findings, Figure 4 demonstrates the Kaplan–Meier curves of the GSE44001 dataset, which further suggested a significantly prolonged survival time in the low-risk patients compared to the high-risk patients (Figure 4A). The AUC at 60 months was 0.59, and the relation between gene expression levels and the risk score was also consistent with that in the TCGA dataset (Figure 4B and C). Therefore, the 6-gene signature has robust prognostic capabilities.

Figure 4 (A) Kaplan–Meier survival curve of patients stratified according to the 6-gene signature in independent validation dataset. (B) ROC curve and AUC of the 6-gene signature. (C) Risk score, survival period and gene expression level.

Independence Evaluation of the Prognostic Signature

To confirm the independence of the 6-gene signature in clinical application, univariate and multivariate Cox regression analyses were conducted with the TCGA training dataset, the TCGA test dataset and the GSE44001 dataset. Variates, including age, grade, pathologic T stage, pathologic N stage, pathologic M stage, tumor stage and 6-gene risk group, were analyzed systemically, and the results are shown in Table 3 (Table S5). Univariate Cox regression analysis verified that high-risk classification, pathologic T3/T4, and tumor stage III/IV were significantly related to survival in the TCGA training dataset; multivariate Cox regression analysis verified that high-risk classification (HR = 2.89, 95% CI = 1.20–6.94, p = 0.018), pathologic T3/T4, and tumor stage III/IV were independent prognostic factors. Moreover, high-risk classification, pathologic T4 stage and pathologic N1 stage were significantly related to survival in the TCGA test dataset, and high-risk classification (HR = 1.444, 95% CI = 1.05–1.97, p = 0.020), pathologic T4 stage and pathologic N1 stage were clinically independent prognostic markers. For the GSE44001 data, high-risk classification was significantly related to survival and was a clinically independent prognostic marker (HR = 1.40, 95% CI = 1.09–1.78, p = 0.007). We have included prognostic-related clinical features T and Stage, and built a nomogram such as Figure S2 with our RiskScore; it can be observed that nomogram, RiskScore and T staging benefit significantly Higher than the extreme curve, in which Nomogram is higher than RiskScore, T stage is similar to RiskScore, which indicates that nomogram has good reliability.

Table 3 Univariate and Multivariate Cox Regression Analysis to Verify the Independence of 6-Gene Signature in Test Datasets

To compare the 6-gene model with existing signatures, we selected several recently reported cervical cancer-related prognostic signatures, such as the 6-gene signature of Chen et al¹¹, the 4-gene signature of mei et al¹², and Pan et al¹³. We used the same method to calculate the risk score of each patient in the TCGA data set, and compared the RiskScore with the 1, 3, and 5-year ROC of these three models, such as Figure S3A–D. It can be observed that the RiskScore is 1, 3 and 5 years have the highest AUC. In addition, comparing the C-index of these four models shows that RiskScore has the highest C-index such as Figure S3E. DCA analysis shows that RiskScore has the highest benefit curve (Figure S3F).

Evaluation of the Prognostic Signature in CSCC by IHC

Table S6 summarizes the characteristics of all patients, including age, disease stage, and tumor grade. The expression of the six genes was evaluated in CSCC by IHC. The differences in CCL17, DES, DPAGT1, FURIN, SLC19A3 and SLC22A3 expression between CC tissues and adjacent normal tissues are shown in Figure 5. The expression of DPAGT1 (53.77 ± 1.84 vs 21.00 ± 3.49, p < 0.001), FURIN (48.55 ± 1.56 vs 28.83 ± 2.39, p < 0.001), SLC19A3 (52.45 ± 1.77 vs 29.67 ± 4.48, p < 0.001) and SLC22A3 (47.16 ± 1.44 vs 25.33 ± 3.32, p < 0.001) was significantly higher in CC tissues than in adjacent cancer tissues. Conversely, CCL17 (26.73 ± 1.65 vs 49.50 ± 3.79, p < 0.001) and DES (30.55 ± 1.50 vs 54.00 ± 2.90, p < 0.01) showed lower expression levels in CC tissues. The correlations between the expression of these genes and CC prognosis are shown in Figure 6. These data reveal that high expression of DPAGT1 (overall survival [OS], HR = 2.469, 95% CI 1.170–4.586, p = 0.016; progression-free survival [PFS], HR = 2.126, 95% CI 1.097–3.851, p = 0.025), FURIN (OS, HR = 2.291, 95% CI 1.158–4.589, p = 0.018; PFS, HR = 1.924, 95% CI 1.035–3.683, p = 0.040), SLC19A3 (OS, HR = 1.908, 95% CI 0.979–4.178, p = 0.058; PFS, HR = 1.993, 95% CI 1.093–4.184, p = 0.027) and SLC22A3 (OS, HR = 2.272, 95% CI 1.147–4.544, p = 0.019; PFS, HR = 2.113, 95% CI 1.135–4.030, p = 0.019) and low expression of CCL17 (OS, HR = 2.396, 95% CI 1.252–5.180, p = 0.010; PFS, HR = 2.102, 95% CI 1.132–4.157, p = 0.020) and DES (OS, HR = 2.179, 95% CI 1.079–4.231, p = 0.030; PFS, HR = 1.981, 95% CI 1.049–3.687, p = 0.036) is associated with a poor prognosis in patients with CC.

Figure 5 IHC for CCL17, DES, DPAGT1, FURIN, SLC19A3 and SLC22A3. CC-adjacent tissue (N=30) and CC tissue (N=130) samples. CC-adjacent tissue samples with weak immunostaining scores for DPAGT1 (G), FURIN (J), SLC19A3 (M), and SLC22A3 (P) and strong immunostaining scores for CCL17 (A) and DES (D). CC samples with weak and strong immunostaining scores for CCL17 (B and C), DES (E and F), DPAGT1 (H and I), FURIN (K and L), SLC19A3 (N and O) and SLC22A3 (Q and R). The expression of the CCL17, DES, DPAGT1, FURIN, SLC19A3 and SLC22A3 genes is depicted in S slides (X 100). **p < 0.01, ***p < 0.001.

Figure 6 OS and PFS curves for CC (N=130) according to CCL17 (A and B), DES (C and D), DPAGT1 (E and F), FURIN (G and H), SLC19A3 (I and J) and SLC22A3 (K and L) gene expression status (low or high). The expression status for each sample was classified according to the median gene expression value.

Discussion

There were estimated 569,800 new CC cases and 311,400 CC-related deaths worldwide in 2018.¹⁴ Nearly 90% of CC-related deaths occur in developing parts of the world. In several Western countries, where screening programs have long been established, CC incidence rates have decreased by as much as 65% over the past 40 years. However, in contrast to the favorable overall trends, the CC incidence rate is reported to be rising in Uganda and in some countries of Eastern Europe (Estonia, Lithuania, and Bulgaria).¹⁵ Although increasingly preventable through the use of vaccines and curable through early cytological detection, women with the advanced or recurrent disease face a dismal prognosis with a potentially considerable risk of morbidity and mortality.¹⁶ CC treatment and prognostic prediction remain a major challenge. Heterogeneity is evident in the data regarding OS, PFS and local recurrence.¹⁷ If tumor behavior can be reliably predicted at the initial diagnosis, the prognosis of CSCC would probably be greatly improved. Therefore, it is critical to investigate the molecular mechanism of CSCC and to identify novel biomarkers.

Some gene signatures, such as multigene signature panels (MSPs), are currently being used in the clinic to evaluate breast cancer recurrence,¹⁸ which indicates that screening gene prognostic biomarkers by gene expression profiling is a reliable molecular identification method. Lin et al constructed a prognostic model for CC patients undergoing radiotherapy and chemotherapy with stage IB-IV disease based on magnetic resonance imaging (MRI), the whole-tumor apparent diffusion coefficient (ADC), and HPV genotyping.¹⁹ Mao et al developed a 26-IncRNA signature to predict the OS of CSCC patients.²⁰ In our study, we downloaded RNA-seq data, CNV data and MAF data for the CC samples from the TCGA database and divided them into a training dataset (N=146) and a test dataset (N=147). By analyzing the CNV data and mutation data of the training set, potential genes related to the occurrence and development of CC were obtained, and the transcriptome data were used to screen for biomarkers related to the prognosis of CC. Subsequently, the two sets of results were integrated. The prognostic signature model was determined by applying the random survival forest algorithm to the overlapping results. A 6-gene model was constructed, and the test dataset and external dataset (GSE44001) further verified that it was robust. The 6-gene signature has a higher AUC and fewer genes than the existing signature, which might be conducive to its clinical application.

In our 6-gene signature, DPAGT1, FURIN, SLC19A3 and SLC22A3 were verified to be risk factors, while CCL17 and DES were verified to be protective factors. Some previous studies have revealed the association between the 6 genes we identified and the pathogenesis of some other cancers. It was reported that CCL17 plays a critical role in preventing tumorigenesis due to its ability to engage CCR4+CD8+ T cells.^21–23 Okada et al transfected an RGD fiber-mutant adenoviral vector into the tumors of melanoma-bearing mice and found that recombinant AdRGD-CCL17 induced CD8+ cytotoxic T lymphocyte (CTL) aggregation in the tumor tissue, which effectively inhibited tumor growth.²⁴ In addition, the Naoko team also used tumor-bearing animal experiments to confirm that injecting recombinant CCL17 adenovirus into the tumors of mice with colorectal cancer induces the infiltration of activated CD8+ CTLs and exerts significant antitumor immune function.²⁵ Zhan et al also inoculated CCL17-loaded drug-eluting stents into pancreatic cancer-bearing mice, and locally released CCL17 recruited CCR4+CD8+ T cells to pancreatic tumor tissues, thereby inhibiting the occurrence and metastasis of pancreatic tumors.²⁶ Interestingly, many studies have shown that CCL17 can also stimulate cancer cell proliferation and migration.^27,28 However, in the present study, the expression of CCL17 was evaluated in CC tissue, and the results showed that lower expression of CCL17 was associated with a poor prognosis in CC patients.^29,30 DES encodes a muscle-specific class III intermediate filament that plays a crucial role in maintaining the structure of sarcomeres, interconnecting Z-disks and forming myofibrils.³¹ A study revealed that the expression of DES was significantly downregulated in gallbladder carcinogenesis according to its correlation with promoter hypermethylation, which indicates its potential as a candidate biomarker for gallbladder carcinogenesis.³² DES was also demonstrated to inhibit telomerase activity in prostate cancer cells.³³ Moreover, previous studies also reported that the DES is a prognostic predictor and therapeutic target for colorectal cancer.^34,35 However, there are few studies on DES in CC. In the present study, we demonstrated that DES was significantly downregulated in CC, which was associated with a poor prognosis. DPAGT1 is an upstream regulator of E-cadherin N-glycosylation status and adherens junction composition, and dysregulation of DPAGT1 may cause disturbances in intercellular adhesion in oral cancer.^36,37 In 2019, data demonstrated that overexpression of DPAGT1 was significantly associated with a poor chemotherapeutic response in hepatitis B virus-positive hepatocellular carcinoma patients.³⁸ FURIN encodes a member of the subtilisin-like proprotein convertase family, which includes proteases that process protein and peptide precursor trafficking through regulated or constitutive branches of the secretory pathway. Studies have revealed that FURIN is associated with the prognosis of pancreatic cancer,³⁹ laryngeal cancer,⁴⁰ and lung cancer.⁴¹ More importantly, Zhao et al recently proposed a five-mRNA signature to predict recurrence and survival in CC that also included DPAGT1 and FURIN.⁴² However, the difference is that they only analyzed gene expression in CC, and we performed an integrated analysis of mRNA, CNV and gene mutation data to identify changes that were significantly related to HPV integration. In the present study, we also evaluated the expression of signature genes in more CC tissue by IHC. SLC19A3 is an ubiquitously expressed transmembrane thiamine transporter. It was demonstrated that downregulation of SLC19A3 can inhibit breast cancer metastasis.⁴³ Another study revealed that aberrant SLC19A3 promoter hypermethylation in plasma may be a novel biomarker for breast and gastric cancer diagnosis.⁴⁴ SLC22A3 is an integral plasma membrane protein located in a cluster on chromosome 6. Xin-Guan et al reported that the expression of SLC22A3 was associated with the early development and progression of familial esophageal cancer.⁴⁵ The rs2504938 SNP in SLC22A3 significantly influenced the OS of pancreatic cancer patients. Moreover, SLC22A3 was found to play an important role in colorectal cancer development and progression.⁴⁶

We integrated multiomics data, including transcriptome, CNV data and mutation data, and eventually constructed a 6-gene signature for CC prognosis. The signature was verified to be robust and stable in different datasets. Furthermore, the signature was shown to be optimally independent, even considering the influence of many clinical factors. However, the current study also has several limitations. First, the present study is based on bioinformatics analysis and verified the prognostic of single gene; however, the signature cannot be verified by experimental verification at present; secondly, the value of the 6-gene signature for therapeutic strategies has not been shown in patients’ cohorts, such as immunotherapeutic and chemotherapeutic prediction; Finally, the signature was constructed by the FPKM data, which is widely used in bioinformatics analysis; therefore, further study should strive to transform the signature into clinically applicable models. Therefore, further experimental verification studies are necessary.

Conclusion

In the present study, we constructed a 6-gene signature for CC prognosis based on multiomics data. The predictive value of the signature was further tested in a training dataset and an external independent validation dataset by applying Kaplan–Meier survival curves. The signature shows strong robustness and clinical independence. Therefore, the 6-gene signature can be used as a prognostic assessment tool for CC patients.