Samples from 425 Chinese patients with lung cancer between January 2019 and June 2020 were retrospectively obtained. The inclusion criteria were as follow: 1) lung cancer, 2) performed with targeted panel sequencing and immunohistochemistry (IHC) assay for PD-L1. The clinical information of the patients is shown in Supplementary Table S1. Genomic and transcriptome data of 969 NSCLC patients from TCGA were used to analyze the difference in gene expression, enriched pathway and cell proportions across the different statuses of HLA LOH. (Thorsson et al., 2018).

Genomic and clinical data from 89 NSCLC patients treated with ICIs were obtained to analyze the predictive efficacy of HLA LOH, TMB, PD-L1 and CD8⁺ T cells on immunotherapy, and the comprehensive predictive model was constructed based on these data (Anagnostou et al., 2020). The classification criteria in the original article was used to define patients as TMB high and low. For PD-L1, the proportion of PD-L1+ cells higher than 1% was defined as PD-L1 high, and the others was defined as PD-L1 low. The cut-off values for CD8⁺ T high and CD8⁺ T low were defined as the median values of CD8⁺ T cells.

The Dako PD-L1 IHC 22C3 pharmDx assay was used to detect PD-L1 protein expression in formalin-fixed paraffin-embedded (FFPE) tumor tissue slides according to the manufacturer’s recommendations (Kulangara et al., 2019; Lantuejoul et al., 2020). Briefly, the FFPE slides were heated at 65°C for 30 min, and then dewaxed, rehydrated, and fixed. Next, the slides were incubated with anti-PD-L1 antibody (clone 22C3). After three thorough washes, the slides were incubated with horseradish peroxidase (HRP)-conjugated secondary antibody. After incubation, the tyramide signal amplification (TSA) was added to the slides at the dilution of 1:100. Finally, the slides were stained with 4′-6′-diamidino-2-phenylindole (DAPI) at the dilution of 1:10. Images were captured with a Vectra 3.0 pathology imaging system microscope (PerkinElmer Inc.). A three-tiered grading system was applied to evaluate the proportion of PD-L1 expression in tumor cells: “negative” (Tumor Proportion Score (TPS) < 1%), “intermediate” (1% ≤ TPS <50%), and “high” (TPS ≥50%).

Genomic DNA was isolated from tumor biopsies and matched peripheral-blood samples using the GeneRead DNA FFPE Kit (Qiagen, 180,134). All sample capture libraries were prepared using the YuceOne Plus v2.2 (YuceBio, Shenzhen, China). The target capture specific probe is hybridized with genomic DNA to enrich the DNA fragments of the target genomic region and then sequenced using high-throughput sequencing technology. SOAPnuke (v1.5.6) was used to filter the original reads to remove low-quality reads with the unknown base “N” greater than 10% (Chen et al., 2018). BWA (v0.7.12) was used to align the clean reads with the human reference genome (UCSC GRCh37/hg19) (Li and Durbin, 2009). Somatic mutations were detected with the VarScan (Version 2.4) (Koboldt et al., 2012). Possible false-positive mutations were filtered using Bcftools (1.14) with the parameter set as follow: “basicfilter = """' (STRLEN (REF) > 50 || STRLEN (ALT) > 50) || INFO/STATUS!∼"Somatic"'""" hotspotfilter = """'INFO/HOTSPOT! = "." && ((INFO/SOR! = 0 && INFO/SOR<3) || INFO/VD < 5 || INFO/AF<0.007 || INFO/SSF>0.05)'""" fpdbfilter = """'INFO/HOTSPOT = "." && ((INFO/FPDB! = "0" && INFO/FPDB! = ".") || (INFO/GERMLINE! = "0" && INFO/GERMLINE! = "."))'""" normalfilter = """'INFO/HOTSPOT = "." && ((INFO/GERMLINE)! = "." || (FORMAT/PMEAN [0]<20)||((INFO/SOR! = 0 && INFO/SOR<5) || INFO/AF<0.02 || INFO/SSF>0.01)||(INFO/AF<0.05 && FORMAT/MQ [0]<50)||(FORMAT/MQ [0]<30)||(INFO/AF<0.05 && FORMAT/QUAL [0]<30) || ((INFO/MSI>10||(INFO/MSILEN>1 && INFO/MSI>4)) && INFO/AF<0.3)||(type! = "snp” && INFO/MSI>3 && ((INFO/MSILEN=(strlen (REF)-1))||(INFO/MSILEN=(strlen (ALT [0])-1))) && INFO/AF<0.1) || (FORMAT/NM [0]>2 && FORMAT/MQ [0]<60 && INFO/AF<0.2) || (FORMAT/NM [0]>3 && (FORMAT/MQ [0]<55||FORMAT/NM (Hellmann et al., 2018a)>3)) || (FORMAT/DP [0]<30 || FORMAT/DP (Hellmann et al., 2018a)<30)|| INFO/VD < 10 || (FORMAT/BIAS [0:0] = "2" && FORMAT/BIAS [0:1] = "1") || (FORMAT/SBF [0] < 0.05 && FORMAT/VD [0]<50) || ((INFO/SOR! = 0 && INFO/SOR<10) && FORMAT/MQ [0]<60))' """“. Identified somatic mutations were annotated with SnpEff (Version 4.3) (Cingolani et al., 2012). Somatic copy number alterations (SCNAs) were analyzed using Allele-Specific Copy number Analysis of Tumors (ASCAT) (v3.1.0) with default parameters (Van Loo et al., 2010).

