Prediction of prognosis and immunotherapy efficacy based on metabolic landscape in lung adenocarcinoma by bulk, single-cell RNA sequencing and Mendelian randomization analyses

Heterogenous metabolic profiles of LUAD

The overall design of our study was shown in the flow chart (Figure 1). To investigate the metabolic reprogramming in LUAD, we extracted 84 metabolic pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The 84 metabolic pathways contained 12 kinds of metabolism, including 14 carbohydrate metabolic pathways, 13 lipid metabolic pathways, 12 cofactors and vitamins metabolic pathway, 13 amino acid metabolic pathways, 14 glycan biosynthesis and metabolic pathways, 2 biosynthesis pathways of other secondary metabolites, 3 energy metabolic pathways, 1 genetic information processing pathway, 6 other amino acids metabolic pathways, 1 terpenoids and polyketides metabolic pathway, 2 nucleotide metabolic pathway and 3 xenobiotics biodegradation and metabolic pathway. We first scored each pathway in all samples using ssGSEA and characterized the metabolic heterogeneity between LUAD and normal tissues (Figure 2A). It was demonstrated that a total of 68 (80.95%) metabolic pathways were dramatically dysregulated (23 upregulated pathways and 45 downregulated pathways in LUAD compared to normal lung tissues) (P < 0.05). The dysfunctional pathways encompassed the three main kinds of metabolism including carbohydrate, lipid and amino acids metabolism. Next, to investigate the intratumor heterogeneity of metabolism, we classified TCGA-LUAD samples into two clusters based on 84 metabolic pathways scores by unsupervised consensus clustering (Figure 2B). PCA and metabolic heatmap displayed that the two clusters had obviously heterogenous metabolic characteristics (Figure 2C, 2D). The cluster A seemed to be a “cold” metabolic subtype but the cluster B seemed to a “hot” metabolic subtype. To further distinguish the metabolic variation in various TME cells, we compared metabolism scores among endothelial cells, fibroblast, malignant cells and pan-immune cells by scRNA-seq data. The results suggested that malignant cells had a significantly higher metabolic level than the other three cells (Figure 2E). The above results indicated that metabolic heterogeneity was common in LUAD and might play a crucial role in the initiation and progression of LUAD.

Figure 1. The flow chart of the study.

Figure 2. The metabolic heterogeneity of lung adenocarcinoma (LUAD). (A) The differences of 84 metabolic pathways scores between LUAD and normal tissues in TCGA. (B) An unsupervised consensus clustering according to 84 metabolic pathways scores in TCGA-LUAD samples. (C) Principal Component Analysis of cluster A and B of TCGA-LUAD. (D) The differences of 84 metabolic pathways scores between cluster A and cluster B. Pathological stage, sex, age, and survival status used as patients’ annotation. (E) The differences of 84 metabolic pathways scores among different cells by single-cell RNA sequencing (scRNA-seq) data. *P < 0.05, **P < 0.01, ***P < 0.001.

Development of a prognostic model based on MPPS and exploration of its clinical relevance

Different sequencing platforms commonly possess different sequencing depth. Consequently, sequencing data from different sequencing platforms had different numbers of genes and significantly different expression levels. When gene expression data from different platforms were utilized, data standardization and scaling and intersecting gene from different platforms are needed, which will cause loss of some genetic information. The metabolic pathway pair model can reduce the effects of some gene deletion on prediction and eliminate the shortcomings of data standardization and scaling in gene expression data processing and effectively avoid the interference caused by the sequencing platform. To develop a prognostic model based on MPPS, we firstly paired the 84 pathways and 525 pathway pairs were obtained after removing those pathway pairs that proportion of 0 or 1 was more than 80% or less than 20%. Subsequently, we conducted univariable Cox regression analysis on these pathway pairs and selected 106 meaningful metabolic pathway pairs for Least absolute shrinkage and selection operator (LASSO) regression analysis (Figure 3A). LASSO regression analysis yielded 33 metabolic pathway pairs with nonzero LASSO coefficients according to the optimal λ value (Figure 3B, 3C). Multivariable Cox regression analysis was further performed to identify prognostic metabolic pathway pairs based on Akaike information criterion value and 19 metabolic pathway pairs were finally obtained (Figure 3D). MPPS was calculated using value of 19 metabolic pathway pairs weighted by their multivariable Cox regression coefficients and stratified LUAD patients into high- and low-MPPS groups according to the optimal cut-off point determined by the “survminer” package. PCA analysis showed that LUAD patients could be divided into distinctive groups according to MPPS (Figure 3E). The heatmap showed obvious discrepancy of 19 metabolic pathway pairs between the high- and low-MPPS group (Supplementary Figure 1A). Patients with high-MPPS scores had significantly shortened overall survival (OS) and progression-free survival (PFS) in the TCGA-LUAD training cohort and six GEO validation cohorts (all P < 0.05). The GEO merge cohort integrating the six GEO cohorts also showed the same trend (P < 0.05) (Figure 3F). The risk plot of MPPS indicated that as MPPS increased, OS time decreased while mortality rose (Supplementary Figure 1B).

Figure 3. Construction of a metabolic pathway-pair score (MPPS)-based prognostic model. (A) The univariable Cox regression analysis of metabolic pathway pairs. (B, C) Determination of the number of metabolic pathway pairs by the LASSO regression analysis. (D) The multivariable Cox regression analysis of metabolic pathway pairs. (E) Principal Component Analysis of the high- and low-risk groups. (F) The Kaplan-Meier analysis of the high- and low-risk groups in TCGA-LUAD cohort and GEO validation cohorts. OS: overall survival; PFS: progression-free survival. The survival analysis was tested by log-rank test.

To determine the correlation of MPPS and clinical traits, we compared the differences in MPPS among different clinical subgroups based on age, sex, survival status and pathological stage. Patients in alive, stage I, stage T1 and stage N0 subgroups had lower MPPS compared to the other subgroups (P < 0.05), while there was no significant difference of MPPS in age, sex and M stage subgroups (Supplementary Figure 2A–2G). The Sankey diagram illustrated the distribution and correspondence of LUAD patients in MPPS groups, survival status, age, sex, and pathological stage (Supplementary Figure 2H). In addition, MPPS also showed robust performance on predicting prognosis in different clinical subgroups, including age, sex, TNM stage (P < 0.05) (Supplementary Figure 2I–2P).

