In this study, the RNA-seq data and clinical information of 528 patients with HNSCC were obtained from The Cancer Genome Atlas (TCGA). A total of 45 patients were excluded owing to missing survival time (n = 2), OS time < 30 days (n = 8), non-primary solid tumor (n = 30) or missing clinical variables (n = 5). Among the remaining 483 patients, CECT data of 211 patients were obtained within 2 weeks before surgery and corresponding clinical information were obtained from The Cancer Imaging Archive (TCIA). A total of 72 patients were excluded because of poor image quality (n = 60) and non-intersection with bioinformatic data (n = 12). After exclusion, 139 patients were enrolled in the study ( Figure 1 ). Genomic data for HNSCC and normal tissues were obtained from the University of California Santa Cruz Xena platform.

Kaplan–Meier (KM) plotting analysis was used to evaluate the prognostic value of IL1B in HNSCC using the “survival” and “survminer” packages in R (version 4.0.4). The cut-off value of IL1B expression in HNSCC tumor tissues was determined as 3.383 using the surv_cutpoint function of the “survminer” package in R. Accordingly, 483 patients were separated into two groups: those with high IL1B expression (n = 191) and those with low IL1B expression (n = 292). The log-rank test was used to estimate the survival probability, median survival times, and 95% confidence intervals (CIs) of the two groups. Because the KM curves of the IL1B high and low expression groups crossed, conditional landmark KM analyses of survival percentage were performed at 12, 24, 36, 48, 60, 72, 84 and 96 months to compare the survival curves stratified by IL1B expression. Landmark KM curves were plotted using the “jskm” package.

Through univariate Cox regression analysis, the prognostic value of the clinical and pathological features of HNSCC was identified. Since the listed features were proven to be common and important variables associated with OS in the clinic and in the literature, they were all included in the multivariate Cox analyses. Cox regression analyses were conducted with “survival” and “forestplot” packages in R.

Subsequently, subgroup analysis was utilized using “cmprsk”, “survival” and “forestplot” packages in R to evaluate the effect of IL1B in each subgroup. Interactions between IL1B and prespecified subgroups were evaluated by including age, gender, radiotherapy, chemotherapy, primary tumor site, HPV status, perineural invasion, T stage, N stage, M stage, and tumor grade in the Cox regression model.

By using the “clusterProfiler” package in R, we analyzed the enriched terms in the Gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) to get insight into the IL1B-related pathways. Differentially expressed genes were obtained using the “limma” package in R, with the threshold set at |log2FC| >1 and adjP < 0.05, by comparing the high and low IL1B expression groups. We employed the molecular function (MF), biological process (BP), and cellular component (CC) annotation databases for GO enrichment. The lollipop figures show the top 30 GO terms and top 30 KEGG pathways. Furthermore, 483 RNA-seq data were uploaded to the CIBERSORTx database to characterize the immune cell composition (27). Spearman’s rank correlation analysis was performed to examine the relationship between IL1B expression and immunocyte infiltration.

A professional radiologist performed region of interest (ROI) segmentation on all 139 CECT interpretations using the Medical Imaging Interaction Toolkit without referring to any clinical information. To avoid artificial and subjective deviations, another radiologist randomly selected 20 patients for the test and retest analysis and multiclinician segmentation. The radiomics features retrieved from the ROI segmentation were automatically extracted using PyRadiomics, an open-source program with 107 radiomics features, to assess the first-order statistical, morphological, and texture features. Table S3 displays all the radiomic feature categories. Prior to feature extraction, the data were standardized using Z-score normalization. The intraclass correlation coefficient (ICC) was used to measure the repeatability of both radiologists delineating the ROI. ICC > 0.8 represents high consistency, 0.51-0.79 represents medium consistency, and < 0.5 represents poor consistency. ICC < 0.8 features were omitted from further analysis.

