m6A‐ and immune‐related lncRNA signature confers robust predictive power for immune efficacy in lung squamous cell carcinoma

Immune checkpoint blockade has revolutionized immunotherapy of lung squamous cell carcinoma (LUSC), but the effective rate is less than 30% attributing to no effective and reliable method for predicting immune response. However, the clinical value of N6‐methyladenosine (m6A) and immune‐related lncRNA (mirlncRNA) concerning immune efficacy remains unknown. First, we identified a packet of specific mirlncRNA. Based on this, patients with LUSC were optimally divided into mirlncRNA clusters A, B, and C. The mirlncRNA cluster A was categorized as an immune‐inflamed phenotype distinguished by infiltration of numerous immune cells such as tumor‐infiltrating lymphocyte cells and highlighted by pathways such as regulation of myeloid leukocyte differentiation, while clusters B and C were found to correspond to immune‐desert and immune‐excluded phenotypes, respectively. Furthermore, the immune‐inflamed phenotype was shown to possess the highest immune infiltration, the lowest chromatin accessibility, survival rates, half inhibitory concentration (IC50), and the best immune efficacy. Finally, risk scores derived from the mirlncRNA signature helped identify subgroups of patients who could significantly benefit from immunotherapy. Encouragingly, in the population with poor response to the three targeted drugs, the immune response of patients with low drug sensitivity is significantly improved, indicating the vitality of combined therapy. The mirlncRNA signature not only identifies molecular typing and distinguishes chromatin accessibility, but also further highlights the immune efficacy and drug sensitivity, which might contribute to developing a new strategy for immunotherapy‐based individualized treatment.

drug sensitivity is significantly improved, indicating the vitality of combined therapy. The mirlncRNA signature not only identifies molecular typing and distinguishes chromatin accessibility, but also further highlights the immune efficacy and drug sensitivity, which might contribute to developing a new strategy for immunotherapy-based individualized treatment.

K E Y W O R D S
chromatin accessibility, immune efficacy, long-chain non-coding RNA, lung squamous cell carcinoma, N6-methyladenosine (m6A) methylation INTRODUCTION Lung cancer is the leading cause of cancer death in 2020. 1 Most lung cancers are non-small-cell lung cancer (NSCLC), 2 of which lung squamous cell carcinoma (LUSC) accounts for 20%-30%. 3 A number of LUSC patients have metastatic diseases, in which several distinct patterns of co-mutations are enriched at diagnosis, so that the overall survival (OS) rate remains poor. 4 Although outstanding progress has been made in treating LUSC using immunotherapies, and the combination of targeted therapy and immunotherap, 5 there are many problems as before. Targeted and immunotherapy therapies have changed the standard of care for several malignancies, and the prospects for combining them proved to be encouraging. 6 Due to the high tumor heterogeneity in LUSC patients, the outcomes of immunotherapy generally differ. 7 However, the molecular mechanisms underlying LUSC remain unclear. Therefore, it is crucial to identify potential prognostic biomarkers and therapeutic targets that may enhance immune efficacy in combating LUSC.
Previous studies have shown that long non-coding RNAs (lncRNAs) are promising therapeutic targets and innovative biomarkers for cancer. 8 Prognostic-related lncRNAs can be used to identify potential targets for immunotherapy. 9 Therefore, it is feasible to use lncR-NAs to evaluate LUSC. Research showed that a risk assessment model was established using five prognosticrelated lncRNAs for LUSC. 10 However, no analysis of the tumor microenvironment (TME) or immune checkpoints was conducted. Considering that growing evidence indicates that a molecular feature of TME has blossomed into a promising candidate during cancer formation and progression, 11 immunotherapy targeting immune checkpoints (ICs) has made significant strides in improving the survival of patients. 12 We believe that further research is required.
Another study established that the relationship between non-coding RNAs and N6-methyladenosine (m6A) modifications provides a new direction for investigating the potential regulatory mechanisms of gene expression in malignant tumors. 13 m6A methylation is orchestrated by demethylases, binding proteins, and methyltransferases, also designated as "erasers," "readers," and "writers," and refers to estimate TME infiltration characterization to provide an insight into effectual immunotherapy for cancer. 14 Nevertheless, the full influence of m6A regulators on lncRNA dysregulation remains unclear. Consequently, analyzing the relationship between lncRNAs modified by m6A and the development of LUSC may help to identify biomarkers that may allow for precise molecular subtyping to guide immune therapy. Research has indicated that biomarkers derived from two specific gene sets can help improve the accuracy of the model. 15 Some researchers have found that the expression profiles of immune genes can be used to reveal the underlying structures of the immune landscape within tumors, the intrinsic properties of which may be incorporated when developing biomarkers. 16 Therefore, we combined m6A regulators, immune genes, and lncRNAs to build a better signature.
In this study, we first developed an m6A-related lncR-NAs (mirlncRNA) signature based on a group of specific lncRNAs associated with m6A regulators and differentially expressed immune genes (DEI genes), and then conducted unsupervised clustering to classify LUSC patients into distinct clusters. Subsequently, we analyzed chromatin accessibility, estimated immune infiltration, and biological function, and explored immune efficacy and drug sensitivity. Based on this, we further evaluated the correspondence between the different clusters and the three recognized tumor immune subtypes. Finally, we fitted the risk score used for risk stratification and construction of the nomogram model to predict the survival probability.

