Research Paper Volume 14, Issue 22 pp 9020—9036

Molecular subtypes identified by pyroptosis-related genes are associated with tumor microenvironment cell infiltration in colon cancer

Yiting Ling1, , Yinda Wang2, , Chenxi Cao1, , Lianzhong Feng2, , Binzhong Zhang2, , Senjuan Li1, ,

  • 1 Department of Anorectal Surgery, Jiaxing Second Hospital, Jiaxing 314000, Zhejiang, P.R. China
  • 2 Department of Gastrointestinal Surgery, Jiaxing Second Hospital, Jiaxing 314000, Zhejiang, P.R. China

Received: June 3, 2022       Accepted: November 2, 2022       Published: November 16, 2022      

https://doi.org/10.18632/aging.204379
How to Cite

Copyright: © 2022 Ling et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

The important role of pyroptosis in tumor progression has been well characterized in recent years. However, little is known about the impact of tumor pyroptosis characteristics on patient prognosis and tumor microenvironment (TME) as well as efficacy of immunotherapy. In this study, we successfully classified colon cancer samples into three pyroptosis characterizations with different prognosis and TME cell infiltration patterns based on the expression of pyroptosis-related genes. Cluster 2, with the characterizations of immunosuppression, was classified as immune-desert cell infiltration patterns. Cluster 3, with the patterns of immune-inflamed cell infiltration, had the feature of an activated innate and adaptive immunity and significant prolonged survival. The activation of stromal pathways including EMT, angiogenesis and TGF-β in cluster 1 may mediate the impaired immune penetration of this cluster, which was classified as immune-excluded cell infiltration patterns. Our results demonstrated the PyroSig signature was a robust and independent biomarker for predicting patient prognosis. Patients with low PyroSig signature was confirmed to be correlated with treatment advantages and significant prolonged survival in two anti-checkpoint immunotherapy cohorts. This study identified three pyroptosis-related subtypes with distinct molecular features, clinical and microenvironment cell infiltration patterns in colon cancer, which could promote individualized immunotherapy for colon cancer.

Introduction

As one of the most common tumor types in digestive system, the incidence rate and mortality of colon cancer are increasing year by year [13]. Although in recent years, the comprehensive treatment mainly including surgical resection, chemotherapy, and radiotherapy has enriched the treatment methods of colon cancer, the prognosis of patients has not been significantly improved, and the 5-year survival rate is 40% - 60% [4, 5]. The amazing benefits of immunotherapy represented by PD-1/PD-L1 immune checkpoint inhibitors on a variety of solid tumors have provided the new direction and strategy for the treatment of colon cancer [6, 7]. Unfortunately, anti-PD-1/PD-L1 therapy did not show the desired therapeutic effect in patients with colon cancer. Only 10% - 20% of patients could benefit from the anti-PD-1/PD-L1 therapy, which was far from meeting the clinical needs [8, 9]. Significant individual heterogeneity especially the existed intrinsic and adaptive immune resistance in tumor microenvironment (TME) was significantly associated with treatment failure. Therefore, it is urgent to develop new biomarkers to predict the prognosis of patients and determine the treatment regimens. Increasing evidence indicated that the pyroptosis played a crucial role in the initiation and progression of colon cancer via multiple biological pathways. Pyroptosis, a distinct form of programmed cell death distinguished from apoptosis, is characterized by cell swelling and rupture, releasing inflammatory cellular contents, which triggers a robust inflammatory response [10, 11]. A variety of cytokines and stress-related signaling molecules are activated concomitantly with the process of pyroptosis, which promotes immune cell infiltration as well as inflammatory responses. During pyroptosis, activated caspase-1 promotes the production of proinflammatory cytokines such as IL-18 an IL-1β7 to regulate the tumor immune microenvironment. Furthermore, since pyroptosis is an innate immune mechanism, it could likewise suppress tumor progression [1215]. However, the TME cell infiltrating patterns mediated by distinct pyroptosis characteristics in colon cancer remains unknown. Exploring the impact of distinct pyroptosis characteristics on immune cells of the TME will help to increase the understanding of the complexity and heterogeneity of the TME and provide novel directions and strategies for personalized immunotherapy for colon cancer patients. In the present study, we performed genomic analyses to comprehensively assess the pyroptosis characterizations in colon cancer based on the expression of pyroptosis-related genes, and correlated these characterizations with TME immune cell infiltration patterns. Three pyroptosis characterizations with distinct prognostic features and TME cell infiltration patterns were successfully identified in colon cancer based on more than 1000 cases. The constructed pyroptosis-related score signature (PyroSig) could serve as an independent biomarker to predict patient prognosis and efficacy of anti-PD-1/PD-L1 immunotherapy.

Results

Genomic variation of pyroptosis-related genes in colon cancer

Based on the published studies, we in total extracted 31 pyroptosis-related genes from the GEO datasets. We used the STRING to establish the protein-protein interaction network (PPI), and reveal widespread protein interactions among these genes (Figure 1A). The GO enrichment analysis was used to reveal the biological functions of these genes, and we found these genes were remarkably related to pyroptosis and immune regulation pathways including pyroptosis, nterleukin-1 production and cytokine production involved in immune response (Figure 1B). only 114 of 399 samples, accounting for 28.57%, exhibited at least one mutation type among these genes (Figure 1C). We also observed an obvious mutation co-occurring relationship existed among these genes (Figure 1D). In addition, the CNV alteration analyses showed that these pyroptosis-related genes had a relatively low CNV alteration frequency in colon cancer. GSDMB, GSDMA, GSDMC, PLCG1, PYCARD, IL6, AIM2 and GSDME primarily focused on copy number amplification frequencies, whereas CASP9, CASP3, IL18, ELANE, GPX4 and CASP6 had copy number deletions (Figure 1E). For the mRNA expression, a significant difference in expression was observed between tumor and normal tissues (Figure 1F). Considering the difference in expression, we performed survival analyses for these pyroptosis-related genes, and the univariate Cox regression model indicated these genes were significantly associated with prognosis of patients, of these, NOD1, PRKACA and NLRP1 were the protective factors for colon cancer (Figure 1G).

