Differences in actionable genomic alterations between brain metastases and non‑brain metastases in patients with non‑small cell lung cancer.

Brain metastases (BM) have been closely associated with increased morbidity and poor survival outcomes in patients with non‑small cell lung cancer (NSCLC). Excluding risk factors in histological subtypes, genomic alterations, including epidermal growth factor receptor mutations and anaplastic lymphoma kinase (ALK) rearrangements have been also regarded as greater risk factors for BM in the aspect of molecular subtypes. In the present study, 69 tumor tissues and 51 peripheral blood samples from patients with NSCLC were analyzed using a hybridization capture‑based next‑generation sequencing (NGS) panel, including 95 known cancer genes. Among the 90 patients with stage IV NSCLC, 26 cases suffered from BM and 64 cases did not. In total, 174 somatic mutations in 35 mutated genes were identified, and 12 of these genes were concurrently present in the BM group and the non‑BM group. Importantly, five mutated genes including ALK, cytidine deaminase (CDA), SMAD family member 4 (SMAD4), superoxide dismutase 2 (SOD2) and Von Hippel‑Lindau tumor suppressor (VHL) genes were uniquely detected in the BM group, and they were enriched in the Hippo signaling pathway, pyrimidine metabolism and pantothenate and co‑enzyme A (CoA) biosynthesis, as demonstrated using Kyoto Encyclopedia of Genes and Genomes enrichment analysis. RNA polymerase II transcription regulator complex and promyelocytic leukemia nuclear body were the top functional categories according to the Gene Ontology enrichment analysis in the BM group and non‑BM group, respectively. Furthermore, 43.33% (13/30) of mutated genes were detected by both tumor tissue deoxyribonucleic acid (DNA) and plasma‑derived circulating tumor DNA (ctDNA) in the non‑BM group, while this percentage was only limited to 29.41% (5/17) in the BM group. To summarize, significant differences in somatic mutations, somatic interactions, key signaling pathways, functional biological information, and clinical actionability for the therapy of targeted agents were founded between the BM group and the non‑BM group, and ctDNA analysis may by applied as a more credible alternative for genomic profiling in patients with advanced NSCLC without BM, due to its higher consistency for genomic profiling between ctDNA analysis and tissue DNA analysis.

Introduction

Brain metastases (BM) can be detected in 10-20% of patients with non-small cell lung cancer (NSCLC) at the initial diagnosis and the percentage extends up to 50% during disease progression (1-3). The development of BM leads to an evident incidence of neurocognitive and functional deficits in patients with NSCLC, which adversely affects the quality of life (4). BM is closely related to a poor prognosis, and the median overall survival (OS) of patients with NSCLC with BM is limited only up to 4-6 months (5,6). However, the molecular mechanisms underlying the high morbidity, poor prognosis, and high mortality rate of patients with NSCLC with BM remain largely unknown.

In terms of histological subtypes, BM occurs most frequently in patients with lung adenocarcinoma (54 and 58.6%), as compared to non-adenocarcinoma subtypes (large cell carcinoma, 17.7%; and squamous cell carcinoma, 9.9%) (7,8). Moreover, patients with NSCLC harboring epidermal growth factor receptor (EGFR) mutations or anaplastic lymphoma kinase (ALK) rearrangements are at a greater risk of developing BM, in the aspect of molecular subtypes (9,10). Pathogenic driver mutations have been reported to be vital for therapeutic decision-making, since targeted therapies significantly improve survival outcomes for the majority of patients (11). Specifically, a subgroup of patients with BM responds well to EGFR tyrosine kinase inhibitors (EGFR-TKIs) (12). However, actionable genomic alterations are not currently available for recurrent or progressive diseases, due to the unsafety of the invasive biopsy procedure and/or the inaccessibility of tumor sites, particularly for patients with BM. Currently, knowledge about pathogenic driver mutations is most frequently obtained from primary or metastatic tumor tissue at easily accessible sites; however, a divergent mutation landscape from these sites has been observed (13,14). Thus, safe, convenient and replaceable approaches are essential to acquiring genomic characterization for patients with NSCLC with BM.