2.2
Identification of mirlncRNA signature for LUSC and unsupervised clustering mrlncRNAs were obtained using correlation analysis between 23 m6A regulators and 14,080 lncRNAs with the following parameters: R = .4, p = .001. Univariate Cox regression was used to determine the six mrlncRNAs associated with prognosis in LUSC.
Unsupervised cluster analysis was applied to identify different clusters according to the enrichment scores of the 29 immune cell subtypes, and the patients were classified for subsequent processing. The above steps were performed using the "ConsensuClusterPlus" package, and 50 repetitions were performed to ensure the stability of the classification. Then, 183 DEI genes were screened between the clusters using the "limma" package. The criteria for screening DEI genes were logFC = 1 and FDR = 0.05. m6A-and immune-related lncRNAs (mirlncRNAs) were identified according to the correlation analysis of six mrlncRNAs and 183 DEI genes with the following parameters: R = .35, p = .001. Different mirlncRNA clusters emerged according to the expression of the five mirlncRNAs by unsupervised clustering analysis.

Analyses of ATAC-seq (assay for targeting accessible chromatin with high-throughput sequencing) data and somatic mutations
We obtained ATAC-seq data for LUSC from the TCGA database. 17 The "CHIPseeker" package was used to conduct the annotation of 3 kb of data upstream and downstream from the peak center, and the ATAC-seq data were classified into three groups according to the mirlncRNA clusters. Heatmaps and average profiles were also generated using the "plotHeatmap" and "plotProfile" scripts in DeepTools. 18 The somatic mutation data of 494 LUSC patients from distinct mirlncRNA clusters were ana-lyzed, summarized, and visualized using the "maftools" package.

Immune cell infiltration analysis
The single-sample gene-set enrichment analysis (ssGSEA) algorithm was used to quantify the relative abundance of 29 immune signatures in LUSC. 19 Immune cell abundance identifier (ImmuCellAI), an online tool that can estimate tumor immune infiltration and predict immunotherapy response powerfully and uniquely, 20 was used to evaluate the abundance of 24 immune cells. Estimation of stromal and immune cells in malignant tumors using expression data (ESTIMATE) is a method for deducing the proportion of immune and stromal cells in tumor samples using gene expression data. 21 Using ESTIMATE, stromal, immune, and ESTIMATE scores were calculated to reveal dissimilar immune infiltrating levels in each LUSC patient.