A total of 139 patients were randomly assigned to the training and validation cohorts in a ratio of 6:4. The mRMR method was used to prioritize feature significance to reduce overfitting of the data and locate the best linked features. The mutual information (MI) to class labels was maximized, and MI with other characteristics was minimized to rank the input features. The RFE method was then applied to the top 20 features for high IL1B expression in the training cohort. During each iteration, the feature that contributed the least to enhancing the performance of the model was eliminated from the selection process. Feature ranking was determined using the number of iterations in which the feature was eliminated.

The features were then centered and scaled using the caret preprocess function and analyzed using the GBM method, which is a machine learning approach suitable for regression and classification issues. The “gbm” function in the “caret” package in R was used to produce a projected class probability score and construct a GBM model. Additionally, we employed a leave-one-out-based cross-validation approach to train and test the GBM model to minimize bias in feature selection. The GBM model was trained using all other samples in the dataset and was evaluated on the remaining sample that was not part of the training cycle. The top performing model was maintained and used. For each sample, the model generated a radiomics score (rad-score). The expression of IL1B was higher when the rad-score was higher.

The radiomics model was developed from the training cohort and estimated in the validation cohort. Receiver operating characteristic (ROC) curves were used to evaluate the diagnostic ability of the model. The evaluation indicators included the area under the ROC curve (AUC), accuracy (ACC), and specificity (SPEC). The validation cohort was independently validated to evaluate prediction accuracy. The calibration degree of the radiomics model was examined by performing the Hosmer–Lemeshow goodness-of-fit test and drawing the calibration curve. To evaluate whether the model is affected by data imbalances, the precision recall (PR) curve was generated using the “modEvA”package in R. To evaluate clinical utility, decision curve analysis (DCA) was performed using the "rmda" package in R. The Delong test was used to compare the diagnostic efficiency of each model.

The rad-score for each patient was computed using a radiomics model. Rad-scores were compared between the IL1B high- and low-expression group in the training and validation cohorts using the Wilcoxon test. A recent study demonstrated that IL1B pro-motes tumor metastatic outgrowth by affecting the epithelial-to-mesenchymal transition (EMT) of tumor cells (28). In this study, we focused on the potential association among rad-scores, IL1B expression, and EMT-related gene expression using Spearman correlation analysis with the “stats” package in R.

The “survminer” package in R was used to assess the rad-score and OS time of the 139 patients involved in the research. Using the greatest Youden index, the cut-off value for the rad-score was found to be 0.445. Accordingly, the patients were separated into two groups: those with high rad-scores (n = 69) and those with low rad-scores (n = 70). The KM curve was utilized to calculate survival probability and median survival time. The log-rank test was used to compare the survival probabilities of patients in the two groups.

SPSS (version 25.0) and R software (version 4.0.4) were used for all statistical analyses. Continuous and categorical variables are shown as medians. The optimal cut-off value of IL1B expression in HNSCC with continuous values was determined using the “survminer” package. The cut-off value of the rad-score was identified as the highest Youden index of the ROC curve. Categorical variables were compared using the Chi-square test or Fisher’s exact test. Continuous variables were compared utilizing the t-test or Mann–Whitney U-test. Statistical significance was defined as a two-tailed P-value < 0.05.

Four hundred and eighty-three patients with HNSCC from the TCGA database were included in this study. The cut-off value of IL1B expression was determined as 3.383 using the “survminer” package in R. Accordingly, the patients were separated into two groups: those with high IL1B expression (n = 191, IL1B expression > 3.383) and those with low IL1B expression (n = 292, IL1B expression < 3.383). Table 1 shows the baseline patient characteristics. There were no statistically significant differences in age, gender, chemotherapy, primary tumor site, perineural invasion, T stage, N stage, M stage, and tumor grade between the IL1B high- and low-expression groups (P > 0.05), except for HPV status (P = 0.012) and radiotherapy (P = 0.04).

