Construction and validation of a PANoptosis-related lncRNA signature for predicting prognosis and targeted drug response in thyroid cancer

Functional analysis of PANoptosis-related genes

Based on previous literature and databases, we identified 35 PANoptosis-related genes (PRGs) (Table S1). The selected PRGs were characterized by Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment utilizing the “clusterprofiler” package of R software (version 4.0.5; R Core Team, 2021). The results were visualized and plotted with the help of “enrichplot” and “ggplot2” packages.

Identification of PANoptosis-related lncRNAs

Gene expression data of 571 TC samples (512 tumor and 59 normal) and matching clinical data of 507 TC patients (Patients with incomplete tracking data were not included) were acquired from The Cancer Genome Atlas (TCGA) databank (https://portal.gdc.cancer.gov/). The expression data is available at https://www.jianguoyun.com/p/DeqNJTAQ9pbeCxjIxY0FIAA, and clinical data are listed in Table S2. To identify PRlncRNAs, Pearson correlation was used to determine the correlation between lncRNAs and PRGs in the present study. Correlation coefficient |R²| > 0.4 and p < 0.001 were used as inclusion criteria for PRlncRNAs. The differentially expressed LncRNAs in TC were obtained using the “limma” package.

Establishment of the prognostic model

We randomly distributed all TC patients in the TCGA dataset into the training set (n = 256) and the test set (n = 256) in a 1:1 ratio, and the training set was used for model construction. Data from the test set and entire set were utilized to estimate the performance of the prognostic model. The “Survival” package was applied to conduct univariate/multivariate Cox analysis for PRlncRNAs with p-value ≤ 0.05. Subsequently, we determine the optimal prediction model by means of the least absolute choice operator (LASSO) (Simon et al., 2011). The risk scores of individual patients were computed by incorporating the expression level of lncRNAs in the model and the regression coefficient in the multivariate analysis. The formula is as follows: Σcoefficient of (PRlncRNAi) × expression of (PRlncRNAi).

Validation of model

The samples were categorized as high- or low-risk according to their average risk scores. Kaplan–Meier method was used for survival analysis of high-risk and low-risk groups. The “pheatmap” package was used to map the risk profile, survival status, and lncRNA expression heat map of the groups. The “survivalROC” and “timeROC” packages were used to plot the Receiver Operating Characteristic (ROC) curve and the Area Under the Curve (AUC) to determine the level of accuracy of the signature. Univariate and multivariate Cox regression analyses were performed to verify whether the predictive signature was an independent prognostic variable for Overall Survival (OS). Using the “rms” R package (Iasonos et al., 2008), we combined the OS of TC patients with age, gender, and tumor stage to create a nomogram. A calibration curve was used to test the agreement between the actual results and model predictions according to the Hosmer–Lemeshow test. To compare the accuracy, we calculated the consistency index (c-index) and restricted mean survival (RMS) of the published models (Guo et al., 2022; Qin et al., 2022; Ding et al., 2022; Liu et al., 2022; Shan et al., 2022) and the model we developed.

Gene enrichment analysis

To explore the molecular and biological differences between these two groups, the KEGG and HALLMARK gene sets referenced by the molecular signature database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) were analyzed by the “clusterProfiler” package in the TCGA cohort (p < 0.05, FDR < 0.25). Single-sample GSEA (ssGSEA) was conducted for significant gene sets using the “GSVA” package (Hänzelmann, Castelo & Guinney, 2013).

Estimation of tumor-infiltrating immune cells

