TCR Repertoire Diversity of Peripheral PD-1⁺CD8⁺ T Cells Predicts Clinical Outcomes after Immunotherapy in Patients with Non–Small Cell Lung Cancer

Immune checkpoint blockades (ICB), such as programmed cell death 1 (PD-1) and programmed death ligand 1 (PD-L1) antibodies, have shown clinical efficacy in many cancer types (1–4). However, although biomarkers can personalize the delivery of ICB therapy, only a subset of patients benefit from ICB therapy. Despite advances in identifying biomarkers, such as PD-L1 expression (1, 5), tumor mutational burden (TMB; refs. 6, 7), immune-regulated gene expression scores, and tumor-infiltrating lymphocytes (8), few patients can provide tumor tissue samples of enough quality and quantity to support the necessary biomarker tests. Here, we explore noninvasive blood-based approaches to identify patients likely to respond to ICB.

The TMB measured via circulating tumor DNA (ctDNA) in the blood is associated with response to ICBs and subsequent progression-free survival (PFS; ref. 9). A high TMB indicates elevated production and release of tumor neoantigens, which improves tumor cell recognition and eradication by cytotoxic T cells invigorated after ICB (10). Activation of the host immune response against cancer cells includes recognition of neoantigen peptides by clonally proliferating T-cell receptors (TCR; ref. 11). Thus, TCR-based biomarkers might predict the response to ICBs. Sequencing of the TCR repertoire can be conveniently performed with samples of peripheral blood, presenting a noninvasive approach to predict response to ICB. Studies of the peripheral TCR repertoire in various cancer types have indicated that the TCR repertoire could serve as either a predictive or a prognostic biomarker for cytotoxic T-lymphocyte–associated protein 4 (12, 13) or PD-1 inhibitors (14, 15). However, these conclusions remain controversial (16, 17) in other studies that used isolated peripheral blood mononuclear cells (PBMC) for TCR sequencing to predict the response to ICBs. The utility of cancer-specific TCR-based biomarkers for ICB prediction needs further analysis.

The double-positive PD-1 and CD8 (PD-1⁺ CD8⁺) T cells are a subpopulation of T lymphocytes that are targeted by PD-1 blockade; this subpopulation contains exhausted CD8⁺ T cells whose inhibition is mediated by the PD-1 pathway (18). Blocking the PD-1 pathway can partially reinvigorate T cells and lead to the proliferation of PD-1⁺ CD8⁺ T cells. The early proliferation of peripheral PD-1⁺ CD8⁺ T cells is associated with positive clinical outcomes after anti–PD-1 therapy (19, 20). A study on melanoma revealed that preexisting PD-1⁺ CD8⁺ T cells were functional cytotoxic T cells that targeted cancer and experienced exhaustion, which was induced by an immune checkpoint pathway: these cells could be detected in the peripheral blood (21). These results suggested that the TCR repertoire profiling of isolated peripheral PD-1⁺ CD8⁺ T cells might serve as a biomarker for predicting the response and prognosis associated with anti–PD-1 therapy.

In this context, we sequenced the complementarity-determining region 3 (CDR3) in isolated peripheral PD-1⁺ CD8⁺ TCR variable β (V-β) to investigate its value for predicting the response to anti–PD-1/PD-L1 therapy in patients with non–small cell lung cancer (NSCLC). We aimed to develop a noninvasive approach to facilitate the risk–benefit stratification of patients receiving ICB treatments.

This study involved two independent patient cohorts (A and B). To examine the ability of TCR diversity to predict the response to anti–PD-1/PD-L1 therapy, 25 patients were enrolled in the training set (cohort A, from March 14th, 2017, to November 21st, 2017). To confirm these findings, an additional 15 patients were enrolled in the validation set (cohort B, from November 23th, 2017, to May 2nd, 2018; Table 1). A total of 40 of the 51 patients with stage IIIB–IV NSCLC who were receiving on-study anti–PD-1/PD-L1 therapy were recruited (Fig. 1). Pre- (baseline, before the first administration of ICBs) and post-ICB (the timepoint of the first imaging evaluation, approximately 4–6 weeks after the first administration of ICBs) peripheral blood samples were prospectively collected, and PD-1⁺ CD8⁺ T cells were isolated through flow cytometry for TCR sequencing. CtDNA was extracted for mutation analysis, and the diversity and clonality of the TCR repertoire were examined for biomarker exploration.