Liquid biopsy is regarded as a promising minimally-invasive approach to obtaining tumor cells, as opposed to invasive tissue biopsies in molecular diagnostics in recent years (15). Cell-free deoxyribonucleic acid (cfDNA), which has been previously detected in blood and body fluids, is released from cells when they undergo necrosis, apoptosis and lysis (16,17). Circulating tumor DNA (ctDNA), the tumor-derived fraction of cfDNA, possesses the potential ability to represent the whole tumor burden across different metastatic sites, partly circumventing tumor spatial heterogeneity issues (18). Since the reliability of a single-tumor tissue biopsy for the obtainment of the whole mutation landscape, the limitations of personalized medicine approaches are a dilemma. Multiregional biopsy analysis has already been proposed by Gerlinger et al (19), in order to profile a more complete mutation landscape and predict the therapeutic outcome. Additionally, several studies have suggested that a dynamic sampling of somatic mutations from ctDNA analysis may represent a larger clonal hierarchy (20-23), rendering the realization of therapeutic decisions and tracking therapeutic outcomes safer and more convenient, even across different metastatic sites.

In the present study, targeted next-generation sequencing (NGS) of tumor tissues and peripheral blood samples from patients with NSCLC was conducted, using a hybridization capture-based panel consisting of 95 known cancer genes. The present study aimed to identify the characteristic of genomic alternations in patients with advanced NSCLC with or without BM. The novelty of the present study was, to the best of our knowledge, systematical comparisons of the differences in somatic mutations, somatic interactions, key signaling pathways, functional biological information and clinical actionability for the therapy of targeted agents between the BM group and the non-BM group. The findings presented herein may possibly aid towards the elucidation of the underlying molecular mechanisms underlying the initiation or progression of BM in NSCLC. Furthermore, it was also investigated whether ctDNA analysis may serve as a more credible alternative for genomic profiling in patients with advanced NSCLC with BM than in patients without BM.

Patients and methods

Patients and sample collection

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).

DNA extraction and quality control

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.

Library preparation and hybridization capture

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.

Illumina sequencing

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.).

Bioinformatics analysis

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).

Statistical analysis

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.

Results

Patient characteristics

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).

Table I

Characteristics of patients with advanced NSCLC according to brain metastatic progression.

Clinical characteristics	No. of patients	Brain metastases		P-value
		Yes	No
Total sample	90	26	64	0.6378
Tissue	53	14	39
ctDNA	37	12	25
Age, years				0.9074
≤50	18	5	13
>50	72	21	51
Sex				0.8135
Male	54	15	39
Female	36	11	25
Histology				0.999
Adenocarcinoma	80	24	56
Squamous carcinoma	9	2	7
Adenosquamous carcinoma	1	0	1
Family history of lung cancer				0.999
Yes	1	0	1
No	89	26	63
History of pulmonary infection				0.8931
Yes	44	13	31
No	46	13	33
History of smoking				0.6771
Once	12	2	10
Now	26	8	18
Never	52	16	36
History of alcohol consumption				0.9312
Once	8	2	6
Now	19	6	13
Never	63	18	45
Pre-existing metabolic disease				0.1556
Yes	19	3	16
No	58	20	38
Unknown	13	3	10
Vascular invasion				0.7495
Yes	76	23	53
No	14	3	11
Nerve invasion				0.5002
Yes	11	2	9
No	79	24	55
CYFRA21-1 at baseline				0.5084
Normal	13	6	7
Elevated	26	9	17
Unknown	51	11	40
CEA at baseline				0.603
Normal	24	6	18
Elevated	60	20	40
Unknown	6	0	6
NSE at baseline				0.7294
Normal	28	12	16
Elevated	12	4	8
Unknown	50	10	40
SCC at baseline				0.1427
Normal	24	11	13
Elevated	12	2	10
Unknown	54	13	41

