Identification of immune‐related lncRNA panel for predicting immune checkpoint blockade and prognosis in head and neck squamous cell carcinoma

Abstract Purpose Immunotherapy is changing head and neck squamous cell carcinoma (HNSCC) treatment pattern. According to the Chinese Society of Clinical Oncology (CSCO) guidelines, immunotherapy has been deemed as first‐line recommendation for recurrent/metastatic HNSCC, marking that advanced HNSCC has officially entered the era of immunotherapy. Long non‐coding RNAs (lncRNAs) impact every step of cancer immunity. Therefore, reliable immune‐lncRNAs able to accurately predict the immune landscape and survival of HNSCC are crucial to clinical management. Methods In the current study, we downloaded the transcriptomic and clinical data of HNSCC from The Cancer Genome Atlas and identified differentially expressed immune‐related lncRNAs (DEir‐lncRNAs). Further then, Cox and least absolute shrinkage and selection operator (LASSO) regression analyses were performed to identify proper DEir‐lncRNAs to construct optimal risk model. Low‐risk and high‐risk groups were classified based on the optimal cut‐off value generated by the areas under curve for receiver operating characteristic curves (AUC), and Kaplan–Meier survival curves were utilized to validate the prediction model. We then evaluated the model based on the clinical factors, immune cell infiltration, and chemotherapeutic and immunotherapeutic efficacy between two groups. Results In our study, we identified 256 Deir‐lncRNAs in HNSCC. A total of 18 Deir‐lncRNA pairs (consisting of 35 Deir‐lncRNAs) were used to construct a risk model significantly associated with survival of HNSCC. Cox proportional hazard regression analysis confirmed that our risk model could be served as an independent prognostic indicator. Besides, HNSCC patients with low‐risk score significantly enriched of CD8+ T cell, and corelated with high chemosensitivity and immunotherapeutic sensitivity. Conclusion Our risk model could be served as a promising clinical prediction indicator, effective discoursing of the immune cell infiltration of HNSCC patients, and distinguishing patients who could benefit from chemotherapy and immunotherapy.


| INTRODUC TI ON
Head and neck squamous cell carcinoma (HNSCC) mainly arises in the oral cavity, oropharynx, hypopharynx, and larynx. HNSCC was the seventh most common cancer worldwide. 1 Multiple pathogenic factors (including heavy use of tobacco, alcohol, and chronic infection of human papilloma virus) have been proven in association with the initial of HNSCC. 2,3 Smoking could increase the risk of HNSCC from 5-to 25-fold, and the risk increases with the quantity of cigarette. 3,4 The consumption of alcohol also independently doubles the risk of HNSCC. 5 Cancer is characterized by the accumulation of various of genetic alterations, which results in the acquisition of ten biological capabilities during the multistep development of human tumors. 6 Subsequent molecular studies have revealed that these genetic or epigenetic alterations often result in mutated cellular antibodies or biomarkers and aberrantly expressed genes. [7][8][9] These abnormal expressed molecules on the surface of cancer cells could be recognized by CD8 + T cell, distinguishing them from their normal counterparts and even inhibiting the outgrowth of cancer cells. 10 Therefore, in the past decades, researchers were dedicated to efficiently activate cancer immunity system to fight cancer cells, which is called "immunotherapy." However, the cancer immunity cycle does not perform optimally in cancer patient; sometimes, rarely protective efforts, even promoting assists by T cell were found during the progression of human cancers. [11][12][13] Either tumor antibodies may not be detectable, or dendritic cells and T cells may classify these antigens as self thereby creating low ratio of cytotoxic T lymphocytes. [14][15][16] Besides, negative regulators to T cell responses are also responsible for the failure of immune protection. 17 The goal of cancer immunotherapy is to initiate or reinstate the cycle of cancer immunity to exert anti-tumor effects. The cancer immunotherapies must be carefully configured to overcome these negative feedback mechanisms. Indeed, immunotherapy exhibits marvelous efficacy in clinical trials 18 ; however, only a minority of patients have incredible response to immunotherapy in clinical practice. Therefore, a deep understanding of immune system and identification of these sorts of HNSCC patients will be helpful to develop effective therapeutic strategies for cancer immunotherapy and to improve the overall survival of HNSCC. In the present study, we aim to comprehensively explore immune-related lncRNA expression profile and construct a model for predicting outcome, drug sensitivity, and adapting implement clinical strategies in HNSCC.

| Data download and pretreatment
The online database TCGA program was aimed to molecularly characterize human cancers. In our study, a series of transcriptome RNA-

