Skip to main content

Integration of QTL and comprehensive analysis in the circulating inflammatory cytokines for pan-cancer

Abstract

Background

Chemokines and cytokines are components of the tumor microenvironment and also influence tumorigenesis and its composition. However, whether they genetically proxy tumorigenesis is unclear. For causal inferences, eQTL and pQTL were used to determine the role of chemokines and cytokines in pan-cancer. The impact on the tumor immune microenvironment was also explored.

Methods

This study leveraged summary statistics from respective genome-wide association studies (GWAS) of 109 cytokines and chemokines in 18 types of solid tumors. Single nucleotide polymorphisms (SNPs) robustly associated with the cytokines and chemokines, located in or close to their coding gene (cis), were used as instrumental variables. A two-sample MR design was employed, followed by comprehensive sensitivity analyses to validate the robustness of results. The impact on immune infiltration was investigated using the TIMER and TISIDB websites. Survival analysis was conducted using the K-M plotter and TIMER 2.0 websites. The TISCH and GEO databases were used to carry out scRNA cell analysis.Analyzing relevant proteins using the STRING database and conducting enrichment pathways for GO analysis of the identified proteins.

Results

The results of the inverse-variance weighted (IVW) method using cis-protein QTL (cis-pQTL) instruments showed the causal effects of TNF in reducing the risk of squamous cell lung cancer (LUSC) and HGF in reducing the risk of head and neck cancer (HNSC).The results were consistent with the eQTL. HGF was associated with better overall survival (OS) in HNSC, regardless of the types of cells enriched. However, high expression of the ligand MET for HGF leads to a decrease in overall survival in LUSC. TNF was related to poor OS in LUSC with no significant impact. However, in CD8 + T cell-enriched, eosinophil-enriched, macrophage-enriched, and NK cell-deficient types of LUSC, high expression of TNF leads to a poor prognosis, and there is statistical significance. The results showed a significant positive correlation between TNF and most immune cell infiltration, immunomodulator and chemokine in LUSC. HGF is positively correlated with the majority of immune cells except CD56 + cells, as well as some immune regulatory factors and chemotactic factors. According to single-cell sequencing results, HGF is mainly secreted by fibroblasts and myofibroblasts in HNSC, while in LUSC, it is primarily secreted by macrophages and CD8 + T cells secrete TNF. The GO/KEGG analysis suggests that proteins related to HGF are mainly involved in regulating peptidyl-tyrosine phosphorylation and positive regulation of the MAPK cascade. Proteins related to TNF are primarily associated with the regulation of I-kappaB kinase/NF-kappaB signaling and cytokine-mediated signaling pathway.

Conclusions

HGF is primarily secreted by fibroblasts in HNSC and may have a protective effect on the occurrence and prognosis of HNSC. These effects are independent of immune cell influence, and this role may not necessarily be mediated through the HGF/MET pathway. On the other hand, TNF in LUSC is mainly secreted by immune cells like CD8 + T cell, and it may have a protective effect on the occurrence of LUSC. However, it’s impact on the prognosis of LUSC through the immune microenvironment may have a different effect.

Plain Language Summary

Chemokines and cytokines are not only components of the tumor microenvironment but also affect tumorigenesis and the composition of the tumor microenvironment. However, whether they genetically proxy tumorigenesis is unclear. For causal inferences, eQTL and pQTL were used to define the role of chemokines and cytokines in pan-cancer. The impact on the tumor immune microenvironment was also explored. This study leveraged the summary statistic from respective genome wide association study (GWAS) of 109 cytokines and chemokines to 18 types of solid tumor. Single nucleotide polymorphisms (SNPs) robustly associated with the cytokines and chemokines, located in or close to their coding gene (cis), were used as instrumental variables. A two-sample MR design was employed, followed by comprehensive sensitivity analyses to validate the robustness of results. The results showed HGF is primarily secreted by fibroblasts in HNSC, and it may have a protective effect on the occurrence and prognosis of HNSC. These effects are independent of immune cell influence, and this role may not be mediated through the HGF/MET pathway. On the other hand, TNF in LUSC is mainly secreted by immune cells like CD8 + T cell, and it may have a protective effect on the occurrence of LUSC. However, it’s impact on the prognosis of LUSC through the immune microenvironment may have a different effect.

Peer Review reports

Introduction

The aggregation of biological evidence has solidified inflammation as a recognized characteristic of cancer. The concept suggests that a condition of mild inflammation could elevate mutation rates and enhance the growth of mutated cells by providing trophic signals [1]. Apart from potential direct promotion of cell growth, activated inflammatory cells have the capability to induce reactive oxygen species and the buildup of reactive nitrogen intermediates in adjacent cells [2]. These mechanisms may cause damage to DNA and its protein products, either directly or indirectly, thereby exerting tumor-promoting effects [3]. On the other hand, immune surveillance plays a role in anti-tumor effects and has toxicity effects on tumors, exhibiting anti-cancer effects in specific circumstances [4]. The causal relationship between specific immune inflammatory factors and tumor occurrence is still unclear, as well as their role in the subsequent development of tumors. Observational findings indicate that specific inflammatory markers in circulation have been correlated with cancer development in prospective cohort studies [5]. For instance, higher pre-diagnostic levels of interleukin 1 alpha (IL-1a), IL-8, and tumor necrosis factor alpha (TNF-A) have been associated with an elevated risk of ovarian cancer, while concentrations of serum amyloid A, soluble tumor necrosis factor receptor-2 (sTNF-RII), and monokine induced by gamma interferon (MIG) have shown positive associations with lung cancer risk [6, 7]. If these associations prove to be causal, intervening in or preventing inflammation pathways could be a strategy to mitigate cancer risk.

