Identification of fatty acid-related subtypes, the establishment of a prognostic signature, and immune infiltration characteristics in lung adenocarcinoma

Materials and Methods

Data source and preprocessing

The data of 539 patients diagnosed with LUAD were retrieved from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/), including RNA sequencing transcriptome (TPM format), mutation, copy number variation (CNV) and clinical data. The microarray data (GSE31210, GSE68465 and GSE72094) of three LUAD patient cohorts were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). Patients with incomplete survival information were excluded from our cohort. The batch effect caused by non-biotechnology deviation among different data sets was reduced by using ‘SVA’ R package’s ‘combat’ algorithm [17]. In total, 503 patients were used from the TCGA database as training sets, and 1066 patients from the GEO database were included in our study as external validation sets. The baseline clinical information of patients was presented in Supplementary Table 1. The procedure for this study is exemplified in Figure 1.

Figure 1. Schematic representation of the present study.

Unsupervised cluster analysis

Gene sets related to FA metabolic processes were downloaded from the molecular marker database (MsigDB, https://www.GSEA-msigdb.org/GSEa/msigdb/genesets.jsp) as FARGs. We identified 38 overlapping FARGs (Supplementary Table 2) in all datasets and extracted the expression level of FARGs from each case for further analysis. Using unsupervised cluster analysis, different FA subtypes were determined according to 38 FARGs, and all patients were classified. The ‘ConsensusClusterPlus’ R package was employed to perform these analyses and it was repeated 1000 times to ensure the stability of clustering. The optimal number of clusters was determined according to the consistent clustering algorithm.

Gene set variation analysis (GSVA) and single sample gene set enrichment analysis (ssGSEA)

In order to study the differences in biological processes leading to FA characteristic patterns, the ‘GSVA’ R package was used to perform GSVA [18, 19]. Additionally, GSVA was also performed by downloading ‘c2.cp.kegg.v6.2.symbols’ from MsigDB. Calibration p<0.05 was considered as the statistical significance of the ‘limma’ software package among different subgroups. Simultaneously, the tumor-infiltrating immune cells in each sample were quantified using ssGSEA algorithm of the ‘GSVA’ R package. In addition, the distribution of LUAD patients according to the mRNA expression level of FARGs was shown by using principal component analysis (PCA).

Determination and annotation of differentially expressed genes (DEGs) in different FA metabolic patterns

‘limma’ package was used to obtain DEGs among patients with different FA subtypes and the significance criteria for selecting DEGs were set as false discovery rate (FDR) <0.05 and | log2 fold change (FC) | ≥ 1. Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) tools were then used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathway analysis to conclude the potential function of these DEGs. FDR <0.05 was considered to be statistically significant.

Construction and validation of FA metabolism score

Univariate Cox analysis was used to distinguish DEGs related to prognosis. The generation of gene clusters and evaluation of their stability was done by the unsupervised clustering algorithm. Furthermore, the LASSO method in the ‘survminer’ package was used to deal with them to avoid overfitting and deleting those closely related genes and the minimum penalty item (λ) was selected by utilizing five-fold cross-validation. The FA-related prognostic marker was established for LUAD patients, and the calculation formula for FA score was as follows: FA score= e^{sum(UBE2S expression* 0.150928887713532+ HMMR expression* 0.153839775933459+ TMPRSS11E expression* 0.0892095657835886+CHIT1 expression *-0.0625439224900231 +KRT6A expression* 0.125096336613015)}. Two groups of patients were created following the median FA score: the high FA score group and the low FA score group. The Kaplan-Meier (K-M) survival curve was drawn, and the differences between various FA score groups were evaluated with a log-rank test. The ‘time-dependent receiver operating characteristic (time-ROC)’ package was used to construct the time-ROC curve, which was then validated in three GEO cohorts. The additional prognostic value of the FA score was further confirmed by stratified analysis, and univariate and multivariate cox analyses were both used to identify the independent prognostic indicators of LUAD.

Construction of nomogram

According to the results of univariate and multivariate cox analysis, the ‘rma’ R package was used to construct a comprehensive nomogram of independent factors to quantitatively evaluate the prognosis of FA score [20]. The nomogram was tested for its accuracy by constructing a calibration curve of 1, 3 and 5 years. The predictive ability of the nomogram was evaluated by using Concordance index (C-index) curve and the time-ROC curve. The nomogram’s prognostic value and FA score were compared according to the C-index [21]. In addition, the net benefit of the nomogram was measured by decision curve analysis (DCA).

Analysis of genome variation

In order to explore the somatic mutation of FA score, the waterfall diagram was drawn with the ‘maftools’ R package to show the mutation of different FA score groups in LUAD patients. The TMB value reflecting the total mutation variable of each LUAD patient was calculated by non-synonymous mutation, and 38MB was used as the estimated value of exon size [22, 23]. In addition, ‘maftools’ R package was employed for the analysis of the significantly mutated genes and the interaction of gene mutations between high and low FA score groups. In these two analyses, only genes with more than 30 mutations in at least one group were considered. The GISTIC2.0 algorithm was employed to study the changes in somatic copy number between two different FA score groups. The fraction genome altered (FGA), the fraction of genome gained (FGG) and the fraction of the genome lost (FGL) values of each LUAD sample was determined [24]. The gene position on the chromosome was located using ‘RCirco’ R package.

Analysis of tumor immune characteristics and pathway enrichment of FA score

The evaluation of tumor immune characteristics includes two aspects: (1) The expression level of immune checkpoints [25]; (2) The score and degree of infiltration of infiltrative immunity and stromal cells calculated by Estimation of Stromal and Immune cells in malignant Tumor tissues using Expression data (ESTIMATE) [26] and ssGSEA algorithm [27]. GSEA was employed to enrich potential pathways related to FA score, whereas GSVA was utilized to verify the tumor pathway differences between high and low FA score groups. FDR values <0.05 were considered significant enrichment.

Prediction of immunotherapy response and chemosensitivity