Landscape of pyroptosis related genes in colon cancer. (A) The protein-protein interactions (PPI) network between pyroptosis genes using STRING database. (B) GO functional enrichment analyses for pyroptosis-related genes. (C) The mutation landscape of pyroptosis genes in TCGA-COAD cohort. (D) The mutation co-occurrence and exclusion analyses for these genes. (E) The copy number variation frequency of pyroptosis genes. (F) Expression of pyroptosis genes in tumor and normal samples based on TCGA-COAD cohort. (G) Survival analyses for the pyroptosis genes using univariate Cox regression model.

Figure 1. Landscape of pyroptosis related genes in colon cancer. (A) The protein-protein interactions (PPI) network between pyroptosis genes using STRING database. (B) GO functional enrichment analyses for pyroptosis-related genes. (C) The mutation landscape of pyroptosis genes in TCGA-COAD cohort. (D) The mutation co-occurrence and exclusion analyses for these genes. (E) The copy number variation frequency of pyroptosis genes. (F) Expression of pyroptosis genes in tumor and normal samples based on TCGA-COAD cohort. (G) Survival analyses for the pyroptosis genes using univariate Cox regression model.

Identification of distinct pyroptosis characterization clusters

Considering the detailed clinical annotations in GSE39582 cohorts, we used this cohort for the subsequent analyses. The consensus clustering was utilized to reveal the pyroptosis characterization in colon cancer and classify patients into distinct molecular subtypes. Based on the expression of entire pyroptosis-related genes, we successfully classified all samples into three different pyroptosis characterizations, which was referred to as cluster 1 subtype, cluster 2 subtype and cluster 3 subtype, respectively (Figure 2A). PCA methods presented three completely disjoint populations among three clusters (Figure 2B). Among three subtypes, patients classified as cluster 2 subtype experienced a poor survival, while patients with cluster 3 subtype showed a particularly prominent prolonged survival (Figure 2C). The significant distinction in all pyroptosis-related gene expression was observed among the three clusters (Figure 2D, 2E). In these three pyroptosis characterizations of colon cancer, we found CASP4, GSDMC, IL1B, IL6, NLRC4, NLRP1, NLRP3 and TNF was significantly up-regulated in cluster 1 subtype compared to other two subtypes. This implied that pyroptosis in this cluster was mainly regulated by CASP4, IL1B, IL6, NLRC4, NLRP1, NLRP3 and TNF. The pyroptosis characterization characterized by high expression of this set of genes mediated the immune-excluded phenotype of tumor microenvironment, which led to a relatively poor prognosis. Similarly, cluster 2 subtype presented an upregulation expression of GPX4, GSDMC, NOD1 and PLCG1, whose pyroptosis characterization, regulated by this gene set, led to an immune-desert microenvironment. And cluster 3 subtype showed a high expression of AIM2, CASP1, CASP3, CASP5, CASP6, CASP8, GSDMB, IL18 and TIRAP, whose pyroptosis characterization, regulated by this gene set, resulted in an improved prognosis and immune-inflamed tumor microenvironment (Figure 2D, 2E).

Identification of distinct pyroptosis characterization clusters in colon cancer. (A) Consensus matrices of pyroptosis related genes, for k=3. (B) Principal component analysis for the three clusters based on the pyroptosis gene expression and revealed three entirely disjoint populations in the meta cohort. (C) Survival analyses for the three pyroptosis characterizations including cluster 1, cluster 2 and cluster 3 in the meta cohort. (D) The hierarchical clustering of pyroptosis genes among three pyroptosis characterization clusters. (E) Difference in pyroptosis gene expression between the three clusters.

Figure 2. Identification of distinct pyroptosis characterization clusters in colon cancer. (A) Consensus matrices of pyroptosis related genes, for k=3. (B) Principal component analysis for the three clusters based on the pyroptosis gene expression and revealed three entirely disjoint populations in the meta cohort. (C) Survival analyses for the three pyroptosis characterizations including cluster 1, cluster 2 and cluster 3 in the meta cohort. (D) The hierarchical clustering of pyroptosis genes among three pyroptosis characterization clusters. (E) Difference in pyroptosis gene expression between the three clusters.

TME cell infiltration patterns under three pyroptosis characterizations

In order to evaluate the correlation between TME cell infiltration patterns and pyroptosis characterizations, we calculated the infiltration abundance of immune cells in each cluster. We found a high level of immune cell infiltration in the tumor microenvironment of cluster 1, sharing the similar immune cell infiltration abundance with the cluster 3 (Figure 3A, 3B). However, despite the high immune level of cluster 1, we did not observe that patients with cluster 1 showed significant survival advantages compared with the other two subtypes (Figure 2C). We then applied the GSVA enrichment analysis to explore this phenomenon by revealing the biological behaviors under these distinct three clusters. As shown in Figure 3C, 3D, we noticed that the immune related signaling pathway was significantly activated in cluster 3 subtype such as cytokine-cytokine receptor interaction, Toll like receptor signaling pathway, T cell receptor signaling pathway, and Nod like receptor signaling pathway (Figure 3D). The activity of stromal-related and carcinogenic pathways was remarkably up-regulated in the cluster 1 subtype including ECM receptor interaction and focal adhesion (Figure 3C). While in cluster 2 subtype, we found the immunosuppressive pathways were significantly activated, which was consistent with the results of immune cell infiltration.

TME cell infiltration patterns in distinct pyroptosis characterizations. (A) The abundance of 28 TME cell infiltration among three pyroptosis characterization clusters visualized by heatmap. (B) Differences of 28 TME cell infiltration abundance between three pyroptosis characterization clusters. (C, D) GSVA enrichment showing the activation states of biological pathways in distinct clusters. (C) cluster 1 vs cluster 3; (D) cluster 2 vs cluster 3. (E) Difference in the immune-activation related gene expression among three clusters. (F) Difference in the TGFβ-EMT pathway-related gene expression among three pyroptosis characterization clusters. (G) Difference in the immune-checkpoint related gene expression among three clusters.