NSCLC, non-small cell lung cancer; ctDNA, circulating tumor DNA; CYFRA21-1, cytokeratin 19 fragment; CEA, carcinoembryonic antigen; NSE, neuron specific enolase; SCC, squamous cell carcinoma antigen.

Landscape of somatic mutations

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).

Figure 1

The landscape of somatic mutation in patients with NSCLC from (A) tumor tissue DNA (n=69) and (B) ctDNA (n=51). (C) Venn diagram of mutated genes derived from tissue DNA analysis or ctDNA analysis. Mutated genes are ranked by mutation frequency. The bars on the right represent the mutation frequencies of each gene. TMB (mutations per Mb) is demonstrated in the upper panel. Patients were arranged along the x-axis. NSCLC, non-small cell lung cancer; ctDNA, circulating tumor DNA; TMB, tumor mutation burden.

Differences in somatic mutations between the BM and non-BM group

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.

Figure 2

(A) coBarplot diagram of genomic subtyping with somatic mutations and (B) coOncoplot diagram of mutated genes in advanced NSCLC patients with (left panel) and without (right panel) BM. NSCLC, non-small cell lung cancer; BM, brain metastases.

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.

Figure 3

The protein structure and mutational proportion of (A) KRAS, (B) RB1, (C) TP53, (D) CD3EAP, (E) EGFR, and (F) CTNNB1 in patients with advanced NSCLC with (top panels) and without BM (bottom panels). Protein domains are marked in different colors. Lollipops represent the locations of protein-altering variants. Square brackets indicate the proportion of patient-harbored non-synonymous mutations in each group. CD3EAP, DNA-directed RNA polymerase I subunit RPA34; CTNNB1, catenin beta-1; NSCLC, non-small cell lung cancer; BM, brain metastases.

Differences in somatic interactions between the BM and 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).

Figure 4

Spectrum of co-occurring and mutually exclusive genomic alterations in patients with advanced NSCLC (A) with and (B) without BM. NSCLC, non-small cell lung cancer; BM, brain metastases.

Key signaling pathways and biological functions analysis of somatic mutations in the BM and the non-BM group

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)].

Figure 5

(A and B) Signaling pathway analysis by KEGG and (D and E) functional terms by GO enrichment analysis in patients with advanced NSCLC (A and D) with and (B and E) without BM. The value represents the number of mutated genes enriched in these functional terms. Venn diagram of (C) signaling pathways or (F) functional terms in the BM group and the non-BM group. KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; NSCLC, non-small cell lung cancer; BM, brain metastases.

Table II

Signaling pathways enriched by KEGG analysis in the BM group and the non-BM group.

