Machine Learning-derived Multi-omics Prognostic Signature of Pyroptosis-related lncRNA with Regard to ZKSCAN2-DT and Tumor Immune Infiltration in Colorectal Cancer

2.1. Data Processing and Pyroptosis-related lncRNA

The TCGA database yielded a total of 473 cases of colon adenocarcinoma. Patients with COAD who lacked complete clinical or survival data (n=21) were removed. Also, the following study contained 452 samples altogether (Table 1). According to their Ensembl IDs, genes were briefly divided into protein-coding and non-coding genes [22]. From GENCODE, the GRCH37 long non-coding RNA annotation file was retrieved. From previously published literature, 52 genes associated with pyroptosis were found [23]. To investigate the relationship between 52 pyroptosis-related genes and 14086 lncRNAs, we used Pearson correlation analysis. In the end, 1730 lncRNAs associated with pyroptosis were verified with correlation coefficients |Cor| > 0.4 and P-0.01.

2.2. Identification of Prognostic Pyroptosis-related Hub lncRNA

We analyzed the 1730 pyroptosis-related lncRNAs in order to investigate their predictive potential. The first step was to identify 277 differential pyroptosis-related lncRNAs. The pyroptosis-related lncRNA was discovered using univariate Cox analysis, along with the LASSO regression. Additionally, the multivariate Cox analysis was also applied. Finally, according to the regression coefficient of multivariate Cox analysis and lncRNA expression, a novel predictive pyroptosis-associated lncRNA risk signature was developed using a total of 7 pyroptosis-related lncRNAs. Risk score = β1 × Expression_lncRNA1+β2 × Expression_lncRNA2+…+ βn × Expression_lncRNAn [24].

2.3. Confirmation of Pyroptosis-related lncRNA Signature

Patients' clinical data were gathered and integrated based on the appropriate risk score. Furthermore, patients were allocated into two groups, namely, high-risk subpopulation, and low-risk subpopulation, according to the optimal risk score. Next, survival analysis comparing two subpopulation was applied via Kaplan-Meier curve. Additionally, to evaluate the predictive ability of the model, the 1-, 3-, and 5-year ROC curves were presented [25]. Furthermore, using Cox regression analysis, both univariate and multivariate, we combined gender, and age, along with stage, TNM, and risk scores, to assess their prognostic markers' accuracy with forest plots.

2.4. Development of ROC and Nomogram

To evaluate the prognostic ability of risk score and other clinical features (including age, gender, stage and TNM), the ROC curves of 1/3/5-year were constructed and plotted to evaluate the AUC value [25]. Additionally, we developed a nomogram including risk scores, age, gender, stage, and TNM to statistically assess the effectiveness of our signature's 1/3/5-year prognostic capability. Furthermore, the capability of the nomogram was detected by plotting a 1/3/5-year calibration curve.

2.5. Role of Pyroptosis-related lncRNA Signature in Clinicopathological Features

We compared the risk score in several clinical feature stratifications based on the clinical characteristics of the patients (i.e., age, gender, pathologic stage, AJCC-TNM.). In order to examine if our signature played a significant influence in the onset and development of COAD, a survival curve was applied to assess the overall survival status of various stratifications.

2.6. Tumor Immune Infiltration and Immune Checkpoint Inhibitors

Based on previous research, Cibersort was utilized to calculate the number of immune cells in each sample [26]. Besides, we conducted ssGSEA analysis to evaluate the function of infiltrating immune cells [27]. Furthermore, TIDE was aimed at assessing the immune responses of each patient [28]. From previously published literature, immune checkpoint genes were found that may be relevant to immune checkpoint therapy [29-32]. In order to investigate the possible use of our model in COAD immune checkpoint therapy, we examined whether our model is associated with 47 immune checkpoint genes.

2.7. Evaluation of Pyroptosis-related lncRNA Expression in Immune Cells

Our study utilized the Tumor Immune Single Cell Hub Database (TISCH) [33], a comprehensive scRNA-seq atlas that was constructed by collecting data from both the Gene Expression Omnibus and ArrayExpress databases. From TISCH, we conducted research to evaluate the expression levels of pyroptosis-related lncRNAs, specifically within immune cells.

2.8. Cell Culture and Q-PCR