This study was approved by the ethics committees of the National Cancer Center/Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College (NCC2016JZ-03 and NCC2018–092). All enrolled patients provided written informed consent.

PBMCs were isolated from fresh peripheral blood (20 mL) containing an anticoagulant by density gradient centrifugation using Lymphoprep (Progen). Based on FACS analysis (BD FACSAriaTM), PD-1⁺ CD8⁺ T cells were purified from the PBMCs. The following human-protein–specific flow cytometry antibodies were purchased from eBioscience: CD3-efluor 450 (clone: OKT3), CD8-APC (clone: RPA-T8), and CD279 (PD-1)-PE (MIH4).

DNA from PD-1⁺ CD8⁺ T cells was isolated using a commercially available kit (QIAamp DNA Mini Kit, QIAGEN; catalog: 51306). The total DNA yield was greater than or equal to 1 μg. The UV absorption ratios of wavelengths 260/280 and 260/230 were greater than or equal to 1.8 and 2, respectively.

Flow cytometry was used for a PD-1 antibody competition assay. PBMCs were isolated, collected, and then incubated (1 × 10⁷ cells) with increasing concentrations of therapeutic PD-1 antibodies (0, 1, and 100 μg/mL) for 30 minutes. The PD-1⁺ ratios of CD279 (PD-1)-PE (MIH4) in CD8⁺ T cells in PBMCs that were treated or untreated with a therapeutic PD-1 antibody were compared.

CDR3 in the TCRβ chain (TRB) was inclusively and semiquantitatively amplified by multiplex PCR; the multiplex amplification included both first and second rounds of PCR (PCR1 and PCR2). The primer sequences have been filed as a part of a Chinese patent (CN105087789A; Supplementary Table S1). During the first round of PCR (PCR1), 10 cycles were used to amplify CDR3 sequences using specific primers against each V and J gene. In the second round of PCR, PCR was performed using universal primers. For PCR1, template DNA (600 ng) was amplified after adding 2× QIAGEN Multiplex PCR Master Mix (25 μL), 5× Q solution (5 μL), the forward primer set pool (1 μL), and the reverse primer set pool (1 μL) to form a reaction system by using a Multiplex PCR Kit (QIAGEN). Then, PCR was performed with 1 cycle of 95°C for 15 minutes, 10 cycles of denaturation at 94°C for 30 seconds, and 10 cycles of annealing at 60°C for 90 seconds and extension for 30 seconds at 72°C. After a final extension for 5 minutes at 72°C, the system was cooled to 4°C. The target fragment of the multiplex PCR products was purified with magnetic beads (Agencourt no. A63882, Beckman). All PCR1 products were used as templates for the second step of amplification after adding pooled primers (2 μL), Phusion master mix prepared using a Phusion High-Fidelity PCR Kit (25 μL; New England Biolabs, America), and nuclease-free water to reach a total volume of 50 μL. The reactions were then transferred to a thermal cycler that carried out the following program: one cycle at 98°C for 1 minute; 25 cycles of denaturation at 98°C for 20 seconds, annealing at 65°C for 30 seconds, and extension at 72°C for 30 seconds; and a final extension at 72°C for 5 minutes. The samples were then held at 4°C. Size selection was performed by agarose gel electrophoresis (400 mA/100 V, 2 hours), and the targeted fragments (between 200 and 350 bp) were retrieved and purified by a QIAquick Gel Purification Kit (QIAGEN). The paired-end sequencing of these samples was carried out with a read length of 151 bp using an Illumina HiSeq 3000 platform (22).

Raw sequencing data were processed and analyzed as follows. (i) Undesired sequences that did not contain the primers were filtered using Cutadapt https://cutadapt.readthedocs.org/, (ii) reads were merged to obtain contigs using Pear (https://cme.h-its.org/exelixis/web/software/pear/doc.html), (iii) sequences were aligned to reference TRB V/(D)/J gene sequences (http://www.imgt.org) using MiXCR to determine the TRB V/(D)/J gene segment in each contig, (iv) the CDR3 region was identified based on the conserved sequence of the CDR3 region, and (v) CDR3 species were clustered to eliminate sequencing errors according to the base quality and sequence similarity.

The immune repertoire was characterized by examining the diversity and clonality. The diversity of the TCR repertoire was calculated based on the Shannon–Wiener index (Shannon index), which is a function of both the relative number of clonotypes present and the relative abundance or distribution of each clonotype (23). The Shannon index is calculated as follows. In the Shannon index, ni is the clonal size of the ith clonotype (i.e., the number of copies of a specific clonotype), S is the number of different clonotypes, and N is the total number of TCR/B-cell receptor sequences analyzed.

