- 1Department of Urology, The Affiliated Hospital of Qingdao University, Qingdao, China
- 2Department of Nutrition, The Affiliated Hospital of Qingdao University, Qingdao, China
This study aimed to identify critical cell cycle-related genes (CCRGs) in prostate cancer (PRAD) and to evaluate the clinical prognostic value of the gene panel selected. Gene set variation analysis (GSVA) of dysregulated genes between PRAD and normal tissues demonstrated that the cell cycle-related pathways played vital roles in PRAD. Patients were classified into four clusters, which were associated with recurrence-free survival (RFS). Moreover, 200 prognostic-related genes were selected using the Kaplan–Meier (KM) survival analysis and univariable Cox regression. The prognostic CCRG risk score was constructed using random forest survival and multivariate regression Cox methods, and their efficiency was validated in Memorial Sloan Kettering Cancer Center (MSKCC) and GSE70770. We identified nine survival-related genes: CCNL2, CDCA5, KAT2A, CHTF18, SPC24, EME2, CDK5RAP3, CDC20, and PTTG1. Based on the median risk score, the patients were divided into two groups. Then the functional enrichment analyses, mutational profiles, immune components, estimated half-maximal inhibitory concentration (IC50), and candidate drugs were screened of these two groups. In addition, the characteristics of nine hub CCRGs were explored in Oncomine, cBioPortal, and the Human Protein Atlas (HPA) datasets. Finally, the expression profiles of these hub CCRGs were validated in RWPE-1 and three PRAD cell lines (PC-3, C4-2, and DU-145). In conclusion, our study systematically explored the role of CCRGs in PRAD and constructed a risk model that can predict the clinical prognosis and immunotherapeutic benefits.
Introduction
Prostate cancer (PRAD) is the most common carcinoma in men (1) and a significant global health concern (2). The main treatment for PRAD is radical prostatectomy (3). Although most PRAD patients could benefit from this treatment, nearly 27%–53% of patients who have undergone this procedure will progress into advanced PRAD and castration-resistant prostate cancer (CRPC) (4, 5). Therefore, timely diagnosis of PRAD and exploring the detailed mechanisms involved in PRAD are critical for improving the prognosis of patients with PRAD.
The cell cycle is one of the vital biological processes in the organism (6). The cell cycle regulates the process of cell division and duplication of genetic materials (7), which is highly associated with the growth and proliferation of cancer cells. Increasing numbers of studies showed that certain genes and drugs serve as potential cycle regulators. For instance, in prostate cancer cells, upregulation of PHLDA3 inhibited cell proliferation by inducing cell cycle arrest at G1 via a decrease in AKT phosphorylation and activation of Wnt/β-catenin (8). The decreased expression of SMARCC1 dramatically accelerated prostate cancer cell proliferation by enhancing cell cycle progression (9). Platycodin D could promote sorafenib-induced apoptosis and cell cycle arrest in prostate cancer cells (10) and BTT-3033 attenuated prostate cancer cell viability and proliferation by cell cycle arrest (11). Thus, cell cycle arrest may be used as a novel therapeutic strategy. However, such research is a lengthy process whose results do not immediately translate into clinical practice.
In this study, we firstly identified the molecular hallmarks in normal samples and PRAD tissues using the Gene Set Enrichment Analysis (GSEA) method. The results showed that cell cycle-related pathways, such as DNA_REPAIR, E2F_TARGETs, G2M_CHECKPOINT, MYC _TARGETS_V1, and MYC_TARGETS_V2, were enriched in the PRAD tissues. Considering the critical roles of cell cycle-related genes (CCRGs) in the initiation and progression of cancers, then we hypothesized that CCRGs may provide a novel insight into the treatment and prognosis of PRAD. Therefore, there is an urgent need to explore the roles of CCRGs in the prognosis and treatment of PRAD patients.
Results
Functional Pathway Screening Using Gene Set Variation Analysis
The clinical information of 551 subjects, including 499 PRAD patients and 52 healthy volunteers, was downloaded from UCSC Xena. Based on The Cancer Genome Atlas (TCGA)-PRAD cohort data, gene set variation analysis (GSVA) results showed the CCRG sets, such as HALLMARK_DNA_REPAIR, HALLMARK_E2F_TARGETs, HALLMARK_G2M_CHECKPOINT, and HALLMARK_MYC_TARGETS_V1 (Figure 1A) were enriched in PRAD patients.
Figure 1 Identification of cell cycle-related DEGs and consensus clustering in TCGA-PRAD cohort. (A) The hallmark enrichment of prostate cancer compared with normal tissues. (B) The volcano plot of cell cycle-related DEGs between high- and low-risk groups (FDR < 0.05 and logFC > 0.5). (C) Consensus matrices for a solution with 4 clusters based on DEGs. (D) Kaplan–Meier curves show the prognostic value of the four patterns. DEGs, differentially expressed genes; TCGA-PRAD, The Cancer Genome Atlas—Prostate Cancer; FDR, false discovery rate.
Cluster Analysis Based on Cell Cycle-Related Genes
The limma (12) R package was used to detect differentially expressed CCRGs. The volcano map of the differentially expressed CCRGs showed that there were 79 upregulated genes and 131 downregulated genes (Figure 1B). TCGA-PRAD cohort could be divided into four clusters based on differentially expressed CCRGs (Figure 1C). Moreover, patients had significant differences among these four clusters (p < 0.039, Figure 1D).
Construction and Validation of Cell Cycle-Related Gene Prognostic Model
As for the 210 differentially expressed CCRGs, univariate Cox regression analysis showed that 57 CCRGs were significantly associated with recurrence-free survival (RFS) in PRAD (Table 1). Thus, these 57 CCRGs served as input to construct a random survival forest survival model. The out-of-bag (OOB) prediction error estimator indicated that the forest prediction error tended to be steady when the number of trees was nearly 400 (Figure 2A). During the hub gene selection process, the top 15 ranked genes in both minimal depth and VIMP were chosen for further model construction (Figure 2B). Finally, CCNL2, CDCA5, KAT2A, CHTF18, SPC24, EME2, CDK5RAP3, CDC20, and PTTG1 were included in the prognosis model construction, and the risk score of each patient was calculated based on the following formula with coefficients showed in Figure 2C: Risk score = (0.6097 * ExpCCNL2) + (0.5850 * ExpCDCA5) + (0.0006 * ExpKAT2A) + (−0.0005 * ExpCHTF18) + (0.1037 * ExpSPC24) + (0.0283 * ExpEME2) + (−0.1527 * ExpCDK5RAP3) + (0.2168 * ExpCDC20) + (−0.1638 * ExpPTTG1). Based on the median value of the risk score, TCGA-PRAD patients were divided into high and low risk group. Principal component analysis (PCA) indicated that the two risk groups were distributed in two directions (Figure 2D). The Kaplan–Meier (KM) curve showed that the high-risk group patients had poorer RFS than the low-risk group (Figure 2E, p < 0.001). The prognosis performance of the risk score for RFS was assessed by time-dependent receiver operating characteristic (ROC) curves, with the area under the curve (AUC) for 1, 3, and 5 years being 0.786, 0.739, and 0.679, respectively (Figure 2F). Moreover, the results of two independent cohorts showed that the high-risk group was associated with worse RFS (Memorial Sloan Kettering Cancer Center (MSKCC), p = 0.041, and GSE70770, p < 0.001) (Figures 3A, D), which were consistent with the results in TCGA-PRAD cohort. The AUCs in MSKCC were 0.771, 0.72, and 0.691 for 1, 3, and 5 years, respectively (Figure 3B). In the GSE70770 cohort, the AUCs for risk scores at 1, 3, and 5 years were 0.671, 0.712, and 0.770, respectively (Figure 3E). The mRNA expression profiles of nine hub CCRGs were differently expressed in the different risk groups in MSKCC and GSE70770 (Figures 3C, F). In addition, the univariate and multivariate Cox regression analyses showed that the risk score was an independent prognostic predictor for RFS in TCGA-PRAD and MSKCC (Table 2).
Figure 2 Gene selection and risk prognostic model of CCRGs based on TCGA cohort. (A) Estimation of the OOB prediction error rate based on the random forest. (B) The top 15 genes according to both minimal depth and variable importance. (C) The coefficient for genes of risk prognostic model. (D) The distribution of risk scores based on PCA. The high-risk group was annotated by yellow and the low-risk group by blue. (E) Kaplan–Meier curves of RFS stratified by risk score. (F) Time-dependent ROC curve analysis for risk score. CCRGs, cell cycle-related genes; TCGA, The Cancer Genome Atlas; OOB, out of bag; PCA, principal component analysis; RFS, recurrence-free survival; ROC, receiver operating characteristic.
Figure 3 Validation of the 9-gene signature in the MSKCC and GSE70770 cohort. (A–C) The distribution of RFS, RFS status, KM curves, and ROC curves in the MSKCC cohort. (D–F) The distribution of RFS, RFS status, KM curves, and ROC curves in the GSE70770 cohort. MSKCC, Memorial Sloan Kettering Cancer Center; RFS, recurrence-free survival; KM, Kaplan–Meier; ROC, receiver operating characteristic. ns, no significance. *p < 0.05, ***p < 0.001.
Differences in Genomic Alteration Profiles, Copy Number Variation, and Tumor Mutational Burden Between Cell Cycle-Related Gene Risk Groups
Considering that genetic alterations were associated with the prognostic outcome of many cancers (13, 14), we compared the genomes of high- and low-risk groups. The top 10 genes with the highest mutation frequency in TCGA-PRAD cohort, high-risk group, and low-risk group are shown in Figure 4A. The mutation rates of TP53 and FOXA1 were 17% and 9% in the high-risk group, but 6% and 4% in the low-risk group, respectively, which indicated that CCRGs were probably associated with TP53 and FOXA1 pathways. For the copy number variation (CNV) status, the results demonstrated that the high-risk group had high burden of amplification (pArm-Amp = 0.001, pFocal-Amp < 0.001) and deletion (pArm-del < 0.001, pFocal-del < 0.001) at both the arm and focal levels (Figure 4B). Furthermore, tumor mutational burden (TMB) was significantly higher in the high-risk group than low-risk group (p < 0.001, Figure 4C).
Figure 4 Mutational landscapes between high- and low-risk groups. (A) The distribution of frequently mutated genes in total, high-risk, and low-risk patients. (B) Arm-level, focal-level copy number amplification and deletion. (C) Tumor mutational burden difference.
Functional Enrichment Between High-Risk and Low-Risk Groups
Considering the prognostic risk model was constructed based on CCRGs, which were associated with cell proliferation, then mRNAsi and MKI67 were analyzed. The results revealed reduced mRNAsi in the low-risk group (p < 0.001, Figure 5A), and the risk score was positively associated with MKI67 (R = 0.590, p < 0.001, Figure 5B). Next we investigated the potential functions in the high- and low-risk groups using GSEA and GSVA methods. Hallmark analysis showed that cell cycle-related pathways were enriched in the high-risk group (Figure 5C). Additionally, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis using GSVA showed that cell cycle-related pathways such as DNA_REPLICATION, MISMATCH_REPAIR, and CELL_CYCLE were also enriched in the high-risk group (Figure 5D). In summary, the results of mRNAsi, MKI67, and functional enrichment all demonstrated that the high-risk group was associated with cell cycle-related pathways, meaning that the prognostic signature could be represented by the CCRGs.
Figure 5 The differences of mRNAsi, MKI67, and potential biological pathways of high- and low-risk groups. (A) mRNAsi. (B) The correlations between MKI67 and risk score. (C) The hallmark enrichment based on the GSEA method. (D) The KEGG enrichment calculated by GSVA method. MSI, microsatellite instability; GSEA, Gene Set Enrichment Analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSVA, gene set variation analysis.
Evaluation of Immune Cell Infiltration Characterization in Tumor Microenvironment
The prognosis of tumor patients was significantly associated with the tumor microenvironment, especially with immune cells (15–17). Therefore, we hypothesized that the distribution of immune cells and expression of immune checkpoint genes would be significantly different in the two risk groups. We evaluated 24 immune cells by CIBERSORT method in normal and tumor samples. Low abundance immune cells were excluded; thus, only 21 immune cell types were assessed (Figure 6A). Of these, the high-risk group had a higher proportion of infiltration by B-cell memory cells, Macrophage M0, Macrophage M2, T-cell follicular helper, and T-cell regulatory (Tregs). The relationships between risk score and common immune checkpoint genes showed that many checkpoint genes were more highly expressed in the high-risk group, including PDL2 (PFCDL1G2), CD48, CD44, and CD200, while TNFRSF4, TNFRSF14, TNFRSF18, TNFRSF25, NRP1, LAG3, and CTLA4 were reduced (Figure 6B). Moreover, several members of the nine hub genes demonstrated significantly positive relationships with B-cell memory, Macrophage M2, T-cell follicular helper, and T-cell regulatory (Tregs) infiltration (Figure 6C).
Figure 6 TME immune cell infiltration landscapes of different risk groups. (A) Differences of 24 TME infiltration cells based on CIBERSORT algorithm. (B) The mRNA expression profiles of common immune checkpoint genes. (C) The correlations between 9-risk genes and TME infiltration cell type. Red, positive; blue, negative. (D) The correlation between 9-risk genes and immune checkpoint molecules. (E) The correlation between REC, CTLA4, and PDCD1 (PD1). TME, tumor microenvironment. ns, no significance. *p < 0.05, **p < 0.01, ***p < 0.001.
We also found that the nine hub genes were positively associated with PDCD1 (PD-1) and CTLA4 (Figure 6D). Of these, CDC20 and PTTG1 were relatively significant (Figure 6E).
Screening of Immunotherapeutic Benefits, Estimated IC50, and Candidate Drugs
To explore the potential performance of the risk signature on immunotherapeutic benefits, we analyzed the IMvigor210 cohort, who received anti-PD-L1 immunotherapy. Notably, a high response rate of anti-PD-L1 therapy was associated with a higher risk core (p < 0.001, Figure 7A). Moreover, the high-risk group had a more favorable survival rate (p = 0.04, Figure 7B) and objective response to anti-PD-L1 than the low-risk group (p < 0.001, Figure 7C). Based on the Genomics of Drug Sensitivity in Cancer (GDSC) dataset, the IC50 value of 138 compounds showed that the low-risk group might be more sensitive to 51 compounds, such as dasatinib, DMOG, and MG.132; and the patients in the high-risk group were likely sensitive to 30 drugs, such as thapsigargin, bleomycin, and vinblastine (false discovery rate (FDR) < 0.05, Figure 7D). In addition, the CMap dataset was utilized to predict the candidate drugs for the risk signature. The results showed that thiostrepton, GW8150, phenoxybenzamine, chrysin, camptothecin, menadione, DL-thiorphan, and sanguinarine were the most potential target drugs due to their enrichment scores (<−0.90) and p < 0.05 (Table 3), of which the 2D conformers are displayed in Figure 7E.
Figure 7 The roles of risk scores in the prediction of immuno-/chemotherapeutic benefits and candidate drugs. (A) Risk scores in high- and low-risk groups with different anti-PD-1 clinical response statuses. p < 0.001. (B) KM curves for high- and low-risk groups in the IMvigor210 cohort. Log-rank test, p = 0.032. (C) Rate of clinical response rate to anti-PDL1 immunotherapy in high- and low-risk groups in the IMvigor210 cohort. (D) Estimated IC50 for 138 compounds based on GDSC dataset. (E) 2D conformer of six significant candidate drugs. KM, Kaplan–Meier; GDSC, Genomics of Drug Sensitivity in Cancer.
The Characteristics of Hub Genes in the Oncomine and cBioPortal Datasets
We used Oncomine dataset to search the expression profiles of nine hub genes whose expression was increased (p < 0.05) in several types of cancer (Figure 8A), especially in PRAD patients, consistently with TCGA-PRAD cohort (Figure 8B). Moreover, the correlation among these hub genes was very high; for example, the correlation coefficients between PTTG1 and SPC24 was 0.87 (p < 0.05, Figure 8D). In addition, the genetic alteration status based on cBioPortal showed that the genetic alteration rates of these hub genes were less than 5% (Figure 8C), and several items belonged to the CNV change. The correlation between CNV and mRNA expression of EME2 was 0.44 (Figure 8E), while that of the others was less than 0.3. The methylation levels of nine genes (expect CCNL2, as data were absent from cBioPortal) were negatively associated with mRNA expression (p < 0.05, Figure 8E).
Figure 8 The mRNA expression patterns, genomic alterations, and methylation of nine hub genes. (A) The overview of nine hub genes in the Oncomine database. (B) The mRNA expression profiles of nine hub genes in TCGA-PRAD. (C) The genetic alterations of nine hub genes based on cBioPortal. (D) The correlations between these nine hub genes. (E) The correlations between mRNA expression and methylation, and copy number variation of nine hub genes. TCGA-PRAD, The Cancer Genome Atlas-Prostate Cancer. NA, not avaiable. *p < 0.05, **p < 0.01, ***p < 0.001.
Verification of Hub Genes Based on the Human Protein Atlas Datasetand RT-qPCR
To verify the reliability of these hub genes, we detected the protein levels from the Human Protein Atlas (HPA) website in normal samples and PRAD tissues. The results showed that eight proteins (CCNL2, CDCA5, CDC20, CDK4RAP3, EME2, KAT2A, PTTG1, and SPC24; note that CHTF18 was absent) were significantly dysregulated in PRAD tissues compared with normal prostate tissues (Figure 9A). To further confirm the expression levels based on bioinformatics analysis, the mRNA expression profiles of these hub genes were detected by RT-qPCR from normal prostatic epithelial cells (RWPE-1) and prostate cancer cells (PC-3, C4-2, and DU-145). The results showed that all nine genes were significantly upregulated in PC-3 and DU-145 compared with RWPE-1 (Figure 9B), which were in accordance with the contents from TCGA and Oncomine datasets.
Figure 9 The IHC expression pattern based on HPA dataset and mRNA levels by qRT-PCR of CCNL2, CDCA5, CDC20, CDK4RAP3, EME2, KAT2A, PTTG1, REC8, and SPC24. (A) The IHC results of nine hub genes in prostate cancer and normal tissues based on HPA. (B) The mRNA level of nine hub genes in normal prostatic epithelial cell (RWPE-1) and prostate cancer cell lines. IHC, immunohistochemistry; HPA, Human Protein Atlas. ns, no significance. *p < 0.05, **p < 0.01, ***p < 0.001.
Discussion
The cell cycle is an essentially biological process (18, 19). Under normal conditions, cells proliferate only in response to mitotic signals, which are more significant for normal organisms, while the proliferation of cancer cells is out of control (7). This suggested that the proliferation of cancer cells is associated with the dysregulation of proliferation-related signals, which could be controlled by the cell cycle. Several previous studies showed that dysregulated genes could give rise to the expression of key factors involved in cancer cell cycle. Therefore, we hypothesized that CCRGs have excellent performance in PRAD.
Recent studies have shown that high throughput RNA sequences and microarray profiles have been utilized to develop signatures for the outcome events of several clinical diseases (20). In our study, we firstly used the GSVA method to explore the differential hallmark pathways, and results showed that CCRGs were highly enriched in PRAD patients (Figure 1A). Based on the differential CCRGs, TCGA-PRAD patients could be classified into four clusters, which were significantly associated with RFS (Figures 1C, D). With regard to risk model construction, the random survival method was utilized to select the hub CCRGs based on variable importance (VIMP) and minimal depth (Figures 2A, B). Then nine hub genes were included in the next step, and multivariate Cox regression was applied to calculate risk scoring. The performances on the prediction of prognostic ability were validated in the MSKCC and GSE70770 datasets (Figure 3, Table 2). Moreover, hallmark and KEGG pathways (Figures 5D, E) showed that the high-risk group was enriched in cell cycle-related pathways. As genomic instability takes critical roles in the development of cancers (21, 22), then genetic alterations, such as mutation status, CNV load, and TMB, were analyzed in TCGA-PRAD cohort. The results showed that the high-risk group was significantly associated with high mutation rates, especially for TP53 and FOXA1 (Figure 4A), CNV load (Figure 4B), and TMB (Figure 4C), which may help tailor personalized treatment.
In recent years, numerous previous studies indicated that cell cycle gene signatures have the potential for evaluating immune cell infiltration, immune evasion, and immune responses (23–25). However, the relationship between cell cycle-related signatures and tumor immune situation in PRAD was not explored before. Therefore, we compared the immune cell abundance based on CIBERSORT and common checkpoint genes in TCGA-PRAD cohort. The high-risk group had inflammatory infiltrates with a higher proportion of M2 macrophages and T-cell regulatory (Tregs) (Figure 6A), which were associated with immune evasion (26, 27). Considering the essential roles of immune checkpoint genes in immunotherapy response, immune checkpoint genes were differently distributed between low- and high-risk groups, such as LAG3, CTLA4, NRP1, and CD276 (Figure 6B). Increased evidence demonstrated that single genes could reshape the tumor immune microenvironment and immune cell infiltration (28, 29). We revealed obviously positive correlations between these nine hub genes and the expression of B-cell memory cells, Macrophage M2, T-cell follicular helper, and T-cell regulatory (Figure 6C), consistent with the results of immune cell infiltration and risk groups. We also noted that most of these nine hub genes were positively correlated with PDCD1 and CTLA4 (Figure 6D), suggesting that these genes could mediate immune evasion and response to immunotherapy. In the IMvigor210 cohort with anti-PD-L1 immunotherapy, we found that patients with high-risk scores might benefit from anti-PD1 immunotherapy (Figures 7A–C).
Additionally, the analysis of drug sensitivity based on the GDSC dataset demonstrated that the CCRG risk signature might be useful for therapeutic applications. Several drugs, such as dasatinib, MG.132, lapatinib, and docetaxel, responded differently between the low- and high-risk groups (Figure 7D), which suggested that CCRGs influenced drug response to chemotherapy and targeted treatment. Fortunately, we found 8 drugs with p < 0.001 and enrichment less than −0.9 (Table 3), including GW-8510, phenoxybenzamine, thiostrepton, chrysin, camptothecin menadione, DL-thiorphan, and sanguinarine (Figure 7E, Table 3). GW-8510, an inhibitor of CDK2, has a similar effect to gemcitabine in inhibiting pancreatic cancer cells (30, 31). Phenoxybenzamine, an alpha blocker, has been used to inhibit histone deacetylases (32, 33) in human cancer cells. Thiostrepton, a natural antibiotic produced by bacteria, could induce upregulation of several heat shock proteins in various human cancer cells (34) and inhibit cancer stem cell growth (35). Chrysin could inhibit cancer growth via induction of apoptosis, alteration of the cell cycle, and inhibition of angiogenesis without causing any toxicity to normal cells (36, 37). Camptothecin (38), menadione (39, 40), DL-thiorphan (41), and sanguinarine (42) are also reported to have a strong relationship with cancer therapy. Moreover, many previous studies have been reported to have anticancer effects in prostate cancer cells, such as chrysin (43), phenoxybenzamine, thiostrepton (44), sanguinarine (45), camptothecin (46), menadione (47), and sanguinarine (48, 49). Finally, the mRNA expression profiles of nine hub genes were explored in TCGA-PRAD and validated in Oncomine, the HPA dataset, and prostate cancer cells based on RT-qPCR.
Although our study have comprehensive analyzed the CCRGs in PRAD, there were still some limitations. Firstly, we constructed and validated the CCRGs risk model only including public datasets, so more real-world data of PRAD are needed to verify the model’s prognostic performance. Secondly, we only validated the mRNA expression profiles of nine hub CCRGs in prostate cancer cells; more experimental studies are still necessary to confirm for clinical application, and more underlying mechanisms of these genes should be explored for further research.
In conclusion, this preliminary research of CCRGs in PRAD patients has profiled the expression levels and genetic alterations of CCRGs in PRAD, which may open up the development of novel drugs against PRAD. The prognostic model based on nine hub CCRGs was strongly correlated with high CNV load, TMB load, mRNAsi, infiltration of different types of immune cells, and chemo-/immunotherapy response, which may provide novel ideas for PRAD with patients chemo-/immunotherapy response.
Materials and Methods
Prostate Cancer Datasets and Samples
We used TCGA-PRAD cohort data, which includes mRNA expression data, clinicopathological features, and RFS data from the UCSC Xena (https://xenabrowser.net/datapages/) dataset. The raw read counts of RNA-seq were transformed into transcripts per kilobase million (TPM) values. The DNA methylation information and genetic mutations were collected from cBioPortal (50).
The external validation cohort, including the MSKCC (GSE21032) and GSE70770, were log2 normalized microarray matrix downloaded from http://cbio.mskcc.org/cancergenomics/prostate/data/ and Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) dataset.
Gene Set Variation Analysis, Gene Set Enrichment Analysis, Difference Analysis, and Consensus Clustering
The hallmark enrichment between PRAD and normal tissues and the KEGG enrichment for high and low risk were performed using the GSVA method. The hallmark enrichment between the high- and low-risk groups was generated by the GSEA method. Then, a panel of 1,875 CCRGs were recognized from MSigDB (51), as previously described (7). The differentially expressed CCRGs were calculated using the limma package in R with FDR < 0.05 and |logFC| > 0.5. Next, a consensus clustering algorithm was utilized to evaluate the prognostic ability of CCRGs using the ConsensusClusterPlus (52) R package.
Calculation of the Cell Cycle-Related Gene Score for Prostate Cancer Patients
In TCGA-PRAD cohort, differentially expressed CCRGs were processed to univariate Cox regression analysis. The random forest algorithm was applied to select the candidate genes. Only genes included in the top 15 lists for both minimal depth and VIMP were selected. Then, PCA was also performed in R. Finally, the candidate genes were brought into the multivariate Cox regression analysis to build the prognosis model. The risk score of the CCRGs for each PRAD patient was computed based on the mRNA expression of selected CCRGs weighted by the multivariate Cox regression coefficient. Based on the median CCRG score, the patients were divided into high- and low-risk groups. Next, PCA was performed to explore the distribution of different risk groups. Finally, the prognostic values of CCRG score were evaluated using KM curve and the 1-, 3-, and 5-year ROC curves. A heatmap plot was utilized to display the expression profiles for the different groups.
As for the MSKCC and GSE70770 cohorts, the KM curve, ROC curve, and heatmap plot were also drawn to validate the prognostic performance.
Identification of Somatic Alteration, Copy Number Variation, Tumor Mutational Burden, Tumor Stemness, and Immune Cell Infiltration in Different Groups
Somatic mutations data were downloaded from TCGA GDC (https://portal.gdc.cancer.gov/) using “TCGABiolinks” (53) R package. Copy number alterations data were calculated by GISTIC2.0 from GDAC Firehose (https://gdac.broadinstitute.org). The total number of mutations per megabyte of tumor tissue (TMB) was generated by the number of nonsynonymous mutations per million. Tumor stemness was previously study calculated by using the one-class logistic regression (OCLR) method (54). The CIBERSORT (55) method was processed to infer the proportion of immune cell infiltration in different tumor groups. The connectivity Map (CMap) (56) was utilized to screen potential candidate drugs.
Screening of Immunotherapy, Chemotherapy Responses, and Potential Drugs
The IMvigor210 dataset was downloaded from IMvigor210CoreBiologies (57) to evaluate the predictive power of the CCRG scores. The estimated IC50 of 138 compounds in the GDSC (58) website was calculated using the “pRRophetic” (59) R package.
The Characteristic of Nine Hub Cell Cycle-Related Genes
The mRNA expression levels of these nine hub CCRGs among pan-cancer were retrieved from Oncomine (https://www.oncomine.org/) with a threshold p-value <0.05 (60). The mRNA expression profiles of nine hub CCRGs in TCGA-PRAD were downloaded from UCSC Xena (61). The mutation alteration data, CNV, and DNA methylation data were downloaded from cBioPortal (http://www.cbioportal.org). The immunohistochemistry (IHC) results of nine hub genes were collected from the HPA (https://www.proteinatlas.org/) (62).
Detection of the mRNA Levels of Hub Genes Using Real-Time Quantitative PCR
The human RWPE-1 cell lines were cultured in K-SFM (Biotecnómica, Porto, Portugal) medium, and human epithelial PRAD cell lines (PC-3, C4-2, and DU145) were cultured with DMEM 1640 medium. Total RNA, cDNA, and RT-qPCR were operated according to the manufacturer’s protocol. Independent experiments were performed in triplicate, and GAPDH served as the internal control. The primers are as follows:
CCNL2 gene 5′-GTACTCCGGGGTGCTCATC-3′ (sense) and 5′-GAGGTCGGTCTCTGTGTCG-3′ (antisense).
CDCA5 gene 5′-GAGGTCCCAGCGGAAATCAG-3′ (sense) and 5′-TCTTTAAGACGATGGGCTTTCTG-3′ (antisense).
KAT2A gene 5′-GCAAGGCCAATGAAACCTGTA-3′ (sense) and 5′-TCCAAGTGGGATACGTGGTCA-3′ (antisense).
SPC24 gene 5′-GCCTTCCGCGACATAGAGG-3′ (sense) and 5′-CCTGCTCCTTCGCATTGAGA-3′ (antisense).
EME2 gene 5′-CGCCGTTACCAAGGCTCTC-3′ (sense) and 5′-GCTGACCCGACTGAACTGC-3′ (antisense).
CHTF18 gene 5′-GAGCCGACTGACGGTCAAG-3′ (sense) and 5′-CGGTTGGTGAAGTCATCACTG-3′ (antisense).
CDK5RAP3 gene 5′-GAGTCTGGTGCTGACGATCC-3′ (sense) and 5′-TGTGAAGAGTATCGGCCAAAAAT-3′ (antisense).
CDC20 gene 5′-GCACAGTTCGCGTTCGAGA-3′ (sense) and 5′-CTGGATTTGCCAGGAGTTCGG-3′ (antisense).
PTTG1 gene 5′-ACCCGTGTGGTTGCTAAGG-3′ (sense) and 5′-ACGTGGTGTTGAAACTTGAGAT-3′ (antisense).
Statistical Analysis
All statistical analyses were executed on the R platform (Version 4.1.0). A significant difference between the two groups was assessed using the Wilcoxon rank test, and among three or more groups, it was calculated using one-way ANOVA and Kruskal–Wallis tests. In addition, KM curves were produced, and a log-rank test was used, together with hazard ratios (HR) with 95% CI using univariate and multivariate Cox regression analyses, as well as correlation tests with Spearman’s method. All statistical p-values were two-sided, with p < 0.05 statistically significant.
Data Availability Statement
The data about TCGA, GEO, CMap and HPA dataset are publicly available, other original contributions presented in the study are included in the article. Further inquiries can be directed to the corresponding author.
Author Contributions
ZL and RP: conceptualization, supervision, methodology, and writing—review and editing. WL and YL: data curation, formal analysis, and writing—original draft. ZL, RP, and YL: validation, visualization, and investigation. All authors contributed to the article and approved the submitted version.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1. Zhang E, He J, Zhang H, Shan L, Wu H, Zhang M, et al. Immune-Related Gene-Based Novel Subtypes to Establish a Model Predicting the Risk of Prostate Cancer. Front Genet (2020) 11:595657. doi: 10.3389/fgene.2020.595657
2. Li J, Du H, Chen W, Qiu M, He P, Ma Z. Identification of Potential Autophagy-Associated lncRNA in Prostate Cancer. Aging (Albany N Y) (2021) 13:13153–65. doi: 10.18632/aging.202997
3. Taneja SS. Re: Radical Prostatectomy or Watchful Waiting in Prostate Cancer-29-Year Follow-Up. J Urol (2019) 202:210–1. doi: 10.1097/JU.0000000000000323
4. Van den Broeck T, van den Bergh RCN, Arfi N, Gross T, Moris L, Briers E, et al. Prognostic Value of Biochemical Recurrence Following Treatment With Curative Intent for Prostate Cancer: A Systematic Review. Eur Urol (2019) 75:967–87. doi: 10.1016/j.eururo.2018.10.011
5. Liu H, Gao L, Xie T, Li J, Zhai TS, Xu Y. Identification and Validation of a Prognostic Signature for Prostate Cancer Based on Ferroptosis-Related Genes. Front Oncol (2021) 11:623313. doi: 10.3389/fonc.2021.623313
7. Zhang Z, Ji M, Li J, Wu Q, Huang Y, He G, et al. Molecular Classification Based on Prognostic and Cell Cycle-Associated Genes in Patients With Colon Cancer. Front Oncol (2021) 11:636591. doi: 10.3389/fonc.2021.636591
8. Ma S, Quan P, Yu C, Fan X, Yang S, Jia W, et al. PHLDA3 Exerts an Antitumor Function in Prostate Cancer by Down-Regulating Wnt/beta-Catenin Pathway via Inhibition of Akt. Biochem Biophys Res Commun (2021) 571:66–73. doi: 10.1016/j.bbrc.2021.07.038
9. Xiao ZM, Lv DJ, Yu YZ, Wang C, Xie T, Wang T, et al. SMARCC1 Suppresses Tumor Progression by Inhibiting the PI3K/AKT Signaling Pathway in Prostate Cancer. Front Cell Dev Biol (2021) 9:678967. doi: 10.3389/fcell.2021.678967
10. Lu Z, Song W, Zhang Y, Wu C, Zhu M, Wang H, et al. Combined Anti-Cancer Effects of Platycodin D and Sorafenib on Androgen-Independent and PTEN-Deficient Prostate Cancer. Front Oncol (2021) 11:648985. doi: 10.3389/fonc.2021.648985
11. Salemi Z, Azizi R, Fallahian F, Aghaei M. Integrin Alpha2beta1 Inhibition Attenuates Prostate Cancer Cell Proliferation by Cell Cycle Arrest, Promoting Apoptosis and Reducing Epithelial-Mesenchymal Transition. J Cell Physiol (2021) 236:4954–65. doi: 10.1002/jcp.30202
12. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res (2015) 43:e47. doi: 10.1093/nar/gkv007
13. Aoki K, Nakamura H, Suzuki H, Matsuo K, Kataoka K, Shimamura T, et al. Prognostic Relevance of Genetic Alterations in Diffuse Lower-Grade Gliomas. Neuro Oncol (2018) 20:66–77. doi: 10.1093/neuonc/nox132
14. Fraser M, Sabelnykova VY, Yamaguchi TN, Heisler LE, Livingstone J, Huang V, et al. Genomic Hallmarks of Localized, non-Indolent Prostate Cancer. Nature (2017) 541:359–64. doi: 10.1038/nature20788
15. Galluzzi L, Humeau J, Buque A, Zitvogel L, Kroemer G. Immunostimulation With Chemotherapy in the Era of Immune Checkpoint Inhibitors. Nat Rev Clin Oncol (2020) 17:725–41. doi: 10.1038/s41571-020-0413-z
16. Hinshaw DC, Shevde LA. The Tumor Microenvironment Innately Modulates Cancer Progression. Cancer Res (2019) 79:4557–66. doi: 10.1158/0008-5472.CAN-18-3962
17. Bader JE, Voss K, Rathmell JC. Targeting Metabolism to Improve the Tumor Microenvironment for Cancer Immunotherapy. Mol Cell (2020) 78:1019–33. doi: 10.1016/j.molcel.2020.05.034
18. Zhao Z, Dong Q, Liu X, Wei L, Liu L, Li Y, et al. Dynamic Transcriptome Profiling in DNA Damage-Induced Cellular Senescence and Transient Cell-Cycle Arrest. Genomics (2020) 112:1309–17. doi: 10.1016/j.ygeno.2019.07.020
19. Kamal KY, van Loon J, Medina FJ, Herranz R. Differential Transcriptional Profile Through Cell Cycle Progression in Arabidopsis Cultures Under Simulated Microgravity. Genomics (2019) 111:1956–65. doi: 10.1016/j.ygeno.2019.01.007
20. Yin H, Zhang C, Gou X, He W, Gan D. Identification of a 13mrna Signature for Predicting Disease Progression and Prognosis in Patients With Bladder Cancer. Oncol Rep (2020) 43:379–94. doi: 10.3892/or.2019.7429
21. Peng B, Li H, Na R, Lu T, Li Y, Zhao J, et al. Identification of a Novel Prognostic Signature of Genome Instability-Related LncRNAs in Early Stage Lung Adenocarcinoma. Front Cell Dev Biol (2021) 9:706454. doi: 10.3389/fcell.2021.706454
22. Tubbs A, Nussenzweig A. Endogenous DNA Damage as a Source of Genomic Instability in Cancer. Cell (2017) 168:644–56. doi: 10.1016/j.cell.2017.01.002
23. Shi R, Bao X, Rogowski P, Schafer C, Schmidt-Hegemann NS, Unger K, et al. Establishment and Validation of an Individualized Cell Cycle Process-Related Gene Signature to Predict Cancer-Specific Survival in Patients With Bladder Cancer. Cancers (Basel) (2020) 12(5):1146. doi: 10.3390/cancers12051146
24. Oshi M, Takahashi H, Tokumaru Y, Yan L, Rashid OM, Matsuyama R, et al. G2M Cell Cycle Pathway Score as a Prognostic Biomarker of Metastasis in Estrogen Receptor (ER)-Positive Breast Cancer. Int J Mol Sci (2020) 21(8):2921. doi: 10.3390/ijms21082921
25. Petroni G, Formenti SC, Chen-Kiang S, Galluzzi L. Immunomodulation by Anticancer Cell Cycle Inhibitors. Nat Rev Immunol (2020) 20:669–79. doi: 10.1038/s41577-020-0300-y
26. Nishikawa H, Koyama S. Mechanisms of Regulatory T Cell Infiltration in Tumors: Implications for Innovative Immune Precision Therapies. J Immunother Cancer (2021) 9(7):e002591. doi: 10.1136/jitc-2021-002591
27. Chen J, Cao X, Li B, Zhao Z, Chen S, Lai SWT, et al. Warburg Effect Is a Cancer Immune Evasion Mechanism Against Macrophage Immunosurveillance. Front Immunol (2020) 11:621757. doi: 10.3389/fimmu.2020.621757
28. Zhou L, Jiang Y, Liu X, Li L, Yang X, Dong C, et al. Promotion of Tumor-Associated Macrophages Infiltration by Elevated Neddylation Pathway via NF-kappaB-CCL2 Signaling in Lung Cancer. Oncogene (2019) 38:5792–804. doi: 10.1038/s41388-019-0840-4
29. Del Prete A, Sozio F, Schioppa T, Ponzetta A, Vermi W, Calza S, et al. The Atypical Receptor CCRL2 Is Essential for Lung Cancer Immune Surveillance. Cancer Immunol Res (2019) 7:1775–88. doi: 10.1158/2326-6066.CIR-19-0168
30. Chen P, Wu JN, Shu Y, Jiang HG, Zhao XH, Qian H, et al. Gemcitabine Resistance Mediated by Ribonucleotide Reductase M2 in Lung Squamous Cell Carcinoma is Reversed by GW8510 Through Autophagy Induction. Clin Sci (Lond) (2018) 132:1417–33. doi: 10.1042/CS20180010
31. Gencalp D, Atay S, Aydin HH, Ak H. Synergistic Interactions Between GW8510 and Gemcitabine in an In Vitro Model of Pancreatic Cancer. Anticancer Agents Med Chem (2021) 16:2204–15. doi: 10.2174/1871520621999210104201553
32. Inchiosa MA Jr. Anti-Tumor Activity of Phenoxybenzamine and its Inhibition of Histone Deacetylases. PLoS One (2018) 13:e0198514. doi: 10.1371/journal.pone.0198514
33. Lin XB, Jiang L, Ding MH, Chen ZH, Bao Y, Chen Y, et al. Anti-Tumor Activity of Phenoxybenzamine Hydrochloride on Malignant Glioma Cells. Tumour Biol (2016) 37:2901–8. doi: 10.1007/s13277-015-4102-y
34. Sandu C, Ngounou Wetie AG, Darie CC, Steller H. Thiostrepton, a Natural Compound That Triggers Heat Shock Response and Apoptosis in Human Cancer Cells: A Proteomics Investigation. Adv Exp Med Biol (2014) 806:443–51. doi: 10.1007/978-3-319-06068-2_21
35. Huang TH, Wu ATH, Cheng TS, Lin KT, Lai CJ, Hsieh HW, et al. In Silico Identification of Thiostrepton as an Inhibitor of Cancer Stem Cell Growth and an Enhancer for Chemotherapy in non-Small-Cell Lung Cancer. J Cell Mol Med (2019) 23:8184–95. doi: 10.1111/jcmm.14689
36. Kasala ER, Bodduluru LN, Madana RM, Attira KV, Gogoi R, Barua CC. Chemopreventive and Therapeutic Potential of Chrysin in Cancer: Mechanistic Perspectives. Toxicol Lett (2015) 233:214–25. doi: 10.1016/j.toxlet.2015.01.008
37. Salama AAA, Allam RM. Promising Targets of Chrysin and Daidzein in Colorectal Cancer: Amphiregulin, CXCL1, and MMP-9. Eur J Pharmacol (2021) 892:173763. doi: 10.1016/j.ejphar.2020.173763
38. Venditto VJ, Simanek EE. Cancer Therapies Utilizing the Camptothecins: A Review of the In Vivo Literature. Mol Pharm (2010) 7:307–49. doi: 10.1021/mp900243b
39. Kishore C, Sundaram S, Karunagaran D. Vitamin K3 (Menadione) Suppresses Epithelial-Mesenchymal-Transition and Wnt Signaling Pathway in Human Colorectal Cancer Cells. Chem Biol Interact (2019) 309:108725. doi: 10.1016/j.cbi.2019.108725
40. Semkova S, Zhelev Z, Miller T, Sugaya K, Aoki I, Higashi T, et al. Menadione/Ascorbate Induces Overproduction of Mitochondrial Superoxide and Impairs Mitochondrial Function in Cancer: Comparative Study on Cancer and Normal Cells of the Same Origin. Anticancer Res (2020) 40:1963–72. doi: 10.21873/anticanres.14151
41. Mizerska-Kowalska M, Kreczko-Kurzawa J, Zdzisinska B, Czerwonka A, Slawinska-Brych A, Mackiewicz Z, et al. Neutral Endopeptidase (NEP) Inhibitors - Thiorphan, Sialorphin, and its Derivatives Exert Anti-Proliferative Activity Towards Colorectal Cancer Cells In Vitro. Chem Biol Interact (2019) 307:105–15. doi: 10.1016/j.cbi.2019.04.033
42. Zhang H, Zhang J, Venkat PS, Gu C, Meng Y. Sanguinarine Exhibits Potent Efficacy Against Cervical Cancer Cells Through Inhibiting the STAT3 Pathway In Vitro and In Vivo. Cancer Manage Res (2019) 11:7557–66. doi: 10.2147/CMAR.S212744
43. Ryu S, Lim W, Bazer FW, Song G. Chrysin Induces Death of Prostate Cancer Cells by Inducing ROS and ER Stress. J Cell Physiol (2017) 232:3786–97. doi: 10.1002/jcp.25861
44. Pandit B, Gartel AL. New Potential Anti-Cancer Agents Synergize With Bortezomib and ABT-737 Against Prostate Cancer. Prostate (2010) 70:825–33. doi: 10.1002/pros.21116
45. Samarghandian S, Afshari JT, Davoodi S. Chrysin Reduces Proliferation and Induces Apoptosis in the Human Prostate Cancer Cell Line Pc-3. Clinics (Sao Paulo Brazil) (2011) 66:1073–9. doi: 10.1590/S1807-59322011000600026
46. Yuan X, Liu L, Wang W, Gao YR, Zhang D, Jia TT, et al. Development of (G3-C12)-Mediated Camptothecin Polymeric Prodrug Targeting to Galectin-3 Receptor Against Androgen-Independent Prostate Cancer. Int J Pharm (2020) 580:119123. doi: 10.1016/j.ijpharm.2020.119123
47. Tomasetti M, Nocchi L, Neuzil J, Goodwin J, Nguyen M, Dong L, et al. Alpha-Tocopheryl Succinate Inhibits Autophagic Survival of Prostate Cancer Cells Induced by Vitamin K3 and Ascorbate to Trigger Cell Death. PLoS One (2012) 7:e52263. doi: 10.1371/journal.pone.0052263
48. Sun M, Liu C, Nadiminty N, Lou W, Zhu Y, Yang J, et al. Inhibition of Stat3 Activation by Sanguinarine Suppresses Prostate Cancer Cell Growth and Invasion. Prostate (2012) 72:82–9. doi: 10.1002/pros.21409
49. Rahman A, Pallichankandy S, Thayyullathil F, Galadari S. Critical Role of H2O2 in Mediating Sanguinarine-Induced Apoptosis in Prostate Cancer Cells via Facilitating Ceramide Generation, ERK1/2 Phosphorylation, and Par-4 Cleavage. Free Radic Biol Med (2019) 134:527–44. doi: 10.1016/j.freeradbiomed.2019.01.039
50. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, et al. Integrative Analysis of Complex Cancer Genomics and Clinical Profiles Using the Cbioportal. Sci Signal (2013) 6:pl1. doi: 10.1126/scisignal.2004088
51. Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) Hallmark Gene Set Collection. Cell Syst (2015) 1:417–25. doi: 10.1016/j.cels.2015.12.004
52. Wilkerson MD, Hayes DN. ConsensusClusterPlus: A Class Discovery Tool With Confidence Assessments and Item Tracking. Bioinformatics (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170
53. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: An R/Bioconductor Package for Integrative Analysis of TCGA Data. Nucleic Acids Res (2016) 44:e71. doi: 10.1093/nar/gkv1507
54. Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, et al. Machine Learning Identifies Stemness Features Associated With Oncogenic Dedifferentiation. Cell (2018) 173:338–54.e315. doi: 10.1016/j.cell.2018.03.034
55. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12:453–7. doi: 10.1038/nmeth.3337
56. Lamb J, Crawford ED, Peck D, Modell JW, Blat IC, Wrobel MJ, et al. The Connectivity Map: Using Gene-Expression Signatures to Connect Small Molecules, Genes, and Disease. Science (2006) 313:1929–35. doi: 10.1126/science.1132939
57. Necchi A, Joseph RW, Loriot Y, Hoffman-Censits J, Perez-Gracia JL, Petrylak DP, et al. Atezolizumab in Platinum-Treated Locally Advanced or Metastatic Urothelial Carcinoma: Post-Progression Outcomes From the Phase II IMvigor210 Study. Ann Oncol (2017) 28:3044–50. doi: 10.1093/annonc/mdx518
58. Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, et al. Genomics of Drug Sensitivity in Cancer (GDSC): A Resource for Therapeutic Biomarker Discovery in Cancer Cells. Nucleic Acids Res (2013) 41:D955–61. doi: 10.1093/nar/gks1111
59. Geeleher P, Cox N, Huang RS. Prrophetic: An R Package for Prediction of Clinical Chemotherapeutic Response From Tumor Gene Expression Levels. PLoS One (2014) 9:e107468. doi: 10.1371/journal.pone.0107468
60. Rhodes DR, Yu J, Shanker K, Deshpande N, Varambally R, Ghosh D, et al. ONCOMINE: A Cancer Microarray Database and Integrated Data-Mining Platform. Neoplasia (2004) 6:1–6. doi: 10.1016/S1476-5586(04)80047-2
61. Goldman MJ, Craft B, Hastie M, Repecka K, McDade F, Kamath A, et al. Visualizing and Interpreting Cancer Genomics Data via the Xena Platform. Nat Biotechnol (2020) 38:675–8. doi: 10.1038/s41587-020-0546-8
Keywords: prostate cancer, GSVA, cell cycle, recurrence-free survival, immunotherapy
Citation: Liu Z, Pan R, Li W and Li Y (2022) Comprehensive Analysis of Cell Cycle-Related Genes in Patients With Prostate Cancer. Front. Oncol. 11:796795. doi: 10.3389/fonc.2021.796795
Received: 17 October 2021; Accepted: 06 December 2021;
Published: 11 January 2022.
Edited by:
Yubing Zhou, First Affiliated Hospital of Zhengzhou University, ChinaReviewed by:
Lei Gao, Second Hospital of Hebei Medical University, ChinaLifeng Zhang, Changzhou No. 2 People’s Hospital, China
Copyright © 2022 Liu, Pan, Li and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yanjiang Li, bHlqMjAwMTM1M0AxNjMuY29t