NCM460 (human colonic epithelial cell) and three human COAD cell lines (SW480, RKO, and HCT116 cells) were cultured in DMEM medium containing 10% fetal bovine serum (FBS). Quantitative real-time polymerase chain reaction was performed on four different cell lines. TRIzol reagent extracted total RNA according to the instruction manual. Next, total RNA was transferred into cDNA using RNA reverse transcription kit. qRT-PCR was performed on cDNA in the applied Biosystems 7500 Fast Real-Time PCR System to determine the ZKSCAN2-DT level. All samples were made in triplicate, and β-actin functioned as a control. Additionally, primer sequences used in PCR were listed as follows: ZKSCAN2-DT, forward TCCAGCTTGGTTGACAAAGTGAGAC and reverse, CCTCCTCGCC TTGCTCTTAATGC; β-actin, Forward, ATCGTGCGT GACATTAAGGAGAAG and Reverse, AGGAAGGAAGG CTGGAAGAGTG.

2.9. Statistical Analysis

Using R software and the Perl language, all results were examined. To determine the statistical significance of two variables, the Student's t-test was used. Statistical significance was considered as *P<0.05, **P<0.01, ***P<0.001.

3.1. Seven Prognostic Pyroptosis-related Hub lncRNA were Identified

In addition, 52 genes relevant to pyroptosis were compiled from published literature (Table S1 ^{(3.1MB, pdf)} ). To identify the pyroptosis-associated lncRNAs, we used Person correlation analysis. Finally, 1730 lncRNAs associated with pyroptosis were discovered (Table S2 ^{(3.1MB, pdf)} ).

Then, based on the “limma” package, 277 distinct pyroptosis-related lncRNAs were discovered (Table S3 ^{(3.1MB, pdf)} ). On the 277 variables, we performed LASSO regression analysis, along with Cox analysis, both univariate and multivariate (Fig. 1A and B). Ultimately, it was determined that seven pyroptosis-related lncRNAs were biological prognostic markers in COAD patients. All 7 of the pyroptosis-related lncRNAs were shown in Fig. (1C) and Table 2 to be risk factors for the prognosis of COAD. Additionally, TCGA-COAD transcription data revealed that these pyroptosis-related lncRNAs were significantly expressed differently in COAD patients (Fig. 1D-1J).

3.2. Confirmation of a Pyroptosis-related lncRNA Signature

According to Fig. (2A), the high-risk subpopulation had a more miserable prognosis outcome than the low-risk group. COAD patients in the low-risk category possessed a greater prognosis result and decreased death, according to the distribution of risk scores (Fig. 2B and 2C). We split the group of COAD patients into high-risk subpopulation and low-risk subpopulation according to the optimal risk score value for each COAD patient (Fig. 2D). The ROC curves were examined with the aim of assessing the diagnostic functions of the risk score (Fig. 2E, AUC1-year=0.647, AUC3-year=0.715, AUC5-year=0.714). Additionally, univariate Cox analysis showed that risk score might be a standalone prognostic factor in COAD patients when paired with other clinical parameters (Fig. 2F, HR=1.110, 95% CI: 1.069-1.157). Similarly, the HR value of multivariate Cox regression analysis was 1.112 (Fig. 2G, HR=1.112, 95% CI: 1.069-1.157).

3.3. Predictive Analysis of the Pyroptosis-related lncRNA Signature with ROC Curve and Nomogram

Gender, age, tumor stage, and TNM were gathered as potential prognostic markers in order to confirm the predictive performance of our signature in comparison to other clinical parameters. When we compared our risk signature to other parameters, we discovered that it had a strong and reliable capacity for predicting the overall survival status of COAD atients (Fig. 3A-3C). Additionally, we statistically use the nomogram to demonstrate the model's ability to predict changes in age, gender, stage, and TNM (Fig. 3D). Our risk model displayed independent prognostic prediction potential for 1/3/5-year OS in COAD patients, as demonstrated by the calibration curve (Fig. 3E-3G). All things considered, the risk signature we developed had a strong and precise predictive capacity to direct COAD prediction.

3.4. Association with other Clinicopathological Features

Age and gender, however, were not yet statistically significant (Fig. 4A and 4B). In addition, we sought to look at how our signature related to other clinicopathological traits in TCGA-COAD. Student's t-test results showed that stage and TNM significantly differed in how the risk score was displayed. In other words, the risk score steadily grew as the stage, or TNM advanced (Fig. 4C-F). Additionally, we conducted a stratified analysis to see if a particular discovered lncRNA may predict OS in different COAD subgroups. The COAD patient with a poor prognosis in various categories (e.g., age65, age>65, female, stage, TNM, etc.) is shown in Fig. (4G-4R).

3.5. TME Features and Immune Infiltration in TCGA-COAD