The clinical response of patients to immune checkpoint inhibitors (ICIs) was predicted using Tumor Immune Dysfunction and Exclusion (TIDE) algorithm [28]. A high TIDE score indicated a worse response to immunotherapy. The calculation process of immunophenoscore (IPS) was described in the previous article [29], where we observed a higher score indicated a better effect of immune checkpoint inhibition treatment. We can obtain IPS from the cancer immune Atlas website (TCIA, https://tcia.at/home). The sensitivity of patients to chemotherapeutic drugs was predicted using CTRP2.0 and PRISM databases, containing the data of area under the curve (AUC) of drug sensitivity, which is used as the standard to measure the drug sensitivity by these two databases. A lower value of the AUC indicates a higher sensitivity to the treatment [30]. Additionally, for the evaluation of the FA score’s potential to predict the response to immunotherapy, an independent anti-PD-L1 immunotherapy cohort (IMvigor210) was included in this study [31].

Acquisition of tissue samples and cell culture

Fresh cancer and adjacent tissue samples from 30 LUAD patients, who were not treated before surgery, were collected from the Fujian Provincial Hospital. The consent forms were signed by all patients after the specimen extraction, and the ethical approval of the Fujian Provincial Hospital (Ethics Approval Number K2022-05-019) was taken for our research plan. Each specimen was placed in a centrifuge tube containing RNA preservation solution and stored at -80° C. All cell lines were bought from the cell bank of the Chinese Academy of Sciences. BEAS-2B was cultured in high sugar DMEM (Sigma-Aldrich, USA) containing 10% FBS (GIBCO, USA). RPMI-1640 medium (Sigma-Aldrich, USA), containing FBS was used to culture A549, HCC827, and BEAS-2B cells in a humid environment containing with 5% CO2 at 37° C.

Extraction of total RNA and analysis of real-time quantitative polymerase chain reaction (RT-qPCR)

The Trizol reagent (Invitrogen, USA) was used for total RNA extraction from tissues and cell lines following the manufacturer’s instructions. cDNA synthesis was carried out by RNA reverse transcription using PrimeScript RT kit (Promega Corporation, Madison, WI, USA). RT-qPCR (Roche, Germany) was performed on the Roche LightCycler480 II system according to the instructions of Promega SYBR-Green PCR Master Mix (Promega Corporation, Madison, WI, USA). The determination was carried out on a 96-wells plate, and each sample had three duplicate holes. GAPDH was used as the internal reference gene, and the relative expression level of mRNA was calculated using 2^-ΔΔCT method. The primer sequence of RT-qPCR is shown in Supplementary Table 3.

Western blot analysis

The Western blot analysis was conducted following the established protocol [32]. Tissues and cells underwent homogenization in radioimmunoprecipitation assay buffer (RIPA, Solarbio, China, R0010), supplemented with 1% protease inhibitor cocktail, and the supernatant concentration was ascertained utilizing the bicinchoninic acid assay. Target proteins were loaded onto a sodium dodecyl sulfate-polyacrylamide gel, followed by protein transfer to a polyvinylidene fluoride membrane. The membranes were then incubated overnight at 4° C with the following primary antibodies: KRT6A (ABclonal, China, 1:1000), TMPRSS11E (Thermo Fisher Scientific, USA, 1:1000), HMMR (ABclonal, China, 1:1000), UBE2S (ABclonal, China, 1:1000) CHIT1 (ABclonal, China, 1:1000) and GAPDH (ABclonal, China, 1:1000) antibodies after a 1.5 hour blocking period. Upon washing with Tris-Buffered Saline with EDTA and Tween 20, the membranes were exposed to horseradish peroxidase-conjugated secondary antibodies for 1 hour at ambient temperature. The blots were developed using an enhanced chemiluminescence substrate (ECL reagents) and subsequently analyzed by Image Lab image analysis software.

Cell transfection

The HanBio Company (Fuzhou, Fujian, China) provided KRT6A small interfering RNAs (siRNA-KRT6A) and non-target small interfering RNA (siRNA-control). Two siRNA and si-control sequences are presented in Supplementary Table 4. Following the manufacturer’s instructions, transfection was carried out using HanBio RNA-specific transfection reagent (Invitrogen, Carlsbad, CA, USA) in Opti-MEM medium (Gibco, Rockville, MD, USA). mRNA expression levels after 48 hours and protein expression levels after 72 hours were carried out to assess the efficiency of cell transfection. After stable transcription, the cells were collected for the other cell-related experiments.

CCK8 (cell counting kit-8) assay

The cell proliferation ability after transfection was studied by CCK8. Briefly, the KRT6A-knockdown cells in the logarithmic phase of growth were seeded into 96-well plates with 2000 cells per well. After the cells adhered, 10μl CCK8 reagent (Cellcook, Guangzhou, China) was added to the into wells and OD (optical density) was measured after incubation at 37° C and 5% CO2 for 2 hours. The detection time was set as 0, 24, 48 and 72 hours respectively. The absorbance at 450nm for each well was read using a microplate reader.

Clone formation assay

The transfected cell suspension was seeded in a 6-well plate at 1000 cells/well. Cells in 6-well plates were then grown in a complete medium. The cells were cultured for 10 days, the cells were fixed with 4% paraformaldehyde for 30 minutes and stained with 0.1% crystal violet for 30 minutes. The formed cell colonies were observed and photographed with a microscope. The number of colonies formed was counted by ImageJ.

5-Ethynyl-20-deoxyuridine (EdU) assay

Employing an EdU Kit (Beyotime, China, C0071s), the EdU assay was conducted. Approximately 60,000 cells per well were distributed in 24-well plates and subsequently cultured with EdU reagent at a 1:1000 dilution for 2 hours on the subsequent day. The cells were then fixed utilizing 4% paraformaldehyde and stained with fluorescent dye and Hoechst. Post-staining, the cells were imaged and enumerated under a fluorescence microscope. The EdU-positive cell index was ascertained by calculating the proportion of positive cells relative to the total cell count.

Transwell assay

To determine the migratory capacity of cells, 100μl of the cell suspension (50000 cells) with serum-free medium were seeded into the upper chamber of a transwell plate. Next, 500μl complete 15% FBS-containing medium was added to the lower compartments. After incubation for 24 hours, the cells were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet for 30 minutes respectively. Then a cotton swab was used on the cells in the upper chambers. Eventually, the cells on the lower surface of the membrane were observed and photographed by the microscope, and the results were measured and recorded using ImageJ. For the invasion assay, the upper chamber was coated with 100μl of 10% Matrigel. The rest of the procedure was the same as migration assay.

Statistical analysis

R software (version 4.2.1), GraphPad Prism (version 9.0) and SangerBox web tool (http://sangerbox.com/) were used to perform all the statistical analyses. Some standard tests which included Student’s T test, Wilcoxon rank sum test, chi-square test, and Fisher exact test, were employed to determine the differences in variables among different groups. The p-value for multivariate multiple comparisons underwent adjustment employing the Benjamin-Hochberg method, and all tests utilized two-tailed p-values [33]. Moreover, the statistical significance was set to p <0.05.

Results

Genetic and transcriptional changes of FARGs in LUAD

A total of 38 FARGs were enrolled in this study and the incidence of somatic mutations analysis in these 38 FARGs revealed that the overall frequency of mutation was relatively low in the LUAD cohort (Figure 2A). Only 171 of 503 patients had mutations (Figure 2A), among which ACSL1 had the highest mutation frequency (3%), while 14 FARGs had no mutation. Furthermore, the somatic CNV in these FARGs was studied and it was found that CNV was common in all 38 FARGs (Figure 2B). ALDH9A1 and ACOX1 exhibited a significant increase in CNV, whereas ACADM displayed a decrease in CNV. The location of CNV changes in FARGs on their respective chromosomes is illustrated in Figure 2C. Additionally, the level of expression of all FARGs in tumor and non-tumor tissues were statistically different (Figure 2D). In univariate Cox analysis, 18 FARGs related to prognosis were detected. Figure 2E illustrates the network of FARGs interactions, prognostic values, and regulator connections in patients with LUAD.

Figure 2. Genetic and transcriptional alterations of FARGs in TCGA-LUAD patients. (A) Mutation frequencies of 38 FARGs. (B) CNV alteration frequency of 38 FARGs. (C) Location of the CNV alteration of 38 FARGs on chromosomes. (D) Expression of 38 FARGs between normal and tumor tissues. (E) The interactions among FARGs in LUAD are visually represented in the diagram. The interconnecting lines between FARGs signify their interdependence, with the line width representing the strength of the correlation between FARGs. Negative correlations are illustrated in green, while positive correlations are denoted in pink. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Identification of FA subtypes

A total of 1569 samples from four independent LUAD cohorts were retained with complete survival information. Through unsupervised clustering of 38 FARGs expressions, the whole cohort was divided into two subtypes: subtype A and B (Supplementary Figure 1 and Supplementary Table 5). Through PCA analysis, significant differences in FARGs transcriptional profiles were revealed between the two subtypes (Figure 3A). Survival analysis revealed that patients with subtype A had a significantly smaller OS than those with subtype B (Figure 3B). Moreover, the comparison of the clinical characteristics of different subtypes of LUAD indicated significant differences in the FARGs expression and clinical pathological characteristics, and the subtype of patients was related to age and T stage (Figure 3C). The KEGG enrichment pathway and the immune cell composition of tumor microenvironment (TME) between the two subtypes were analyzed to understand the association between intrinsic biological traits and various clinical phenotypes. GSVA revealed that subtype A was enriched in the cell process, while subtype B was mainly related to the metabolism of various substances (Figure 3D and Supplementary Table 6). According to the results of the ssGSEA analysis, subtype B had considerably higher levels of infiltration of B cells, DC cells, eosinophils, T cells, and Th cells than subtype A, while subtype A had significantly higher levels of infiltration of NK cells, Treg cells, and Th2 cells than subtype B (Figure 3E and Supplementary Table 7).

Figure 3. Consistent clustering for FARGs clusters, biological processes, and characteristics of immune infiltration cells of each cluster. (A) PCA plot based on the FARGs. (B) Kaplan–Meier survival curves for the different FARGs clusters. (C) Heatmap for different clinicopathologic features and expression levels of FARGs between two clusters. (D) GSVA of biological pathways between two clusters. (E) Different expressions of immune infiltration cells in each FARGs cluster. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Gene cluster identification based on DEGs

To investigate the potential biological behavior inherent in each FA pattern, we identified 222 DEGs related to FARGs subtypes. Subsequently, we conducted a meticulous functional enrichment analysis to gain a deeper understanding of these DEGs (Supplementary Table 8). These DEGs were then enriched and analyzed. The FARGs were significantly enriched in the process of metabolism and cell cycle (Figure 4A, 4B and Supplementary Table 9), indicating that FA plays a significant role in cell cycle and metabolism. The Cox regression analysis was then used to detect the relationship between these FARGs and OS in LUAD patients, and 197 FARGs with strong prognostic values were identified to help with further analysis of the core characteristics of FA (Supplementary Table 10). To verify this conclusion, we undertook an additional unsupervised clustering step on a set of 197 FARGs, which resulted in the identification of three distinct gene clusters for FARGs (Supplementary Figure 2 and Supplementary Table 11). Obviously, significant differences in the mRNA expression of FARGs were observed between gene clusters in the metadata set (Figure 4C). The PCA map displays the relative distances of the three gene clusters using current prognostic DEGs (Figure 4D). Similar to FARGs subtyping, gene clusters can also be utilized to differentiate immune-infiltrating cells and OS in all patients. (Figure 4E, 4F).

Figure 4. Identification of gene clusters based on DEGs. (A, B) GO and KEGG enrichment analyses of DEGs among two FARGs clusters. (C) Expression of FARGs between gene clusters A, B and C. (D) PCA plot based on the DEGs. (E) Kaplan–Meier curves for the different gene clusters. (F) Different expressions of immune infiltration cells in each gene cluster. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Construction and validation of prognostic FARGs score

In the establishment of a model for quantification of each patient, five of the 197 prognostic DEGs were maintained by applying multivariate Cox regression analysis and the LASSO regression model. These were used to establish FARGs scores, named ‘FA score’ (Supplementary Figure 3 and Supplementary Table 12). Furthermore, the value of the FA score was determined by the prediction of the prognosis of patients. The patients were divided into two groups: the high FA score group and the low FA score group according to the median score. Figure 5A illustrates the distribution of patients in two FA subtypes, three gene phenotypes and two FA score groups. There were differences detected in relation to expressions of the great majority of FARGs in the high and low FA score groups (Figure 5B). Besides, significant differences in FA score among FA clusters and gene clusters were identified (Figure 5C, 5D). We observed that the low FA score group had a significant survival advantage over the group with a high FA score (Figure 5E). To showcase the comprehensive significance of the FA score, its validation across supplementary cohorts produced consistent results (Supplementary Figure 4A–4C and Supplementary Table 12). Additionally, with the use of time-ROC analysis, it was further verified that the FA score is a good indicator for predicting the prognosis of LUAD patients (Figure 5F and Supplementary Figure 4D–4F). In order to test the good applicability of the FA score, a stratified analysis of LUAD cancer patients was performed based on clinical and pathological data. The comparison between patients with high and low FA scores revealed that individuals with high FA scores had a shorter OS period at most levels including age over 60 years (Supplementary Figure 5). However, it may be because of the small sample size, and no statistical difference in OS of M1 patients of the high and low FA score groups. The correlation analysis between FA score and several clinical features showed that the survival status, M, N, T, age and stage of LUAD patients were notably associated with FA score (Figure 5G) and these findings indicate that FA score is a highly reliable marker. As a result, the FA score was included in Cox regression analysis as an effective index and clinical characteristics, and we observed that FA score, N, and T are independent variables that impact the prognosis of LUAD patients (Table 1).

Figure 5. Construction of FA score. (A) Alluvial diagram of clusters distributions in groups with different FA score and survival outcomes. (B) Expression of FARGs between high and low FA score groups. (C) Differences in FA score between FARGs clusters. (D) Differences in FA score between gene clusters. (E) Kaplan–Meier curves for patients with high and low FA score in the TCGA-LUAD cohort. (F) Time-dependent receiver operating characteristic curve of FA score for predicting the prognosis of the LUAD patients in the TCGA-LUAD. (G) The circular pie chart for the proportion difference of clinical indices. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Table 1. Univariate and multivariate Cox analysis of the clinicopathological features and FA score with OS.

Characteristics	Univariate Cox		Multivariate Cox
	HR(95%CI)	P value	HR(95%CI)	P value
Stage	1.977(1.586-2.463)	< 0.001	1.256(0.883-1.786)	0.205
M	1.727(1.18-2.527)	0.005	1.184(0.75-1.87)	0.469
N	1.942(1.575-2.394)	< 0.001	1.584(1.193-2.104)	0.001
T	1.816(1.386-2.38)	< 0.001	1.495(1.067-2.094)	0.02
Age	1.038(0.822-1.31)	0.754
Sex	1.041(0.847-1.28)	0.7
FA score	0.579(0.467-0.717)	< 0.001	0.695(0.535-0.903)	0.006
Significant value is given in bold.

Construction of nomogram for predicting the survival rate

Considering the inconvenience of using FA score to predict OS in LUAD patients clinically, a nomogram was established featuring the FA score and clinicopathological parameters to predict the OS rate of patients (Figure 6A). Predictive factors included independent factors N, T, and FA score, and the calibration chart illustrated that the proposed nomogram had a similar performance in the TCGA queue compared with the ideal model (Figure 6B). The time-ROC and C-index curves both suggested that the nomogram had the best prediction effect (Figure 6C, 6D). In addition, DCA disclosed that the nomogram garnered greater net advantages in predicting 1, 3, and 5-year prognostic outcomes compared to univariate analysis (Figure 6E, 6F). More importantly, these results were validated in the independent cohort GSE68465 (Supplementary Figure 6).

Figure 6. Development of a nomogram by integrating the FA score and clinicopathological characters in TCGA-LUAD cohort. (A) Nomogram for predicting the 1-, 3-, 5-years OS. (B) Calibration curve of the nomogram for predicting the 1-, 3-, and 5-years OS. (C) ROC curve for predicting the different years’ OS. (D–F) Decision curves showing the comparison of net benefits of the nomogram, N, T and FA score for 1-, 3-, and 5-years OS.

Correlation between genomic changes and FA score

The distribution of somatic variation in LUAD driver genes was evaluated across groups with high and low FA scores. The top 20 driver genes with the highest frequency of change were analyzed and evaluated using the ‘maftools’ tool. (Figure 7A, 7B). These findings provide a new direction for studying the composition of tumor ICIs and the mechanism of gene mutation in immune checkpoint blockade (ICB) treatment. When comparing the mutation frequency between samples from high and low FA score groups, more somatic mutations, including non-synonymous mutations and synonymous mutations, were observed in the high FA score group (Figure 7C–7E). The TCGA cohort mutation annotation file analysis revealed a substantial difference in the mutation frequency of 37 genes, including TP53 and UBR4, and between the two FA score groups, a significant co-occurrence of mutations in these genes was observed (Figure 7F, 7G). Additionally, the patients’ TMB was considerably higher in the high FA score group than that in the low FA score group (Figure 7H and Supplementary Table 13). According to the optimal cut-off value of TMB, patients were divided into two groups. K-M survival analysis revealed that the patients in the high TMB group had a substantially ideal OS rate (Figure 7I). Subsequently, the TMB score and FA score were combined to conduct a stratified survival analysis, which showed that the high and low TMB subgroups classified in accordance with the FA score displayed significant survival differences (Figure 7J). Further analysis was performed of the GISTIC score and copy number gain/loss frequency of high and low FA score groups and the results illustrated that the frequency of copy number increase/decrease in the high FA score group was greater than that in low FA score group (Figure 7K). Furthermore, the differences in FGA, FGG, and FGL between subtypes were evaluated and we observed greater values of the FGA, FGG and FGL in the high FA score group than those in the low FA score group (Figure 7L and Supplementary Table 14). This result indicated that the increase in copy number/loss frequency might be a factor leading to the higher FA score in LUAD patients.

Figure 7. Integrated comparisons of somatic mutation and CNVs between FA score groups in the TCGA cohort. (A, B) Waterfall plots showing the mutation information of the top 20 genes with the highest mutation frequency in high and low FA score groups. (C–E) Association between all mutation counts, synonymous mutation counts, nonsynonymous mutation counts, and FA score and their distribution in the low and high FA score groups. (F) Differentially mutated genes between high and low FA score groups are displayed as a forest plot. (G) Interaction effect of genes mutating differentially in patients in the low and the high FA score groups. (H) Distribution of TMB in the low and the high FA score groups. (I) Kaplan–Meier curves for the OS of the high-TMB and low-TMB groups. (J) Kaplan–Meier curves for patients stratified by both TMB and FA score. (K) Gene fragments profiles with amplification (red) and deletion (green) among the two groups. (L) Comparison of the fraction of the genome altered, lost, and gained between the two groups. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Evaluation of TME between high FA score and low FA score patients

The effects of immunotherapy and its consequences on the cancer cells are determined by the state of TME. The correlation between FA score and cancer immune cycle activity was analyzed and the activity of some anti-cancer immune responses, such as Release of cancer cell antigens, Priming and activation, CD8 T cell recruiting, Neutrophil recruiting and MDSC recruiting, was observed to be substantially different in the high and low FA score groups (Figure 8A and Supplementary Table 15). Moreover, the association between FA score and immune cell abundance was evaluated using ssGSEA algorithm. The scatter diagram illustrates the majority of immune cells were vastly infiltrated in the low FA score group (Figure 8B). A low FA score was closely related to a high immune score, stromal score and ESTIMATE score, while a high FA score was related to high tumor purity (Figure 8C and Supplementary Table 16). Additionally, TME score exhibits a negative correlation with TMB (Figure 8D). The FA score was found to be positively correlated with the expression of many immunological checkpoints as well as the enrichment score of gene features relevant to immunotherapy response (Figure 8E and Supplementary Tables 17, 18). GSVA was conducted in the high FA score and low FA score groups to investigate the cancer marker pathways related to FA score. Compared with the low FA score group, the 17 landmark pathways in the high FA score group were substantially increased (Figure 8F and Supplementary Table 19). Moreover, GSEA corroborated that 11 oncogenic pathways were upregulated in the high FA score group, the majority of which were linked to the well-established carcinogenic pathways (Figure 8G–8J and Supplementary Table 20).

Figure 8. The relationship between immune characters of TME and FA score in the whole cohort. (A) Differences in activities of the cancer immunity cycles between high and low FA score groups. (B) Different expressions of immune infiltration cells in each FA score group. (C) Violin plots for the immune score, stromal score, ESTIMATE score, and tumor purity in the low and high FA score groups. (D) Spearman correlation analysis between TMB and immune score, stromal score, and ESTIMATE score. (E) The correlations between FA score, immunotherapy-predicted pathways and immune checkpoints. (F) The difference in the hallmark gene sets between FA score groups. (G–J) The GSEA results for the 11 overlapping upregulated hallmark pathways in terms of the high FA score groups. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Application of FA score in predicting chemotherapy and immunotherapy

IPS files downloaded from TCIA were utilized to determine whether FA score can predict the response of LUAD patients to immunotherapy. An elevated IPS value was discerned in the low FA score cohort, signifying an enhanced immunotherapeutic response among the patients (Figure 9A and Supplementary Table 21). The TIDE algorithm was also used to predict the immunotherapeutic effect of ICB which suggested that the TIDE score of the high FA score group was higher, implying that patients with high FA score may display a poor immunotherapy response (Figure 9B and Supplementary Table 22). An immunotherapy cohort (IMvigor210) was introduced to further investigate if an FA score can predict response to immunotherapy. A better survival rate was observed in patients with low FA scores, and they showed a higher objective remission rate than those with high FA score (Figure 9C–9E). While immunotherapy remains a preeminent modality for cancer management, chemotherapy has historically served as a crucial postoperative intervention. Therefore, it is necessary to forecast potential LUAD therapeutic alternatives with the CTRP and PRISM repositories. According to the results, drugs such as leptomycin B and paclitaxel were predicted to be better options for patients with high FA scores. Moreover, it was observed that out of the seven chemotherapy drugs, the value of AUC of gemcitabine was the lowest, which means that LUAD patients may have good therapeutic sensitivity to gemcitabine (Supplementary Figure 7).

Figure 9. Application of the FA score for immunotherapy prediction in the high and low FA score groups in the TCGA-LUAD cohort. (A) IPS is used to predict the responsiveness to CTLA-4 and PD-1. (B) Distribution of TIDE scores in the whole cohort. (C) Kaplan–Meier curves for patients with high and low FA score in the IMvigor210 cohort. (D) The rate of response to immunotherapy for patients with high and low FA score in the IMvigor210 cohort. (E) The distribution of FA score in different patient statuses in the IMvigor210 cohort. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

The expression level of genes in FA score in LUAD cell lines and tissues

UBE2S, HMMR, TMPRSS11E and KRT6A were found to be up-regulated in tumor tissues relative to surrounding non-cancerous and normal tissues in the screening cohort (Supplementary Figure 8). The RT-qPCR and WB for LUAD cell lines and patient tissues were used to validate the expression levels of 4 genes in LUAD. The expression of UBE2S, HMMR and KRT6A expression in tumor tissues was substantially higher than that in adjacent non-cancerous tissues (Figure 10A). For cell lines, the expression of UBE2S, HMMR, TMPRSS11E and CHIT1 was upregulated, whereas KRT6A was downregulated in LUAD cell lines compared with cell line BEAS-2B (Figure 10B).

Figure 10. The expression validation of the genes in the FA score by RT-qPCR and WB. (A) Expression levels of 5 genes in LUAD and matched paracancerous tissues. (B) Expression levels of 5 genes in the normal lung epithelial cell line BEAS-2B and two LUAD cell lines (HCC827, A549). ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. T, tumor tissues; P, paracancerous tissues.

Knockdown of KRT6A inhibits LUAD cell proliferation, migration and invasion

Considering the important role of FA score in LUAD, the independent prognostic genes may have a greater impact on the biological function of LUAD cells. We chose KRT6A, exhibiting the most robust prognostic association, to further substantiate our conjecture. To explore the role of KRT6A in LUAD in vitro, siRNAs were introduced into A549 and HCC827 cells to suppress KRT6A expression, as evidenced by RT-qPCR (Figure 11A). The CCK-8 assay revealed that suppression of KRT6A led to a reduction in the total cell viability of both A549 and HCC827 cells (Figure 11B). The results of the clone formation and EDU assays corroborated the results of CCK8 (Figure 11C, 11D). In addition, a significant decrease in the LUAD cell invasion was also shown by the transwell assay (Figure 11E).

Figure 11. KRT6A silencing inhibits LUAD cell proliferation, migration and invasion (A549 and HCC827). (A) RT-qPCR was used to verify the efficiency of KRT6A knockdown in LUAD cell lines. (B) CCK-8, (C) clone formation and (D) EDU assays were used to test the proliferation of LUAD cell lines. (E) Transwell assay was applied to measure cell migration and invasion of LUAD cell lines. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Author Contributions

Yi Huang and Yi Zhang conceived and designed the study. Yi Huang and Yuzhi Wang performed the collection and assembly of data. Yuzhi Wang and Xiaoxiao Huang analyzed the data. Yunfei Xu and Yuzhi Wang performed experiments and wrote the manuscript. Xiaoxiao Huang, Hong Fan and Zelin Qi revised the article critically. All authors read and approved the final manuscript.

Acknowledgments

We thank Bullet Edits Limited for the linguistic editing and proofreading of the manuscript.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Ethical Statement and Consent

Ethics approval was sought and approved by the Ethics Committee of Fujian Provincial Hospital (Ethics Approval Number K2022-05-019). The patients/participants provided their written informed consent to participate in this study.

Funding

This study was supported in part by Medical Vertical Project of Fujian Province (Grant No. 2020CXB001) to Yi Huang, Joint fund of science and technology innovation of Fujian province (Grant No. 2021Y9024) to Yi Huang, Key Project of Natural Science Foundation of Fujian province (Grant No. 2022J02048) to Yi Huang.

References

• 1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, Bray F. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021; 71:209–49. https://doi.org/10.3322/caac.21660 [PubMed]
• 2. Bade BC, Dela Cruz CS. Lung Cancer 2020: Epidemiology, Etiology, and Prevention. Clin Chest Med. 2020; 41:1–24. https://doi.org/10.1016/j.ccm.2019.10.001 [PubMed]
• 3. Gridelli C, Rossi A, Carbone DP, Guarize J, Karachaliou N, Mok T, Petrella F, Spaggiari L, Rosell R. Non-small-cell lung cancer. Nat Rev Dis Primers. 2015; 1:15009. https://doi.org/10.1038/nrdp.2015.9 [PubMed]
• 4. Yamanashi K, Okumura N, Takahashi A, Nakashima T, Matsuoka T. Surgical and survival outcomes of lung cancer patients with intratumoral lung abscesses. J Cardiothorac Surg. 2017; 12:44. https://doi.org/10.1186/s13019-017-0607-3 [PubMed]
• 5. Schenk EL, Patil T, Pacheco J, Bunn PA Jr. 2020 Innovation-Based Optimism for Lung Cancer Outcomes. Oncologist. 2021; 26:e454–72. https://doi.org/10.1002/onco.13590 [PubMed]
• 6. Swinnen JV, Brusselmans K, Verhoeven G. Increased lipogenesis in cancer cells: new players, novel targets. Curr Opin Clin Nutr Metab Care. 2006; 9:358–65. https://doi.org/10.1097/01.mco.0000232894.28674.30 [PubMed]
• 7. DeBerardinis RJ, Thompson CB. Cellular metabolism and disease: what do metabolic outliers teach us? Cell. 2012; 148:1132–44. https://doi.org/10.1016/j.cell.2012.02.032 [PubMed]
• 8. Santos CR, Schulze A. Lipid metabolism in cancer. FEBS J. 2012; 279:2610–23. https://doi.org/10.1111/j.1742-4658.2012.08644.x [PubMed]
• 9. Medes G, Thomas A, Weinhouse S. Metabolism of neoplastic tissue. IV. A study of lipid synthesis in neoplastic tissue slices in vitro. Cancer Res. 1953; 13:27–9. [PubMed]
• 10. Wen YA, Xing X, Harris JW, Zaytseva YY, Mitov MI, Napier DL, Weiss HL, Mark Evers B, Gao T. Adipocytes activate mitochondrial fatty acid oxidation and autophagy to promote tumor growth in colon cancer. Cell Death Dis. 2017; 8:e2593. https://doi.org/10.1038/cddis.2017.21 [PubMed]
• 11. Shi J, Fu H, Jia Z, He K, Fu L, Wang W. High Expression of CPT1A Predicts Adverse Outcomes: A Potential Therapeutic Target for Acute Myeloid Leukemia. EBioMedicine. 2016; 14:55–64. https://doi.org/10.1016/j.ebiom.2016.11.025 [PubMed]
• 12. Shao H, Mohamed EM, Xu GG, Waters M, Jing K, Ma Y, Zhang Y, Spiegel S, Idowu MO, Fang X. Carnitine palmitoyltransferase 1A functions to repress FoxO transcription factors to allow cell cycle progression in ovarian cancer. Oncotarget. 2016; 7:3832–46. https://doi.org/10.18632/oncotarget.6757 [PubMed]
• 13. Wang YY, Attané C, Milhas D, Dirat B, Dauvillier S, Guerard A, Gilhodes J, Lazar I, Alet N, Laurent V, Le Gonidec S, Biard D, Hervé C, et al. Mammary adipocytes stimulate breast cancer invasion through metabolic remodeling of tumor cells. JCI Insight. 2017; 2:e87489. https://doi.org/10.1172/jci.insight.87489 [PubMed]
• 14. Talmadge JE, Fidler IJ. AACR centennial series: the biology of cancer metastasis: historical perspective. Cancer Res. 2010; 70:5649–69. https://doi.org/10.1158/0008-5472.CAN-10-1040 [PubMed]
• 15. Schoors S, Bruning U, Missiaen R, Queiroz KC, Borgers G, Elia I, Zecchin A, Cantelmo AR, Christen S, Goveia J, Heggermont W, Goddé L, Vinckier S, et al. Fatty acid carbon is essential for dNTP synthesis in endothelial cells. Nature. 2015; 520:192–7. https://doi.org/10.1038/nature14362 [PubMed]
• 16. Wong BW, Wang X, Zecchin A, Thienpont B, Cornelissen I, Kalucka J, García-Caballero M, Missiaen R, Huang H, Brüning U, Blacher S, Vinckier S, Goveia J, et al. The role of fatty acid β-oxidation in lymphangiogenesis. Nature. 2017; 542:49–54. https://doi.org/10.1038/nature21028 [PubMed]
• 17. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007; 8:118–27. https://doi.org/10.1093/biostatistics/kxj037 [PubMed]
• 18. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
• 19. Ferreira MR, Santos GA, Biagi CA, Silva Junior WA, Zambuzzi WF. GSVA score reveals molecular signatures from transcriptomes for biomaterials comparison. J Biomed Mater Res A. 2021; 109:1004–14. https://doi.org/10.1002/jbm.a.37090 [PubMed]
• 20. Liu TT, Li R, Huo C, Li JP, Yao J, Ji XL, Qu YQ. Identification of CDK2-Related Immune Forecast Model and ceRNA in Lung Adenocarcinoma, a Pan-Cancer Analysis. Front Cell Dev Biol. 2021; 9:682002. https://doi.org/10.3389/fcell.2021.682002 [PubMed]
• 21. Eng KH, Schiller E, Morrell K. On representing the prognostic value of continuous gene expression biomarkers with the restricted mean survival curve. Oncotarget. 2015; 6:36308–18. https://doi.org/10.18632/oncotarget.6121 [PubMed]
• 22. Chalmers ZR, Connelly CF, Fabrizio D, Gay L, Ali SM, Ennis R, Schrock A, Campbell B, Shlien A, Chmielecki J, Huang F, He Y, Sun J, et al. Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. 2017; 9:34. https://doi.org/10.1186/s13073-017-0424-2 [PubMed]
• 23. Sun J, Shi R, Zhang X, Fang D, Rauch J, Lu S, Wang X, Käsmann L, Ma J, Belka C, Su C, Li M. Characterization of immune landscape in papillary thyroid cancer reveals distinct tumor immunogenicity and implications for immunotherapy. Oncoimmunology. 2021; 10:e1964189. https://doi.org/10.1080/2162402X.2021.1964189 [PubMed]
• 24. Moore LE, Jaeger E, Nickerson ML, Brennan P, De Vries S, Roy R, Toro J, Li H, Karami S, Lenz P, Zaridze D, Janout V, Bencko V, et al. Genomic copy number alterations in clear cell renal carcinoma: associations with case characteristics and mechanisms of VHL gene inactivation. Oncogenesis. 2012; 1:e14. https://doi.org/10.1038/oncsis.2012.14 [PubMed]
• 25. Johnston RJ, Su LJ, Pinckney J, Critton D, Boyer E, Krishnakumar A, Corbett M, Rankin AL, Dibella R, Campbell L, Martin GH, Lemar H, Cayton T, et al. VISTA is an acidic pH-selective ligand for PSGL-1. Nature. 2019; 574:565–70. https://doi.org/10.1038/s41586-019-1674-5 [PubMed]
• 26. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
• 27. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005; 102:15545–50. https://doi.org/10.1073/pnas.0506580102 [PubMed]
• 28. Lu X, Jiang L, Zhang L, Zhu Y, Hu W, Wang J, Ruan X, Xu Z, Meng X, Gao J, Su X, Yan F. Immune Signature-Based Subtypes of Cervical Squamous Cell Carcinoma Tightly Associated with Human Papillomavirus Type 16 Expression, Molecular Features, and Clinical Outcome. Neoplasia. 2019; 21:591–601. https://doi.org/10.1016/j.neo.2019.04.003 [PubMed]
• 29. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017; 18:248–62. https://doi.org/10.1016/j.celrep.2016.12.019 [PubMed]
• 30. Yang C, Huang X, Li Y, Chen J, Lv Y, Dai S. Prognosis and personalized treatment prediction in TP53-mutant hepatocellular carcinoma: an in silico strategy towards precision oncology. Brief Bioinform. 2021; 22:bbaa164. https://doi.org/10.1093/bib/bbaa164 [PubMed]
• 31. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE II, Koeppen H, Astarita JL, Cubas R, Jhunjhunwala S, Banchereau R, Yang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018; 554:544–8. https://doi.org/10.1038/nature25501 [PubMed]
• 32. Bu J, Qiao X, He Y, Liu J. Colonic electrical stimulation improves colonic transit in rotenone-induced Parkinson’s disease model through affecting enteric neurons. Life Sci. 2019; 231:116581. https://doi.org/10.1016/j.lfs.2019.116581 [PubMed]
• 33. Thissen D, Steinberg L, Kuang D. Quick and Easy Implementation of the Benjamini-Hochberg Procedure for Controlling the False Positive Rate in Multiple Comparisons. Journal of Educational and Behavioral Statistics. 2002; 27:77–83. https://doi.org/10.3102/10769986027001077
• 34. Hirschey MD, DeBerardinis RJ, Diehl AME, Drew JE, Frezza C, Green MF, Jones LW, Ko YH, Le A, Lea MA, Locasale JW, Longo VD, Lyssiotis CA, et al, and Target Validation Team. Dysregulated metabolism contributes to oncogenesis. Semin Cancer Biol. 2015; 35:S129–50. https://doi.org/10.1016/j.semcancer.2015.10.002 [PubMed]
• 35. Ma Y, Temkin SM, Hawkridge AM, Guo C, Wang W, Wang XY, Fang X. Fatty acid oxidation: An emerging facet of metabolic transformation in cancer. Cancer Lett. 2018; 435:92–100. https://doi.org/10.1016/j.canlet.2018.08.006 [PubMed]
• 36. Rodríguez-Enríquez S, Hernández-Esquivel L, Marín-Hernández A, El Hafidi M, Gallardo-Pérez JC, Hernández-Reséndiz I, Rodríguez-Zavala JS, Pacheco-Velázquez SC, Moreno-Sánchez R. Mitochondrial free fatty acid β-oxidation supports oxidative phosphorylation and proliferation in cancer cells. Int J Biochem Cell Biol. 2015; 65:209–21. https://doi.org/10.1016/j.biocel.2015.06.010 [PubMed]
• 37. Peng G, Chi H, Gao X, Zhang J, Song G, Xie X, Su K, Song B, Yang J, Gu T, Li Y, Xu K, Li H, et al. Identification and validation of neurotrophic factor-related genes signature in HNSCC to predict survival and immune landscapes. Front Genet. 2022; 13:1010044. https://doi.org/10.3389/fgene.2022.1010044 [PubMed]
• 38. Zhao S, Chi H, Yang Q, Chen S, Wu C, Lai G, Xu K, Su K, Luo H, Peng G, Xia Z, Cheng C, Lu P. Identification and validation of neurotrophic factor-related gene signatures in glioblastoma and Parkinson’s disease. Front Immunol. 2023; 14:1090040. https://doi.org/10.3389/fimmu.2023.1090040 [PubMed]
• 39. Shen Y, Chi H, Xu K, Li Y, Yin X, Chen S, Yang Q, He M, Zhu G, Li X. A Novel Classification Model for Lower-Grade Glioma Patients Based on Pyroptosis-Related Genes. Brain Sci. 2022; 12:700. https://doi.org/10.3390/brainsci12060700 [PubMed]
• 40. Liu Z, Xu L. UBE2S promotes the proliferation and survival of human lung adenocarcinoma cells. BMB Rep. 2018; 51:642–7. https://doi.org/10.5483/BMBRep.2018.51.12.138 [PubMed]
• 41. Ho JY, Lu HY, Cheng HH, Kuo YC, Lee YA, Cheng CH. UBE2S activates NF-κB signaling by binding with IκBα and promotes metastasis of lung adenocarcinoma cells. Cell Oncol (Dordr). 2021; 44:1325–38. https://doi.org/10.1007/s13402-021-00639-4 [PubMed]
• 42. Jiang X, Tang L, Yuan Y, Wang J, Zhang D, Qian K, Cho WC, Duan L. NcRNA-Mediated High Expression of HMMR as a Prognostic Biomarker Correlated With Cell Proliferation and Cell Migration in Lung Adenocarcinoma. Front Oncol. 2022; 12:846536. https://doi.org/10.3389/fonc.2022.846536 [PubMed]
• 43. Yang B, Zhang W, Zhang M, Wang X, Peng S, Zhang R. KRT6A Promotes EMT and Cancer Stem Cell Transformation in Lung Adenocarcinoma. Technol Cancer Res Treat. 2020. https://doi.org/10.1177/1533033820921248 [PubMed]
• 44. Erin N, Grahovac J, Brozovic A, Efferth T. Tumor microenvironment and epithelial mesenchymal transition as targets to overcome tumor multidrug resistance. Drug Resist Updat. 2020; 53:100715. https://doi.org/10.1016/j.drup.2020.100715 [PubMed]
• 45. Hong W, Yuan H, Gu Y, Liu M, Ji Y, Huang Z, Yang J, Ma L. Immune-related prognosis biomarkers associated with osteosarcoma microenvironment. Cancer Cell Int. 2020; 20:83. https://doi.org/10.1186/s12935-020-1165-7 [PubMed]
• 46. Ayers M, Lunceford J, Nebozhyn M, Murphy E, Loboda A, Kaufman DR, Albright A, Cheng JD, Kang SP, Shankaran V, Piha-Paul SA, Yearley J, Seiwert TY, et al. IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest. 2017; 127:2930–40. https://doi.org/10.1172/JCI91190 [PubMed]
• 47. Pabani A, Butts CA. Current landscape of immunotherapy for the treatment of metastatic non-small-cell lung cancer. Curr Oncol. 2018 (Suppl 1); 25:S94–102. https://doi.org/10.3747/co.25.3750 [PubMed]
• 48. Niemira M, Collin F, Szalkowska A, Bielska A, Chwialkowska K, Reszec J, Niklinski J, Kwasniewski M, Kretowski A. Molecular Signature of Subtypes of Non-Small-Cell Lung Cancer by Large-Scale Transcriptional Profiling: Identification of Key Modules and Genes by Weighted Gene Co-Expression Network Analysis (WGCNA). Cancers (Basel). 2019; 12:37. https://doi.org/10.3390/cancers12010037 [PubMed]
• 49. Yu Y, Zeng D, Ou Q, Liu S, Li A, Chen Y, Lin D, Gao Q, Zhou H, Liao W, Yao H. Association of Survival and Immune-Related Biomarkers With Immunotherapy in Patients With Non-Small Cell Lung Cancer: A Meta-analysis and Individual Patient-Level Analysis. JAMA Netw Open. 2019; 2:e196879. https://doi.org/10.1001/jamanetworkopen.2019.6879 [PubMed]
• 50. Liu D. Cancer biomarkers for targeted therapy. Biomark Res. 2019; 7:25. https://doi.org/10.1186/s40364-019-0178-7 [PubMed]
• 51. Wang X, Wang F, Zhong M, Yarden Y, Fu L. The biomarkers of hyperprogressive disease in PD-1/PD-L1 blockage therapy. Mol Cancer. 2020; 19:81. https://doi.org/10.1186/s12943-020-01200-x [PubMed]
• 52. Sanchez-Vega F, Mina M, Armenia J, Chatila WK, Luna A, La KC, Dimitriadoy S, Liu DL, Kantheti HS, Saghafinia S, Chakravarty D, Daian F, Gao Q, et al, and Cancer Genome Atlas Research Network. Oncogenic Signaling Pathways in The Cancer Genome Atlas. Cell. 2018; 173:321–37.e10. https://doi.org/10.1016/j.cell.2018.03.035 [PubMed]
• 53. Hellmann MD, Ciuleanu TE, Pluzanski A, Lee JS, Otterson GA, Audigier-Valette C, Minenza E, Linardou H, Burgers S, Salman P, Borghaei H, Ramalingam SS, Brahmer J, et al. Nivolumab plus Ipilimumab in Lung Cancer with a High Tumor Mutational Burden. N Engl J Med. 2018; 378:2093–104. https://doi.org/10.1056/NEJMoa1801946 [PubMed]
• 54. Hellmann MD, Callahan MK, Awad MM, Calvo E, Ascierto PA, Atmaca A, Rizvi NA, Hirsch FR, Selvaggi G, Szustakowski JD, Sasson A, Golhar R, Vitazka P, et al. Tumor Mutational Burden and Efficacy of Nivolumab Monotherapy and in Combination with Ipilimumab in Small-Cell Lung Cancer. Cancer Cell. 2019; 35:329. https://doi.org/10.1016/j.ccell.2019.01.011 [PubMed]
• 55. Lee WC, Reuben A, Hu X, McGranahan N, Chen R, Jalali A, Negrao MV, Hubert SM, Tang C, Wu CC, Lucas AS, Roh W, Suda K, et al. Multiomics profiling of primary lung cancers and distant metastases reveals immunosuppression as a common characteristic of tumor cells with metastatic plasticity. Genome Biol. 2020; 21:271. https://doi.org/10.1186/s13059-020-02175-0 [PubMed]
• 56. Davoli T, Uno H, Wooten EC, Elledge SJ. Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science. 2017; 355:eaaf8399. https://doi.org/10.1126/science.aaf8399 [PubMed]
• 57. Jamal-Hanjani M, Wilson GA, McGranahan N, Birkbak NJ, Watkins TBK, Veeriah S, Shafi S, Johnson DH, Mitter R, Rosenthal R, Salm M, Horswell S, Escudero M, et al, and TRACERx Consortium. Tracking the Evolution of Non-Small-Cell Lung Cancer. N Engl J Med. 2017; 376:2109–21. https://doi.org/10.1056/NEJMoa1616288 [PubMed]
• 58. Dyer BA, Feng CH, Eskander R, Sharabi AB, Mell LK, McHale M, Mayadev JS. Current Status of Clinical Trials for Cervical and Uterine Cancer Using Immunotherapy Combined With Radiation. Int J Radiat Oncol Biol Phys. 2021; 109:396–412. https://doi.org/10.1016/j.ijrobp.2020.09.016 [PubMed]
• 59. Cromie MM, Gao W. Epigallocatechin-3-gallate enhances the therapeutic effects of leptomycin B on human lung cancer a549 cells. Oxid Med Cell Longev. 2015; 2015:217304. https://doi.org/10.1155/2015/217304 [PubMed]
• 60. Özdaş S, Özdaş T. Crm1 knockdown by specific small interfering RNA reduces cell proliferation and induces apoptosis in head and neck cancer cell lines. Turk J Biol. 2018; 42:132–43. https://doi.org/10.3906/biy-1711-8 [PubMed]
• 61. Fabi F, Adam P, Vincent K, Demontigny F, Parent S, Joncas FH, Asselin E. Inhibition of CRM1 activity sensitizes endometrial and ovarian cell lines to TRAIL-induced cell death. Cell Commun Signal. 2018; 16:39. https://doi.org/10.1186/s12964-018-0252-z [PubMed]
• 62. Turner JG, Dawson J, Cubitt CL, Baz R, Sullivan DM. Inhibition of CRM1-dependent nuclear export sensitizes malignant cells to cytotoxic and targeted agents. Semin Cancer Biol. 2014; 27:62–73. https://doi.org/10.1016/j.semcancer.2014.03.001 [PubMed]
• 63. Zasadil LM, Andersen KA, Yeum D, Rocque GB, Wilke LG, Tevaarwerk AJ, Raines RT, Burkard ME, Weaver BA. Cytotoxicity of paclitaxel in breast cancer is due to chromosome missegregation on multipolar spindles. Sci Transl Med. 2014; 6:229ra43. https://doi.org/10.1126/scitranslmed.3007965 [PubMed]
• 64. Blair HA, Deeks ED. Albumin-Bound Paclitaxel: A Review in Non-Small Cell Lung Cancer. Drugs. 2015; 75:2017–24. https://doi.org/10.1007/s40265-015-0484-9 [PubMed]
• 65. Plunkett W, Huang P, Xu YZ, Heinemann V, Grunewald R, Gandhi V. Gemcitabine: metabolism, mechanisms of action, and self-potentiation. Semin Oncol. 1995; 22:3–10. [PubMed]
• 66. Hayashi H, Kurata T, Nakagawa K. Gemcitabine: efficacy in the treatment of advanced stage nonsquamous non-small cell lung cancer. Clin Med Insights Oncol. 2011; 5:177–84. https://doi.org/10.4137/CMO.S6252 [PubMed]
• 67. Jordan MA. Mechanism of action of antitumor drugs that interact with microtubules and tubulin. Curr Med Chem Anticancer Agents. 2002; 2:1–17. https://doi.org/10.2174/1568011023354290 [PubMed]
• 68. Asai N, Ohkuni Y, Kaneko N, Yamaguchi E, Kubo A. Relapsed small cell lung cancer: treatment options and latest developments. Ther Adv Med Oncol. 2014; 6:69–82. https://doi.org/10.1177/1758834013517413 [PubMed]