Clonality is defined as 1 - (Shannon index)/ln(# of productive unique sequences) (24). A maximally diverse population is associated with a clonal score of 0, and a perfectly monoclonal population is associated with a clonality score of 1. A change in clonality was defined as [clonality(posttreatment) - clonality(baseline)]/clonality(baseline).

Peripheral blood samples were centrifuged within 4 hours to separate the plasma from the peripheral blood cells. Plasma was separated by centrifugation, which included two steps. The first centrifugation was at 1,600 × g for 10 minutes, followed by 16,000 × g for 10 minutes. Circulating cell-free DNA (cfDNA) was isolated from plasma samples using a QIAamp Circulating Nucleic Acid Kit (QIAGEN). Genomic DNA and tissue DNA were extracted using a QIAamp DNA DNeasy Blood & Tissue Kit (QIAGEN). All DNA extraction procedures were performed according to the manufacturer's instructions.

Genomic alterations (mutations, insertions, deletions, and amplifications) were detected in ctDNA extracted from plasma samples using a broad-targeted next-generation sequencing–based 1,021-gene panel, which included most prevalent tumor-related genes in human cancers (Supplementary Table S2). After the isolation of ctDNA by hybrid capture, the assay was performed using molecular barcoding and proprietary bioinformatics algorithms with massively parallel sequencing on an Illumina HiSeq 3000 platform in a Clinical Laboratory Improvement Amendments/College of American Pathologists–accredited laboratory (Geneplus). Variants in plasma ctDNA were assessed in samples from 19 patients that were collected before therapy and after two cycles of treatment with anti–PD-1/PD-L1 immunotherapy. The maximum variant allele frequency (MVAF) was used to examine the changes in ctDNA amounts in each sample (25).

Radiographic imaging was acquired by the indicated approaches, such as CT or MRI for tumor response assessment, which was evaluated by both the investigator and an independent radiologist. Baseline tumor assessments were performed within 1 to 28 days prior to the initiation of anti–PD-1/PD-L1 treatment. Subsequent assessments were performed every 4 to 8 weeks until there was disease progression. The efficacy evaluation was performed by using RECIST version 1.1, which included complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD). Disease control (DC) indicated that patients had achieved DC (i.e., CR, PR, or SD according to RECIST version 1.1). PFS was defined as the time from the start of anti–PD-1/PD-L1 treatment until either objective disease progression (assessed by an investigator using RECIST version 1.1) or death from any cause. Overall survival (OS) was defined as the time from the start of anti–PD-1/PD-L1 treatment until death from any cause. We censored patients who had not progressed, who had died at the time of statistical analysis, or who were lost to follow-up but had not died at the time of their last evaluation.

The optimal thresholds for T-cell diversity were calculated by ROC analysis. The Youden index (J = sensitivity+specificity-1), a measure of overall diagnostic effectiveness, was used in conjunction with ROC analysis. Independent predictive variables for PFS were identified by multivariate Cox regression analysis using a backward selection procedure. The association between the response and high versus low PD-1⁺ CD8⁺ T-cell diversity was assessed by the Fisher exact test. The Kaplan–Meier curve analysis and the log-rank test were used to compare the PFS and OS. The frequency of the top TCRβ variable (TRBV) clone in different subgroups was compared by the Kruskal–Wallis test and the Steel–Dwass test. All reported P values are two-tailed. For all analyses, P < 0.05 was considered statistically significant unless otherwise specified. IBM SPSS software (version 23.0), GraphPad Prism (version 6.01), and R 3.4.2 were used for the statistical analysis.

A total of 40 of 51 patients with stage IIIB–IV NSCLCs who were receiving anti–PD-1/PD-L1 therapy provided pre-ICB peripheral PD-1⁺ CD8⁺ T-cell samples for subsequent TCR sequencing and were ultimately recruited for analysis. All patients received anti–PD-1/PD-L1 therapy as a second-line or higher treatment. The demographic and clinical characteristics of these patients are summarized in Table 1. In the discovery set (cohort A, n = 25), 10 patients obtained a PR, and 6 had SD. In the validation set (cohort B, n = 15), 3 showed a PR, and 4 had SD. When both patient groups were pooled, the patients received a median of 6 cycles (range, 2–31 cycles) of anti–PD-1/PD-L1 treatment. The best response was a PR, which was observed in 13 patients (13/40, 32.5%); SD was observed in 10 patients (10/40, 25%), and PD was observed in 17 patients (17/40, 42.5%). The median follow-up time was 13.1 months (range, 3.2–22.9 months). There were no significant differences between cohorts A and B in terms of the baseline characteristics, including age, gender, the Eastern Cooperative Oncology Group (ECOG) performance score, smoking status, pathologic type, disease stage, and treatment line (Table 1).

