Qian Huang
Xingyu Peng
Qingqing Li
Jinfeng Zhu
Ju Xue4
Hua Jiang- 1Department of General Practice, Shanghai East Hospital, Tongji University School of Medicine, Shanghai, China
- 2Department of Geriatrics, Shanghai East Hospital, Tongji University School of Medicine, Shanghai, China
- 3Department of General Surgery, The Second Affiliated Hospital of Nanchang University, Nanchang, China
- 4Department of Pathology, Hunan Provincial People’s Hospital, The First Affiliated Hospital of Hunan Normal University, Changsha, China
Background: Treatment of cancer with pyroptosis is an emerging strategy. Molecular subtypes based on pyroptosis-related genes(PRGs) seem to be considered more conducive to individualized therapy. It is meaningful to construct a pyroptosis molecular subtypes-related prognostic signature (PMSRPS) to predict the overall survival (OS) of patients with pancreatic adenocarcinoma(PAAD) and guide treatment.
Methods: Based on the transcriptome data of 23 PRGs, consensus clustering was applied to divide the TCGA and GSE102238 combined cohort into three PRGclusters. Prognosis-related differentially expressed genes(DEGs) among PRGclusters were subjected to LASSO Cox regression analysis to determine a PMSRPS. External cohort and in vitro experiments were conducted to verify this PMSRPS. The CIBERSORT algorithm, the ESTIMATE algorithm and the Immunophenoscore (IPS) were used to analyze the infiltrating abundance of immune cells, the tumor microenvironment (TME), and the response to immunotherapy, respectively. Wilcoxon analysis was used to compare tumor mutational burden (TMB) and RNA stemness scores (RNAss) between groups. RT-qPCR and in vitro functional experiments were used for evaluating the expression and function of SFTA2.
Results: Based on three PRGclusters, 828 DEGs were obtained and a PMSRPS was subsequently constructed. In internal and external validation, patients in the high-risk group had significantly lower OS than those in the low-risk group and PMSRPS was confirmed to be an independent prognostic risk factor for patients with PAAD with good predictive performance. Immune cell infiltration abundance and TME scores indicate patients in the high-risk group have typical immunosuppressive microenvironment characteristics. Analysis of IPS suggests patients in the high-risk group responded better to novel immune checkpoint inhibitors (ICIs) than PD1/CTLA4. The high-risk group had higher TMB and RNAss. In addition, 10 potential small-molecule compounds were screened out. Finally, we found that the mRNA expression of SFTA2 gene with the highest risk coefficient in PMSRPS was significantly higher in PAAD than in paracancerous tissues, and knockdown of it significantly delayed the progression of PAAD.
Conclusions: PMSRPS can well predict the prognosis, TME and immunotherapy response of patients with PAAD, identify potential drugs, and provide treatment guidance based on individual needs.
1 Introduction
Pancreatic adenocarcinoma (PAAD) is one of the most lethal malignancies, with a 5-year survival rate of only 11% after diagnosis, and is expected to become the second leading cause of cancer death in the United States over the next few decades (1, 2). PAAD is difficult to detect in the early stage, and metastatic spread often occurs when it is detected. It is estimated that only about 15-20% of patients with PAAD currently have the chance to have a radical cure through surgery, and the majority of patients’ prognoses are improved by adjuvant systemic chemotherapy (1). Despite this, patients with PAAD who receive a combination of surgery, chemotherapy, and radiotherapy benefit only in a small percentage of cases (3). In addition, four major immune defects also prevent the majority of patients with PAAD from responding optimally to emerging immunotherapies, and they are the absence of effective intratumoral T cells, heterogeneous dense stroma, immunosuppressive tumor microenvironment, and a lack of tumor-killing T cells (4). In short, the poor efficacy and prognosis of patients with PAAD make it particularly important to construct a new prognostic signature and guide individual treatment plans.
Pyroptosis is a new form of programmed cell death, which manifests as swelling of cells until their membrane ruptures, activating a strong inflammatory response by releasing cellular contents. It is biochemically characterized by inflammasome formation, activation of the caspase family and gasdermin, formation of membrane pores, and release of numerous proinflammatory cytokines, such as IL-1β and IL-18 (5–8). Pyroptosis-induced inflammatory responses can protect the host from microbial infection through classical or non-classical pathways (9). In addition to infectious diseases, more and more studies have recently confirmed the role of pyroptosis in malignant tumors’ s occurrence and development (10, 11). Chen et al. reported that the combination of ruthenium (II) polypyridyl complex Δ-Ru1 and Taxol enhance the anti-cancer effect on Taxol-resistant cancer cells through Caspase-1/GSDMD-mediated pyroptosis (12). Zhang et al. showed that injecting intratumorally with DM-αKG significantly inhibited tumor growth and metastasis through caspase-8-mediated GSDMC-dependent pyroptosis (13). Moreover, commonly used chemotherapeutic drugs, such as cisplatin and paclitaxel, are also effective at inhibiting tumor proliferation and metastasis by inducing pyroptosis (14, 15).
In recent years, gene expression signatures based on pyroptosis-related genes(PRGs) have been reported to predict the prognosis of many cancers, including bladder cancer (16), hepatocellular carcinoma (17), uveal melanoma (18), PAAD (19) among others. However, identification of molecular subtypes on the basis of gene expression seems to be a promising new approach as it helps to rapidly identify cancer features and derive the most appropriate treatment strategies (20). A high quality review also confirmed that the classification of pancreatic tumors based on their molecular characteristics is important to improve the accuracy of clinical treatment decisions (21).
The intervention of pyroptosis identifies a new area of research for the prognosis and treatment of patients with PAAD. Therefore, we established a novel prognostic signature based on the pyroptosis molecular subtype to predict OS and better guide treatment in patients with PAAD.
2 Materials and methods
2.1 Collection of data
The RNA expression data and clinical data of four independent PAAD cohorts(TCGA-PAAD, n = 182 [178 tumors, 4 normal]; GSE102238, n = 100 [50 tumors, 50 normal]; GSE57495, n=63[63 tumors]; ICGC-PACA-CA, n = 262 [262 tumors]) were downloaded from the following public databases: The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/), the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) and the International Cancer Genomics Consortium (ICGC) (https://dcc.icgc.org/). Cell lines gene expression matrices for primary PAAD were obtained from the Cancer Cell Line Encyclopedia(CCLE) dataset (https://portals.broadinstitute.org/ccle/about). TCGA and GSE102238 cohorts are used to build the signature, and the ICGC and GSE57495 cohorts are used to externally verify the signature. To merge TCGA with the GSE102238 cohort, we proceeded as follows. First, we converted the fragments per kilobase of transcript per million value of the RNA expression data of the TCGA cohort to transcripts per kilobase million value. Then, the intersection genes of TCGA and GSE102238 cohorts were extracted, and the corresponding expression data were matched. Next, we removed normal samples of the TCGA and GSE102238 cohorts, converted the TCGA TPM data to log2(TPM+1) data, and normalized the data from the GSE102238 cohort. Finally, we merged the processed TCGA with the GSE102238 cohort into a cohort, and used the ComBat function of the “SVA” package to remove batch effects on the merged data(n=227) (22). To reduce the probability of non-cancer death, we excluded patients with PAAD with a survival time of<30 days.
2.2 Establishment of molecular subtypes based on PRGs expression and screening of differentially expressed genes
52 PRGs (Supplementary Table 1) were from previous researches (23–25). The TCGA and GSE102238 cohorts were combined into one cohort (n=227), and the expression data of 23 PRGs were obtained after matching with 52 PRGs. k-means clustering algorithm was used to obtain the different molecular subtypes associated with PRGs expression. Every molecular subtype was named pyroptosis-related gene cluster (PRGcluster). Single sample gene set enrichment analysis algorithm (ssGSEA) and Normalized enrichment score (NES) were used to quantify the infiltrating abundance of 23 immune cells in different PRGclusters (26). Principal component analysis (PCA) can determine whether the three PRGclusters can be separated. Gene set variation analysis (GSVA) enrichment analysis was used to discover the underlying biological functions between PRGclusters. |log2FC|>0.585 and false discovery rate (FDR)<0.05 were considered as the criteria for screening differentially expressed genes (DEGs) among the three PRGclusters (27).
2.3 Gene enrichment analysis
The Gene Ontology(GO) and Kyoto Encyclopedia of Genes and Genomes(KEGG) gene enrichment analysis was used to explore the relevant cytological functions and pathways of DEGs. q-value<0.05 was the cutoff criterion for determining whether a gene is significantly enriched.
2.4 Establishment and identification of prognostic signature
Univariate Cox regression analysis helped screen prognosis-related differentially expressed genes (PRDEGs). The TCGA+GSE102238 cohort composed of transcriptome data of PRDEGs and the corresponding survival information were randomly divided into the training cohort and testing cohort. The least absolute shrinkage and selection operator (LASSO) Cox regression analysis was performed on PRDEGs in the training cohort to construct a refined prognostic signature (28). Here is the formula for the risk score:
“Coef(Xi)”, “exp(Xi)”, and “n” represent the gene coefficient, expression level, and number of genes, respectively. The median risk score of the training cohort divided patients in the training cohort and validation cohorts into the low- and high-risk groups. Kaplan-Meier survival curves and time-dependent ROC were used to evaluate the prognostic predictive performance of pyroptosis molecular subtypes-related prognostic signature (PMSRPS). To analyze PMSRPS’ independent prognostic performance, univariate and multivariate Cox regression analyses were performed.
2.5 Immune and mutant landscapes between different risk groups
The Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm was used to quantify immune cell infiltration. Based on the Estimation of STromal and Immune cells in Malignant Tumours using the Expression data (ESTIMATE) algorithm (29), we calculated the stromal score, immune score, ESTIMATE score (sum of Stromal and immune scores) for each sample to quantify tumors Microenvironment. Potential immune checkpoint molecules refer to published papers (30, 31) (n=40, Supplementary Table 2). Immunophenoscore (IPS) was downloaded from The Cancer Immune Atlas (TCIA) (https://www.tcia.at/home) to compare responses to immune checkpoint inhibitors(ICIs) across different risk groups. TMB, the number of gene mutations, and the type of gene mutation were derived from the somatic mutation data(n=182) of the TCGA database. The gene copy number (n=185) was downloaded from the specific website(https://xena.ucsc.edu) to observe the copy number variation of PRGs. The RNA stemness scores (RNAss) were downloaded from the Pan-Cancer Atlas Hub (https://pancanatlas.xenahubs.net) (32).
2.6 Chemotherapeutic drug sensitivity and Identification of small-molecule compounds
Half maximal inhibitory concentration (IC50) was using for predicting the sensitivity of chemotherapy drugs in the high- and low-risk groups. |log2FC|>0.585 and FDR<0.05 were the screening criteria for DEGs between high- and low-risk groups. After uploading the up-regulated genome (log2FC>0) and down-regulated genome (log2FC<0) to the L1000FWD website (https://maayanlab.cloud/L1000FWD), a table including drug name, similarity score, and q-value among others will be obtained. Furthermore, 2D and 3D images of small molecule compounds were obtained from the PubChem website (https://pubchem.ncbi.nlm.nih.gov/).
2.7 Cell culture and transfection
The human PAAD cell line CFPAC-1 was derived from Procell (CL-0059, Wuhan, China) and was cultured in a special medium (Procell, CM-0059) in a 37°C, 5% CO2 incubator. CFPAC-1 cells have been identified by short tandem repeat (STR) sequences. The SFTA2-specific siRNA and negative control used in this study were from GenePharma (Shanghai, China). The sequences of the siRNAs are shown in Supplementary Table 3. CFPAC-1 cells were transfected using Lipofectamine 3000 reagent (Invitrogen, Waltham, USA) according to standard guidelines.
2.8 Samples and real-time quantitative PCR
The tumor tissue and paired adjacent tissue of PAAD diagnosed after operation were collected and stored at -80°C. Ethical approval was obtained from the Medical Research Ethics Committee of the Second Affiliated Hospital of Nanchang University. Consistent with a previous study (33), total RNA was extracted from tissues and transfected cells, respectively, and reverse transcription and real-time PCR were performed. The following primers used in the experiments are also shown in Supplementary Table 3.
2.9 In vitro functional experiments
2.9.1 Cell counting kit-8
Cells in each group were seeded into 96-well plates at 5 × 103 cells per well. After the cells were attached, adding 10 μl of CCK-8 reagent (GlpBio, GK10001, USA) to each well at 0, 24, 48 and 72 hours, respectively. After culturing for two hours, the absorbance value of each well was detected using a microplate reader at a wavelength of 450 nm.
2.9.2 5-ethynyl-2’-deoxyuridine assay
Cells in each group were seeded into 96-well plates at 1.5 × 104 cells per well. The next day, EdU staining was performed using an EdU kit (RiboBio, C10310-2, China) according to the manufacturer’s instructions. The percentage of EdU positive cells was calculated by the following formula: EdU positive rate=number of EdU positive cells/number of DAPI positive cells×100%.
2.9.3 Transwell migration assay
Cells in each group were seeded into transwell chambers (Corning, 3422, USA) at 5 × 104 cells per well. After 24 h, cells were fixed with 4% formaldehyde solution for 20 min, and then stained with 0.4% crystal violet solution for 20 min. The cells inside the cell were wiped with a moist cotton swab and allowed to dry, and the migrated cells outside the transwell chambers were counted using an inverted microscope.
2.10 Statistical analysis
This study used R software version 4.1.1 to analyze data. Log-rank test was used for survival analysis. Spearman correlation was used to analyze the correlation of risk scores with immune cell infiltration abundance, TMB, and ssRNA. Wilcoxon rank-sum test and chi-square test were used to judge the significant difference(P<0.05) between two or more groups.
3 Results
3.1 Effects of 52 PRGs on genetic variation, expression level and prognosis of PAAD
Figure 1 shows the overall process of the study. Because of 52 PRGs obtained in this study were derived from ovarian cancer, esophageal adenocarcinoma and glioma studies (23–25), whether they have an effect on PAAD is uncertain. Therefore, we decided to conduct preliminary verification of these PRGs. Somatic mutation and gene expression data have been reported to help infer cancer progression (34). Beside gene mutation, about 95% of patients with PAAD have other genetic changes such as gene amplification and deletion (35). Therefore, we focused on analyzing the genetic variation and expression of PRGs in PAAD. The waterfall plot (Figure 2A) showed the frequency and type of somatic mutation of 52 PRGs in TCGA-PAAD. We found that among 158 samples, there are 90 mutations with a frequency of 56.96%, and TP53 demonstrated the highest mutation frequency, and most samples had the nonsense mutation(a type of point mutation in which a change in a base mutates a codon representing an amino acid into a stop codon, resulting in premature termination of peptide synthesis). The copy number variation (CNV) is an important basis and source of genetic variation (36). Figure 2B revealed that CNV changes occurred in all 52 PRGs. Of these, 27 PRGs(51.92%) were dominated by copy loss, 18 PRGs(34.62%) were dominated by copy gain, and the remaining PRGs had the same copy loss and copy gain frequency. The CNV of these genes on the chromosome is shown in Figure 2C. Next, we compared the expression of 52 PRGs in PAAD and normal pancreatic tissues by combining TCGA and the Genotype-Tissue Expression project (GTEx)(https://xenabrowser.net/datapages/) database. The results showed that 36 PRGs were significantly overexpressed and only 11 underexpressed in PAAD (Figure 2D). We also explored the prognostic value of PRGs matched with survival information and found that the expression of the majority of PRGs (16/23) correlated with the prognosis of patients with PAAD (Supplementary Figure 1). In total, the genetic variation, expression level and prognostic changes presented by 52 PRGs in PAAD indicated that they may play a role in the occurrence and development of PAAD, and it is convincing to use them for further analysis.
Figure 1 Rough flow chart of this study.
Figure 2 The genetic variation landscape of 52 pyroptosis-related genes (PRGs). (A) Somatic mutation frequency and type of PRGs. (B) Copy number variation frequency map of PRGs. (C) Copy number circle plot. (D) Differential expression of PRGs in the tumor and normal samples(Combined TCGA and GTEx databases). ***P< 0.001.
3.2 Identification of pyroptosis molecular subtypes and extraction of DEGs
Prevalent intratumoral and intertumoral heterogeneity has been revealed an important cause of poor prognosis of PAAD. Therefore, in recent years, increasing evidence supports the subtyping of pancreatic tumors based on their molecular characteristics to improve the accuracy of clinical decision-making on treatment (21). Among them, transcriptome subtyping has been used by more and more people because of its unbiased classification that is robust and reproducible. Therefore, we conducted consensus clustering and typing of 23 PRGs that matched transcriptome data and found that the intergroup correlations were the lowest and the intragroup correlations were the highest when clustering variable k = 3, indicating that the joint cohort (n = 227, merged by TCGA and GSE102238) could be well clustered into three subtypes namely PRGclusters A, B, and C according to the expression of the 23 PRGs (Figure 3A). The results of PCA analysis also confirmed the rationality of this classification (Figure 3B). The results of the Kaplan-Meier survival analysis showed that the OS of PRGcluster C was significantly better than that of PRGclusters A and B (Figure 3C). The results of ssGSEA analysis showed that PRGcluster B had the highest infiltration abundance in most immune cells, which not only included anti-tumor immune cells, but also pro-tumor immunosuppressor cells and immune cells with function of anti-tumor and pro-tumor, followed by PRGcluster A and PRGcluster C (P<0.05, Figure 3D). Only the infiltrate abundance of Type.17.T.helper. cell and CD56dim.natural. killer(NK) did not differ significantly among the three PRGclusters. Next, we tried to find the correlation between the prognosis of different PRGclusters and immune cell infiltration and found that PRGcluster C had the least infiltration in some immunosuppressive cells that promoted tumors: MDSC, immature B cell, mast cell, plasmacytoid. dendritic cell, regulatory T cell(Tregs) and Type.2.T.helper.cell, which seem to help understand why PRGcluster C has the best prognostic effect. However, PRGcluster B has a lower prognosis than PRGcluster C despite having the highest invasion abundance among most antitumor immune cells. We speculate that this is the result of the checks and balances and interactions between different functional immune cells in PRGclsuter. Together these results suggest a complex association between immune cell infiltration and prognosis of different pyroptosis molecular subtypes of PAAD. The heatmap presents the clinicopathological features of the three PRGclusters (Figure 3E). In addition, we performed GSVA analysis on three PRGclusters and found that PRGcluster A was highly expressed in apoptosis, linoleic acid metabolism, and P53 signaling pathway; PRGcluster B was highly expressed in infection, immunity, pyroptosis, and various tumors; and PRGcluster C was up-regulated in olfactory transduction and glycine serine and threonine metabolism (Figure 3F). These results indicated that different types of pyroptosis had significantly different effects on the prognosis, immune activity, and biological function of patients with PAAD. To highlight these differences, a differential analysis of the expression levels of three PRGclusters was performed, and a total of 828 DEGs were obtained (Figure 3G).
Figure 3 Molecular subtypes based on the expression of 23 PRGs and identification of DEGs. (A) 227 PAAD samples were divided into three PRGclusters using the Euclidean algorithm. (B) Principal component analysis (PCA) among three PRGclusters. (C) Kaplan-Meier survival analysis among three PRGclusters. (D) Comparison of immune cell infiltration abundance among three PRGclusters through ssGSEA analysis. (E) Heatmap showing clinicopathological information of three PRGclusters. (F) GSVA analyses of any two PRGclusters. (G) The Venn map of 828 differentially expressed genes (DEGs). (H) A sample of 276 DEGs associated with prognosis was divided into three geneClusters using the Euclidean algorithm. (I) Kaplan-Meier survival analysis among three geneClusters. (J) Differential expression of PRGs among three geneClusters. *P< 0.05, **P< 0.01, ***P< 0.001. ns, no statistical difference.
3.3 The potential biological function of DEGs
To better understand the potential function of 828 DEGs extracted, we performed GO and KEGG analyses. The GO enrichment results showed that DEGs were significantly enriched in leukocyte mediated immunity and lymphocyte mediated immunity (biological process), immunoglobulin complex, and external side of the plasma membrane (cellular component), as well as antigen binding and immunoglobulin receptor binding (molecular function) (Figure 4A). The KEGG enrichment results revealed that DEGs prominently enriched in staphylococcus aureus infection and hematopoietic cell lineage among others (Figure 4B). Therefore, DEGs may have biological functions that are relevant to cellular immunity and pyroptosis.
Figure 4 Potential biological functions of 828 DEGs. (A) The results of GO enrichment analysis. (B) The results of KEGG enrichment analysis. q-value<0.05 was the cutoff criterion for significant gene enrichment.
3.4 Verification of stability of pyroptosis molecular subtypes and establishment of PMSRPS
Next, we did a univariate Cox regression analysis for the DEGs combined with survival information and got 276 PRDEGs(Supplementary Table 4). Prior to grouping and modelling PRDEGs, it was realized that PMSRPS were established based on pyroptosis molecular subtypes(i.e. PRGclusters), failure to confirm the stability of subtypes inevitably affects the credibility of the signature, prompting us to perform another consensus clustering and typing of PRDEGs that match transcriptional data (37). The optimal cluster number supported the existence of three distinct and robust geneClusters in patients with PAAD (Figure 3H). Among these three geneClusters, the significant difference in OS was strikingly consistent with the result of PRGclusters (Figure 3I). Also, the infiltration abundance of 23 immune cells was highly in accordance with the differences among the PRGclusters (Figure 3J). More interestingly, we found that the majority of the geneClusters’ patients with PAAD came from the corresponding PRGclusters. In the survival analysis, the proportions of PRGclusters A, B and C in geneClusters A, B and C were 76.69%, 67.12% and 91.67%, respectively. In the immunoinfiltration analysis, the proportions of PRGclusters A, B and C in geneClusters A, B and C was 77.54%, 67.53% and 91.67%, respectively. This may provide a reasonable explanation for the high similarity in the survival and immune cell infiltration’ results between two clusters. All in all, these results suggest that pyroptosis molecular subtypes we identified are robust and reliable.
276 PRDEGs carrying transcriptomic data were subsequently matched with survival information to get a TCGA+GSE102238 cohort(n=218). The cohort was equally divided into a training cohort (n = 109) and a testing cohort (n = 109), and a significant difference was not found in proportion of clinical data and database sources between the two cohorts (Table 1). Next, we performed a LASSO Cox regression analysis on these PRDEGs in the training cohort (Figure 5A) and obtained a PMSRPS with the best fitting effect. The signature was composed of SFTA2, NCAM1, and SPRR1B. According to the risk coefficients of these three genes, we got a formula: Risk score=expression (SFTA2)×0.118+expression(NCAM1)×(-0.256)+expression(SPRR1B)×0.109. The median risk score of the training cohort, 1.053, divided the cohort of patients with PAAD into the high- and low-risk groups. The Venn diagram showed the relationship between clusters and risk scores, as well as survival status (Figure 5B). By comparison, we observed that both PRGcluster A and geneCluster A patients had significantly higher risk scores than other clusters (P<0.001, Figures 5C, D). In addition, the expression of 16 PRGs was also significantly different between the high- and low-risk groups (P<0.05, Figure 5E).
Table 1 Comparison of clinical data between the training and testing cohorts.
Figure 5 The establishment of PMSRPS and its relationship to molecular subtypes and PRGs. (A) Lasso regression analysis. (B) The Venn plot shows the relationship between molecular subtypes and risk score as well as survival state. (C) Distribution of risk scores among three PRGclusters. (D) Differences in risk scores among three geneClusters. (E) Differences in PRGs expression between the high- and low-risk groups. *P< 0.05, **P< 0.01, ***P< 0.001.
3.5 Internal validation of PMSRPS
The Kaplan-Meier survival curve showed that in the training cohort, the OS of patients with PAAD in the high-risk group was worse than that in the lower-risk group (P=0.001, Figure 6A). Also, the survival rate was obviously decreased in patients with PAAD with increasing risk scores (Figure 6A). Time-dependent ROC curves showed that the area under the curves (AUCs) for predicting 1-, 3-, and 5-year survival were 0.732, 0.684, and 0.606, respectively, in the training cohort (Figure 6A). As expected, we observed highly consistent results in the testing and TCGA+GSE102238 cohorts (all P<0.05, Figures 6B, C). This indicated that PMSRPS had a good predictive performance for the prognosis of patients with PAAD.
Figure 6 Internal cohorts validation of the predicted performance of PMSRP. (A–C) Comparison of overall survival (OS) between the high- and low-risk groups, and time-dependent ROC curves of risk scores in the training (A), the testing (B), and the TCGA+GSE102238 cohorts (C). (D) Univariate and multivariate Cox regression analyses for independent prognostic performance assessment of PMSRPS. (E) Relationship between risk score and T stage. (F–K)The differences in OS between the high- and low-risk groups when patients with age<=65 (F) or >65 years (G), the gender of female (H) or male (I), and the tumor T stage of T1-2 (J) or T3-4 (K).
Subsequently, we verified the independence, clinical correlation and applicability of PMSRPS in the TCGA+GSE102238 cohorts. The results showed that risk scores could independently predict poor outcomes in patients with PAAD (P< 0.001, Figure 6D). That is, for all patients with PAAD, the higher the risk score, the worse the prognosis. There was also a significant relationship between risk scores and T stage, and we found that the risk scores of T3-4 were significantly higher than those of T1-2(P=0.029, Figure 6E), which also seemed to help explain the propensity of high risk scores for poor prognosis. In addition, results of applicability analysis showed that the survival rate was lower in the high-risk group than in the low-risk group when patients with age<=65 years (P= 0.001, Figure 6F) or >65 years(P= 0.082, Figure 6G), the gender of female or male (all P< 0.01, Figures 6H, I), and the tumor T stage of T1-2 or T3-4 (all P< 0.05, Figures 6J, K). This showed that PMSRPS could accurately predict and distinguish the prognosis of high- and low-risk groups at different ages, different genders, and different T stages. Taken together, the above results suggested that PMSRPS was an independent predictor of poor prognosis in patients with PAAD, and has general applicability.
3.6 External cohort verification and performance comparison of PMSRPS
To confirm the robustness of PMSRPS, we calculated the risk scores of patients with PAAD in the ICGC and GSE57495 cohorts using the same formula, and divided patients into the high- and low-risk groups based on median risk score(1.053) (Figures 7A, F). Likewise, we observed that the number of survivors decreased with increasing risk scores (Figures 7A, F), and patients in the high-risk group had significantly lower OS than those in the low-risk group (all P< 0.01, Figures 7B, G). Time-dependent ROC curves showed that AUCs for PMSRPS to predict 1, 3 and 5-year survival rates for patients in ICGC cohort were 0.602, 0.610 and 0.738, respectively, while AUCs for predicting 1, 3 and 5-year survival rates for patients in GSE57495 cohort were all >0.7 (Figures 7C, H). Cox regression analysis revealed that risk score were independent prognostic factors in patients with PAAD (all P< 0.01, Figures 7D, E, I, J). Also, we verified the clinical relevance of PMSRPS in the ICGC cohort and found that the expression of the risk gene SFTA2 was higher in patients with PAAD aged 65 years and younger than in patients older than 65 years (Figure 7K). Moreover, the risk gene SPRR1B was highly expressed in the disease state of tumor recurrence/progression (Figure 7L).
Figure 7 External cohort validation of the predicted performance of PMSRPS and performance comparison. (A, F) Risk grouping, survival status, and risk heatmap between the high- and low-risk groups in the ICGC (A) and GSE57495 (F) cohorts. (B, G) Comparison of OS between the high- and low-risk groups in the ICGC (B) and GSE57495 (G) cohorts. (C, H) Time-dependent ROC curves of PMSRPS in the ICGC (C) and GSE57495 (H) cohorts. (D, E, I, J) Results of independent prognostic analysis of PMSRPS in the ICGC (D, E) and GSE57495 (I, J) cohorts. (K, L) Results of external clinical correlation analysis for the high- and low-risk groups in the ICGC cohort. (M–O) Comparison of the performance of PMSRPS with Bai-Sig, Xiao-Sig and Wu-Sig in predicting 1- (M), 3- (N), 5-year (O) OS in patients with PAAD.
Furthermore, to highlight the advantages of PMSRPS, we compared the performance of PMSRPS with three PAAD prognostic signatures recently published. The first is a pyroptosis-related risk signature constructed by Bai et al. (38). The second is a DNA-methylation-driven genes based prognostic signature (GPRC5A, SOWAHC, S100A14, ARNTL2) created by Xiao et al(39). The third is the m6A-related RNA signature(AP005233.2, AC092171.3, AC010175.1, CASC8, TP53TG1, SNAI3.AS1, FLRT1, AC022098.1, DCST1.AS1) constructed by Wu et al (40). The results showed that the AUCs of PMSRPS (0.676, 0.631, 0.600) for predicting the 1-, 3-, and 5-year overall survival of patients with PAAD were higher than those of Bai-Sig (0.609, 0.594, 0.542), Xiao-Sig (0.623, 0.614, 0.557) and Wu-Sig (0.667, 0.630, 0.551) (Figures 7M–O). This reveals that the predictive performance of PMSRPS is better than some existing prognostic signatures.
3.7 Prediction of tumor immune microenvironment and immunotherapeutic response in patients with PAAD by PMSRPS
The TIME of dynamic evolution during tumor progression is an important histological feature of PAAD, which makes us eager to know whether there are TIME differences between high - and low-risk groups with different survival outcomes. We first analyzed the infiltration of immune cells. Results of the this study suggested that the abundance of Tregs, NK cells activated, Macrophages M0, and Mast cells activated infiltration was distinctly higher in the high-risk group than in the low-risk group (all P< 0.01, Figure 8A), and with increasing risk scores, their infiltration abundance increased prominently (all P< 0.01, Figure 8B). However, the abundance of T cells CD8, T cells CD4 memory resting, B cells naïve, Monocytes, and Mast cells resting infiltration in the high-risk group was significantly lower than that in the low-risk group (all P< 0.05, Figure 8A), and they were less abundant with greater risk scores (all P< 0.01, Figure 8B). In addition to infiltrating immune cells, infiltrating stromal cells are also vital components of tumor TIME, which together perturb the tumor signal and play an important role in cancer biology (29). Therefore, we used the ESTIMATE method to assess the distribution of these components between the high- and low-risk groups, and found that the stromal score, immune score, and estimated score were significantly lower in the high-risk group than in the low-risk group (all P<0.05, Figure 8C), which meant PMSRPS could accurately predict and distinguish TIME in the different risk groups, and that the high-risk group had higher tumor purity and more typical immunosuppressive “cold” tumor microenvironment. Furthermore, researchers found the expression of immune checkpoint molecules on a significant percentage of tumor cells, and they have been shown to promote epithelial-mesenchymal transformation, resistance to apoptosis and antitumor drugs, and propensity to spread and metastasize (41). In this study, we observed that about half (19/40) of the immune checkpoint molecules had significant expression differences between the two groups. Among them, a total of 12 genes, including PDCD1, BTLA, and CD28, were significantly highly expressed in the low-risk group, while 7 genes, including TNFRSF14, LGALS9, and CD276, were significantly highly expressed in the high-risk group (all P< 0.05, Figure 8D). To overcome the immunosuppressive power of checkpoint molecules, a new immunotherapy called ICIs works by blocking these immunosuppressive molecules and reactivating effector T cells to specifically kill tumor cells. Currently, the main ICIs used in clinical practice are CTLA4 inhibitors, PD-1 inhibitors and PD-L1 inhibitors. We evaluated the response of patients with PAAD in different risk groups to ICIs by IPS scores and found that the IPS (PD-1/PD-L1/PD-L2(-) and CTLA-4(-)) scores of patients in the high-risk group were significantly higher in the lower-risk group(P< 0.01, Figure 8E), rather than IPS-PD1/PDL1/PDL2 blocker score, IPS-CTLA4 blocker score, and IPS-CTLA4 and PD1/PDL1/PDL2 blocker score (Supplementary Figure 2). This suggested that patients with PAAD in the high-risk group might be more suitable for novel ICIs therapy rather than PD1/CTLA4 immunotherapy. Besides PMSRPS, we were surprised to find that the model gene had good predictive performance for the existing immunotherapy cohort via a web server for Comprehensive Analysis on Multi-Omics of Immunotherapy in Pan-cancerimmune checkpoint inhibitors (CAMOIP) (http://www.camoip.net/). The results showed that the highly expressed prognostic risk gene SPRR1B had significantly worse OS (HR: 5.5, P=0.001, Figure 8F) in the melanoma immunotherapy cohort of Hugo. et al. (42). And the highly expressed prognostic protective gene NCAM1 had significantly better OS (HR: 0.35, P=0.022, Figure 8G) in the melanoma immunotherapy cohort of Auslander et al. (43). Taken together, the above results manifested that PMSRPS might be a reliable tool for predicting immune activity and immunotherapy response in patients with PAAD.
Figure 8 Prediction of TIME and immunotherapy response by PMSRPS. (A) Comparison of immune cell infiltration in the high- and low-risk groups. (B) Correlation of risk score with infiltrating abundance of 10 immune cells. (C) Comparison of tumor microenvironment scores in the high- and low-risk groups. (D) Differential expression of immune checkpoint molecules in the high- and low-risk groups. (E) Comparison of IPS-PD1/PDL1/PDL2(-) and CTLA4(-) scores between the high- and low-risk groups. (F, G) The predictive role of SPRR1B (F) and NCAM1 (G) genes in PMSRPS on existing immunotherapy cohorts was explored via a web server for Comprehensive Analysis on Multi-Omics of Immunotherapy in Pan-cancerimmune checkpoint inhibitors (CAMOIP). *P< 0.05, **P< 0.01, ***P< 0.001.
3.8 Prediction of PMSRPS on other measures related to immunotherapy response
Although ICIs is considered a breakthrough in cancer treatment, the benefit to patients is limited, so it makes sense to understand the underlying indicators that influence response to ICIs treatment. Cancer stemness was reported to be significantly positively associated with ICIs resistance in cancer (44). TMB predicts clinical response to ICIs in a variety of tumors and is associated with improved survival after receiving ICIs (45). Genomic stability was the third proven predictor of ICIs response (46). Therefore, we explored these indicators in the high- and low-risk groups with different immunotherapy responses. As shown in Figures 9A, B, the TMB and RNAss of the high-risk group were visibly higher than those of the low-risk group (all P< 0.05), and TMB and RNAss had distinctly increased with the increase of the risk score (all P ≤ 0.01, Figures 9C, D). In addition, the genomic mutation rate in the high-risk group was significantly higher than that in the low-risk group (90.48% [76/84] vs 59.26% [32/54], P<0.001, Figures 9E, F). Among them, gene mutation rates of the top two genes, KRAS and TP53, increased more significantly from the low- to the high-risk groups (KRAS: 71% [60/84] vs. 30% [16/54], P<0.001; TP53: 67% [56/84] vs 37% [20/54], P=0.001). These results suggest that high TMB, cancer stemness, and genomic instability are characteristic of high-risk patients compared with low-risk patients, which can help evaluate and predict immunotherapy effects in patients with different risk groups from multiple perspectives.
Figure 9 Prediction of PMSRPS on other measures related to immunotherapy response. (A, B) Comparison of tumor mutational burden(TMB) (A) and RNA stemness scores (RNAss) (B) between the high- and low-risk groups. (C, D) Correlation of risk score with TMB (C) and RNAss (D). (E, F) Changes in gene mutation frequency and mutation type from the low-risk group (E) to the high-risk group (F).
3.9 Sensitivity analysis of chemotherapeutic drugs and screening of small molecule drugs
Immunotherapy has important therapeutic value for some malignancies, but in most phase I and II trials, It has failed to show good clinical efficacy in patients with PAAD when used alone, unless combined with chemotherapy (47). This suggests that it is necessary to simultaneously screen for more effective chemotherapy drugs. The results of chemotherapeutic drug sensitivity analysis in this study revealed that the high-risk group was more sensitive to commonly used chemotherapy drugs for PAAD (Paclitaxel, Doxorubicin, Docetaxel among others) than the low-risk group (all P< 0.05, Figure 10). Moreover, PMSRPS can predict chemotherapy responses in different risk populations to other cancers (Bicalutamide and Lenalidomide among others). (all P< 0.001, Figure 10).
Figure 10 Comparison of sensitivities to commonly used chemotherapeutics for PAAD and other cancers in the high- and low-risk groups.
In addition, to expand and search for potential small-molecule drugs for PAAD, we extracted DEGs between the high- and low-risk groups. Next, the up-regulated and down-regulated DEGs were uploaded separately to the L1000FWD website and matched with small molecule therapies. This study presented the 10 most important potential small molecule drugs, followed by similarity scores, q-values, and mechanism of action(Table 2). Among them, the top three drugs with negative similarity scores were tamoxifen, ZM-241385, and BRD-A24021119, which were expected to negatively regulate the expression of DEGs. 2D and 3D images of these three drugs could be seen in Supplementary Figure 3. These potential small-molecule drugs might reverse the high expression of high-risk group genes and guide the development of PAAD-targeted drugs. However, the exact efficacy of these drugs needs to be further confirmed in future prospective studies.
Table 2 Potential small molecule drugs screened for PAAD treatment.
3.10 Expression level and functional verification of SFTA2 gene
Finally, to further verify the reliability of the PMSRPS, we performed expression and in vitro functional experiments on the SFTA2 gene with the highest correlation coefficient. As shown in Figure 11A, we found that SFTA2 was highly expressed in PAAD tissues compared with the corresponding paracancerous tissues (P<0.01). Further, we selected PAAD CFPAC-1 cell lines with high expression through CCLE database for SFTA2 gene interference (Figure 11B). Figure 11C shows that si-SFTA2#1 and si-SFTA2#2 had more obvious inhibitory effects, which was used for further functional experiments. Through CCK-8 experiments, we found that interfering with the expression of SFTA2 inhibited the cell viability of PAAD cells (Figure 11D). By EdU staining, we found that compared with the control group, the proliferation rate of the SFTA2-inhibited group was significantly reduced (Figure 11E). Through transwell migration experiments, we found that interfering with SFTA2 expression would inhibit the migration of PAAD cells (Figure 11F).
Figure 11 Expression and function of the highest-risk SFTA2 gene in PMSRPS. (A) qRT-PCR analysis of SFTA2 mRNA levels in randomly selected 16 pairs of PAAD tissues and corresponding adjacent tissues. (B) The Cancer Cell Lines Encyclopedia (CCLE) database was searched for PAAD CFPAC-1 cell lines suitable for SFTA2 intervention. (C) Validation of knockdown efficiency of SFTA2 interfering fragments in CFPAC-1 cells by qRT-PCR. (D) CCK-8 assay was used to assess the effect of SFTA2 silencing on CFPAC-1 cell viability. (E) EdU staining was used to assess the effect of SFTA2 silencing on the proliferative capacity of CFPAC-1 cells. (F) Transwell experiments were used to evaluate the effect of SFTA2 silencing on the migration ability of CFPAC-1 cells. *P< 0.05, **P< 0.01, ***P< 0.001.
4 Discussion
The overall prognosis for patients with PAAD has long been disappointing (48). Authoritative cancer statistics in 2022 showed that PAAD has become the third leading cause of cancer-related death in both men and women (2). The poor prognosis of patients with PAAD is associated with the limited therapeutic response, and one possible reason for the poor therapeutic response is the ability of PAAD cells to avoid induction of death (49). Cell death is an important physiological process for maintaining tissue homeostasis, which can be divided into accidental cell death and regulated cell death (RCD) according to its occurrence rate and potential control (49, 50). RCD with clear mechanism of effect can be further subclassified into apoptotic and non-apoptotic subcategories (49). Since apoptosis is strongly resisted by cancer cells (51), targeting non-apoptotic cell death is considered a more promising therapeutic approach (49). Pyroptosis is a non-apoptotic death type characterized by the release of a large number of inflammatory factors, and its induced activation can produce powerful antitumor activity (52, 53). This is illustrated by the fact that wang et al. revealed that pyroptosis of less than 15% of tumour cells was sufficient to clear the entire 4T1 mammary tumour graftI by building a bioorthogonal chemical system (54). Moreover, Zhang et al’ s study showed that GSDME expression could not only enhance the phagocytosis of tumor-related macrophages, but also strengthen the number and function of NK cells and CD8+T cells infiltrated by tumors (55). In studies related to PAAD, Cui et al. demonstrated that MST1 inhibits pancreatic cancer progression via ROS-induced pyroptosis (56). Peng et al. showed that ICy OH produced by ICy Q could damage mitochondrial membranes, induce intracellular inflammatory responses, and selectively induce pancreatic cancer cell death via the cell pyroptosis pathway (a series of enzymatic reactions leading to the production of fragments of pyroptosis protein GSDMD-N) (57). Considering the promise of immunotherapy, researchers recently tested the immunotherapy effect of pyroptosis using the membrane targeted photosensitizer TBD-3C. The results showed that pyroptosis-aroused immunological responses could transfer the immunosuppressive “cold” tumor microenvironment(TME) into an immunogenic “hot” TME, which not only inhibited the growth of primary PAAD, but also attacked distant tumors (58). In brief, a new strategy to induce pyroptosis may provide more effective cancer treatment options (59). However, estimates of disease status relying solely on changes in gene expression are unstable (60). As a result, we sought to establish a prognostic signature through the analysis of pyroptosis molecular subtypes to aid in the accurate prediction of PAAD prognosis and guidance for individualized treatment.
We performed consensus clustering and typing of patients with PAAD based on the expression levels of 23 PRGs and obtained 828 DEGs. Through GO and KEGG enrichment analysis, we found that the potential biological functions of these DEGs were related to the occurrence of cellular immunity and pyroptosis. It has been reported that immune cells macrophages and neutrophils can interact to drive the pathogenesis of PAAD (61). Among them, M2 and a small fraction of M1 cells are not only unable to phagocytose tumor cells, but they also migrate to other tissues and organs without being killed(62). The immature c-Kit+ neutrophil subsets can generate ROS by participating in oxidative mitochondrial metabolism, thereby inhibiting the immune function of CD4+ T cells and promoting tumor progression (63). Furthermore, overactivation of the NLRP3 inflammasome may promote the development of hematopoietic malignancies (64).
Subsequently, we refined Lasso Cox regression analysis on the extracted 276 prognostic DEGs in the training cohort to obtain PMSRPS. Internal and external validation results showed that PMSRPS was an independent and effective prognostic tool for PAAD. Compared with some recently published PAAD prognostic signatures (38–40), PMSRPS had better predictive power in predicting 1-, 3-, and 5-year survival in patients with PAAD. PMSRPS consists of two risk prognostic genes (SFTA2 and SPRR1B) and one protective prognostic gene (NCAM1). SFTA2 is the gene with the highest risk factor in PMSRPS, and although it has been reported to be associated with poor prognosis of PDAC, experimental validation evidence is lacking (65). Therefore, this study verified for the first time that the expression of SFTA2 in PAAD tissues was higher than that in paired paracancer tissues in our own samples. It was further confirmed that the vitality, growth rate and migration ability of PAAD cells were significantly decreased after the SFTA2 expression level was knocked down. In addition, SPRR1B was a prognostic or diagnostic biomarker for various malignancies such as PAAD, lung adenocarcinoma, metastatic cutaneous melanoma, and oral squamous cell carcinoma among others (66–69). NCAM1 was confirmed to be up-regulated mediated by miR-141-3p and inhibited ameloblastoma cell migration (70).
Then, we used PMSRPS to predict the TIME of different risk groups and found that the high-risk group presented a typical tumor immunosuppressive microenvironment. This was reflected in less T cells CD8 and T cells CD4 memory resting as well as more Tregs and macrophage M0 in the high-risk group. Studies have shown that higher CD8 T cell infiltration in PDAC is strongly associated with long-term survival (71). In contrast, Tregs are important factors in maintaining immune tolerance, not only suppressing effector cells within the tumor but also restricting antitumor immune responses by interacting with stroma, vasculature, and lymphatic vessels (72). Also, M0 macrophages can differentiate into M2 macrophages in the presence of M-CSF, IL-4, or IL-10 and promote immune escape through the high expression of PD-L1, IL-10, or TGFβ (73). In addition, mast cells are key regulators of inflammation and immunosuppression (74). Besides suppressing anti-tumor immunity by releasing anti-inflammatory cytokines such as IL-10 and TGFβ, they can provide oxygen for tumor growth by regulating angiogenesis (75, 76). Excitingly, in this study, PMSRPS was confirmed to predict the expression of multiple immune checkpoint molecules in different risk groups and showed that the high-risk group had higher IPS (PD-1/PD-1/PD-1/PD-1/PD-1) and CTLA-4(-)) scores compared with the low-risk group. This suggests that PMSRPS may be a potential biomarker for predicting response to novel immunotherapies in PAAD.
Cancer stem cells with self-renewal ability can induce tumor metastasis and recurrence and are contribute to drug resistance (77, 78). Tumor stemness leading to immunotherapy resistance can be attributed to its restriction of the antitumor immune response by inhibiting type I IFN signaling (79). The results of this study showed that the RNAss of patients in the high-risk group were notably higher than those in the low-risk group, and the RNAss were distinctly positively correlated with the risk score, indicating that the high-risk score might lead to enhanced tumor stem cell characteristics (26), which may explain the poorer prognosis in the high-risk group. Interestingly, TMB was also elevated in the high-risk group and positively correlated with tumor stemness, whereas high TMB was associated with better treatment response and better prognosis. One possible reason for our analysis is that as TMB increases, so does the neoantigens, one or more of which are more likely to be immunogenic and trigger a T cell response that enhances the antitumor response(80). In summary, we can take a more comprehensive view of the immunotherapeutic effect of PAAD.
Finally, we screened several potential small-molecule drugs for patients with PAAD, with tamoxifen in the top spot. As a selective estrogen receptor modulator (SERM), tamoxifen is clinically used for hormone receptor-positive breast cancer. However, in recent years, tamoxifen may inhibit peritoneal mesothelium-mesenchymal transition, as well as peritoneal mesenchymal cell migration, interstitial fibrosis and neovascularization due to its ability to quiescent the peritoneal stroma (81). This may aid in the immunotherapy of PAADs with desmoplastic stroma. Another SERM, bazedoxifene, was shown to have the function of inhibiting STAT3 phosphorylation and STAT3 DNA binding, inducing apoptosis, and suppressing tumor growth in PAAD cells with persistent IL6/GP130/STAT3 signaling (82). In addition, a non-selective adenosine receptor antagonist caffeine, and its analog CGS 15943 has been reported to block the proliferation in HCC and PDAC cell lines by inhibiting the PI3K/Akt pathway (83). However, the exact function of above small molecule drugs remains unclear and needs to be further confirmed by future prospective studies.
So far, pyroptosis-related prognostic models of PAAD have mushroomed. Similar to these studies, our study also carried out methods such as survival analysis, ROC and external cohort validation of PMSRPS, but it still has its advantages. First, after knowing the importance of subtyping based on molecular characteristics of cancer for improving clinical treatment decision-making, we did not directly use PRGs to establish a prognostic signature as many previous studies had done (84–88), but instead identified three stable pyroptosis molecular subtypes based on PRGs expression levels. By analysis, we found that different PRGclusters had different survival outcomes and immune cell infiltration manifestations. To highlight this difference, we extracted the DEGs between PRGclusters to establish a PMSRPS, which made the modeling process more rigorous. Secondly, there are only three signature genes in this study, which is significantly less than other studies (mostly between 5 and 8), which helps simplify the calculation of risk scores and the risk grouping process for patients with PAAD, thus reducing cost estimates and enabling faster prognosis assessment. Finally, we screened new potential small-molecule drugs for patients with PAAD based on DEGs in the high- and low-risk groups, which provides a reference for PAAD to develop new treatment options. For all this, there are still some shortcomings in this study. One is that despite the combination of the TCGA and GSE102238 cohorts in this study, the PAAD sample size used to construct the PMSRPS is still small, mainly due to the limited number of patients with PAAD with complete clinical information we obtained from current public databases. Another shortcoming is that the transcriptional data for the PRGs included in the merged cohort (TCGA+GSE102238) were incomplete. This situation prevents us from comparing the prediction performance of PMSRPS with other existing PRG models in the TCGA+GSE102238 cohort. In short, a large prospective multicenter study is needed in the future to further verify the exact effect of PMSRPS on the survival and treatment of patients with PAAD.
5 Conclusion
In conclusion, we identified a novel signature, PMSRPS for patients with PAAD. This signature is a good predictor of prognosis, immune microenvironment, immunotherapy effect, genomic instability and tumor stemness in patients with PAAD.
Data availability statement
Publicly available datasets were analyzed in this study. This data can be found here: The TCGA database can be found here: https://portal.gdc.cancer.gov/, the GEO database https://www.ncbi.nlm.nih.gov/geo/, the ICGC database https://dcc.icgc.org/, the GTEx database https://xenabrowser.net/datapages/, and the original data for experimental validation are available through corresponding authors.
Ethics statement
The studies involving human participants were reviewed and approved by The Second Affiliated Hospital of Nanchang University Medical Research Ethics Committee. The patients/participants provided their written informed consent to participate in this study.
Author contributions
QH conducted the formal analysis and wrote the original manuscript. XP and JZ collected the samples and conducted the experiments. JX took the raw data from the public database and did the collation. QH and JZ completed the drawing. QL participated in the major revision of the manuscript. HJ guided design, reviewed the manuscript and provided funding acquisition. All authors contributed to the article and approved the submitted version.
Funding
This research was funded by Natural Science Foundation of China (Number: 82103094).
Acknowledgments
We thank GEO, TCGA, ICGC and GTEx databases for providing invaluable data for statistical analyses.
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.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1111494/full#supplementary-material
References
1. Mizrahi JD, Surana R, Valle JW, Shroff RT. Pancreatic cancer. Lancet (2020) 395(10242):2008–20. doi: 10.1016/S0140-6736(20)30974-0
2. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics 2022. CA Cancer J Clin (2022) 72(1):7–33. doi: 10.3322/caac.21708
3. Narayanan S, Vicent S, Ponz-Sarvise M. PDAC as an immune evasive disease: Can 3D model systems aid to tackle this clinical problem? Front Cell Dev Biol (2021) 9:787249. doi: 10.3389/fcell.2021.787249
4. Upadhrasta S, Zheng L. Strategies in developing immunotherapy for pancreatic cancer: Recognizing and correcting multiple immune "Defects" in the tumor microenvironment. J Clin Med (2019) 8(9):1472. doi: 10.3390/jcm8091472
5. Aglietti RA, Estevez A, Gupta A, Ramirez MG, Liu PS, Kayagaki N, et al. GsdmD p30 elicited by caspase-11 during pyroptosis forms pores in membranes. Proc Natl Acad Sci U.S.A. (2016) 113(28):7858–63. doi: 10.1073/pnas.1607769113
6. Shi J, Gao W, Shao F. Pyroptosis: Gasdermin-mediated programmed necrotic cell death. Trends Biochem Sci (2017) 42(4):245–54. doi: 10.1016/j.tibs.2016.10.004
7. Christgen S, Place DE, Kanneganti TD. Toward targeting inflammasomes: insights into their regulation and activation. Cell Res (2020) 30(4):315–27. doi: 10.1038/s41422-020-0295-8
8. Christgen S, Tweedell RE, Kanneganti TD. Programming inflammatory cell death for therapy. Pharmacol Ther (2022) 232:108010. doi: 10.1016/j.pharmthera.2021.108010
9. Burdette BE, Esparza AN, Zhu H, Wang S. Gasdermin d in pyroptosis. Acta Pharm Sin B (2021) 11(9):2768–82. doi: 10.1016/j.apsb.2021.02.006
10. Hou J, Hsu JM, Hung MC. Molecular mechanisms and functions of pyroptosis in inflammation and antitumor immunity. Mol Cell (2021) 81(22):4579–90. doi: 10.1016/j.molcel.2021.09.003
11. Li T, Zheng G, Li B, Tang L. Pyroptosis: A promising therapeutic target for noninfectious diseases. Cell Prolif (2021) 54(11):e13137. doi: 10.1111/cpr.13137
12. Chen D, Guo S, Tang X, Rong Y, Bo H, Shen H, et al. Combination of ruthenium (II) polypyridyl complex delta-Ru1 and taxol enhances the anti-cancer effect on taxol-resistant cancer cells through caspase-1/GSDMD-mediated pyroptosis. J Inorg Biochem (2022) 230:111749. doi: 10.1016/j.jinorgbio.2022.111749
13. Zhang JY, Zhou B, Sun RY, Ai YL, Cheng K, Li FN, et al. The metabolite alpha-KG induces GSDMC-dependent pyroptosis through death receptor 6-activated caspase-8. Cell Res (2021) 31(9):980–97. doi: 10.1038/s41422-021-00506-9
14. Wang Y, Gao W, Shi X, Ding J, Liu W, He H, et al. Chemotherapy drugs induce pyroptosis through caspase-3 cleavage of a gasdermin. Nature (2017) 547(7661):99–103. doi: 10.1038/nature22393
15. Zhang CC, Li CG, Wang YF, Xu LH, He XH, Zeng QZ, et al. Chemotherapeutic paclitaxel and cisplatin differentially induce pyroptosis in A549 lung cancer cells via caspase-3/GSDME activation. Apoptosis (2019) 24(3-4):312–25. doi: 10.1007/s10495-019-01515-1
16. Lu H, Wu J, Liang L, Wang X, Cai H. Identifying a novel defined pyroptosis-associated long noncoding RNA signature contributes to predicting prognosis and tumor microenvironment of bladder cancer. Front Immunol (2022) 13:803355. doi: 10.3389/fimmu.2022.803355
17. Deng M, Sun S, Zhao R, Guan R, Zhang Z, Li S, et al. The pyroptosis-related gene signature predicts prognosis and indicates immune activity in hepatocellular carcinoma. Mol Med (2022) 28(1):16. doi: 10.1186/s10020-022-00445-0
18. Zhang F, Deng Y, Wang D, Wang S. Construction and validation of a pyroptosis-related gene signature associated with the tumor microenvironment in uveal melanoma. Sci Rep (2022) 12(1):1640. doi: 10.1038/s41598-022-05599-9
19. Chen Y, Liu Y, Wang M. Identification of a pyroptosis-related gene signature and effect of silencing the CHMP4C and CASP4 in pancreatic adenocarcinoma. Int J Gen Med (2022) 15:3199–213. doi: 10.2147/IJGM.S353849
20. Shao W, Yang Z, Fu Y, Zheng L, Liu F, Chai L, et al. The pyroptosis-related signature predicts prognosis and indicates immune microenvironment infiltration in gastric cancer. Front Cell Dev Biol (2021) 9:676485. doi: 10.3389/fcell.2021.676485
21. Huang X, Zhang G, Liang T. Subtyping for pancreatic cancer precision therapy. Trends Pharmacol Sci (2022) 43(6):482–94. doi: 10.1016/j.tips.2022.03.005
22. Huo J, Wu L, Zang Y. Development and validation of a CTNNB1-associated metabolic prognostic model for hepatocellular carcinoma. J Cell Mol Med (2021) 25(2):1151–65. doi: 10.1111/jcmm.16181
23. Ye Y, Dai Q, Qi H. A novel defined pyroptosis-related gene signature for predicting the prognosis of ovarian cancer. Cell Death Discovery (2021) 7(1):71. doi: 10.1038/s41420-021-00451-x
24. Zeng R, Huang S, Qiu X, Zhuo Z, Wu H, Jiang L, et al. Predicting the prognosis of esophageal adenocarcinoma by a pyroptosis-related gene signature. Front Pharmacol (2021) 12:767187. doi: 10.3389/fphar.2021.767187
25. Zhang M, Cheng Y, Xue Z, Sun Q, Zhang J. A novel pyroptosis-related gene signature predicts the prognosis of glioma through immune infiltration. BMC Cancer (2021) 21(1):1311. doi: 10.1186/s12885-021-09046-2
26. Huo J, Cai J, Guan G, Liu H, Wu L. A ferroptosis and pyroptosis molecular subtype-related signature applicable for prognosis and immune microenvironment estimation in hepatocellular carcinoma. Front Cell Dev Biol (2021) 9:761839. doi: 10.3389/fcell.2021.761839
27. Xia C, Chen L, Sun W, Yan R, Xia M, Wang Y, et al. Total saponins from Paris forrestii (Takht) h. li. show the anticancer and RNA expression regulating effects on prostate cancer cells. BioMed Pharmacother (2020) 121:109674. doi: 10.1016/j.biopha.2019.109674
28. Motamedi F, Perez-Sanchez H, Mehridehnavi A, Fassihi A, Ghasemi F. Accelerating big data analysis through LASSO-random forest algorithm in QSAR studies. Bioinformatics (2022) 38(2):469–75. doi: 10.1093/bioinformatics/btab659
29. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612
30. Campbell KS. Mystery checkpoint revealed: KIR3DL3 finally found a ligand in HHLA2. Cancer Immunol Res (2021) 9(2):128. doi: 10.1158/2326-6066.CIR-20-0996
31. Lu L, Hu Y, Wang C, Jiang F, Wu C. Methylation and expression of the exercise-related TLR1 gene is associated with low grade glioma prognosis and outcome. Front Mol Biosci (2021) 8:747933. doi: 10.3389/fmolb.2021.747933
32. Li N, Zhan X. Integrated genomic analysis of proteasome alterations across 11,057 patients with 33 cancer types: clinically relevant outcomes in framework of 3P medicine. EPMA J (2021) 12(4):605–27. doi: 10.1007/s13167-021-00256-z
33. Zhu J, Huang Q, Peng X, Luo C, Liu S, Liu Z, et al. Identification of LncRNA prognostic signature associated with genomic instability in pancreatic adenocarcinoma. Front Oncol (2022) 12:799475. doi: 10.3389/fonc.2022.799475
34. Fleck JL, Pavel AB, Cassandras CG. Integrating mutation and gene expression cross-sectional data to infer cancer progression. BMC Syst Biol (2016) 10:12. doi: 10.1186/s12918-016-0255-6
35. Yu T, Tan H, Liu C, Nie W, Wang Y, Zhou K, et al. Integratively genomic analysis reveals the prognostic and immunological characteristics of pyroptosis and ferroptosis in pancreatic cancer for precision immunotherapy. Front Cell Dev Biol (2022) 10:826879. doi: 10.3389/fcell.2022.826879
36. Iafrate AJ, Feuk L, Rivera MN, Listewnik ML, Donahoe PK, Qi Y, et al. Detection of large-scale variation in the human genome. Nat Genet (2004) 36(9):949–51. doi: 10.1038/ng1416
37. Li L, Gao H, Wang D, Jiang H, Wang H, Yu J, et al. Metabolism-relevant molecular classification identifies tumor immune microenvironment characterization and immunotherapeutic effect in cervical cancer. Front Mol Biosci (2021) 8:624951. doi: 10.3389/fmolb.2021.624951
38. Bai Z, Xu F, Feng X, Wu Y, Lv J, Shi Y, et al. Pyroptosis regulators exert crucial functions in prognosis, progression and immune microenvironment of pancreatic adenocarcinoma: a bioinformatic and in vitro research. Bioengineered (2022) 13(1):1717–35. doi: 10.1080/21655979.2021.2019873
39. Xiao M, Liang X, Yan Z, Chen J, Zhu Y, Xie Y, et al. A DNA-Methylation-Driven genes based prognostic signature reveals immune microenvironment in pancreatic cancer. Front Immunol (2022) 13:803962. doi: 10.3389/fimmu.2022.803962
40. Wu Q, Chen L, Miao D, Jin Y, Zhu Z. Prognostic signature based on m6A-related lncRNAs to predict overall survival in pancreatic ductal adenocarcinoma. Sci Rep (2022) 12(1):3079. doi: 10.1038/s41598-022-07112-8
41. Marcucci F, Rumio C, Corti A. Tumor cell-associated immune checkpoint molecules - drivers of malignancy and stemness. Biochim Biophys Acta Rev Cancer (2017) 1868(2):571–83. doi: 10.1016/j.bbcan.2017.10.006
42. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, et al. Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell (2017) 168(3):542. doi: 10.1016/j.cell.2017.01.010
43. Auslander N, Zhang G, Lee JS, Frederick DT, Miao B, Moll T, et al. Robust prediction of response to immune checkpoint blockade therapy in metastatic melanoma. Nat Med (2018) 24(10):1545–9. doi: 10.1038/s41591-018-0157-9
44. Zhang Z, Wang ZX, Chen YX, Wu HX, Yin L, Zhao Q, et al. Integrated analysis of single-cell and bulk RNA sequencing data reveals a pan-cancer stemness signature predicting immunotherapy response. Genome Med (2022) 14(1):45. doi: 10.1186/s13073-022-01050-w
45. Samstein RM, Lee CH, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet (2019) 51(2):202–6. doi: 10.1038/s41588-018-0312-8
46. Fuca G, Corti F, Ambrosini M, Intini R, Salati M, Fenocchio E, et al. Prognostic impact of early tumor shrinkage and depth of response in patients with microsatellite instability-high metastatic colorectal cancer receiving immune checkpoint inhibitors. J Immunother Cancer (2021) 9(4):e002501. doi: 10.1136/jitc-2021-002501
47. Schizas D, Charalampakis N, Kole C, Economopoulou P, Koustas E, Gkotsis E, et al. Immunotherapy for pancreatic cancer: A 2020 update. Cancer Treat Rev (2020) 86:102016. doi: 10.1016/j.ctrv.2020.102016
48. Kleeff J, Korc M, Apte M, La Vecchia C, Johnson CD, Biankin AV, et al. Pancreatic cancer. Nat Rev Dis Primers (2016) 2:16022. doi: 10.1038/nrdp.2016.22
49. Chen X, Zeh HJ, Kang R, Kroemer G, Tang D. Cell death in pancreatic cancer: from pathogenesis to therapy. Nat Rev Gastroenterol Hepatol (2021) 18(11):804–23. doi: 10.1038/s41575-021-00486-6
50. Galluzzi L, Vitale I, Aaronson SA, Abrams JM, Adam D, Agostinis P, et al. Molecular mechanisms of cell death: recommendations of the nomenclature committee on cell death 2018. Cell Death Differ (2018) 25(3):486–541. doi: 10.1038/s41418-017-0012-4
51. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell (2011) 144(5):646–74. doi: 10.1016/j.cell.2011.02.013
52. Xia X, Wang X, Cheng Z, Qin W, Lei L, Jiang J, et al. The role of pyroptosis in cancer: pro-cancer or pro-"host"? Cell Death Dis (2019) 10(9):650. doi: 10.1038/s41419-019-1883-8
53. Ding B, Sheng J, Zheng P, Li C, Li D, Cheng Z, et al. Biodegradable upconversion nanoparticles induce pyroptosis for cancer immunotherapy. Nano Lett (2021) 21(19):8281–9. doi: 10.1021/acs.nanolett.1c02790
54. Wang Q, Wang Y, Ding J, Wang C, Zhou X, Gao W, et al. A bioorthogonal system reveals antitumour immune function of pyroptosis. Nature (2020) 579(7799):421–6. doi: 10.1038/s41586-020-2079-1
55. Zhang Z, Zhang Y, Xia S, Kong Q, Li S, Liu X, et al. Gasdermin e suppresses tumour growth by activating anti-tumour immunity. Nature (2020) 579(7799):415–20. doi: 10.1038/s41586-020-2071-9
56. Cui J, Zhou Z, Yang H, Jiao F, Li N, Gao Y, et al. MST1 suppresses pancreatic cancer progression via ROS-induced pyroptosis. Mol Cancer Res (2019) 17(6):1316–25. doi: 10.1158/1541-7786.MCR-18-0910
57. Peng B, Chen G, Li Y, Zhang H, Shen J, Hou JT, et al. NQO-1 enzyme-activated NIR theranostic agent for pancreatic cancer. Anal Chem (2022) 94(32):11159–67. doi: 10.1021/acs.analchem.2c01189
58. Wang M, Wu M, Liu X, Shao S, Huang J, Liu B, et al. Pyroptosis remodeling tumor microenvironment to enhance pancreatic cancer immunotherapy driven by membrane anchoring photosensitizer. Adv Sci (Weinh) (2022) 9(29):e2202914. doi: 10.1002/advs.202202914
59. Nagarajan K, Soundarapandian K, Thorne RF, Li D, Li D. Activation of pyroptotic cell death pathways in cancer: An alternative therapeutic approach. Transl Oncol (2019) 12(7):925–31. doi: 10.1016/j.tranon.2019.04.010
60. Wei R, Zhang H, Cao J, Qin D, Li S, Deng W. Sample-specific perturbation of gene interactions identifies pancreatic cancer subtypes. Int J Mol Sci (2022) 23(9):4792. doi: 10.3390/ijms23094792
61. Pratt HG, Steinberger KJ, Mihalik NE, Ott S, Whalley T, Szomolay B, et al. Macrophage and neutrophil interactions in the pancreatic tumor microenvironment drive the pathogenesis of pancreatic cancer. Cancers (Basel) (2021) 14(1):194. doi: 10.3390/cancers14010194
62. Zhou J, Tang Z, Gao S, Li C, Feng Y, Zhou X. Tumor-associated macrophages: Recent insights and therapies. Front Oncol (2020) 10:188. doi: 10.3389/fonc.2020.00188
63. Rice CM, Davies LC, Subleski JJ, Maio N, Gonzalez-Cotto M, Andrews C, et al. Tumour-elicited neutrophils engage mitochondrial metabolism to circumvent nutrient limitations and maintain immune suppression. Nat Commun (2018) 9(1):5099. doi: 10.1038/s41467-018-07505-2
64. Ratajczak MZ, Kucia M. The Nlrp3 inflammasome - the evolving story of its positive and negative effects on hematopoiesis. Curr Opin Hematol (2021) 28(4):251–61. doi: 10.1097/MOH.0000000000000658
65. Atay S. Integrated transcriptome meta-analysis of pancreatic ductal adenocarcinoma and matched adjacent pancreatic tissues. PeerJ (2020) 8:e10141. doi: 10.7717/peerj.10141
66. Liu Y, Zhu D, Xing H, Hou Y, Sun Y. A 6gene risk score system constructed for predicting the clinical prognosis of pancreatic adenocarcinoma patients. Oncol Rep (2019) 41(3):1521–30. doi: 10.3892/or.2019.6979
67. Sheng Z, Han W, Huang B, Shen G. Screening and identification of potential prognostic biomarkers in metastatic skin cutaneous melanoma by bioinformatics analysis. J Cell Mol Med (2020) 24(19):11613–8. doi: 10.1111/jcmm.15822
68. Sasahira T, Kurihara-Shimomura M, Shimomura H, Bosserhoff AK, Kirita T. Identification of oral squamous cell carcinoma markers MUC2 and SPRR1B downstream of TANGO. J Cancer Res Clin Oncol (2021) 147(6):1659–72. doi: 10.1007/s00432-021-03568-9
69. Zhang Z, Shi R, Xu S, Li Y, Zhang H, Liu M, et al. Identification of small proline-rich protein 1B (SPRR1B) as a prognostically predictive biomarker for lung adenocarcinoma by integrative bioinformatic analysis. Thorac Cancer (2021) 12(6):796–806. doi: 10.1111/1759-7714.13836
70. Guan G, Niu X, Qiao X, Wang X, Liu J, Zhong M. Upregulation of neural cell adhesion molecule 1 (NCAM1) by hsa-miR-141-3p suppresses ameloblastoma cell migration. Med Sci Monit (2020) 26:e923491. doi: 10.12659/MSM.923491
71. Balachandran VP, Luksza M, Zhao JN, Makarov V, Moral JA, Remark R, et al. Identification of unique neoantigen qualities in long-term survivors of pancreatic cancer. Nature (2017) 551(7681):512–6. doi: 10.1038/nature24462
72. Scott EN, Gocher AM, Workman CJ, Vignali DAA. Regulatory T cells: Barriers of immune infiltration into the tumor microenvironment. Front Immunol (2021) 12:702726. doi: 10.3389/fimmu.2021.702726
73. Roulleaux Dugage M, Nassif EF, Italiano A, Bahleda R. Improving immunotherapy efficacy in soft-tissue sarcomas: A biomarker driven and histotype tailored review. Front Immunol (2021) 12:775761. doi: 10.3389/fimmu.2021.775761
74. Huang B, Lei Z, Zhang GM, Li D, Song C, Li B, et al. SCF-mediated mast cell infiltration and activation exacerbate the inflammation and immunosuppression in tumor microenvironment. Blood (2008) 112(4):1269–79. doi: 10.1182/blood-2008-03-147033
75. Varricchi G, Galdiero MR, Loffredo S, Marone G, Iannone R, Marone G, et al. Are mast cells MASTers in cancer? Front Immunol (2017) 8:424. doi: 10.3389/fimmu.2017.00424
76. Xiong Y, Liu L, Xia Y, Qi Y, Chen Y, Chen L, et al. Tumor infiltrating mast cells determine oncogenic HIF-2alpha-conferred immune evasion in clear cell renal cell carcinoma. Cancer Immunol Immunother (2019) 68(5):731–41. doi: 10.1007/s00262-019-02314-y
77. Lu Y, Zhou X, Liu Z, Wang W, Li F, Fu W. Characteristic analysis of featured genes associated with stemness indices in colorectal cancer. Front Mol Biosci (2020) 7:563922. doi: 10.3389/fmolb.2020.563922
78. Yi L, Huang P, Zou X, Guo L, Gu Y, Wen C, et al. Integrative stemness characteristics associated with prognosis and the immune microenvironment in esophageal cancer. Pharmacol Res (2020) 161:105144. doi: 10.1016/j.phrs.2020.105144
79. Miranda A, Hamilton PT, Zhang AW, Pattnaik S, Becht E, Mezheyeuski A, et al. Cancer stemness, intratumoral heterogeneity, and immune response across cancers. Proc Natl Acad Sci U.S.A. (2019) 116(18):9020–9. doi: 10.1073/pnas.1818210116
80. Addeo A, Friedlaender A, Banna GL, Weiss GJ. TMB or not TMB as a biomarker: That is the question. Crit Rev Oncol Hematol (2021) 163:103374. doi: 10.1016/j.critrevonc.2021.103374
81. Wilson RB, Archid R, Reymond MA. Reprogramming of mesothelial-mesenchymal transition in chronic peritoneal diseases by estrogen receptor modulation and TGF-beta1 inhibition. Int J Mol Sci (2020) 21(11):4158. doi: 10.3390/ijms21114158
82. Wu X, Cao Y, Xiao H, Li C, Lin J. Bazedoxifene as a novel GP130 inhibitor for pancreatic cancer therapy. Mol Cancer Ther (2016) 15(11):2609–19. doi: 10.1158/1535-7163.MCT-15-0921
83. Edling CE, Selvaggi F, Ghonaim R, Maffucci T, Falasca M. Caffeine and the analog CGS 15943 inhibit cancer cell growth by targeting the phosphoinositide 3-kinase/Akt pathway. Cancer Biol Ther (2014) 15(5):524–32. doi: 10.4161/cbt.28018
84. Tao S, Tian L, Wang X, Shou Y. A pyroptosis-related gene signature for prognosis and immune microenvironment of pancreatic cancer. Front Genet (2022) 13:817919. doi: 10.3389/fgene.2022.817919
85. Xie W, Li X, Yang C, Li J, Shen G, Chen H, et al. The pyroptosis-related gene prognostic index associated with tumor immune infiltration for pancreatic cancer. Int J Mol Sci (2022) 23(11):6178. doi: 10.3390/ijms23116178
86. Xu X, Liang JH, Li JH, Xu QC, Yin XY. Values of a novel pyroptosis-related genetic signature in predicting outcome and immune status of pancreatic ductal adenocarcinoma. Gastroenterol Rep (Oxf) (2022) 10:goac051. doi: 10.1093/gastro/goac051
87. Yan C, Niu Y, Li F, Zhao W, Ma L. System analysis based on the pyroptosis-related genes identifies GSDMC as a novel therapy target for pancreatic adenocarcinoma. J Transl Med (2022) 20(1):455. doi: 10.1186/s12967-022-03632-z
Keywords: pancreatic cancer, pyroptosis molecular subtypes, prognostic signature, tumor immune microenvironment, tumor stemness, chemotherapeutic drug sensitivity, small molecule compounds
Citation: Huang Q, Peng X, Li Q, Zhu J, Xue J and Jiang H (2023) Construction and comprehensive analysis of a novel prognostic signature associated with pyroptosis molecular subtypes in patients with pancreatic adenocarcinoma. Front. Immunol. 14:1111494. doi: 10.3389/fimmu.2023.1111494
Received: 29 November 2022; Accepted: 23 January 2023;
Published: 03 February 2023.
Edited by:
Yutian Zou, Sun Yat-sen University Cancer Center (SYSUCC), ChinaReviewed by:
Jiang Deng, Institute of Health Service and Transfusion Medicine, ChinaLiu Le Ping, Central South University, China
Copyright © 2023 Huang, Peng, Li, Zhu, Xue and Jiang. 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: Hua Jiang, huajiang2013@tongji.edu.cn
†These authors have contributed equally to this work