HLA-I typing of tumors and matched normal samples was performed as previously described (Yi et al., 2021). In short, the patched opitype software was used for HLA typing (Szolek et al., 2014). Then loss of heterozygosity in human leukocyte antigen (LOHHLA) algorithm was performed to identify HLA LOH in tumor samples (McGranahan et al., 2017). It is classified as LOH if the following two conditions are met: 1) A copy number <0.5 and 2) allelic imbalance is determined with p < 0.01 using the paired Student’s t-Test between the two distributions. According to the proportion of LOH in HLA alleles, LOH status was further divided into three grades: negative (all HLA alleles had no LOH), low (LOH/HLA alleles are 1/6, 2/6, or 1/5), high (LOH/HLA alleles are 1/3, 1/4, ≥3/6 or ≥2/5).

TMB was determined as the number of non-synonymous mutations per megabase of the genome examined. All nonsynonymous mutations and their upstream and downstream nucleotide sequences of 10 amino acids were selected and translated into 21-mer peptides. The translated peptide was used to produce 8- to 11-mer peptides with a sliding window approach. The binding affinity of peptide and major histocompatibility complex (MHC) class I was predicted by NetMHCpan (version 3.0) (Hoof et al., 2009). If the predicted half-maximum inhibitory concentration (IC50) binding affinity was not greater than 500 nM, the peptide was selected. Multiple selected peptides generated by the same mutation were counted as a neoantigen. Tumor neoantigen burden (TNB) was measured as the number of such peptides per megabase of the genome examined. The chromosomal instability (CIN) score was used to estimate copy number burden, which was calculated as follows: the ploidy of the sample was generated by the ASCAT algorithm (v3.1.0). For each of the 22 autosomal chromosomes, the percentage of gained and lost genomic material was calculated relative to the ploidy of the sample. The CIN score of a sample was defined as the average of this percentage value over the 22 autosomal chromosomes (Braun et al., 2020). Gene mutation in pathways analysis was compared with the previously reported gene list (Sanchez-Vega et al., 2018).

The normalized counts matrix was used to identify differentially expressed genes (DEGs) with edgeR package (Chen et al., 2016). DEGs with |log2FoldChange| > 1 and a p-value < 0.05 were considered as significant DEGs. The Matplotlib (v3.4.3) package was used to generate volcano plots. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, Canonical Pathways, WikiPathways, GO Biological Processes and Reactome Gene Sets enrichment were performed with Metascape web tool with a p-value < 0.05 (Zhou et al., 2019). The gene expression matrix was used to estimate the level of 64 cells with xCell (Aran et al., 2017).

The decision-tree was evaluated to construct the comprehensive predictive model with biomarkers of CD8⁺ T cells, TMB, PD-L1 and HLA LOH to predict the efficacy of ICIs in NSCLC patients. The predictive efficacy of the model was analyzed with the leave-one-out method. Different depth of the tree, from 3 to 20, were tested to calculate the F1 score of the comprehensive predictive model. The most optimal depth of decision tree was determined based on the highest calculated F1 score.

