Immune-related genes and gene sets for predicting the response to anti-programmed death 1 therapy in patients with primary or metastatic non-small cell lung cancer.

Although antibodies targeting the immune checkpoint protein programmed death-1 (PD-1) exert therapeutic effects in patients with primary or metastatic non-small cell lung cancer (NSCLC), the majority of patients exhibit partial or complete resistance to anti-PD1 treatment. Thus, the aim of the present study was to identify reliable biomarkers for predicting the response to anti-PD-1 therapy. The present study analyzed tumor specimens isolated from 24 patients (13 with primary and 11 with metastatic NSCLC) prior to treatment with approved PD1-targeting antibodies. The expression profile of 395 immune-related genes was examined using RNA immune-oncology panel sequencing. The results demonstrated that six immune-related differently expressed genes (DEGs), including HLA-F-AS1, NCF1, RORC, DMBT1, KLRF1 and IL-18, and five DEGs, including HLA-A, HLA-DPA1, TNFSF18, IFI6 and PTK7, may be used as single biomarkers for predicting the efficacy of anti-PD-1 treatment in patients with primary and with metastatic NSCLC, respectively. In addition, two DEG sets comprising either six (HLA-F-AS1, NCF1, RORC, DMBT1, KLRF and IL-18) or two (HLA-A and TNFSF18) DEGs as potential combination biomarkers for predicting the efficacy of anti-PD-1 therapy in patients with NSCLC. Patients with a calculated expression level of the DEG sets >6.501 (primary NSCLC) or >6.741 (metastatic NSCLC) may benefit from the anti-PD-1 therapy. Overall, these findings provided a basis for the identification of additional biomarkers for predicting the response to anti-PD-1 treatment. Copyright: © Chen et al.

Introduction

According to global cancer statistics in 2018, lung cancer is the most commonly diagnosed cancer and the leading cause of cancer mortality, accounting for 18.4% of total cancer deaths (1). Non-small cell lung cancer (NSCLC) accounts for 80–85% of all cases of lung cancer (2). The optimal therapy for improving survival rates of NSCLC patients depends on an early diagnosis at stage I and II when surgical resection remains a feasible and consistent option (3). However, surgery is not an effective therapy for patients with an advanced or metastatic NSCLC; instead, chemotherapy, radiotherapy and targeted therapy (e.g. targeting epidermal growth factor receptor and vascular endothelial growth factor) can be used alone or in combination (4,5). Monoclonal antibodies blocking immunological checkpoints, also termed immune checkpoint inhibitors (ICIs), have become promising therapeutic options for patients diagnosed with advanced NSCLC with or without metastasis (6). Among them, anti-programmed death-1 (PD-1) monoclonal antibodies have been widely used for treating advanced NSCLC (7–10). Antibodies targeting PD-1 prevent tumor cells from escaping immune-mediated destruction (11).

The PD-1 receptor is a vital immune checkpoint molecule expressed on activated T cells that mediates immunosuppression (12). Binding of PD-1 to its ligands (PD-L1 and PD-L2) on cancer cells suppresses T cells, resulting in evasion of the immune response (12). Anti-PD-1 antibodies bind to PD-1 receptors to disrupt the inhibition of T cells by tumor cells, enhancing the antitumor effects of the immune system (12). Previous clinical trials, including CheckMate 017 (13) and 057 (14), KEYNOTE-010(9) and KN-024 (10), have demonstrated that compared with docetaxel, anti-PD-1 antibodies significantly prolonged the overall survival and had a favorable safety profile in patients with advanced NSCLC. Despite the prominent therapeutic effects of anti-PD-1 monoclonal antibodies on patients with advanced NSCLC, only a limited fraction of patients benefit from this immunotherapeutic agent (15). Therefore, developing reliable methods to predict the efficacy of anti-PD-1 monoclonal antibody in treating patients with advanced NSCLC may provide economic relief for non-responding patients and save time to start other types of therapy.

