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

Immunotherapy has been a remarkable clinical advancement in cancer treatment, but only a few patients benefit from it. Metabolic reprogramming is tightly associated with immunotherapy efficacy and clinical outcomes. However, comprehensively analyzing their relationship is still lacking in lung adenocarcinoma (LUAD). Herein, we evaluated 84 metabolic pathways in TCGA-LUAD by ssGSEA. A matrix of metabolic pathway pairs was generated and a metabolic pathway-pair score (MPPS) model was established by univariable, LASSO, multivariable Cox regression analyses. The differences of metabolic reprogramming, tumor microenvironment (TME), tumor mutation burden and drug sensitivity in different MPPS groups were further explored. WGCNA and 117 machine learning algorithms were performed to identify MPPS-related genes. Single-cell RNA sequencing and in vitro experiments were used to explore the role of C1QTNF6 on TME. The results showed MPPS model accurately predicted prognosis and immunotherapy efficacy of LUAD patients regardless of sequencing platforms. High-MPPS group had worse prognosis, immunotherapy efficacy and lower immune cells infiltration, immune-related genes expression and cancer-immunity cycle scores than low-MPPS group. Seven MPPS-related genes were identified, of which C1QTNF6 was mainly expressed in fibroblasts. High C1QTNF6 expression in fibroblasts was associated with more infiltration of M2 macrophage, Treg cells and less infiltration of NK cells, memory CD8+ T cells. In vitro experiments validated silencing C1QTNF6 in fibroblasts could inhibit M2 macrophage polarization and migration. The study depicted the metabolic landscape of LUAD and constructed a MPPS model to accurately predict prognosis and immunotherapy efficacy. C1QTNF6 was a promising target to regulate M2 macrophage polarization and migration.


INTRODUCTION
Lung cancer remains a leading cause of cancer-related death worldwide [1].Lung adenocarcinoma (LUAD) is the most common histological subtype, accounting for 40% lung cancer cases [2].Although great progress has been made in LUAD treatment, the five-year survival rate of patients remains dismal [3].Immunotherapy has AGING led to striking clinical improvements while not all cancer patients can benefit from immunotherapy due to heterogeneity and adaptive evolution of tumor.Only about one third of patients acquire durable alleviation from it [4].To give patients more personalized medicine, it is essential to reveal the mechanism underlying distinct immunotherapy responses and develop signatures to predict prognosis and immunotherapy efficacy.
Recent studies revealed that oncogenic transformation induces a well-characterized metabolic phenotype in tumor cells, which in turn affects tumor microenvironment (TME) [5].As a new hallmark of malignant tumors, metabolic reprogramming improves malignant cells adaptation to meet bioenergetic, biosynthetic, redox balance demands and immune evasion, thus providing a selective advantage during tumorigenesis [5].Aerobic glycolysis (the Warburg effect) is a special metabolic pattern that tumor cells consume glucose and produces lactate even when oxygen is sufficient.Aerobic glycolysis not only provides enough ATP but also numerous precursor metabolites for lipids, amino acids, and nucleotides biosynthesis to support rapid proliferation [6].Dysregulated lipid metabolism is another prominent metabolic alteration in cancer [7].Under energy stressful conditions, tumor cells can harness lipid hydrolyzation to generate ATP and second messengers including diacylglycerol, arachidonic acid, lysophosphatidic acid, and phosphatidic acid to activate oncogenic signaling pathways [7][8][9].Other metabolic pathways such as amino acids metabolism, one carbon metabolism, purine and pyrimidine metabolism, are also dysregulated in tumor cells due to mutation of oncogenes, tumor suppressor genes or metabolic enzymes [10,11].Increasing evidence has suggested that tumor metabolic heterogeneity is greatly associated with TME status and immunotherapy [12][13][14].Several studies have suggested that glycolysis of tumor cells restricts glucose utility of tumor-infiltrating lymphocytes, thereby inducing T cells exhaustion and immune escape [15].Glutamine deprivation inhibits the transformation of CD4 + T cells to inflammatory subtypes, production and secretion of pro-inflammatory cytokines (IL-1, IL-6, and TNF) by macrophages, and promotes the apoptosis of immune cells [16][17][18].Large amount of lactic acid produced by tumor cells increases the acidity of TME and impairs the anti-tumor function of T cells and natural killer (NK) cells [19,20].Consequently, comprehensively depicting tumor metabolic landscape is promising to predict the prognosis and immunotherapy response of cancer patients and develop new treatment strategies.
Based on some metabolic features, many metabolic signature-based prognostic models have been established and acquired good predicting performance [14,21,22].
However, most models are based on a single metabolic pathway, lacking comprehensive exploration for tumor metabolism.Moreover, most models are constructed by focusing on exact gene expression and commonly unapplicable in another cohort sequenced by different platforms.To overcome the above shortcomings, we comprehensively assessed 84 metabolic pathways from 12 kinds of metabolism in LUAD by single sample gene set enrichment analysis (ssGSEA), developed and validated a metabolic pathway-pair score (MPPS) model to accurately predict the prognosis and immunotherapy efficacy of LUAD patients regardless of sequencing platforms.The model performed better than 51 published signatures of LUAD and was applicable to pan-cancers.The distinct metabolic features, TME between high-and low-MPPS groups were depicted.Weighted gene co-expression network analysis (WGCNA) and 117 machine learning algorithm combinations were performed and identified 7 MPPS-related genes, of which C1QTNF6 was mainly expressed in fibroblast.C1QTNF6 expression in fibroblast is positively related to fibroblasts, M2 macrophages and Treg cells infiltration but negatively related to memory CD8 + T cells and NK cells infiltration.Silencing C1QTNF6 expression in fibroblast impaired M2 macrophage polarization and migration in vitro assays.Meanwhile, Mendelian randomization (MR) also indicated that C1QTNF6 was cause of lung cancer onset.The MPPS model overcomes the obstacle of sequencing data from different platforms and is promising to guide LUAD patients' selection for immunotherapy.

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 AGING 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 panimmune 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.

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).
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).
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).

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.
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 upregulated and 19 metabolic pathways were downregulated (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 cancerimmunity 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.
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).

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 AGING 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 (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.AGING 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 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 lowrisk 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.

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 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).

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.2and 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 learningbased 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.