| Immune-related differential expressed lncRNAs acquisition
Firstly, immune-related genes were obtained from the ImmPort database (http://www.immpo rt.org) and was used to screen immunerelated lncRNAs (ir-lncRNAs) by a co-expression strategy. Then, correlation analysis was performed between all lncRNAs of HNSCC samples and known immune-related genes by Pearson's correlation.
The cut-off criteria of adjusted p-value (adj. p-value) was set as 0.001 and the criterion of correlation coefficients was set as more than 0.4. Subsequently, we applied limma package of R software to identify the differentially expressed ir-lncRNAs. The thresholds were set as log fold change (FC) >1 along with false discovery rate (FDR) <0.05. Besides, to analyze the function of ir-lncRNAs, gene ontology (GO) was performed using the online database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david.ncifc rf.gov/, version 6.8).

| Pairing DEir-lncRNAs
The Cancer Genome Atlas (TCGA) is a comprehensive database that includes multi-layered cancer genome profiles. Large-scale collection of data inevitably generates batch effects introduced by K E Y W O R D S drug sensitivity, head and neck squamous cell carcinoma, immune-related long non-coding RNAs, immunity cell infiltration, immunotherapy differences in processing at various stages from sample collection to data generation. To eliminate these factors, we cyclically paired these differentially expressed ir-lncRNAs (DEir-lncRNAs). We constructed a new 0-or-1 matrix by comparison of all the DEir-lncRNAs expression in each patient of TCGA. We defined the DEir-lncRNA pair as 1 if the expression level of lncRNA A is higher than lncRNA B; otherwise, this pair was defined as 0. This new matrix was used as further analyses.

| Survival-related DEir-lncRNAs pair
Paired DEir-lncRNAs associated with clinical outcomes in HNSCC patients were identified as survival-related DEir-lncRNA pairs.
Survival-related DEir-lncRNA pairs were selected by univariate Cox regression applying survival package of R software. These survivalrelated Deir-lncRNA pairs were specified for subsequent research. were obtained. Then, the highest point of risk score of maximum AUC were set to the cut-off criteria to classify HNSCC patients into high or low risk of group.

| Validation of the constructed risk model
To validate the constructed risk model, we performed Kaplan-Meier analysis to show the survival difference of high and low risk group. Subsequently, we also re-assessed the relationship between the model and clinicopathological characteristics. Then, to confirm whether the model can be used as an independent clinical prognostic biomarker, we performed the univariate and multivariate Cox regression analyses, and the results were represented by forest map.

| Immune cell abundance analysis
We collected the currently acknowledged methods to calculate the immune infiltration statues from CIBERSORT. 19 The results were shown in a lollipop diagram, which was finished by violet package of R software.

| Prediction of clinical application of the risk model
Recent studies have predicted the response to immune checkpoint inhibitors (ICIs) of cancers based on their genomic features, and the immunophenoscore (IPS) was correlated with responses to ICIs and heatmap (1B) plots of DEir-lncRNAs of HNSCC immunotherapy. 20,21 The score of IPS-CTLA-4 and IPS-PD1 stands for the potential for response rates of ICIs in HNSCC. Besides, we calculated a half inhibitory centration (IC 50 ) of common chemotherapeutic drugs, such as cisplatin, docetaxel, gemcitabine, and paclitaxel. The difference in the IC50 between the high-and low-risk groups were assessed, and the results were shown as box drawings by pRRophetic package of R. 22 Then, we also analyzed the relationship between the model and the expression of ICI-related genes.
Subsequently, the immune-related mRNAs (ir-mRNAs) of TCGA were achieved by the intersection of mRNA matrix of TCGA and a list of well-known immune-related genes from ImmPort database.
Then, the immune-related lncRNAs (ir-lncRNAs) were obtained by Pearson's analysis between lncRNA matrix and ir-mRNAs of HNSCC from TCGA. Eventually, we obtained a total of 804 ir-lncRNAs.
Subsequently, we identified differentially expressed ir-lncRNAs (DEir-lncRNAs) between 501 HNSCC tumor tissues and 44 normal tissues, the results showing 256 DEir-lncRNAs with FDR-adjusted p < 0.05 and absolute log fold change (FC) >1 (Table S1). These DEir-lncRNAs were represented by heatmap and volcano plots in Figure 1. The functions of lncRNAs are thought to reflected by their associated mRNAs. Therefore, to exploit the underlying functions of immune-related lncRNA signature, gene ontology (GO) analyses was carried out for differentially co-expressed mRNAs. The results of the GO analysis are shown in Figure S1; many immunity-associated processes were identified, including immune response, innate immune response, adaptive immune response, and T and B cell receptor signaling pathway.