In the discovery set (cohort A, n = 25), a positive correlation was observed between the pre-ICB peripheral PD-1⁺ CD8⁺ TCR diversity and a superior response (DC group vs. PD group) to anti–PD-1/PD-L1 therapy [median, 3.72 (3.31–4.24) vs. 2.45 (2.21–2.84), P < 0.001; Fig. 2A]. Based on an optimal threshold of 3.14 for the maximum ROC curve value, TCR diversity was used to dichotomize patients into “high” and “low” diversity subsets (Fig. 2B). Patients with “high” diversity had a significantly longer PFS than those with “low” diversity [6.4 months, 95% confidence interval (CI), 3.8–9.2 vs. 2.5 months, 95% CI, 1.9–3.2; HR, 0.39; 95% CI, 0.17–0.94; log-rank P = 0.021; Fig. 2C]. In the independent validation set (cohort B, n = 15), patients who had been previously classified as having “high” diversity had a higher DC rate (85.7%; 95% CI, 77.4%–96.8% vs. 12.5%; 95% CI, 3.2%–21.4%; P < 0.001) and longer PFS (7.3 months, 95% CI, 3.6–11.1 vs. 2.6 months, 95% CI, 1.6–3.6; HR, 0.26; 95% CI, 0.08–0.86; log-rank P = 0.002; Fig. 2D) than those in the “low” diversity subgroup. Similarly, patients who had been previously classified as having “high” diversity exhibited a prolonged OS, regardless of whether they were in the discovery set (not reached vs. 6.9 months, 95% CI, 3.1–10.7; HR, 0.31; 95% CI, 0.09–0.98; log-rank P = 0.029) or the validation set (not reached vs. 13.8 months, 95% CI, 9.2–18.4; HR, 0.16; 95% CI, 0.03–0.82; log-rank P = 0.058; Fig. 2E and F).

Based on the results described above, the two cohorts (A and B) were merged for further analysis. When stratified by TCR diversity, similar trends for response and PFS were observed in the merged dataset (Supplementary Fig. S1). The pre-ICB PD-1⁺ CD8⁺ TCR diversity exhibited an optimal Youden's index of 0.81, with a sensitivity of 0.87 and a specificity of 0.94, for stratifying the clinical response to ICBs (DC group vs. PD group; Fig. 2G). Patients with “high” TCR diversity also had a significantly prolonged OS after ICB compared with those with “low” diversity (not reached vs. 10.1 months, 95% CI, 6.0–14.2; HR, 0.31; 95% CI, 0.11–0.83; log-rank P = 0.013; Supplementary Fig. S2). In the multivariate Cox regression models, which included TCR diversity, pathologic type, gender, age, ECOG score, and the kind of immunotherapy, only TCR diversity was identified as an independent factor for both PFS (HR, 0.57; 95% CI, 0.37–0.78; P = 0.004; Fig. 3A; Supplementary Table S3) and OS (HR, 0.42; 95% CI, 0.02–0.89; P = 0.033; Fig. 3B). To further validate the predictive performance of PD-1⁺ CD8⁺ TCR diversity for immunotherapy, a control set was designed that analyzed the pre-ICB TCR diversity of PBMCs from 18 patients (Supplementary Table S4). This analysis which failed to show any association with clinical outcomes (response and PFS) associated with ICB treatment (Fig. 4A and B). Thus, isolated PD-1⁺ CD8⁺ T cells provide the basis for predictive biomarkers with clinical utility.

