Single-cell RNA profiling reveals an immunosuppressive microenvironment in EGFR double-mutant non-small cell lung cancer

Introduction

Lung cancer has the highest mortality rate among all cancer types, with non-small-cell lung cancer (NSCLC) accounting for approximately 85%. The EGFR and KRAS oncogenes are the two most common mutated genes in NSCLC (1). Targeted therapies using tyrosine kinase inhibitors (TKIs) have proven effective in NSCLC. Patients with progressive NSCLC often acquire additional gene mutations such as EGFR C797S after TKI failure (2,3). Previous studies have revealed that NSCLC patients with co-occurring genetic alterations experience notably shorter median progression-free survival (PFS) and overall survival than patients with single mutations (3-7). However, no standard treatment guidelines currently work for this co-driver gene mutant subgroup, which is difficult to study due to the low incidence of co-mutation and its poor prognosis. A deeper understanding of its mutational landscape and clinical implications is critical for guiding effective therapeutic strategies in NSCLC.

In 2015, Cai et al. found spatial intratumoral heterogeneity in NSCLC through bulk DNA detection of micro dissected lesions, and they observed that EGFR mutations did not always coexist with ALK fusion in different tissue areas (8). More recently, Chen et al. revealed genetic colony evolution through single-cell DNA sequencing in NSCLC, in which co-mutated cells were present in a low proportion, and new sub-clones with mutated genes appeared after treatment failure (9). However, the functional implications and mechanistic differences driven by concomitant EGFR double-mutant (DM) and EGFR single-mutant (SM) status remain unclear.

EGFR-mutant tumors are known to exhibit reduced CD8⁺ T cell infiltration and impaired cytotoxic function, resulting in a more suppressive immune microenvironment than that in EGFR wild-type (10). A single-cell RNA sequencing (scRNA-seq) study further confirmed that EGFR-mutant NSCLC is deficient in CD8⁺ T cells and that EGFR-mutant cancer cells have less contact with T cells than EGFR wild-type (11). It is reported that NSCLC harboring KRAS or TP53 mutations often have an immunosuppressive microenvironment (12). However, there is a notable lack of research investigating the tumor microenvironment (TME) in EGFR compound mutations at the single-cell level.

Conventional scRNA-seq techniques are unable to reliably identify and isolate specific mutated cancer cells, resulting in a limited understanding of their unique transcriptomic features. To address this limitation, we employed scFocuSCOPE, a scRNA-seq platform capable of detecting targeted driver gene mutations within individual cells. This method enables the simultaneous identification of mutations in targeted genomic regions and the comprehensive capture of whole-transcriptome expression profiles. By applying scFocuSCOPE to an EGFR-mutant NSCLC cohort, we were able to uncover distinct TME differences, offering new insights into the clinical implications of co-occurred driver mutations. We present this article in accordance with the MDAR reporting checklist (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2025-708/rc).

Methods

Study participants

Patients were diagnosed with NSCLC at the Guangdong Lung Cancer Institute, China, and fresh tissue specimens were collected through puncture biopsy or surgery from primary or metastatic tumor sites between March 2020 and September 2021. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Research Ethics Committee of Guangdong Provincial People’s Hospital [No. GDREC2019217H (R1)] and informed consent was obtained from all individual participants.

Definition of genetic alteration groups

All samples underwent DNA sequencing and were confirmed to harbor one or more driver gene alterations (Table S1). Samples were categorized into DM or SM groups primarily based on genetic status detected by targeted scRNA-seq profiling whenever available (Table S2). For a few samples without targeted scRNA-seq data (due to insufficient tissues or technical limitations), classification was determined as SM group.

The DM cohort in this study encompasses two distinct genetic scenarios. EGFR compound mutations: two or more oncogenic mutations within the EGFR gene (e.g., L858R + T790M). EGFR with concomitant driver alterations: an EGFR driver alteration co-occurring with one or more mutations in other driver genes (Table S1).

DNA extraction and next-generation sequencing

DNA was extracted from fresh-frozen tumor tissues and white blood cells, checked for quality, fragmented, and used for library construction with a 31-gene panel (Table S3, Berry Oncology Co., Ltd., Fuzhou, China). Libraries were sequenced on an Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) (13), and the resulting reads were aligned to the GRCh38 (hg38) human reference genome. Following targeted therapy for solid tumors, only patients with targeted genetic mutations were enrolled for scRNA-seq (Table S1).

Single-cell dissociation and two complementary DNA (cDNA) library constructions

Fresh tissues were preserved, washed, dissociated, and processed for cell viability analyses (14). Each cell was tagged with a barcode, unique molecular identifiers (UMIs), and probes targeting specific messenger RNA (mRNA) sequences. The captured mRNAs were reverse transcribed and amplified into two cDNA libraries, one for scRNA-seq library construction and the other for target scRNA-seq, following the protocol of the FocuSCOPE^® Single Cell Multiomics Lung Cancer Druggable Mutation Analysis Kit (Singleron, Cologne, Germany).

scRNA-seq raw reads analysis and quality control

The raw reads were processed to generate gene expression profiles using CeleScope^®. Barcodes and UMIs were extracted from the R1 reads and corrected. Adapter sequences and poly-A tails were trimmed from R2 reads using Cutadapt v3.7. After preprocessing, the R2 reads were aligned against the hg38 transcriptome using STAR (15) v2.6.1b. Uniquely mapped reads were assigned to exons using FeatureCounts (v2.0.1). Viral-track was used to filter background contamination (16). For each variant, Otsu’s method (17) was used to determine the supporting read threshold T.

scRNA-seq clustering and cell type annotation

Cells were filtered based on genes counted as <200 or in the top 2%, which were subsequently removed, excluding the top 2% UMIs. Cells with more than 20% mitochondrial content and genes expressed in fewer than five cells were also excluded. After filtering, 130,684 cells were retained for downstream analysis, with an average of 1,778 genes and 6,123 UMIs per cell.

Cell type identification, dimension reduction, and clustering analysis were performed using the Seurat program (18). After normalization and scaling the data, the top 2,000 variable genes were selected using Find Variable Features for principal component analysis. Cells were divided into 46 clusters by Find Clusters using the top 20 principal components and a resolution parameter set at 1.2. The uniform manifold approximation and projection (UMAP) algorithm was applied to visualize the cells in a two-dimensional space. Differentially expressed genes (DEGs) of the clusters were identified using the Find Markers function. The major cell types were subjected to a second round of unsupervised clustering. The clustering procedure started with an unfiltered expression matrix. Subsequently, doublets and low-quality cells were removed, and the resolution was set to 0.8 (18). The cell type identity of each cluster was determined based on the expression of canonical markers found in the DEGs using the SynEcoSys database.

Single-cell copy number analysis

Copy number variations in epithelial cells were assessed using the infer CNV R package (v1.8.0) (19). Genes expressed in more than three cells were sorted based on their loci on each chromosome. The relative expression values were centered at 1, with a 1.5 standard deviation from the residual-normalized expression values used as the ceiling. A slide window size of 101 genes was used to smooth the relative expression on each chromosome, removing the effect of gene-specific expression. A random subset of fibroblasts, dendritic cells (DCs), and granulocyte-monocyte progenitors was used as a reference, and the unsupervised dimensional reduction and sub-clustering of epithelial cells were measured.

DEGs, pathway enrichment, and trajectory analysis

Genes expressed in more than 10% of the cells in a cluster with a fold change >0.25 were selected as DEGs using Seurat v3.1.2 Find Markers based on the Wilcoxon rank-sum test with the default parameters. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted using the cluster Profiler R package v4.0.2. Pathways with an adjusted p (p_adj) value less than 0.05 were considered significantly enriched. The enrichment score was defined as log₁₀ (corrected P value). The pathway was upregulated in the DM group (positive score) or in the SM group (negative score). Cellular trajectory reconstruction analysis using gene counts and expression (Cyto TRACE) (20) was used to predict differentiation states from the scRNA-seq data. Monocle 2 (21) was used to construct the trajectory, and 2,000 highly variable genes were selected for dimensionality reduction using Seurat v3.1.2 Find Variable Features. Dimension reduction was performed using the DDR Tree.