Expression, function, prognosis analyses of MPPSrelated 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  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 MPPSrelated 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).

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 panimmune cells (Figure 13B).By comparing multiple single-cell datasets, the similar expression trend was obtained (Figure 13C).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.
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 FBScontaining medium (1:1) for 48h.Compared to the control, the group with C1QTNF6 silencing had significantly decreased PD-L1 and M2 macrophagerelated 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 was obtained after harmonization.The funnel plot and leave-one-out sensitivity analysis showed that there was 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)  14D).

DISCUSSION
LUAD is a highly aggressive malignancy with an unfavorable prognosis and average 5-year survival rate of only 15% [24].With the rapid development of immunotherapy, it has shown great potential in the treatment of cancer.However, only about one third of patients can benefit from immunotherapy due to heterogeneity and adaptive evolution of tumor cells [4].To advance precision medicine, it is necessary to stratify cancer patients into distinct groups according to their prognosis and immunotherapy response before treatment.With the advances in sequencing technology, more and more gene expression-based prognostic models have been constructed to predict the prognosis and immunotherapy response of cancer patients [21,25,26].Unfortunately, the most models have not robust performance in other cohorts due to sequencing data from different platforms.To overcome this obstacle, pairing multiple markers to construct a prognostic model was put forward creatively.Metabolic reprogramming has been identified as a new hallmark of cancer and tightly associated with clinical outcomes and immunotherapy efficacy.Tumor cells reprogram their metabolism to compete for nutrients with other cells in TME, deal with oxidative stress, and reshape an immunosuppressive TME to evade the immune system [27,28].Comprehensively depicting the metabolic profile of LUAD is promising to predict the survival and immunotherapy efficacy of LUAD patients.
In this study, we assessed 84 metabolic pathways involved in 12 kinds of metabolism in LUAD by ssGSEA and analyzed the metabolic heterogeneity of LUAD.Then, we paired the 84 pathways and identified 19 metabolic pathway pairs by univariable, LASSO, multivariable Cox regression analysis.Using the 19 metabolic pathway pairs, we established a MPPS system and stratified LUAD patients into the high-and low-MPPS group.The high-MPPS group was characterized by high galactose metabolism, FA elongation, pyrimidine metabolism, cysteine and methionine metabolism, one carbon pool by folate and aminoacyl-tRNA biosynthesis.The low-MPPS group was characterized by dominant caffeine metabolism, valine, leucine and isoleucine degradation, selenocompound metabolism and arachidonic acid metabolism.Galactose is another important carbohydrate and involved in glycosylation, energy storage and pentose phosphate pathway directly or indirectly [29].Many tumors preferentially use glycolysis for survival and proliferation and have metabolic vulnerability to galactose.It has been reported that tumor cells with Akt activation will be induced cell death in galactose culture [30].Thus, LUAD with high galactose metabolism may be more adaptative for various energy substances.FA biosynthesis includes de novo synthesis and FA elongation.Elongation of very long-chain fatty acid (ELOVL) family enzymes are responsible for catalyzing FA elongation.Disruption of FA elongation by silencing ELOVL5 can suppress proliferation and invasion of renal cell carcinoma [31].Moreover, VLCFA deficiency results in a marked decrease in ceramides as well as downstream glucosylceramides and sphingomyelins, which impairs mitochondrial morphology and renders cancer cells sensitive to oxidative stress and cell death [32].Pyrimidine metabolism, one carbon pool by folate and aminoacyl-tRNA biosynthesis are tightly associated to nucleotides biosynthesis and translation, which are indispensable for rapid proliferation of malignant cells.UBE2T-mediated Akt K63 ubiquitination and Akt/βcatenin activation accelerate hepatocellular carcinoma development by increasing pyrimidine metabolism [33].Combination of pyrimidine synthesis inhibitors and other anti-tumor drugs is promising to kill tumor cells [34].Nucleotide synthesis and DNA methylation are highly dependent on one carbon pool by folate, which supports vital events for growth and survival [35].Methionine and cysteine, two of the most representative sulfur amino acids, play a crucial role in protein structure, metabolism, immunity, and especially, oxidation.They are extremely sensitive to almost all forms of reactive oxygen species and protect cells from oxidative stress damage [36].Dietary restriction of methionine and cysteine will alter the energetic metabolism and enhance the sensitivity of gliomas to ferroptosis [37].These metabolic pathways are highly elevated in high-MPPS group and may shape a refractory phenotype.Conversely, many anti-tumor metabolic pathways were elevated in the low-MPPS group.Caffeine can enhance anti-tumor activity of anti-PD-1 monoclonal antibody by increasing the infiltration of CD4 + and CD8 + T lymphocytes and decreasing the infiltration of Treg cells [38].
The branched-chain amino acids (BCAAs) (valine, leucine, and isoleucine) are essential amino acids that play important roles in metabolic regulation.The accumulation of BCAAs can activate mTOR signaling pathway to promote tumor proliferation [39].Thus, degradation of BCAAs may be harmful to tumor progression.Se compounds have been demonstrated as anticancer agents in vivo and in vitro experiments.They can prevent oncogene activation and cancer cell differentiation through scavenging of ROS, tumorpromoting eicosanoids and inducing tumor suppressor genes expression [40].Arachidonic acid metabolism is a double-edged sword in tumor initiation and progression.On the one hand, arachidonic acid can inhibit M2 macrophage polarization and enhance ferroptosis sensitivity to suppress tumor progression [41,42].On the other hand, it promotes stromal cell-mediated immunosuppression in NSCLC [43].
Consequently, the MPPS system divided LUAD patients into distinct metabolic reprogramming subgroups well.
The Increasing evidence demonstrates that metabolic reprogramming in TME affects anti-tumor immunity.
For example, targeting glutamine metabolism increased antitumor immunity in mouse models by upregulating mitochondrial metabolism of cytotoxic T lymphocytes in NSCLC [44,45].Treg cells rely on oxidative phosphorylation and FA oxidation to support their survival and differentiation [46].Lipid metabolic reprogramming can prevent effector T cells senescence and enhance immunotherapy efficacy [47].These also reveal that deeply understanding and depicting metabolic heterogeneity can favor immunotherapy.However, up to now, there is still a lack of comprehensive depiction of heterogeneous metabolic landscape in TME.The evaluation of 84 metabolic pathways in LUAD revealed the metabolic heterogeneity of LUAD in this study.
Considering the tight association of metabolism and immunotherapy, we wondered whether the LUAD patients with different MPPS had different responses to immunotherapy.Using seven independent immunotherapy cohorts, we found that the patients with low-MPPS scores commonly had higher immunotherapy response rates than those with high-MPPS scores.To further explore the alteration of immune cells, molecules and function, it was revealed that more immune cells infiltration, immune-related genes expression, and immune function activation were in the low-MPPS group, such as 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.The low-MPPS group also had higher cancer-immunity cycle scores in cancer antigen presentation, priming and activation, CD4 + T cell, dendritic cell, B cell, Th17 cell recruiting, and immune cells tumor infiltration.These results implied that LUAD with the low-MPPS score was inclined to be a "hot" TME and LUAD with the high-MPPS score was a "cold" TME.TIDE score also validated the conclusion and T cell dysfunction was higher in the low-MPPS group than the high-MPPS AGING group.With the increase of MPPS, the inflamed TME was transformed to the excluded TME.TMB is emerging as another indicator of immunotherapy except for PD-L1 expression.The high-MPPS group had higher TMB compared to the low-MPPS group.The LUAD patients with high TMB and low MPPS had the best prognosis and those with low TMB and high MPPS had the worst prognosis.Consequently, the bad prognosis of the high-MPPS group is not likely due to TMB.By identifying the sensitivity of 137 chemotherapy drugs, multiple drugs sensitive to the high-or low-MPPS group were determined, which may be helpful to guide precision medicine of LUAD patients.
Although targeting cancer metabolism to improve immunotherapy efficacy is highly promising, the crosstalk of metabolic pathways between tumor cells and immune cells in TME lead to disruption of normal metabolic pathways in immune cells by strategies to inhibit/alter cancer metabolism [48].Thus, it is critical to target the specific metabolic pathways to kill tumors without interfering with or even promoting anti-tumor immunity.To identify such pathways, we analyzed the metabolic pathway pair in different kinds of cells by scRNA-seq data.Interestingly, the average value of cysteine and methionine metabolism/ganglio series biosynthesis is significantly elevated in malignant cells than the other cells including immune cells, fibroblasts, endothelial cells.Many previous studies have reported that tumor cells are highly dependent on cysteine and methionine metabolism than normal cells and they are promising targetable weaknesses of cancer cells [49].Ganglio series biosynthesis are also tightly related to some malignant phenotypes such as metastasis [50].As a result, this metabolic pathway pair may be promising to be a metabolic target in LUAD therapy.
The previous studies mostly choose the modeling algorithms to identify the hub genes based on their knowledge limitations and preferences.To overcome this shortcoming, we firstly identified the MPPS-related hub gene module by WGCNA and then, integrated 117 machine learning algorithms to further recognize the prognostic signature.Finally, seven genes were identified, in which C1QTNF6 caught our attention due to its specific expression in fibroblast.Some studies have suggested that silencing C1QTNF6 in LUAD cells can suppress the proliferation, migration and invasion of LUAD cells [51].C1QTNF6 is a prognostic indicator for poor survival across many cancers including LUAD and one of the most relative genes of TAM [52,53].However, there is still little knowledge about the function of C1QTNF6 in tumors.
By analyzing multiple scRNA-seq datasets, we found C1QTNF6 expression was mainly focused on fibroblast and its expression in fibroblast was positively related to the infiltration of M2 macrophages, Treg cells, and negatively related to the infiltration of memory CD8 + T cells, NK cells.Moreover, there existed strong interaction between M2 macrophages and fibroblast by intercellular communication analysis.In vitro experiments also validated that the CM from fibroblast C1QTNF6-/-would promote the transformation of M0 into M1 but not M2 macrophage, decrease PD-L1 expression, and reduce M2 macrophage migration.Hallmarks enrichment analysis showed that NF-κB signaling pathway, glycolysis, IL6/JAK/STAT3 signaling pathway, TGF-β signaling pathway, Wnt/β-catenin signaling pathway, and Hedgehog signaling pathway were enriched in high C1QTNF6 expression group, which were reported to participate in M2 macrophage polarization.Inhibition of autophagic degradation of RELA will rescue activity of NF-κB signaling pathway and shape the phenotype of hepatoma-polarized M2 macrophages [54].Activation of IL6/JAK/STAT3 signaling pathway in macrophages can promote M2 polarization and PD-L1 expression [55,56].A large amount of lactate produced by glycolysis induces M2 macrophage polarization and promotes the invasion of pituitary adenoma [57].Mesenchymal stem cells can induce M2 polarization phenotype via secreting TGF-β to activate Akt/FoxO1 pathway in LPSstimulated macrophages [58].It is also reported that crosstalk between hepatic tumor cells and macrophages by Wnt/β-catenin signaling pathway can promote M2 polarization [59].FOXM1 can induce M2 polarization through SEMA3C/NRP2/Hedgehog signaling [60].
The results indicated that C1QTNF6 may be tightly associated with M2 polarization.Lin et al. reported that after silencing C1QTNF6, the enrichment of cytokinecytokine receptor interaction pathways was reduced in LUAD cell by RNA sequencing, which indicated that C1QTNF6 may participate in cytokine-cytokine receptor interaction pathways directly or indirectly [61].Consequently, C1QTNF6 expression in fibroblast may promote M2 macrophage polarization and migration by regulating cytokine-cytokine receptor interaction.
Mendelian randomization (MR) is as a valuable tool for inferring causal relationships between exposure and outcome by leveraging Genome wide association study (GWAS) data.The result of MR suggested that C1QTNF6 expression had the increased risk of lung cancer although there was no evidence of colocalization.The MR result was consistent with the expression and prognosis of C1QTNF6 in LUAD.
There are still some limitations in our study.Although we identified two distinct metabolic subtypes with significantly different prognosis and immunotherapy efficacy, some immunotherapy cohorts were from the studies about urothelial cancer or melanoma and more AGING LUAD-related immunotherapy cohorts are needed to validate our conclusion.The drug sensitivity needs further validation by IC50 assays.Although we identified the potential metabolic pathways associated with prognosis, the targetable molecules for the pathways remain to be explored.The underlying mechanism that C1QTNF6 regulated M2 macrophage polarization and migration remains to be elucidated.Moreover, the relationships of C1QTNF6 and the other immune cells need further exploration.The conclusion of MR needs further experimental validation.The above insufficient will be the focus of our future study.