Through a competition inhibition experiment, we showed that the postimmunotherapy PD-1⁺ CD8⁺ T-cell sorting that required MIH4 binding to PD-1 was not disturbed by the previous binding of PD-1 inhibitors (Supplementary Fig. S3). We next examined the early dynamic changes in the TCR clonality of isolated PD-1⁺ CD8⁺ T cells after ICB treatment. In 19 patients with matched pre- and post-ICB peripheral PD-1⁺ CD8⁺ TCR sequencing data, TCR clonality displayed a trend of consistent change that was positively associated with tumor shrinkage (Fig. 5A) and was negatively associated with the MVAF of mutations in the major clones (Fig. 5B). Patients with primary PD or an initial PR underwent decremental or incremental evolution in TCR clonality, respectively (Fig. 5C). Patients with increased TCR clonality demonstrated a superior PFS (7.3 months, 95% CI, 3.3–11.4 vs. 2.6 months, 95% CI, 2.3–2.9; HR, 0.28; 95% CI, 0.11–0.74; log-rank P = 0.002; Fig. 5D) and OS (not reached vs. 7.5 months, 95% CI, 1.8–13.2; HR, 0.23; 95% CI, 0.07–0.79; log-rank P = 0.034; Fig. 5E) compared with those with decreased TCR clonality.

Two patients who had pseudo-PD (no. 16, PR; and no. 3, SD with 25% tumor shrinkage; Supplementary Fig. S4) showed a similar fold change in PD-1⁺ CD8⁺ TCR clonality before versus after anti–PD-1/PD-L1 treatment compared with patients with a PR (mean, 1.54 ± 0.24 vs. 1.51 ± 0.21). Patients with both pseudo-PD and those with a PR exhibited a significantly higher fold change in the TCR clonality than patients with PD (mean, 0.78 ± 0.19; P = 0.016 for pseudo-PD compared with PD and 0.003 for PR compared with PD; Supplementary Fig. S5).

To further demonstrate the similarity between pseudo-PD and PR, the distribution and dynamic changes in TRBV, representing patterns of TCR V-β families, for each patient were analyzed during ICB treatment (Fig. 6). Generally, the frequencies of dominant TRBV clones (rank 1) after ICB were significantly higher than those before ICB treatment in patients with a PR (P = 0.029) and pseudo-PD (P = 0.033) but not in those with SD and PD (P = 0.174; Fig. 6A and B). The dominant TRBV clones in patients with a PR and pseudo-PD obviously expanded from pre-ICB to post-ICB (Fig. 6C and D). For example, in patient no. 16 (with pseudo-PD), the post–ICB-dominant TRBV clone (TRBV28) represented 48% of all PD-1⁺ CD8⁺ TCRs, which was more abundant than the pre–ICB-dominant TRBV clone (TRBV20-1, 19%). The TRBV7-8 clone of this patient expanded from 4% at pre-ICB to 47% at post-ICB.

Finally, we calculated the overlap rate of the top ten clones in paired pre- and post-ICB samples. The mean overlap rates were 57.4% and 64.3% for the 4 patients with a PR and the 2 patients with pseudo-PD, respectively, and were 26.8% and 39.8% in patients with PD and SD, respectively. In addition, the overlap rate in the DC group was significantly higher than that in the PD group (55.7% vs. 26.8%, P = 0.001).

Our analysis revealed that peripheral PD-1⁺ CD8⁺ TCR diversity could predict clinical benefits of anti–PD-1/PD-L1 therapy. The TCR repertoire evolved during immunotherapy, and the early changes in TCR clonality were associated with the clinical outcomes to anti–PD-1/PD-L1 therapy in NSCLC. These results suggested that the peripheral PD-1⁺ CD8⁺ TCR repertoire could be a predictive biomarker for ICB.

The spectrum of TCR epitopes responsible for tumor neoantigen recognition is diverse due to the random formation of neoantigens, which are derived from millions of individualized genetic alterations (26). Thus, even tumors with similar histologic origins can carry diverse “genomic landscapes” (27). The TCR diversity may reflect the probability of neoantigen recognition. The PD-1⁺ CD8⁺ phenotype is a marker of exhausted T cells and has the highest possibility including neoantigen-specific cytotoxic T cells (28). The host immune system needs to maintain a diversified TCR repertoire to recognize the variety of tumor neoantigens (11). Hence, a PD-1⁺ CD8⁺ TCR repertoire with a higher diversity provides more opportunities for tumor neoantigen recognition and indicates that there is a higher proportion of exhausted T cells that can be reinvigorated by ICB, subsequently leading to improved immunologic responses (29).

Results from previous studies on TCR diversity for ICB prediction differed (16, 17), perhaps due to the existence in those studies of large numbers of nontumor antigen–specific TCRs that diluted the tumor neoantigen–specific TCRs. This speculation was supported by our TCR analysis in PBMCs versus specific PD-1⁺ CD8⁺ T cells, which showed that TCR repertoire sequencing in nonisolated PBMCs did not correlate with the clinical outcomes of anti–PD-1/PD-L1 therapy. Thus, isolation of specific T cells may help evaluation of TCR diversity as an immunotherapy biomarker.