Figure 3. TME cell infiltration patterns in distinct pyroptosis characterizations. (A) The abundance of 28 TME cell infiltration among three pyroptosis characterization clusters visualized by heatmap. (B) Differences of 28 TME cell infiltration abundance between three pyroptosis characterization clusters. (C, D) GSVA enrichment showing the activation states of biological pathways in distinct clusters. (C) cluster 1 vs cluster 3; (D) cluster 2 vs cluster 3. (E) Difference in the immune-activation related gene expression among three clusters. (F) Difference in the TGFβ-EMT pathway-related gene expression among three pyroptosis characterization clusters. (G) Difference in the immune-checkpoint related gene expression among three clusters.

We analyzed the expression of chemokines and cytokines that characterize these three pyroptosis characterization clusters to further explore their characteristics of TME immunoregulation (Figure 3E3G) [16]. The results indicated that the mRNA expression of TGFb/EMT pathway was significantly upregulated in cluster 1 subtypes (Figure 3F), and the cluster 1 and cluster 3 subtypes demonstrated increased expression of mRNAs associated with immune-activation transcripts (Figure 3E). In addition, the cluster 3 subtypes showed an upregulated expression of immune checkpoint molecules, which may mediate the immune escape of this subtype (Figure 3G). Based on the above results, we were surprised to find that the three pyroptosis characterizations of colon cancer were consistent with the three immunophenotypes of the tumor microenvironment, which all shared similar immune cell infiltration patterns. The immune cell infiltrate characteristic of cluster 1 subtype was consistent with the immune-excluded phenotype, which was characterized by marked stromal activation of the tumor microenvironment. The cluster 2 subtype was consistent with an immune-desert phenotype, with a microenvironment with little immune cell infiltration. Whereas the tumor microenvironment of cluster 3 subtype phenocopies a relatively high level of infiltration of innate and adaptive immune cells, which was consistent with an immune-inflamed phenotype. We then used the gene signature to further evaluate the classification accuracy. In order to further validate the stability of these findings, we used the ssGSEA to estimate the activity of stromal pathways through stromal signature, and found the activation of angiogenesis pathways, pan-fibroblast transforming growth factor beta response (Pan-F-TBRS), and epithelial-mesenchymal transition (EMT) (Figure 4A). The ESTIMATE algorithm showed the stromal score was significantly higher in cluster 1 subtype compared with other two subtypes (Figure 4B). Subsequent analysis revealed that the majority of patients with a molecular phenotype of CIN presented with the pyroptosis characterizations of cluster 2, whereas the majority of patients with a dMMR molecular subtype presented with the pyroptosis characterizations of cluster 1 (Figure 4C, 4F). The above analyses also demonstrated that TME cell infiltration patterns under the three pyroptosis characterizations were in accordance with three immunophenotypes of tumor including immune-inflamed, immune-excluded and immune-desert phenotypes, indicating the accuracy of our classification.

Construction of pyroptosis related score signature (PyroSig). (A) Variations between three distinct pyroptosis characterization clusters in pathways with stroma activation. (B) ESTIMATE algorithm analyses revealing the overall TME stromal score among three clusters. (C) The proportion of molecular subtypes in the three clusters. (D) The venn diagram showing 1219 overlap DEGs between three clusters. (E) GO functional enrichment analyses for overlap DEGs. (F) The changes of molecular subtypes, pyroptosis characterization clusters and PyroSig, visualized by alluvial diagram. (G) Spearman correlation between the known signatures and PyroSig values. (H) Differences in PyroSig signature across three distinct pyroptosis characterization clusters. (I) Difference of stromal signature between low and high PyroSig groups. (J) Differences in PyroSig across distinct molecular subtypes. (K) Kaplan-Meier curves showing the survival difference between the low and high PyroSig groups. (L) Survival analyses of four subgroups, where patients were stratified according to adjuvant chemotherapy.

Figure 4. Construction of pyroptosis related score signature (PyroSig). (A) Variations between three distinct pyroptosis characterization clusters in pathways with stroma activation. (B) ESTIMATE algorithm analyses revealing the overall TME stromal score among three clusters. (C) The proportion of molecular subtypes in the three clusters. (D) The venn diagram showing 1219 overlap DEGs between three clusters. (E) GO functional enrichment analyses for overlap DEGs. (F) The changes of molecular subtypes, pyroptosis characterization clusters and PyroSig, visualized by alluvial diagram. (G) Spearman correlation between the known signatures and PyroSig values. (H) Differences in PyroSig signature across three distinct pyroptosis characterization clusters. (I) Difference of stromal signature between low and high PyroSig groups. (J) Differences in PyroSig across distinct molecular subtypes. (K) Kaplan-Meier curves showing the survival difference between the low and high PyroSig groups. (L) Survival analyses of four subgroups, where patients were stratified according to adjuvant chemotherapy.

Construction of pyroptosis related signature (PyroSig)

To further explore the potential biological characteristics of the three pyroptosis characterizations, we investigated the transcriptome differences between the three clusters and a total of 1219 DEGs were determined (Figure 4D). We analyzed differentially expressed genes (DEGs) common to the three subtypes to uncover biological pathways that differ among subtypes. A total of 1219 DEGs were identified, and we found that these DEGs were similarly enriched for immune-related signaling pathways, confirming again that pyroptosis characterizations in colon cancer were significantly correlated with TME anti-tumor immune response (Figure 4E). We established the PyroSig signature by the LASSO analysis to further evaluate the role of pyroptosis characterizations in patient survival and TME cell-infiltrating patterns. The coefficient of each signature gene was summarized in Supplementary Table 2 and Supplementary Figure 1. We used the MaxStat algorithm to classify patients into low and high PyroSig groups. The flow of samples including molecular subtypes, pyroptosis clusters and PyroSig was presented by the alluvial diagram (Figure 4F). The PyroSig signature was negatively correlated with DNA damage repair and was positively correlated with stromal activity (Figure 4G). The cluster 1 and 2 pyroptosis characterization displayed the highest median PyroSig, while the cluster 2 showed the lowest median PyroSig (Figure 4H). The stroma in patients with high PyroSig was remarkably activated compared with that in patients with low PyroSig (Figure 4I). The CSC molecular subtype exhibited a highest PyroSig compared to the other three molecular subtypes (Figure 4J). Patients with low Pyrosig scores were associated with significantly prolonged survival (Figure 4K). Furthermore, we found that the Pyrosig signature could similarly be used to predict survival in patients who received adjuvant chemotherapy. Patients with lower Pyrosig were consistently associated with improved survival regardless of whether or not patients received adjuvant chemotherapy. (Figure 4L).

