DOI: https://doi.org/10.21203/rs.3.rs-1973956/v1
Background: The prognosis for advanced and recurrent metastatic cervical cancer is poor. Angiogenesis plays a vital role in tumor development and the tumor microenvironment (TME).
Methods: We performed a consensus clustering analysis of RNA-seq data based on ARG for CESC patients downloaded from TCGA. Then we analyzed the characteristics, prognosis, and immune infiltration status among the subtypes. Then we constructed predictive models and ARGscore. And we explored the relationship between ARG scores and prognosis, TME, and immunotherapy correlation.
Results: We found that most ARG expression was upregulated in CESC compared to normal samples and clarified the mutation of ARG in CESC. We divided the 290 CESC patients into 2 ARG clusters by consensus clustering. We observed significant differences in their survival and immune infiltration status. Subsequently, ARGscore that can predict prognosis was established. We found that the high-risk group predicted a poorer prognosis. We have verified that ARG scores have good accuracy. At the same time, we confirmed that ARG scores were closely related to TME. A reliable nomogram was developed to facilitate the clinical competence of ARG scores. In addition, we explored the relationship between ARG scores and TMB and found no correlation. However, the prognosis of the high-TMB group was better than that of the low-TMB group, and the ARGscore could offset the advantage. The TIDE score validated the possibility that ARG scores predict immunotherapy.
Conclusion: With this study, we obtained an ARG score based on the ARG established to assess the TME status and prognostic risk of patients and provide a basis for immunotherapy.
Cervical cancer (CC) is one of the most common cancers in women worldwide, and it occurs mainly in low- and middle-income countries [1]. The World Health Organization (WHO) launched a global strategy project that includes human papillomavirus (HPV) vaccination and cervical screening to accelerate the elimination of CC [2]. With the widespread availability of CC screening and vaccines, initial progress has been made in actions to prevent CC. However, the incidence of CC is still on the rise in some countries, and the prognosis for advanced and recurrent metastatic CC is unfavorable, with a 5-year recurrence rate higher than 20%[3]. Studies have shown that even in resource-rich areas, about 30% of CC patients still recur or even die [1]. The late-stage stages and recurrence of CC are the leading cause of death in patients with CC. Treatment of patients with advanced and recurrent CC can be tricky due to prior radiation or chemotherapy history.
In the TME, tumor cells form an intricate network of abnormal blood vessels by secreting large amounts of pro-angiogenic factors (e.g., VEGF, PDFG-B, TGF-β) that promote tumor appreciation and metastasis. The formation of the vascular network also hinders the pernicious effect of immune cells of TME on tumor cells[4], and the activity of immune cells in TME determines the effectiveness of immunotherapy. So there is a close correlation between immune response and angiogenesis [5].
Immunotherapy has been an emerging therapeutic tool in recent years for treating many malignancies. Immune checkpoint inhibitors (ICIs) enhance tumor killing by targeting PD-1, PD-L1, and CTLA-4 receptors, thereby activating T-cell activity [6]. Advanced CC can benefit from immunotherapy with promising antitumor activity [7]. Commonly used predictive markers (e.g., PD-L1) have low sensitivity when evaluating immunotherapy efficacy in CC yet [8]. There are no specific predictive biomarkers for the clinical efficacy of immunotherapy in CC.
Recent studies have shown the efficacy of immunotherapy combined with targeted therapy in advanced and recurrent CC [9, 10]. Therefore, targeted therapy and immunotherapy are the future treatment directions for CC. The development of various targeted drugs, especially anti-angiogenic therapy, has a promising future in treating CC. However, the current systemic treatment plan for recurrent metastatic CC is still incomplete.
Multiple lines of evidence suggest a close relationship between angiogenesis, tumorigenesis, and progression. Nevertheless, there are few basic studies on the relationship between angiogenesis and CC. A few clinical trial studies on anti-angiogenesis for advanced CC have shown promising benefits of anti-angiogenesis for both OS and PFS in CC [11–13]. Given the critical role of ARG in tumors, this study uses ARG for risk scoring and grouping based on gene expression profiles from TCGA. It develops a prognosis-related ARG prediction signature with a nomogram. At the same time, CC is a highly heterogeneous disease. The staging by TNM and FIGO alone can no longer meet the needs of its risk assessment. Thus, it is crucial to explore biomarker that is important for the treatment and prognostic assessment of CC. To establish a multidimensional signature to identify high-risk patients for personalized treatment of CC patients.
RNA sequencing, somatic mutation, and corresponding clinical data were downloaded for Cervical endocervical adenocarcinoma and squamous cell carcinoma (CESC) from the UCSC Xena database (http://xena.ucsc.edu/) in TCGA TARGET GTEx. Cases with no survival data were excluded. A total of 290 cases of tumor tissue and 10 cases of normal tissue were included. Gene expression and clinical data of GSE44001 were downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/) as a validation set. 181 angiogenesis-related genes were obtained from the GeneCard database (https://www.genecards.org/) (Table S1).
To observe the difference in ARG expression between tumor and normal tissues in CESC, we obtained differentially expressed genes (DEG) by R package "DESeq2" (|log2FoldChange|>2, P < 0.05). Then, the interactions between ARG gene proteins were viewed by the STRING database (https://cn.string-db.org/), and hub genes were found by Cytoscape. In addition, we used the "maftools" package to check the mutation status of ARG in CESC patients.
Consistency clustering can distinguish CESC into several subtypes based on ARG so that new ARG subtypes can be discovered and a comparative analysis of ARG subtypes can be performed. The "Consensus ClusteringPlus" package implements this process. To ensure the stability of the classification, we repeated this analysis 1000 times. Differentially expressed genes (DEGs) were identified by differential analysis for gene set variation analysis (GSEA) to determine the biological functional differences between ARGs. GSEA can make a straightforward enrichment analysis of the functional gene dataset composed of genes with the complete data obtained from sequencing and microarrays. We applied the MSIgDB dataset for the analysis.
We investigated the relationship between clinical characteristics and survival analysis to determine the clinical significance of the consistent clustering results. Clinical parameters include stage, grade, T, N, and M. Kalpan-Meier analysis was performed with the "survival" and "survminer" packages to assess the differences in OS between subgroups.
ESTIMATE is a tool for predicting tumor purity and for predicting infiltrating stromal/immune cells in tumor tissue using gene expression data. Based on this algorithm, we can evaluate the substrate score, immune score, and ESTIMATE score of the ARG cluster. CIBERSORT is a tool for deconvolution of the expression matrix of human immune cell subtypes based on the principle of linear support vector regression, allowing the calculation of the levels of 22 immune cell subtypes per patient to understand TME. The infiltrated fraction of immune cells was identified by single-sample gene enrichment analysis (ssGSEA). Moreover, we analyzed the expression of PD-1, PD-L1, and CTLA-4 between clustered subgroups.
To determine the DEG of the ARG cluster, we used the "DESeqs" package (|log2FoldChange| > 2, p < 0.05) for the analysis of variance. GO, and KEGG are databases of gene-related functions stored based on different taxonomic ideas. The GO database classifies target genes by cellular component (CC), molecular function (MF), and bioengineering (BP). The KEGG database is one of the pathway-related databases. GO and KEGG enrichment analyses were performed using the "clusterProfiler" package based on these DEGs. We clustered patients by DEG using the "Consensus ClusteringPlus" package to explore the biological activity between ARG clusters further. Then we performed a survival analysis of the DEG cluster. Meanwhile, we observed the expression of ARG between DEG clustering subgroups.
The ARGscore was used to assess angiogenesis in patients with CESC quantitatively. We performed a univariable Cox analysis on DEGs of the ARG cluster and selected genes with survival significance for the next step of the analysis. We build predictive models by Least absolute shrinkage and selection operator (LASSO) Cox regression to obtain the genes and correlation coefficients for building the models. ARGscore = Expression of gene [1] * corresponding correlation coefficient [1] + Expression of gene [2] * corresponding correlation coefficient [2] + Expression of gene [n] * corresponding correlation coefficient [n].
We calculated the respective risk scores for the ARG cluster and DEG cluster. Then we compared them, as well as compared the survival analysis of the different risk groups. The accuracy of ARGscore in predicting patient survival at 1, 3, and 5 years was assessed by timeROC.
We investigated the correlation between ARGscore and clinical parameters. To determine whether the ARGscore was an independent prognostic factor, we performed univariate and multifactorial Cox analyses for all clinical parameters. We integrated clinical parameters based on COX regression analysis. We scored according to the contribution of each clinical parameter to 1-year, 3-year, and 5-year OS and calculated Patients' OS at 1, 3, and 5 years. A nomogram was constructed for risk prognosis.
We compared the correlation between the levels of immune cells and ARGscore between different risk groups by the results of CIBERSORT. Also, we compared the correlation between ARGscore-related genes and immune cells. The expression of immune checkpoints in high and low-risk groups was also observed.
We investigated the relationship between ARGscore and TMB by calculating the TMB of patients from somatic mutation data. In addition, we used the TIDE web tool ( http://tide.dfci.harvard.edu/) to predict the potential response to immunotherapy in high- and low-risk groups.
R software (version 4.1.3) and its corresponding software packages were used to process, analyze, and present the data. P < 0.05 was considered to be significant.
We determined the expression levels of 181 ARGs in TCGA TARGET GETx-CESC in tumor and normal samples. DEG was selected between tumor tissue and normal tissue. Then the DEG and ARG were taken to intersect to obtain the DEG associated with angiogenesis between the two samples. As shown in Fig. 1A, most DEGs were abundantly expressed in tumor tissues. To reveal DEG interactions, we established a protein-protein interaction (PPI) analysis through the STRING website (Fig. 1B). Further, 10 hub genes were found by Cytoscape. They are STS3, CASP3, MAPK1, CREB1, MAPK14, MAPK14, FGF2, NFKB1, GLI1, MKI67,and MKI67.
We next further explored the 181 ARG somatic mutations in CESC. Somatic mutations were found in 137 of 287 samples (44.7%). MKI67 has the highest mutation rate. It was followed by TP53 and FBN1 (Fig. 1C)
To explore the relationship between angiogenesis and CESC, we selected 290 tumor patients with TCGA-CESC. The prognostic value of 181 ARGs was shown by univariate COX regression (Table S2).
To further clarify the relationship between ARG expression and CESC, we classified CESC into Cluster 1 (n = 194) and Cluster 2 (n = 96) by consensus clustering analysis based on ARG expression levels. PCA analysis showed a significant difference in gene expression between ARG clusters (Fig. 2A). Furthermore, by comparing the OS time of ARG clusters, we found a significant difference in the survival of the ARG cluster (Fig. 2B). After analyzing the clinical information of ARG clusters, we found no statistical difference between the two clusters in TNM stage, grade, and stage (Table S3).
According to the results of the GSEA analysis, the DEG between ARG clusters was mainly enriched in cell movement-related pathways and therefore associated with tumor metastasis (Table S4). To determine the relationship between TME in ARG and CESC, we investigated the level of infiltration of immune cell subpopulations in both groups by CIBERSORT and ssGSEA (Fig. 3A-B). Substantial differences in the degree of enrichment of most immune cells between ARG clusters. The enrichment levels of CD8 T cells, CD4 memory T cells, follicular helper T cells, gamma delta phenotype T cells, NK cells, M1-type macrophages, dendritic cells, Th1 and Th2, were significantly higher in Cluster 1 than in Cluster 2. CD8 central type memory T cells, mast cells, eosinophils, natural killer cells, and Th17 showed the opposite. It suggests that ARG clusters have different immune-related. Meanwhile, we explored the expression levels of three vital immune checkpoints (ICPs) between ARG clusters. These three ICPs are the drug inhibitor targets of choice for current clinical trials. We found that PD-1 was highly expressed in Cluster 1, while PD-L1 was highly expressed in cluster 2 (Fig. 3C).
TME scores can assess the abundance of immune and matrix elements in TME. We further implemented ESTIMATE's algorithm to evaluate TME scores between different clusters, including matrix scores, immune scores, and ESTIMATE scores. The study results showed that the TME score of cluster 2 was significantly higher than that of cluster 1 (Fig. 3D). The results of the immune evaluation showed the existence of different immune statuses in ARG clusters. CD8 T cell expression was significantly different. ICPs activate the immune response by regulating T cells in the immune response. We speculate that in CESC, different clusters may respond differently to using the ICIs process.
To investigate the potential biological activity between ARG clusters, we obtained 476 DEGs between the two clusters. Functional enrichment analysis was performed on DEG (Fig. 4A). These DEGs are mainly focused on biological processes such as angiogenesis and cell motility. KEGG pathway analysis yields associations with angiogenesis, EGFR resistance pathways, cell motility, etc. This result is similar to the GO enrichment results (Fig. 4B and Table S5). We then conduct a univariate COX analysis. The presence of 52 genes was determined to be survival significance (P < 0.05). To investigate specific regulatory mechanisms, we classified CESC patients into 2 DEG clusters according to DEG by consensus clustering (Fig. 4D). Kaplan-Meier analysis showed a shorter survival time for DEG cluster A (P = 0.033) (Fig. 4C). Figure 4E shows that most ARGs show substantial differences in expression between the two DEG clusters. Thus the DEG cluster verified the differences between the ARG clusters.
ARGscore are built from the DEG of the ARG cluster. We performed LASSO Cox on 52 DEGs to establish the best prediction model. Finally, we obtained 12 genes (Figure S1). The ARGscore were as follows: Risk score=(-0.0487*CFAP73) + (0.0041*TSPAN8) + (0.1323*TLL1) + (0.1427*IL1B) + (-0.3913*C11ORF16) +(0.0890*BCO0) +(0.3722*ATOH1) +(0.6696*PCDHAC2) +(1.4462*CCDC175) +(0.0198*PAPPA) +(0.5123*SYT16) +(0.0238*SLC35F3).
By comparing the respective risk scores of ARG and DEG clusters, we found that cluster 2 had a higher score than the ARG cluster (Fig. 5A). The fraction of Cluster B in the DEG cluster is higher (Fig. 5B). We considered the association of a high ARGscore with reduced survival in conjunction with survival analysis. In addition, the Kaplan-Meier analysis showed better survival in low-risk patients (Fig. 5C). The AUCs at 1, 3, and 5 years were 0.815,0.753 and 0.744, respectively (Fig. 5D). Figure 5E-F also shows that as the ARGscore increases, the time to OS decreases, and mortality increases.
To assess the predictive robustness of the ARGscore, we calculated the risk scores for the validation set. The ARGscore was divided into high and low-risk groups (Figure S2). Kaplan-Meier analysis reveals that survival is lower in the high-risk group compared to the low-risk group. Predictions of 1-, 3-, and 5-year survival demonstrated that the ARGscore still had a respectable AUC score. It means that the ARGscore has a good performance in assessing the prognosis of CESC patients.
To determine the relationship between ARGscore and clinicopathological features, we discussed the relationship between ARGscore and stage, grade, T, N, M, and survival status. We found that the higher the T stage, the higher the risk score (Fig. 6A-J). Furthermore, we performed univariate Cox and multivariate Cox analyses to explore the prognostic independence of multiple clinical factors. The results showed significant differences in risk scores (Fig. 6K-J), which were independent prognostic factors.
Since risk scores are highly correlated with patient prognosis, we built a nomogram combining clinical parameters (Fig. 6M). This nomogram was used to assess patients' OS at 1, 3, and 5 years. The calibration curve of the nomogram shows the essential accuracy between the actual observed and predicted values (Figure S3). In addition, we found that this predictive model with multiple clinical factors had a more significant net benefit in predicting prognosis.
We explored the correlation between ARGscore and immune cells based on the results of the CIBERSORT analysis. Figure 7A shows that ARGscore was positively correlated with NK cells, mast cells, B cells, and eosinophils. It is inversely related to helper follicular T cells, CD8 T cells, CD4 T cells, M1 macrophages, dendritic cells, and neutrophils. We then investigated the correlation between the 12 genes that established the predicted ARGscore and immune cells. Most immune cells were found to be closely associated with the selected genes (Fig. 7B). In addition, we selected 22 genes associated with ICPs and assessed their expression between high and low-risk groups. The study showed that the expression of most ICPs was significantly different between the two groups (Fig. 7C). It indicated that the ARGscore has a guiding role in immunotherapy.
Numerous studies have demonstrated that TMB is a valuable predictor of tumor immune response. High TMB can benefit from immune checkpoint inhibitors [14]. We compared the expression of TMB in high and low-risk groups. The results showed no significant differences and correlation between ARGscore and TMB (Fig. 8A-B). However, we divided the patients into a high TMB group and a low TMB group (TMB = 2.64 as the cut point). The high TMB group was discovered to have a better prognosis than the low TMB group (Fig. 8C). We then combined TMB and ARGscores for survival analysis of CESC patients. The ARGscore eliminated the prognostic benefit in the high TMB group (Fig. 8D).
Immunotherapy is of great clinical value in the treatment of tumors. Nevertheless, only some patients are responders to immunotherapy. TIDE algorithm is a Method using induction of T cell dysfunction in tumors with high infiltrating cytotoxic T cells (CTL) and prevention of T cell infiltration in tumors with low infiltrating cytotoxic T cells. It can mimic the tumor immune escape mechanism to predict the response to tumor immunotherapy[15]. We calculated TIDE scores to assess the immune response of CESC patients. The findings displayed a higher TIDE score in the low-risk group than in the high-risk group (Fig. 8E). It was indicated that the high-risk group responds less to immunotherapy than the low-risk group. The low-risk group is more likely to benefit from immunotherapy.
In TME, a variety of stromal cells, such as fibroblasts and immune cells, are also present in addition to tumor cells. Together they constructed the TME [16]. Among them, immune cells play an essential dual anti-tumor and pro-tumor role in TME. Thus the immune status presented by TME also provides food for thought for immunotherapy of tumors. Angiogenesis is one of the characteristics of TME. Haphazard tumor vasculature inhibits T cell entry into TME and suppresses its activity [5]. Angiogenic factors (e.g., VEGF) can also induce CD8 + T cell failure [17] and loss of movement. Th2 released from M2-type macrophages and Tregs promote aberrant and abnormal functional angiogenesis [18]. Cytokines released from inflammatory immune cells in TME also promote tumor angiogenesis, leading to tumor cell appreciation and metastasis [19]. There is a close relationship between CC and tumor neoangiogenesis. Some studies have demonstrated that HPV mediates tumor angiogenesis [20]. This study designed the CESC predictive model based on ARG and explored its immune function and TME status to be clinically instructive.
We observed the expression and mutation of ARGs in the TCGA-CESC cohort. Most ARG-associated DEGs were upregulated in tumor tissues and were associated with prognosis. It implies that increased ARG expression may be related to tumorigenesis and proliferation. According to ARG, we innovatively propose to cluster the CESC into cluster 1 and cluster 2 according to ARG. We identified two clusters that showed significant differences in survival and immune status of patients. Cluster 1 is enriched in CD8 + T cells, natural killer cells, dendritic cells, and M1-type macrophages. It was shown that granzyme and perforin secreted by natural killer cells and Th1 secreted by M1-type macrophages promote tumor cell killing by immune cells in TME, thereby inhibiting tumor appreciation [21, 22]. In addition, PD-1 expression was higher in Cluster 1, while PD-L1 expression was lower than in cluster 2. PD-1 receptor expression is highly regulated in activated T cells. PD-L1 receptors are mainly expressed on antigen-presenting cells (APCs) and tumor cells [23]. Upregulation of PD-L1 expression in TME by tumor cells and APC is a tumor immune escape strategy, thus reducing T cell activity. It is suggested that there may be more activated T cells in cluster 1, which are more capable of killing tumor cells. And cluster 2 may be more sensitive to anti-PD-L1 immunotherapy.
Subsequently, by their DEGs, we clustered the ARG clusters again into two DEG clusters (cluster A and B). The results showed that the DEG clusters differed significantly in survival and angiogenesis. We established an ARGscore based on the DEG of the ARG cluster for predicting the prognosis of patients with CESC. Six genes that set the ARGscore are involved in tumor development. TSPAN8 is a prognostic marker. It is highly expressed in various cancer tissues and mediates cancer cell proliferation and metastasis [24]. In particular, TSPAN8 is a potential target for immunotherapy of certain tumors [25]. IL1B is a pro-inflammatory cytokine that has been associated with tumorigenesis. It was demonstrated that IL1B is overexpressed in CC and predicts a poor prognosis [26]. PAPPA is a protease that plays a vital role in pregnancy. But it is an oncogene in areas other than pregnancy [27]. In TME of some tumors, increased expression of PAPPA promotes tumor appreciation and inhibits the pernicious effect of CD8 + T cells [28]. ATOH1 is essential for differentiating many cells, and various studies have demonstrated that ATOH1 deficiency enhances tumor formation and progression, especially in colorectal cancer [29, 30]. SLC35F3 was shown to be used as a new marker for the detection of colorectal cancer in a recent study [31]. SYT16 is less frequently reported in tumors. It has been reported that SYT16 is a prognostic marker for glioma and is closely associated with immune infiltration [32]. We further evaluated the respective risk scores of the ARG cluster and the DEG cluster. Cluster 2 (ARG cluster) and cluster A (DEG cluster) were found to have a higher risk score for a worse prognosis. Patients in the high-risk group had lower survival rates than those in the low-risk group. It was implied that high ARGscore were associated with poor prognosis, which would make the ARGscore a predictive marker of prognosis.
To further validate the relationship between ARGscore and prognosis, we made a observation about the relationship between ARGscore and clinical parameters and survival. The ARGscore was found to be an independent prognostic biomarker that can be used to assess the survival outcome of CESC patients. Based on this, we created a clinical prediction model. The ARGscore has been validated to have good accuracy and can accurately predict prognosis. Of particular interest is the integrated prospect that our study reveals the interaction between the angiogenic features of CESC and the immune phenotype. Angiogenesis and immune escape are complementary. They are usually present at the same time. There is growing evidence that angiogenic factors inhibit the ability of immune cells to fight tumors through a range of mechanisms [33]. Combined with previous evidence, our study adds to the evidence for a complex interaction between angiogenesis in CESC and immune function.
The development and application of ICI bring hope to many patients with advanced tumors. But the clinical application also faces the problems of drug resistance and limited audience. Anti-vascular targeting combined with ICI therapy not only inhibits blood vessels necessary for tumor growth and metastasis but also recodes the tumor immune microenvironment (TIME). The efficacy of anti-vascular targeting combined with ICI therapy has been validated in several clinical trials [34]. In this study, we investigated the differences in the expression of immune checkpoint genes between high and low-risk groups. It indicates that CESC patients in different risk groups have other immunotherapy effects. The TIDE algorithm proves this point as well. It implies that the ARGscore is instructive for immunotherapy.
Currently, the main clinical predictor of ICI treatment effectiveness is the expression level of PD-L1 measured by immunohistochemistry (IHC). However, PD-L1 expression is less predictive in advanced CC. In addition to PD-L1, microsatellite instability (MSI) and defective mismatch repair (dMMR) were shown to be equally predictive of ICI. High microsatellite instability (MSI-H) is more sensitive to ICI treatment Studies have shown that MSI-H status is deficient in patients with CC [35, 36]. TMB has recently been considered a new marker for predicting clinical response to ICI [37]. There was no difference in TMB expression in this study's high and low-risk groups. It also did not correlate with risk scores. However, grouping TMB at a cut point of 2.64mut/Mb showed that the high TMB group had lower OS than the low TMB group. TMB combined with ARGscore also made a difference to OS. Huang et al. [38] displayed that TMB showed sensitivity to immunotherapy at a cut point of 5mut/Mb and correlated with OS. This study could not evaluate the immune efficacy of the high and low TMB groups due to data limitations. Nevertheless, the analysis of the results of TMB and ARGscore on OS suggests that TMB can predict the response and prognosis of immunotherapy by combining ARGscore. It is recommended that TMB can expect the answer and prognosis of immunotherapy by combining ARGscore. The development of new checkpoint immunotherapy response biomarkers is essential to improve treatment efficacy. The ARGscore can be used as a predictive marker for the effect of CESC immunotherapy in this study.
Existing studies of scoring systems and prognostic scores for CESC have focused on immunohistology and genetic alterations. In contrast, the ARGscore developed in our research is a promising breakthrough in regulating TME by CESC in terms of angiogenesis. ARGscore provides new insights into the immune diversity of CESC from an angiogenic perspective. Meanwhile, it emphasizes that ARGscore can predict the sensitivity of immunotherapy. Considering ARGscore when selecting comprehensive anticancer therapy may improve patient outcomes. However, to maximize the efficacy of immunotherapy, we should incorporate more clinical and tumor microenvironmental factors.
In conclusion, the current research demonstrates that ARG influences the TME of CC. Further we confirmed that the ARG-based risk score model could predicts prognosis and sensitivity to immunotherapy. The study provides a basis for future research on the use of immunotherapy and clinical trials.
Data Avilability Statement
The datasets analyzed during the current study are available in the UCSC Xena database (http://xena.ucsc.edu/)，GEO (https://www.ncbi.nlm.nih.gov/geo/)，and GeneCard database (https://www.genecards.org/).
Acknowledgments
Thanks to XIANTAO Academic website (https://www.xiantao.love/)and ChiPlot (https://www.chiplot.online/)for visualizing the data.
Funding
This research received no external funding.
Author Information
Authors and Affiliations
Department of Oncology, First Teaching Hospital of Tianjin University of Traditional Chinese Medicine, Tianjin, 300380, China.
Zixin Li, Ying Zhang, Jiaqiao Pei, Huixin Chen& Yingying Huang
National Clinical Research Center for Chinese Medicine Acupuncture and Moxibustion, Tianjin, China
Zixin Li, Ying Zhang, Jiaqiao Pei, Huixin Chen& Yingying Huang
Graduate School of Tianjin University of Traditional Chinese Medicine, Tianjin, China.
Zixin Li, Jiaqiao Pei, Zhe Xu, Huixin Chen& Yingying Huang
School of Integrative Medicine, Tianjin University of Traditional Chinese Medicine, Tianjin, 301617, China.
Zhe xu
Contributions
Design, Z.X.L. and Y.Z.; Data Curation, Z.X.L. and Z.X; Writing, Z.X.L. and J.Q.P.; Inspection, H.X.C. and Y.Y.H. All authors have read and agreed to the published version of the manuscript.
Corresponding authors
Correspondence to Ying Zhang.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for Publication
Not applicable.
Institutional Review Board Statement
Not applicable.
Conflicts of Interest
The authors declare no conflict of interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Zixin Li and Ying Zhang contributed equally to this work and share first authorship.