IL1B was significantly overexpressed in HNSCC tissues than in healthy tissues (P < 0.001), and the median difference between the two groups was 1.338 ( Figure 2A ). According to the KM survival curves, there was a statistical critical level of difference between the IL1B high- and low-expression groups (P = 0.074, Figure 2B ). Owing to the crossovers on the KM curve, we conducted a landmark KM analysis. In the 0-36 months, the IL1B low expression group had a considerably better survival rate than the IL1B high expression group (P = 0.027, Figure 2C ). The survival rate from 48-78 months was not significantly different between the groups (P > 0.05). From 78-120 months in the late stage, the survival rate was higher in the IL1B low expression group (P = 0.044, Figure 2D ). In conclusion, survival analysis shows that high IL1B expression in patients with HNSCC in early (< 36 months) and late stages (> 78 months) indicates a poor survival probability.

We found that high IL1B expression was a risk factor for OS in univariate Cox analysis with a critical level of statistical significance (High vs. Low) (HR = 1.289, 95% CI:0.975-1.706, P = 0.075, Figure 3A ). After multivariate adjustment, high IL1B proved to be a risk factor for OS (High vs. Low) (HR = 1.56, 95%CI:1.161-2.096, P = 0.003, Figure 3B ). “HR = 1.56” indicates that high IL1B expression significantly increased the risk of death in patients with HNSCC. Differentially expressed IL1B exhibited a prognostic value in some clinicopathological subgroups ( Figure 3C ), including patients underwent radiotherapy (High vs. Low) (HR = 1.87, 95% CI: 1.19-2.937, P = 0.007), patients underwent chemotherapy (High vs. Low) (HR = 2.514, 95% CI: 1.533-4.122, P < 0.001), primary tumor in the larynx (High vs. Low) (HR = 1.845, 95% CI: 1.033-3.296, P = 0.038), tumors in N1/N2/N3/NX stages (High vs. Low) (HR = 1.47, 95% CI: 1.06-2.038, P = 0.021), and tumors in GE/G4/GX grade (High vs. Low) (HR = 1.795, 95% CI: 1.072-3.004, P = 0.026). The above analyses indicate that increased IL1B expression in patients with HNSCC predicts a potentially poor prognosis.

To investigate the potential mechanism underlying the role of IL1B in HNSCC, functional enrichment analysis was performed on DEGs between the IL1B high- and low-expression groups using the “clusterProfile” package in R. In the BP, CC, and MF categories, the most significantly enriched GO terms relative to the IL1B expression groups were response to xenobiotic stimulus, cell-substrate junction, and signaling receptor activator activity, respectively ( Figure 4A ). The most enriched KEGG pathways were herpes simplex virus 1 infection, pathways of neurodegeneration-multiple diseases, HPV infection, and MAPK signaling pathway, et al. ( Figure 4B ). Owing to the significance of IL1B, the association between immunocyte infiltration and IL1B expression was further investigated using the CIBERSORTx website, as shown in Figure 4C .

Among the 483 patients with RNA-Seq data in TCGA, 139 patients with corresponding CECT data were acquired from TCIA and used to develop a radiomics model for IL1B expression. The workflow is shown in Figure 5 . Two professional radiologists extracted 107 imaging features for each patient. The reproducibility of the radiomics features was assessed using the ICC ( Table S1 ). The mean ICC value of the imaging features was 0.941, and 97 radiomics features with an ICC value ≥ 0.80 (90.7% of all features) were included in in the follow-up feature selection.

A total of 139 patients were randomly assigned to the training cohort (n = 84) and validation cohort (n = 55) at a 6:4 ratio. Regarding clinical characteristics and IL1B expression, no significant differences were observed between the training and validation cohorts (P > 0.05 for all, Table S2 ). All 97 extracted features (ICC ≥ 0.80) were included in the mRMR algorithm, which selected 20 optimal features for further analysis ( Table S3 ). Then, the 20 features were imported into the RFE algorithm, and the top three features were adopted as the best features. The three features were shape_Sphericity, glszm_SmallAreaEmphasis, and firstorder_Kurtosis.

