In-depth study of pyroptosis-related genes and immune infiltration in colon cancer

Data acquisition and pre-processed

Expression profiling data for CC (The Cancer Genome Atlas-Colon Adenocarcinoma (TCGA-COAD)) were downloaded from the UCSC Xena dataset. Normalised counts were converted to transcripts per million (TPM) values by selecting ‘Count’ as the data type. In total, 521 human colon tumor transcriptomic datasets from TCGA, comprising 480 tumours and 41 normal samples, were included in the analysis. Patient data included age, sex, and TNM stage, among others. Subsequently, the transcriptome data were matched with patient ID extracted from the clinical information. After excluding incomplete survival information and unmatched transcriptome samples, data from 431 patients were included. Genome annotation data were obtained from GENCODE (https://www.gencodegenes.org/), data types were converted based on Ensembl ID, and relevant coding genes were extracted.

Cluster analysis

Based on the literature, a total of 33 pyroptosis genes were considered pyroptosis gene sets (Fu & Song, 2021; Man & Kanneganti, 2015). Unsupervised clustering does not rely on predefined labels, allowing the revelation of unknown patterns and subgroups, thereby providing novel biological hypotheses and perspectives. Therefore, to better identify potential biological subgroups and uncover the internal structure of data, unsupervised cluster analysis was performed based on these genes and the expression data from TCGA-COAD using the ‘ConsensusClusterPlus’ R package. The consensus clustering algorithm was also employed to determine the optimal number of clusters by analysing 100 iterations, thus ensuring classification stability. Principal component analysis (PCA) is a multivariate statistical analysis tool for unsupervised learning. The proportion of ambiguous clustering (PAC) is an important index for assessing the clarity of group division, where a lower PAC score indicates a clearer classification result and a more distinct sample allocation. Therefore, in this study, the consensus clustering method was combined with PAC to determine the optimum number of clusters (k). Based on the consensus clustering results, the samples were divided into subgroups A and B. Subsequently, the gene expression of PRGs was compared between these two subgroups.

Differential analysis and PPI analysis

To identify PRGs, the R package DESeq2 (Ritchie et al., 2015) was used to identify DEGs across different pyroptosis subtypes. This differential expression was displayed in a volcano plot drawn using the ggplot2 package, and a heatmap was constructed using the pheatmap package. The selection conditions were —log2(fold change) —>1.0 and p <0.05. To clarify their interactions, PPI analyses of the above mRNAs were performed using STRING (Szklarczyk et al., 2019) version 11, which is a biological database widely used for PPI analysis.

Functional enrichment analysis and gene set variation analysis

As a common method for large-scale functional enrichment research, Gene Ontology (GO) enrichment analyses (Ashburner et al., 2000) include biological processes, molecular function, and cellular components. Kyoto Encyclopedia of Genes and Genomes (KEGG) (Ogata et al., 1999) is an important gene database containing data on genomes, pathways, and drugs. Therefore, GO and KEGG analyses were performed on DEGs using the clusterProfiler package (Yu et al., 2012) in R, with p < 0.05 as the significance judgement criterion. Based on the ‘c2.cp.kegg. v7.0. symbols’ gene sets, a gene set variation analysis was conducted using gene set variation analysis (GSVA) (Hänzelmann, Castelo & Guinney, 2013). Single-set gene set enrichment analysis (ssGSEA) was used to calculate the score of KEGG pathways according to the gene expression matrix for each sample separately, and the results were displayed as a heatmap.

Differential gene expression analysis

The chromosome location information and TPM expression values of HES4, TNNC1, RNF208, ADAM12, DRD4, POLR2L, RHOD, WDR24, ZNF771, and PRMT1 were extracted from TCGA-COAD transcriptome data. The Rcircos package was used to generate a chromosome ring map to illustrate the chromosome distribution of the 10 genes, and the Wilcoxon rank-sum test was used to compare the expression differences of these genes in the tumor and normal groups and subtypes A and B. In addition, the ‘circlize’ package was used to generate correlation circles among 10 genes, with gene correlations assessed using Spearman’s correlation test.

Immune infiltration analysis

Deconvolution analysis was performed on the transcriptome expression matrix using CIBERSORT to estimate the composition and abundance of immune cells in mixed cells. Gene expression matrix data were uploaded to CIBERSORT, and the samples were filtered based on p <0.05 to construct an immune infiltration matrix. The distribution of the 22 immune cells infiltrating each sample was displayed in a histogram using the ggplot2 package in R.

Immunohistochemistry

This study evaluated the expression of TNNC1 in CC tissue samples using immunohistochemistry (IHC). CC tissue samples were obtained from volunteers who underwent CC surgery at The Second Hospital of Dalian Medical University (Dalian, Liaoning Province, China) between January 2020 and December 2022. All participants were informed in writing about the purpose of the study, the guarantee of anonymity, and the methods of data collection and storage. Each participant signed an informed consent form. This study was approved by the ethics committee of our hospital (Licence number: 2022-036).

Paraffin sections of human colonic tissues and CC tissues were dewaxed and rehydrated according to standard procedures. Antigen retrieval was then performed in EDTA buffers. The sections were pretreated with 3% hydrogen peroxide for 10 min and rinsed three times in phosphate-buffered saline (PBS), with 3 min per rinse. Finally, the sections were incubated overnight with anti-TNNC1 antibody (1:400, 21652-1-AP; Proteintech, Rosemount, IL, USA) at 4 °C. The next day, the protocol of a universal two-step test kit (PV-9000, Beijing Zhongshan Jinqiao Biological) was followed for the subsequent steps. Briefly, after washing to remove the primary antibody, a reaction enhancer was added to the sections, which were then incubated at room temperature for 20 min and rinsed with PBS three times (3 min, each). Subsequently, the sections were incubated with a secondary antibody (enzyme-labelled goat anti-mouse/rabbit IgG polymer) for 20 min, and colour development was achieved using nickel-diaminobenzidine (DAB). The results were assessed based on the IHC score; total IHC scores were calculated by multiplying the ‘staining intensity score’ (0: negative, 1: weak, 2: moderate, and 3: strong) by the ‘positive rate’ (0–300%) (Li et al., 2020). For TNNC1 expression level analysis, an expression level ≤125% was classified as the low-level group, while a score >125% was classified as the high-level group.

Cell culture and validation of the expression level of TNNC1

Various CC cell lines were culture, and the tumor cells with the highest TNNC1 expression were selected to construct stably transfected cell lines. All cells were cultured in a 37 °C (5% carbon dioxide) incubator. CCD841CoN (Meisen, Zhejiang, China), SW480, and SW620 colon cancer cells (Procell, Wuhan, China) were maintained in DMEM supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin (PS). HT29 and HCT 116 cells were purchased from Procell (Wuhan, China) and grown in McCoy’s 5A medium supplemented with 10% FBS and 1% PS, respectively. The expression levels of TNNC1 mRNA were compared between normal human colon epithelial cell lines and human CC cell lines using real-time quantitative PCR. TransZol Up (ET111-01-V2; TransGen Biotech, Beijing, China) was used to extract total RNA from the cell lines (normal colon epithelial cells and four colorectal cancer cell lines). The A260/A280 ratio of the extracted RNA ranged from 1.8 to 2.0. Total RNA concentration was measured using the Cytology 3 Imaging Reader, and an appropriate amount of RNase-free water was added to dilute the RNA to 1 µg/µL. The PerfectStart^® Uni RT & qPCR kit (AUQ-01; TransGen Biotech) was used to synthesise cDNA strands by adding 1 µL of total RNA to the reaction mixture. qPCR was then conducted to detect the differential TNNC1 expression in cells using the Applied Biosystems Step One system (Applied Biosystems, Waltham, MA. USA). The relative expression of TNNC1 was determined using the 2^−ΔΔCT method. All consumables, including gun tips, centrifuge tubes, and organic solvents used throughout the experimental process, were free of RNase enzyme contamination.

The following primer sequences were used: TNNC1, forward primer: 5′-CTACAAGGCTGCGGTAGAGC-3′, and reverse primer: 5′-CCCAGCACGAAGATGTCGAA-3′; for GAPDH, forward primer: 5′-GACAGTCAGCCGCATCTTCT-3′, and reverse primer: 5′-GCGCCCAATACGACCAAATC-3′.

Statistical analysis

All bioinformatics analyses were carried out with R software version 4.0.2. Statistical analyses for IHC and qRT-PCR were conducted using GraphPad Prism 8. The significance of the variables in the two groups was assessed using the Chi-square test and F-test. Non-normally distributed variables were compared using the Mann–Whitney U test (Wilcoxon rank-sum test), while the independent Student’s t-test was used to assess the statistical significance of normally distributed variables.

Identification of pyroptosis subtypes

We conducted an unsupervised clustering analysis to identify two stable subgroups (Fig. 1). We observed consensus clustering from k = 2 to k = 4 in the sample clustering heatmap and tracing curves (Figs. 1A and 1B). When the unsupervised K-mean was 2, the heatmap indicated clear sample clustering, with the highest e discrimination between groups. Moreover, the consistency matrix demonstrated the strongest consensus between the groups, with no significant cross-group confounding. In addition, PCA analysis revealed that when k = 2, groups A and B were well separated in principal component space, indicating distinct gene expression signatures (Fig. 1C). The PAC score was significantly lower at k = 2 than at other k-means (such as k = 3 and k = 4), further supporting the validity of selecting k = 2. We then compared the expression of PRGs between the different subgroups and observed significant differences between subgroups A and B. Patients in the two subgroups exhibited different pyroptosis states (Fig. 1D).

Differential genes and interaction network of differential pyroptosis subgroups

To assess the transcriptomic differences among the multiple regulatory modes of pyroptosis, we performed a differential analysis for each mRNA between the different subgroups, revealing 445 pyroptosis-associated DEGs (Figs. 2A and 2B). Subsequently, we constructed a PPI network for these DEGs, setting the minimum required interaction score to 0.7 to ensure high confidence. This filtering yielded 441 nodes (mRNA) and 204 edges (interaction relationships), with a mean node degree of 0.925. The interactions observed exceeded random distribution, indicating that these proteins, as a group, are at least partially biologically connected. Using k-means clustering, we divided the entire protein network into three clusters: clusters 1, 2, and 3 comprising 144, 115, and 182 proteins, respectively. This result suggests that these pyroptosis-related proteins may participate in pyroptosis via three distinct pathways (Fig. 2C).

GO and KEGG enrichment analyses of pyroptosis-related DEGs

We then performed GO analysis and KEGG pathway enrichment analyses to analyse the functions and pathways associated with these DEGs. GO analysis (Figs. 3A–3C; Table S1) revealed that these DEGs participated in several biological processes such as developmental processes and hormone regulation. They are also associated with cellular components such as the extracellular region, extracellular space, and integral components of the plasma membrane. The DEGs were also enriched in various molecular functions, including receptor regulator activity, neuropeptide receptor activity, structural molecule activity, and G protein-coupled peptide receptor activity. The KEGG pathway analysis (Fig. 3D; Table S2) revealed that these DEGs were enriched in several pathways, including the complement and coagulation cascades, neuroactive ligand–receptor interaction, PPAR, and Wnt pathways.

Gene set variation analysis of pyroptosis-related DEGs

We conducted a KEGG pathway analysis of the mRNA expression profiles of the different pyroptosis subtypes using GSVA. After obtaining the KEGG scores, we screened for the differential pathways between the subtypes using the limma package. GSVA analysis revealed differences in several KEGG pathways between the two subtypes; a total of 50 significantly enriched KEGG pathways are presented in the corresponding heatmaps in Fig. 4, including the transforming growth factor β (TGF-β) signaling pathway and autophagy regulation.

Differential gene expression analysis

We further screened candidate genes closely related to the prognosis of patients with COAD through survival and differential expression analyses. By combining these findings with existing literature on genes related to pyroptosis, we finally identified 10 genes: HES4, TINNC1, RNF208, ADAM12, DRD4, POLR2L, RHOD, WDR24, ZNF771, and PRMT1. These genes are closely related to pyroptosis processes and are significantly differentially expressed in COAD. Using a circos diagram, we illustrated the chromosomal distribution of these 10 mRNAs (Fig. 5A). Moreover, the expression levels of these 10 mRNAs were significantly different between normal and tumor tissues; specifically, their expression levels were elevated in tumor tissues (Fig. 5B). As shown in Fig. 5C, we analysed the correlation between their expression and various other genes using Spearman’s correlation. We observed significant associations between WDR24 and DRD4 (rho = 0.05, p < 0.05) and between WDR24 and PRMT1 (rho = 0.66, p < 0.05). Therefore, we focused on the expression of these mRNAs in different pyroptosis subgroups and found that all these mRNAs were differentially expressed in various subgroups. Specifically, except for ADAM12, the expression of all mRNAs was elevated in subgroup B, suggesting that these mRNAs may play a role in the regulation of pyroptosis (Fig. 5D).

Immune infiltration

We analysed the composition of immune cell infiltrates in the tumor microenvironment (TME) using the CIBERSORT algorithm. Figure 6A shows the immune infiltration profile of various cells within the TME of COAD, while Fig. 6B shows the correlation between immune cell scores. According to the immune infiltration score, we performed Spearman’s correlation analysis between the expression levels of HES4, TINNC 1, RNF208, ADAM12, DRD4, POLR2L, RHOD, WDR24, ZNF771, and PRMT1 and immune cells. We found that the expression of several mRNAs was correlated with immune cell infiltration (Figs. 6C–6L).

Immunohistochemical analysis of TNNC1

In our initial screening and analysis, TNNC1 exhibited significant differential expression in CC and was closely associated with the pyroptosis process. Additionally, TNNC1 was identified as one of the candidate genes closely related to the prognosis of patients with CC through survival and differential expression analyses. Therefore, we conducted an IHC analysis to verify TNNC1 protein expression in normal colon tissues and CC tissues. In normal colon tissues, we observed that the intestinal mucosa remained intact and the gland arrangement was regular. A low level of TNNC1 expression was observed in paracancer tissues. However, unlike the normal colonic tissue, the arrangement of glands in CC tissue was disordered, and TNNC1 protein staining was positive in endothelial cells (Figs. 7A and 7B).

Expression of TNNC1 in CC cell lines

TNNC1 expression levels were assessed in four CC cell lines by qRT-PCR. As shown in Fig. 7C, compared to the normal colon epithelial cell line (CCD841CoN), TNNC1 expression levels were higher in colon cancer cell lines, including HCT116, HT29, and SW620. Among the CC cell lines, the expression level of TNNC1 was the highest in SW480.