In total, 120 patients with NSCLC from the Department of Pathology at the Affiliated 3201 Hospital of Xi'an Jiaotong University were recruited between May, 2017 and December, 2020. The pathological diagnosis was verified by three pulmonary pathologists, based on the 4th edition of the World Health Organization Classification of Lung Tumors (24). Tumors with histological components other than NSCLC were excluded. In total, one fresh tumor tissue, 68 formalin-fixed and paraffin-embedded (FFPE) tumor specimens and 51 peripheral blood specimens were collected for next-generation sequencing (NGS) analysis. Written informed consent was acquired from all participating individuals. The study was performed according to the Code of Ethics of the World Medical Association (Declaration of Helsinki) (25), and it was approved by the Medical Ethics Committee of Affiliated 3201 Hospital of Xi'an Jiaotong University. The corresponding Institutional Review Board (IRB) number was No.008(2017).

Genomic DNA (gDNA) from fresh tumor tissues, FFPE tumor specimens and plasma was extracted using the QIAamp DNA Mini kit (Qiagen GmbH), GeneRead DNA FFPE kit (Qiagen GmbH) and the HiPure Circulating DNA Midi kit C (Magen Biotechnology Co., Ltd.), respectively. The Qubit^® 3.0 Fluorometer (Invitrogen; Thermo Fisher Scientific, Inc.) and NanoDrop ND-1000 (Thermo Scientific, Inc.) were used for the quantity and purity evaluation of the gDNA. The fragmentation status was evaluated using the Agilent 2100 Bioanalyzer instrument (Agilent Technologies, Inc.) with the High Sensitivity DNA Reagent (Agilent Technologies, Inc.) to produce a DNA integrity number. Additionally, the step of quality control (QC) was also performed to evaluate 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, sized at 200-400 base pairs (Table SI). The Agilent 2100 Bioanalyzer instrument (Agilent Technologies, Inc.) 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 Corporation) by each amplicon's yield ratio.

In total, 300 ng gDNA from each sample was mechanically fragmented using an E220 focused ultrasonicator Covaris (Covaris, LLC.). The targeted size of the DNA fragment was between 150 and 200 bp. Subsequently, 10-100 ng DNA was used for library construction with the KAPA library preparation kit (Kapa Biosystems Inc.; Roche Diagnostics), which was constructed with end-repair, A-tailing and adapter ligation without additional fragmentation, according to the manufacturer's instructions. Finally, the NGS libraries were captured using the xGen Lockdown Probe pool (Integrated DNA Technologies, Inc.), and the captured DNA fragments were amplified for 13 cycles of PCR, using 1X KAPA HiFi Hot Start Ready Mix (Kapa Biosystems Inc.; Roche Diagnostics). The thermal cycling conditions used were as follows: Initial denaturation for 45 sec at 98°C followed by 13 cycles of 98°C for 15 sec, 60°C for 30 sec, and 72°C for 30 sec. The final extension was performed at 72°C for 60 sec.

Following QC and quantification using the Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.) and Qubit^® 3.0 Fluorometer (Invitrogen; Thermo Fisher Scientific, Inc.), the NGS libraries were sequenced on an Illumina NextSeq CN500 platform with a medium flux chip (NextSeq CN500 Mid Output v2 kit; Illumina Inc.).

Clean data were obtained following filtering low-quality reads, which includes reads with adapter sequences and reads with length <36 bp. All filtered reads were aligned to the human genome (University of California Santa Cruz ID: hg19) using the Burrows-Wheeler-Alignment Tool (BWA v.0.7.12; Wellcome Trust Sanger Institute) (26). Subsequently, the Picard and Genome Analysis Toolkit (GATK v.3.2; Broad Institute) method was used for duplicate removal, local realignment, and base quality score recalibration, and it was also adopted for generating quality statistics, including mapped reads, mean mapping quality, and mean coverage. Finally, the VarDict (v.1.6.0; GitHub, Inc.) was adopted for the identification of single nucleotide variation (SNV) and Insertion/Deletion (InDel) (27).