At present, identification of patients with NSCLC for anti-PD-1 therapy is mainly based on the visual assessment of PD-1 expression levels in tumor tissue specimens by immunohistochemistry (IHC) (16). However, PD-1 IHC data interpretation is subjective and may be inconsistent due to a varied response to anti-PD-1 therapy in patients (17). Although a number of potential indicators for predicting response to anti-PD-1 therapy have been identified, including mutational burden (18), gut microbiome (19) and tumor-infiltrating lymphocyte abundance (20), these indicators lack sufficient sensitivity and specificity (21). Based on next-generation sequencing, gene expression profiling allows simultaneous assessment of a large number of genes and has been well used to developing response signatures for a number of types of tumor, including NSCLC (22). The sensitivity of RNA immune-oncology (IO) panel sequencing is >20-fold higher compared with that of whole transcriptome sequencing; specifically, RNA IO panel sequencing can detect lowly expressed coding genes with high repetition and identify differential expression with >2-fold changes (23). Thus, RNA IO panel sequencing provides a sensitive and accurate approach for identifying biomarkers and feature genes in clinical studies (24).

Previous studies have reported that the expression levels of PD-L1 are markedly different between primary tumors and nodal metastases in patients with advanced NSCLC (25,26). In addition, the effects of anti-PD-1 antibodies differ between patients with primary and metastatic NSCLC (8). The present study aimed to investigate the roles of robust biomarkers in determining whether patients with primary or metastatic advanced NSCLC may benefit from anti-PD-1 antibody treatment.

Materials and methods

Study design

The present study was an observational study, which was conducted in accordance with the Declaration of Helsinki and approved by the Hunan Provincial Tumor Hospital (Changsha, China; approval no. 2019 fast review of scientific research [08]). All patients enrolled in the present study met the following criteria: i) They were confirmed to be affected by advanced NSCLC; ii) they received anti-PD-1 monoclonal antibody therapy as the second-line treatment; iii) they were available for a 3-year follow-up period; and iv) they signed informed consent forms. The exclusion criteria were as follows: i) Other types of tumors; ii) history of myocardial infarction, unstable angina, cerebral apoplexy or uncontrollable arrhythmias; iii) pregnancy or lactation period; iv) history of mental disorders; and v) poor compliance with the study protocol. According to these criteria, 24 patients with NSCLC were selected and enrolled in the study, including 13 patients with primary and 11 with metastatic carcinoma. Primary tumor or metastatic lymph node samples were collected from each patient by needle biopsy prior to the start of monoclonal antibody therapy and stored following formalin fixation and embedding in paraffin at room temperature. All patients were evaluated by examination of the samples and computed tomography or magnetic resonance imaging using the Response Evaluation Criteria in Solid Tumors version 1.1 (RECIST 1.1)(27). Responses to treatment were assessed every 6 weeks by computed tomography or magnetic resonance imaging using RECIST 1.1 and confirmed by a subsequent evaluation ≥4 weeks from the start of treatment. Highly selective humanized monoclonal IgG4 antibodies against PD-1/PD-L1, including Nivolumab, IBI308 and Duravalumab, were administered by intravenous infusion every 2 weeks according to the treatment regimen prescribed by each patient's primary care physician. Based on the response to the anti-PD-1 treatment (28), the primary carcinoma group was subdivided into the responding (n=7) and the non-responding (n=6) groups. Similarly, the metastatic carcinoma group was classified into the responding (n=5) and the non-responding (n=6) groups. Overall survival was defined as the time between the start of treatment until death. Progression-free survival defined as the time between the start of treatment, that a patient lives with NSCLC but it does not get worse.

RNA extraction