Correlations between HLA LOH and clinical characteristics were analyzed using the Fisher’s exact test for categorical variables. Kruskal–Wallis rank-sum tests were used for comparisons of continuous variables across multiple groups. The t-test was used to compare differences between groups. Receiver operating characteristic (ROC) curve analysis was performed to measure the discriminatory ability of HLA LOH, TMB, PD-L1 and CD8⁺ T cells. Kaplan-Meier curve, Log-rank test, and cox-regression were used to assess the association of HLA LOH, TMB, PD-L1 and CD8⁺ T cells with patient survival. p < 0.05 was considered significant. Statistical analyses were performed in the R (v3.6.1) and Python (v3.8.8). The packages in Python were used in the study including lifelines (v0.27.0) and scipy (v1.7.3).

To examine tumor genomic characteristics and their potential association with different HLA LOH statuses, samples from 425 Chinese lung cancer patients were retrospectively collected and subjected to targeted panel sequencing and PD-L1 expression detection. According to the different statuses of HLA class I alleles, which are encoded by three genes (HLA-A, HLA-B, and HLA-C), 425 patients were divided into three groups: HLA LOH negative group, HLA LOH low group and HLA LOH high group. Compared with the HLA LOH negative and HLA LOH low groups, there was a significant gender discrepancy in the HLA LOH high group [68% (Male) versus 32% (Female)] (Supplementary Table S2).

The genomic landscapes with different statuses of HLA LOH were shown in Figure 1A. The top 10 genes were listed in order of mutated frequency: EGFR, TP53, KRAS, PIK3CA, KEAP1, ZFHX4, RBM10, FAT1, CDKN2A, and CTNNB1 (Figure 1A), which had a slight difference with the TCGA database (Cancer Genome Atlas Research, 2014). In the TCGA dataset, the top five genes with the highest mutation frequency were TP53 (46%), KRAS (33%), KEAP1 (17%), SKT11 (17%) and EGFR (14%). Interestingly, a high frequency of BCL6 amplification was found in the HLA LOH subgroup (the proportion in the HLA LOH negative, HLA LOH low and HLA LOH high groups were 4%, 17%, and 15%, respectively). BCL6 is related to the stress response in breast cancer, lung cancer, glioma, and other solid tumors. Overexpression of BCL6 is associated with tumor immune surveillance and drug resistance during the process of chemotherapy and radiotherapy (Fernando et al., 2019).

Among 425 lung cancer patients, the proportions of HLA LOH negative, HLA LOH low and HLA LOH high were 76.24% (324/425), 9.88% (42/425) and 13.88% (59/425), respectively. (Figures 1B,C and Supplementary Table S1). The frequency of HLA LOH in lung adenocarcinomas (LUAD) was lower than in lung squamous cell carcinomas (LUSC) (Figure 1D), which was consistent with previous studies (Zhao et al., 2021).

Furthermore, the association of genomic features and PD-L1 expression with the different statuses of HLA LOH were analyzed. It was found that TMB was significantly higher in the HLA LOH high group than in the other groups (p < 0.0001) (Figure 1E, Supplementary Table S2). However, there was no significant difference in TNB among these groups (Figure 1F, Supplementary Table S2). In addition, the CIN score was significantly different between the HLA LOH high and HLA LOH negative groups (p = 0.0033; Figure 1G). Consistent with the previous reports, the HLA LOH was related to PD-L1 expression among the three groups (p = 0.07642) (McGranahan et al., 2017; Montesion et al., 2021) (Figure 1H, Supplementary Table S2).

To further investigate the association of gene alteration with different statuses of HLA LOH in Chinese lung cancer patients, a comparison analysis between HLA LOH high and HLA LOH negative groups was conducted. It was found that the mutation frequency of 25 genes was significantly higher in the HLA LOH high group than in the HLA LOH negative group (Figure 2A), such as the TP53 mutations (85.7% versus 50.9%; p < 0.0001)), FAT1 mutations (16.1% versus 6.5%; p = 0.032), RB1 mutations (16.1% versus 5.2%; p = 0.0098), TRRAP mutations (14.3% versus 4%; p = 0.0026), ERBB2 mutations (12.5% versus 4.9%; p = 0.019) and BCL6 amplification (12.5% versus 4%; p = 0.0026) (Figure 2B). Conversely, EGFR mutation was associated with HLA LOH negative group when compared with HLA LOH high group (50.3% versus 37.5%, p = 0.0677). Other lung cancer-related genes, such as APC, KRAS, and BRAF, did not show significant differences among the groups (Figure 2B). Furthermore, the mutation in oncogenic signaling pathways were analyzed. It was found that several pathways were significantly different among the groups, such as the TP53 pathway, cell cycle pathway, Hippo pathway, WNT pathway and the RTK-RAS pathway (Figure 2C).