ID	Description	Adjusted P-value	geneID
BM group
hsa05226	Gastric cancer	1.85E-12	TP53/TERT/KRAS/RB1/EGFR/SMAD4/CTNNB1/APC/NRAS
hsa05225	Hepatocellular carcinoma	5.53E-12	TP53/TERT/KRAS/RB1/EGFR/SMAD4/CTNNB1/APC/NRAS
hsa05210	Colorectal cancer	1.21E-10	TP53/KRAS/EGFR/SMAD4/CTNNB1/APC/NRAS
hsa05213	Endometrial cancer	7.60E-10	TP53/KRAS/EGFR/CTNNB1/APC/NRAS
hsa05223	Non-small cell lung cancer	2.89E-09	TP53/KRAS/RB1/EGFR/ALK/NRAS
hsa05224	Breast cancer	5.41E-09	TP53/KRAS/RB1/EGFR/CTNNB1/APC/NRAS
hsa05219	Bladder cancer	1.06E-08	TP53/KRAS/RB1/EGFR/NRAS
hsa05215	Prostate cancer	1.78E-08	TP53/KRAS/RB1/EGFR/CTNNB1/NRAS
hsa01522	Endocrine resistance	1.89E-08	TP53/KRAS/RB1/EGFR/CYP2D6/NRAS
hsa05165	Human papillomavirus infection	6.64E-08	TP53/TERT/KRAS/RB1/EGFR/CTNNB1/APC/NRAS
Non-BM group
hsa05230	Central carbon metabolism in cancer	1.47E-20	PDGFRA/NRAS/TP53/PIK3CA/KRAS/EGFR/HR AS/MTOR/KIT/FGFR3/NTRK1/IDH1/ERBB2
hsa05219	Bladder cancer	2.47E-19	NRAS/TP53/RB1/KRAS/EGFR/HRAS/CDKN2A/FGFR3/SRC/BRAF/ERBB2
hsa01522	Endocrine resistance	1.58E-18	CYP2D6/NRAS/TP53/PIK3CA/RB1/KRAS/EGFR/HRAS/CDKN2A/MTOR/SRC/BRAF/ERBB2
hsa01521	EGFR tyrosine kinase inhibitor resistance	7.84E-18	PDGFRA/NRAS/PIK3CA/KRAS/EGFR/HRAS/JAK2/MTOR/FGFR3/SRC/BRAF/ERBB2
hsa05215	Prostate cancer	1.06E-16	PDGFRA/NRAS/TP53/PIK3CA/RB1/KRAS/EGFR/HRAS/CTNNB1/MTOR/BRAF/ERBB2
hsa05214	Glioma	3.59E-16	PDGFRA/NRAS/TP53/PIK3CA/RB1/KRAS/EGFR/HRAS/CDKN2A/MTOR/BRAF
hsa05224	Breast cancer	3.76E-16	NRAS/TP53/PIK3CA/RB1/KRAS/EGFR/HRAS/CTNNB1/APC/MTOR/KIT/BRAF/ERBB2
hsa05226	Gastric cancer	4.50E-16	NRAS/TP53/PIK3CA/RB1/KRAS/EGFR/HRAS/CTNNB1/APC/MTOR/BRAF/TERT/ERBB2
hsa05213	Endometrial cancer	1.77E-15	NRAS/TP53/PIK3CA/KRAS/EGFR/HRAS/CTNNB1/APC/BRAF/ERBB2
hsa05225	Hepatocellular carcinoma	2.21E-15	NRAS/TP53/PIK3CA/RB1/KRAS/EGFR/HRAS/CDKN2A/CTNNB1/APC/MTOR/BRAF/TERT

Table III

Functional terms enriched by GO enrichment analysis in the BM and the non-BM group.

ID	Description	Adjusted P-value	geneID
BM group
GO:0090575	RNA polymerase II transcription regulator complex	0.0010	TP53/RB1/SMAD4/CTNNB1
GO:0030877	Beta-catenin destruction complex	0.0017	CTNNB1/APC
GO:1990909	Wnt signalosome	0.0017	CTNNB1/APC
GO:0016605	PML body	0.0023	TP53/TERT/RB1
GO:0016342	Catenin complex	0.0060	CTNNB1/APC
GO:0005667	Transcription regulator complex	0.0060	TP53/RB1/SMAD4/CTNNB1
GO:0019897	Extrinsic component of plasma membrane	0.0060	KRAS/CTNNB1/APC
GO:0061695	Transferase complex, transferring phosphorus-containing groups	0.0163	TP53/TERT/RB1
GO:0016328	Lateral plasma membrane	0.0173	CTNNB1/APC
GO:0019898	Extrinsic component of membrane	0.0225	KRAS/CTNNB1/APC
Non-BM group
GO:0016605	PML body	0.0018	TP53/RB1/MTOR/TERT
GO:0005925	Focal adhesion	0.0018	KRAS/EGFR/JAK2/CTNNB1/FGFR3/SRC
GO:0030055	Cell-substrate junction	0.0018	KRAS/EGFR/JAK2/CTNNB1/FGFR3/SRC
GO:0009925	Basal plasma membrane	0.0021	ERBB4/EGFR/ERBB2
GO:0019898	Extrinsic component of membrane	0.0021	PIK3CA/KRAS/CTNNB1/APC/SRC
GO:0045121	Membrane raft	0.0021	KRAS/EGFR/JAK2/CTNNB1/SRC
GO:0098857	Membrane microdomain	0.0021	KRAS/EGFR/JAK2/CTNNB1/SRC
GO:0019897	Extrinsic component of plasma membrane	0.0021	KRAS/CTNNB1/APC/SRC
GO:0098589	Membrane region	0.0021	KRAS/EGFR/JAK2/CTNNB1/SRC
GO:0030877	Beta-catenin destruction complex	0.0021	CTNNB1/APC