RNA extraction was performed using the MagMAX-96 Total RNA Isolation kit (cat. no. AM1830; Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. Samples were sonicated and incubated with proteinase K for 15 min at 56°C for protein digestion, followed by 1 h at 80°C for the disruption of nucleic acid-protein cross-links. The digested samples were centrifuged at 12,000 × g for 5 min at 4°C, and the collected supernatant was treated with DNase I. Subsequently, B1 buffer and 100% ethanol were added into the supernatant to produce the binding solution, which was transferred to a 0.45-µm cellulose acetate microcentrifuge spin column and concentrated. Following washing with a washing buffer, the column was eluted with an elution buffer to collect the RNA solution. The RNA integrity number (RIN) in the RNA solution was measured using the 2100 Bioanalyzer Instrument (Agilent Technologies, Inc.). All samples had a RIN value >7. Finally, the RNA solution was subjected to spectrophotometric analysis for determining the A260/A280 and A260/A230 ratios. The A260/A280 ratio ranged between 1.9 and 2.1, and the A260/A230 ratio was >2.0.

Reverse transcription and library construction

The RNA was reverse transcribed into cDNA using the High-Capacity cDNA Reverse Transcription kit (Applied Biosystems; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. The targeted cDNA was amplified with the multiplex immune response primer pool (included in Oncomine IO Panel; Genecast Biotechnology Co., Ltd.) targeting 395 genes. Following amplification, the amplicons were partially digested using a FuPa Reagent (10 min at 50°C, 10 min at 55°C, 20 min at 60°C and 60 min at 10°C). Barcode adapters were subsequently ligated to the partially digested amplicons (30 min at 22°C, 10 min at 72°C and 60 min at 10°C). The barcode-tagged amplicons were purified and amplified. The amplified products were then dissolved in the low EDTA TE buffer for preparing the RNA library. The library concentration was determined using NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Inc.), and the quality of libraries was analyzed by measuring the gene length in the library using the 2100 Bioanalyzer Instrument.

Sequencing

Qualified RNA libraries were quantified, and the concentrations of different libraries were presented in molar concentration according to the average fragment length. All the libraries were pooled based on fragment lengths and the used chip, and sequenced on the Ion Torrent S5 platform. The sequencing kits used included Ion 540 Chef Reagents (cat. no. A27758), Ion S5 chef solutions (cat. no. A27754), Ion S5 chef supplies (cat. no. A27755), I Ion S5 sequencing reagents (cat. no. A27768) and Ion S5 Sequencing solutions (cat. no. A27767). The direction of sequencing was single-end, and the loading concentration was determined by Qubit. The nucleotide length and loading concentration of the final library are presented in Table I.

Table I.

Type of sequencing and loading concentration of the final library.

Sample no.	Nucleotide length, bp	Loading concentration of the final library, pM
201600596	204	45,177.3
201526397	202	66,084.3
201509607	202	18,998.1
201517774	203	54,449.1
201518131	201	56,903.4
201600702	198	40,450.5
201723630	201	59,812.2
201614107	202	25,179.3
201811713	200	46,904.4
201726177	199	60,266.7
201726775	201	5,9085
201619238	201	32,269.5
201613253	199	57,630.6
201729934	200	53,358.3
201729044	200	66,811.5
20180193	204	44,541.0
201817287	201	51,176.7
201814975	195	57,994.2
201811594	200	50,813.1
201800133	200	48,904.2
201801858	200	56,085.3
201703656	204	67,902.3
201715487	200	42,177.6
201816274	200	41,450.4

RNA IO panel sequencing and data normalization

The RNA IO panel sequencing produced 1–2 million reads per sample. Data normalization and processing were performed as previously described (29). Briefly, 10 HK genes were used as endogenous controls. The absolute readout of each HK gene was compared against a predetermined HK reads per million (RPM) profile. The baseline HK RPM profile was established by measuring the average RPM of 10 replicates of GM12878 cell line samples across various sequencing runs. The fold-change ratio for each HK gene was calculated as follows: Ratio of HK = absolute read count of HK/RPM profile of HK. The median value of all HK ratios was then used as the normalization ratio for each sample: Normalization ratio = Median of all HK ratios. The normalized RPM (nRPM) of all genes of each sample was calculated as: nRPM (sample S, gene G) = Absolute read count (sample S, gene G)/Normalization ratio (sample S).

Data analysis

The sequencing data were subjected to quality control according to the following standards: Mapped reads ≥200,000; number of detected HK genes ≥6; and valid reads (on-target ratio) ≥67%. The sequencing quality control data of all patients are presented in Table II. The gene expression levels in the qualified samples were quantified and normalized as follows: Firstly, the RPM value of HK genes in each sample was divided by the standard HK gene value to generate a raw value, and the values of all samples were calculated accordingly; secondly, the median of all the values of samples was presented as the normalized value; finally, the RPM of each sample was divided by the normalized value to produce the corresponding nRPM. An R/Bioconductor software package ‘limma’ (30) running on the R 3.5.3 software (31) was used to perform differential expression analysis according to the log2(nRPM+1) value. The DEGs were determined using the ‘limma’ package with a false discovery rate of 0.05 and absolute fold-change ≥2. The enrichment analysis of DEGs was performed using the ‘ClusterProfiler’ package (32), including Gene Ontology (GO; molecular function, biological process and cellular component; http://geneontology.org), Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG; http://www.genome.jp/kegg) database and Reactome Pathway Database (https://reactome.org), to determine the primary functions of the DEGs as well as the associated metabolic and signaling pathways. Subsequently, differential analysis of gene sets was performed using the gene set variation analysis (GSVA) package for R (33). A receiver operating characteristic (ROC) curve of the DEG set was used to calculate the cut-off values.

Table II.

The sequencing quality control result of all the patients.

Sample no.	Mapped reads, n	On target reads, n	Valid reads	Detected genes, n	HK genes, n	QC result
201600596	2,250,344	1,928,995	86%	362	10	Pass
201526397	2,072,309	1,830,263	88%	361	10	Pass
201509607	2,215,655	2,040,175	92%	353	10	Pass
201517774	2,538,679	2,367,826	93%	380	10	Pass
201518131	2,710,910	2,464,217	91%	376	10	Pass
201600702	2,793,598	2,578,212	92%	376	10	Pass
201723630	2,068,023	1,888,519	91%	366	10	Pass
201614107	2,270,331	2,074,856	91%	358	10	Pass
201811713	2,498,210	2,260,380	90%	369	10	Pass
201726177	2,367,111	2,152,414	91%	360	10	Pass
201726775	1,085,151	950,267	88%	357	10	Pass
201619238	2,049,446	1,837,943	90%	350	10	Pass
201613253	2,444,755	2,256,998	92%	361	10	Pass
201729934	2,597,198	2,342,932	90%	360	10	Pass
201729044	2,107,692	1,897,977	90%	365	10	Pass
20180193	2,414,549	2,192,169	91%	369	10	Pass
201817287	2,358,980	2,077,554	88%	364	10	Pass
201814975	2,240,242	2,056,542	92%	373	10	Pass
201811594	2,371,493	2,164,936	91%	356	10	Pass
201800133	2,427,701	2,202,896	91%	365	10	Pass
201801858	2,461,863	2,186,381	89%	364	10	Pass
201703656	2,116,778	1,929,443	91%	368	10	Pass
201715487	2,359,914	2,189,056	93%	369	10	Pass
201816274	2,339,855	2,138,627	91%	361	10	Pass

HK, housekeeping; QC, quality control.

Statistical analysis

Data are presented as the mean ± standard deviation. Data were analyzed by SPSS 20.0 (IBM Corp.) and R 3.5.3 software. The expression level of the DEG set for each patient was calculated by the mean of log2(nRPM+1) values. The mean value was calculated by summing the four upregulated gene expression values, subtracting the two downregulated gene expression values and by dividing the result by 6. The survival analysis for various groups was performed using the Kaplan-Meier survival analysis with the log-rank test. Fisher's exact test was used for the analysis of patient characteristics. Mann-Whitney U test was used to assess the DEG expression levels in the responding and non-responding groups. P<0.05 was considered to indicate a statistically significant difference.

Results

Associations between clinicopathological features and the efficacy of anti-PD-1 monoclonal antibody treatment

The present study performed the RNA IO panel sequencing on samples form 24 patients with NSCLC prior to anti-PD-1 therapy. Following therapy, all 24 patients were divided into the responding (partial response or stable disease; n=7 primary and 5 metastatic cancer cases) and non-responding (progressive disease; n=6 primary and 6 metastatic cancer cases) groups based on a previously reported response pattern (28). As illustrated in Table III, there were no significant differences in the age, sex, pathological diagnosis, clinical stage, therapeutic regimens and history of smoking between the responding and non-responding groups (P>0.05). These results suggested that the assessed clinicopathological features were not associated with the efficacy of anti-PD-1 therapy in these patients.

Table III.

Clinicopathological features of all patients with NSCLC.

	Primary NSCLC (n=13)			metastatic NSCLC (n=11)
Characteristic	Responding (n=7)	Non-responding (n=6)	P-value	Responding (n=5)	Non-responding (n=6)	P-value
Age, years	54.25	52.79		58.73	53.83
Sex, n
Male	7	5	0.462	5	5	>0.999
Female	0	1		0	1
Disease stage, n
III B	1	0	>0.999	0	0	NA
IV	6	6		5	6
Pathological diagnosis, n
Adenocarcinoma	4	2	0.592	3	4	0.592
Squamous cell carcinoma	3	4		2	2
Therapeutic regimen, n
Duravalumab	4	1	0.266	3	1	0.437
IBI308	0	1		1	1
Nivolumab	3	4		1	4
History of smoking, n
No	6	2	0.103	2	5	0.242
Yes	1	4		3	1

NSCLC, non-small cell lung cancer.

Identification of DEGs between the responding and non-responding groups of patients with NSCLC

The gene expression levels in tumor tissues were analyzed and compared between the responding and non-responding groups of patients with NSCLC patients. As demonstrated in Fig. 1A and B, six genes were identified as significant DEGs between the responding and non-responding groups of patients with primary NSCLC; among them, four were upregulated and two were downregulated in the responding group compared with the non-responding group. Furthermore, the heatmap of DEGs revealed the downregulation of the levels of HLA-F antisense RNA 1 (HLA-F-AS1) and neutrophil cytosolic factor 1 (NCF1), involved in autoimmune diseases, as well as the upregulation of the levels of transcription factor RAR-related orphan receptor C (RORC), deleted in malignant brain tumors 1 (DMBT1), involved in the interaction between the tumor cells and immune system, killer cell lectin-like receptor F1 (KLRF1) and interleukin-18 (IL-18) in the responding group compared with those in the non-responding group (Fig. 1C).

Figure 1.

Effective gene expression map of tumor tissues from the responding and non-responding groups. (A) Volcano plot of DEGs in patients with primary NSCLC. (B) Box plot of DEG expression levels in patients with primary NSCLC. (C) The expression levels of DEGs associated with various pathways in patients with primary NSCLC. (D) volcano map in patients with metastatic NSCLC. (E) Box plot of DEG expression levels in patients with metastatic NSCLC. (F) The expression levels of DEGs associated with various pathways in patients with metastatic NSCLC. NSCLC, non-small cell lung cancer; DEG, differentially expressed gene; HLA-F-AS1, HLA-F antisense RNA 1; NCF1, neutrophil cytosolic factor 1; RORC, RAR-related orphan receptor C; DMBT1, deleted in malignant brain tumors 1; KLRF1, killer cell lectin-like receptor F1; IL-18, interleukin-18; HLA-A, major histocompatibility complex class IA; HLA-DPA1, major histocompatibility complex class II DP α1; TNFSF18, tumor necrosis factor ligand superfamily member 18; IFI6, interferon α-inducible protein 6; PTK7, inactive tyrosine-protein kinase 7.

The gene expression analysis also identified major histocompatibility complex class IA (HLA-A), major histocompatibility complex class II DP α1 (HLA-DPA1), tumor necrosis factor ligand superfamily member 18 (TNFSF18), interferon α-inducible protein 6 (IFI6) and inactive tyrosine-protein kinase 7 (PTK7) as DEGs between the responding and non-responding patients with metastatic NSCLC (Fig. 1D and E). As demonstrated in the heatmap in Fig. 1F, the expression levels all five DEGs were upregulated in the responding group compared with those in the non-responding group.

Enrichment of DEGs in patients with primary or metastatic NSCLC

Functional enrichment analysis of the six DEGs in patients with primary NSCLC was next performed. As presented in Fig. 2, the six DEGs were significantly enriched in the ‘phosphatidylinositol 3-kinase signaling’, ‘T-helper cell differentiation’, ‘phosphatidylinositol-mediated signaling’ and ‘antigen processing-cross presentation’ signaling pathways. Similarly, the significantly enriched signaling pathways of the five DEGs in patients with metastatic NSCLC mainly included ‘interferon signaling’, ‘cytokine signaling in immune system’, ‘antigen processing and presentation of exogenous peptide antigen’, ‘antigen binding’ and ‘interferon-gamma-mediated pathway’ (Fig. 3).

Figure 2.

Enrichment analysis of six differentially expressed genes from patients with primary non-small cell lung cancer. (A) Gene Ontology analysis, including biological process, cell component and molecular function. (B) Kyoto Encyclopedia of Genes and Genomes pathway analysis. (C) Reactome enrichment analysis.

Figure 3.

Enrichment analysis of five DEGs in patients with metastatic NSCLC. (A) Gene Ontology analysis, including biological process, cellular component and molecular function. (B) Kyoto Encyclopedia of Genes and Genomes analysis. (C) Reactome enrichment analysis.

When all 395 analyzed genes in patients with primary or metastatic NSCLC were subjected to cluster analysis of 35 immune-related signaling/functional pathways using GSVA, the results demonstrated that there were no significant differences in any of these pathways between the responding and non-responding groups (adjusted P>0.05; Fig. 4).

Figure 4.

Associations between DEGs and PFS. (A) Associations between DEGs and PFS in patients with primary NSCLC. (B) Associations between DEGs and PFS in patients with metastatic NSCLC. DEGs, differentially expressed genes; PFS, progression-free survival; NSCLC, non-small cell lung cancer.

DEG set-based prediction of the efficacy of anti-PD-1 therapy in patients with primary NSCLC

Cluster analysis of the six DEGs in 13 patients with primary NSCLC was further performed (Fig. 5A). Although the results of one patient deviated from the others, the remaining patients had the correct attribution and could be distinguished, suggesting that the gene set composed by these six DEGs could fairly distinguish the patients between the responding and non-responding groups. Therefore, the present study sought to determine whether the six DEGs may form a DEG set for predicting the effect of anti-PD-1 treatment in these patients. As illustrated in Fig. 5B, the expression levels of the DEG set in the responding group were significantly higher compared with those in the non-responding group (P=0.011). Subsequently, the ROC curve was used to determine the cut-off value indicating the efficacy of anti-PD-1 treatment. The expression level of the DEG set exhibited statistical significance in predicting the efficacy of anti-PD-1 treatment [area under the curve (AUC) =0.881; P=0.022; sensitivity, 100%; specificity, 83.3%; Fig. 5C). Based on the calculated cut-off value, patients with an expression level of the DEG set >6.501 benefited from anti-PD-1 therapy.

Figure 5.

DEG set-based prediction of the efficacy of anti-PD-1 treatment in patients with primary NSCLC. (A) Cluster analysis of the six DEGs identified in patients with primary NSCLC. (B) The expression levels of the DEG set in the responding group were significantly higher compared with those in the non-responding group. (C) Receiver operating characteristic curve was used to predict the efficacy of anti-PD-1 treatment. Patients with a DEG set expression value >6.501 benefited from anti-PD-1 therapy. PD-1, programmed death-1; NSCLC, non-small cell lung cancer; DEG, differentially expressed gene; AUC, area under the curve.

DEG set-based prediction of the efficacy of anti-PD-1 treatment in patients with metastatic NSCLC

Cluster analysis of the five DEGs in patients with metastatic NSCLC was subsequently performed (Fig. 6A). The results demonstrated that the set of five DEGs did not fully distinguish between patients in the responding and non-responding groups. Therefore, cluster analysis was performed on subsets of the 5 DEGs. As presented in Fig. 6B, the DEG set comprising HLA-A and TNFSF18 effectively distinguished the two groups of patients. The expression level of the gene set was calculated by the mean of log2(nRPM+1) values. Notably, the responding group displayed a significantly higher expression levels of the DEG set compared with those in the non-responding group (P=0.004; Fig. 6C). As indicated by the cut-off value (Fig. 6D), the expression level of the gene set had statistical significance in predicting the efficacy of anti-PD-1 monoclonal antibody treatment (AUC=0.967; P=0.011; sensitivity, 100%; specificity, 83.3%); patients with an expression level of the DEG set >6.741 benefited from anti-PD-1 antibody therapy.

Figure 6.

DEG set-based prediction of the efficacy of anti-PD-1 treatment in patients with metastatic NSCLC. (A) Cluster analysis of the five DEGs identified in patients with metastatic NSCLC. (B) Further cluster analysis of the gene set for distinguishing patients with metastatic NSCLC. (C) The expression levels of the DEG set in the responding group were significantly higher compared with those in the non-responding group. (D) Receiver operating characteristic curve was used to predict the efficacy of anti-PD-1 treatment. Patients with a DEG set expression value >6.741 benefited from anti-PD-1 therapy. PD-1, programmed death-1; NSCLC, non-small cell lung cancer; DEG, differentially expressed gene; AUC, area under the curve.

Associations between DEGs and progression-free survival (PFS)

Lastly, the present study conducted survival analysis on 13 patients with primary carcinoma and 10 patients with metastatic carcinoma. As demonstrated in Fig. 7A, in the primary carcinoma group, patients with a longer PFS exhibited higher expression levels of DMBT1, KLRF1, RORC and the 6-gene set compared with those in patients with a shorter PFS. Similarly, in the metastatic carcinoma group, patients with a longer PFS displayed higher expression levels of HLA-A, TNFSF18 and the 2-gene set compared with those in patients with a shorter PFS (Fig. 7B).

Figure 7.

Survival analysis of patients with primary and metastatic carcinoma. (A) Among patients with primary carcinoma, high expression levels of DMBT1, KLRF1, RORC and the 6-gene set were observed in patients with longer PFS. (B) In metastatic carcinoma, patients with a longer PFS displayed high expression levels of HLA-A, TNFSF18 and the 2-gene set. PFS, progression-free survival; OS, overall survival; HLA-F-AS1, HLA-F antisense RNA 1; NCF1, neutrophil cytosolic factor 1; RORC, RAR-related orphan receptor C; DMBT1, deleted in malignant brain tumors 1; KLRF1, killer cell lectin-like receptor F1; IL-18, interleukin-18; HLA-A, major histocompatibility complex class IA; TNFSF18, tumor necrosis factor ligand superfamily member 18.

Discussion

Since surgery is ineffective for patients with advanced NSCLC, chemotherapy remains the preferred treatment option (34). Anti-PD-1 monoclonal antibody, an immune checkpoint inhibitor, has provided a breakthrough in the treatment of patients with advanced NSCLC (35). The role of anti-PD-1 antibodies in the first- and second-line treatment of NSCLC or in the local adjuvant therapy of NSCLC has been established in previous studies (36–38). However, the effectiveness of anti-PD-1 therapy on tumors is limited, and rapid growth of tumors is observed in a number of patients (35). The high cost and difficulty in predicting the efficacy have become the bottleneck for the promotion of anti-PD-1 therapy (39). Thus, there is an urgent need to identify efficient and precise biomarkers for screening patients with NSCLC that may respond to anti-PD-1 antibody therapy. Analysis of the functional mechanism of anti-PD-1 monoclonal antibody suggests that immune-related genes and pathways serve an important role in its antitumor effect (40). In the present study, RNA IO panel sequencing was used to examine the expression levels of 395 genes associated with immune pathways in patients with primary or metastatic NSCLC prior to the standard anti-PD-1 antibody therapy (41). Literature review and data analysis revealed that five immune-related genes and two gene sets may potentially be used for predicting the therapeutic effects of PD-1 inhibitors in patients with primary or metastatic NSCLC.

In the present study, among patients with primary NSCLC, the responding group exhibited lower expression levels of HLA-F-AS1 and NCF1 compared with those in the non-responding group. HLA-F-AS1 is a long non-coding RNA that is significantly downregulated in human lung adenocarcinoma tissues compared with matched adjacent non-tumor tissues (42). The NCF1 protein is an essential component of the phagocytic NADPH oxidase type 2, which is involved in autoimmune inflammatory disorders (43). Kelkka et al (44) have reported that mice lacking NCF1 developed markedly fewer Lewis lung carcinoma tumors compared with those in the wild-type controls. Consistently, the results of the present study demonstrated that patients with primary NSCLC with a longer PFS exhibited higher expression levels of HLAF-AS1 and NCF1 compared with those in patients with a shorter PFS. Thus, low levels of HLA-F-AS1 and NCF1 may be biomarkers for predicting response of patients with primary NSCLC to anti-PD-1 therapy. In addition, low expression levels of HLA-F-AS1 may indicate improved efficacy of anti-PD-1 treatment (45,46). DMBT1 has been proposed as a candidate tumor suppressor (45,46). DMBT1 is highly expressed in normal lung tissues, but is present at low levels in lung cancer cell lines and primary NSCLC tissues (45). In the present study, among patients with primary NSCLC, the responding group exhibited higher levels of DMBT1 compared with those in the non-responding group, whereas increased expression levels of DMBT1 were present in patients with a longer PFS compared with those in patients with a shorter PFS. Although DMBT1 is lowly expressed in patients with NSCLC, its relatively high expression levels may potentially be used as an index for predicting the efficacy of anti-PD-1 treatment in patients with primary NSCLC.

Among patients with metastatic NSCLC in the present study, the responding group presented with significantly higher levels of HLA-A and TNFSF18 compared with those in the non-responding group. HLA-A belongs to the HLA class I antigens and serves a crucial role in presenting tumor cell immunogenic polypeptide to T cells as well as promoting the antitumor effects of cytotoxic T lymphocytes (47,48). However, HLA-A levels are markedly downregulated in the majority of primary NSCLC tumors and all metastatic lymph nodes compared with those in normal lung tissues (49). TNFSF18, also termed glucocorticoid-induced TNFR-related protein (GITRL), participates in the functioning of effector and regulatory T cells, which is important for the development of immune responses (50). Upregulation of GITRL has been demonstrated to improve antitumor immunity in murine Lewis lung carcinoma (51,52). In addition, in the present study, patients with metastatic NSCLC with a longer PFS presented with higher expression levels of HLA-A and TNFSF18 compared with those in patients with a shorter PFS. Therefore, patients with metastatic NSCLC with high expression levels of HLA-A and TNFSF18 may benefit from anti-PD-1 treatment, suggesting that HLA-A and TNFSF18 may be potential biomarkers for predicting the efficacy of anti-PD-1 therapy in patients with metastatic NSCLC. PTK7 is a member of the receptor protein tyrosine kinase family (53). Studies have demonstrated that PTK7 is highly expressed in tumor tissues of patients with primary lung adenocarcinoma, and inhibition of PTK7 reduces the number of tumor-initiating cells and induces tumor regression (53,54). By contrast, one study has reported that the mRNA and protein expression levels of PTK7 are downregulated in human lung squamous cell carcinoma compared with those in normal lung tissues, and overexpression of PTK7 in lung cancer cells inhibits cell proliferation, invasion and migration (55). Thus, it remains to be determined whether PTK7 is associated with the development of NSCLC or the response to anti-PD-1 treatment.

Single-gene predictive biomarkers are usually considered unsatisfactory in terms of accuracy and precision. In recent years, an increasing number of studies have demonstrated that biomarkers consisting of gene sets (multiple DEGs) are more accurate compared with single-gene biomarkers (56,57). Li et al (58) have established a 4-gene set biomarker that predicts early relapse in advanced epithelial ovarian cancer after initial platinum-paclitaxel chemotherapy with an accuracy ~65.5%. In addition, a minimal driver gene set has been developed to predict bone metastasis in breast cancer (59). Another study has proposed that an immune gene-set based signature may serve as a promising biomarker for estimating overall survival of patients with ovarian cancer (60). The results of the present study demonstrated a gene set comprising six DEGs (HLA-F-AS1, NCF1, RORC, DMBT1, KLRF and IL-18) may be used for predicting the efficacy of anti-PD-1 therapy in patients with primary NSCLC; specifically, patients with a calculated expression level of the DEG set >6.501 may benefit from anti-PD-1 therapy. In addition, a DEG set comprising two DEGs (HLA-A and TNFSF18) may be applied to predict the efficacy of anti-PD-1 therapy in patients with metastatic NSCLC. Patients with an expression level of the gene set >6.741 may benefit from anti-PD-1 monoclonal antibody treatment.

The present study had certain limitations due to the small sample size. In addition, there were no overlapping DEGs or gene sets observed for both primary and metastatic NSCLC in the present study. In two previous studies (61,62), patients with primary and metastatic cancer also exhibited inconsistent gene expression alterations; this problem should be addressed in depth in future studies.

In summary, the present study conducted RNA IO panel sequencing to identify potential biomarkers for predicting the response to anti-PD-1 therapy in patients with primary or metastatic NSCLC. The results of the present study demonstrated that five immune-related DEGs and two DEG sets may be used, respectively, as single and combination biomarkers for the prediction of treatment efficacy. Although these results provided a basis for identification of additional biomarkers to predict the response to anti-PD-1 treatment, they need to be verified in further studies.