Analysis of therapeutic efficacy
The expression of immune checkpoints (ICs) was analyzed to explore differences in immune efficacy. The IPS has been scientifically proven to predict patient responses to immunotherapy. 22 Therefore, the IPS of patients with LUSC was analyzed using the "limma" package to predict the efficacy of immunotherapy. The R package "pRRopheticto" was utilized to investigate different drug responses among distinct mirlncRNA clusters. Half inhibitory centration (IC50) was used to measure drug sensitivity. 2D structures of the drugs were downloaded from https://pubchem.ncbi.nlm.nih.gov/. The 2D structures were transformed into 3D structures using Chem3D software.

Functional enrichment analysis
Using the R package "limma," differentially expressed mRNAs (DEmRNAs) were identified between every two mirlncRNA clusters. The criteria for screening DEImRNAs were logFC = 1 and FDR = 0.01. Kyoko Encyclopaedia of Genes and Genomes (KEGG) pathway enrichment and gene ontology (GO) analysis were executed by using packages "org.Hs.eg.db," "ggplot2," and "enrichplot." Gene set enrichment analysis (GSEA) is a calculation method used to determine whether a transcendentally defined set of genes is statistically significant, accordant differences between two biological states, 23 which was applied between every two mirlncRNA clusters, processed by GSEA software (version 4.2.1).

Construction of prognosis risk model
Optimal prognostic mirlncRNAs were scr eened by the least absolute shrinkage and selection operator (LASSO) method using the R package "glmnet." 24,25 Subsequently, the risk score was calculated using the following formula: Here, n is the number of mirlncRNAs, exp i is the expression level of lncRNA i , and a i is the regression coefficient of lncRNA i .

Statistical analysis
R (version 4.4.1) and GSEA (version 4.2.1) were used in our study. Correlation analysis was conducted using Spearman analysis, and the landscapes of the co-expression networks were plotted using the R package "igraph." Survival curves were constructed using the Kaplan-Meier method, and the log-rank test was used to identify the significance of differences. Heatmaps were created using the R package "pheatmap." Statistically significant p-values were defined as values less than .05. 26

Construction of mirlncRNA signature for LUSC
To explore the relationship between 14,080 lncRNAs and 23 m6A regulators, a co-expression network was established. Nineteen m6A regulators and 323 lncRNAs were incorporated according to the following parameters: R = .4, p = .001 ( Figure 1A). A univariate Cox regression analysis was performed to select six mrlncRNAs that were significantly associated with prognosis (p = .05, Figure 1B). Unsupervised clustering of the 494 LUSC samples was then conducted using the "ConsensusClusterPlus" package. Finally, two distinct clusters were identified, which were termed as m6A clusters 1 and 2, k = 2 chosen as the beat clustering number in accordance with the heatmap, the cumulative distribution function value, and delta area ( Figure 1C and Figure S1A,B). Significant differences were observed in the transcriptional profile of mrlncRNAs between the two dissimilar m6A clusters ( Figure 1D). Patients with m6A cluster 2 (n = 125) had a better prognosis, while patients of m6A cluster 1 (369 patients) had poorer prognosis ( Figure 1E). The expression of the six mrlncRNAs in normal adjacent tissues was significantly lower than that in LUSC tissues ( Figure S1C).
To identify the mirlncRNA signature, common differentially expressed immune genes (DEI genes) were identified. For the DEI genes being selected well and truly, unsupervised clustering analysis, based on the enrichment scores of 29 immune signatures calculated by ssGSEA, was performed to classify patients into three immuno-related subtypes according to the heatmap and delta area, named clusters A, B, and C ( Figure S2A,B). Based on the heatmap of 29 immune signatures and three immunotypingrelated biomarkers (CD274, TGFβ-1, and stromal cells), three clusters were found to have significantly distinct immune infiltration characteristics ( Figure 2A). To exhibit the results more intuitively, a different comparison of immune cell infiltration was performed. The expression of TGFβ-1 in cluster A, being in close proximity to that in cluster C, was the highest and cluster B the lowest ( Figure 2B,C). As shown in Figure S2C, cluster A was highly enriched in stromal activation pathways, including the TGF-β signaling pathway. Cluster C expressed enrichment pathways involved in fully activated immunity, including cytokine-cytokine receptor interactions ( Figure S2D).
The differences in immune genes among three clusters were analyzed using the "limma" package, three groups of DEI genes were screened out between cluster A, B, and C, respectively (logFC = 1, FDR = 0.05). Then, 183 common DEI genes were identified after intersection analysis ( Figure 2D). Finally, the co-expression relationship between the six mirlncRNAs and 183 DEI genes is shown in Figure 2E. Five mirlncRNAs and eight DEI genes were incorporated according to the following parameters: R = .35, p = .001. The five mirlncRNAs can be used for the clinical outcome of patients with LUSC, which were named as mirlncRNA signatures and were fitted as risk scores according to LASSO analysis and Cox regression.
Unsupervised clustering of the 494 tumor samples was conducted based on the expression of the mirl-ncRNA signature. Three was selected as the optimal clustering number based on the heatmap and delta area ( Figure 3A,B). These three clusters were named mirl-ncRNA clusters A, B, and C. Prognostic analysis of the three mirlncRNA clusters revealed that the particularly prominent survival disadvantage was in mirlncRNA cluster A, which had the lowest expression of mirlncRNAs ( Figure 3C). The risk score decreased from one cluster to the next in the following order: mirlncRNA clusters A, C,  Figure 3D). It can be found that the mirlncRNA transcriptional profiles were significantly different among normal adjacent tissues and three diverse mirlncRNA clusters ( Figure 3E). The distribution of riskscore in the patient population was exhibited in Figure 3F. Attribute changes in individual patients were visualized using the alluvial diagram ( Figure 3G). Similar results were obtained in the validation cohort ( Figure S3A-F).