Evaluation of the MPPS model

To investigate the accuracy of the MPPS model, ROC analysis was conducted and showed good performance in both training cohort and validation cohort (1-, 3-, 5-year AUC: 0.755, 0.781, 0.785 in OS of TCGA-LUAD; 0.664, 0.689, 0.671 in PFS of TCGA-LUAD; 0.599, 0.658, 0.695 in OS of GSE3141; 0.832, 0.705, 0.682 in OS of GSE13213; 0.699, 0.641, 0.633 in OS of GSE30219; 0.789, 0.661, 0.678 in OS of GSE31210; 0.671, 0.624, 0.693 in OS of GSE50081; 0.641, 0.663, 0.684 in OS of GSE72094; 0.672, 0.633, 0.66 in OS of GEO merge cohort) (Figure 4A). With the developments in next-generation sequencing, a considerable number of prognostic models were developed. To compare the performance of MPSS with other signatures, we retrieved 51 published signatures of LUAD including 15 lncRNA signatures and 36 mRNA signatures. These signatures encompassed various biological processes, such as autophagy, immune response, ferroptosis, stemness, epithelial-mesenchymal transition (EMT), hypoxia, ageing, methylation et al. The results suggested that MPSS had highest 1-, 3-, 5-year AUC of TCGA-LUAD cohort (Figure 4B). MPPS also displayed higher C-index in TCGA-LUAD cohort than almost all models (Figure 4C).

Figure 4. Evaluation of MPPS model. (A) The 1-, 3-, 5-year receiver operating characteristic (ROC) curves of MPPS model in TCGA-LUAD cohort and GEO validation cohorts. (B) The 1-, 3-, 5-year area under curves (AUC) of MPPS model and 51 published prognostic models of LUAD. (C) C-index of MPPS model and 51 published prognostic models of LUAD.

To further investigate the performance of MPPS on predicting prognosis of other tumors, we performed survival analyses of patients in the high- and low-MPPS groups involving 32 types of tumors in TCGA other than LUAD. Patients in the high-MPPS group had significantly worse OS than low-MPPS group in all 32 tumors (P < 0.05) (Figure 5A). MPPS also displayed high 1-, 3-, 5-year AUC in predicting prognosis of 32 tumors (Figure 5B–5D).

Figure 5. Evaluation of MPPS model in pan-cancer. (A) The Kaplan-Meier analysis of the high- and low-risk groups across 32 tumors in TCGA database except LUAD. The survival analysis was tested by log-rank test. (B–D) The 1-, 3-, 5-year AUC of MPPS model across 32 tumors in TCGA database except LUAD.

Pathway enrichment and function annotation of the high- and low-MPPS groups

To explore the underlying mechanism of survival variation in different MPPS groups, we first performed pathway enrichment and function annotation of the high- and low-MPPS groups by KEGG and Gene Ontology (GO) analyses. The high-MPPS group had higher metabolic levels in pentose phosphate pathway, pyrimidine, cysteine and methionine metabolism. Multiple proliferation-related pathways including DNA replication, cell cycle, mismatch repair were enriched in the high-MPPS group (P < 0.05). On the contrary, immune response pathways such as B cell receptor signaling pathway, JAK/STAT signaling pathway, T cell receptor signaling pathway and cytokine-cytokine receptor interaction were more enriched in the low-MPPS group (Figure 6A). GO function annotation also demonstrated that DNA replication and translation were more associated with high MPPS, and immune cell development, maturation, activation and response were more related to low MPPS (P < 0.05) (Figure 6B). Additionally, multiple oncogenic pathways including hypoxia, epithelial-mesenchymal transition, DNA damage response, glycolysis, unfolded protein response and mTOC1 signaling et al. were significantly enriched in the high-MPPS group (Supplementary Figure 3A). The above results all indicated that LUAD with high MPPS possessed higher malignancy.

Figure 6. Pathway enrichment and function annotation of the high- and low-risk groups. (A) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses between the high- and low-risk groups. (B) Gene Ontology (GO) between the high- and low-risk groups. (C) The differences of 31 metabolic pathways scores between the high- and low-risk groups. (D) The bubble diagram was drawn by the average of 19 metabolic pathway pairs in pan-immune cells, endothelial cells, fibroblasts, and malignant cells. (E) The differences of genes in 19 metabolic pathway pairs between the high- and low-risk groups. Pathological stage, sex, age, and survival status used as patients’ annotation. (F) The protein-protein interaction network by STRING website and the software Cytoscape v3.9.1. *P < 0.05, **P < 0.01, ***P < 0.001.

To characterize metabolic reprogramming involved in MPPS, we firstly analyzed the variation of 31 pathways in MPPS model between LUAD and normal lung tissues. A total of 26 pathways were dysregulated in LUAD, in which 7 metabolic pathways were up-regulated and 19 metabolic pathways were down-regulated (Supplementary Figure 3B). Further analyses revealed that 22 pathways were significantly variant between the high- and low-MPPS groups. The high-MPPS group had higher metabolic scores in the galactose metabolism, fatty acid (FA) elongation, pyrimidine metabolism, cysteine and methionine metabolism, one carbon pool by folate and aminoacyl-tRNA biosynthesis et al. indicating that biomass synthesis and proliferation were more active in the high-MPPS group. The low-MPPS group had higher metabolic scores in the caffeine metabolism, valine, leucine and isoleucine degradation, selenocompound metabolism and arachidonic acid metabolism, etc. (Figure 6C). The correlation of MPPS and 31 metabolic pathways was shown in Supplementary Figure 3C. By calculating the averages of 19 metabolic pathways pairs in each kind of TME cell, we found the lung malignant cells had higher level of cysteine and methionine metabolism/ganglio series biosynthesis than the immune cells, endothelial cells and fibroblasts (Figure 6D). The further differential expression analysis displayed that there existed obviously differential expression in 31 metabolic pathway genes between the high- and low-MPPS groups and the protein-protein interaction analysis demonstrated that the 76 DEGs had complex regulatory network (Figure 6E, 6F).