Characteristics of PyroSig in tumor somatic mutation

We used TCGA-COAD cohort as the validation set. Consistent with the GEO cohort, patients with low PyroSig also experienced an improved overall survival compared to those with high PyroSig (Figure 5A). Compared with CIN and invasive molecular subtype, patients with MSI/CIMP exhibited a relatively higher PyroSig (Figure 5B). However, we found MSI and MSS phenotype did not show an obvious difference (Figure 5C). We then explore the somatic mutation landscape of the high and low PyroSig group to reveal the correlation between tumor mutation burden and PyroSig (Figure 5D5F). The tumor mutation burden (TMB) in tumors with low PyroSig were similar with those with high PyroSig (Figure 5F). We plotted the frequency of mutations as well as the types of mutations in the top 30 mutated genes between high and low PyroSig group using waterfall plots (Figure 5D, 5E). Based on the multivariate Cox regression model consisting of sex, age, tumor location and MMR status, we confirmed that the PyroSig was an independent and robust prognostic biomarker to predict the outcomes of patients with colon cancer (Figure 5G). Additionally, the predictive performance of PyroSig signature to evaluate the one, three and five year survival rates, confirmed by the ROC curves, reached 0.741, 0.765 and 0.741, whose reliability was far superior to the traditional pathological evaluation (Figure 5H).

Characteristics of PyroSig in tumor somatic mutation. (A) Kaplan-Meier curves showing the survival analyses of high and low PyroSig groups in TCGA-COAD cohort. (B) Differences in PyroSig between distinct TCGA colon cancer molecular subtypes. (C) Differences in PyroSig between different microsatellite status. (D, E) The waterfall plot showing the differences of TMB landscape between low and high PyroSig groups. (D) High PyroSig group. (E) Low PyroSig group. (F) Differences in tumor burden mutation between low and high PyroSig groups. (G) Multivariate cox regression analysis for PyroSig in predicting patient’s survival. (H) ROC curves showing the predictive values of PyroSig in prognosis.

Figure 5. Characteristics of PyroSig in tumor somatic mutation. (A) Kaplan-Meier curves showing the survival analyses of high and low PyroSig groups in TCGA-COAD cohort. (B) Differences in PyroSig between distinct TCGA colon cancer molecular subtypes. (C) Differences in PyroSig between different microsatellite status. (D, E) The waterfall plot showing the differences of TMB landscape between low and high PyroSig groups. (D) High PyroSig group. (E) Low PyroSig group. (F) Differences in tumor burden mutation between low and high PyroSig groups. (G) Multivariate cox regression analysis for PyroSig in predicting patient’s survival. (H) ROC curves showing the predictive values of PyroSig in prognosis.

Predictive performance of PyroSig in anti-PD-1/PD-L1 immunotherapy

To further determine the prognostic predictive value of the PyroSig signature, we generalized it to other colon cancer cohorts. Based on the GSE17536, GSE37892 and GSE38832 collected from GEO database, we found that patients with lower PyroSig were consistently associated with improved survival compared with those with higher PyroSig, which further confirmed the potential of PyroSig signature as an independent prognostic biomarker in colon cancer (Figure 6A6C). The successful application of immunotherapy on multiple solid tumors, represented by immune checkpoint inhibitors, has provided a new strategy for the comprehensive treatment of colon cancer. In this study, we collected two immunotherapy cohorts to evaluate the role of PyroSig signature in predicting the efficacy of anti-PD-1/PD-L1 immunotherapy. In the IMvigor210 cohort, which investigated the efficacy of anti-PD-L1 regimens, we found the patients with low PyroSig presented a prominent improved clinical response and survival than patients with high PyroSig (Figure 6D6F). Consistent with anti-PD-L1 regimens, PyroSig signature could be also used to predict the outcomes of patients treated with anti-PD-1 regimens (Figure 6G, 6H). The above results showed a potential predictive value of PyroSig in anti-PD-1/PD-L1 immunotherapy.

Role of PyroSig in predicting efficacy of immunotherapy. (A) Kaplan-Meier curves showing the survival analyses of high and low PyroSig groups in GSE17536 cohort. (B) Kaplan-Meier curves showing the survival analyses of high and low PyroSig groups in GSE37892 cohort. (C) Kaplan-Meier curves showing the survival analyses of high and low PyroSig groups in GSE38832 cohort. (D) Kaplan-Meier curves displaying the survival difference of high and low PyroSig groups in IMvigor210 cohort. (E) The ratio of clinical response types in high PyroSig and low PyroSig groups in the IMvigor210 cohort when treated with anti-PD-1 immunotherapy. (F) Differences in PyroSig between different clinical response types in the IMvigor210 cohort. (G) Survival analyses for PyroSig in GSE78220 anti-PD-1 immunotherapy cohort. (H) The ratio of clinical response types in high PyroSig and low PyroSig groups in the GSE78220 cohort when treated with anti-PD-1 immunotherapy. (I) Differences in PyroSig between different clinical response types in the GSE78220 cohort.

Figure 6. Role of PyroSig in predicting efficacy of immunotherapy. (A) Kaplan-Meier curves showing the survival analyses of high and low PyroSig groups in GSE17536 cohort. (B) Kaplan-Meier curves showing the survival analyses of high and low PyroSig groups in GSE37892 cohort. (C) Kaplan-Meier curves showing the survival analyses of high and low PyroSig groups in GSE38832 cohort. (D) Kaplan-Meier curves displaying the survival difference of high and low PyroSig groups in IMvigor210 cohort. (E) The ratio of clinical response types in high PyroSig and low PyroSig groups in the IMvigor210 cohort when treated with anti-PD-1 immunotherapy. (F) Differences in PyroSig between different clinical response types in the IMvigor210 cohort. (G) Survival analyses for PyroSig in GSE78220 anti-PD-1 immunotherapy cohort. (H) The ratio of clinical response types in high PyroSig and low PyroSig groups in the GSE78220 cohort when treated with anti-PD-1 immunotherapy. (I) Differences in PyroSig between different clinical response types in the GSE78220 cohort.