We trained these three features using the GBM algorithm and the leave-one-out-based cross-validation approach to obtain an ideal machine learning model for predicting IL1B expression in HNSCC tissues. Figure 6A outlines the importance of each feature in the final radiomics model ( Table S4 ). Figures 6B, C depicts the performance of the IL1B expression prediction. In the training cohort, the rad-score of the radiomics model was substantially greater in the group with high IL1B expression than that in the group with low IL1B expression (P = 1.404*10^-8, Figure 6D ). This trend was further confirmed in the validation cohort (P = 0.01; Figure 6E ).

The radiomics model performed well in the IL1B grading evaluation based on ROC curve analysis. The AUC of the radiomics model for the training group was 0.861 (95%CI:0.782-0.94, Figure 7A ), with ACC: 0.81, SPEC: 0.822. The AUC of the radiomics signature for the validation group was 0.703 ((95%CI:0.565-0.842, Figure 7B ), with ACC: 0.618 and SPEC: 0.667. The AUC of the PR curve in the training cohort was 0.803 ( Figure 7C ), whereas the AUC of the PR curve in the validation cohort was 0.627 ( Figure 7D ). The P-values of the H-L goodness of fit test were greater than 0.05 (training cohort: P = 0.393, Figure 7E ; validation cohort: P = 0.231, Figure 7F ), indicating a good agreement between the evaluated IL1B expression by the calibration curve of the model and the actual IL1B expression. Figure 6 depicts the radiomics model’s DCA curve. The DCA in the training cohort revealed that if patients’ high-risk threshold is between 15%-100%, utilizing the radiomics model to forecast IL1B expression offers greater value than “ALL” and “NONE” strategies ( Figure 7G ). The DCA curves for the validation cohort also demonstrated that the radiomics model adds greater benefit than “ALL” and “NONE” techniques in the 30%-85% high-risk threshold ( Figure 7H ), indicating high clinical utility. The Delong test demonstrated that the AUCs of the training cohort were not markedly different from those of the validation cohort (P > 0.05) in this study.

Because EMT is a common cellular phenotype that occurs only in tumors, we further explored the correlation among rad-score, IL1B expression, and EMT-related genes ( Table S5 ) expression in HNSCC tissue. The rad-score was significantly positively correlated with IL1B expression (cor = 0.472, P = 4.490*10^-9). The correlation between the rad-score and EMT-related genes was similar to that between IL1B expression and EMT-related genes. IL1B expression was positively and significantly correlated with CNN1, CXCL6, CXCL8, GADD45A, GLIPR1, INHBA, LAMA3, LAMC2, PLAUR, SERPINE1, THBS1, TNFAIP3, TNFRSF12A, and TPM4 (P < 0.05), and the rad-score was also positively and significantly correlated with these genes. IL1B expression was negatively and significantly correlated with ABI3BP, BGN, COMP, GPC1, ITGB5, MGP, SGCD, and SGCG (P < 0.05), whereas the rad-score had a similar relationship with EMT-related genes ( Figure 8 ).

The IL1B expression level was strongly associated with the prognosis of HNSCC and was well-predicted by the radiomics model. As a result, we used KM curves to investigate the association between the rad-score and OS in patients with HNSCC. The optimal rad-score cutoff value was determined as 0.445. Accordingly, the 139 patients were divided into a high rad-score group (n = 69) and a low rad-score group (n = 70). According to the log-rank test, the OS of the two groups was significantly different (P = 0.041, Figure 9 ); the median OS of the high and low rad-score groups was 61.267 and 96.667 months, respectively. The survival rates of the high and low rad-score groups at 65.7 months were 35.859% (95%CI:20.667%-62.194%) and 68.668% (95%CI:55.584%-84.832%), respectively. The corresponding survival rates at 96.5 months were 28.682% (95%CI:14.187%-57.985%) and 35.606% (95%CI:14.501%-87.424%), respectively. KM survival analysis indicated that a higher rad-score resulted in worse OS.