Analyses of chromatin accessibility and tumor mutation between three mirlncRNA clusters
To further characterize the mirlncRNA signature, genomewide chromatin accessibility of samples from three mirl-ncRNA clusters was investigated using ATAC-seq. The ATAC-seq signal in mirlncRNA cluster A was lower than in the other two clusters ( Figure 4A). The scaled reads increased from one cluster to the next, in this order: mirl-ncRNA cluster A, C, and B ( Figure 4B-D, p < .001). The tumor mutational burden (TMB) was the highest in mirl-ncRNA cluster B and the lowest in mirlncRNA cluster A ( Figure 4E). The "maftools" package was used to analyze LUSC patients' somatic mutation profiles. TP53, TTN, and CSMD3 were the three most common mutated genes, and TP53 was the top one ( Figure 4F-H). It was worth noting that the higher the expression of the mirlncRNA signature, the higher the TMB.

Analyses of immune infiltration, immune efficacy, and drug sensitivity between three mirlncRNA clusters
To investigate the immune infiltration characteristics in distinct mirlncRNA clusters, the enrichment scores of 29 immune signatures were calculated using the ssGSEA algorithm. Surprisingly, in the training cohort, differential analysis showed that immune infiltration of all immune signatures was the highest in mirlncRNA cluster A and the lowest in mirlncRNA cluster B ( Figure 5A). However, patients with mirlncRNA cluster A did not show a matching survival advantage and had the worst prognosis. The same analysis was performed on 24 cells, whose infiltrations were calculated using the ImmuCellAI algorithm. In the training cohort, the infiltration of most ImmuCellAI cells was highest in mirlncRNA cluster A, second highest in mirlncRNA cluster C, and lowest in mirlncRNA cluster B ( Figure 5B). In addition, the ESTIMATE, stromal, and immune scores were the highest in mirlncRNA cluster A and the lowest in mirlncRNA cluster B according to ESTIMATE analysis in the training cohort ( Figure 5C-E).
Most results were consistent with those of the training cohort, except for two ImmuCellAI cells (CD4 naïve cells and effector memory cells) with no significant dissimilarities in the validation cohort ( Figure S4A-E). Spearman's correlation analysis showed that the enrichment scores of the 29 immune signatures were positively correlated with risk scores ( Figure S5A).
To appraise the predictive ability of the mirlncRNA signature for immunotherapy efficacy, we first tracked the correlation between the mirlncRNA signature and ICs. Then, a negative correlation was found, indicating that the higher the expression of the mirlncRNA signature, the lower the expression of ICs ( Figure S5B,C). As shown in Figure 6A-E, the expression of HAVCR2, CD48, TIGIT, CTLA4, and LGALS9 was significantly higher than that of the other two clusters. Similar results for ICs were found in the validation cohort ( Figure S6A-E).
Immunotherapy using checkpoint-blocking antibodies targeting CTLA-4 and PD-1/PD-L1 has emerged as a promising approach for the treatment of various malignancies. 27 Thus, the IPS file of patients was used to forecast the clinical reaction to immunotherapy. Intriguingly, in the training cohort, we found that patients in mirlncRNA cluster A with the highest risk score and the worst prognosis had the best immunotherapy response when using anti-PD-1 monotherapy (nivolumab, p < 0.01, Figure 6F), anti-CTLA-4 monotherapy (ipilimumab, p < .05, Figure 6G), and combined immunotherapy (p < 0.001, Figure 6H). The results of patients receiving combined immunotherapy were verified in the validation cohort (p < .001, Figure S6F).
To analyze the responses to chemotherapy among the three mirlncRNA clusters further and comprehensively, the "pRRophetic" package was used to estimate the chemotherapeutic response. In the training cohort, we found that patients in the mirlncRNA cluster A had the lowest IC50 of the six chemotherapeutic agents, with the highest risk score and worst prognosis ( Figure 6I,L,M,P,Q,T). The 3D structure tomographs of these drugs are shown in Figure 6J,K,N,O,R,S. Similar results were obtained for the validation cohort ( Figure S6G-T). Risk score between low-and high-risk groups were significantly distinct ( Figure S7A-B), and    people with higher risk score survived worse ( Figure  S7C) and had higher IC 50 of AZD3759, BI-2536 and Gefitinib ( Figure S7D,G,J). These patients with LUSC were divided into low-and high-IC 50 groups according to the median of patients' IC50 of each agent. In low-risk group, patients with higher IC 50 of BI-2536 had higher IPS, while AZD3759 and Gefitinib not ( Figure S7E,H,K). In high-risk group, those with higher IC 50 of every agent had better immunotherapy response ( Figure S7F,I,L).