The ANNOVAR software tool (v. 20210202; https://annovar.openbioinformatics.org/en/latest/) was used for annotating somatic variants (28). The candidates of somatic variants were identified by the following filter conditions: i) Removal of the variants coverage depth (VDP) <10; ii) removal of the 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 Exome Aggregation Consortium (ExAC) (https://ncbiinsights.ncbi.nlm.nih.gov/tag/exac/); iii) retainment of variant sites with MAF≥0.001 and <0.1 in the 1,000 Genomes databases with COSMIC evidence (http://cancer.sanger.ac.uk/cosmic); iv) retainment of variations in the exon or splicing region (10 bp upstream and downstream of splicing sites); v) remove synonymous mutations; vi) remove unknown variant classification; and vii) removal of the functional benign variant sites, predicted by PolyPhen-2 (Polymorphism Phenotyping v2; http://genetics.bwh.harvard.edu/pph2/) (29) or MutationTaster (MutationTaster2020; https://www.muta-tiontaster.org/) (30). Additionally, the association between the identified somatic mutations and their clinical significance was established by OncoKB Precision Oncology Database (http://oncokb.org/). Kyoto Encyclopedia Of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analysis were used to explore the biological consequences of mutant genes by using the cluster Profiler package (http://bioconductor.org/packages/release/bioc/html/clusterProfiler.Html) in R software (R 4.0.3, R Core Team; https://www.RProject.org) (31).

The maftools package in R software (R 4.0.3, R Core Team;https://www.R-Project.org) was used to create somatic mutation landscapes, co-Barplots, co-Oncoplots, lollipop plots, and spectrums of co-occurring and mutually exclusive genomic alterations. Fisher's exact test was used to evaluate the statistical differences in categorical variables between the BM group and the non-BM group using R software (R 4.0.3, R Core Team; https://www.RProject.org). Statistical analyses were two-sided, and P<0.05 was considered to indicate a statistically significant difference.

In the present retrospective study, the numbers of patients with different stages of NSCLC were as follows: I, II, III and IV were as follows: i) Stage I, 6 patients (aged 61-78 years; median, 76 years); ii) stage II, 3 patients (aged 57-71 years; median, 58 years); iii) stage III, 21 patients (aged 43-77 years; median, 61 years); and iv) stage IV, 90 patients (aged 30-84 years; median, 63 years). Patients with stage IV NSCLC (n=90) were divided into the BM group (n=26, aged 43-84 years) and non-BM group (n=64, aged 30-84 years) by clinically detectable metastatic lesions. There were no significant differences in the demographic and clinical characteristics between the BM and the non-BM group (Table I).

To delineate the somatic mutation landscape, the somatic mutations from the tumor tissue DNA of 69 patients were first analyzed by applying an NGS panel of 95 known cancer genes (Table SII). The present study mainly focused on protein-altering variants, based on the annotation of somatic SNVs and InDels. A total of 157 somatic variants of 32 mutated genes were detected in 62 out of 69 (89.86%) tumor tissues (Fig. 1A and Table SIII). To explore the feasibility of genomic profiling using peripheral blood samples, the same NGS panel was used to detect 51 plasma-derived ctDNA. In total, 77 somatic variants of 30 mutated genes were identified in 38 out of 51 (74.51%) peripheral blood samples, which was lower than that from tissue DNA (Fig. 1B and Table SIII). Furthermore, the percentages of tissue DNA-specific mutated genes and ctDNA-specific mutated genes were 28.57% (12/42) and 23.81% (10/42), respectively, whereas the consistency rate of mutated genes was 47.62% (20/42) in NSCLC (Fig. 1C and Table SIII).

Patients with NSCLC with BM have been shown to be significantly associated with an increased mortality rate (4). However, the underlying molecular mechanisms of the initiation and progression of BM have not yet been fully elucidated in NSCLC. In the present study, to explore this issue, the genomic characterizations of patients both from the BM and the non-BM group were first depicted using an NGS panel. In total, 47 somatic variants of 17 mutated genes were identified in 23 out of 26 (88.46%) patients with BM (Figs. 2 and S1A, and Table SIV), and 127 somatic variants of 30 mutated genes were observed in 55 of 64 (85.94%) non-BM patients (Figs. 2 and S1C, and Table SIV). To further discover the molecular differences between the BM and the non-BM group, a Venn diagram was plotted (Fig. S1B). Among the total number of 35 mutant genes, ALK, cytidine deaminase (CDA), SMAD family member 4 (SMAD4), superoxide dismutase 2 (SOD2), and Von Hippel-Lindau tumor suppressor (VHL) were uniquely present in the BM group, whereas 18 mutant genes [ATM serine/threonine kinase (ATM), cyclin-dependent kinase inhibitor 2A (CDKN2A), Erb-B2 receptor tyrosine kinase 4 (ERBB4), HRAS, PIK3CA, etc.] only emerged in the non-BM group. In total, 12 mutant genes [APC regulator of WNT signaling pathway (APC), DNA-directed RNA polymerase I subunit RPA34 (CD3EAP), catenin beta-1 (CTNNB1), cytochrome P450 2D6 (CYP2D6), dihydropyrimidine dehydrogenase (DPYD), EGFR, Kirsten rat sarcoma (KRAS), NRAS, TP53, telomerase reverse transcriptase (TERT), thiopurine S-methyltransferase (TPMT) and RB transcriptional co-repressor 1 (RB1)] were concurrently present in both groups.

Moreover, the frequencies and sites of the above concurrently mutated genes were also investigated. Of note, apart from CD3EAP, DPYD and TPMT, the frequencies and/or sites of the concurrently mutated genes were different (Fig. 3 and Table SV). In relation to the KRAS gene, marked differences were observed in the mutation sites and frequencies (p.L19F vs. p.G12D, p.G12V, p.Q61H, and p.A146T, 3.85 vs. 12.5%) of these two groups (Fig. 3A and Table SV). Although about half of the patients harbored TP53 mutations in the two groups, the types and sites of amino acid alterations were different (Fig. 3C and Table SV). The frameshift insertion of TP53 was exclusively identified in the BM group, and in-frame deletion and nonsense mutation were only detected in the non-BM group.