In addition, neoantigen-specific T cells were expected to experience clonal expansion after responding to ICB treatments, which can be reflected by TCR clonality (30). In the present study, patients with increased PD-1⁺ CD8⁺ TCR clonality had improved survival outcomes compared with those with decreased TCR clonality. PD-1⁺ CD8⁺ TCR clonality displayed a consistent trend of change that was positively correlated with the tumor burden index (tumor shrinkage and MVAF). These results indicated that PD-1⁺ CD8⁺ T cells underwent clonal expansion after responding to anti–PD-1/PD-L1 therapy. Our conclusion was partially supported by a study by Kim and colleagues (19) that described the proliferative response of peripheral blood PD-1⁺ CD8⁺ T cells after the first week of treatment, which was correlated with an improved response and survival prognosis to anti–PD-1 therapy. These results and those of our study suggest that the dynamic monitoring of functional T cells or of their TCR repertoire may be a useful complementary and potential surrogate biomarker for the early identification of patients benefiting from ICB when a pre-ICB predictor is not available.

Pseudo-PD often poses difficulties for clinicians when determining whether to continue with ICB delivery (31). In our study, patients who were identified as having pseudo-PD showed a fold change in PD-1⁺ CD8⁺ TCR clonality and in the expansion of dominant TCR clones that was similar to patients who experienced a PR but not to patients who had PD. Thus, monitoring the TCR repertoire aids ICB prediction. To further clarify the details of TCR clonal evolution after anti–PD-1/PD-L1 therapy, we analyzed the TRBV overlap rate to compare the similarity of TCR clones before and after anti–PD-1/PD-L1 treatment. We found that better clinical outcomes were accompanied by greater PD-1⁺ CD8⁺ TCR repertoire similarities and with a consistent change in TCR clonality during treatment, suggesting that the improved anticancer immunity resulted from preexisting functional T-cell clones that were maintained and increased during immunotherapy. Although the small sample size limited the reliability of our conclusions, our results suggest a potential approach for the early identification of pseudo-PD, which should be further validated in larger cohorts.

Another limitation of our study is that the TCR sequences identified by current technology could not exactly represent the functional T cells that are specific to an individualized tumor-specific neoantigen, even in the isolated PD-1⁺ CD8⁺ T cells. Assays including epitope mapping or antigen discovery could be more accurately interpreted if paired blood and tumor tissue samples were included in the analysis. Moreover, some confounding variables, such as PD-L1 expression and the TMB, were not included in our analysis due to the limited number of available samples, which might have led to a potential bias. Finally, these findings are exploratory and need to be confirmed in additional, larger cohorts.

In conclusion, the sequencing of the peripheral PD-1⁺ CD8⁺ TCR repertoire suggests a noninvasive approach to stratify patients with NSCLC treated with anti–PD-1/PD-L1 treatments.

No potential conflicts of interest were disclosed.

Conception and design: Z. Wang, J. Wang

Development of methodology: J. Han, Y. Wang, S. Chen, Z. Yao

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J. Han, J. Duan, R. Wan, X. Wang, S. Chen, D. Wang, K. Fei, Z. Yao, S. Wang, Z. Lu, Z. Wang, J. Wang

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Han, J. Duan, S. Chen, D. Wang, S. Wang, Z. Wang, J. Wang

Writing, review, and/or revision of the manuscript: J. Han, J. Duan, H. Bai, D. Wang, S. Wang, Z. Lu, Z. Wang, J. Wang

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J. Han, J. Duan, Y. Wang, Y. Tian, Z. Wang

Study supervision: J. Wang

The authors thank all patients who were involved in this study. The authors also thank the support from Geneplus-Beijing.

This work was supported by the National Natural Sciences Foundation Key Program (81630071); CAMS Innovation Fund for Medical Sciences (CIFMS 2016-I2 M-3-008 and 2017-I2 M-1-005); Aiyou Foundation (KY201701); Ministry of Education Innovation Team Development Project (IRT-17R10); CAMS Key Lab of Translational Research on Lung Cancer (2018PT31035); China National Natural Sciences Foundation (81871889); Beijing Natural Science Foundation (7172045); Beijing Novel Program Grants, cross-cooperation (Z181100006218130); and Non-profit Central Research Institute Fund of the Chinese Academy of Medical Sciences (2018RC320009).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.