We further explored whether the risk signature we established was associated with tumor microenvironment characteristics and immune infiltration in TCGA-COAD. The ESTIMATE algorithm exhibited that the high-risk group of COAD patients was significantly linked to a lower immune score and elevated tumor purity (Fig. 5A-5C). On the contrary, it was significantly positively associated with B cells naïve and Tregs (Fig. 5C and 5D). Besides, we found that the signature as constructed was a notably negative correlation with neutrophils, mast cells resting, mast cells activated, eosinophils, and T cells CD4 activated (Fig. 5E-5I). Cibersort showed the tumor immune infiltration statutes between high- and low-groups (Fig. 5J). Our ssGSEA results exhibited the differential distribution of 29 immune cell functions in two subgroups (Fig. 5K). The outcomes provided strong evidence that the signature, as established, was significantly associated with TME features and immune infiltration.

3.6. Correlation between the Pyroptosis-related lncRNA Model and all Crucial Immune Checkpoint Inhibitors

We subsequently explored the following six hub immune checkpoint inhibitors related genes: PDCD1LG2, PDCD1, HAVCR2, CTLA4, CD274, and IDO1 [25-27]. The results expressed that the signature we constructed was significantly negatively correlated with the levels of CD274 and IDO1, implying that our model served as a vital function in asses of ICI treatment response in COAD (Fig. 6A and 6B). In addition, TIDE showed that COAD patients with high-risk scores had lower microsatellite instability (MSI) and stronger immune exclusion, further indicating that our risk model could perfectly predict the responsiveness to ICI therapy (Fig. 6C and 6D). We analyzed the correlation between our pyroptosis-related lncRNA model and immune checkpoint inhibitors (ICI) related targets to identify the potential functions of our model on the immunotherapy of COAD patients (Fig. 6E). By and last, we explored and analyzed all crucial immune checkpoint inhibitors-related genes in high- and low-risk groups, and some of them (i.e., CTLA4, TNFRSF9, ICOS, etc.) were markedly differentially expressed (Fig. 6F).

3.7. The Clinical Significance of ZKSCN2-DT in vitro and COAD Study

ZKSCAN2-DT level was one of the most overexpressed of seven prognostic pyroptosis-related lncRNAs in COAD. Herein, we further explored the ZKSCAN2-DT expression in COAD using Q-PCR. We analyzed the ZKSCAN2-DT expression in three COAD cell lines (SW480, RKO, and HCT116 cells) and normal human colonic epithelial cells (NCM460). Fig. (7A) presented that ZKSCAN2-DT was significantly overexpressed in the COAD cell line compared with normal human colonic epithelial cells. Furthermore, paired analysis was employed, and the results showed that ZKSCAN2-DT expression was significantly upregulated in COAD tissues (Fig. 7B). To evaluate the latent role of ZKSCAN2-DT in COAD, the expression matrix of TCGA-COAD was employed for analysis. As shown in Fig. (7C), patients with higher ZKSCAN2-DT had a worse survival prognosis. Combined with the pathological stage, we could find that ZKSCAN2-DT expression was markedly differentially expressed in different pathological stages. That was, the worse the pathological stage, the higher the ZKSCAN2-DT expression (Fig. 7D). Finally, the correlation between ZKSCAN2-DT and 49 ICI-related genes was analyzed. The finding presented that the higher ZKSCAN2-DT expression possessed a stronger level of ICI-related genes, suggesting that ZKSCAN2-DT might play an important role in COAD immunotherapy (Fig. 7E).

3.8. Role of ZKSCN2-DT in Immune Status and TME

To clarify the relationship between ZKSCAN2-DT and TME characteristics in COAD, patients were further divided into the high ZKSCAN2-DT subgroup and the low ZKSCAN2-DT subgroup based on the median expression of ZKSCAN2-DT. ESTIMATE results showed that the subgroup with lower ZKSCAN2-DT expression had significantly higher stromal score, immune score, along with ESTIMATE score compared with the high-ZKSCAN2-DT group, indicating more immune cells and stroma cells in the low expression subgroup (Fig. 8A-8D). CIBERSORT results revealed that B cells, T cells, and M1 macrophages showed a significantly elevated level in the low-ZKSCAN2-DT subtype (Fig. 8E). ssGSEA analysis results showed that the immune infiltration fraction of most immune cells (i.e., APC co-inhibition, B cells, CD⁸⁺Tcells, iDCs, Th1, Type-1 IFN Response, etc.) decreased significantly with the increase of ZKSCAN2-DT expression in TCGA-COAD (Fig. 8F). We used GSE166555 of the TISCH database to evaluate pyroptosis-related lncRNA expression in tumor immune cells. As shown in Fig. (9), the relatively higher expression level of AL161729.4 was observed in malignant cells, suggesting high AL161729.4 expression in malignant colorectal cancer cells. LINC02381 showed higher expression levels in fibroblasts, and AC007128.1 showed less expression in malignant colorectal cancer cells.