CONCLUSIONS
Based on 84 metabolic pathways, we constructed a MPPS model to accurately predict the prognosis and immunotherapy efficacy of LUAD patients.Targeting C1QTNF6, a MPPS-related gene, is promising to suppress M2 macrophage polarization and migration.

Data collection and processing
Gene Provincial Hospital.The prognostic information was also followed up.

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 AGING 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: 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 highand 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, AGING 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.compgenomics.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.

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 13acetate (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 nonmigrated 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 −8 ); (ii) showed independent association [linkage disequilibrium (LD) clumping r 2 < 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.

Figure 1 .
Figure 1.The flow chart of the study.

Figure 2 .
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.

Figure 3 .
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.

Figure 5 .
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.

Figure 6 .
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.

Figure 7 .
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.

Figure 8 .
Figure 8. Prediction of immunotherapy by MPPS model.Survival analysis (A) and response to anti-PD-L1 therapy (B) between the

Figure 9 .
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.

Figure 13 .
Figure 13.Effect of C1QTNF6 on infiltrating immune cells of TME.(A) The heatmap of 7 genes expression in different cells by scRNAseq.(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.

Figure 14 .
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.

Table 1 . The results of univariable and multivariable Cox regression analyses.
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.AGING (approval number: SWYX: NO. 2023-283).The requirement for written informed consent was waived.All procedures involving human participants were performed in accordance with the 1975 Helsinki declaration and its later amendments.