Evaluation of TME and immunotherapeutic benefits in the high- and low-MPPS groups

Considering many immune-related pathways were enriched in the low-MPPS group, we evaluated the TME components between the high- and low-MPPS groups by estimate algorithm. The low-MPPS group had higher stromal score, immune score and ESTIMATE score than high-MPPS group (Figure 7A). The GSVA enrichment analysis was performed to evaluate immune filtrating cells and immunologic functions. The results showed that low-MPPS group had higher immune cells infiltration and immunologic functions activation in total including the infiltration of activated B cells, activated CD8+ T cells, activated dendritic cells, eosinophil and macrophage and the immune checkpoint, HLA, T cell co-inhibition or stimulation, type II IFN response (Figure 7B, 7C). A total of 29 immune checkpoint genes were differentially expressed between the high- and low-MPPS groups, in which 27 immune checkpoint genes, accounting for 93.1% were highly expressed in the low-MPPS group (Figure 7D). In addition, a total of 16 (66.7%) HLA genes expression were altered and they were all upregulated in the low-MPPS group (Figure 7E). The cancer–immunity cycle elucidates antitumor immune responses and offers an opportunity to understand the interactions between cancer and its immune system [23]. The low-MPPS group had higher cancer–immunity cycle scores in cancer antigens presentation, priming and activation, CD4⁺ T cell, dendritic cell, B cell, Th17 cell recruiting, and immune cells tumor infiltration but lower scores in cancer antigens release and eosinophil recruiting than the high-MPPS group (Figure 7F). To further investigate the correlation between MPPS and immunotherapy efficacy, we calculated the TIDE score. The results suggested that the low-MPPS group had higher T cell dysfunction score than the high-MPPS group (Supplementary Figure 4A). LUAD with high-MPPS score was more inclined to immune-desert or excluded phenotype and LUAD with low-MPPS score was more inclined to immune-inflamed phenotype (Supplementary Figure 4B). Higher immunophenoscore (IPS) was also exhibited by patients in the low-MPPS group compared with those in the high-MPPS group (Supplementary Figure 4C). The above results indicated that patients in the low-MPPS group may be more sensitive to immunotherapy.

Figure 7. TME landscapes of the high- and low-risk groups. (A) The differences of stromal score, immune score and ESTIMATE score between the high- and low-risk groups. (B) The differences of infiltrating immune cells between the high- and low-risk groups by single sample gene set enrichment analysis (ssGSEA). (C) The differences of immune function between the high- and low-risk groups by ssGSEA. (D) The differences of immune checkpoint genes expression between the high- and low-risk groups. (E) The differences of HLA-related genes expression between the high- and low-risk groups. (F) The differences of cancer-immunity cycle scores between the high- and low-risk groups. *P < 0.05, **P < 0.01, ***P < 0.001.

To further validate the speculation, seven independent immunotherapy cohorts in the published literatures were used to validate immunotherapy efficacy and prognosis including advanced urothelial cancer treated with atezolizumab, an anti-PD-L1 antibody, melanoma treated with anti-CTLA4 and anti-PD-1 therapy, metastatic melanoma treated with anti-CTLA4 therapy, non-small cell lung cancer (NSCLC) treated with nivolumab or pembrolizumab, an anti-PD-1 antibody, NSCLC treated with anti-PD-1/PD-L1 antibody, melanoma treated with ACT, Melanoma treated with anti-PD-1 antibody. The results showed that the low-MPPS group had significant survival advantage and higher immune response rate compared to the high-MPPS group in all validation cohorts (Figure 8A–8N). The response to anti-PD-1 and anti-CTLA4 therapy was calculated using the TIDE website based on TCGA cohort. Patients in the low-MPPS group were more likely to be responders and benefit from immunotherapy (Figure 8O, 8P).

Figure 8. Prediction of immunotherapy by MPPS model. Survival analysis (A) and response to anti-PD-L1 therapy (B) between the high- and low-risk groups in advanced urothelial cancer (IMvigor210 cohort). Survival analysis (C) and response to anti-CTLA4 and anti-PD1 therapy (D) between the high- and low-risk groups in melanoma (GSE91061). Survival analysis (E) and response to anti-CTLA4 therapy (F) between the high- and low-risk groups in metastatic melanoma. Survival analysis (G) and response to anti-PD1 therapy (H) between the high- and low-risk groups in NSCLC (GSE126044). Survival analysis (I) and response to anti-PD-1/PD-L1 therapy (J) between the high- and low-risk groups in NSCLC (GSE135222). Survival analysis (K) and response to adoptive T cell therapy (L) between the high- and low-risk groups in melanoma. Survival analysis (M) and response to anti-PD-1 therapy (N) between the high- and low-risk groups in melanoma (GSE78220). (O) Difference of responder between low- and high-risk group of LUAD in TCGA. (P) Difference of benefits between low- and high-risk group of LUAD in TCGA.

TMB and drug sensitivity analysis

To explore the correlation of MPPS and tumor mutation, Spearman correlation analysis was performed and significant positive correlation was found between MPPS and TMB (R=0.17, P=0.00013) (Figure 9A). LUAD patients in the high-MPPS group had higher TMB than those in the low-MPPS group (P=0.0038) (Figure 9B). By Kaplan-Meier analysis, we found LUAD patients with low-MPPS score and high TMB had the best survival advantages and LUAD patients with high-MPPS score and low TMB had the worst prognosis (P<0.001) (Figure 9C). The distribution of somatic mutations in the high- and low-MPPS groups was investigated in the TCGA-LUAD cohort. Patients in the high-MPPS group displayed significantly higher frequencies of somatic mutations compared with those in low-MPPS group (93.33% vs 85.81%), especially in TP53 (52% vs 39%), TTN (51% vs 35%), MUC16 (43% vs 36%), RYR2 (37% vs 32%), CSMD3 (41% vs 29%) and LRP1B (34% vs 25%). Moreover, missense mutation and multi-hit were the main mutation type in both high- and low-MPPS groups (Figure 9D, 9E). To further explore the clinical utility of MPPS in precision medicine, we assessed the sensitivity of 137 chemotherapeutic or targeted therapy drugs in different MPPS groups (Figure 9F). The results showed that the patients in the high-MPPS groups had lower IC50 values of 67 drugs, indicating sensitivity. Alternatively, the patients in the low-MPPS group were sensitive to 26 drugs. Together, the results may provide a standard of reference for treatment stratification of patients with LUAD.

