We recruited 395 NSCLC patients from the Department of Oncology in Affiliated Hospital of Chengde Medical University (Chengde, China) and the Department of Oncology in First Affiliated Hospital of Yangtze University (Jingzhou, China) between February 2017 and November 2019. The pathological diagnosis was verified by independent pulmonary pathologists based on the 4th edition of the World Health Organization Classification of Lung Tumors (41), and tumors with histological components other than NSCLC were excluded. Eight fresh tumor tissues, 213 formalin-fixed and paraffin-embedded (FFPE) tumor specimens, and 174 plasma were collected for NGS analysis. For the group of targeted therapy in patients with advanced-stage (IIIB and IV) NSCLC, 69 tumor tissues and 68 plasma samples were collected before the initiation of any therapies. Clinical characteristics of 395 NSCLC patients are shown in Table 1 .

Genomic DNA (gDNA) from the fresh tumor tissues, FFPE tumor specimens, and plasma was extracted by QIAamp DNA Mini Kit (Qiagen, Hilden, Germany), GeneRead DNA FFPE Kit (Qiagen, Hilden, Germany), and HiPure Circulating DNA Midi Kit C (Magen, Guangzhou, China), respectively. Qubit^® 3.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) and NanoDrop ND-1000 (Thermo Scientific, Wilmington, DE, USA) were used to assess the quantity and purity of gDNA. Fragmentation status was evaluated by the Agilent 2100 Bioanalyzer instrument (Agilent Technologies) using the High-Sensitivity DNA Reagent (Agilent Technologies, Santa Clara, CA, USA) to produce a DNA integrity number (DIN). Additionally, the step of quality control (QC) was also performed to assess FFPE DNA integrity by multiplex polymerase chain reaction (PCR). In brief, gDNA (30 ng) was amplified using three different primers of the glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene, the size set of which was 200–400 base pairs. An Agilent 2100 bioanalyzer instrument (Agilent Technologies) was used to determine the concentration of multiplex PCR products. The fragmentation of gDNA from FFPE was estimated by the average yield ratio (AYR) value, which was calculated by dividing the yield ratio of reference DNA (Promega Madison, WI, USA) into each amplicon.

Three hundred nanograms of gDNA from each sample determined by Qubit quantification was mechanically fragmented via an E220 focused ultrasonicator Covaris (Covaris, Woburn, MA, USA). The targeted size of the DNA fragment was between 150 and 200 bp. Then, 10–100 ng DNA was used for library construction by the KAPA library preparation kit (Kapa Biosystems Inc., Wilmington, MA, USA), which was constructed with end repair, A-tailing, and adapter ligation without additional fragmentation following the manufacturer’s instructions. Finally, the NGS libraries were captured using the xGen Lockdown Probe pool (IDT Technologies), and the captured DNA fragments were amplified with 12–13 cycles of PCR using 1× KAPA HiFi Hot Start Ready Mix.

After QC and quantification by Agilent 2100 Bioanalyzer (Agilent Technologies) and Qubit^® 3.0 Fluorometer (Invitrogen; Thermo Fisher Scientific, Inc.), the NGS libraries were sequenced on an Illumina NextSeq CN500 platform (Illumina Inc, San Diego, CA, USA) Medium flux chip.

Clean data were obtained following filtering adapter, low quality reads, and reads with length <36 bp, and they were aligned to reference human genome (University of California Santa Cruz ID: hg19) using the Burrows–Wheeler Aligner v. 0.7.12. Subsequently, the Picard and Genome Analysis Toolkit (GATK v.3.2) method was used for duplicate removal, local realignment, and base quality score recalibration, and it also generated the quality statistics, such as mapped reads, mean mapping quality, and mean coverage. Finally, VarDict was adopted for the identification of SNV and InDel.

