Identification and functional analysis of lipid metabolism-related DEGs
A total of 217 lipid metabolism-related DEGs were identified from the TCGA-LUAD cohort. A volcano plot was constructed to reveal significant DEGs (Fig. 1A), and a heatmap was created to show the hierarchical clustering analysis of the DEGs (Fig. 1B). For the overall understanding of 217 lipid metabolism-related DEGs, GO terms and KEGG pathway enrichment analysis were conducted using the clusterProfiler package, while canonical pathways analysis was performed by IPA. The results of KEGG pathway enrichment showed that DEGs were significantly enriched in arachidonic acid metabolism, metabolism of xenobiotics by cytochrome P450, glycerophospholipid metabolism, and steroid hormone biosynthesis. In contrast, GO terms analysis showed that genes were significantly enriched in the fatty acid metabolic process, glycerolipid metabolic process, fatty acid dericative metabolic process, and steroid metabolic process (Fig. 1C). The genes in each KEGG pathway and GO term are presented in additional file 2. IPA identified significant canonical networks associated with the DEGs. IPA showed that the top canonical pathways associated with common DEGs including eicosanoid signaling, FXR/RXR activation, and atherosclerosis signaling (Fig. 1D). Combining the results of the three functional analyses showed that DEGs mainly overlapped in glycerophospholipid and steroid metabolism. Furthermore, non-overlapping pathways provided additional information indicating further exploration of the role of lipid metabolism in LUAD.
Interaction network construction and cytoHubba analysis
Lipid metabolism-related DEGs were analyzed by the STRING tool. Ultimately, an interaction network with 216 nodes and 1140 edges was established and visualized in Cytoscape (Fig. 2). According to 12 ranked methods in cytoHubba, 6 hub genes were identified by the overlap of the top 10 genes (Additional file 3). Moreover, these genes were related to Insulin (INS), Lipoprotein Lipase (LPL), Hematopoietic Prostaglandin D Synthase (HPGDS), Diacylglycerol O-Acyltransferase 1 (DGAT1), UDP Glucuronosyltransferase Family 1 Member A6 (UGT1A6), and Cytochrome P450 Family 2 Subfamily C Member 9 (CYP2C9).
The expression level analysis of hub genes
DEG results of hub genes are presented in Table 1. The data showed that CYP2C9, UGT1A6, INS, and DGAT1 were upregulated, while HPGDS and LPL were downregulated in TCGA-LUAD tissues compared to normal tissues. To verify the expression results of hub genes, GEPIA and ONCOMINE databases were used. In GEPIA databases, HPGDS and LPL were significantly downregulated in LUAD samples (Fig. S1). In addition, correlation analysis showed that LPL/DGAT1 (r = 0.15; P < 0.01), UGT1A6/HPGDS (r = -0.11; P = 0.02), and HPGDS/DGAT1 (r = -0.09; P < 0.05) were significantly correlated (additional file 4). Meta-analysis of 6 six hub genes of lung cancer was performed by ONCOMINE databases, and showed that UGT1A6 and DGAT1 were upregulated, while HPGDS and LPL were downregulated (Fig. S2).
Table 1. DEG results of hub genes.
|
logFC
|
P
|
FDR
|
CYP2C9
|
1.027946
|
0.001961
|
0.004128
|
UGT1A6
|
3.382161
|
4.80E-31
|
6.74E-30
|
INS
|
1.773607
|
2.42E-07
|
8.55E-07
|
DGAT1
|
1.041643
|
4.93E-14
|
3.24E-13
|
HPGDS
|
-1.19395
|
6.18E-18
|
5.18E-17
|
LPL
|
-1.96376
|
5.62E-83
|
2.01E-81
|
Survival analysis of hub genes
In this study, the relationship between mRNA expression of hub genes and clinical outcome was examined using the Kaplan-Meier plotter. The results showed that high expression of CYP2C9 [HR = 1.50 (1.19–1.90), P < 0.01], UGT1A6 [HR = 1.61 (1.26–2.06), P < 0.01], and INS [HR = 1.46 (1.15–1.85), P < 0.01], and low expression of DGAT1 [HR = 0.78 (0.62–0.98), P = 0.04], HPGDS [HR = 0.58 (0.45–0.73), P < 0.01], and LPL [HR = 0.54 (0.43–0.69), P < 0.01], were associated with a worse OS for 719 LUAD patients (Fig. 3).
Prediction model based on survival-related hub genes and validation
Based on the Cox regression model, a nomogram was built to predict the prognosis of TCGA-LUAD patients, using mRNA expression of the six survival-related hub genes (Fig. 4A). The concordance index of the nomogram was 0.61. Subsequently, the prognosis score of each patient was calculated, which showed that patients in the high-risk score group had a worse OS of 3 years [HR = 1.88 (1.09 - 3.25), P = 0.02] (Fig. 4B). A total of 486 patients with a complete information, including gender, tumor stage, age, and smorking status were included for the multivariate Cox regression analysis. Except for HPGDS and LPL, the HR of CYP2C9, DGAT1, UGT1A6, and INS was not significant. In addition, the risk score calculated from the six-gene signature was an independent prognostic factor (Fig. S3). The model was validated and demonstrated that patients in the high-risk score group had a worse OS [HR = 1.91 (1.02 - 3.50), P = 0.04] (Fig. 4C).