Cell-cell interaction analysis

Cell-cell interactions between cell types were predicted based on known ligand receptor pairs using Cellphone DB (v2.1.7). The permutation number for calculating the null distribution of the average ligand-receptor pair expression in randomized cell identities was set to 1,000. Predicted interaction pairs with a P<0.05 and an average log expression >0.1 were considered significant.

Cancer cell signature analysis

The DEG profile was evaluated using the gene counting state (GCS). The correlation coefficient with GCS in Cyto TRACE was set at >0.6, resulting in 281 genes being filtered out, including H2AFZ, KRT18, SOD1, and KPM as “stem cell” genes; EIF1, EIF4A1, and EIF3F in mRNA translations, which were crucial for tumor cell plasticity; and Bcl11b and Atoh8 in cellular plasticity (Table S4). We then used the hallmarks of cancer cells as another signature by incorporating gene sets from the Molecular Signatures Database (MSigDB) [GSEA | MSigDB | Browse Gene Sets (gsea-msigdb.org)]. The expression of each gene set was calculated using the UCell signature score R package UCell v1.1.0 (https://github.com/carmonalab/UCell).

Histological analyses

Formalin-fixed paraffin-embedded NSCLC sections (4 µm) were prepared for hematoxylin and eosin staining and multicolored antibody detection. The immunohistochemistry (mIHC) panel includes antibodies targeting CD8 (SP16; ZSGB-BIO), programmed death-1 (PD-1) (UMAB199; ZSGB-BIO), programmed death ligand-1 (PD-L1) (SP263; Ventana), FOXP3 (236A/E7; Abcam), CD68 (PG-M1; ZSGB-BIO), and PanCK (AE1/AE3; ZSGB-BIO). mIHC was performed using the Opal 7-Color Automation IHC Kit (Akoya Biosciences, Marlborough, MA, USA) on a Leica Bond RXm fully automated immunostainer. The slides were scanned using the PerkinElmer Vectra Polaris system, which generates a single unmixed whole slide scan [pixel resolution: 0.50 µm (20×)] of up to seven colors. Stained slides were digitized using a multispectral slide imaging platform (Phenochart, Akoya Biosciences, Marlborough, MA, USA). After tissue segmentation and digital cell phenotyping, the density and positivity rate of each biomarker were semi-automatically assessed using the inForm Advanced Image Analysis software (inForm 2.5.1; Akoya Biosciences).

Plasmid construction

pLVX-EGFR L858R was constructed using nested polymerase chain reaction (PCR), and the EGFR L858R mutation was introduced using the following primers:

• P1F: TAGAGCTAGCGAATTCGCCACCATGCGACCCTCCGGG;
• P1R: AGTTTGGCCCGCCCAAAATCTGTGATC;
• P2F: ATTTTGGGCGGGCCAAACTGCTGGGTG;
• P2R: GATGATGATGCTCGAGTCATGCTCCAATAAATTCACTG.

The PCR products were cloned into the pLVX-puro backbone using the EcoR-I and Xho-I enzymes. The plasmid pLVX-ERBB2_S310F was constructed using the following primers:

• P1F: TCTAGAGCTAGCGAATTCGCCACCATGGAGCTGGCGG;
• P1R: CAGACGAGGGTGCAGAATCCCACGTCCGTAG;
• P2F: CTACGGACGTGGGATTCTGCACCCTCGTCTG;
• P2R: ATGATGATGATGCTCGAGTCACACTGGCACGTCCAGAC.

PCDH-EF1-Hygro_ERBB2 S310F, for the construction of DM cell lines, was obtained by cloning ERBB2 S310F from PLVX-ERBB2 S310F into the PCDH-EF1-Hygro backbone using the Nhel and NotI restriction sites.

Stable cell line establishment

H1299 cells were used to construct SM cell lines containing EGFR_L858R or ERBB2_S310F via lentiviral infection. Lentiviruses were produced using the pLVX-EGFR L858R and pLVX-ERBB2 S310F vectors with the packaging vectors pSPAX2 and pMD2.G in LentiX-293T (Takara, San Jose, CA, USA) cells, as previously described (PMID: 38412092). Transfection was performed using the PEI-Max reagent (Polysciences, Warrington, PA, USA), according to the manufacturer’s instructions. Four hours after transfection, the medium was replaced with fresh Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (FBS). The viral supernatant was collected and used to infect HT1299 cells. The stable integrants were selected using 1 µg/mL puromycin (Sigma-Aldrich Chemical Co., St. Louis, MO, USA) for two weeks and maintained in a complete DMEM medium at 0.25 µg/mL.

HT1299 cells with EGFR_L858R or ERBB2_S310F DMs were constructed based on HT1299 EGFR_L858R cells infected with virus packaging from the PCDH-ERBB2 S310F vector using the method described above. Cells were screened by culturing in a complete medium containing 200 µg/mL hygromycin (Sigma) for two weeks and maintained in a culture medium with 50 µg/mL hygromycin. All stable cells were cultured in a drug-free medium for 48 h before use.

Cell counting kit 8 (CCK8) cytotoxicity assay

Approximately 4,000 cells per well were seeded in 96-well plates for 24 h. The cells were then treated with a series of drug concentrations diluted in the corresponding culture medium. Each measurement was repeated at least three times. After incubation with the indicated drugs for 72 h, 10 µL of the CCK-8 reagent (Med Chem Express, Monmouth Junction, NJ, USA) was added to each well and incubated for a further 2 h. The optical density (OD) at 450 nm was measured using a multifunctional microplate reader (Infinite M200 Pro, Tecan, Männedorf, Switzerland). The percentage of each concentration that accounted for the control was presented as cell viability. The IC50values were calculated using SPSS.

Colony formation

HT1299, HT1299 EGFR_L858R, HT1299 ERBB2_S310F, and HT1299 EGFR_L858R/ERBB2_S310F cells were seeded in six-well plates at 1×10³ cells/well, and on the next day, they were treated with the indicated drugs. The drug concentrations used were as follows: 0.6 µL for dimethyl sulfoxide (DMSO), 0.1 µM for pyrotinib, 0.4 µg/mL for DS-8201a, 3 µM for efitinib, 0.1 µM for afatinib, and 1 µM for osimertinib. After 14 days, the colonies were fixed, stained with crystal violet, and counted. The area of colonies per well was measured using ImageJ software.

Western blot

To compare the protein expression in HT1299, HT1299 EGFR_L858R, HT1299 ERBB2_S310F, HT1299 EGFR_L858R, and ERBB2_S310F cells, cell lysates were prepared in NP-40 buffer [20 mM Tris pH 7.4, 150 mM NaCl, 5 mM EDTA, 0.5% NP-40, 1 mM DTT, supplemented with EDTA-free protease inhibitor tablet (Roche, Basel, Switzerland)], and immunoblotting was performed as described previously (PMID: 38412092). Anti-EGFR, ERBB2, ERK1/2, and AKT antibodies were purchased from Abcam (Cambridge, UK), and phospho-specific antibodies were purchased from BD Biosciences (San Diego, CA, USA). Pan-PARP antibodies were obtained from Cell Signaling Technology (Beverly, MA, USA). The actin antibody was purchased from Santa Cruz Biotechnology (Santa Cruz, CA, USA).

RNA sequencing

Total RNA was assessed for purity and concentration using Nano Drop 2000, and RNA integrity and quantity were evaluated using an Agilent 2100/4200 system. mRNA was isolated from the total RNA, fragmented, and converted into cDNA. Following blunt-end conversion and adaptor ligation, the cDNA library was enriched by PCR. The concentration of the final library was checked using a Qubit^® fluorometer and quantitative PCR (qPCR), and size distribution was confirmed by agarose gel electrophoresis. The prepared libraries were pooled and sequenced on an Illumina platform using PE150for lncRNA-Seq. Clean reads were aligned to the GRCh38 (hg38) human reference genome using Hisat2, and gene expression levels were quantified using Feature Count. DEGs were identified using Edge R, with significant genes defined as those with |log₂ (fold change)|>1 and Q value <0.05, after multiple testing corrections.

Statistical analysis

Statistical analyses were performed using R (v4.0.3) software. Gene expression or signatures between two groups of cells were compared using an unpaired two-tailed Student’s t-test. Cell distributions of paired samples were compared using a paired two-tailed Wilcoxon rank-sum test. For group comparisons, an unpaired two-sided Wilcoxon rank-sum test was used, unless otherwise specified. Statistical significance was set at P<0.05.

Results

Single-cell capturing and sample grouping

To overcome the limitations of conventional scRNA-seq in capturing the small population of mutated cancer cells, we employed scFocuSCOPE, a targeted scRNA-seq technique designed to detect targeted gene mutations at the single-cell level. The detailed workflow is presented in the Methods section and illustrated in Figure 1A and Figure S1A. Two cDNA libraries were constructed in the process: a common scRNA-seq library for the TME study and a specific driver gene library to identify the genetic mutation status in every cancer cell.

Figure 1 Single-cell capture and sample grouping. (A) Pipelines of the conventional single-cell and targeted mutated single-cell RNA sequencing technology (scFocuSCOPE). These two cDNA libraries were obtained and combined for the subsequent data analysis. (B) UMAP of 130,684 color-coded cells from distinct cell types, including cancer cells, lymphocytes, MPs, ECs, pDCs, AT1, and AT2. (C) UMAP of all cells from two groups: DM and SM driver genes. (D) Patients (n=7) with double gene mutations at both DNA and single-cell RNA levels. DM-P1 and DM-P4 showed only the first sample. Driver gene mutation status at the DNA level (left) and the Venn diagram listing the mutated cancer cell fractions (right) are presented. (E) Boxplots showing the fraction of each cell type in the DM and SM groups, analyzed using an unpaired two-sided Wilcoxon rank-sum test. *, P<0.05. AT1, alveolar type I; AT2, alveolar type II; DM, double-mutant; ECs, endothelial cells; MPs, mononuclear phagocytes; pDCs, plasmacytoid dendritic cells; SM, single-mutant; UMAP, uniform manifold approximation and projection.

Although the low incidence of co-mutation, 25 NSCLC samples were finally collected and all were identified with EGFR alterations (Table S2 by DNA sequencing (Berry Oncology, China) (Figure S1B). Although 17 samples were detected with two or more driver gene mutations based on DNA sequencing, only nine samples from seven patients had more than two types of genetic alterations at the scRNA level. In total, 130,684 qualified cells from these samples were annotated into 12 clusters based on their genetic signatures (Figure 1B,1C, and Figure S1C). To better understand the characteristics of cancer cells at the single-cell level, we classified the samples primarily based on their genetic status by the targeted scRNA-seq results (Table S2, groups determined in the “Methods” section). Nine samples from seven patients with NSCLC harboring EGFR compound or concomitant driver gene variations were categorized as the DM group (Figure 1D). Sixteen samples from 15 patients were categorized as the SM group. Four samples from four patients with DM were treatment-naïve, whereas the others were treatment-resistant. For the genetic evolutionary study, we combined the data from treatment-naïve and treatment-resistant samples (9) (Figure S1B).

Epithelial cells were classified into cancer and normal cells (alveolar type I and alveolar type II, AT1/2; Figure 1B). Most cancer cells in the DM samples were non-mutated, and double-mutated cells were a minority among all mutated cells (Figure 1D), consistent with previous research on intratumoral heterogeneity (8). No cells with co-mutations were identified in the DM sample, DM-P7-Lu (Figure 1D). Additionally, significantly more lymphocytes and fewer mononuclear phagocytes (MPs) were observed in the SM group than in the DM group (Figure 1E).

Molecular characterization of cancer hallmarks in SM and DM tumor cells

All cancer cells from both DM and SM groups were subdivided into 14 clusters (Figure 2A), revealing remarkable discrepancies and heterogeneity among the tumor cells (Figure 2B and Figure S2A). To elucidate the characteristics of SM and DM cancer cells, we analyzed U Cell scores of cancer hallmarks (22-24) (Figure 2C,2D). These cancer cell hallmarks include sustained angiogenesis (ANGIOGENESIS); programmed cell death (APOPTOSIS); genome instability and mutation (DNA_REPAIR); tissue invasion and metastasis (EMT); self-sufficiency in growth signals (PROLIFERATION), reprogramming of energy metabolism (GLYCOLYSIS), tumor-promoting inflammation (INFLAMMATION), and plasticity (22,25-29) (Table S5).

Figure 2 Analysis of oncogenic signatures in single- and double-mutant cancer cells. (A) UMAP of 14 color-coded cancer cell clusters labeled as c1 to c14. (B) UMAP visualization of cancer cells from DM (burgundy) and SM (light denim) patients. (C) Octagons represented the statistical scores of tumor cell hallmarks for DM and SM groups, including ANGIOGENESIS, APOPTOSIS, DNA_REPAIR, EMT, PROLIFERATION, GLYCOLYSIS, INFLAMMATION, and PLASTICITY as mentioned in the “Methods” section. (D) Heatmap displaying each cancer cell with the U cell scores for cancer cell hallmarks. (E) U cell scores of cancer cell hallmark-relevant genes were compared between SM and DM groups. ***, P<0.001. (F) Trajectories of cancer cells in the DM and SM groups; the direction of the arrow represented the estimated direction of evolution. (G) KEGG analysis showing the differential genes in cancer cells between two groups. (H) PD-L1 expression (red dots) levels at each cancer cell based on (B). (I) Communications among all distinct cell types were analyzed and represented as the minus of potential interaction numbers between the DM and SM groups. All differences were significant. DM, double-mutant; KEGG, Kyoto Encyclopedia of Genes and Genomes; SM, single-mutant; UMAP, uniform manifold approximation and projection.

Two groups exhibited significantly different cancer hallmarks (Figure 2E, Figure S2B,S2C). Most SM cancer cells exhibited higher enrichment of plasticity, proliferation, and DNA repair, whereas majority of DM cancer cells (cluster 8–10, 12, and 14) were characterized by angiogenesis, EMT, and inflammation (Figure 2D, Figure S2D). The DEGs set in cancer cells were determined using the U Cell score and compared with the Cyto TRACE score, which revealed similar differences between the DM and SM cohorts (Figure S2B,S2C,S2E). Consistent with their plasticity, SM cancer cells showed a more uniform distribution in cell velocity analysis (Figure 2F, Figure S2F). From the cell trajectory direction, more DM cancer cells clustered at the endpoints and were excluded from the starting point, whereas an almost equal number of SM cancer cells were scattered at every phase and preferred clustering at the beginning (Figure 2F). SM cancer cells were more probably in the early stages of evolution than DM cancer cells, aligning with the hallmark of plasticity observed in SM cancer cells. To corroborate the cancer hallmarks, KEGG analyses revealed significantly higher active pathways of EMT-related PI3K-ATK and focal adhesion in the DM cohort than SM cohort (24,30-32) (Figure 2G, Figure S2G). Conversely, the oxidative phosphorylation and viral carcinogenesis pathways associated with plasticity were enriched in the SM group compared to those in the DM group (Figure 2G). Alterations in driver genes have been reported to upregulate the expression of PD-L1 (33,34), which may contribute to immune escape. More DM cancer cells expressed PD-L1than SM cancer cells (Figure 2H, Figure S2H). Tumor cells in the DM group interacted more with adjacent normal cells, such as AT1, mural cells, and endothelial cells (ECs), and had less contact with fibroblasts than those in the SM group (Figure 2I, Figure S2I).

“Desert” immune microenvironment of DM cohort

The number of immune cells is reported to be significantly reduced in patients harboring EGFR SM (11). In this study, we attempted to determine the immune microenvironment in patients with DM. Significantly fewer lymphocytes were present in the DM group (Figure 1D). To determine the composition and aggressiveness of T cells, lymphocytes were reclustered into three major types: T cells, natural killer (NK) cells, and B cells (Figure 3A). Subsequently, they were subdivided into seven, five, and four sub-clusters, respectively (Figure 3A, Figure S3A). In total, 1,539 T cells were obtained from the DM (n=9) and 15,388 T cells from the SM (n=15) groups (Figure 3A), which limited tissue availability and fewer samples in DM group might have contributed to this large discrepancy in T cell counts. Fewer naïve T cells and more proliferating T cells were observed in the DM cohort (Figure 3B). T cells from the DM cohort were located almost exclusively in the terminal stage, whereas SM T cells were more evenly distributed across all trajectory directions (Figure 3C). Proliferating T cells in the SM group were concentrated at both endpoints (Figure 3C and Figure S3B). To further validate the killing ability of T cells, we computed the activation- and exhaustion-related genes of T cells; activation-related genes included GNLY, GZMB, NKG7, GZMA, GZMK, and PRF1, and exhaustion-related genes included CTLA4, LAG3, TIGIT, HAVCR2, PDCD1, and TOX (35,36) (Figure 3D). The expression levels of T cell activity-related genes were often lower in the DM group than those in the SM group, particularly in proliferating T cells. Conversely, almost all exhaustion-related genes showed significantly higher expression levels in the DM cohort than in the SM group. Overall, lower activation and more exhausted behaviors of T cells in the DM cohort probably suggested weakened T cell functions in the DM group for tumor cell killing. Focal adhesion, PI3K-ATK, and gap junction signaling pathway enrichment of T cells in the DM group indicated their high potential to promote tumor invasion and metastasis by the DEGs and KEGG analysis (Figure S3C,S3D), similar to the characteristics of DM cancer cells. Herpes simplex infection and leishmaniasis pathways were abundant in the SM group, suggesting the activated status of T cells (Figure S3D). The concept of the “desert” immune TME in the DM group was further confirmed by lower pathway enrichment for T cell receptor signaling, NK-mediated cytotoxicity, and antigen processing and presentation (Figure 3E).

Figure 3 Immune “desert” of the double-mutant cohort. (A) T cells, NK cells, and B cells were subdivided into seven, five, and four subtypes, respectively, from lymphocytes. We identified naïve T cells, exhausted T cells (CD8 Tex and CD4 Tex), cytotoxic T cells (CD8 Teff), regulatory T cells (Treg), proliferating T cells, and NKT cells from T cells. NK1 to NK5 were from NK cells. Memory B cells (B mem), naïve B cells (Naïve B), plasma B cells, and a cluster of highly expressing TPT1 cells (B-TPT1) were identified from B cells. DM group-derived lymphocytes were colored in burgundy, and SM cohort in light denim in the UMAP. (B) Boxplot showing the percentage of each T cell subcluster from the DM (burgundy) or SM (light denim) group compared with a paired two-tailed Wilcoxon rank-sum test. *, P<0.05. (C) Trajectories of each T cell subtype from the DM or SM group with distinct colors and the direction of the arrow represented the estimated direction of evolution. (D) Activation (left) and exhaustion (right) signature genes of every T cell subtype were analyzed and compared using an unpaired two-tailed Student’s t-test. The solid line indicated a significant difference between the DM and SM groups (P<0.05). (E) KEGG analysis presented the pathway enrichment of differential genes in proliferating T cells between the two groups. (F) Average expression levels of genes associated with proliferation in proliferating T cells. The dot size represented the percentage of cells expressing the targeted genes. Color implies expression levels. (G) Cell communications between cancer cells (DM, left; SM, right) and T cell subtypes. Square size was in accordance with the frequency of cell interaction. DM, double-mutant; KEGG, Kyoto Encyclopedia of Genes and Genomes; SM, single-mutant; UMAP, uniform manifold approximation and projection.

The characteristics of proliferating T cells were analyzed by cell cycle- or proliferation-associated genes, such as MKI67, CDK126, and TOP2A (37,38). The average expression levels of these genes indicated better regrowth ability to proliferate T cells in the DM cohort (Figure 3F), as confirmed by the DEGs and KEGG analysis (Figure 3E, Figure S3E). To some extent, the higher proportion and better functional proliferating T cells in the DM cohort indicated the exhaustion of other T cell subtypes.

Cancer cells exhibited probable ligand-receptor interaction states in both groups; however, the SM cohort exhibited stronger interactions with proliferating T cells and weaker interactions with regulatory T cells (Tregs) (Figure 3G, Figure S3F). Overall, fewer quantities and poor functional T cells in the DM group could contribute to a “desert” immune environment with poor outcomes.

Although the proportions of B cell sub-clusters between the SM and DM groups were not significantly different (Figure S4A,S4B), ERBB2 and B cell receptor signaling pathways were more active in the DM group than those in the SM group (Figure S4C,S4D). The NK cell cluster, a relatively small population, was subdivided into five clusters based on regularly associated marker genes (Figure S4E). Distinct immune response modes of NK cells were observed in the two cohorts, as verified by KEGG enrichment analysis (Figure S4F). The DM group showed a slightly higher proportion of NK1 and NK2 clusters (Figure S4G). NK cells in the DM group were in a more active physiological state because of the higher enrichment of the FoxO and MAPK signaling pathways, whereas Fc gamma R-mediated phagocytosis was less active (Figure S4F,S4H). Higher GZMH expression in the SM group suggested a stronger antitumor ability of NK cells than that in the DM cohort (Figure S4I).

Immune-suppressive macrophages in the DM cohort

In addition to T cells, we explored differences in other immune cells between the two groups, such as macrophages, mature DCs, type 2 conventional dendritic cells (cDC2), osteoclasts, master cells, and granulocyte-macrophage progenitor cells (Figure 4A). Macrophages were divided into nine subtypes (Figure 4A) based on distinct gene set annotations (Figure 4B and Figure S5A-S5C). The others were sub-clustered into mature DCs, cDC2, and plasmacytoid dendritic cells (pDCs). Cell trajectory analysis of macrophages was similar to that observed in cancer and T cells (Figure 4C). The direction of the cell trajectory indicated the exhaustion of macrophages in the DM group and the potential differentiation of SM cells into DM cells. Therefore, development-related pathways, focal adhesion, and the FoxO signaling pathway were more upregulated in the SM group than those in the DM group (Figure 4D). NEAT1 and NEAT2 (MALAT1), cell proliferation-related genes, were highly expressed in the Mac-NEAT1 subcluster (Figure 4B and Figure S5D). The SM group had a remarkably higher proportion of cells in the Mac-NEAT1 cluster than that in the DM group (Figure 4E). Other macrophage clusters exhibited similar proportions in both cohorts (Figure 4E). NEAT 1/2 involving pathways, such as regulation of the actin cytoskeleton and MAPK signaling, was more active in the SM group than in the DM group (Figure 4F).

Figure 4 Immune-suppressive macrophages in the double-mutant cohort. (A) Macrophages were sub-clustered into nine types, and DC into three types, as exhibited in the UMAP. DM group originated macrophages were colored burgundy, and those from the SM group were shown in light denim, with the same for DC. (B) Signature genes of each T cell subcluster showed distinct expression levels and cell proportions. (C) Cell trajectories of each macrophage subtype from the DM (left) or SM (right) group. Each dot represented one cell. (D) KEGG enrichment pathways in macrophages were listed between DM and SM groups. (E) The fraction of macrophage subtypes (%) was measured using an unpaired two-sided Wilcoxon rank-sum test. **, P<0.01. (F) KEGG enrichment in MAC-NEAT1 cells was shown between the two groups. DC, dendritic cells; DM, double-mutant; KEGG, Kyoto Encyclopedia of Genes and Genomes; SM, single-mutant; UMAP, uniform manifold approximation and projection.

Communication among the different macrophage subtypes and cancer cells was similar (Figure 4F and Figure S5E,S5F), indicating similar crosstalk between macrophages from the DM and SM groups and cancer cells, including Mac-NEAT1. The cell proportion of mature DCs, cDC2, and pDCs between the two cohorts had no significant differences (Figure S5G). Overall, the smaller proportion of macrophages and less active Mac-NEAT1 subcluster in the DM cohort than that in the SM group suggested a slightly more immunosuppressive environment.

Sequential TME detection revealed clone evolution

To further explore clonal evolution and the impact of treatment on the immune microenvironment, we conducted related analyses in patients from the DM group. Two patients from the DM group (DM-P1 and DM-P4) both provided two serial scRNA-seq samples (Figure 5A,5B). Prior to targeted therapy and DNA sequencing, patient DM-P1 was diagnosed with NSCLC with EGFR 19del and EGFR T790M (Table S2). Subsequently, the EGFR C797S mutation was acquired after osimertinib treatment (DM-P1-Li). Two time-sequential specimens (DM-P4-LN and DM-P4-SN) from DM-P4 had identical concurrent EGFR L858R and ERBB2 S310F mutations (Figure 1C,5A,5B). After TKI treatment, the proportion of cancer cells decreased in both patients (Figure S6A,S6B). Consistent with previous observations of more MPs and fewer lymphocytes in the DM group than those in the SM cohort, a high proportion of MPs and an extremely low percentage of lymphocytes were identified in the cell annotation regardless of treatment (Figure S6A,S6B). After TKI treatment, gene enrichment downstream of the EGFR and bypass signaling pathways (PI3K-Akt and Ras signaling pathways) was enhanced, suggesting that the invasion, metastasis, and inflammation of tumors were enhanced after TKI failure (Figure 5C,5D, and Figure S6C,S6D). As aforementioned, among the eight characteristics of cancer cells, the hallmarks of apoptosis, DNA damage repair, EMT, inflammation, and plasticity were significantly enhanced in both patients after targeted therapy (Figure 5E,5F).

Figure 5 Changes of the immune microenvironment before and after treatment in two patients with double mutation. (A) Timeline of patient DM-P1 (up). Tissue sampling for scRNA-seq was indicated before treatment with tyrosine kinase inhibitors (Prior TKIs, DM-P1-Lu) and post-TKI treatment (Post TKIs, DM-P1-Li). Fishbone displayed the gene mutational dynamics at the single-cell RNA level. Mutations were shown in different colors. (B) Timeline of patient DM-P4 (up). Two tissue samples were obtained for scRNA-seq, before TKI treatment (Prior, DM-P4-LN) and post-TKI treatment (Post TKIs, DM-P4-SN). Three samples were obtained for mIHC, prior TKIs, post-TKIs, and Post ICIs. (C,D) Differential gene KEGG enrichment of the two patients between Prior TKIs and Post TKIs was analyzed. (E,F) Dynamics of cancer cell hallmark tendency of DM-P1 (E) and DM-P4 (F) mentioned during TKI treatment as described above were analyzed using an unpaired two-sided Wilcoxon rank-sum test. Except for the hallmark of glycolysis in DM-P4, comparisons of other cancer cell hallmarks showed significant differences. (G) Three samples (Prior TKIs, Post TKIs, and Post ICIs) from DM-P4 underwent mIHC analysis. The vertical Y-axis represented the proportion of positive cells (%) per field. PanCK (malignant cells, red), CD8 (CD8⁺ T cell, cyan), CD68 (macrophage, orange), FOXP3 (regulatory T cell, yellow), PD-1 (magenta), PD-L1 (green), and DAPI (blue) were shown with single or total colors. The cell numbers per area with specific fluorescence between three samples were compared. *, P≤0.05. Scale bar 200 µm. Anl, anlotinib; AP, pemetrexed + cisplatin; Atezo, atezolizumab; Bev, bevacizumab; Car, carboplatin; DAPI, 4',6-diamidino-2-phenylindole dihydrochloride; Erl, erlotinib; Gef, gefitinib; IMRT, intensity-modulated radiation therapy; KEGG, Kyoto Encyclopedia of Genes and Genomes; Osi, osimertinib; Pem, pemetrexed; Post ICIs, post-immune checkpoint inhibitors; Pry, pyrotinib; TKI, tyrosine kinase inhibitor.

To further explore the dynamics of the immune TME during treatment, samples from one patient, DM-P4 prior to TKI, post-TKIs, and post-immune checkpoint inhibitor (ICI) treatment were analyzed using multiplex fluorescence mIHC (Figure 5G). The number of CD8⁺ T cells decreased significantly after targeted therapy, while increasing after immunotherapy (n=1) (Figure 5G). Conversely, the number of macrophages increased significantly after targeted therapy, but reduced after immunotherapy. The number of PD-L1-expressing cells decreased remarkably after immunotherapy. However, no significant difference was observed before and after targeted therapy (Figure 5G, Figure S6E). This suggested an ongoing immune escape after immunotherapy.

Differentiation identification between the DM and SM groups in vitro

As demonstrated in the DM-P4 patient harboring EGFR L858R and ERBB2 S310F mutations, dynamic changes in the TME were observed during treatments (Figure 5B,5G). To validate these single-cell level findings, we conducted in vitro experiments using engineered cell lines. Specifically, EGFR L858R and ERBB2 S310F mutations were introduced into wild-type H1299 lung cancer cells to generate corresponding mutant cell lines for further functional analysis (Figure S7A,S7B). Consistent with the elevated PD-L1 levels in the DM (Figure 2H), the expression of PD-L1 was significantly higher in the DM cell line than that in the SM line (Figure 6A,6B). Bulk RNA sequencing of these cells revealed significantly higher levels of cell growth-related genes (IGF2) and matrix-interacting proteins (SPP1) in DM cells than in SM cells (Figure 6B-6D).In the ERBB2 mutant line, an increase in phosphorylated EGFR (pEGFR) protein levels was observed (Figure 6B), probably because of the formation of EGFR-ERBB2 heterodimers, as suggested by previous studies (39). The downstream RTK pathway proteins, pAKT and pERK1/2, were more elevated in the DM cells than those in SM cells (Figure 6B), consistent with prior scRNA-seq analysis in cancer cells (Figure S2G). Pathways involved in viral process regulation and vasculature development were more highly activated in the EGFR-mutant cells than those in the DM cells (Figure 6E, Figure S7C-S7L), with viral carcinogenesis activation similar to that observed in patients with SM, as indicated by scRNA-seq analysis (Figure 2G). Almost no difference was observed in the effects of EGFR TKIs between EGFR SM and DM cell lines. However, afatinib was more effective than gefitinib and osimertinib in killing DM cells (Figure 6F,6G). Upon treatment with TKIs, the DM cell line exhibited slightly higher colony formation than that of the SM EGFR or ERBB2 cell lines, suggesting a marginally worse response to the drugs (Figure 6H,6I). Notably, a few colonies formed in the ERBB2 mutant cell line treated with the antibody-drug conjugate, DS-8201a (trastuzumab deruxtecan), highlighting the potent cytotoxicity of DS-8201a. Treatment with either EGFR or ERBB2 TKIs alone inhibited colony formation in the DM cell only mildly. Interestingly, significant reductions in colony numbers were achieved when EGFR and ERBB2 TKIs were combined. This indicated that DM cells are probably more sensitive to dual-targeted drugs than to single-targeted drugs (Figure 6H). However, EGFR and the downstream RTK pathway proteins had no significant differences between single TKI and combined TKIs (Figure S8).

Figure 6 Differences between single and double mutations in the H1299 lung cancer cell lines. (A) Western blot results showed EGFR, ERBB2, downstream signaling proteins ERK1/2, AKT, and corresponding activated proteins in the four cell lines. (B) Heatmap of differential gene expression levels from RNA sequencing results for the H1299 wild-type, ERBB2m, EGFRm, and ERBB2/EGFR DM cell lines. (C,D) Volcano plots showed DEGs in the DM cell line compared to those in the ERBB2m or EGFRm cell lines, respectively. (E) GSEA and GO enrichment analysis of DEGs in the DM cell line versus the EGFR SM cell line. (F,G) IC50 curves for targeted drugs against ERBB2 (F) and EGFR. (E) In different mutant cell lines: DM (burgundy), ERBB2m (lake-green), and EGFRm (sky-blue). Except for the concentration of DS-8201a in µg/mL, all others were in µM. (H) Colony formation assay was performed for four cell lines after targeted drug treatment; DMSO served as a control. (I) The bar graph showed the statistical significance of the colony formation assay results. *, P≤0.05. A, afatinib; D, DS-8201a; DEGs, differentially expressed genes; DM, double-mutant; EGFRm, epidermal growth factor receptor single-mutant; ERBB2m, erbb2 receptor tyrosine kinase single-mutant; G, gefitinib; GO, gene ontology; GSEA, gene set enrichment analysis; O, osimertinib; P, pyrotinib.

Discussion

In this study, we used scRNA sequencing to explore the microenvironment and heterogeneity in double and single EGFR-mutated NSCLC. While substantial heterogeneity exists, most cancer cellular clusters in DM NSCLC exhibited more invasive tumor characteristics. Furthermore, the overall microenvironment in DM tumors was more immunosuppressive than in SM.

Previously, mutated cancer cells in solid tumors were difficult to capture and analyze at the single-cell level because of the low proportion of mutated cancer cells and technical limitations. The key technological advancement in this study, scFocuSCOPE, not only provides comprehensive whole-transcriptome profiling but also significantly enhances the accuracy and efficiency of identifying mutation-bearing cancer cells. Interestingly, a subset of cases showed inconsistent mutation profiles between targeted scRNA-seq and DNA sequencing, which may be attributed to a range of biological and technical factors. To our knowledge, this is the first study to explore the origin and evolution of gene mutations in cancer cells using single-cell transcriptomics.

Patients with double-driver gene mutations or EGFR compound mutations have been reported to have poorer outcomes, including reduced treatment responses and shorter survival (4,7,40-43). Several studies have examined differences between EGFR DM and SM (8,9,44). Zhao et al. demonstrated that patients with NSCLC harboring a single EGFR 19Del had longer PFS than that in patients with two or more EGFR mutations, particularly in cases with a highly active p53 signaling pathway (45). While the TME of DM NSCLC remains largely unexplored. Cai et al. proposed that concurrent gene mutations co-exist within some cancer cells (8). Chen et al. also found DNA-level evidence in osimertinib-treated patients with NSCLC (9), showing only a few cells showing co-occurring mutation (8,9). Our study aimed to further elucidate clonal genetic origin and evolution at the single-cell level. The transcriptional landscape analysis presented here is the first demonstration that EGFR DM coexists within a single cancer cell in treatment-naïve patients, confirming the oncogenic hypothesis of EGFR DM. In all DM samples, the majority of cancer cells were non-mutated, and double-mutated cells constituted a minority among all mutated cells, highlighting the large heterogeneity of cancer cells in the DM group. This suggests that DM cancer cells may arise and develop earlier and more separately than commonly believed, which may affect treatment response and clinical outcomes.

To explore the mechanism underlying poor outcomes in patients with DM, we examined the distinct characteristics of cancer cells in the DM and SM groups. Cancer cells in the SM group showed more potent plasticity and higher expression of plasticity-related genes than those in the DM group. Cell trajectory analysis suggested that SM cells may evolve and differentiate into DM cells through sub-clonal progression. Malignant epithelial cells in the DM group exhibited higher angiogenesis, tissue invasion and metastasis, which may explain the poor response to single TKI treatment. However, dual-targeted treatment might overcome this difficulty in patients with DM NSCLC, as the response to combined targeted drugs is significantly better than that to single TKI drugs in DM cell lines. However, it is important to note that our in vitro drug experiments were conducted only in cell lines harboring the EGFR L858R and ERBB2 S310F double mutations. Therefore, the therapeutic responses observed may be narrowly applicable to this specific genetic context and may not be generalizable to other types of double mutations. Ongoing clinical observations in lung cancer have demonstrated the efficacy of dual TKIs (46-50). Further preclinical and clinical studies are warranted to explore drug efficacy in other DM subtypes.

It is well established that EGFR-mutant NSCLCs are generally considered “cold tumors” with poor responses to immunotherapy. Previous scRNA-seq studies showed distinct evolutionary trajectories in EGFR-mutant NSCLC compared to EGFR wild-type and found significantly fewer immune T cells, indicating a poor immune environment and weak immunotherapy response (11). Several clinical trials have demonstrated that combining immunotherapy with targeted therapy leads to worse clinical outcomes compared to targeted therapy alone (51-53). While EGFR SM is known to shape a poor immune environment (11), our findings suggest that DM cancer cells have an even stronger effect in creating an immune “desert”. A limited number of T cells in the DM group appeared to be more exhausted and less activated than those in the SM samples, particularly fewer CD8⁺ Teff cells (Figure 1D). This desert microenvironment may contribute to the poor outcome of patients with DM. It is important to note that this observation is based on a single patient. Therefore, it requires further validation in larger cohorts.

A previous study using cell lines demonstrated that high levels of EGFR activation and phosphorylation significantly upregulated PD-L1 expression (34). Targeted drugs that reduced EGFR phosphorylation levels consistently downregulated PD-L1 expression in vitro (34). However, other studies have reported that in SM patients, targeted therapy does not significantly affect the proportion of PD-L1-expressing cells (54,55). In our study, the DM group exhibited significantly higher PD-L1 expression compared to the SM cohort. Interestingly, we also found that immunotherapy further reduced PD-L1 expression in one DM patient (Figure 5G), providing a potential mechanistic explanation for the limited efficacy of ICIs in EGFR-mutant NSCLC. Nonetheless, this observation should be interpreted with caution, as it is based on mIHC results from only one DM patient, and further validation in larger, independent cohorts is necessary to confirm these findings.

In our study, higher EGFR mRNA expression and more activated EGFR-related downstream pathways were observed in the DM cell line than in the SM cell line. Correspondingly, slightly higher PD-L1 expression was observed in the DM cells. Co-culture experiments with peripheral blood mononuclear cells (PBMCs) in a previous study revealed that cell lines with high PD-L1 expression elevated T cell apoptosis levels in PBMCs (34). Based on these observations, we hypothesize that the aggressive behavior of DM cancer cells, combined with the suppressive immune TME, contributes to poor clinical outcomes. However, whether DM cancer cells have a similar effect and the mechanism of immune regulation remains unclear. Further studies are required to address these questions.

There are a few limitations in this study. First, scFocuSCOPE sequencing could only detect the targeted driver gene mutations in specific regions. Second, all samples were processed using scRNA-seq; however, a few patients with SM did not have sufficient samples for scFocuSCOPE sequencing. Third, owing to the low incidence of double mutations, the recruited DM cohort was smaller than the SM cohort. The sample size of the DM cohort (n=7) limits the statistical power of certain comparisons. While we observed consistent and pronounced trends in the TME between DM and SM groups, future studies with larger cohorts are needed to confirm these findings with greater statistical confidence. Fourth, the patients did not have identical stages or clinical treatment regimens, and comparing the differences in clinical outcomes between the two groups was difficult. Fifth, we acknowledge the limitations of the limited longitudinal cohort (n=2). Subsequent studies in genetically matched patient groups under uniform treatment would help establish broader validity. Sixth, our DM cohort contained heterogeneous mutation patterns, each with potentially distinct biological properties. Larger studies with mutation-stratified cohorts are needed to determine whether these findings extend to all EGFR DM populations. At last, we adopted only one H1299 cell line, which may not fully recapitulate the biology of patient-derived models or primary tumors. Future validation in patient-derived organoids or xenograft models and different cell lines with varied compound driver gene mutations would strengthen these observations.

In summary, we employed conventional and targeted single-cell transcriptional profiles, revealing the substantial heterogeneity among malignant cancer cells, and the most cellular subpopulations displayed a more suppressive immune microenvironment in patients with DM NSCLC. These findings not only enhance our understanding of the molecular and cellular complexities underlying DM NSCLC but also have broad implications for the field. Further detailed studies with large cohorts are needed to validate the conclusions and related mechanisms.

Conclusions

Using targeted driver gene-capturing scRNA-seq, we identified a higher metastatic potential in DM non-small cell lung cancer patients compared to SM patients, characterized by an “immune desert” of fewer, dysfunctional immune cells.

Acknowledgments

All authors appreciate the helps from the staff in Guangdong Lung Cancer Institute, Guangdong Provincial People’s Hospital, and Guangdong Academy of Medical Sciences.

Footnote

Reporting Checklist: The authors have completed the MDAR reporting checklist. Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2025-708/rc

Data Sharing Statement: Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2025-708/dss

Peer Review File: Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2025-708/prf

Funding: This study was funded by the High-Level Hospital Construction Project (grant No. DFJH201809, to J.J.Y.), the Beijing Xisike Clinical Oncology Research Foundation (grant No. Y-2022METAZZD-0110, to J.J.Y.), the National Natural Science Foundation of China (grant No. 81972164, to J.J.Y.), the Provincial Natural Science Foundation of Guangdong Province, China (grant No. 2019A1515010931, to J.J.Y.), the Natural Science Foundation of Hunan Province, China (grant No. 2021JJ30941, to J.L.), the Postdoctoral Foundation of Guangdong Provincial People’s Hospital, China (grant No. BY012022043, to J.Y.), and the Key Lab System Project of Guangdong Science and Technology Department-Guangdong Provincial Key Laboratory of Translational Medicine in Lung Cancer, China (grant No. 2017B030314120, to Yi-Long Wu).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2025-708/coif). H.R.L. and W.W. are employed by Berry Oncology Corporation. The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Research Ethics Committee of Guangdong Provincial People’s Hospital [No. GDREC2019217H (R1)] and informed consent was obtained from all individual participants.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Skoulidis F, Heymach JV. Co-occurring genomic alterations in non-small-cell lung cancer biology and therapy. Nat Rev Cancer 2019;19:495-509. [Crossref] [PubMed]
2. Yang Z, Yang N, Ou Q, et al. Investigating Novel Resistance Mechanisms to Third-Generation EGFR Tyrosine Kinase Inhibitor Osimertinib in Non-Small Cell Lung Cancer Patients. Clin Cancer Res 2018;24:3097-107. [Crossref] [PubMed]
3. Choudhury NJ, Marra A, Sui JSY, et al. Molecular Biomarkers of Disease Outcomes and Mechanisms of Acquired Resistance to First-Line Osimertinib in Advanced EGFR-Mutant Lung Cancers. J Thorac Oncol 2023;18:463-75. [Crossref] [PubMed]
4. Lara-Mejía L, Cardona AF, Mas L, et al. Impact of Concurrent Genomic Alterations on Clinical Outcomes in Patients With ALK-Rearranged NSCLC. J Thorac Oncol 2024;19:119-29. [Crossref] [PubMed]
5. Barnet MB, O'Toole S, Horvath LG, et al. EGFR-Co-Mutated Advanced NSCLC and Response to EGFR Tyrosine Kinase Inhibitors. J Thorac Oncol 2017;12:585-90. [Crossref] [PubMed]
6. Kobayashi S, Canepa HM, Bailey AS, et al. Compound EGFR mutations and response to EGFR tyrosine kinase inhibitors. J Thorac Oncol 2013;8:45-51. [Crossref] [PubMed]
7. Li W, Bai R, Guo H, et al. Epidermal growth factor receptor compound and concomitant mutations: advances in precision treatment strategies. Chin Med J (Engl) 2023;136:2776-86. [Crossref] [PubMed]
8. Cai W, Lin D, Wu C, et al. Intratumoral Heterogeneity of ALK-Rearranged and ALK/EGFR Coaltered Lung Adenocarcinoma. J Clin Oncol 2015;33:3701-9. [Crossref] [PubMed]
9. Chen J, Facchinetti F, Braye F, et al. Single-cell DNA-seq depicts clonal evolution of multiple driver alterations in osimertinib-resistant patients. Ann Oncol 2022;33:434-44. [Crossref] [PubMed]
10. Qiao M, Zhou F, Liu X, et al. Interleukin-10 induces expression of CD39 on CD8+T cells to potentiate anti-PD1 efficacy in EGFR-mutated non-small cell lung cancer. J Immunother Cancer 2022;10:e005436. [Crossref] [PubMed]
11. Yang L, He YT, Dong S, et al. Single-cell transcriptome analysis revealed a suppressive tumor immune microenvironment in EGFR mutant lung adenocarcinoma. J Immunother Cancer 2022;10:e003534. [Crossref] [PubMed]
12. Jiang B, Hu L, Dong D, et al. TP53 or CDKN2A/B covariation in ALK/RET/ROS1-rearranged NSCLC is associated with a high TMB, tumor immunosuppressive microenvironment and poor prognosis. J Cancer Res Clin Oncol 2023;149:10041-52. [Crossref] [PubMed]
13. Lv W, Liang L, Chen X, et al. Noninvasive Prenatal Testing of Methylmalonic Acidemia cblC Type Using the cSMART Assay for MMACHC Gene Mutations. Front Genet 2021;12:750719. [Crossref] [PubMed]
14. Wu F, Fan J, He Y, et al. Single-cell profiling of tumor heterogeneity and the microenvironment in advanced non-small cell lung cancer. Nat Commun 2021;12:2540. [Crossref] [PubMed]
15. Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 2013;29:15-21. [Crossref] [PubMed]
16. Bost P, Giladi A, Liu Y, et al. Host-Viral Infection Maps Reveal Signatures of Severe COVID-19 Patients. Cell 2020;181:1475-1488.e12. [Crossref] [PubMed]
17. Otsu N. A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man, and Cybernetics 1979;9:62-6.
18. Satija R, Farrell JA, Gennert D, et al. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol 2015;33:495-502. [Crossref] [PubMed]
19. Bassez A, Vos H, Van Dyck L, et al. A single-cell map of intratumoral changes during anti-PD1 treatment of patients with breast cancer. Nat Med 2021;27:820-32. [Crossref] [PubMed]
20. Gulati GS, Sikandar SS, Wesche DJ, et al. Single-cell transcriptional diversity is a hallmark of developmental potential. Science 2020;367:405-11. [Crossref] [PubMed]
21. Qiu X, Hill A, Packer J, et al. Single-cell mRNA quantification and differential analysis with Census. Nat Methods 2017;14:309-15. [Crossref] [PubMed]
22. Drake JM, Paull EO, Graham NA, et al. Phosphoproteome Integration Reveals Patient-Specific Networks in Prostate Cancer. Cell 2016;166:1041-54. [Crossref] [PubMed]
23. Gong J, Li R, Chen Y, et al. HCC subtypes based on the activity changes of immunologic and hallmark gene sets in tumor and nontumor tissues. Brief Bioinform 2021;22:bbaa427. [Crossref] [PubMed]
24. Cooperberg MR, Erho N, Chan JM, et al. The Diverse Genomic Landscape of Clinically Low-risk Prostate Cancer. Eur Urol 2018;74:444-52. [Crossref] [PubMed]
25. Liberzon A, Birger C, Thorvaldsdóttir H, et al. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 2015;1:417-25. [Crossref] [PubMed]
26. de Mey S, Dufait I, De Ridder M. Radioresistance of Human Cancers: Clinical Implications of Genetic Expression Signatures. Front Oncol 2021;11:761901. [Crossref] [PubMed]
27. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell 2011;144:646-74. [Crossref] [PubMed]
28. Hanahan D. Hallmarks of Cancer: New Dimensions. Cancer Discov 2022;12:31-46. [Crossref] [PubMed]
29. Hanahan D, Weinberg RA. The hallmarks of cancer. Cell 2000;100:57-70. [Crossref] [PubMed]
30. Yang J, Qin G, Luo M, et al. Reciprocal positive regulation between Cx26 and PI3K/Akt pathway confers acquired gefitinib resistance in NSCLC cells via GJIC-independent induction of EMT. Cell Death Dis 2015;6:e1829. [Crossref] [PubMed]
31. Yuan Y, Liao H, Pu Q, et al. miR-410 induces both epithelial-mesenchymal transition and radioresistance through activation of the PI3K/mTOR pathway in non-small cell lung cancer. Signal Transduct Target Ther 2020;5:85. [Crossref] [PubMed]
32. Li M, Wang Y, Li M, et al. Integrins as attractive targets for cancer therapeutics. Acta Pharm Sin B 2021;11:2726-37. [Crossref] [PubMed]
33. Hong S, Chen N, Fang W, et al. Upregulation of PD-L1 by EML4-ALK fusion protein mediates the immune escape in ALK positive NSCLC: Implication for optional anti-PD-1/PD-L1 immune therapy for ALK-TKIs sensitive and resistant NSCLC patients. Oncoimmunology 2016;5:e1094598. [Crossref] [PubMed]
34. Chen N, Fang W, Zhan J, et al. Upregulation of PD-L1 by EGFR Activation Mediates the Immune Escape in EGFR-Driven NSCLC: Implication for Optional Immune Targeted Therapy for NSCLC Patients with EGFR Mutation. J Thorac Oncol 2015;10:910-23. [Crossref] [PubMed]
35. Jiang VC, Hao D, Jain P, et al. TIGIT is the central player in T-cell suppression associated with CAR T-cell relapse in mantle cell lymphoma. Mol Cancer 2022;21:185. [Crossref] [PubMed]
36. Zheng L, Qin S, Si W, et al. Pan-cancer single-cell landscape of tumor-infiltrating T cells. Science 2021;374:abe6474. [Crossref] [PubMed]
37. Pabla S, Conroy JM, Nesline MK, et al. Proliferative potential and resistance to immune checkpoint blockade in lung cancer patients. J Immunother Cancer 2019;7:27. [Crossref] [PubMed]
38. Lu X, Pang Y, Cao H, et al. Integrated Screens Identify CDK1 as a Therapeutic Target in Advanced Gastrointestinal Stromal Tumors. Cancer Res 2021;81:2481-94. [Crossref] [PubMed]
39. Tsutsumi H, Iwama E, Ibusuki R, et al. Mutant forms of EGFR promote HER2 trafficking through efficient formation of HER2-EGFR heterodimers. Lung Cancer 2023;175:101-11. [Crossref] [PubMed]
40. Vokes NI, Chambers E, Nguyen T, et al. Concurrent TP53 Mutations Facilitate Resistance Evolution in EGFR-Mutant Lung Adenocarcinoma. J Thorac Oncol 2022;17:779-92. [Crossref] [PubMed]
41. Stockhammer P, Grant M, Wurtz A, et al. Co-Occurring Alterations in Multiple Tumor Suppressor Genes Are Associated With Worse Outcomes in Patients With EGFR-Mutant Lung Cancer. J Thorac Oncol 2024;19:240-51. [Crossref] [PubMed]
42. Yu X, Zhang X, Zhang Z, et al. First-generation EGFR tyrosine kinase inhibitor therapy in 106 patients with compound EGFR-mutated lung cancer: a single institution's clinical practice experience. Cancer Commun (Lond) 2018;38:51. [Crossref] [PubMed]
43. Lin YT, Tsai TH, Wu SG, et al. Complex EGFR mutations with secondary T790M mutation confer shorter osimertinib progression-free survival and overall survival in advanced non-small cell lung cancer. Lung Cancer 2020;145:1-9. [Crossref] [PubMed]
44. Blakely CM, Watkins TBK, Wu W, et al. Evolution and clinical impact of co-occurring genetic alterations in advanced-stage EGFR-mutant lung cancers. Nat Genet 2017;49:1693-704. [Crossref] [PubMed]
45. Zhao W, Song A, Xu Y, et al. Rare mutation-dominant compound EGFR-positive NSCLC is associated with enriched kinase domain-resided variants of uncertain significance and poor clinical outcomes. BMC Med 2023;21:73. [Crossref] [PubMed]
46. Wu YL, Cheng Y, Zhou J, et al. Tepotinib plus gefitinib in patients with EGFR-mutant non-small-cell lung cancer with MET overexpression or MET amplification and acquired resistance to previous EGFR inhibitor (INSIGHT study): an open-label, phase 1b/2, multicentre, randomised trial. Lancet Respir Med 2020;8:1132-43. [Crossref] [PubMed]
47. Wu YL, Zhang L, Kim DW, et al. Phase Ib/II Study of Capmatinib (INC280) Plus Gefitinib After Failure of Epidermal Growth Factor Receptor (EGFR) Inhibitor Therapy in Patients With EGFR-Mutated, MET Factor-Dysregulated Non-Small-Cell Lung Cancer. J Clin Oncol 2018;36:3101-9. [Crossref] [PubMed]
48. Yang JJ, Fang J, Shu YQ, et al. A phase Ib study of the highly selective MET-TKI savolitinib plus gefitinib in patients with EGFR-mutated, MET-amplified advanced non-small-cell lung cancer. Invest New Drugs 2021;39:477-87. [Crossref] [PubMed]
49. Fang MM, Cheng JT, Chen YQ, et al. Molecular features and clinical outcomes of EGFR-mutated, MET-amplified non-small-cell lung cancer after resistance to dual-targeted therapy. Ther Adv Med Oncol 2024;16:17588359241234504. [Crossref] [PubMed]
50. Yu HA, Goldberg SB, Le X, et al. Biomarker-Directed Phase II Platform Study in Patients With EGFR Sensitizing Mutation-Positive Advanced/Metastatic Non-Small Cell Lung Cancer Whose Disease Has Progressed on First-Line Osimertinib Therapy (ORCHARD). Clin Lung Cancer 2021;22:601-6. [Crossref] [PubMed]
51. Oxnard GR, Yang JC, Yu H, et al. TATTON: a multi-arm, phase Ib trial of osimertinib combined with selumetinib, savolitinib, or durvalumab in EGFR-mutant lung cancer. Ann Oncol 2020;31:507-16. [Crossref] [PubMed]
52. Yang JC, Shepherd FA, Kim DW, et al. Osimertinib Plus Durvalumab versus Osimertinib Monotherapy in EGFR T790M-Positive NSCLC following Previous EGFR TKI Therapy: CAURAL Brief Report. J Thorac Oncol 2019;14:933-9. [Crossref] [PubMed]
53. Atagi S. EGFR TKIs and Immune Checkpoint Inhibitors: Is This an Optimal Combination? J Thorac Oncol 2018;13:1245-7. [Crossref] [PubMed]
54. Kim TJ, Hong SA, Kim O, et al. Changes in PD-L1 expression according to tumor infiltrating lymphocytes of acquired EGFR-TKI resistant EGFR-mutant non-small-cell lung cancer. Oncotarget 2017;8:107630-9. [Crossref] [PubMed]
55. Fujimoto D, Uehara K, Sato Y, et al. Alteration of PD-L1 expression and its prognostic impact after concurrent chemoradiation therapy in non-small cell lung cancer patients. Sci Rep 2017;7:11373. [Crossref] [PubMed]