Clinical actionability for the therapy of targeted agents

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).

Figure 6

(A) The clinical actionability of somatic mutations annotated according to OncoKB. The highest level of actionable alterations in (B) patients with advanced NSCLC (n=90), (C) the BM group (n=26), and (D) the non-BM group (n=64). Distribution of actionable alterations in (E) patients with advanced NSCLC, (F) the BM group, and (G) the non-BM group. Distribution of alteration types in (H) patients with advanced NSCLC, (I) the BM group, and (J) the non-BM group. NSCLC, non-small cell lung cancer; BM, brain metastases.

ctDNA analysis has a higher consistency of genomic profiling in the non-BM group as compared with that in the BM group

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).

Figure 7

Somatic mutation landscapes of patients with advanced NSCLC with (A and B) or without (C and D) BM from (A and C) tumor tissue DNA and (B and D) ctDNA. (E) Venn diagrams of mutated genes derived from tumor tissue DNA analysis and ctDNA analysis in the (top panel) BM group and (bottom panel) the non-BM group.

Discussion

Exploring the genomic alterations is crucial for clinical management in NSCLC patients with BM. Although dynamic mutation landscapes have been reported, systematic comparisons of genomic characteristics between the BM and the non-BM groups remain limited. In the present study, 174 somatic mutations of 35 mutated genes were identified in 90 patients with stage IV NSCLC using an NGS panel of 95 known cancer genes. Significant differences between the BM and the non-BM group were detected in somatic mutations, somatic interactions, key signaling pathways, functional biological information and clinical actionability for the therapy of targeted agents. Finally, it was also observed that ctDNA analysis presented with a higher consistency for genomic profiling of the non-BM than that of the BM group, indicating that ctDNA analysis may serve as a more credible alternative for genomic profiling in advanced NSCLC patients without BM.

In the present study, 17 mutated genes and 30 mutated genes were identified in the BM and the non-BM group, respectively. Among these genes, five genes, including ALK, CDA, SMAD4, SOD2 and VHL were uniquely present in patients with BM. ALK is a tyrosine kinase and its constitutively activated mutation renders ALK a formidable cancer driver gene (33,34). BM occurs frequently in tumors harboring ALK rearrangements (10), and its clinical significance has been considered to be critical. Several ALK inhibitors have been reported to demonstrate conspicuous activity in brain metastatic patients with crizotinib-resistant ALK-positive NSCLC, including second-generation (brigatinib and alectinib) (35-37), and third-generation therapeutics (lorlatinib) (38,39). Additionally, ALK and VHL (P=0.0435) are exclusively co-occurring genes in the BM group, which was reported in Chinese patients with NSCLC for the first time, to the best of our knowledge. More notably, as the first generation of the blood-brain barrier (BBB)-penetrating TKIs, AZD3759 can activate a p53-SMAD4 positive feedback loop and lead to apoptosis in hepatoma cells (40), offering a promising future approach for the treatment of brain metastatic NSCLC patients by AZD3759 (41). Collectively, the data of the present study may contribute to an improved comprehension of the underlying molecular mechanisms of patients with NSCLC with BM and may provide prospective therapeutic targets for this specific subgroup.