Functional annotation between three mirlncRNA clusters
To further analyze the benefits of the mirlncRNA signature, functional enrichment analysis was performed using KEGG pathway enrichment and GO analysis. In total, three groups of DEmRNAs were identified between mirlncRNA clusters A, B, and C (FDR = 0.01, logFC = 1).
As expected, DEmRNAs between mirlncRNA clusters A, B, and C were enriched in functional modules linked with immunity such as immune receptor activity, respectively (Figure 7A, D, G). Furthermore, signaling pathways enriched by KEGG analysis based on DEmRNAs between mirlncRNA clusters A, B, and C are exhibited in Figure 7B, E, H, respectively.
GSEA was performed to explore the biological functions of the three mirlncRNA clusters. As expected, when mirlncRNA cluster B was used as the control group, mirl-ncRNA cluster A was significantly enriched in biological processes related to regulatory T cells, myeloid suppressors, and cytokines, such as negative regulation of myeloid leukocyte differentiation ( Figure 7C). When mirlncRNA cluster C was used as the control group, mirlncRNA cluster A was significantly enriched in the regulation of myeloid leukocyte-mediated immunity, and so forth ( Figure 7F). When mirlncRNA cluster B was used as the control group, mirlncRNA cluster C showed enrichment in biological processes, such as vascular endothelial growth factor production ( Figure 7I). However, mirlncRNA cluster B lacked immune biological processes ( Figure 7C,I).