Discussion

Although the diagnosis and treatment of colon cancer have made some progress in recent years with the development of technology, it has not significantly reduced the incidence rate and mortality of colon cancer. In addition, for advanced colon cancer, existing treatment methods are still limited. Therefore, the early diagnosis and treatment as well as the mechanism of recurrence monitoring are particularly important in colon cancer [17]. Due to the individual heterogeneity of colon cancer, some clinical indicators, such as age, stage, image characteristics and blood biomarkers, have limited effect in predicting prognosis and evaluating treatment plans. In recent years, with the wide application of large-scale sequencing technology, evaluating the gene expression level and mutation characteristics has become one of the important means for prognosis monitoring of many solid tumors. However, a single gene expression level is easily affected by various factors in vivo, and its prediction accuracy is often poor. It has become an important means to improve the prediction efficiency to identify novel molecular subtypes and build a multi gene model by combining a group of gene signature with machine learning [18]. Pyroptosis plays a dual role of promoting or anti-tumor in mediating inflammation and tumor progression [1921]. On the one hand, inflammasomes have the function of eliminating microorganisms and maintaining the integrity of intestinal epithelium, which can prevent tumor attraction. In contrast, inflammasomes can also stimulate the production of protective factors of cnacer cells to help cancer cells escape immune killing, in which, immunosuppressive factors such as Il-18 and il-1 β accumulate in the tumor microenvironment, impairing the function of natural killer immune cells and mediating the immunosuppressive microenvironment [22]. Yang et al. reported that cisplatin exerted antitumor roles in triple-Negative Breast Cancer through promoting MEG3/NLRP3/caspase-1/GSDMD pathway to mediate pyroptosis [23]. Zhang et al. revealed that GSDME silencing remarkably suppressed cisplatin induced secondary necrosis/pyroptosis, but could not inhibit paclitaxel induced secondary necrosis/pyroptosis [24]. At the same time, the antineoplastic features of pyroptosis have been widely demonstrated, possibly due to the protective effect of inflammasomes on the gastrointestinal epithelium. However, the TME cell-infiltrating patterns mediated by different pyroptosis characteristics in colon cancer remains unknown. In this study, we comprehensively evaluated the molecular characteristics of genes related to pyroptosis to construct new classifiers and prognostic signature, and associated them with immune cell infiltration in the tumor microenvironment of colon cancer, so as to guide more appropriate treatment strategies.

Here, we integrated multi-omics data on colon cancer from the TCGA database to further explore the genomic signatures of cell-coke-death-related genes in colon cancer tissues. Although the copy number variation and mutation frequency of these genes are relatively low, the expression levels of these genes are significantly different between colon cancer tissues and normal tissues. Through consensus clustering of these genes, we revealed three distinct cellular pyroptosis characterizations, implying that colon cancer patients could be divided into three distinct molecular subtypes based on pyroptosis-related genes. GSVA and GO enrichment analysis showed that these three subtypes were highly correlated with immune related signaling pathways. The evaluation of TME immune cell infiltration confirmed that these three pyroptosis related molecular subtypes had significantly different TME cell infiltration patterns, and were consistent with the three immune phenotypes of tumors. In general, innate immune cells and stroma were significantly activated in cluster 1 subtypes, which was consistent with the immune-excluded phenotype. Cluster 2 subtype presented an inhibitory microenvironment with relatively few innate and adaptive immune cells, so it was classified as immune-desert phenotype. Contrary to the former two, the microenvironment in cluster 3 subtype had a large number of immune cells infiltrating, and the immune related signal pathway also showed an activated phenotype, which was consistent with the immune-inflamed phenotype. Although both cluster 1 and cluster 3 subtypes belong to “cold tumors”, their mechanisms were not identical. The stroma in the TME was significantly activated in cluster 1 subtype, mediating tumor immune escape. There were abundant immune cells in both immune-excluded and immune-inflamed tumors. However, unlike the immune-inflamed tumors, the activated stroma retained immune cells around the nests of tumor cells, limiting immune cell infiltration into the parenchyma of the tumor. The interaction between the stroma and the immune cells made the immune cells appear to be present inside the tumor. Although PD-1/PD-L1 inhibitors could stimulate T cell activation and proliferation around the stroma, they could not stimulate infiltration, thus limiting the clinical response rate of this subtype to immunotherapy [2529]. In addition, activation of the classical oncogenic pathway in cluster 1 subtype further reduced the level of immune cell infiltration [30]. Our results suggested that impaired immune permeability of the cluster 1 subtype may be associated with activation of EMT, TGF-β, and angiogenic pathways.

At present, it has become a hot topic to transform “cold tumor” into “hot tumor” to increase immune infiltration of microenvironment [3133]. In this study, we identified three molecular subtypes closely related to the immune phenotype through pyroptosis. This suggested that pyroptosis may be a factor that could not be ignored in mediating the complexity of immune infiltration in microenvironment. Regulating pyroptosis-related genes to change the pyroptosis characteristics may be a potential strategy to reshape the microenvironment. We used LASSO algorithm to construct PyroSig signature to further characterize the role of pyroptosis in prognosis and TME cell infiltration patterns. Cluster 3 subtypes were characterized by low PyroSig signature, while cluster 1 and cluster 2 subtypes were characterized by high PyroSig signature. Subsequent analysis revealed that PyroSig signature could be used as an independent biomarker to predict the prognosis of colon cancer patients. We also found that pyroptosis could be involved in the patient's resistance to immune checkpoint immunotherapy. In patients receiving anti-PD-1 and anti-PD-L1 immunotherapy, low PyroSig was closely associated with enhanced clinical response and significantly prolonged survival. Although our results indicated the pyroptosis-related signature could be associated with efficiency of immunotherapy, due to the lack of single-cell sequencing data related to pyroptosis and immunotherapy, we were unable to further determine which cell subset of pyroptosis mediated the efficacy of immunotherapy. Determining which cell subset of pyroptosis played a crucial role in immunotherapy could help to further guide more effective immunotherapeutic strategies. This was a potential limitation of our study.