Additionally, five genes exclusively identified in the BM group were distributed in the Hippo signaling pathway, pyrimidine metabolism (PyM), and pantothenate and CoA biosynthesis, according to KEGG enrichment analyses. As a key mediator in the Hippo signaling pathway, Yes-associated protein (YAP) has been founded to facilitate drug resistance and tumorigenesis in NSCLC (42-45). Furthermore, YAP has been reported to play a crucial role for the promotion of BM in lung adenocarcinoma patients, and the inhibition of YAP by shRNA may significantly suppress migration and invasion abilities of metastatic NSCLC cell lines H2030-BrM3 in a murine model (46). Combined with the current results, these findings may provide prospective therapeutic approaches by modulating the members or mediators of the Hippo signaling pathway in brain metastatic patients with NSCLC. As a complex enzymatic network, the main function of PyM is to integrate de novo nucleotide synthesis, nucleoside salvage, and catalytic degradation of pyrimidines. In cancer cells, the de novo nucleotide synthesis pathway continuously provides deoxyribonucleoside triphosphates (dNTPs) to sustain uncontrolled proliferation, being different from normal cell de novo nucleotide synthesis (47,48). Until recently, PyM has been mainly implicated in the differentiation of leukemic cells; however, little is known about its roles in the differentiation of solid tumors (49). To the best of our knowledge, the present finding is the first report on PyM as an exclusive signaling pathway in Chinese patients with NSCLC with BM. However, further studies are required in order to discern whether PyM signaling pathway plays a critical role in the initiation and/or progression of BM in NSCLC.

More importantly, it was also demonstrated that ctDNA analysis is more feasible as an alternative for somatic mutation landscapes of non-BM patients than that of BM patients by the higher consistency between ctDNA analysis and tumor DNA analysis (43.33 vs. 29.41%). A possible explanation for the discrepancies between ctDNA analysis and tumor DNA analysis may be the inhibition of tumor cell release into the bloodstream by the BBB in patients with BM (50,51). Thus, the basic detection rate of genomic alterations derived from peripheral blood ctDNA in the BM group has been reported to be lower than that in the non-BM group (52,53). In a similar study, Aldea et al (53) demonstrated that ctDNA was positive in 52% of isolated central nervous system progression (iCNS) vs. 84% in extra-CNS only (noCNS), which was accompanied by a lower detection rate of pathogenic driver mutations (37 vs. 77%) and resistance alterations (6 vs. 45%). However, it cannot be overlooked that liquid biopsy is a potent method with which to improve the identification of actionable biomarkers when tumor tissue is unavailable (52,54,55). In the present study, plasma-derived ctDNA analysis improved the detection rate of EGFR actionable mutations by a 15.39% (4/26) increase in the BM group, and four patients had the opportunity to receive targeted therapies (erlotinib/erlotinib + ramucirumab/afatinib/gefitinib/osimertinib/dacomitinib) and/or participate in clinical trials. Consequently, the data of the present study are consistent with those of previous studies in which the identification of actionable mutations is growing in advanced NSCLC patients with the aid of plasma-derived ctDNA (56,57).

In the present study, one of the main limitations is that NGS data were obtained from tumor tissue DNA or plasma-derived ctDNA without the simultaneous analysis of matched normal tissue to delete the germline mutations. Thus, for a single sample the analysis was not complete; however, it may be considered adequate for the acquisition of actionable genomic alterations for the application of guided clinical treatment based on the suitable filter conditions (please see the 'Materials and methods' section, 'Bioinformatics analysis') (30,58-65). However, it cannot be disregarded that either multiregional biopsies or more than one type of biopsies may be costlier than the analysis of a single tumor sample without the inclusion of matched normal tissue. Another main limitation is the absence of a genomic profile, derived from brain tumor tissue DNA due to difficulties in the acquisition of brain tissue samples from NSCLC patients with BM. Additionally, although all patients were recruited for a prospective study, the collection of NGS data was retrospective. Lastly, further multiple-institution research with larger sample sizes is required to validate the present conclusions.

In conclusion, the somatic mutation landscapes of NSCLC with and without BM were compared, and significant differences in somatic mutations, somatic interactions, key signaling pathways, functional biological information, and clinical actionability for the therapy of targeted agents were observed between the BM group and the non-BM group. Moreover, plasma-derived ctDNA analysis may be a more reliable alternative for genomic profiling of advanced patients without BM, based on the higher consistency between ctDNA analysis and tumor DNA analysis in NSCLC.