DISCUSSION
In the present study, we identified a packet of specific lncR-NAs associated with m6A regulators and DEI genes. Based on the expression of the mirlncRNA signature, patients with LUSC were optimally divided into mirlncRNA clusters A, B, and C. The mirlncRNA cluster A was categorized as an immune-inflamed phenotype distinguished by infiltration of numerous immune cells such as TILs and highlighted by pathways such as regulation of myeloid leuko-cyte differentiation, while clusters B and C were found to correspond to immune-desert and immune-excluded phenotypes, respectively. Furthermore, the immune-inflamed phenotype (cluster A) was shown to possess the highest immune infiltration, the lowest ATAC-seq signal, survival rates, IC50, and TMB, as well as the best immune efficacy. Finally, quantified risk scores derived from the mirlncRNA signature helped identify the subgroups of patients with LUSC who could significantly benefit from immunotherapy. Taken together, the mirlncRNA signature was shown to contribute to the identification of immune subsets and prediction of immune response in patients with LUSC ( Figure 8 and Figure S7). (Figure 8 and Figure S8). Currently, many dysregulated lncRNAs have been detected, such as JPX and PKMYT1AR, which are associated with tumorigenesis and tumor progression. 28 The expression of lncRNAs is associated with the occurrence, metastasis, and prognosis of LUSC. 29 Inspired by the view that biomarkers derived from two specific gene sets improve the accuracy of the model, 15 we combined m6A regulators, immune genes, and lncRNAs to improve the accuracy of the prognostic signature. After LASSO analysis and Cox regression analysis, we identified five mirlncR-NAs (SNHG21, AP001469.3, AL122125.1, AC138035.1, and AL096701.3) as prognostic signatures (Figures 1 and 2). In our study, all five mirlncRNAs were upregulated in LUSC tissue ( Figure S1C), and the lowest expression patients had the highest risk score ( Figure S3D) and the poorest prognosis ( Figure 3C). Eventually, these five mirlncRNAs were regarded as protective factors (HR < 1) whose lowest expression was associated with the worst survival. Therefore, our mirlncRNA signature represents a great innovation for predicting clinical outcomes.
Open chromatin-accessible regions have been found to contain cis-regulatory elements and possibly modulate gene activity. 30 A study has pointed out that in NSCLC cells, the transcriptional activity of migration-associated genes is downregulated using the Feiyanning formula treatment, which prolongs patient survival. 31 Therefore, we hypothesized that low ATAC-seq signal is associated with a better immune response. In our study, patients in mirlncRNA cluster A had the lowest ATAC-seq signal ( Figure 4A-D) and the best immunotherapeutic response ( Figure 6F-H).
TMB is considered a representation, in a tumor exome, of non-synonymous coding mutations that have been proven to differ between patients. 32 In our study, TP53, TTN, and CSMD3 were the three genes with the top three mutation rates in LUSC patients ( Figure 4F-H), which were related to the occurrence, treatment, and prognosis of lung squamous cell carcinoma and other diseases. 33 Previous studies concluded that patients with high TMB were associated with increased sensitivity and clinical benefits of anti-PD-1/PD-L1 therapy. 34 However, in our study, patients with the lowest TMB had the best immune efficacy. Some scholars have suggested that high TMB does not guarantee good immune efficacy from immune checkpoint inhibitor (ICI) therapy 35 and TMB may not be invariably related to ICI reactivity. 36 Our findings confirmed these conclusions. Moreover, the definition of TMB cutoffs for patient stratification remains uncertain. Foun-dationOne calls high TMB when greater than or equal to 20 mutations/Mb are detected; TMB-Intermediate: 6-19 mutations/Mb; and TMB-Low: ≤5 mutations/Mb. 37 The mean TMB values of the three mirlncRNA clusters in our study were all below six ( Figure 4E). If the criteria mentioned above were followed, none of the clusters in our study would have high TMB. Therefore, although the TMB level was statistically different among the three mirl-ncRNA clusters, the difference between the TMB of the three mirlncRNA clusters may not be enough to cause a difference in immune efficacy. This may be why some scholars believe that TMB needs to be further refined before it can be broadly adopted across tumor types and patient populations. 38 Based on the mirlncRNA signature, we revealed three distinct clusters that presented significantly dissimilar TME cell infiltration characteristics ( Figure 5A,B and Figure S4A,B) and were in accordance with the three phenotypes mentioned above. The first is the immunedesert phenotype, which can reportedly be the result of immunological ignorance, the induction of tolerance, or a lack of appropriate T-cell priming or activation. 39 In our study, mirlncRNA cluster B was categorized as an immune-desert phenotype owing to the lack of signaling pathways related to immunology (Figures 7C  and 8I) and had the lowest immune infiltration. Regarding immune-excluded phenotypes, some studies have pointed out that immune-excluded tumors may be associated with particular stromal-based inhibition, specific vascular factors or barriers, or a particular chemokine state. 40 Therefore, mirlncRNA cluster C was categorized as an immune-excluded phenotype due to extracellular F I G U R E 8 Schematic diagram. We identified a packet of specific lncRNAs associated with m6A regulators and DEI genes. Based on the expression of the mirlncRNA signature, patients with LUSC were optimally divided into mirlncRNA clusters A, B, and C. The mirlncRNA cluster A was categorized as an immune-inflamed phenotype distinguished by infiltration of numerous immune cells such as TILs and highlighted by pathways such as regulation of myeloid leukocyte differentiation, while clusters B and C were found to correspond to immune-desert and immune-excluded phenotypes, respectively. Furthermore, the immune-inflamed phenotype (cluster A) was shown to matrix, chemokines, and vascular factors [70], such as positive regulation of chemokine production, regulation of extracellular matrix disassembly, and vascular endothelial growth factor production, among others ( Figure 7I). And last but not least was the immune-inflamed phenotype that presented a colossal amount of immune cell infiltration in TME, 41 including suppressor B cells, cancer-associated fibroblasts, myeloid-derived suppressor cells, as well as immune-inhibitory regulatory T cells. 40 Therefore, mirl-ncRNA cluster A was categorized as an immune-inflamed phenotype distinguished by the infiltration of numerous immune cells such as TILs, B cells, and immune-inhibitory regulatory T cells, and distinguished by pathways such as regulation of myeloid leukocyte differentiation and positive regulation of regulatory T-cell differentiation ( Figure 7C,F). Some studies showed patients characterized by immune-inflamed phenotype had the best prognosis, 14 while the others did not. 42 In this study, mirlncRNA cluster A, classified as an immune-inflamed phenotype, had the poorest prognosis. The possible causes of this result are analyzed below. First, considering the upregulation of ICs in patients in mirlncRNA cluster A, we inferred that ICs mediated immune escape so that the antitumor effect of immune cells was inhibited. Second, numerous studies showed that some immune cells are involved in tumor progression in certain situations such as neutrophils, 43 immunosuppressive M2 macrophages, 44 natural killer cells, 45 T helper 17 cells, 46 T cells, and natural killer T cells. 47 Furthermore, crosstalk between proximal immune cells and cancer cells eventually gives rise to an environment that is a major contributor to tumor growth and metastasis. 47 In addition to the above conclusions, considering that the mechanisms that dictate the balance between cancer-promoting and cancer-inhibitory inflammation within TME remain unclear, we believe that current research on TME is not thorough enough. Therefore, it was not surprising that mirlncRNA cluster A, the immuneinflamed phenotype, had the greatest immune cell infiltration, but the poorest prognosis ( Figures 3C and 5). In our study, mirlncRNA cluster A, which was categorized as an immune-inflamed phenotype, had the highest immune infiltration, and the expression of ICs, such as CTLA-4, was also the highest in mirlncRNA cluster A ( Figure 6A-E). Likewise, another study suggested that the upregulation of ICs was positively correlated with high lymphocyte infiltration. 48 In other words, cancer-immune phenotypes classified using our mirlncRNA signature may be able to predict IC expression. Furthermore, it has also been proven that the expression of CTLA-4 was significantly elevated in the cohort that showed clinical benefit compared with the one that showed no clinical benefit, indicating that the expression of ICs may correlate with immunotherapy response. 49 That is, patients in mirlncRNA cluster A with the highest immune infiltration and ICs expression should respond well to immunotherapy. Interestingly, our research showed this result ( Figure 6F-H). In conclusion, our mirlncRNA signature provides predictive advantages for precision immunotherapy for LUSC.
Although surgical resection is the capital treatment for non-metastatic lung cancers, it is only suitable for a small proportion of patients with NSCLC and is limited by the number and location of the lesions. 50 Patients who are not eligible for surgery universally need chemotherapy to enhance their quality of life and prolong their survival time. 51 We found that patients in mirlncRNA cluster A, with the lowest IC50, were the most sensitive to some potential chemotherapeutic agents ( Figure 6I-T). Furthermore, compared with traditional chemotherapy, ICI treatment has resulted in long-term survival in some patients with advanced NSCLC. 52 In accordance with the tumor immunoediting hypothesis, the immune system can promote or select tumor cells with reduced immunogenicity, thereby providing developing tumors with a mechanism to escape immune responses that kill tumors. 53 This may upregulate the expression of immunosuppressive molecules in tumors, such as CTLA-4 and PD-1. 54 Considering that the expression of ICs was significantly different in the three mirlncRNA clusters ( Figure 6A-E), we speculated that different groups of patients would respond differently to immunotherapy. We detected that patients in mirlncRNA cluster A were most likely to respond to anti-CTLA-4 or anti-PD-1 therapy ( Figure 6F-H). Therefore, we believe that the prognosis of these patients may improve after immunosuppressant or more specific chemotherapy drugs are used. We found that patients with high-risk score responded significantly less well to targeted therapy drugs. Those patients with high IC 50 of drugs responded better to immunotherapy while those with low IC 50 not ( Figure S7). In other words, for those high-risk patients with poor prognosis, although the efficacy of targeted therapy drugs was unsatisfactory, the efficacy of immunotherapy was positive, which meant the risk model may contribute to the clinical treatment strategy of the combination of targeted therapy drugs and immunotherapy agents. [55][56][57] possess the highest immune infiltration, the lowest ATAC-seq signal, survival rates, IC50, and TMB, as well as the best immune efficacy. Finally, quantified risk scores derived from the mirlncRNA signature helped identify the subgroups of patients with LUSC who could significantly benefit from immunotherapy. Taken together, the mirlncRNA signature was shown to contribute to the identification of immune subsets and prediction of immune response in patients with LUSC.
Despite our positive findings, this research has some limitations. First, the results of the study were obtained using bioinformatics methods based on TCGA database. Although some findings have been verified internally, the relevant results still need to be verified experimentally. Second, the data used to assess the efficacy of immunotherapy were acquired from TCGA dataset as opposed to an actual immunotherapy cohort. Third, a larger sample size is required to enhance the reliability of immune signatures.
In summary, the mirlncRNA signature based on five miraculous lncRNAs can not only distinguish molecular typing and identify the corresponding tumor immune subtypes and their chromatin accessibility but also further highlight the immune efficacy and drug sensitivity of tumor immune subtypes. The mirlncRNA signature may provide new ideas and methods for identifying cancer immunotyping and may contribute to developing a new strategy for immunotherapy-based individualized treatment of patients with LUSC.

A U T H O R C O N T R I B U T I O N S
Zihan Ding and Haijian Zhang have access to all the data in this study and take responsibility for the integrity and accuracy of the data analyses. Haijian Zhang, Yang Liu, Changhong Cheng, and Yi Han: study concept and design. Zihan Ding and Qinhua Huang drafted the manuscript. Zihan Ding, Lining Song and Chen Zhang: statistical analysis. Haijian Zhang, Xiaohong Cui, and Yan Wang supervised the study. All authors have read and approved the final manuscript.

A C K N O W L E D G M E N T S
The authors acknowledge The Cancer Genome Atlas (TCGA) database for the convenience of this research. All of the study sponsors had no role in the study design, collection, analysis, and interpretation of data.

C O N F L I C T O F I N T E R E S T S TAT E M E N T
The authors declare that there is no conflicts of interest.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The datasets used in the current study and related scripts of bioinformatics analysis are available from the corresponding author upon reasonable request. The link to the code is as follows: https://github.com/hjz005/ZihanDing20220813 O R C I D Haijian Zhang https://orcid.org/0000-0003-2516-4215