We obtained transcription profiling data and clinical data for LUAD patients from the TCGA database (https://portal.gdc.cancer.gov/) (535 LUAD samples and 59 normal LUAD samples) and the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE11969) (GSE11969 and GPL7015) (94 LUAD samples). The clinical characteristics of the TCGA database are shown in Supplementary Table 1 . The gene IDs of the samples were converted from the human gene annotation file to the corresponding gene symbols and averaged if multiple probes targeted the same gene ID. TCGA-LUAD is the test set, and GEO-LUAD is the train set.

In combination with previous studies (9), we obtained three gene sets (KEGG fatty acid metabolism pathway, Reactome fatty acid metabolism genes, and Hallmark fatty acid metabolism genes) from Molecular Signature Database v7.2 (MSigDB). The screening yielded 309 genes related to fatty acid metabolism ( Supplementary Table 2 ).

First, the “limma” R package was used to perform differential analysis to screen for differential genes related to fatty acid metabolism, and genes with LogFC < 0.585 and FDR < 0.05 were considered statistically significant. The “clusterprofiler” R package was used to enrich the GO and KEGG pathways of differential genes to determine their main biological features and cellular functional pathways. Differences were statistically significant when p-values and corrected p-values <0.05. Finally, the results of the enrichment analysis were visualized using the “ggplot 2” and “goplot” R packages. After removing patients with a survival time of less than 30 days, sequencing data of differentially expressed genes associated with fatty acid metabolism in the samples were combined with survival data, and genes associated with prognosis were screened from those associated with fatty acid metabolism by univariate Cox regression analysis based on the train set, and the cutoff point is set to p-value < 0.05. The correlation between the mutation frequency of genes in the train set samples and the mutated genes was analyzed using the “maftools” R package. The “Glmnet” R package was used for genes associated with the prognosis of fatty acid metabolism in LUAD. Based on the TCGA database, a prognostic risk score model for predicting OS in LUAD samples was developed using Least Absolute Shrinkage Selection Operator (LASSO) Cox regression analysis. The risk score formula was as follows.

The “Coef” represents non-zero regression coefficients calculated using the LASSO Cox regression analysis (Supplementary Table 3), and “ExpGene” is the expression values of genes from the prognostic risk score model.

LUAD patients were divided into high-risk and low-risk groups according to the median fatty acid risk score of the TCGA-LUAD cohort sample, and the K-M method was used to analyze whether the high-risk and low-risk groups differed in terms of survival prognosis. The feasibility of the model was further validated using the GEO database. Nomograms, PCA, and ROC were used to ensure the accuracy of the model.

PCA of gene expression profiles of fatty acid risk scoring models was performed using the “limma” R package for both train and test sets, including the expression profiles of differentially expressed genes associated with fatty acid metabolism in the train set. Fatty acid risk scoring models were constructed and the results were visualized using the “ggplot2” R package.

Clinical information (stage) and fatty acid risk score of patients with LUAD were combined after excluding survival time of less than 30 days and missing data. Independent prognostic indicators were screened using univariate and multivariate Cox regression (p-value < 0.05).

To further investigate the overall survival (OS) of individual LUAD patients, a predictive model based on independent clinical parameters was developed using the “nomogram” R package. ROC curves and calibration plots were used to measure the ability of nomogram to predict prognosis.

Gene mutation data were downloaded from the TCGA database to calculate the tumor mutation burden (TMB) of LUAD patients. The “ggpubr” R package was used to explore whether there was an association between patient risk score and the frequency of tumor mutations in target genes. The immune cell infiltration file was downloaded from Timer2.0 (http://timer.cistrome.org/) to estimate the relationship between immune cell infiltration and risk score (for all TCGA tumors), using the “limma” and “pheatmap” R package for difference analysis, and the results were visualized. We further analyzed whether there were differences in immune-related functions between the high- and low-risk groups using the “GSVA” and “GSEABase” packages. The “PRRophetic” R package was used to predict the semi-inhibitory concentrations of cisplatin, gemcitabine, and paclitaxel in each sample, which indicates the effectiveness of a substance to inhibit a specific biological or biochemical function. The TIDE online database (http://tide.dfci.harvard.edu/) was used to predict the effect of immunotherapy in the high-risk and low-risk groups, and differences were considered statistically significant when the p-value < 0.05.

After screening for differential genes between high- and low-risk scoring groups, PPI network data [generated interaction score >0.90 (medium confidence)] were plotted online using the STRING online database (https://cn.string-db.org/). The 10 most pooled hub genes were screened using the cubHubba plugin of Cytoscape software (version: 3.9.1). Based on these 10 hub genes, the target genes differentially expressed in tumor tissues and normal tissues (cutoff value log2FC > 1, p-value < 0.05) and associated with postnatal survival were screened by the GEPIA (http://gepia.cancer-pku.cn) online database. The CIBERSORTx (https://cibersortx.stanford.edu/index.php online database was used to analyze the infiltration of 22 tumor-infiltrating lymphocyte-associated target genes in the microenvironment of high-risk and low-risk LUAD patients for significance ranking analysis 1,000 times, and “reshape2” and “ggpubr” R packages were used to visualize the difference results. The “limma” and “ggpubr” R packages were used to analyze the relationship between target gene expression and clinical characteristics (stage, T, N, M, age, and gender) of LUAD patients. Finally, the correlations between target genes were evaluated based on the GEPIA (http://gepia.cancer-pku.cn) online database and the Spearman test.

By comparing the expression levels of fatty acid metabolism-related genes in tumor and normal samples in the TCGA database, 126 genes were screened in the TCGA-LUAD cohort (p-value < 0.05, FDR < 0.585), of which 79 genes were upregulated in tumor tissue samples and 47 genes were downregulated in tumor tissue samples. Heat maps and volcano maps of differentially expressed genes in normal and tumor samples are shown in Figures 1A, B . We know from GO enrichment analysis that among the biological processes, fatty acid metabolic processes, long-chain fatty acid metabolic processes, and fatty acid biosynthesis processes are highly enriched terms ( Figure 1C ), and Figure 1D shows the 72 most significantly enriched genes and enriched pathways. The results of KEGG enrichment analysis of genes and genomes in Kyoto showed that Fatty acid metabolism, Fatty acid degradation, Fatty acid elongation, and Fatty acid biosynthesis were all highly enriched KEGG items ( Figure 1E ), and Figure 1F shows the 55 most significantly enriched genes and enriched pathways. These results suggest that fatty acid metabolism-related genes in LUAD are clustered in biological pathways related to fatty acid metabolism and are closely related to fatty acid metabolism.

The TCGA-LUAD cohort was used as a train set. A total of 126 cases of fatty acid metabolism-related genes were screened from the TCGA-LUAD cohort. After excluding data from patients with a survival time of less than 30 days and null values, 29 fatty acid metabolism-related genes associated with patient survival were screened using univariate Cox analysis ( Figure 2A ). The somatic mutation profiles of the 29 prognosis-related genes showed a mutation frequency of 17.29% in the 561 LUAD samples (a total of 97 cases were mutated) ( Figure 2B ), with ADH1B and CYP4B1 having the highest mutation frequencies, followed by ELOVL6, MDH2, SMS, LTA4H, ENO3, ALDOA, ALOX15, ALOX15B, MAOA, LDHA, CA4, ELOVL2, DPEP2, CEL, HSD17B4, INMT, AOC3, and ACOXL, and no mutations in other genes. Further analysis revealed a mutational positive relationship between ACOXL and HSD17B4, ADH1B; CYP4B1 and AOC3, and PTGR1 and LDHA ( Figure 2C ). These 29 genes were further incorporated into the LASSO logistic regression algorithm based on the TCGA-LUAD cohort. A total of 14 genes for constructing fatty acid risk score models were obtained by LASSO Cox regression analysis ( Figures 2D, E ), namely, ALDH2, HACD1, ELOVL2, ENO3, CEL, CA4, CYP2U1, LDHA, ALOX5AP, SMS, ALDOA, CYP4B1, DPEP2, and ELOVL6.

Using this risk score model, the LUAD samples (low and high risk) were completely distinguished ( Figures 2F, G ), while the GEO-LUAD cohort was similarly distinguished as a test set ( Figure 2H ). In this fatty acid risk model, the median risk score of the TCGA-LUAD cohort was used as the cutoff value, and 490 TCGA-LUAD patients were classified into a high-risk group (n = 245) and a low-risk group (n = 245), and 94 GEO-LUAD patients were classified into a high-risk group (n = 51) and a low-risk group (n = 43). In both the test and train groups, the low-risk group had a better clinical prognosis (p-value < 0.05) ( Figures 2I, J ). Univariate prognostic COX analyses of the train and test sets showed that stage and risk score were independent prognostic factors, and multivariate COX analyses also showed this result (p-value < 0.05) ( Figures 2K–N ). Using the ROC curve estimation model to evaluate the reliability of the fatty acid risk score, the area under the curve (AUC) was 0.727, 0.715, and 0.677 for the train sets at 1, 3, and 5 years, respectively ( Figure 2O ). The AUC was 0.802, 0.656, and 0.653 for the test sets at 1, 3, and 5 years, respectively ( Figure 2Q ); the train set risk model has the largest area under the ROC curve ( Figure 2P ) and the test set risk model. The ROC curve was second only to stage, which proved the most reliable model of fatty acid metabolism risk ( Figure 2R ).

Prognostic factors such as risk score, stage, age, and sex were included in the nomogram to assess the predictive ability of individual OS, and to calculate the prediction effect of the nomogram on 1-, 3-, and 5-year OS in patients with LUAD ( Figure 3A ). The nomogram was applied to patients in the test set (GEO database), the effect was good ( Figure 3C ), and the calibration charts were close to the ideal curve ( Figures 3B, D ). In addition, ROC curve analysis is carried out to verify the practicability of the nomogram, and the AUC of train set and test set risk, nomogram, age, gender, and the stage is calculated ( Figures 3E,F ). Nomogram AUC is the largest, indicating that the nomogram has the best prediction effect. Further univariate and multivariate Cox analysis of the train set data showed that the classification and risk score were independent prognostic factors (p-value < 0.001) ( Figures 3G, H ). In conclusion, the prognostic ability of this fatty acid prognostic model has been validated from several perspectives.

Immune cell infiltration showed immunological differences between the high-risk and low-risk groups, with significantly increased abundance of Macrophage M0, Macrophage M1, T cell CD4+ Th1, T cell CD4+ Th2, and T cell CD4+ memory activated in the high-risk group and B cell memory, Macrophage M2, T cell CD4+ central memory, T cell CD4+ effector memory, T cell CD8+, Myeloid dendritic cell resting, Myeloid dendritic cell activated, and T cell regulatory significantly increased in the low-risk group ( Figure 4A ). Immune function analysis ( Figure 4B ) showed that HLA and Type_II_IFN_Response immune-related functions were active in the low-risk group and MHC_class_I immune-related functions were active in the high-risk group. TIDE scores were lower in the high-risk group than in the low-risk group ( Figure 4C ), which proved that the high-risk group had a better effect of immunotherapy. In the risk score and chemical drug sensitivity analysis, we found that the risk score was negatively correlated with cisplatin, gemcitabine, and paclitaxel chemotherapy drug sensitivity. The IC50 of patients in the high-risk group with the lower value ( Figures 4D–I ) proves that patients in the high-risk group were more sensitive to chemotherapy drugs than those in the low-risk group. In conclusion, the quantification of fatty acid metabolic risk scores is of great importance in patients with LUAD, not only to assess the prognosis of immunotherapy, but also to evaluate the effect of chemotherapy, which may be a new biomarker to change the outcome of treatment in patients with LUAD.

Differential gene interactions between low- and high-risk groups were analyzed using the STRING online database. The differential gene PPI network is shown in Figure 5A . Hub genes were identified from differential genes using the Cytoscape plugin cytoHubba. We chose a total of 10 genes in the network ( Figure 5B ) and ranked CDK1, BUB1, CCNA2, CCNB1, CDC20, BUB1B, CCNB2, DLGAP5, TPX2, and TTK utilizing the degree method Figure 5C ). Through differential analysis and survival analysis, three hub genes were screened as target genes. We performed a differential analysis based on TCGA and GTEx databases using GEPIA and found that BUB1B ( Figure 5D ), CCNB1 ( Figure 5E ), and TTK ( Figure 5F ) were significantly overexpressed in LUAD tissues compared with normal samples (p-value < 0.05). BUB1B ( Figures 5G, H ), CCNB1 ( Figures 5I, J ), and TTK ( Figures 5K, L ) were significantly associated with survival prognosis (OS and RFS) (p-value < 0.01), and the higher the gene expression, the worse the prognosis.

The results of immune cell infiltration showed that the upregulated immune cells in the gene BUB1B high expression group were T cells CD8, T cells CD4 memory activated, NK cells resting, Macrophages M0, and Macrophages M1, and the downregulated immune cells were B cells naive, T cells CD4 memory resting, T cells regulatory (Tregs), NK cells activated, Monocytes, Dendritic cells resting, and Mast cells resting ( Figure 5M ); the upregulated immune cells in the CCNB1 high expression group were T cells CD8, T cells CD4 memory activated, T cells follicular helper, NK cells resting, Macrophages M0, and Macrophages M1, and the downregulated immune cells were B cells memory, T cells CD4 memory resting, Monocytes, Dendritic cells resting, and Mast cells resting ( Figure 5N ); the upregulated immune cells in the TTK high expression group were T cells CD8, T cells CD4 memory activated, Macrophages M0, Macrophages M1, and Dendritic cells activated, and downregulated immune cells were B cells memory, Plasma cells, T cells CD4 memory resting, T cells regulatory (Tregs), NK cells activated, Monocytes, Dendritic cells resting, and Mast cells resting ( Figure 5O ). We noted that the infiltration of activated immune cells was abundant and active in the high expression group, which was suitable for immunotherapy. Clinical correlation analysis showed that in the BUB1B, CCNB1, and TTK genes, stage I patients had significantly lower expression levels than stage II, III, and IV patients, and higher gene expression levels in men and patients younger than 65 years (p-value < 0.05) ( Figures 5P–R ). There was no significant difference between BUB1B, CCNB1 mutant, and wild-type risk scores ( Figures 5S, T ), and there was a significant difference between TTK mutant and wild-type risk scores (p-value < 0.05) ( Figure 5U ). Based on the GEPIA (http://gepia.cancer-pku.cn) online database, correlations between genes were analyzed using the Spearman test. The three genes BUB1B, CCNB1, and TTK were positively correlated with each other with statistically significant differences ( Figures 5V–X ) (p-value < 0.001).