While there has been extensive discussion on the potential therapeutic benefits of reducing chemokine levels to inhibit tumorigenesis, observational studies linking specific circulating inflammatory cytokine concentrations to cancer risk are limited, often characterized by relatively small sample sizes, and susceptible to potential biases such as unmeasured confounding and reverse causation. To address these challenges and strengthen the evidence for a potential causal role of inflammation in cancer risk, we employed Mendelian Randomization (MR). In MR, germline genetic variants serve as instrumental variables to proxy tumorigenesis exposure to a particular factor of interest, in this case, circulating cytokines, chemokines, growth factors, and interferons (collectively referred to as cytokines). Our study utilized MR to capture the typical cytokine concentration experience over the tumorigenesis, focusing on usual levels rather than expression variations resulting from epigenetic alterations. We leveraged genetic variants robustly associated with circulating cytokines to estimate the association between genetically proxied inflammatory cytokine concentrations and the risk of pan-cancer. Additionally, we utilized outcome data from well-established consortia, publicly available for analysis. Furthermore, we conducted a comprehensive examination and assessed the potential implications of circulating inflammatory cytokines on cancer prognosis and their association with immune infiltration.

Methods

Study design

The flowchart of the current work is outlined in Fig. 1. To begin, we selected genetic variants as instrumental variables (IVs) for circulating cytokines. Next, we collected the latest summary statistics from the corresponding genome-wide association studies (GWASs) for pan-cancers. Then, we performed a two-sample MR study with the IVW method as our primary analysis [8]. And a series of sensitivity analyses were conducted, including MR-Egger regression, weighted-median, contamination mixture (ConMix), and MR-Pleiotropy Residual Sum and Outlier (MR-PRESSO) [9, 10]. We identified inflammatory factors with positive eQTL and their associated tumor types. Subsequently, we selected SNPs representative protein-level qQTL as IVs for MR analysis and related validation in eQTL-positive tumors. Finally, we conducted a comprehensive analysis of cytokines that are positive for both eQTL and qQTL, along with their corresponding cancer types, including prognostic analysis, immune infiltration analysis, the function and enrichment of interacting protein analysis, and scRNA cell analysis.

Fig. 1
figure 1

Flowchart of the analyses performed

Data sources and IVs selection

To investigate circulating inflammatory cytokines as exposures, we obtained IVs for eQTL from the latest GWAS meta-analysis across three independent sources [11] [specifically, the SCALLOP consortium [12], the INTERVAL study [13], and the Northern Finland Birth Cohort 1966 [14]], encompassing a total of 31,112 individuals of European ancestry. The original paperc [11] provides detailed information on SNP selection methods and meta-analysis procedures. Additionally, we incorporated pQTL data from the integrated version [15] of five previously published GWAS datasets [16,17,18,19,20] for our primary analysis, identifying 738 cis-pQTLs associated with 734 plasma proteins. Supplement Table 1 presents the characteristics of 18 types of cancer as the outcome data utilized in this study.

MR relies on three key assumptions for accurate causal inferences (1): IVs must be genuinely associated with exposures (2), independent of confounding factors, and (3) exclusively influence outcomes through exposures, without involvement in other pathways. To bolster confidence in the fulfillment of these assumptions, particularly the first and the third, we opted for cis quantitative trait loci (QTLs) as IVs to enhance instrument strength. This includes cis-protein QTLs (cis-pQTLs, within a gene range of ± 500 kb) and cis-expression QTLs (cis-eQTLs, within a gene range of ± 500 kb).

Cis-QTLs, situated at or in proximity to the gene of origin, inherently exhibit a stronger correlation with gene expression and subsequently protein concentrations compared to other variants. In our primary analysis, we employed cis-pQTLs as the IVs. For supplementary analysis, we utilized cis-eQTLs, as these may capture the effects of pQTLs through gene expression, although not all pQTLs are represented by eQTLs [21]. Specifically, both cis-pQTLs and cis-eQTLs were associated with circulating cytokine levels and gene expression aggregated across tissues.

It’s important to note that, for a balanced number and strength of instrumental variables and potentially informative results, we selected a relatively lenient threshold of 1E-4, deviating from the traditional GWAS threshold of 5E-8. Palindromic variants with a minor allele frequency (MAF) greater than 0.05 were excluded, and clumping was performed with a pairwise linkage disequilibrium (LD) cutoff of r2 < 0.1. Harmonization of QTL alleles between the exposure and outcome was conducted to ensure proper alignment of effects.