| Identification of survival-related DEir-lncRNAs pair
We cyclically paired DEir-lncRNAs pairs through comparison of all

| Validation of the risk model
According to the cut-off point confirmed above, 195 patients were grouped into high-risk group and 304 into low-risk group. Kaplan-Meier analysis showed that patients in the low-risk group exhibited a longer survival time than those in the high-risk group (p < 0.0001) ( Figure 4). Subsequently, we performed stratified analysis between high-and low-risk group by clinical characteristics (such as gender, age, TNM stage, and grade), and subgroup analysis results robustly supported that better overall survival of patients in low-risk group than in high-risk group ( Figure S2). In order to avoid basis from clinical features, univariate and multivariate Cox proportional hazard regression analysis were performed to evaluate the independent prognostic value of risk model. The results showed that our risk model, clinical stage, and age could be served as independent prognostic indicators ( Figure 5).
Then, we assessed the relationship between the risk score of HNSCC and clinicopathological characteristics (TNM stage, age, grade). The results showed that patients with advanced T stage more likely gathered in the high-risk group, indicating that higher scores of the 18 DEir-lncRNA pairs might be significantly associated with the progression of HNSCC ( Figure S3).

| Immune cell abundance analysis based on risk model
We analyzed the infiltration difference of immune cells between the low-and the high-risk group classified by the risk model. The immune infiltration estimations for TCGA expression profiles was calculated by CIBERSORT (https://ciber sort.stanf ord.edu/). 23 Then, the intersection of the infiltration estimation file and our risk score file were acquired for further Pearson's correlation analysis. The association between immune cell infiltration difference and different risk group was presented by vioplot ( Figure 6). Patients in lowrisk group showed higher infiltration levels of B cell, CD8 T cells,

| Prediction drug sensitivity based on risk model
Considering the positive association between two important ICI biomarkers, we assess the innate sensitivity or resistance to anti-PD-1 therapy. Our results showed that HNSCC patients with low-risk score represented significantly higher IPS-CTLA4 and IPS-PD1 (immunophenoscore, IPS), indicating the potential of ICI application for HNSCC ( Figure 8). Besides, we also identified associations between risk and the efficacy of common chemotherapeutics (including cisplatin, docetaxel, gemcitabine, and paclitaxel) in HNSCC (Figure 9). The results indicated that HNSCC patients with low risk represented higher sensitivity to docetaxel and gemcitabine.

| DISCUSS ION
HBSCCs are a heterogeneous group of cancers originated from oral cavity, pharynx, and larynx, which are characterized by different molecular subgroups and clinical features. The primary treatment of early stage of HNSCC is radical radiotherapy or surgery, and the treatment option of advanced stage of HNSCC is radiotherapy combination with chemotherapy. However, a substantial proportion of HNSCC patients die from their disease after such treatment.
Besides, these extensive surgery and high-dose chemotherapies are also mutilating for HNSCC patients. We analyzed the expression profile of lncRNA in HNSCC for the following reasons: (1) the amount of ncRNAs in the whole genome is 98%, which were once thought to be "garbage sequence";  Table S1) were consist of an amount of well-studied lncRNAs associated with human cancers, such as LINC01614 in breast cancer, lung cancer, and gastric cancer. 30,31 Changes in the composition, physical properties, and spatial con-  A stroma-related lncRNA signature (SLS) has been reported to predict adjuvant chemotherapy benefit in patients with colon cancer. 43 LncRNAs are also participated in the acquired resistance to chemotherapy, 44,45 and targeting to these lncRNAs can reverse drug resistance and enhance the sensitivity of cancer cells to chemotherapy. 46 Several studies have indicated that lncRNAs not only are involved in the typical hallmarks of cancer but also are closely correlated with the regulation of immune response. 47,48 Given the crucial role of lncRNAs in cancer, we also assessed the sensitivity to ICIs, such as pembrolizumab and nivolumab; these patients represented favor outcome, which rationalize the above phenomenon. And our prediction of drug sensitivity also showed that HNSCC patients with low-risk score represented significantly higher sensitivity to PD-1 and CTLA4.
In summary, we identified an immune-related lncRNAs in HNSCC, and constructed a risk model for prediction prognosis for patients with HNSCC and might help in distinguishing patients who could benefit from anti-tumor immunotherapy.

AUTH O R CO NTR I B UTI O N S
All authors contributed to the planning and design of the study. CZ had full access to all of the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.

CO N FLI C T O F I NTE R E S T
The authors declare that there are no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The datasets analyzed during the current study are available in the TCGA database (https://portal.gdc.cancer.gov/).