Figure 9. Tumor mutation burden (TMB) and drug sensitivity analysis. (A) The correlation of MPPS and TMB in TCGA-LUAD samples. (B) The differences of TMB between the high- and low-risk groups. (C) The Kaplan–Meier curves show OS differences stratified by TMB and MPPS. Visualization of the top 20 gene mutations in high-risk group (D) and low-risk group (E). (F) The sensitivity of 117 drugs between the high- and low-risk groups. P < 0.05 was considered as statistical significance.

Establishment of a MPPS-based nomogram

To construct a MPPS-based nomogram for convenient use, we analyzed the prognostic value of MPPS, age, sex, pathological stage and treatment type by univariable and multivariable Cox regression analyses (Table 1). In the univariable Cox regression analysis, age (HR (95%CI): 1.43 (1.056-1.936), P=0.021), stage (HR (95%CI): stage II 2.302 (1.590-3.332), P<0.001; stage III/IV 3.295 (2.309-4.702), P<0.001) and MPPS (HR (95%CI): 1.042 (1.036-1.049), P<0.001) were significantly related to LUAD prognosis. After adjusted by multivariable Cox regression analysis, MPPS, age and pathological stage were identified as independent prognostic factors and used to construct a prognostic nomogram (HR (95%CI): age 1.45 (1.065-1.975), P = 0.018; stage II 2.371 (1.637-3.434), P < 0.001; stage III/IV 2.537 (1.756-3.667), P < 0.001; MPPS 1.039 (1.032-1.046), P < 0.001). The prognostic nomogram made quantitative predictions of the 1-, 3-, and 5-year OS probabilities in patients with LUAD (Figure 10A). The calibration curves exhibited a high consistency between the predicted and actual 1-, 3-, and 5-year OS (Figure 10B). The ROC curves displayed the nomogram had higher AUC values than the single predictor such as MPPS, age, sex, stage, treatment type (1-, 3-, 5-year AUC: 0.793, 0.821, 0.82) (Figure 10C–10E).

Figure 10. Development and evaluation of a prognostic nomogram. (A) Nomogram composed of MPPS, age, stage to predict 1-, 3-, 5-year OS probability. (B) Calibration curves of 1-, 3-, 5-year OS by nomogram. 1- (C), 3- (D), 5- (E) year ROC curves of MPPS, nomogram, age, sex, stage, and treatment type.

Identification of a MPPS-related gene signature by WGCNA and machine learning

To identify MPPS-related modules, WGCNA analysis was performed and 21 modules were identified. 231 genes with gene significance (GS)>0.25, module membership (MM)>0.2 and P < 0.05 were considered as hub MPPS-related genes. Therefore, the hub genes in cyan, tan and turquoise modules met the criterion (Figure 11A). Intersecting with GEO genes from 6 GEO cohorts and DEGs between TCGA-LUAD and normal tissues, 104 hub MPPS-related genes were identified for subsequent analysis (Figure 11B). Based on the expression profiles of 104 hub MPPS-related genes, univariable Cox analysis identified 82 prognostic genes. These 82 genes were subjected to our machine learning-based integrative procedure to develop a consensus MPPS-related gene signature. In the TCGA-LUAD dataset, we fitted 117 kinds of prediction models via the LOOCV framework and further calculated the C-index of each model across 6 GEO validation datasets (Figure 11C). Finally, the 7-gene signature composed of ECT2, ANLN, SLC2A1, LDHA, GAPDH, C1QTNF6 and KRT8 identified by a combination of Lasso regression and survival-SVM had the highest mean C-index in the 6 validation cohorts (Figure 11D, 11E). A gene-based risk score for each patient was calculated by the survival-SVM algorithm and divided patients into the high- and low-risk group according to the optimal cutoff value determined by the “survminer” package. To validate the prognostic value of the gene signature, we performed Kaplan-Meier analysis. The patients in the high-risk group had significantly dismal OS and PFS compared to the low-risk group in the TCGA-LUAD training cohort and six GEO validation cohorts (all P < 0.05) (Figure 11F). The GEO merge cohort also showed the same trend (P < 0.05). In addition, ROC analysis measured the discrimination of the gene signature, with 1-, 3-, 5-year AUCs of 0.697, 0.704, 0.626 in OS of TCGA-LUAD; 0.643, 0.615, 0.556 in PFS of TCGA-LUAD; 0.763, 0.771, 0.686 in OS of GSE3141; 0.861, 0.676, 0.692 in OS of GSE13213; 0.673, 0.752, 0.777 in OS of GSE30219; 0.777, 0.727, 0.757 in OS of GSE31210; 0.758, 0.718, 0.698 in OS of GSE50081; 0.684, 0.643, 0.662 in OS of GSE72094; 0.716, 0.689, 0.699 in OS of GEO merge cohort (Supplementary Figure 5A). To further verify the predicting performance of the gene signature in the clinical practice, we next evaluated the mRNA expression of the 7 genes in a clinical cohort of 42 LUAD patients by qRT-PCR. The Kaplan-Meier analysis showed the low-risk group had better prognosis than the high-risk group (Supplementary Figure 5B). The model had high accuracy in predicting OS with 1-, 3-, 5-year AUCs of 0.763, 0.725, 0.762 in the clinical practice (Supplementary Figure 5C). Together, the MPPS-related gene signature had robust performance in prognostic prediction of LUAD.

Figure 11. Identification of MPPS-related genes by WGCNA and machine learning. (A) Correlation analysis between module eigengenes and MPPS by WGCNA. (B) The intersection of WGCNA hub genes, TCGA-DEGs, and GEO genes. (C) The C-index of 117 machine learning algorithm combinations via LOOCV framework across all validation datasets. (D, E) Determination of the number of MPPS-related genes by the LASSO regression analysis. (F) The Kaplan-Meier analysis of the high- and low-gene risk scores groups stratified by LASSO and survival-SVM in the training and validation cohorts.