To explore the TME characteristics of patients with different HLA LOH statuses, genomic and transcriptional data of 969 NSCLC patients were obtained from TCGA. The HLA LOH was defined as the copy number of HLA class I or II-related genes less than 2. The differentially expressed genes (DEGs) between HLA LOH and HLA LOH negative groups were analyzed. There were 1081 genes significantly upregulated and 497 genes significantly downregulated in the HLA LOH groups when compared with HLA LOH negative groups (Figure 3A). Downregulated genes are associated with immunity, which suggests the pre-existing immune recognition of the tumor in HLA LOH negative patients (Figure 3B). To further verify this phenomenon, abundance of immune cells was evaluated by xCell. It was found that the proportion of B cells, CD8⁺ T cells, DC and M1 macrophages as well as the immune score and microenvironment score were significantly higher in HLA LOH negative patients when compared with HLA LOH patients. (Figure 3C). These results suggested that a pre-existing immunologically active TME may be related to HLA LOH negative patients.

To investigate the predictive efficacy of HLA LOH on the treatment of ICIs, genomic and clinical data of 89 NSCLC patients treated with ICIs were obtained (Supplementary Table S3). As shown in Figures 4A,B, the median overall survival (mOS) and median progression-free survival (mPFS) of patients with HLA LOH were 17 and 3 months, respectively. For patients with HLA LOH negative, the mOS was not yet reached, and the mPFS was 11 months. Patients with HLA LOH negative had a longer survival period than those with HLA LOH, which indicated that HLA LOH has a negative effect on the prognosis and may represent a mechanism of tolerance to immunotherapy. Compared with HLA LOH, it was found that TMB, PD-L1 and CD8⁺ T cell were insufficient to identify responders (Supplementary Figures S1A–F). Furthermore, multivariate cox regression was analyzed. And the result showed that only HLA LOH was significantly associated with the prognosis (Figure 4C).

TMB and PD-L1 expression were higher in HLA LOH high group (Figures 1E–H). To further investigate the prognostic association of HLA LOH with TMB and PD-L1 expression, additional stratification of HLA LOH for TMB or PD-L1 was analyzed. As shown in Figures 4D–G, the combination of HLA LOH with TMB or PD-L1 could further distinguish responders.

Compared with LUAD, the frequency of HLA LOH in LUSC was higher, which suggested that HLA LOH may have different effects on the efficacy of ICIs in LUAD and LUSC (Zhao et al., 2021). Therefore, the predictive effect of HLA LOH on immunotherapy in NSCLC with different pathological classifications was further analyzed. As shown in Figures 4H–K, HLA LOH was significantly correlated with the survival period in LUSC patients.

To further evaluate the predictive ability of HLA LOH for ICIs response, the ROC was plotted. As shown in Figure 5A, the area under the curve (AUC) value of HLA LOH was 0.42, indicating that HLA LOH was insufficient to distinguish responders from non-responders who were subjected to immunotherapy. Then other reported biomarkers, such as TMB, PD-L1 expression and CD8⁺ T cells were analyzed; however, none of the AUC values for these biomarkers exceeded 0.7, which suggested that a single biomarker was not effective enough to precisely distinguish patients who would benefit from the treatment of immunotherapy. In order to further improve the predictive efficacy, a comprehensive predictive model containing multiple factors related to anti-tumor immunity (HLA LOH, TMB, PD-L1 and CD8⁺ T cells) was constructed. As shown in Supplementary Figures S2A,B, the model had the highest F1 score when the decision tree depth was three. Therefore, the parameter was used to construct the comprehensive predictive model. Furthermore, the responders predicted by the comprehensive predictive model (pDCB) had a significantly longer survival period than non-responders predicted by the comprehensive predictive model (pNDB) (Figures 5B,C), and the predictive ability was better than the single biomarker (Figures 4A,B, Supplementary Figures S1A–F).