Conclusions

This study identified three pyroptosis characterizations with distinct clinical, molecular characteristics and TME cell infiltration patterns, and constructed PyroSig signature based on these pyroptosis characterizations, which could be served as a robust and independent biomarker for predicting patient outcomes and efficacy of immunotherapy. It might help to promote individualized immunotherapy for colon cancer from the perspective of pyroptosis characterizations.

Materials and Methods

Sample datasets collection and processing

After systematic search of GEO and TCGA database, we collected a total of 5 colon cancer cohorts including the TCGA-COAD, GSE39582, GSE17536, GSE37892, GSE38832, which contain detailed clinical information [3437]. All cohorts from the GEO database were based on the Affymetrix platforms, so we used the affy R package for background correction and normalization [38]. For TCGA-COAD cohort, we downloaded the FPKM value of the original gene expression and converted FPKM into transcripts per kilobase million (TPM) values. We merged the GEO cohorts and used the sva R package to perform batch correction. The detailed clinical information of all the included cohorts was summarized in Supplementary Table 1.

Identification of pyroptosis characterizations in colon cancer

The pyroptosis-related genes were derived from published studies. In order to identify distinct pyroptosis characterization clusters in colon cancer, we used ConsensuClusterPlus R package to execute consensus clustering based on the mRNA expression of these pyroptosis-related genes. Consensus clustering is repeated 1000 times to ensure the stability of molecular classification [39].

Functional annotation analysis

The Gene Set Variation Analysis (GSVA) and GO enrichment analysis was used to execute functional annotation for uncovering the signaling pathways involved in these pyroptosis characterizations. The “c2.cp.kegg.v6.2.symbols” gene set was downloaded from GSEA database [4042]. The limma package was utilized to identify differentially expressed genes (DEGs) between three clusters. The DEGs with P value < 0.05 was considered as significant [43]. We then conducted GO enrichment analysis to reveal the biological pathways associated with these DEGs [44].

Inference of TME cell abundance

Multiple algorithms have now been developed for assessing the abundance of various immune cell infiltrates of the tumor microenvironment based on gene expression. In this study, we used single sample gene set enrichment analysis (ssGSEA) to measure the abundance of microenvironmental immune cell infiltration under three pyroptosis characterizations. We obtained the gene sets with 28 immune cells from a previous study [45, 46].

Establishment of pyroptosis related signature (PyroSig)

Considering that pyroptosis characterizations played an important role in the prognosis and microenvironment of colon cancer, we constructed pyroptosis related signature. First, a univariate Cox regression model was used to explore the correlation between patient prognosis and the expression of these DEGs. The least absolute shrinkage and selection operator regression method (LASSO) was then performed for the expression of the DEGs with the prognosis P value less than 0.001 [47]. The PyroSig signature was defined as follows:

PyroSig=i=1nCoefiExpri

where i is the expression of prognosis-related genes.

Statistical analysis

We used the Kaplan-Meier method to plot survival curves and used the log-rank tests for significance test based on the survminer R package. In addition, a univariate and multivariate Cox regression models were used to reveal the prognostic value of PyroSig signature as a continuous variable. The MaxStat R package was utilized to determine the optimal cut-off point of PyroSig to classify patients into low and high groups. We executed the Wilcoxon test to calculate the difference significance between two groups. The One-way ANOVA test and Kruskal-Wallis tests were utilized to calculate the difference significance groups of three or more [48]. All data analyses were handled based on the software R 4.0.5. All statistical P-values were two-sided, with a statistical significance of p < 0.05.

Data availability statement