Expression, function, prognosis analyses of MPPS-related gene signature

To investigate the correlation of MPPS and the gene risk score, Spearman correlation analysis was performed and significantly positive correlation was observed with R=0.6 and P < 2.2e-16 (Figure 12A). Next, the 7 genes all exhibited obviously higher expression in LUAD than normal lung tissue (Figure 12B). Compared to the low-MPPS group, the high-MPPS group had significantly upregulated expression (Figure 12C). To explore the correlation of 7 genes and 31 metabolic pathways, the correlation heatmap was drawn (Figure 12D). The result showed that there was a positive correlation among the 7 genes expression and their expression was positively associated with FA elongation, pyrimidine metabolism, cysteine and methionine metabolism and one carbon pool by folate and negatively associated with valine, leucine and isoleucine degradation, selenocompound metabolism, glycerophospholipid metabolism and arachidonic acid metabolism. Subsequently, we analyzed the correlation between the gene signature and 14 functional states across 18 cancers using CancerSEA data. The results manifested that the signature was positively related to LUAD proliferation, invasion, cell cycle, DNA damage and repair (Figure 12E). To investigate the mechanism underlying dysregulated expression, the CNV analysis was applied. The ECT2, SLC2A1, KRT8, ANLN and C1QTNF6 showed widespread CNV amplification. In contrast, GAPDH and LDHA had prevalent CNV depletion. The locations of CNV alterations of the 7 MPPS-related genes on chromosomes are shown in Figure 12F. Finally, the prognostic value of the 7 MPPS-related genes was analyzed by Kaplan-Meier curve in the TCGA-LUAD cohort (Figure 12G). The upregulation of the 7 genes all indicated worse survival (P<0.001).

Figure 12. Expression, function, prognosis analyses of MPPS-related gene signature. (A) Correlation of MPPS and gene risk score by LASSO and survival-SVM. (B) The differences of the expression of 7 genes between TCGA-LUAD and normal samples. (C) The differences of the expression of 7 genes between the high- and low-MPPS groups. (D) The correlation of 7 genes expression and metabolic pathways. (E) The correlation of 7 genes expression and 14 biological processes by cancerSEA website. (F) The copy number variation frequency and location in chromosomes of the 7 genes. (G) The Kaplan-Meier analysis of the high- and low-expression of the 7 genes.

Effect of C1QTNF6 on infiltrating immune cells of TME

By investigating MPPS-related genes expression in specific cells in TME, we found only C1QTNF6 was highly expressed in fibroblasts and the other six genes were mainly expressed in malignant cells (Figure 13A). Moreover, referring to the published literatures, little is known about the function of C1QTNF6 compared to the other six genes. Consequently, we focused on the function of C1QTNF6. The violin plot showed that the expression of C1QTNF6 was the highest in fibroblasts, followed by endothelial cells, malignant cells and pan-immune cells (Figure 13B). By comparing multiple single-cell datasets, the similar expression trend was obtained (Figure 13C). Next, we analyzed the effect of C1QTNF6 expression in fibroblasts on infiltrating immune cells of TME by scTIME Portal website. GZMK⁺FOS⁺CD8⁺ T cells have been identified as memory T cells and are prevalent in para-carcinoma tissues or normal donors. FCGR3A⁺ NK cell is a classic NK cell cluster and plays important roles in anti-tumor immunity. C1QTNF6 expression in fibroblasts was dramatically negative relation to the infiltration of GZMK⁺FOS⁺CD8⁺ T cells (R=-0.515, P=0.049) and FCGR3A⁺ NK cell (R=-0.485, P=0.012) (Figure 13D, 13E). SPP1⁺ACP5⁺ macrophage and SPP1⁺CLEC5A⁺ macrophage are biased toward an M2 signature. Their infiltrating abundance was significantly positive to C1QTNF6 expression in fibroblasts (SPP1⁺ACP5⁺ macrophage: R=0.699, P<0.001; SPP1⁺CLEC5A⁺ macrophage: R=0.589, P<0.01) (Figure 13F, 13G). CTLA4⁺CD4⁺ Tregs are enriched in tumors and commonly involved in T cell inhibition across patients and tumor types. C1QTNF6 expression in fibroblasts was positive related to CTLA4⁺CD4⁺ Tregs abundance in TME, indicating high C1QTNF6 expression in fibroblasts might result in more CTLA4⁺CD4⁺ Tregs infiltration in TME (R=0.456, P=0.019) (Figure 13H). Interestingly, C1QTNF6 expression in fibroblasts was also positively associated with fibroblast infiltration in TME (R=0.547, P<0.01) (Figure 13I). Subsequently, we analyzed the cell communication between fibroblasts and immune cells. The results suggested that fibroblasts had strong interactions with M2 macrophages including C1QC⁺PLTP⁺ macrophage, SPP1⁺ACP5⁺ macrophage, SPP1⁺CLEC5⁺ macrophage (Supplementary Figure 6A). Furthermore, the ligand-receptor interaction analysis suggested that fibroblasts were very likely to interact with M2 macrophage through CD74-COPA, CD74-APP and CD74-MIF (Supplementary Figure 6B–6D). By analyzing the hallmarks enrichment between high and low C1QTNF6 expression groups, we found multiple pathways related to M2 polarization were enriched in high C1QTNF6 expression groups including NF-κB signaling pathway, glycolysis, IL6/JAK/STAT3 signaling pathway, TGF-β signaling pathway, Wnt/β-catenin signaling pathway, and Hedgehog signaling pathway (Supplementary Figure 6E). Consequently, we hypothesized that C1QTNF6 expression in fibroblasts may affect M2 macrophage polarization or recruitment.