In NSCLC, patients presenting with EGFR mutations have a much higher incidence of BM compared to those without EGFR mutation. EGFR mutations and KRAS mutations are usually mutually exclusive, and KRAS mutations could confer resistance to EGFR-TKIs when they co-exist (32). In the present study, the somatic interactions in the BM group were different from the interactions of the non-BM group. EGFR and KRAS were a mutually exclusive set of genes in the non-BM group (P=0.0183) (Fig. 4B and Table SV), whereas limited mutually exclusive interactions were observed in the BM group (P=0.999) (Fig. 4A and Table SVI). Additionally, obvious differences were also detected in the co-occurring set of genes between the two groups. As demonstrated in Fig. 4, VHL and ALK (P=0.0435), SOD2 and KRAS (P=0.0435), TERT and KRAS (P=0.0435), TERT and SOD2 (P=0.0435) were co-occurring pair of genes in the BM group (Fig. 4A and Table SV), while FGFR3 and CD3EAP (P=0.0141), isocitrate dehydrogenase [NADP(+)] 1 (IDH1) and BRAF (P=0.0182), neurotrophic receptor tyrosine kinase 1 (NTRK1) and KRAS (P=0.0189), PIK3CA and KRAS (P=0.0340), CDKN2A and BRAF (P=0.0364), and interactions between other six co-occurring pair of genes were detected in the BM group (Fig. 4B and Table SVI).

In order to acquire a more incisive understanding of the biological consequences in these two groups, KEGG and GO enrichment analyses were performed. The top 10 KEGG pathways enriched by mutated genes were depicted according to gene count and P-value, and most signaling pathways were cancer-related (Fig. 5A and B, and Table II). Since patients with BM respond well to EGFR-TKIs, the P-value of EGFR-TKI resistance was markedly higher in the BM group, than in the non-BM group (4.50E-04 vs. 7.84E-18) (Tables II and SVII). In GO enrichment analysis, functional categories were most involved in RNA polymerase II transcription regulator complex in the BM group (Fig. 5D and Table III) and promyelocytic leukemia nuclear body in the non-BM group (Fig. 5E and Table III). Furthermore, five mutated genes uniquely present in the BM group were prevailingly distributed in transcription regulator complex, RNA polymerase complex, bicellular tight junction, tight junction, apical junction complex, and so forth (Tables III and SVII). Collectively, a high consistency of altered signaling pathways between these two groups according to KEGG analysis was observed (Fig. 5C), whereas the percentage was decreased in GO analysis-related altered functional terms [Fig. 5F; KEGG analysis, 77.27% (85/110) vs. GO analysis, 30.00% (15/50)].

In order to evaluate the clinical utility of anticipative molecular profiling, all mutations were divided into different levels, according to the evidence of clinical actionability in OncoKB (Fig. 6A). As standard therapeutic biomarkers, a cluster of gene mutations was approved by the FDA. In the present cohort, 47 out of 120 (39.17%) patients possessed at least one actionable alteration. Among the patients with stage IV disease, level_1 accounted for 34.44% (31/90), including missense mutations of BRAF, EGFR, IDH1, PDGFRA and PIK3CA, a nonsense mutation of ATM and an in-frame insertion of EGFR; level_2 accounted for 4.44% (4/90), including missense mutations of NRAS and PIK3CA and an in-frame deletion of PIK3CA; level_3 accounted for 1.11% (1/90), including an in-frame insertion of EGFR; level_4 accounted for 1.11% (1/90), including missense mutations of CDKN2A (Fig. 6B, E and H, and Table SVIII). Additionally, it was also observed that non-BM patients had a slightly higher percentage of actionable alterations than patients with BM, namely 45.31% (29/64) vs. 30.77% (8/26) (Fig. 6D, F, G, I, J, and Table SVIII).

To compare the feasibility of genomic profiling of advanced patients with or without BM using plasma-derived ctDNA, somatic mutations from 14 tumor tissues and 12 peripheral blood samples were analyzed in the BM group by the above NGS panel. A total of 32 somatic variants of 12 mutated genes were identified in 13 out of 14 (92.86%) tumor tissue DNA samples (Fig. 7A and Table SIX), and 15 somatic variants of 10 mutant genes were detected in 10 of 12 (83.33%) plasma-derived ctDNA (Fig. 7B and Table SIX). Meanwhile, eighty-three somatic variants of 22 mutated genes in 36 of 39 (92.31%) tumor tissue DNA were also detected (Fig. 7C and Table SX), as well as 44 somatic variants of 21 mutated genes in 19 out of 25 (76.00%) plasma-derived ctDNA (Fig. 7D and Table SX) in the non-BM group. In summary, 43.33% (13/30) of the mutated genes were detected by both tumor tissue DNA analysis and ctDNA analysis in the non-BM group (Fig. 7E), whereas the percentage was only 29.41% (5/17) in the BM group (Fig. 7E).