Identification and verification of aging-related lncRNAs for prognosis prediction and immune microenvironment in patients with head and neck squamous carcinoma

Introduction

Head and neck squamous cell carcinoma (HNSCC) is a type of malignancy originating from the mucosal epithelium of the lip, oral cavity, larynx, or naso-, oro-, hypo-pharynx. HNSCC has a mortality rate of nearly 50% in the first five years after diagnosis, with over 850,000 new cases and more than 400,000 deaths globally in 2020 [1–3]. Though diverse treatment modalities such as surgery, targeted therapy, chemo-, radio- and immuno-therapy are used in HNSCC patients, their clinical benefits are far from satisfactory [4].

Aging is an irreversible growth arrest condition closely associated with the development of many chronic diseases and cancers [5]. Cellular senescence can be induced by stimuli such as telomere shortening due to extensive replication, oxidative stress, DNA damage, and oncogene overexpression, which are regulated by senescence-related genes [6]. After senescent cells lose proliferative capacity, hindering tumor progression, they remain viable and metabolically active, secreting various cytokines that promote the growth, migration, and invasion of tumor cells [7,8]. The potential prognostic value of aging-related genes in HNSCC and the correlations with inflammation and tumor immunity have recently been explored and confirmed [9].

Long non-coding RNAs (lncRNAs) are a class of non-protein-coding RNAs that are more than 200 nucleotides in length and have specific functions (e.g., splicing, transcriptional and post-transcriptional modulation of mRNA), which are largely linked to tumor development, metastasis, and tumor immunity [10]. The action mechanisms of lncRNAs in aging have not yet been elucidated, and the prognostic impact of lncRNAs associated with aging-related genes in HNSCC remains unknown.

Here, we use bioinformatics to probe into the mechanisms of aging-related lncRNAs (ARLs) in HNSCC, identify some new molecular biomarkers linked with prognosis, and build an effective prognostic model to forecast the survival of HNSCC patients and function as their new treatment target.

Materials and Methods

Data collection

We downloaded the RNA sequencing data and clinical data of HNSCC and normal tissues from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). Aging-related genes were acquired from the Human Ageing Genomic Resources (HAGR; Table S1). Expression data were standardized to fragment per kilobase million. Based on Strawberry Perl, we integrated and processed the data, distinguished between lncRNA and mRNA, extracted complete clinical information, and removed data missing survival status and unknown survival time or <30 days.

Identification of ARLs in HNSCC

We performed co-expression analysis of lncRNAs and AGs in tumor samples (correlation coefficient > 0.5 and p < 0.001) with packages Strawberry Perl and limma R and mapped a network with igraph R. Differential expression of ARLs in cancer versus normal tissue samples was analyzed (Log2 fold change (FC) > 1, false discovery rate (FDR) < 0.05, p < 0.05) [11], and a heatmap was plotted with the R package pheatmap.

Construction of the Prognostic Signature

Based on the clinical data of HNSCC, lncRNAs with significant effects on OS were screened from ARLs (p < 0.05) using univariate Cox analysis. Samples were randomly and evenly separated into a training risk group and a test risk group on the R package caret. In the training group, we conducted the least absolute shrinkage and selection operator (LASSO) Cox regression on the R package glmnet and multivariate Cox analysis to recognize the ARLs for risk models and computed the risk score as follows:

Risk score =Expression of the 1st lncRNA * coefficient+Expression of the 2nd lncRNA * coefficient+Expression of the nth lncRNA * coefficient 

where the coefficient denoted the regression coefficient of the corresponding lncRNA [12]. The overall group was classified by the median risk score into low- and high-risk groups.

Verification of risk signature and establishment of nomogram

We performed the Kaplan–Meier survival curve and log-rank test on the R package survminer to clarify whether the low-risk group has longer OS than the high-risk group. The 1, 3- and 5-year survival receiver’s operating characteristic (ROC) curves and the areas under the curves (AUCs) were conducted on the R package timeROC to assess the prediction efficacy in comparison with other clinical data [13]. Based on the R package rms, a nomogram for forecasting survival was built on basis of the risk score and clinical features, which were detected by both uni-/multi-variate Cox regression. We assessed the accuracy of the nomogram via calibration curve and ROC analysis.

Gene set enrichment analyses (GSEA) and assessment of the immune landscape