The ANNOVAR software tool was used for annotating somatic variants. The candidate somatic variants were identified by the following filter conditions: (i) remove mutations with coverage depth <10×; (ii) remove variant sites with mutant allele frequency (MAF) >0.001 in the 1,000 Genomes databases (1,000 Genomes Project Consortium; https://www.internationalgenome.org/) and remove variant sites with MAF >0.001 in the ExAC (https://ncbiinsights.ncbi.nlm.nih.gov/tag/exac/); (iii) retain variant sites with MAF ≥0.001 and <0.1 in the 1,000 Genomes databases with COSMIC evidence; (iv) retain variations in the exonic or splicing region (10 bp upstream and downstream of splicing sites); (v) remove synonymous mutations; (vi) remove Unknown Variant_Classification; and (vii) the functional benign variant sites predicted by PolyPhen-2 and MutationTaster were removed.

CNVkit was employed to identify copy number variations (CNVs). Focal copy number events were determined with Genomic Identification of Significant Targets in Cancer (GISTIC) 2.0 (threshold: q-value <0.25). The correlation between the identified somatic variants and their clinical significance was established by OncoKB Precision Oncology Database (http://oncokb.org/).

The maftools package in R software (R 3.5.1, R Core Team; https://www.RProject.org) was used to create the mutational landscape, including somatic single nucleotide variations (SNVs) and InDels. The enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) were carried out to explore the biological significance of the candidate mutant genes by ClusterProfiler package (http://bioconductor.org/packages/release/bioc/html/clusterProfiler.Html) (42).

Tumor responses (CR, complete remission; PR, partial response; SD, stable disease; and PD, progressive disease) were assessed by the Response Evaluation Criteria in Solid Tumors (RECIST version 1.1) (43–46). The disease control rate (DCR) was defined as all positive responders (including complete, partial, and stable responders) divided by the total number of response-evaluable patients. The response rate (RR) was calculated by dividing complete and partial responders into the total number of response-evaluable patients. Meanwhile, the clinic pathological characteristics and DCR of targeted therapy were compared between two groups via chi-square tests by R software and GraphPad Prism (v. 7.0; GraphPad Software, La Jolla, CA) software, respectively. Fisher exact test was used to evaluate the statistical differences in categorical variables between two groups. A two-sided p-value <0.05 was used to evaluate statistical significance. The group of fewer than five samples was not adopted to perform statistical analyses. Tumor response assessments were only applied to patients with advanced-stage (IIIB and IV) NSCLC and receiving the first-line targeted therapy.

In this retrospective study, 395 NSCLC patients were enrolled including 127 patients from the Department of Oncology in Affiliated Hospital of Chengde Medical University and 268 patients from the Department of Oncology in First Affiliated Hospital of Yangtze University. Among them, 340 patients were diagnosed as LUAD, 54 patients were diagnosed as LUSC, and 1 patient was diagnosed as large cell carcinoma. In terms of tumor staging, the number of NSCLC patients in stages I–IV was 27, 14, 47, and 255, respectively. One hundred thirty-seven patients with advanced-stage (III and IV) NSCLC were treated with targeted therapy before any other therapies. For other clinical indicators, 273 of 395 cases had lymph node metastasis pathologically, 224 of 395 patients had never smoked, and 227 were male patients. Till the last time of our follow-up, 307 patients were still alive, 47 patients had died, and others were lost to follow-up ( Table 1 ).

To delineate the mutation landscape of NSCLC, we first analyzed somatic mutations from 221 tumor tissue samples by an NGS panel of 95 known cancer genes ( Supplementary Table S1 ). A total of 455 somatic variants in 46 genes were detected in 199 of 221 (90.05%) tumor tissues, including 39 nonsense mutations, 298 missense mutations, 7 frameshift insertions, 3 frameshift deletions, 17 in-frame insertions, 79 in-frame deletions, and 12 splice sites. The top 10 frequent mutated genes were TP53 (54.30%), EGFR (48.87%), RB1 (14.48%), KRAS (8.60%), PIK3CA (8.60%), CD3EAP (7.24%), CTNNB1 (5.43%), ERBB2 (5.43%), APC (4.07%), and BRAF (2.26%) ( Figure 1A ; Supplementary Table S2 ).

To explore the feasibility of genomic profiling by peripheral blood, the NGS panel was also applied to the analysis of 174 plasma samples. A total of 236 somatic variants in 33 genes were identified in 125 of 174 (71.84%) plasma-derived ctDNA samples, including 11 nonsense mutations, 157 missense mutations, 4 frameshift insertions, 1 frameshift deletion, 3 in-frame insertions, 51 in-frame deletions, and 9 splice sites, which was lower than that from tumor tissue DNA ( Figure 1B ; Supplementary Table S2 ). Expectedly, the mutation frequencies of the most 10 common mutated genes were also lower, which were TP53 (32.76%), EGFR (25.86%), CD3EAP (12.64%), RB1 (8.05%), PIK3CA (7.47%), KRAS (5.75%), APC (2.87%), TERT (2.87%), FBXW7 (2.30%), and HRAS (1.72%) ( Figure 1B ; Supplementary Table S2 ). In short, the percentage of ctDNA-specific mutated genes and tissue DNA-specific mutated genes were 6.12% (3/49) and 32.65% (16/49), respectively, while 61.22% (30/49) of mutated genes were shared by ctDNA analysis and tumor tissue DNA analysis ( Figure 1C ; Supplementary Figure S1 , Supplementary Table S2 ).

A total of 213 tumor tissue samples and 139 plasma samples were profiled for somatic copy number alterations (SCNAs) in NSCLC. There were 20 significant peaks of copy number gain in 213 tumor tissue samples including 1q21.3, 3q26.31, 3q26.33(PIK3CA), 5p15.33(TERT), 5p13.3, 7p22.2, 7p15.3, 7p11.2(EGFR), 7q31.2(TP53TG1), 8p11.22, 8q21.13, 8q24.21, 11q13.3, 12q14.3, 13q14.2(RB1), 13q34, 14q13.3, 17q12, 19q13.12, and 20q13.2 ( Figures 2A , left, 2B ; Supplementary Table S3 ), while there were only six significant peaks of copy number gain in 139 plasma samples including 1q23.3, 5p15.33(TERT), 5p14.1, 7p11.2(EGFR), 7q31.2, and 8q22.3(TP53INP1) ( Figures 2C , left, 2D ; Supplementary Table S3 ). Additionally, we also identified 15 significant peaks of copy number loss including 1p36.33, 1q23.1, 4p16.3, 5p15.33(TERT), 6q22.1, 7p22.3, 7q34, 8q24.3, 9q34.3, 10q11.21, 10q26.13, 11p15.5, 13q34, 16p13.3, and 20q13.33 in 213 tissue samples ( Figures 2A , right, 2B ; Supplementary Table S3 ), but only three significant peaks of copy number loss including 6q22.1, 7q34 and 10q26.13 were found in 139 plasma samples ( Figures 2C , right, 2D ; Supplementary Table S3 ). Lastly, SCNAs of LUAD ( Supplementary Figure S2 ) and LUSC ( Supplementary Figure S3 ) were also profiled, which further indicated that the feasibility of genomic profiling using ctDNA analysis was lower than using tumor tissue DNA analysis.

To compare the feasibility of genomic profiling using plasma samples in NSCLC subtypes, ctDNA from 143 LUAD patients and 31 LUSC patients was sequenced by the NGS panel of 95 known cancer genes.

In the group of LUAD, a total of 417 somatic variants in 44 genes were identified in 176 of 197 (89.34%) tumor tissue samples, including 35 nonsense mutations, 271 missense mutations, 7 frameshift insertions, 3 frameshift deletions, 17 in-frame insertions, 75 in-frame deletions, and 9 splice sites ( Figure 3A ; Supplementary Table S4 ). There were 197 somatic variants in 33 genes identified in 102 of 143 (71.33%) ctDNA samples including 8 nonsense mutations, 131 missense mutations, 2 frameshift insertions, 3 in-frame insertions, 46 in-frame deletions, and 7 splice sites ( Figure 3C ; Supplementary Table S4 ). Moreover, the mutation frequencies of target genes detected by ctDNA analysis were lower than that by tumor tissue DNA analysis ( Figures 3A, C ; Supplementary Table S4 ). In summary, the percentage of ctDNA-specific mutated genes and tumor tissue DNA-specific mutated genes were 6.38% (3/47) and 29.79% (14/47), respectively, while 63.83% (30/47) of mutated genes were detected by both ctDNA analysis and tumor tissue DNA analysis ( Figure 3 ; Supplementary Table S4 ).

In the group of LUSC, there were 36 somatic variants in 13 genes were identified in 22 of 23 (95.65%) tissue samples, including 3 nonsense mutations, 26 missense mutations, 4 in-frame deletions, and 3 splice sites ( Figure 3B ; Supplementary Table S4 ). Meanwhile, a total of 39 somatic variants in 9 genes were identified in 23 of 31 (74.19%) plasma samples including 3 nonsense mutations, 26 missense mutations, 5 in-frame deletions, 1 frameshift deletion, 2 frameshift insertions, and 2 splice sites ( Figure 3D ; Supplementary Table S4 ). Additionally, the mutation frequencies of target genes detected by ctDNA analysis were similar to that by tumor tissue DNA analysis ( Figures 3B, D ; Supplementary Table S4 ). Collectively, the percentage of ctDNA-specific mutated genes, tissue DNA-specific mutated genes, and mutated genes detected by both biopsies were 13.33% (2/15), 40% (6/15), and 46.67% (7/15), respectively ( Figure 3 ; Supplementary Table S4 ). All of the results were consistent with previous studies that the genetic heterogeneity of LUSC is more complex than that in LUAD (47, 48). In short, our data indicated that ctDNA analysis is more feasible as an alternative for genomic profiling of LUAD than that of LUSC by comparing the consistency between ctDNA analysis and tumor DNA analysis.

To expound on the biological function of mutated genes in NSCLC, KEGG and GO enrichment analyses were performed. We identified 104 altered signaling pathways in 199 tumor tissues samples and 109 altered signaling pathways in 125 plasma samples by KEGG enrichment analyses ( Figures 4A, B ). All the mutated genes involved in these pathways are shown in Table 2 and Supplementary Table S5 . Meanwhile, a total of 56 altered functional terms were enriched in 199 tumor tissue samples, and 41 altered functional terms were enriched in 125 plasma samples by GO enrichment analyses ( Table 3 ; Supplementary Table S5 ). Figures 4D, E show the top 10 most abundant altered functional terms according to gene counts and p-value ( Table 3 ). Additionally, the mutation types and distributions of all mutated genes involved in the top 10 altered functional terms in 199 tumor tissue samples and 125 plasma samples were shown in 28 and 19 lollipop plots ( Supplementary File 2 ), respectively. Collectively, we observed a high level of consistency of altered signaling pathways and altered functional terms between ctDNA analysis and tumor tissue DNA analysis, which were 83.62% (97/116) and 59.02% (36/61), respectively ( Figure 4C ).

Additionally, Figure 5 shows the top 10 altered signaling pathways and altered functional terms in LUAD ( Figures 5A–F ; Supplementary Table S6 ) and LUSC ( Figures 5G – L ; Supplementary Table S6 ) enriched by KEGG enrichment analyses ( Figures 5A, B, G, H ) and GO enrichment analyses ( Figures 5D, E, J, K ) using tumor tissues ( Figures 5A, D, G, J ) and ctDNA ( Figures 5B, E, H, K ). A high level of consistency of altered signaling pathways between ctDNA analysis and tumor DNA analysis was observed in both LUAD and LUSC, which were 82.91% (97/117) versus 89.01% (81/91) ( Figures 5C, I ). However, we found about 50% consistency of altered functional terms between these two kinds of samples in both LUAD and LUSC, which were 55.07% (38/69) versus 46.51% (20/43) ( Figures 5F, L ). Among these significantly altered functional terms, we focused on the alteration in transmembrane receptor protein tyrosine kinases ( Figures 5D, E, J, K ), which play key roles in a wide range of cellular processes including growth, motility, differentiation, apoptosis, and metabolism (49, 50). As in previous studies, the abnormity of receptor tyrosine kinases (RTKs) signaling is closely related to diseases, especially malignancy (51–54). Constitutive activation, such as gain-of-function mutations, genomic amplification, chromosomal rearrangements, and/or autocrine activation, serves as excellent examples to show the oncogenic properties of RTKs (49, 55).

To assess the clinical utility of anticipative molecular profiling, all mutations were divided into different levels according to the evidence of clinical actionability in OncoKB ( Supplementary Figure S4A ). As standard therapeutic biomarkers, a cluster of gene mutations was approved by the FDA. In our cohort, 48.35% (191/395) of patients possess at least one actionable alteration. In the group of NSCLC, level_1 accounted for 79.21% (343/433), including missense mutation of EGFR, PIK3CA, BRAF, IDH2, IDH1, ATM, and KRAS, in-frame insertion of EGFR, and in-frame deletion of EGFR; level_2 accounted for 16.17% (70/433), including in-frame insertion of ERBB2 and missense mutation of ERBB2, PIK3CA, KRAS, BRAF, and NRAS; level_3 accounted for 1.62% (7/433), including in-frame deletion of EGFR and missense mutation of AKT1; level_4 accounted for 3.00% (13/433), including missense mutation of EGFR, FGFR3, and mTOR, and nonsense mutation of CDKN2A ( Table 4 ; Supplementary Figures S4B, E, H ; Supplementary Table S7 ). However, patients with LUAD have a higher percentage of actionable alterations than patients with LUSC, which were 52.06% (177/340) versus 25.93% (14/54) ( Supplementary Figures S4C, D, F, G, I, J ; Supplementary Table S7 ), indicating that they may benefit more from targeted therapy.

Next, we wanted to explore the impact of genomic characteristics on the clinical outcomes of patients with advanced-stage (III and IV) NSCLC. Based on EGFR mutations, 46 patients were classified as the group of EGFR alterations without TP53 mutations (called EGFR), and 68 patients were classified as the group of EGFR and TP53 co-alterations (called EGFR&TP53). Among them, 13 out of 46 patients in EGFR group and 19 out of 68 patients in the EGFR&TP53 group were first treated with Gefitinib. As shown in Figure 6 , the disease control rate (DCR = CR + PR + SD) in patients with EGFR alterations alone and in patients with EGFR&TP53 co-alterations was 100% (9/9) and 68.75% (11/16), respectively (p = 0.0608) ( Figure 6A ; Table 5 ). Meanwhile, the response rate (RR = CR + PR) was 44.44% (4/9) in patients with EGFR alterations alone, while only 12.5% (2/16) in patients with EGFR&TP53 co-alterations (p = 0.0726) ( Figure 6A ; Table 5 ). Importantly, the distributions of drug-resistant EGFR mutations were also explored, and the percentage of T790M mutation in patients with PR, SD, and PD was 16.67% (1/6), 14.29% (2/14), and 20% (1/5), respectively ( Table 5 ; Supplementary Table S8 ). Additionally, we also sought to investigate whether the TP53 concomitant mutations are associated with the responses to different first-generation EGFR-TKIs. As shown in Figure 6B and Table 5 , the DCR of Gefitinib was similar to that of Icotinib (68.75% vs 69.23%, p > 0.9778).

Co-occurring genetic alterations might play important roles in the mechanisms of tumor responses, and it may be helpful to explain the marked diversity of individual outcomes (32, 33). As previous studies have reported, a potential role of TP53 mutation is associated with poor therapeutic responses (34–36). Figure 7 shows mutations and protein structures of p53 and EGFR in NSCLC patients who were first treated with Gefitinib, Icotinib, or either of them. TP53 mutations were detected in 11 patients who responded well to Gefitinib, including nonsense mutations (numbered BT1812060034LNCTX), and missense mutations (numbered BT1806210029LNCBP, BT1806220195LNCTV, BT1806290337LNCTB, BT1808040195LNCTV, BT1808230036LNCBP, BT1810270069LNCBP, BT1812020489LNCTS, BT1902210293LNCTV, BT1905230269LNCTV, and BT1906210049LNCBP) ( Figure 7A , up). However, TP53 mutation sites and types in five patients with PD were different from those in the above patients, including nonsense mutations (numbered BT1812160473LNCTV and BT1901200062LNCTB), missense mutations (numbered BT1902140231LNCTX and BT1905110010LNCBP), and splice site (numbered BT1805130336LNCTV) ( Figure 7A , down). Furthermore, we also carried out the comparative analysis of p53 mutations for patients treated with Icotinib ( Figure 7B ; Supplementary Table S8 ) because it has similar chemical structures, molecular mechanisms, and clinical curative effects with Gefitinib but cost a much lower price. Our results also indicated that patients who responded differently to Icotinib carried different mutation sites. However, the roles of TP53 mutations on clinical outcomes in NSCLC patients with EGFR-mutation treated with Gefitinib and Icotinib require further elucidation, since the molecular mechanism among them remains largely unknown.