Mendelian randomization analyses

We conducted distinct analyses using two sets of instruments (cis-pQTL and cis-eQTL) to explore the relationships between genetically influenced circulating cytokine concentrations and the risk of each cancer outcome. In cases where a single SNP was available to construct the instrumental variable, we employed the ratio of coefficients method to derive MR estimates, utilizing first-order weights for standard error generation. For situations with multiple SNPs available for a given cytokine’s IVs, MR estimates from individual SNPs within the instrument were combined using the IVW MR method.

To address multiple hypothesis testing, we computed false discovery rate (FDR)-adjusted p values (q values) in the primary IVW MR analyses using the sequential p value approach by Benjamini and Hochberg [22]. A q value not exceeding 10% was considered significant. Effect estimates indicate the increase in cancer risk per standard deviation higher on the natural scale of each cytokine. For instruments showing evidence of association in the MR analyses, we further assessed potential pleiotropic effects using Phenoscanner, a database incorporating genotype-phenotype associations [23].

Immune infiltration analysis

The “GENE” module of the TIMER database was used to assess the extent of immune cell infiltration across 18 types of cancer. Additionally, we employed TISIDB (cis.hku.hk/TISIDB/index.php) to investigate the correlation between the “GENE” and major histocompatibility complexes (MHCs), as well as chemokine receptors.”

Prognostic analysis

We conducted a prognostic analysis of the potential significance of the “gene” in cancer prognosis, considering enrichment in different cell types. This analysis involved assessing overall survival (OS) using both the K-M plotter and the TIMER2.0 database.

The function and enrichment analysis

We conducted gene co-expression analysis on “gene” using LinkedOmics (www.linkedomics.org/login.php), specifically choosing the “HiSeq RNA” platform and the “TCGA_LUAD” cohort for examination. The correlation between “gene” and its co-expression genes was determined using the Pearson test. GO analysis was used to anaylsis function. And pathview analysis of key KEGG pathway.

Single-cell sequencing analysis

In TISHC analysis, the secretion of target genes in the respective cancer types is examined to identify the components responsible for their expression.

Statistical analyses

We employed two distinct sets of instrumental variables (IVs), namely cis-pQTL and cis-eQTL, in our analyses to investigate the associations between genetically predicted circulating cytokine levels and the risk of pan-cancer. In cases where only a single QTL was present, we utilized the classic Wald ratio for MR estimates. Alternatively, for scenarios with multiple QTLs, we employed the random-effects IVW model through a meta-analysis approach to combine Wald ratios and obtain overall MR estimates [24]. The IVW method is generally considered unbiased, given the fulfillment of three aforementioned assumptions.

Furthermore, we conducted sensitivity analyses employing four pleiotropy-robust methods: MR-Egger regression [25], weighted-median method [26], ConMix [27], and MR-PRESSO [28]. These methods offer valid results under different assumptions, enhancing the robustness of our conclusions. Cochran’s Q test in IVW and the intercept from MR-Egger method were used to assess heterogeneity and horizontal pleiotropy, respectively (P-value < 0.05). Information on R2 variance and F-statistics, quantifying IV strength, was provided, and cumulative F-statistic > 10 was considered to avoid weak instrument bias [29].