All data related to this work can be acquired from the Gene-Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) and the GDC portal (https://portal.gdc.cancer.gov/).

Author Contributions

YTL, YDW and CXC were responsible for the design of this study. YTL and YDW performed the integration and analyses of the data. YTL, YDW, LZF and BZZ wrote this manuscript. YTL, YDW, LZF, BZZ and SJL revised the manuscript. All authors approved this manuscript.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Funding

This work was supported by the Science and Technology Planning Project of Jiaxing City (2018AY32005).

References

  • 1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018; 68:394–424. https://doi.org/10.3322/caac.21492 [PubMed]
  • 2. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin. 2021; 71:7–33. https://doi.org/10.3322/caac.21654 [PubMed]
  • 3. Siegel RL, Miller KD, Goding Sauer A, Fedewa SA, Butterly LF, Anderson JC, Cercek A, Smith RA, Jemal A. Colorectal cancer statistics, 2020. CA Cancer J Clin. 2020; 70:145–64. https://doi.org/10.3322/caac.21601 [PubMed]
  • 4. Dienstmann R, Salazar R, Tabernero J. Personalizing colon cancer adjuvant therapy: selecting optimal treatments for individual patients. J Clin Oncol. 2015; 33:1787–96. https://doi.org/10.1200/JCO.2014.60.0213 [PubMed]
  • 5. Jung F, Lee M, Doshi S, Zhao G, Cheung KLT, Chesney T, Guidolin K, Englesakis M, Lukovic J, O’Kane G, Quereshy FA, Chadi SA. Neoadjuvant therapy versus direct to surgery for T4 colon cancer: meta-analysis. Br J Surg. 2021; 109:30–6. https://doi.org/10.1093/bjs/znab382 [PubMed]
  • 6. Janjigian YY, Shitara K, Moehler M, Garrido M, Salman P, Shen L, Wyrwicz L, Yamaguchi K, Skoczylas T, Campos Bragagnoli A, Liu T, Schenker M, Yanez P, et al. First-line nivolumab plus chemotherapy versus chemotherapy alone for advanced gastric, gastro-oesophageal junction, and oesophageal adenocarcinoma (CheckMate 649): a randomised, open-label, phase 3 trial. Lancet. 2021; 398:27–40. https://doi.org/10.1016/S0140-6736(21)00797-2 [PubMed]
  • 7. Reck M, Rodríguez-Abreu D, Robinson AG, Hui R, Csőszi T, Fülöp A, Gottfried M, Peled N, Tafreshi A, Cuffe S, O’Brien M, Rao S, Hotta K, et al, and KEYNOTE-024 Investigators. Pembrolizumab versus Chemotherapy for PD-L1-Positive Non-Small-Cell Lung Cancer. N Engl J Med. 2016; 375:1823–33. https://doi.org/10.1056/NEJMoa1606774 [PubMed]
  • 8. Chalabi M, Fanchi LF, Dijkstra KK, Van den Berg JG, Aalbers AG, Sikorska K, Lopez-Yurda M, Grootscholten C, Beets GL, Snaebjornsson P, Maas M, Mertz M, Veninga V, et al. Neoadjuvant immunotherapy leads to pathological responses in MMR-proficient and MMR-deficient early-stage colon cancers. Nat Med. 2020; 26:566–76. https://doi.org/10.1038/s41591-020-0805-8 [PubMed]
  • 9. Overman MJ, McDermott R, Leach JL, Lonardi S, Lenz HJ, Morse MA, Desai J, Hill A, Axelson M, Moss RA, Goldberg MV, Cao ZA, Ledeine JM, et al. Nivolumab in patients with metastatic DNA mismatch repair-deficient or microsatellite instability-high colorectal cancer (CheckMate 142): an open-label, multicentre, phase 2 study. Lancet Oncol. 2017; 18:1182–91. https://doi.org/10.1016/S1470-2045(17)30422-9 [PubMed]
  • 10. Shi J, Gao W, Shao F. Pyroptosis: Gasdermin-Mediated Programmed Necrotic Cell Death. Trends Biochem Sci. 2017; 42:245–54. https://doi.org/10.1016/j.tibs.2016.10.004 [PubMed]
  • 11. Shi J, Zhao Y, Wang K, Shi X, Wang Y, Huang H, Zhuang Y, Cai T, Wang F, Shao F. Cleavage of GSDMD by inflammatory caspases determines pyroptotic cell death. Nature. 2015; 526:660–5. https://doi.org/10.1038/nature15514 [PubMed]
  • 12. Dupaul-Chicoine J, Yeretssian G, Doiron K, Bergstrom KS, McIntire CR, LeBlanc PM, Meunier C, Turbide C, Gros P, Beauchemin N, Vallance BA, Saleh M. Control of intestinal homeostasis, colitis, and colitis-associated colorectal cancer by the inflammatory caspases. Immunity. 2010; 32:367–78. https://doi.org/10.1016/j.immuni.2010.02.012 [PubMed]
  • 13. Kolb R, Liu GH, Janowski AM, Sutterwala FS, Zhang W. Inflammasomes in cancer: a double-edged sword. Protein Cell. 2014; 5:12–20. https://doi.org/10.1007/s13238-013-0001-4 [PubMed]
  • 14. Wang YY, Liu XL, Zhao R. Induction of Pyroptosis and Its Implications in Cancer Management. Front Oncol. 2019; 9:971. https://doi.org/10.3389/fonc.2019.00971 [PubMed]
  • 15. Yu J, Li S, Qi J, Chen Z, Wu Y, Guo J, Wang K, Sun X, Zheng J. Cleavage of GSDME by caspase-3 determines lobaplatin-induced pyroptosis in colon cancer cells. Cell Death Dis. 2019; 10:193. https://doi.org/10.1038/s41419-019-1441-4 [PubMed]
  • 16. Peng D, Kryczek I, Nagarsheth N, Zhao L, Wei S, Wang W, Sun Y, Zhao E, Vatan L, Szeliga W, Kotarski J, Tarkowski R, Dou Y, et al. Epigenetic silencing of TH1-type chemokines shapes tumour immunity and immunotherapy. Nature. 2015; 527:249–53. https://doi.org/10.1038/nature15520 [PubMed]
  • 17. Tariq K, Ghias K. Colorectal cancer carcinogenesis: a review of mechanisms. Cancer Biol Med. 2016; 13:120–35. https://doi.org/10.28092/j.issn.2095-3941.2015.0103 [PubMed]
  • 18. Araghi M, Soerjomataram I, Jenkins M, Brierley J, Morris E, Bray F, Arnold M. Global trends in colorectal cancer mortality: projections to the year 2035. Int J Cancer. 2019; 144:2992–3000. https://doi.org/10.1002/ijc.32055 [PubMed]
  • 19. Diakos CI, Charles KA, McMillan DC, Clarke SJ. Cancer-related inflammation and treatment effectiveness. Lancet Oncol. 2014; 15:e493–503. https://doi.org/10.1016/S1470-2045(14)70263-3 [PubMed]
  • 20. Henao-Mejia J, Elinav E, Strowig T, Flavell RA. Inflammasomes: far beyond inflammation. Nat Immunol. 2012; 13:321–4. https://doi.org/10.1038/ni.2257 [PubMed]
  • 21. Xia X, Wang X, Cheng Z, Qin W, Lei L, Jiang J, Hu J. The role of pyroptosis in cancer: pro-cancer or pro-“host”? Cell Death Dis. 2019; 10:650. https://doi.org/10.1038/s41419-019-1883-8 [PubMed]
  • 22. Terme M, Ullrich E, Aymeric L, Meinhardt K, Desbois M, Delahaye N, Viaud S, Ryffel B, Yagita H, Kaplanski G, Prévost-Blondel A, Kato M, Schultze JL, et al. IL-18 induces PD-1-dependent immunosuppression in cancer. Cancer Res. 2011; 71:5393–9. https://doi.org/10.1158/0008-5472.CAN-11-0993 [PubMed]
  • 23. Yan H, Luo B, Wu X, Guan F, Yu X, Zhao L, Ke X, Wu J, Yuan J. Cisplatin Induces Pyroptosis via Activation of MEG3/NLRP3/caspase-1/GSDMD Pathway in Triple-Negative Breast Cancer. Int J Biol Sci. 2021; 17:2606–21. https://doi.org/10.7150/ijbs.60292 [PubMed]
  • 24. Zhang CC, Li CG, Wang YF, Xu LH, He XH, Zeng QZ, Zeng CY, Mai FY, Hu B, Ouyang DY. Chemotherapeutic paclitaxel and cisplatin differentially induce pyroptosis in A549 lung cancer cells via caspase-3/GSDME activation. Apoptosis. 2019; 24:312–25. https://doi.org/10.1007/s10495-019-01515-1 [PubMed]
  • 25. Chen DS, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature. 2017; 541:321–30. https://doi.org/10.1038/nature21349 [PubMed]
  • 26. Hegde PS, Karanikas V, Evers S. The Where, the When, and the How of Immune Monitoring for Cancer Immunotherapies in the Era of Checkpoint Inhibition. Clin Cancer Res. 2016; 22:1865–74. https://doi.org/10.1158/1078-0432.CCR-15-1507 [PubMed]
  • 27. Herbst RS, Soria JC, Kowanetz M, Fine GD, Hamid O, Gordon MS, Sosman JA, McDermott DF, Powderly JD, Gettinger SN, Kohrt HE, Horn L, Lawrence DP, et al. Predictive correlates of response to the anti-PD-L1 antibody MPDL3280A in cancer patients. Nature. 2014; 515:563–7. https://doi.org/10.1038/nature14011 [PubMed]
  • 28. Joyce JA, Fearon DT. T cell exclusion, immune privilege, and the tumor microenvironment. Science. 2015; 348:74–80. https://doi.org/10.1126/science.aaa6204 [PubMed]
  • 29. Salmon H, Franciszkiewicz K, Damotte D, Dieu-Nosjean MC, Validire P, Trautmann A, Mami-Chouaib F, Donnadieu E. Matrix architecture defines the preferential localization and migration of T cells into the stroma of human lung tumors. J Clin Invest. 2012; 122:899–910. https://doi.org/10.1172/JCI45817 [PubMed]
  • 30. Wellenstein MD, de Visser KE. Cancer-Cell-Intrinsic Mechanisms Shaping the Tumor Immune Landscape. Immunity. 2018; 48:399–416. https://doi.org/10.1016/j.immuni.2018.03.004 [PubMed]
  • 31. Casey SC, Tong L, Li Y, Do R, Walz S, Fitzgerald KN, Gouw AM, Baylot V, Gütgemann I, Eilers M, Felsher DW. MYC regulates the antitumor immune response through CD47 and PD-L1. Science. 2016; 352:227–31. https://doi.org/10.1126/science.aac9935 [PubMed]
  • 32. Kortlever RM, Sodir NM, Wilson CH, Burkhart DL, Pellegrinet L, Brown Swigart L, Littlewood TD, Evan GI. Myc Cooperates with Ras by Programming Inflammation and Immune Suppression. Cell. 2017; 171:1301–15.e14. https://doi.org/10.1016/j.cell.2017.11.013 [PubMed]
  • 33. Sodir NM, Swigart LB, Karnezis AN, Hanahan D, Evan GI, Soucek L. Endogenous Myc maintains the tumor microenvironment. Genes Dev. 2011; 25:907–16. https://doi.org/10.1101/gad.2038411 [PubMed]
  • 34. Laibe S, Lagarde A, Ferrari A, Monges G, Birnbaum D, Olschwang S, and COL2 Project. A seven-gene signature aggregates a subgroup of stage II colon cancers with stage III. OMICS. 2012; 16:560–5. https://doi.org/10.1089/omi.2012.0039 [PubMed]
  • 35. Marisa L, de Reyniès A, Duval A, Selves J, Gaub MP, Vescovo L, Etienne-Grimaldi MC, Schiappa R, Guenot D, Ayadi M, Kirzin S, Chazal M, Fléjou JF, et al. Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value. PLoS Med. 2013; 10:e1001453. https://doi.org/10.1371/journal.pmed.1001453 [PubMed]
  • 36. Smith JJ, Deane NG, Wu F, Merchant NB, Zhang B, Jiang A, Lu P, Johnson JC, Schmidt C, Bailey CE, Eschrich S, Kis C, Levy S, et al. Experimentally derived metastasis gene expression profile predicts recurrence and death in patients with colon cancer. Gastroenterology. 2010; 138:958–68. https://doi.org/10.1053/j.gastro.2009.11.005 [PubMed]
  • 37. Tripathi MK, Deane NG, Zhu J, An H, Mima S, Wang X, Padmanabhan S, Shi Z, Prodduturi N, Ciombor KK, Chen X, Washington MK, Zhang B, Beauchamp RD. Nuclear factor of activated T-cell activity is associated with metastatic capacity in colon cancer. Cancer Res. 2014; 74:6947–57. https://doi.org/10.1158/0008-5472.CAN-14-1592 [PubMed]
  • 38. Gautier L, Cope L, Bolstad BM, Irizarry RA. affy--analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 2004; 20:307–15. https://doi.org/10.1093/bioinformatics/btg405 [PubMed]
  • 39. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010; 26:1572–3. https://doi.org/10.1093/bioinformatics/btq170 [PubMed]
  • 40. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
  • 41. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005; 102:15545–50. https://doi.org/10.1073/pnas.0506580102 [PubMed]
  • 42. Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci. 2012; 131:281–5. https://doi.org/10.1007/s12064-012-0162-3 [PubMed]
  • 43. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]
  • 44. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–7. https://doi.org/10.1089/omi.2011.0118 [PubMed]
  • 45. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, Schinzel AC, Sandy P, Meylan E, Scholl C, Fröhling S, Chan EM, Sos ML, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. 2009; 462:108–12. https://doi.org/10.1038/nature08460 [PubMed]
  • 46. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017; 18:248–62. https://doi.org/10.1016/j.celrep.2016.12.019 [PubMed]
  • 47. Meng Z, Yuan Q, Zhao J, Wang B, Li S, Offringa R, Jin X, Wu H. The m6A-Related mRNA Signature Predicts the Prognosis of Pancreatic Cancer Patients. Mol Ther Oncolytics. 2020; 17:460–70. https://doi.org/10.1016/j.omto.2020.04.011 [PubMed]
  • 48. Hazra A, Gogtay N. Biostatistics Series Module 3: Comparing Groups: Numerical Variables. Indian J Dermatol. 2016; 61:251–60. https://doi.org/10.4103/0019-5154.182416 [PubMed]