Figure 13. Effect of C1QTNF6 on infiltrating immune cells of TME. (A) The heatmap of 7 genes expression in different cells by scRNA-seq. (B) The differences of C1QTNF6 expression among pan-immune cells, endothelial cells, fibroblasts and malignant cells. (C) The heatmap of C1QTNF6 expression in TME cells by multiple scRNA-seq datasets. The correlation of C1QTNF6 expression in fibroblasts and immune cells infiltration (memory CD8⁺ T cell (D), NK cell (E), M2 macrophages (F, G), Treg cell (H), fibroblast (I)). (J) qRT-PCR was performed to detect the efficiency of C1QTNF6-siRNA transfection. (K) M0 macrophages were stimulated by conditional medium from MRC-5 cells with C1QTNF6 silencing for 48h. qRT-PCR was performed to detect the expression of PD-L1, M1 and M2 markers. (L) Representatives and summary of M2 macrophage migration assays induced with MRC-5 cells with or without C1QTNF6 silencing. The data were presented as the mean±SD; n = 3. *P < 0.05, **P < 0.01, ***P < 0.001.

To validate the hypothesis, we constructed MRC-5 cells with C1QTNF6 silencing by siRNAs. As shown in Figure 13J, the specific cells were successfully established with high silencing efficiency. After 48h, the supernatant was collected, centrifuged and prepared as CM. To detect the effects of silencing C1QTNF6 in MRC-5 cells on macrophages polarization, we cultured M0 macrophages with the mixture of CM and FBS-containing medium (1:1) for 48h. Compared to the control, the group with C1QTNF6 silencing had significantly decreased PD-L1 and M2 macrophage-related genes expression (CD163, CD206). Conversely, M1 macrophage-related genes expression (CD80, CD86) were obviously increased when C1QTNF6 was silenced in MRC-5 (Figure 13K). Furthermore, we successfully induced M0 to M2 macrophage by IL-4 and IL-13 stimulation (Supplementary Figure 6F). The macrophage migration assay showed that silencing C1QTNF6 in MRC-5 cells could reduce M2 macrophage migration in vitro (Figure 13L). These fundings suggested that C1QTNF6 expression in MRC-5 cells promoted M2 macrophage polarization and recruitment.

Cause effect of C1QTNF6 on lung cancer onset

To evaluate the cause effect of C1QTNF6 on lung cancer onset, MR analysis was performed. 340 eligible SNPs were used as instrumental variables for C1QTNF6 and 274 common SNPs were obtained after harmonization. The funnel plot and leave-one-out sensitivity analysis showed that there were no obviously heterogeneous SNPs (Figure 14A and Supplementary Figure 7). MR analysis revealed that C1QTNF6 expression increased the risk of lung cancer (Figure 14B). Except weighted mode, the other four methods all showed the same trend (OR (95%CI): Inverse variance weighting (IVW) 1.029 (1.023-1.035) P < 0.001; MR Egger 1.015 (1.002-1.028) P =0.029; Weighted median 1.017 (1.008-1.027) P < 0.001; Simple mode 1.084 (1.052-1.117) P < 0.001) (Figure 14C). there was no heterogeneity and horizontal pleiotropy (the Cochrane’s Q-value > 0.1; MR PRESSO global test P = 0.236), indicated that the result of MR analysis was credible. Moreover, Steiger filtering further ensured directionality with all P-values of SNP less than 0.05. However, Bayesian co-localization showed that there was no genetic co-localization between C1QTNF6 and lung cancer (PP.H4=4.59e-03) (Figure 14D).

Figure 14. Mendelian randomization analysis of C1QTNF6 and lung cancer. (A) The funnel plot displayed the distribution of instrumental variables for C1QTNF6. (B) Scatter plot showed that C1QTNF6 increased the risk of lung cancer. (C) Forest plot showed the cause effect of C1QTNF6 on lung cancer onset. (D) The co-localization analysis of C1QTNF6 and lung cancer.

Data collection and processing