We performed GSEA (version 4.2.3) using the curated gene set (kegg.v7.5.1 symbols.gmt), screened at p < 0.05 and FDR < 0.25, to determine the different biofunctions and pathways between the two groups [14]. Immune cell infiltration files for HNSCC in the TCGA were obtained from UCSC Xena (https://xena.ucsc.edu/). Pearson correlation analyses between risk score and immune cell infiltration scores were conducted by the R packages “limma”, “scales”, “ggplot2”, “reshape2”, “tidyverse”, “ggpubr” and “ggtext” and we screened out and visualized the results at p < 0.05. The TME scores of HNSCC were calculated using the ESTIMATE algorithm with estimate R package. Then we assessed immune cell infiltration and immune cell function by single-sample GSEA (ssGSEA), and we accessed immune checkpoint activity by comparing differential expression of immune checkpoint genes.

Prediction of drug therapy response

We adopted the Genomics of Cancer Drug Sensitivity in Cancer (https://www.cancerrxgene.org) and R package pRRophetic to forecast the sensitivity to drug therapies as per the half-maximal inhibitory concentration (IC50) in the two groups. All the code files mentioned above can be downloaded and used at the following website: https://github.com/yexian123/ARLs-HNSCC-analysis.git.

Cell culture and cell transfection

We purchased three human nasopharyngeal carcinoma cell lines (HNE1, CNE1, and CNE2) from the Chinese Academy of Sciences (Shanghai, China). All cells were grown in an RPMI-1640 medium (KeyGEN, Nanjing, China) with 10% fetal bovine serum (FBS, PAN-Seratech, Germany) and 1% penicillin/streptomycin in 5% CO2 at 37°C. The pcDNA3.1 and LINC00861-pcDNA3.1 bought from Genescript Biotechnology were transfected into CNE1 and CNE2 cells using Lipotransfectamine 3000 (Thermo). After 48 h, the CNE1 and CNE2 cells were harvested.

Quantitative real-time PCR (qRT-PCR)

Total RNA from each cell line was extracted with an RNA-easy isolation reagent, and reversed with HiScript® III RT SuperMix for qPCR (+gDNA wiper) (both Vazyme Biotech, China) to synthesize complementary DNA (cDNA). The PCR system was made from SYBR Green® Premix Ex Taq™ (Vazyme Biotech, China). Results of individual lncRNAs were standardized to the GAPDH expression. The primer sequences for LINC00861 and GAPDH are listed in Table S2.

Cell counting kit-8 (CCK-8) assay

Cell proliferation was detected using CCK-8 assay (Beyotime, Shanghai, China) according to the instructions. Cells (1000/well) were planted into 96-well plates and grown in an RPMI-1640 medium containing 10% FBS. At the same time each day, 10 μl of a CCK-8 solution was added to each well, and further incubated for 2 h at 37°C. Absorbance at 450 nm was detected using a microplate spectrophotometer (Thermo, USA) and used to estimate the proliferative capacity of CNE1 and CNE2 cells.

EdU assay

We assessed cell proliferation using an EdU assay kit (RiboBio, China). CNE1 and CNE2 cells (5 × 104 cells/well) were inoculated in a confocal laser cuvette, then cultured in a medium with 50 μM EdU (C10310-1, RiboBio) for 2 h and treated as per the instructions. Images were obtained under a fluorescent microscope. The mean proportion of EdU-positive cells in three random fields of view was analyzed.

Senescence-associated β-galactosidase (SA-β-gal) staining

Cellular SA-β-gal activity were assayed according to the manufacturer’s instructions (C0602, Beyotime, Nanjing, China). CNE1 and CNE2 cells were fixed and stained with the fixative solution and staining solution mixture provided in the kit. Observe and photograph the cells under an ordinary light microscope. The number of SA-β-gal-positive cells in 3 to 5 regions of the 6-well plate was randomly counted as a percentage of the total number of cells, and the mean value was calculated.

Statistical analysis

All statistical analysis and graphical visualization were conducted on R software 4.1.2. Categorical data between groups were compared via the chi-square test. Continuous variables in normal distribution (containing risk score and TME scores) were compared between two or more groups via Student’s t-test or one-way ANOVA. Differences in expression levels between immune checkpoint genes and ARLs were investigated via Wilcox test. The significance level was probability less than 0.05.

Results

Identification of aging-related LncRNAs in HNSCC patients

Fig. 1 illustrates the flowchart of the study. Firstly, we downloaded transcript data of 44 normal and 504 cancer samples from TCGA on 2022-04-15. Then 307 aging-related genes (ARGs) were obtained from HAGR (Table S1). After that, we got 156 AGs and 565 ALRs by correlation analysis (correlation coefficient > 0.5 and p < 0.001).

Figure 1: The flow chart of the whole analysis process.

The regulatory network of these lncRNAs with AGs is shown in Fig. 2A. Based on differential expression analysis of the ALRs in cancer versus normal tissue samples (|Log2FC| > 1 and p < 0.05), we acquired 276 significantly different ALRs, of which 244 were up-regulated and 32 were down-regulated in cancer tissues (Table S3, Fig. 2B).

Figure 2: (A) The regulatory network showed 565 lncRNAs associated with 156 AGs. (B) The 276 significantly different ALRs in cancer versus normal tissue samples.

Establishment and verification of the prognostic model

The 501 eligible tumor samples with complete information were equally divided into a train (N = 251) group and a test group (N = 250). In the training group, univariate Cox regression analysis identified 38 ALRs significantly related with overall survival (OS) (p < 0.05) (Figs. 3A and 3B). To avoid over-fitting of the prognostic features, we used LASSO regression of these lncRNAs and isolated sixteen lncRNAs associated with AGs in HNSCC when the first-order value of Log(λ) was the smallest likelihood of bias (Figs. 3C and 3D), where Log(λ) was −3.96. After that, we performed multi-COX regression analysis and finalized nine ARLs used to construct the model. The Sankey diagram (Fig. 3E) demonstrated that these nine lncRNAs were positively associated with the AGs ADCY5, C1QA, CDKN2A, ELN, GSR, KCNA3, TFDP1 and UCP2.

Figure 3: Selecting eligible ARLs in the training group for building the model. (A) ARLs affecting OS extracted by univariate Cox regression analysis. (B) The heatmap showed the expression profiles of 38 prognostic lncRNAs. (C) The 10-fold cross-validation for variable selection in the LASSO model. (D) The LASSO coefficient profile of 16 ARLs, the intersection of the middle dashed line with the upper X-axis is the number of lncRNAs obtained from the LASSO regression, and the intersection with the lower X-axis is the value of Log(λ). (E) The Sankey diagram of 8 aging-related genes and 9 associated lncRNAs.

Risk scores were calculated as follows: AC246787.2 × (−1.08765834515961) + LINC00996 × (−1.00608978255594) +`CDKN2A-DT` × (−0.754633772241435) + LINC01508 × (0.524759230799357) + AL359502.2 × (0.64309513483616) + Z97653.1 × (−1.22717339963094) + AC243960.1 × (3.02676937595734) + LINC00861 × (−4.95297301646684) + LINC02384 × (−1.88205538422933).

With the risk score formula, we divided the train, test, and entire groups into low- and high-risk groups. We compared the distribution of risk scores, survival status, lncRNA expression, and survival time among the risk groups. All results showed significant differences between the low- and high-risk groups (Figs. 4A–4D). Time-dependent receiver operating characteristics (ROC) were used to assess the sensitivity and specificity of the model for prognosis. We also illustrated the ROC results in terms of the area under the ROC curve (AUC); the 1-, 3- and 5-year AUCs were 0.751, 0.779, and 0.801 for the train set, 0.645, 0.681, and 0.713 for the test set and 0.679, 0.684 and 0.695 for the whole group, respectively (Fig. 4E). In the entire set, patients showed the same results for age, gender, T, N, or M staging, and grade (Fig. 4F).

Figure 4: Prognosis value of the 9 ARLs model in the train, test, and all sets. (A) Exhibition of ARLs model based on risk score of the train, test, and entire sets, respectively. (B) Survival time and survival status between low- and high-risk groups in the train, test, and all sets, respectively. (C) The heat map of 9 ARLs expressions in the train, test, and all sets, respectively. (D) Kaplan-Meier survival curve of OS for patients between risk groups in the train, test, and all sets, respectively. (E) The train, test, and all sets’ 1-, 3-, and 5-year ROC curves, respectively. (F) Kaplan-Meier survival analysis for HNSCC patients stratified by age, gender, T, N, or M staging, and grade in all set.

Independent prognostic ability of the signature

In univariate Cox regression analyses for all tumor samples, the hazard ratio (HR) for the risk score was 1.324, with a 95% confidence interval (CI) of 1.193–1.468 (p < 0.001). In addition, age, and stage were also statistically different (Fig. 5A). In the multivariate Cox regression, risk score (HR: 1.379, 95% CI: 1.225–1.553, p < 0.001), age (HR: 1.023, 95% CI: 1.007–1.038; p = 0.003) and stage (HR: 1.610, 95% CI: 1.322–1.960; p < 0.001) were independent prognostic factors of HNSCC (Fig. 5B).

Figure 5: Nomogram and assessment of the risk model. (A) Uni-Cox and (B) multi-Cox analyses of clinical factors and risk score with OS. (C) The nomogram that integrated the age, tumor stage, and risk score predicted the probability of the 1-, 3-, and 5-year OS. (D) The calibration curves for 1-, 3-, and 5-year OS. (E) The 3-year ROC curves of risk score, nomogram total score, and clinical characteristics.

Then based on risk score, age, and stage, we constructed a nomogram to predict the incidence rates of 1-, 3-, and 5-year OS in HNSCC patients (Fig. 5C). The calibration curves demonstrated a high concordance between the actual and nomogram-predicted survival rates of HNSCC patients (Fig. 5D). The sensitivity and specificity of the model for prognosis were evaluated via time-dependent ROC. The risk score (AUC = 0.684) and nomogram score (AUC = 0.715) had better predictive ability than age, gender, grade, and stage in the 3-year ROC of the risk model (Fig. 5E).

Gene set enrichment analyses and investigation of immune status

Firstly, with the GSEA software, we explored the KEGG pathways in low- and high-risk groups to investigate the biofunction differences. Of the top 10 significantly enriched pathways, the high-risk group was enriched in only one pathway related to steroid biosynthesis. In comparison, the low-risk group was enriched in four immune-associated KEGG pathways (all FDR < 0.25) (Fig. 6A). Then we analyzed the relevance of immune cells to the low- and high-risk groups using data from multiple platforms. Results showed the majority of immune cells were related to the low-risk group, which generally agrees with the GSEA results (Fig. 6B). The content of B cells, follicular helper T cells (cTfh), was negatively related to the risk score, suggesting they were more abundant in the low-risk group (Fig. 6C), which had larger immune score and ESTIMATE score (Fig. 6D). Furthermore, the low-risk group had more immune cell infiltration such as B cells, CD8+ T cells, and dendritic cells (DCs), and more prosperous immune functions (Fig. 6E). The above results imply that the low-risk group is under a higher immune infiltration state.

Figure 6: The investigation of tumor immune factors and immunotherapy. (A) GSEA of the top 10 pathways significantly enriched in the low- and high-risk groups. (B) The immune cell bubble of risk groups. (C) The correlation between risk score and immune cells. (D) The comparison of immune-related scores between low- and high-risk groups. (E) The difference between immune cell infiltration and immune functions in risk groups. *p < 0.05; **p < 0.01; ***p < 0.001.

Clinical treatment response analyses

Most immune checkpoints were also better activated in the low-risk group (Fig. 7A), which means we can group HNSCC patients by the risk pattern to select the appropriate checkpoint inhibitor. The low-risk group with higher immune scores had smaller IC50 in 15 immunotherapeutic agents (Fig. 7B) [15–29] and was more sensitive to 18 targeted therapeutics (e.g., Nilotinib) and six chemotherapeutic agents (e.g., Vinblastine) than the high-risk group. All 39 chemicals were shown in Fig. S1. Sixteen drugs had lower IC50 in the high-risk group (Fig. S2), of which three drugs including Bryostatin.1 were associated with immunotherapy (Fig. 7C) [30–32].

Figure 7: Immune checkpoints and drug sensitivity analyses. (A) The difference in immune checkpoint expression in risk groups. (B) Fifteen immunotherapeutic drugs showed lower IC50 values in the low-risk group. (C) Three immunotherapeutic drugs showed lower IC50 values in the high-risk group. *p < 0.05; **p < 0.01; ***p < 0.001.

LINC00861 inhibited the proliferation of nasopharyngeal carcinoma cells

We selected LINC00861 for further analysis since it had the highest correlation with the risk signature. The qRT-PCR indicated that LINC00861 was differentially expressed in three different nasopharyngeal carcinoma cell lines (Fig. 8A). We transfected the LINC00861-PCDNA3.1 plasmid in CNE1 and CNE2 cells to overexpress it and verified the transfection efficiency by qRT-PCR (Fig. 8B). CCK-8, Edu, and SA-β-gal Staining assays showed that LINC00861 overexpression considerably inhibited the growth and increased the senescence of CNE1 and CNE2 cells (Figs. 8C–8E).

Figure 8: Cellular experiments. (A) and (B) Expression level of LINC00861. (C) CCK-8, (D) Edu, and (E) SA-β-gal staining assay in CNE1 and CNE2 cell lines in different treatment groups. *p < 0.05; **p < 0.01; ***p < 0.001.

Discussion

Surgery, radiotherapy, and standard chemotherapy are commonly used for HNSCC. Recently, medicine including cetuximab, an anti-EGFR, and pembrolizumab, an anti-PD-1 therapy, was accepted for late-stage, recurrent, or metastatic HNSCC [33,34]. Regardless of therapeutic progress, the 5-year OS rate of HNSCC is still low (<50%) due to its ease of invasion, metastasis, development of chemo-resistance, and the fact that most cases are diagnosed at late stages [35]. Moreover, the 3-year OS rate is around 80% for HPV+ HNSCC and 55% for HPV- HNSCC [36]. Recently, a variety of novel biomarkers such as PITX2 methylation and C1GALT1 autoantibody that forecast prognosis, and B7-H3 that predicts immunotherapeutic prognosis and reaction in HNSCC were identified gradually [37–39]. However, the prognostic models and biomarkers need to be updated since more new drugs, including ICP inhibitors, are put into clinical practice. Here, a new model based on ARLs was constructed for the first time in HNSCC, which can both predict prognosis and differentiate the TIME.

Aging is closely connected with the occurrence of many chronic diseases and tumors [5]. Cell senescence is a physiological process whereby cells lose the proliferative ability forever. Cell senescence is triggered by various stimuli and the DNA damage response (DDR), leading to stimulation of the p53 and/or p16INK4A pathways, and its important regulatory factors include the cell cycle regulators p53, p21, p16, and pRb proteins [40]. The increased activity of lysosomal senescence-related β-galactosidase (SA-β-Gal) is an important marker of senescent cells [41]. Cell senescence is a double-edged sword in malignant tumors, as it prevents tumor development and promotes HNSCC progression through the secretion of various factors, called senescence-associated secretory phenotype (SASP). Notch, mTOR, and NF-κB pathways are engaged in SASP management [42]. The majority of studies conclude that cellular senescence acts primarily as a cancer suppressor rather than a cancer promoter [43–50]. Silencing peroxiredoxin1(Prx1)-induced cellular senescence blocks malignant conversion and inhibits the formation of oral squamous carcinoma [51]. In nasopharyngeal carcinoma, senescent cells promote tumor cell invasion by increasing the secretion of matrix metalloproteinase 9 (MMP 9) [52]. Bioinformatics analyses reveal the critical function of AGs in predicting the survival prognosis and revealing the TIME landscape of HNSCC [9,53].

Recently, many studies show that lncRNAs are critical in HNSCC occurrence, and can both promote and suppress cancer. More importantly, some lncRNAs can serve as biomarkers and treatment targets for HNSCC and can predict patient prognosis [2,54,55]. However, as the significant member involved in transcriptional and post-transcriptional regulation, the adjusting role of aging by lncRNAs in HNSCC is unknown. Therefore, we analytically validated the link of lncRNAs with aging and the prognosis of HNSCC, aiming to make some contribution to subsequent studies.

By LASSO and multi-cox regression analysis, we identified nine lncRNAs associated with eight AGs in HNSCC, all of which were positively regulated with the genes (Fig. 3E). Three lncRNAs (LINC01508, AL359502.2, AC243960.1) were related with a high risk of patients with poor prognosis, and the remaining 6 lncRNAs (AC246787.2, LINC00996, `CDKN2A-DT`, Z97653.1, LINC00861, LINC02384) were associated with low risk. These nine lncRNAs were enrolled into a prediction model. Through validation, the model was proved to be reliable and stable (Figs. 4, 5A–5E). By comparing with published prognostic models of HNSCC constructed based on ferroptosis-related lncRNAs, we found that these articles had high AUC predictive values for 3-year survival (0.83, 0.811, 0.726, and 0.687, respectively) [56–59], the 3-year AUC in our results was 0.779 for the train set, so collectively it seems that our proposed model shows a good accuracy in predicting prognosis. Reportedly, LINC01508 upexpression inhibited cisplatin tolerance of ovarian cancer cells by suppressing the Hippo-YAP pathway [60]. In cervical cancer cells, LINC00861 inhibits tumor progression and EMT by acting as a ceRNA for miR-513b-5p and modulating the PTEN/AKT/mTOR pathway [61]. The biological functions of the rest lncRNAs used in the modeling are unknown, about which researchers can carry out future experiments.

Aging not only regulates cell proliferation, but also alters the TIME through the secretion of SASP, which contains MMP, growth factors (VEGF and GM-CSF), inflammatory factors (IL-6 and IL-8), and the accumulation of lipofuscin deposits [62]. GSEA implied that the low-risk group was enriched with immune pathways such as the B cell receptor (BCR) and T cell receptor (TCR) pathways, and had higher immune infiltration and immune score, which may be related to their better prognosis (Figs. 6A and 6D). Considering the function of cell senescence in regulating TME, we conducted ssGSEA to explore the immune state in the low- and high-risk groups. The immune cells (aDCs, B cells, CD8+ T cells, DCs, iDCs, mast cells, neutrophils, NK cells, pDCs, T helper cells, Tfh cells, Th1 cells, Th2 cells, TIL, Tregs) and immune roles (APC (or T cell) co-suppression, APC (or T cell) co-stimulation, CCR, checkpoint, cytolytic activity, HLA, proinflammation, type II IFN response) were more active in the low-risk group (Figs. 6B and 6E). ICPs including PD-1, CTLA4, CD200, CD27, CD86, CD80, CD70, and TIGIT were more expressed in the low-risk group (Fig. 7A). These results imply that we can use the risk scores of ARLs to well distinguish the immune infiltration status of HNSCC and which is higher in the low-risk group. SASP factors released by senescent cells recruit innate immune cells (neutrophils, NK cells) and adaptive immune cells (CD8+ T cells) to modulate the removal of senescent tumor cells [63]. Therefore, the low-risk group may have more active cell senescence and thus is can be more cured by immunotherapy.

Radio-, chemo-, immuno-, and targeted therapies can induce cellular senescence and play an important role in TME through SASP [64]. Most HNSCC patients receive radiotherapy, which alters TME and activates the immune response of tumor cells [65,66] and may account for better outcomes in the low-risk group. As for drug sensitivity, the low-risk group was more sensitive to six inhibitors of the PI3K/AKT/mTOR pathway, such as rapamycin, which not only enhanced immunotherapeutic effects [18,22–25,29] but also blocked SASP-induced tumor progression (Fig. 7B). One study found that rapamycin prevented radiation-induced secretion of NFκB-driven pro-inflammatory SASP from prostate cancer cells and fibroblasts, thereby inhibiting tumor progression [67]. As a result, our study can guide clinical treatment and provide evidence for improving future immunotherapy and seeking appropriate target populations. Moreover, LINC00861 was identified as a protective factor against HNSCC progression in vitro (Fig. 8).

However, our study still has some limitations. Firstly, both our training and validation sets were from TCGA, which may cause bias, so the results may be more plausible if external validation was performed. Secondly, we only performed qRT-PCR, CCK-8, Edu, and SA-β-gal staining assays validation of LINC00861 in cell lines, but did not do further functional experiments or tissue validation, so further experiments are still needed. Finally, more clinical follow-up data shall be used to demonstrate the value of our prognostic model.

In conclusion, we built a robust prognostic prediction model for HNSCC using nine ARLs. It has a good predictive value in survival time, immune infiltration, immune checkpoint expression, and sensitivity to multiple drugs. LINC00861 acts as a protective factor against HNSCC progression.

Acknowledgement: We thank the supporters of Jiangsu Province Hospital Core Facility Centre in cellular experimental.

Funding Statement: This work was supported by the National Natural Science Foundation of China (82003228), the Natural Science Foundation of Jiangsu Province (BK20201080), and the Research Project of Clinical Medical Science and Technology Development Fund of Jiangsu University (JLY2021097).

Author Contributions: Xiaoke Di, Xinchen Sun, and Qing Gao were primarily involved in the design of this study. Xiaoke Di, Yuanyuan Sun, and Shu Zhou performed the bioinformatics analysis. Qing Gao and Yuanyuan Sun completed the cellular experimental validation. Qing Gao, Xiaoke Di, Yujing Shi, and Zeyuan Liu completed the manuscript writing. All authors read and approved the final manuscript.

Availability of Data and Materials: The public datasets to advocate our results can be gained from TCGA (https://portal.gdc.cancer.gov).

Ethics Approval: TCGA belongs to public databases. The patients involved in the database have obtained ethical approval. Users can download relevant data for free for research and publish relevant articles. Our study is based on open source data, so there are no ethical issues and other conflicts of interest.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

Figure S1: All 39 chemicals with lower IC50 values in the low-risk group.

Figure S2: All 16 chemicals with lower IC50 values in the high-risk group.