The statistical power of our entire MR analyses was estimated using the non-centrality parameter-based approach by Brion et al. via an online tool (https://shiny.cnsgenomics.com/mRnd/) [30]. Correlation analysis was conducted to illustrate mutual corroboration or complementarity of results under both instrumental strategies.

All statistical analyses were implemented using the “TwoSampleMR” (version 0.5.6), “MendelianRandomization” (version 0.6.0), and “MR-PRESSO” (version 1.0) packages in R (version 4.1.2). To address multiple comparisons for numerous cytokines, we performed stratified false discovery rate (FDR) analysis using the Benjamini-Hochberg procedure for IVW analysis, estimating adjusted P-values separately for pan-cancer [31]. Additionally, we applied aggregated FDR correction to complement our results. Statistical significance was defined with a threshold of adjusted P-value < 0.10.

Result

eQTL analysis of cytokines and cytokines across pan-cancer

To investigate the causal effect of all analyzable cytokines (listed in Supplement Table 2) on the risk of pan-cancer, including gastric cancer and lung cancer, MR tests were conducted. Out of the 109 analyzed cytokines, 51 had cis-eQTL instruments, explaining 0.1–28.9% of the phenotypic variance. Cumulative F-statistics for cis-pQTL instruments of 29 out of 31 cytokines were greater than 10, indicating the robust strength of genetic instruments (see Supplement Table 2). Detailed MR results for the causal relationships of interest are presented in Supplementary Table 3, based on the cis-eQTL instruments. The visualization of all MR results is depicted in Fig. 2.

Using cis-eQTL instruments in IVW and Wald ratio analysis, genetically predicted higher levels of 13 cytokines showed a suggestive association with an increased risk of different solid tumors. Specifically, IL4 had a suggestive negative association with gastric cancer (interleukin 4, b: -0.451, be: 0.128, P: 0.0004), TNF had a suggestive negative association with hepatocellular carcinoma (tumor necrosis factor, b: -1.339, be: 0.437, P: 0.0022), FLT3LG had a suggestive negative association with lung adenocarcinoma (fms-related tyrosine kinase 3 ligand, b: -0.297, be: 0.096, P: 0.0019), CD40LG and IL18 had suggestive positive associations with lung adenocarcinoma (CD40 ligand, b: 0.644, be: 0.158, P: 0.00004, interleukin 18, b: 0.192, be: 0.043, P: 0.000007), CD70 and TNF had suggestive negative associations with squamous cell lung cancer (CD70 molecule, b: -0.292, be: 0.086, P: 0.0007, tumor necrosis factor, b: -0.267, be: 0.045, P: 0.000000003), CSF1 had a positive association with squamous cell lung cancer (colony stimulating factor 1, b: 0.296, be: 0.086, P: 0.00055), XCL1 had a negative association with thyroid cancer (chemokine (C motif) ligand 1, b: -0.747, be: 0.161, P: 0.00000352), HGF had a negative association with head and neck cancer (hepatocyte growth factor, b: -0.005, be: 0.001, P: 0.00133), IL2RA had a negative association with cervical cancer (interleukin 2 receptor, alpha, b: -0.009, be: 0.002, P: 0.0000012), LTB had a positive association with cervical cancer (lymphotoxin beta, b: 0.004, be: 0.001, P: 0.0008), and FLT3LG had a negative association with malignant neoplasm of the skin (fms-related tyrosine kinase 3 ligand, b: -0.006, be: 0.001, P: 0.000188).

Fig. 2
figure 2

In this summary of MR-IVW or MR-Wald ratio results using the cis-eQTL instrument definition, squared tiles represent associations that are nominally significant (p < 0.0027), and an asterisk indicates significance after correcting for multiple comparisons (FDR ≤ 10%). The color scale reflects MR beta estimates, with white tiles indicating associations for which no instrument was available

Validation of the positive eQTL results using qQTL

Utilizing cis-pQTL instruments in the IVW analysis, genetically predicted higher levels of HGF showed a suggestive association with a decreased risk of HNSC (Hepatocyte Growth Factor, odds ratio [OR]: 0.9985, 95% confidence interval [CI]: 0.9975–0.9996, P: 0.006). Additionally, higher levels of TNF were suggestively associated with a decreased risk of LUSC (Tumor Necrosis Factor, odds ratio [OR]: 0.5875, 95% confidence interval [CI]: 0.4523–0.7632, P: 6.74E-5). Moreover, the reverse Mendelian randomization results suggest that when LUSC as the exposure and TGF as the outcome, the IVW results P = 0.676. Similarly, when HNSC is the exposure and HGF is the outcome, the IVW results also P = 0.861, indicating no reverse causal relationship (Supplementary Tables 5 and Supplementary Figs. 12).

As demonstrated in Fig. 3 and A (scatter plot), Fig. 3B (forest plot), Fig. 4A (scatter plot), and Fig. 4B (forest plot), a causal connection was found by the IVW analysis between HGF and HNSC (27 SNPs as IVs, Supplementary Table 4), and between TNF and LUSC (14 SNPs as IVs, Supplementary Table 4). the funnel plot (Figs. 3D and 4D), a revealed no significant heterogeneity in the observed associations. And the leave-one-out analyses (Figs. 3C and 4C) indicated that the observed associations remained stable, even after the removal of any single SNP.

Focusing on associations that withstand multiple comparison correction (stratified FDR < 10%), their corresponding sensitivity analyses indicated similar estimates, although several methods yielded wide confidence intervals due to lower statistical power (refer to Fig. 5). Further examination revealed little evidence of heterogeneity (the majority of Cochran Q statistic P-values > 0.05) or horizontal pleiotropy (the majority of MR-Egger intercept P-values > 0.05).

Fig. 3
figure 3

The casual impact of HGF on HNSC. (A) Scatter plot showing the relationship between HGF and HGF. All four of the approaches used in the current manuscript were shown. (B) A forest plot was utilized to display the MR estimation and 95% CI values for each SNP. At the bottom are the MR-Egger and IVW findings. (C) Leave-one-out tests to determine if a single instrumental variable was responsible for the causal impact. (D) A funnel plot was used to determine if the connection that was seen was accompanied by noticeable heterogeneity

Fig. 4
figure 4

The casual impact of TNF on LUSC. (A) Scatter plot showing the relationship between TNF and LUSC. All four of the approaches used in the current manuscript were shown. (B) A forest plot was utilized to display the MR estimation and 95% CI values for each SNP. At the bottom are the MR-Egger and IVW findings. (C) Leave-one-out tests to determine if a single instrumental variable was responsible for the causal impact. (D) A funnel plot was used to determine if the connection that was seen was accompanied by noticeable heterogeneity

Fig. 5
figure 5

Forest plots of IVW and sensitivity analyses. Shown are the associations that withstand multiple comparison correction (stratified FDR < 10%) based on the (A) cis-pQTL of immune cytokines and cancer species which was cis-eQTL positive and (B) IVW and sensitivity analyses of cis-pQTL and cis-eQTL positive

The prognostic value of TNF in LUSC and HGF in HNSC

We conducted a comprehensive analysis of the prognostic significance of HGF in various types of Head and Neck Squamous Cell (HNSC) cancers using the K-M plotter. As depicted in Fig. 6, heightened expression of HGF correlated with improved Overall Survival (OS) in HNSC (p = 6.8e − 05). This favorable prognosis was consistent across all subtypes of HNSC, indicating that HGF may not enhance HNSC prognosis through CD8 cells, CD4 cells, other immune cells, or mesenchymal stem cells. The ligand MET, which binds to HGF, exhibited elevated expression linked to poor prognosis in HNSC, while the activated protein HGFAC improved HNSC prognosis, independent of immune cell content (Fig. 6E). Considering that HGF doesn’t seem to exert its effects through the HGF-cMET pathway or modulation of the immune microenvironment to enhance HNSC prognosis.

Fig. 6
figure 6

Illustrates the impact of HGF on the prognosis of different types of Head and Neck Squamous Cell (HNSC) tissues and the effect of TNF on the prognosis of various types of Lung Squamous Cell Carcinoma (LUSC) tissues. (A). The overall prognostic influence and the impact on the prognosis of tumors enriched or deficient in mesenchymal stem cells, macrophages, and NK cells. (B). The effect on the prognosis of tumors enriched or deficient in B cells, CD4 + T cells, CD8 + T cells, and regulatory T cells. (C). The influence on the prognosis of tumors enriched or deficient in eosinophils, neutrophils, type 1 helper T cells, and type 2 helper T cells. (D). The impact on the prognosis of tumors with different levels of tumor neoantigens and high or low tumor mutational burden (TMB). (E). The overall prognostic influence of HGFAC and MET on HNSC

In contrast, heightened expression of TNF was associated with poorer OS in Lung Squamous Cell Carcinoma (LUSC), although the difference was not statistically significant (p = 0.052) (Fig. 6). Among 319 HNSC patients with enriched macrophages, TNF significantly correlated with poorer OS (p = 0.021). However, in 173 HNSC patients with decreased macrophages, TNF was associated with better OS, though not statistically significant (p = 0.17). Similar trends were observed for eosinophils. This opposing trend suggests that TNF may negatively impact LUSC prognosis through macrophages and eosinophils. Additionally, in 54 LUSC patients with enriched NK cells, TNF improved prognosis without significant differences (p = 0.12). Conversely, in 438 LUSC patients with decreased NK cells, TNF negatively affected prognosis (P = 0.037). This opposing trend suggests that NK cells may inhibit the negative impact of TNF on LUSC prognosis. The trends were verified in the TIMER2.0 database.

Fig. 7
figure 7

Illustrates the impact on the immune microenvironment: (A). The influence of HGF in pan-cancer on lymphocytes, immunomodulators, and immunoinhibitors. The black frame representing the effects of HGF in HNSC. The last two depicts the effects of HGF in HNSC on GZMA and PRF1. (B). The impact of TNF in pan-cancer on lymphocytes, immunomodulators, and immunoinhibitors. It further showcases the effects of TNF in Lung Squamous Cell Carcinoma (LUSC) on GZMA and PRF1, with the black frame indicating LUSC

The impact on the immune microenvironment

TIMER 2.0 and TISIDB were employed to investigate the potential impact of TNF in Lung Squamous Cell Carcinoma (LUSC) and HGF in Head and Neck Squamous Cell Carcinoma (HNSC) on the tumor immune microenvironment.

Fig. 8
figure 8

T-SNE diagram (A-B). HGF expression profles were demonstrated at single-cell levels from HNSC by T-SNE diagram. (C). MET expression profles were demonstrated at single-cell levels from HNSC by T-SNE diagram. (D-E). TNF expression profles were demonstrated at single-cell levels from LUSC by T-SNE diagram

In LUSC, TNF upregulation was associated with nearly all lymphocytes, immunomodulators, including Immunoinhibitors, Immunostimulators, Major Histocompatibility Complexes (MHCs), as well as chemokines and their receptors. Furthermore, TNF expression showed a positive correlation with GZMA and PRF1 (cor = 0.207, cor = 0.35, P<0.05), representing levels of cytotoxic T lymphocytes (CTL) in HNSC (Fig. 7B).

Similarly, HGF played a similar role in influencing the immune microenvironment of Head and Neck Squamous Cell Carcinoma (HNSC). In HNSC, HGF was negatively correlated with CD56 cells but positively correlated with all other immune cells. Regarding Immunoinhibitors, HGF was negatively correlated only with PVRL2 and TNFB1, while its correlation with KIR2D and CD160 was not strong. It exhibited positive correlations with other Immunoinhibitors. In Immunostimulators, HGF showed negative correlations with PVR, RAET1E, TNFRSF25, TNFSF9, and weak correlations with CD70 and TNFRSF18. Positive correlations were observed with other Immunostimulators. Concerning MHCs, HGF displayed negative correlations with B2M, HLA-A, HLA-B, HLA-C, HLA-E, HLA-F, HLA-G, TAP1, TAP2, TAPBP, and positive correlations with other MHCs. In terms of chemokines, it mainly showed negative correlations with CCL20, CCL26, CCL27, CXCL14, XCL1, and weak correlations with CCL1, CCL15, CCL16, CCL23, CCL25, CCL28, CXCL1, CXCL10, CXCL11, CXCL17. Positive correlations were observed with other chemokines. Regarding chemokine receptors, HGF exhibited weak correlations with CCR3 and CCR9 and positive correlations with other chemokine receptors.And HGF expression showed a positive correlation with GZMA and PRF1 (cor = 0.239, cor = 0.283, P<0.05), representing levels of cytotoxic T lymphocytes (CTL) in LUSC (Fig. 7A).

Single cell analysis of HGF in HNSC and TGF in LUSC

Sequencing the transcriptome of a single cell is a crucial technique used for analyzing the functions of candidate molecules in single cells. HGF expression profles were demonstrated at single-cell levels from HNSC by T-SNE diagram (Fig. 8A-B). HGF exhibited high expression in monocytes/macrophages, myofibroblasts, fibroblasts, endothelial cells, and malignant cells, with minimal expression in immune cells except for a small population of plasma cells. Additionally, the expression of the HGF receptor MET was predominantly concentrated in muscle cells, malignant cells, and endothelial cells (Fig. 8C). TNF profles were demonstrated at single-cell levels from NSCLC by T-SNE diagram (Fig. 8D-E).TNF was expressed in nearly all immune cells, while its expression in malignant cells, endothelial cells, and epithelial cells was relatively low.

Fig. 9
figure 9

Depicts the functional analysis of the interacting proteins of HGF and TNF: (A-B). The top 23 interacting proteins of HGF and TNF. (C-D). GO analysis results of the interacting proteins of HGF. (E-F). GO analysis results of the interacting prote

The function analysis of HGF and TNF

Figure 9A-B illustrates the top 50 interacting proteins with HGF and TNF. The main functions of HGF-interacting proteins include proliferation, survival, protein synthesis, and glycogenesis, while TNF-interacting proteins primarily function in inflammation and proliferation. Furthermore, GO analysis (Biological function) revealed that TNF predominantly participates in NF-κB signaling, cytokine-mediated signaling, and cellular response to tumor necrosis factor, among others (Fig. 9E-F). HGF, on the other hand, is primarily involved in peptidyl-tyrosine phosphorylation, positive regulation of MAPK cascade, regulation of ERK1 and ERK2 cascade, positive regulation of kinase activity, phosphatidylinositol 3-kinase signaling, response to fibroblast growth factor stimulus, and other processes (Fig. 9C-D).

Discussion

In our comprehensive study involving pan-cancer eQTL analysis of 109 cytokines and chemokines, we identified IL4 as having a repressive effect on gastric cancer development. Similarly, TNF demonstrated an inverse relationship with hepatocellular carcinoma, while FLT3LG hinted at a negative link with lung adenocarcinoma. On the other hand, CD40LG and IL18 were suggestively positively associated with lung adenocarcinoma. In the case of LUSC, both CD70 and TNF showed suggestive negative associations, whereas CSF1 was positively associated. Further findings included negative associations of XCL1 with thyroid cancer, HGF with head and neck cancer, and IL2RA with cervical cancer, while LTB showed a positive association with cervical cancer. FLT3LG also negatively correlated with malignant skin neoplasms. These initial observations were validated using pQTL analysis, which confirmed that HGF was associated with a reduced risk of head and neck cancer, and higher TNF levels corresponded to a lower risk of squamous cell lung cancer, aligning with our eQTL findings. Notably, HGF was predominantly expressed in monocytes/macrophages, myofibroblasts, fibroblasts, endothelial cells, and malignant cells, with minimal presence in immune cells apart from a select group of plasma cells in HNSC. Furthermore, the HGF receptor MET was primarily observed in muscle cells, malignant cells, and endothelial cells in HNSC cases. TNF expression was almost universal across immune cells in NSCLC. Survival analysis proposed that HGF may enhance the prognosis for HNSC, a result that contrasts with expectations based on its interaction with the c-MET ligand. This trend remained consistent across eQTL and pQTL analyses. Given the observed prognosis enhancement did not proceed via the HGF-cMet pathway, subsequent KEGG/GO analysis indicated the involvement of HGF-related proteins in regulating cell proliferation, survival, protein synthesis, and glycogen synthesis, among other processes. However, further examination revealed no significant impact of TNF on the survival rate among LUSC patients, indicating TNF may detrimentally affect survival in macrophage, NK cell, or eosinophil-enriched LUSC. Primarily secreted by inflammatory cells, TNF’s influence on prognosis could also stem from its action on inflammatory cells, such as macrophages and NK cells, necessitating deeper investigation into its role and mechanisms.

In summary, our research revealed that HGF has a significant positive role in the development and prognosis of HNSC, as evidenced by both eQTL and pQTL analyses. Notably, while HGF is connected to various elements of the HNSC immune microenvironment, survival analyses indicated that its beneficial impact on prognosis might not operate via the immune cell components or the HGF-Met signaling pathway. Conversely, TNF exhibited a substantial positive influence at the onset of LUSC across eQTL and pQTL levels. Despite its significant role in modifying the immune microenvironment of LUSC, TNF was associated with a detrimental effect on patient survival when CD8 + T cells, eosinophils, NK cells, macrophages, and regulatory T cells were abundant. This suggests that TNF may influence LUSC prognosis through its interactions with these specific cellular components. And the limitation of this study is that it is primarily based on database data and has not yet been validated through experiments and population studies. In the future, we will conduct further research to complete the validation.

Data availability

The data that support the findings of this study are available in public database Open GWAS.

References

  1. Shinko D, Diakos CI, Clarke SJ, Charles KA. Cancer-related systemic inflammation: the challenges and Therapeutic opportunities for Personalized Medicine. Clin Pharmacol Ther. 2017;102(4):599–610. https://doi.org/10.1002/cpt.789.

    Article  PubMed  Google Scholar 

  2. Lin W, Shen P, Song Y, Huang Y, Tu S. Reactive oxygen species in autoimmune cells: function, differentiation, and Metabolism.Front immunol. 2021 Feb 25:12635021. https://doi.org/10.3389/fimmu.2021.635021

  3. Bi Q, Wu JY, Qiu XM, Zhang JD, Sun ZJ, Wang W. Tumor-Associated inflammation: the Tumor-promoting immunity in the early stages of Tumorigenesis. J Immunol Res 2022 Jun 13:20223128933. https://doi.org/10.1155/2022/3128933

  4. Schreiber RD, Old LJ, Smyth MJ. Cancer immunoediting: integrating immunity’s roles in cancer suppression and promotion. Science. 2011;331(6024):1565–70. https://doi.org/10.1126/science.1203486.

  5. Watson J, Salisbury C, Banks J, Whiting P, Hamilton W. Predictive value of inflammatory markers for cancer diagnosis in primary care: a prospective cohort study using electronic health records. Br J Cancer. 2019;120(11):1045–51. https://doi.org/10.1038/s41416-019-0458-x.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Trabert B, Pinto L, Hartge P, Kemp T, Black A, Sherman ME, Brinton LA, Pfeiffer RM, Shiels MS, Chaturvedi AK, Hildesheim A. Wentzensen N.Pre-diagnostic serum levels of inflammation markers and risk of ovarian cancer in the prostate, lung, colorectal and ovarian cancer (PLCO). Screen Trial Gynecol Oncol. 2014;135(2):297–304. https://doi.org/10.1016/j.ygyno.2014.08.025.

    Article  CAS  PubMed  Google Scholar 

  7. Shiels MS, Katki HA, Hildesheim A, Pfeiffer RM, Engels EA, Williams M, Kemp TJ, Caporaso NE, Pinto LA, Chaturvedi AK. Circulating inflammation markers, risk of Lung Cancer, and Utility for Risk Stratification. J Natl Cancer Inst. 2015;107(10):djv199. https://doi.org/10.1093/jnci/djv199.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Burgess S, Dudbridge F, Thompson SG. Combining information on multiple instrumental variables in mendelian randomization: comparison of allele score and summarized data methods. Stat Med. 2016;35(11):1880–906. https://doi.org/10.1002/sim.6835.

    Article  PubMed  Google Scholar 

  9. Burgess S, Thompson SG. Erratum to: interpreting findings from mendelian randomization using the MR-Egger method. Eur J Epidemiol. 2017;32(5):391–2. https://doi.org/10.1007/s10654-017-0276-5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Burgess S, Thompson SG. Interpreting findings from mendelian randomization using the MR-Egger method. Eur J Epidemiol. 2017;32:377–89. https://doi.org/10.1007/s10654-017-0255-x.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Bouras E, Karhunen V, Gill D, Huang J, Haycock PC, Gunter MJ, et al. Circulating inflammatory cytokines and risk of five cancers: a mendelian randomization analysis. BMC Med. 2022;20(1):3. https://doi.org/10.1186/s12916-021-02193-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Folkersen L, Gustafsson S, Wang Q, Hansen DH, Hedman AK, Schork A et al. Genomic and drug target evaluation of 90 cardiovascular proteins in 30,931 individuals.Nat Metab (2020) 2(10):1135–48. https://doi.org/10.1038/s42255-020-00287-2

  13. Sun BB, Maranville JC, Peters JE, Stacey D, Staley JR, Blackshaw J, et al. Genomic atlas of the human plasma proteome. Nature. 2018;558(7708):73–9. https://doi.org/10.1038/s41586-018-0175-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Karhunen V, Gill D, Malik R, Ponsford M, Ahola-Olli A, Papadopoulou A et al. Genetic study of circulating cytokines offers insight into the determinants,cascades and effects of systemic inflammation. medRxiv (2020). https://doi.org/10.1101/2020.10.26.20219477

  15. Thorolfsdottir RB, Jonsdottir AB, Sveinbjornsson G, Aegisdottir HM et al. Variants at the interleukin 1 gene locus and Pericarditis.JAMA cardiol. 2024;9(2):165–72. https://doi.org/10.1001/jamacardio.2023.4820

  16. Folkersen L, Fauman E, Sabater-Lleal M, Strawbridge RJ, et al. Mapping of 79 loci for 83 plasma protein biomarkers in cardiovascular disease. PLoS Genet. 2017;13(4):e1006706. https://doi.org/10.1371/journal.pgen.1006706. eCollection 2017 Apr.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Suhre K, Arnold M, Bhagwat AM, Cotton RJ, Engelke R et al. Connecting genetic risk to disease end points through the human blood plasma proteome.Nat Commun. 2017 Feb 27:814357. https://doi.org/10.1038/ncomms14357

  18. James J, Lee R, Wedow A, Okbay E, Kong O, Maghzian et al. Gene discovery and polygenic prediction from a genome-wide association study of educational attainment in 1.1 million individuals.Nat Genet. 2018;50(8):1112–21. https://doi.org/10.1038/s41588-018-0147-3

  19. Oselli C, Chaffin MD, Weng LC, Aeschbacher S, Ahlberg G, Albert CM, et al. Multi-ethnic genome-wide association study for atrial fibrillation. Nat Genet. 2018;50(9):1225–33. https://doi.org/10.1038/s41588-018-0133-9.

    Article  CAS  Google Scholar 

  20. Yao C, Chen G, Song C, Keefe J, Mendelson M, et al. Genome-wide mapping of plasma protein QTLs identifies putatively causal genes and pathways for cardiovascular disease. Nat Commun. 2018;9(1):3268. https://doi.org/10.1038/s41467-018-05512-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Consortium GT. The GTEx consortium atlas of genetic regulatory effects across human tissues. Science. 2020;369(6509):1318–30. https://doi.org/10.1126/science.aaz1776.

    Article  CAS  Google Scholar 

  22. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B(Methodological). 1995;57(1):289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x.

    Article  Google Scholar 

  23. Staley JR, Blackshaw J, Kamat MA, Ellis S, Surendran P, Sun BB, et al. PhenoScanner: a database of human genotype-phenotype associations.Bioinformatics (Oxford. England). 2016;32(20):3207–9.

    CAS  Google Scholar 

  24. Burgess S, Butterworth A, Thompson SG. Mendelian randomization analysis with multiple genetic variants using summarized data. Genet Epidemiol. 2013;37(7):658–65. https://doi.org/10.1002/gepi.21758.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Burgess S, Thompson SG. Interpreting findings from mendelian randomization using the MR-egger method. Eur J Epidemiol. 2017;32(5):377–89. https://doi.org/10.1007/s10654-017-0255-x.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Bowden J, Davey Smith G, Haycock PC, Burgess S. Consistent estimation in mendelian randomization with some invalid instruments using a weighted median estimator. Genet Epidemiol. 2016;40(4):304–14. https://doi.org/10.1002/gepi.21965.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Burgess S, Foley CN, Allara E, Staley JR, Howson JMM. A robust and efficient method for mendelian randomization with hundreds of genetic variants. Nat Commun. 2020;11(1):376. https://doi.org/10.1038/s41467-019-14156-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Verbanck M, Chen CY, Neale B, Do R. Detection of widespread horizontal pleiotropy in causal relationships inferred from mendelian randomization between complex traits and diseases. Nat Genet. 2018;50(5):693–8. https://doi.org/10.1038/s41588-018-0099-7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Burgess S, Thompson SG, Collaboration CCG. Avoiding bias from weak instruments in mendelian randomization studies. Int J Epidemiol. 2011;40(3):755–64. https://doi.org/10.1093/ije/dyr036.

    Article  PubMed  Google Scholar 

  30. Brion MJ, Shakhbazov K, Visscher PM. Calculating statistical power in mendelian randomization studies. Int J Epidemiol. 2013;42(5):1497–501. https://doi.org/10.1093/ije/dyt179.

    Article  PubMed  Google Scholar 

  31. Sun L, Craiu RV, Paterson AD, Bull SB. Stratified false discovery control for large-scale hypothesis testing with application to genome-wide association studies. Genet Epidemiol. 2006;30(6):519–30. https://doi.org/10.1002/gepi.20164.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This study was supported by the National Natural Science Foundation of China (No.82100539 to H.L).

Author information

Authors and Affiliations

Authors

Contributions

H.L. Writing the first draft of the manuscript, statistical analysis, data curation, writing review, and editing. B.C. Conceptualization, designed the study, supervised, editing, and reviewed, revised the manuscript. All authors approved the final version.

Corresponding author

Correspondence to Bangwei Cao.

Ethics declarations

Ethics approval and consent to participate

The present Mendelian randomization analysis was based on summary data from previous studies that had gained written informed consent and ethics approval. No ethical permit is required for the secondary analysis of summary data.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lin, H., Cao, B. Integration of QTL and comprehensive analysis in the circulating inflammatory cytokines for pan-cancer. BMC Cancer 24, 1007 (2024). https://doi.org/10.1186/s12885-024-12726-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12885-024-12726-4

Keywords