THCA immune cell infiltration assessment data were were firstly estimated by the Tumor Immune Estimation Resource (http://timer.cistrome.org/). Based on the TCGA database expression data, we applied the TIMER (Li et al., 2017), CIBERSORT (Newman et al., 2015), CIBERSORT-ABS, QUANTISEQ, MCPcounter (Dienstmann et al., 2019), xCELL (Aran, Hu & Butte, 2017), and EPIC (Racle & Gfeller, 2020) algorithms to obtain the cell composition for the TCGA-THCA dataset samples by R software (Table S3). This enabled comparison of immune cell subpopulations in the groups. ssGSEA was also performed using the “GSVA” packages in R to compare the immune-cell infiltration in the groups. ImmuneScore, StromalScore, and ESTIMATEScore (ImmuneScore + StromalScore) were quantified for each patient using the “limma” and “estimate” packages in R. Finally, the “ggpubr” package in R was used to map immune checkpoints for risk groups (Charoentong et al., 2017).

Clusters based on prognostic PRlncRNAs

Given the PRlncRNA expression associated with the risk signature, we obtained two subtypes using the “ConsensusClusterPlus” package in R. We then evaluated the differences in immune cell infiltration and immune scores between the two subtypes.

Drug sensitivity prediction

The accuracy of the signature in forecasting the response to common antitumor drugs on TC was also evaluated. To determine the therapeutic response, the half inhibitory concentration (IC50) of each sample was estimated with the help of the “pRRophetic” package (Geeleher, Cox & Huang, 2014). Lower IC50 values indicate higher drug sensitivity. The Wilcoxon test was used to calculate the p-values of differences in drug sensitivity.

Quantitative PCR

Ten paired tumor and para-tumor TC tissue samples were acquired from Qilu Hospital of Shandong University. Our study was approved by the Research Ethics Committee of Qilu Hospital of Shandong University (Approval Number: KYLL-2017(KS)-248), and all patients signed a written informed consent form. RNA was isolated from tissues with the total RNA extraction reagent TRIzol (DP424; TIANGEN Biotech, Beijing, China). The obtained RNA was then reverse-transcribed to cDNA using HiScript II Reverse Transcriptase (Vazyme, Nanjing, China). Quantitative PCR assays were performed using a Bio-Rad C1000 Thermal Cycler with SYBR Green Master Mix (Q711; Vazyme Biotech, Nanjing, China), and expression levels were calculated using the 2^−ΔΔCt method. β-actin was used as a standard control. The primers designed for qRT-PCR were obtained from Sangon Biotech (Shanghai, China). And the sequences are presented in Table S4.

Cell culture and transfection

The human thyroid cell lines TPC-1 (RRID: CVCL_6298) was obtained from the Stem Cell Bank, Chinese Academy of Sciences (Shanghai, China). The PTC cell line BHP10-3 (RRID: CVCL_6278) were supplied by Qilu Hospital of Shandong University (Zhang et al., 2019). TPC-1 and BHP10-3 cells were cultured in Roswell Park Memorial Institute (RPMI 1640) complete medium (Shanghai Institute of Biological Sciences, Shanghai, China) at 37 °C in an incubator with a constant temperature and 5% CO₂. The cells were dispersed in six Well Cell Culture Plates and transfected the following day. The overexpression plasmid was obtained from MiaoLing Plasmid Platform (Wuhan, China). When the cell density reached approximately 70%, 2,000 ng of negative control or plasmid was transfected into cells using Lipofectamine™ 2000 Transfection Reagent (11668019; Invitrogen, Waltham, MA, USA). After 6 h, the medium was replaced with 2 mL of complete medium. After 24 h of transfection, cells were collected for subsequent experiments.

Proliferation and migration assay

After 24 h of transfection, the cells were evenly distributed in 96-well plates, five parallel replicate wells were set for each group, and 10 μL of CCK-8 reagent was added on days 1, 2, and 3. The absorbance was measured at 450 and 630 nm, which represents the cell growth rate. Thyroid cancer cells (5 × 10⁴) transfected with RBPMS-AS1 overexpression plasmid or negative control were introduced to the upper chamber of Transwell plates in 200 μL serum-free culture medium, and 800 μL medium comprising 10% fetal bovine serum was introduced into the lower chamber as an inducer. After 24 h, the migrated cells were washed with PBS, fixed in methanol, stained with 1% crystal violet, and non-migrating cells in the upper chamber were removed with a cotton swab. At least four randomly selected fields were observed and counted under a microscope (Olympus, Tokyo, Japan).

Wound healing assay

The cells were incubated in 6-well plates. When the cells were cultured to 80% confluence, a 200 μL pipette tip was used to create a scratch. The cells were then washed with PBS and placed in serum-deficient medium for an additional 24 h. The wound healing area at 0 and 24 h was recorded using a microscope (Olympus, Tokyo, Japan).

Functional enrichment analysis of PANoptosis-related genes

From the previous literature, we obtained 35 PRGs. Based on GO analysis, PRGs primarily participate in biological processes (BP), such as the positive regulation of cytokine production, positive regulation of response to external stimulus, and response to virus. GO analysis of cellular components (CC) revealed that PRGs were mainly positioned in the membrane raft, membrane microdomain, and inflammasome complex. According to molecular function (MF), PRGs mostly carried out the cytokine receptor binding function (Figs. 1A, 1B). KEGG enrichment analysis showed that the selected PRGs were mainly involved in signaling pathways including the NOD-like receptor signaling pathway, Influenza A-mediated pathways, and necroptosis (Figs. 1C, 1D).

Construction of the PANoptosis-related lncRNAs predictive signature in TC

Using Pearson correlation analysis, we recognized 2,192 PANoptosis-related lncRNAs (PRlncRNAs). The Sankey Diagram (Yue et al., 2022) shows the relationship between PRGs and lncRNAs (Fig. S1A). We then obtained 819 differentially expressed PRlncRNAs between the normal and tumor samples (log2|FC| > 1, FDR < 0.05). The volcano map is shown in Fig. 2A. Univariate Cox analysis were performed to investigate the link between lncRNA expression and survival. The analysis indicated that 24 lncRNAs were significantly correlated with the OS of TC patients (Fig. 2B). Notably, 21 PRlncRNAs were “risk” genes, while only LINC0115, DPP4-DT, and RBPMS-AS1 could be considered “protective” genes. The heatmap visualizes PRlncRNA expression in normal and tumor samples (Fig. 2C). LASSO Cox regression analysis was conducted to select and validate prognostic lncRNAs (Figs. 2D, 2E). Finally, the multivariate Cox regression analysis selected seven PRlncRNAs for the construction of the risk signature. A PANoptosis-related prognostic lncRNA signature for TC was created by merging the regression factors and expressive values of seven lncRNAs: IDI2-AS1, LINC02154, DPP4-DT, RBPMS-AS1, C5orf34-AS1, LINC01705, and GAPLINC (Fig. 1F). The expression of the seven lncRNAs was closely correlated with PRGs (Fig. 2F).

Evaluating the predictive capacity of the risk signature

To examine the efficacy of the model, the patients in the public database were assigned to training and validation cohorts. The basic information of the patients in both sets is listed in Table 1. Based on the risk score, the patients in each set were classified into high-risk and low-risk groups. Patients in the high-risk group had poor survival according to the Kaplan–Meier curve analysis (p < 0.001) (Fig. 3A). Figures 4B and 4C show the survival status, distribution of patients, and the mortality rate increased with risk scores. Both testing and the entire set suggested a poor prognosis in high-risk patients (Figs. 3F–3H, 3K–3M). Furthermore, ROC was used to validate the risk signature. The AUC values of this risk signature were 1.000, 0.986, and 0.882 for 1, 3, and 5 years respectively, suggesting good predictive power (Fig. 3D). In addition, the ROC analysis results had high accuracy for both the testing set and the entire set (Figs. 3I, 3M). Figures 3E, 3J and 3O depict the heat maps of the expression levels of the seven PRlncRNAs in all sets.

Association of clinicopathological factors with PRlncRNAs risk signature

To verify the predictive effect on patients with different clinicopathologic variables, patients with TC were stratified in terms of age, gender, and clinical stage. Patients in the high-risk group had notably shorter OS in all subgroups of age (≤65 years), gender (female or male), TNM stage (T1-2 or T3-4, N0 or N1, M0), and clinical stage (stage–III–IV) (Fig. 4). Moreover, the model demonstrated the highest sensitivity and specificity compared to other clinical variables such as age, gender, grade, and stage in all sets (Figs. 5A–5C). To better assess the independent forecasting capacity of the signature, the Univariate/Multivariate Cox regression analyses were conducted on risk scores and patient clinical data. Univariate Cox regression analysis revealed that age, stage, and risk scores acted as independent prognostic indicators (p < 0.001; Fig. 5D), while multivariate independent prognostic analysis showed that only age and risk scores acted as independent prognostic factors for TC patients (p < 0.001; Fig. 6E). We then incorporated risk scores with clinical factors and plotted a predictive nomogram for patients with TC (Fig. 5F). Combined with the calibration plots, it was observed that the predicted OS was consistent with the actual observations (Fig. 5G). In addition, a comparison of the c-index and the RMS of seven previously reported risk models shows that the developed model has a significant advantage (Figs. 5H, 5I). Therefore, the risk signature showed high performance in predicting OS in patients with TC in the TCGA dataset.

Discovery of the association between tumor‑infiltrating immune cells and risk signature

We explored the underlying differences in signaling pathways between the risk groups by Gene Set Enrichment Analysis (GSEA). GSEA of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways showed that the enriched pathways in the high-risk group involved ABC transporters, ECM-receptor interaction pathway, hedgehog signaling athway, primary bile acid biosynthesis, and tryptophan metabolism (Fig. 6A). Conversely, the accumulated pathways in the group at low risk included DNA replication, Nucleotide excision repair, Ribosome, and Spliceosome (Fig. 6B).

To elucidate the relationship between the tumor immune microenvironment (TME) and the risk signature, we performed an analysis of RNA microarray data to explore the correlation between risk characteristics and immune cells. The bar graph shows the association between risk score and tumor-infiltrating immune cells (Fig. 6C). The XCELL algorithm results showed that the risk score was significantly positively correlated with B cell plasma, plasmacytoid dendritic cells, and T cell regulatory cells (Tregs), while it was positively correlated with T cell CD4⁺ Th1, CD8⁺, and mast cells. The results of the QUANTISEQ algorithm showed that the risk score was significantly positively correlated with regulatory T cells (Tregs), while negatively correlated with macrophage M1, neutrophils, etc. The EPIC algorithm showed that the risk score was positively correlated with cancer-associated fibroblasts and negatively correlated with macrophages and NK cells. To better explore the immune status between the two groups, we evaluated the stromal score, the immune score, and the estimated score (Fig. 6D). Subsequently, we investigated the variations between the risk points and tumor infiltrating immune cells using ssGSEA. The scores for B cells, macrophages, and T helper cells were markedly higher in the group at high risk (Fig. 6E). Furthermore, we observed that there was a considerable difference (p < 0.001) between the groups in immune checkpoints such as IDO2, CD44, CD40, HHLA2, ADOARA2A, and CD80 (Fig. 6F).

Cluster analysis according to prognostic PRlncRNAs

In an effort to evaluate the immune microenvironment and response of different subtypes, cluster analysis was performed to obtain clustered subtypes. On the basis of the seven PRlncRNAs, we divided the patients into two subtypes by applying the ConsensusClusterPlus package (Fig. 7A). As shown in the Sankey diagram (Fig. 7B), most of the low-risk patients were grouped in cluster 1, whereas the majority of the high-risk patients were assigned to cluster 2. Survival analysis showed that patients in subgroup 1 had a superior OS (p < 0.001) (Fig. 7C). Based on the cluster analysis, the relationship between the different clusters and TME were examined. Box plots show that cluster 1 had higher immune, stromal, and ESTIMATE scores than cluster 2 (Fig. 7D). The heat map shows the variation in immune cell infiltration across different tumor clusters (Fig. 7E).

Drug sensitivity in the PANoptosis-related LncRNA signature

We further explored the differences in drug resistance between the risk groups. We investigated the IC50 values of 80 chemotherapy drugs and inhibitors in both groups. The values of six typical drugs are shown in Fig. 8. Low-risk patients had substantially lower IC50 values for Imatinib, Sunitinib, Pazopanib, Mitomycin C, 5-Fluorouracil, and tipifarnib, suggesting that the drugs may be eligible for the treatment of patients in the low-risk group. In contrast, drugs such as Bortezomib and Phenformin showed better effects in the high-risk group. This means that candidate immunotherapeutic candidates for patients can be identified with TC based on their risk scores.

Assessing the expression and biological function of PRLncRNAs

Among the PRLncRNAs, we verified the gene levels of four lncRNAs in pairs of samples obtained from patients with TC in our center. Comparable expression patterns were observed in the clinical samples (Figs. 9A–9D). GAPLINC showed increased expression in tumor tissues (T) compared to in normal tissues (N), whereas RBPMS-AS1, IDI2-AS1, and LINC02154 showed lower expression levels. In addition, we also used a GEO database (GSE33630) to verify the expression differences of PRLncRNAs between thyroid cancer and normal tissues. However, we were able to obtain only the expression of RBPMS-AS1, IDI2-AS1, and GAPLINC because of limited data from the sequencing platform (Fig. 9E). As shown in Fig. S2, the RNA expression of RBPMS-AS1 was successfully upregulated in both cell lines after transfection. Functionally, its role in regulating TC cell proliferation was tested using the Cell Counting Kit-8 (CCK-8) assay (Fig. 9F). Overexpression of RBPMS-AS1 inhibited the proliferation of BHP10-3 cells significantly but did not affect TPC-1 cell growth (Fig. 9E). By wound healing assay and transwell experiments, the overexpression of RBPMS-AS1 decreased the migration ability of the two TC cell lines significantly (Figs. 9G, 9H).