Gene expression data of LUAD and corresponding clinical characteristics were respectively retrieved from The Cancer Genome Atlas (TCGA) (HTseq-fragments per kilobase million, HTseq-FPKM) (https://cancergenome.nih.gov/) and Gene-Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) databases. Patients without prognostic information or survival time = 0 were excluded. Then 500 LUAD cases from TCGA database and 1009 LUAD cases from GEO database were retrieved as the training cohort and validation cohort (GSE3141: 58 cases; GSE13213: 117 cases; GSE30219: 83 cases; GSE31210: 226 cases; GSE50081: 127 cases; GSE72094: 398 cases). The demographic was shown in Supplementary Table 1. The ComBat method from the ‘SVA’ R package was used to remove the batch effects among different GEO datasets. Pan-cancer mRNA expression profiles and prognostic information were obtained from UCSC Xena website (https://xenabrowser.net/datapages/). The somatic mutation and copy number variation (CNV) of TCGA-LUAD were also curated from TCGA database. Eighty-four metabolic pathway gene sets were extracted from the KEGG database. The abbreviations of 84 metabolic pathways were listed in Supplementary Table 2. In addition, 42 LUAD tissues were collected from the Department of Thoracic Surgery, Shandong Provincial Hospital. The prognostic information was also followed up.

Seven immunotherapeutic cohorts were acquired to validate the prediction of immunotherapy efficacy using the MPPS model: advanced urothelial cancer treated with atezolizumab, an anti-PD-L1 antibody (IMvigor210 cohort) [62]; melanoma treated with anti-CTLA4 and anti-PD-1 therapy (GSE91061) [63]; metastatic melanoma treated with anti-CTLA4 therapy [64]; NSCLC treated with nivolumab or pembrolizumab, an anti-PD-1 antibody (GSE126044) [65]; NSCLC treated with anti-PD-1/PD-L1 antibody (GSE135222) [66]; melanoma treated with adoptive T cell therapy (ACT) (GSE100797) [67]; Melanoma treated with anti-PD-1 antibody (GSE78220). The response and benefit of TCGA cohort were calculated based on the Tumor Immune Dysfunction and Exclusion (TIDE) website (http://tide.dfci.harvard.edu/) by integrating TIDE score, interferon gamma (INFG), microsatellite instability (MSI) score, CD274, Merck18, CD8, myeloid-derived suppressor cells (MDSC), cancer associated fibroblast (CAF) and tumor-associated macrophages (TAM) M2.

To generate eQTL instruments for C1QTNF6, genetic variants located within 1000 kb on either side of the coding sequence (in cis) that are robustly associated with gene expression were extracted using eQTLs summary statistics obtained from the eQTLGen Consortium (https://www.eqtlgen.org/cis-eqtls.html). The data were established based on 26,609 blood samples of Europeans [68]. Lung cancer GWAS data (ieu-a-987) were obtained from the Transdisciplinary Research in Cancer of the Lung (TRICL). The GWAS data included 29,863 cases and 55,586 controls from European.

Estimation of metabolic pathways heterogeneity in LUAD

The levels of 84 metabolic pathways were estimated in each sample by ssGSEA. Then the metabolic differences of LUAD and normal samples were analyzed by “limma” package. An unsupervised consensus clustering according to 84 metabolic pathways scores was performed to identify distinct LUAD metabolic subtypes, which were showed by the principal component analysis (PCA). The metabolic profiles of TME cells including endothelial cells, malignant cells, cancer-associated fibroblasts (CAFs) and pan-immune cells were compared using single-cell RNA sequencing data (GSE111907).

Development and evaluation of a MPPS system

Pairwise comparisons of the 84 metabolic pathways’ scores in the training cohort were performed. The algorithm presented a scoring system in which the score of the metabolic pathway-pair was recorded as 1 if the expression level of the first metabolic pathway’s score was higher than that of the second; otherwise, it was recorded as 0, resulting in the construction of a 0 or 1 matrix. A metabolic pathway pair was deleted if the proportion of 0 or 1 was more than 80% or less than 20% of the samples in the training cohort. The abbreviations of the qualified metabolic pathway pairs were listed in Supplementary Table 2. Next, the qualified metabolic pathway pairs were enrolled for univariable, LASSO and multivariable Cox regression analysis to construct a MPPS system. The MPPS was calculated as follows:

$$\begin{array}{l}MPPS={\displaystyle \sum _{i=1}^{n}Coef(metabolicpathwaypairi)}\\ *Value(metabolicpathwaypairi)\end{array}$$

According to the optimal cut-off point of MPPS determined by the ‘survminer’ package, LUAD patients were stratified into the high- and low-MPPS groups. The optimal cut-off points were determined separately in the training and validation cohorts. The survival rates of the high- and low-MPPS groups were compared by Kaplan-Meier method in the training cohort and validation cohorts. Receiver operating characteristic (ROC) curves, generated by the “timeROC” package and “Aalen” weighting method, and C-index were used to detect the accuracy of MPPS. Univariable and multivariable Cox regression analyses were used to detect the prognostic roles of the clinical characteristics and MPPS. The independent prognostic factors were combined to develop a predicting nomogram by R package “rms”. The calibration curve was used to detect the consistency of the nomogram. The optimal cut-off value was calculated separately for each cancer, when we evaluated the performance of MPPS in pan-cancer.

Enrichment analysis and functional annotation

GO and KEGG pathway analyses were performed to investigate the variation in biological processes between high- and low-MPPS groups. The results of GO annotation were displayed by an online tool bioinformatics (https://bioinformatics.com.cn/). The Hallmark gene set was used to explore the distinction in various biological signatures between the high- and low-MPPS groups. The CancerSea website (http://biocc.hrbmu.edu.cn/CancerSEA/) was used to investigate 14 biological processes of multiple genes across various cancers.

Protein-protein interaction (PPI) network

The differentially expressed genes (DEGs) involved in MPPS between high- and low MPPS groups were input in STRING website (https://string-db.org/). PPI network was constructed with a minimum confidence score >0.4 and visualized by the software Cytoscape v3.9.1.

TME landscape analyses

Immune score, stromal score and ESTIMATE score were calculated using the ESTIMATE algorithm [69]. Immune cells infiltration and functions were evaluated by ssGSEA [70]. Expression of various immune checkpoint genes and HLA-related genes was compared between different MPPS scores groups. The cancer-immunity cycle scores of TCGA-LUAD samples were downloaded from the TIP database (http://biocc.hrbmu.edu.cn/TIP/) [71]. The discrepancy of the cancer-immunity cycle scores between the high- and low-MPPS groups was compared. The differences of MPPS scores among immune-inflamed, excluded and desert phenotypes had been analyzed [62].

Tumor mutation burden (TMB) and drug sensitivity analyses

The “maftools” R package was employed to explore the mutation frequency in different MPPS subgroups [72]. Then, the correlation between MPSS and TMB was analyzed. Subsequently, we evaluated the synergistic effect of TMB and MPPS score on prognostic stratification. A total of 137 drugs sensitivity in different MPPS groups was analyzed by R package “pRRophetic” and visualized in the form of parliament plot by Hiplot Pro (https://hiplot.com.cn) [73].

WGCNA

Co-expression gene networks of TCGA-LUAD were constructed using the WGCNA package. The unsigned network was selected. An appropriate soft threshold β was calculated to meet the criteria for the scale-free network. The optimal β was 4. Then, the weighted adjacency matrix was converted into a topological overlap matrix (TOM), and the corresponding dissimilarity was generated (1- TOM). Finally, the dynamic tree cut algorithm was used to identify the modules, and 80 was selected as the minimum number of genes for each module. The modules with the correlation coefficient R > 0.2, P-value < 0.05 were regarded as the key modules, and the genes with GS > 0.25, MM > 0.2 were regarded as key genes. These genes intersected by TCGA DEGs and GEO genes were used as candidate genes for subsequent analysis.

Hub genes identified from machine learning-based integrative approaches

Ten machine learning algorithms and 117 algorithm combinations were utilized to identify hub genes related to MPPS with high accuracy and stability performance on prognostic prediction. The integrative algorithms included random survival forest (RSF), stepwise Cox, elastic network (Enet), Lasso, Ridge, CoxBoost, partial least squares regression for Cox (plsRcox), generalised boosted regression modelling (GBM), supervised principal components (SuperPC), and survival support vector machine (survival-SVM). The procedure was as follows: (a) Univariable Cox regression analysis was used to screened out prognostic genes from the candidate genes; (b) Next, 117 algorithm combinations were performed on the prognostic genes to identified hub genes based on the leave-one-out cross-validation (LOOCV) framework in the TCGA-LUAD cohort; (c) All hub genes derived from 117 algorithm combinations were validated in six independent validation cohorts (GSE13213, GSE31210, GSE3141, GSE30219, GSE50081, GSE72094); (d) The hub genes with the highest average C-index across all validation cohorts were considered optimal.

scRNA-seq analysis

GSE111907 was retrieved to evaluate the metabolic pathway levels and hub genes expression in malignant, pan-immune cells, endothelial and fibroblast cells. The hub genes expression in various cell subtypes of TME was explored by TISCH2 website (http://tisch.comp-genomics.org/). GSE127465 was used to analyze intercellular communication and correlation of hub gene expression and immune cells infiltration using scTIME Portal website (http://sctime.sklehabc.com/unicellular/home).

RNA extracting and real-time PCR

Total RNA was extracted from LUAD frozen tumor tissues and cells using the AG RNAex Pro Reagent (Accurate Biotechnology (Hunan) Co., Ltd., China). The mRNA (500 ng) was converted into cDNA using Evo MMLVRT Master Mix kit (Accurate Biotechnology (Hunan) Co., Ltd., China). Then, cDNA was amplified with SYBR Premix Ex Tap kit (Accurate Biotechnology (Hunan) Co., Ltd., China). The mRNA levels were assayed by qRT-PCR using the Roche LightCycler® 480 system. 2^-ΔΔCt method was used to obtain relative quantitation (RQ) values, with 18S rRNA as endogenous control. The sequences of the primers were listed in Supplementary Table 3.

Cell culture and transfections

THP-1 cell was purchased from the Procell, Wuhan, China. MRC-5 cell was a gift from Shufang Chen (Shandong Provincial Hospital). MRC-5 cell was cultured in Dulbecco’s Modified Eagle Medium (DMEM) (HyClone, USA), and THP-1 cell was cultured in RPMI H1640 (HyClone, USA), supplemented with 10% fetal bovine serum (FBS) (BI, Israel) in a humidified atmosphere of 5% CO2 and 37° C according to protocol. C1QTNF6 siRNAs (Huzhou Hippo Biotechnology Co., Ltd., Zhejiang Province, China) were transfected into cells using jetPRIME (Polyplus-transfection, Illkirch, France) according to the manual. The sequences of C1QTNF6 siRNAs were as follows.

siC1QTNF6#1:

sense (5’-3’): GGAAUUACAAGGAGACGUA(dT)(dT)

antisense (5’-3’): UACGUCUCCUUGUAAUUCC(dT)(dT)

siC1QTNF6#2:

sense (5’-3’): GGGUCUUUGUGAACCUUGA(dT)(dT)

antisense (5’-3’): UCAAGGUUCACAAAGACCC(dT)(dT)

Preparation of conditioned medium (CM)

MRC-5 cells were transfected with C1QTNF6 siRNAs for 48h. The medium was replaced using serum-free medium and cells were cultured for additional 24h. Next, the supernatant was centrifuged at 300×g for 5 min and collected to induce TAM polarization and migration.

Polarization of THP-1 cells

To explore the effects of C1QTNF6 expression of MRC-5 on macrophage polarization, THP-1 cells were induced to M0 macrophage by phorbol 12-myristate 13-acetate (PMA) stimulation for 24h in six-well plates. Then, 2 ml mixture of CM and FBS-containing medium (1:1) was added for 48h. Finally, the total RNA was extracted and the M1- or M2-related markers were detected by qRT-PCR.

Macrophage migration assay

THP-1 cells were induced to M0 state under PMA stimulation. Then, the M0 macrophages were polarized into M2 macrophages via IL4 and IL13 stimulation. 20,000 M2 macrophages were plated in the upper chamber in the serum-free medium. The lower chambers were filled with 600 μl mixture of CM and FBS-containing medium (1:1). After 48h, the non-migrated cells in the upper chambers were removed and the migrated cells were stained with crystal violet for 30 min. The images were observed by microscope and the numbers of migrated cells were calculated by Image J.

Instrumental variable selection and Mendelian randomization

Only SNPs of C1QTNF6 cis-eQTL satisfying the following criteria were included as strong instrumental variables: (i) showed genome-wide significant association (P < 5 × 10 ⁻⁸); (ii) showed independent association [linkage disequilibrium (LD) clumping r² < 0.1; kb=500]; (iii) F-statistic > 10; (iv) not a palindromic SNP. Finally, 320 SNPs were identified as strong instrumental variables for C1QTNF6. For the MR analysis, the IVW method is the primary method. In addition, MR Egger, weighted median, simple mode and weighted mode methods were also used to detect the cause effect by R package “TwoSampleMR”. Leave-one-out sensitivity analysis was performed to evaluate the influence of each SNP on the outcome. Heterogeneity and potential horizontal pleiotropy were assessed by the Cochrane’s Q-value and MR-PRESSO global test. Steiger filtering was used to detect the directionality of the association between C1QTNF6 and lung cancer. Bayesian co-localization analyses were used to assess the probability that two traits share the same causal variant using the ‘coloc’ package (https://github.com/chr1swallace/coloc) with default arguments [68]. All SNPs within 1 Mb up and down stream of the leading SNPs were retrieved for colocalization analysis to analyze the posterior probability of H4 (PP.H4) PP.H4 > 80% was defined as having evidence of co-localization.

Statistical analyses

The statistical analysis of this study was performed using R v4.1.3, GSEA v4.2.3, GraphPad Prism 8 and SPSS v26. For quantitative data, the statistical significance of normally distributed variables was estimated by the Student’s t-test, and non-normally distributed variables were analyzed using the Wilcoxon rank sum test. When comparing between more than two groups, the Kruskal-Wallis test and one-way analysis of variance as non-parametric and parametric methods were made, respectively. Statistical significance was set at P < 0.05 unless otherwise stated. False discovery rate (FDR) was used to adjust P-value.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repositories and accession numbers have been listed in the article.