
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)13769-3
10.1016/j.heliyon.2024.e37738
e37738
Research Article
Predicting biomarkers in laryngeal squamous cell carcinoma based on the cytokine-cytokine receptor interaction pathway
Chen Qingyong a
Wang Dongqing b
Chen Zhipeng b
Lin Liqiang b
Shao Qiang b
Zhang Han c
Li Peng b
Lv Huaiqing lvhuaiqing@126.com
b⁎
a The Second School of Clinical Medicine of Binzhou Medical University, Yan Tai, China
b Department of Otorhinolaryngology, Linyi People's Hospital, Linyi, China
c No.One Clinical Medicine School of Binzhou Medical University, Bing Zhou, China
⁎ Corresponding author. lvhuaiqing@126.com
10 9 2024
30 9 2024
10 9 2024
10 18 e3773826 5 2024
5 9 2024
9 9 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
Objective

To analyze and validate differential genes in the cytokine-cytokine receptor interaction CCRI pathway in laryngeal squamous cell carcinoma (LSCC) using bioinformatics and Mendelian randomization (MR) to find potential biomarkers for LSCC.

Methods

Five sets of LSCC-related gene chips were downloaded from the GEO database, and four sets of combined datasets were randomly selected as the test set and one set as the validation set to screen for differential genes in the CCRI pathway; two-way Mendelian randomization was performed to analyze the causal relationship between cytokine receptor as the exposure factor and LSCC as the outcome variable; and the causal relationship was analyzed by DGIdb, Miranda, miRDB, miRWalk, TargetScan, spongeScan, and TISIDB databases to analyze the relationship between differential genes and drugs, immune cell infiltration, and mRNA-miNA-lncRNA interactions.

Results

A total of 7 differentially expressed genes CD27, CXCL2, CXCL9, INHBA, IL6, CXCL11, and TNFRSF17 were screened for enrichment in the CCRI signaling pathway; MR analysis showed that the CCRI receptor was a risk factor for LSCC (IVW: OR = 1.629, 95 % CI:1.060–2.504, P = 0.026); Seven differential genes were correlated with drugs, immune cells and mRNA-miNA-lncRNA, respectively; the CCRI differential gene expression analysis in the validation set was consistent with the test set results.

Conclusion

This study provided CCRI differential gene expression by bioinformatics, and MR analysis demonstrated that cytokine receptors are risk factors for LSCC, providing new ideas for the pathogenesis and therapeutic targets of LSCC.

Keywords

Mendelian randomization analysis
Laryngeal squamous cell carcinoma
Cytokine-cytokine receptor interaction pathway
Bioinformatics
Gene-drug interactions
CeRNA regulatory network
Immune infiltration
==== Body
pmc1 Introduction

Laryngeal cancer is one of the most common malignancies in otolaryngology-head and neck surgery, with an incidence rate second only to nasopharyngeal carcinoma and nasal cavity cancer [1]. Laryngeal squamous cell carcinoma (LSCC) is the most common histological subtype of laryngeal cancer, accounting for about 85 % ∼ 95 % of cases [2].It is closely associated with smoking and alcohol consumption, and its pathogenesis is not yet fully elucidated [3].Due to the anatomical complexity and specificity of the laryngeal region, tumor growth is concealed and difficult to detect in early stages, with approximately 60 % of patients diagnosed at an advanced stage of cancer [4].Current treatment options for laryngeal cancer primarily involve multidisciplinary approaches, including surgery, radiation therapy, chemotherapy, biological therapy, and immunotherapy; however, the 5-year survival rate and prognosis remain suboptimal [5].Postoperative patients often suffer from impaired swallowing, respiratory, and phonatory functions, significantly affecting their quality of life [6].Therefore, identifying new potential target genes in laryngeal cancer to explore the pathogenesis and guide treatment is of great importance.

Cytokines are a group of biologically active signaling molecules synthesized and secreted by immune cells and some non-immune cells, including interleukins, tumor necrosis factors, colony-stimulating factors, interferons, chemokines, and growth factors [7]. Cytokines exert their biological effects by binding to specific cytokine receptors. Through cytokine-receptor interactions, they induce cell apoptosis, regulate cell development and differentiation, modulate the immune response, mediate inflammatory reactions, and promote tissue repair [8]. Studies have found that there may be common signaling pathways between cytokine receptor systems and oncogenes that stimulate tissue cells to proliferate, differentiate, and even transform into cancer cells [9]. With advances in understanding the biology of cytokines and receptors, cytokine therapy and anti-cytokine therapy are increasingly being successfully applied in clinical settings [10].Mendelian randomization is a causal inference method that has been widely used in recent years, It takes advantage of the randomness of alleles at meiosis and the irreversibility of genetic variations before the onset of a disease and uses the genetic variation as an instrumental variable to infer the causal relationship between the exposure factor and the outcome of the study (disease), which largely reduces the bias and confounding by confounding factors.With the development of bioinformatics technologies and Mendelian randomization analysis methods, analyzing key genes in tumor development and exploring their causal relationships to unravel the mechanisms of tumor pathogenesis and therapeutic targets has become one of the current research hotspots.

This study utilizes bioinformatics combined with multiple relevant databases to screen for differentially expressed genes in the CCRI signaling pathway. Through Mendelian randomization analysis, the causal relationship between cytokine receptor pathways and LSCC is investigated. Based on the pathway genes, potential drugs, miRNAs, and immune cells related to LSCC are predicted, providing new theoretical references for the pathogenesis and targeted therapy of LSCC.

2 Materials and methods

2.1 Data sources and processing

We conducted a search in the GEO database for "laryngeal squamous cell carcinoma," with selection criteria based on the data type being experimental research and the study type being expression microarray. The samples were human tissue specimens. This search yielded five gene expression datasets: GSE201777, GSE2379, GSE29330, GSE84957, and GSE39400. We processed these datasets using the limma package in R. We randomly selected GSE84957 as the validation set, while the remaining four datasets were combined to form the test set. The test set included 36 normal tissue samples adjacent to carcinoma (control group) and 86 tumor tissue samples (experimental group); the validation set included 9 normal tissue samples adjacent to carcinoma (control group) and 9 tumor tissue samples (experimental group). To avoid affecting the reliability and accuracy of the datasets due to differences in experimental conditions, technical platforms, sample sources, and other factors. Firstly, the normalizeBetweenArrays function [11] in the limma package of R software (version 4.3.2) was used to perform log2 conversion and normalization on the five datasets to ensure consistency between the arrays. Next, the ComBat function in the sva package [12]of the R software was utilized to perform batch correction on the merged data sets. Finally, the statistical method of Principal Component Analysis PCA was utilized visualized and evaluated [13]. See Fig. 1A before batch correction and Fig. 1B after batch correction.Fig. 1 A Before batch correction. B After batch correction.

Fig. 1

2.2 Methods

2.2.1 Differential gene screening

Using R packages such as limma, dplyr, pheatmap, and ggplot2, differential expression analysis is performed on the combined dataset from the test sets. Filtering criteria are set to adjusted P-value (P.adj) < 0.05 and |log2 fold change (FC)| ≥ 1.5 to select differentially expressed genes in the test sets. Heatmaps and volcano plots are generated to visualize these genes, with heatmaps showing the top 50 significantly upregulated and downregulated genes. The volcano plots use log2FC as the x-axis and -log10(adjusted P-value) as the y-axis.

2.2.2 GO and KEGG pathway enrichment analysis

Gene Ontology (GO) enrichment analysis categorizes gene functions into three components: Cellular Component (CC), Molecular Function (MF), and Biological Process (BP), which describe the cellular context, molecular functions, and biological processes the gene products are involved in Ref. [14]. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis systematically examines the metabolic pathways and biological functions of gene products in cells [15]. The R package clusterProfiler [16] is used for conducting GO and KEGG enrichment analyses, and results are visualized using ggplot2. Criteria for GO analysis are P < 0.05; the same criteria apply to KEGG analysis.

2.2.3 Bidirectional two-sample MR analysis of CCRI and laryngeal squamous cell carcinoma

Mendelian randomization (MR) is a causal inference method that has been widely used in recent years, which uses genetic variation as an instrumental variable to infer causal relationships between exposure factors and study outcomes.MR takes advantage of the randomness of alleles at meiosis and the irreversibility of genetic variation before the onset of disease, which allows MR to reduce confounders and reverse interference, compensating for the difficulty of traditional observational studies in determining causal relationships between exposure factors and outcome variables. Currently, the causal relationship between CCRI receptors and laryngeal squamous cell carcinoma has not been well studied. Therefore, the present study evaluated the causal relationship between the two by two-sample bidirectional MR analysis, which provides new ideas for early identification, prediction, and targeted therapy of laryngeal squamous cell carcinogenesis.

2.2.3.1 Data sources

The CCRI pathway in KEGG pathway enrichment analysis was selected as the main study in this study. The data related to CCRI receptor and laryngeal squamous cell carcinoma in this study were obtained from the genome-wide association studies (GWAS) database (https://gwas.mrcieu.ac.uk/) [17]. The exposure variable is cytokine-cytokine receptor interactions, with the outcome variable being laryngeal squamous cell carcinoma. Cytokine-cytokine receptor interaction receptor data from Sun et al. published in a large scale study in Nature Genetics, cytokine-cytokine receptor interaction receptor data comprising approximately 3301 European population patients and 10,534,735 single nucleotide polymorphisms (SNPs) loci for association analysis with the data set 'prot-a-1520' [18]. The laryngeal squamous cell carcinoma dataset, 'finn-b-C3_LARYNX', contains approximately 220,000 European population patients and 16,380,466 single nucleotide polymorphisms (SNPs) loci.

2.2.3.2 MR analysis and sensitivity analysis

The analyses use R software and the TwoSampleMR package [19], employing methods such as Inverse Variance Weighted (IVW), MR-Egger regression, Weighted Median (WM), Simple mode, and Weighted mode to explore the causal relationship between cytokine-cytokine receptor interactions and laryngeal squamous cell carcinoma[20,21,22,23]. Further robustness and stability of results are enhanced by Cochran's Q test, MR-Egger regression test, and MR-Egger intercept to detect heterogeneity and pleiotropy. A leave-one-out approach is used to test the stability of the data by gradually removing individual SNPs and calculating the combined effect of the remaining SNPs. Scatter plots, forest plots, and funnel plots are created for visual assessment of pleiotropy in MR analysis.

2.2.3.3 Reverse MR analysis

We used laryngeal squamous cell carcinoma as an exposure factor, CCRI receptor as an outcome variable, and SNPs associated with laryngeal squamous cell carcinoma (an exposure factor) as instrumental variables for reverse MR causal analysis, and the screening criterion was still based on IVW, and MR-Egger regression, weighted median method, simple model method and weighted model method were used as complementary methods, and the reverse MR results were used for sensitivity analysis by Cochran Q test, MR-Egger regression test, MR-Egger-intercept method and leave-one-out test for sensitivity analysis to further analyze whether laryngeal squamous cell carcinoma affects CCRI receptors.

2.2.3.4 Validation of differential gene expression on the CCRI pathway

The expression of differential genes CD27, CXCL2, CXCL9, INHBA, IL6, CXCL11, and TNFRSF17 on the CCRI pathway in laryngeal squamous cell carcinoma was also visualized with the R language dplyr package, and the screening conditions were corrected P.adj <0.05 and |log2FC| ≥ 1.5.

2.2.4 CCRI pathway gene-drug interaction analysis

The Drug-Gene Interaction database (DGIdb) provides information on known and potential drug-gene interactions, including over 40,000 genes, 10,000 drugs, and more than 100,000 drug-gene interaction relations across 42 drug categories and at least 49 interaction types such as inhibitors, activators, cofactors, ligands, and vaccines [24].This study predicted cytokine-cytokine receptor pathway differential gene-drug interactions with the help of the DGIdb database (https://dgidb.genome.wustl.edu/), and seven differential genes, CD27, CXCL2, CXCL9, INHBA, IL6, CXCL11, and TNFRSF17, were placed in the Genes search box, Preset Filters filtered on Default Approved, Antineoplastic, and Immunotherapies, and Advanced Filters filtered on Source Databases (22 of 22), Gene Categories (43 of 43), Interaction Types (43 of 43), and Genes.), Interaction Types (31 of 31), and the filter results were imported into Cytoscape software to construct gene-drug interaction networks for visualization and analysis.

2.2.5 Prediction of CCRI pathway differential gene ceRNA regulatory network

Differential miRNAs for mRNA binding of cytokine-cytokine receptor interaction pathway genes were first predicted using miRanda, miRDB, miRWalk, and TargetScan target gene prediction databases, and the predicted differential miRNAs were taken as intersections with the mRNAs to construct pairs of mRNA-miRNA interactions. Perl was used to process the data, and the screening criteria were positive miRNAs in all four databases. Secondly, the spongeScan database was used to predict lncRNAs interacting with differential miRNAs, and pairs of miRNA-lncRNA interactions were obtained. Cytoscape software will be integrated to get the ceRNA network of mRNA-miNA-lncRNA interactions, and visualized.

2.2.6 CCRI pathway gene Immune Cell Infiltration Analysis

The CIBERSORT algorithm performs deconvolution analysis based on the principle of linear support vector regression to analyze the relative abundance of various immune cell subtypes in complex mixed cell samples [25]. The algorithm CIBERSORT was used to predict the composition of infiltrating immune cells in pathway differential genes, retain data of P < 0.05, and then analyze the differential infiltration levels of CCRI genes and 22 immune cells in laryngeal tumor tissues (test group) and adjacent normal tissues (control group). The correlation analysis and visualization of immune cell infiltration were performed using the R language tidyverse package, ggplot2 package and linkET package.

2.2.7 Validation of differential genes in the CCRI pathway

Expression of 7 differential genes of the CCRI signaling pathway in laryngeal tumor tissues was verified using the dataset GSE84957 as a validation set.The gene expression profile dataset of GSE84957 was obtained from the GEO database, with 9 cases of paracancerous normal tissue samples (control group) and 9 cases of tumor tissue samples (test group). The expression of 7 differential genes of CCRI signaling pathway in the dataset GSE84957 was analyzed using the R software limma, reshape2 and ggpubr packages and Wilcoxon test, and the expression of the 7 differential genes in the test group and the control group was suggested by * suggesting P < 0.05, ** suggesting P < 0.01, and *** suggesting P < 0.001 to indicate the expression of the 7 differential genes of the test group and the control group significance, and box-and-line diagram visualization.

3 Results

3.1 Differential Gene screening

In the test dataset, differential genes were selected using the filtering criteria of adjusted P-value (P.adj) < 0.05 and |log2 fold change (log2FC)| ≥ 1.5. A total of 105 differential genes were identified, including 74 downregulated genes and 31 upregulated genes. According to |log2FC|, the top 50 up-regulated differential genes and the top 50 down-regulated differential genes were represented by a heat map(Fig. 2A), which shows that the differential gene expression was in opposite directions in the test and control groups, i.e., the differential genes were down-regulated in the control group when they were highly expressed in the test group. The volcano plot with log2FC as the horizontal coordinate and -log10(adj.P.Val) as the vertical coordinate can show the magnitude and confidence of the difference in differential gene expression(Fig. 2B).Fig. 2 2A Heatmap of test set differential genes: red indicates upregulation of gene expression, blue indicates downregulation of gene expression; 2B Differential gene volcano plot of the test set: red indicates up-regulation, blue indicates down-regulation, and gray indicates non-significant difference, the horizontal coordinate Log2FC value is the magnitude of the expression difference of the gene in the two samples, and the larger value represents the larger difference, and the vertical coordinate adj.P.Val value is the corrected P value, and the larger value represents the higher confidence level.

Fig. 2

3.2 GO and KEGG analysis

The results of the GO enrichment analysis show that differentially expressed genes are enriched in 60 biological processes (BP) such as epidermal development, humoral immune response, and organization of external encapsulating structures. They are also enriched in 14 cellular components (CC) including gelatinase granules, tertiary granule membranes, and collagen extracellular matrix, as well as in 17 molecular functions (MF) such as cytokine activity, endopeptidase activity, and serine-type hydrolase activity. Fig. 3A Bar Chart; Fig. 3B Circle Graph; Fig. 3C Bubble Chart.Fig. 3 GO Enrichment Analysis: 3A Bar Chart where the length of each bar represents the number of genes enriched in that category, with the color from blue to red reflecting the smaller P value; 3B Circle Grapht, the outermost circle is the id of GO, the 3 colors represent the 3 categories CC, MF and BP; the 2nd circle is the number of genes distributed in GO entries, the 3rd circle is the number of differentially differentiated genes enriched on GO categories, and the innermost circle reflects the Rich Factor value; 3C Bubble Chart, where the size of each bubble reflects the number of genes enriched in that category, with the color from blue to red reflecting the smaller P value.

Fig. 3

The KEGG pathway enrichment analysis results indicate that the differentially expressed genes are mainly enriched in 26 signaling pathways, including cytokine-cytokine receptor interaction, IL-17 signaling pathway, and Staphylococcus aureus infection. On the cytokine-cytokine receptor interaction pathway, 7 differentially expressed genes are enriched: CD27, CXCL2, CXCL9, INHBA, IL6, CXCL11, and TNFRSF17. Refer to bar Fig. 4A and bubble Fig. 4B.Fig. 4 KEGG Pathway Enrichment Analysis:4A is a bar graph, where the length of each bar represents the number of genes enriched in that signaling pathway. The redder the color, the smaller the P-value. 4B is a bubble chart, where the size of each bubble reflects the proportion of genes enriched in that signaling pathway. The redder the color, the smaller the P-value.

Fig. 4

3.3 Cytokine-cytokine receptor interaction and bidirectional two-sample MR analysis results for laryngeal squamous cell carcinoma

3.3.1 Results of two-way MR analysis and sensitivity analysis

This study utilizes IVW (Inverse Variance Weighted) as the "gold standard" for assessing the causal relationship between cytokine-cytokine receptor interactions and laryngeal squamous cell carcinoma. MR analysis indicates a positive causal association between cytokine-cytokine receptor interactions and laryngeal squamous cell carcinoma, identifying these interactions as risk factors for the development of the carcinoma (IVW: OR = 1.629, 95 % CI: 1.060–2.504, P = 0.026). The reverse MR analysis did not show significant correlation (IVW: P = 0.47 > 0.05), as shown in Fig. 5A and B. The Cochran Q test within IVW (Q = 16.87, P = 0.394 > 0.05) and MR-Egger regression test (Q = 15.96, P = 0.385 > 0.05) indicate no heterogeneity in the study results; the Egger-intercept (P = 0.369 > 0.05) suggests no presence of horizontal pleiotropy. The leave-one-out test analysis demonstrates that the link between cytokine-cytokine receptor interactions and laryngeal squamous cell carcinoma is not dominated by single SNPs. The method of exclusion did not reveal any SNPs that significantly impacted the causal estimate, as shown in Fig. 6A. The funnel plot results display a symmetrical distribution of the included SNPs, and the MR analysis does not show pleiotropy or heterogeneity, as illustrated in Fig. 6B. Scatter plots employing five algorithms (IVW, MR-Egger regression, weighted median method, simple mode, and weighted mode) for the regression of metabolite SNPs demonstrate a significant positive causal relationship between cytokine-cytokine receptor interactions and laryngeal squamous cell carcinoma, with stability in the SNPs included in the study, as shown in Fig. 6C.Fig. 5 A is the forest plot for the forward Mendelian randomization (MR) analysis; Fig. 5B is the forest plot for the reverse MR analysis.

Fig. 5

Fig. 6 A shows the forest plot of the leave-one-out test analysis results for cytokine-cytokine receptor interaction receptors in laryngeal squamous cell carcinoma MR analysis; Fig. 6B is the funnel plot for MR analysis of cytokine-cytokine receptor interactions in laryngeal squamous cell carcinoma; Fig. 6C is the scatter plot for the risk of laryngeal squamous cell carcinoma based on cytokine-cytokine receptor interactions.

Fig. 6

3.3.2 Differential gene expression results on the CCRI pathway

On the CCRI pathway, five differential genes (CXCL2, CXCL9, INHBA, IL6, CXCL11) are upregulated in laryngeal squamous cell carcinoma, while two differential genes (CD27, TNFRSF17) are downregulated. These results are consistent with the MR findings, as shown in Fig. 7.Fig. 7 The figure represents a volcano plot displaying the differential gene expression on the cytokine-cytokine receptor interaction pathway. Red indicates up-regulation, blue indicates down-regulation, and gray indicates no significant difference. The x-axis represents the Log2 fold change (Log2FC) values, which quantify the difference in gene expression between two samples; the larger the value, the greater the difference. The y-axis represents the adjusted P-value (adj.P.Val), with smaller values indicating higher reliability.

Fig. 7

3.4 Pathway gene-drug interactions

The results show that there are 20 drugs that interact with IL6, 5 drugs that interact with CXCL2, 1 drug that interacts with CD27, 1 drug that interacts with INHBA, and 1 drug that interacts with TNFRSF17. These selected drugs are all related to the functional regulation of differential genes, and their discovery provides more options for the future prevention and treatment of laryngeal squamous cell carcinoma. Red elliptical nodes represent key genes, and blue rectangular nodes represent targeted drugs. Fig. 8.Fig. 8 Represents the differential gene and drug interaction network of the CCRI signaling pathway. Red indicates pathway genes, and blue represents drug names. Lines indicate the interactions between genes and drugs.

Fig. 8

3.5 Construction of the ceRNA network

Using miRanda, miRDB, miRWalk, and TargetScan databases for predicting target genes, miRNAs interacting with mRNAs associated with the CCRI pathway were identified. The results were intersected with differentially expressed miRNAs, yielding 50 interacting miRNAs. Interactions with differentially expressed miRNAs were also predicted for lncRNAs using the spongeScan database, and the intersecting results of miRNAs and differentially expressed lncRNAs resulted in 32 interacting lncRNAs. miRNAs and mRNAs not involved in these interactions were removed. This led to the integration of an mRNA-miRNA-lncRNA interaction network within the ceRNA framework. In this network, as shown in Fig. 7, key gene mRNAs are represented by red elliptical nodes, miRNAs by green triangular nodes, and lncRNAs by purple diamond-shaped nodes. The network comprises 48 nodes (5 mRNAs, 11 miRNAs, and 32 lncRNA nodes) and 46 edges, each edge representing interactions between mRNA and miRNA or miRNA and lncRNA. Among them, lncRNA LINC00689 is involved in binding with both hsa-miR-877-3p and hsa-miR-671-5p, while RP11-384K6.6 can bind with miRNAs hsa-miR-185-5p and hsa-miR-214-5p.Fig. 9.Fig. 9 Presents the construction of a ceRNA network. The red nodes represent pathway gene mRNAs, the green nodes represent miRNAs, and the purple diamonds represent lncRNAs. Each line indicates the interaction between mRNA and miRNA, as well as between miRNA and lncRNA.

Fig. 9

3.6 Analysis of immune cell infiltration of pathway genes

Immune infiltration analysis based on 22 immune-related genes showed that the laryngeal tumor pathway gene CD27 was positively correlated with 5 immune cells and negatively correlated with 4 immune cells; the gene CXCL2 was positively correlated with 2 immune cells; the gene CXCL9 was positively correlated with 3 immune cells and negatively correlated with 1; the gene INHBA was positively correlated with 3 immune cells and negatively correlated with 1 immune cells; gene IL6 is positively correlated with 2 immune cells; gene CXCL11 is positively correlated with 3 immune cells and negatively correlated with 2 immune cells; gene TNFRSF1 is positively correlated with 1 immune cell and negatively correlated with 2 immune cells. Horizontal and vertical coordinates represent 22 types of immune cells, the lower left graph represents the correlation between cells, blue reflects the negative correlation between immune cells, red reflects the positive correlation between immune cells, the upper right graph shows the correlation between 7 differential genes of cytokine-cytokine receptor interaction pathway and 22 types of immune cells in laryngeal tumor tissues, and positive correlation between differential genes and immune cells is represented by red line segment, negative correlation is represented by green. Positive correlations between differential genes and immune cells are represented by red lines, negative correlations are represented by green, and irrelevant correlations are represented by gray. The thicker the line and the larger the rectangle, the stronger the correlation.Fig. 10A.Fig. 10 A illustrates the interaction graph between pathway genes and immune cells in throat tumor tissues. The lower left graph displays the associations between cells, where blue indicates a negative correlation among immune cells, and red indicates a positive correlation. The upper right section shows the relationships between pathway genes and immune cells; red lines represent positive correlations, green lines indicate negative correlations, and gray lines signify no correlation. Thicker lines and larger rectangles suggest stronger correlations. Fig. 10B is a box plot, with the x-axis labeled with pathway gene names and the y-axis showing gene expression levels. Red represents the experimental group, and blue represents the control group. An asterisk (*) indicates p < 0.05, two asterisks (**) indicate p < 0.01, and three asterisks (***) indicate p < 0.001.

Fig. 10

3.7 Validation of group differences analysis

To further verify the accuracy of the seven differential genes in the CCRI pathway identified from Test Set 7 in laryngeal tumor tissues, differential genes in the CCRI pathway were analyzed from the validation set GSE84957, using R language to construct boxplot visualizations. The results indicate that the expression levels of the genes CD27, CXCL2, CXCL9, INHBA, IL6, CXCL11, and TNFRSF17 in the validation set GSE84957 are consistent with those in Test Sets GSE201777, GSE2379, GSE29330, and GSE39400 (P < 0.05). The genes CD27 and TNFRSF17 show low expression, while CXCL2, CXCL9, INHBA, IL6, and CXCL11 exhibit high expression in laryngeal tumor tissues. Fig. 9 shows the boxplots, with the experimental group in red and the control group in blue.Fig. 10B.

4 Discussion

Laryngeal cancer accounts for about 20 % of malignant tumors in the head and neck region, making it the second most common malignant tumor of the respiratory tract after lung cancer [26]. It is characterized by high incidence, high mortality, and low cure rates. Laryngeal squamous cell carcinoma (LSCC) is the most common histological type of laryngeal cancer, with the incidence in men approximately ten times higher than in women. With advances in treatment methods, the 5-year survival rate for early-stage laryngeal cancer patients reaches 90 %, while it is only about 30%–40 % for late-stage patients [27]. Therefore, early diagnosis and treatment of LSCC are crucial for improving prognosis. Studies have shown that the development of resistance in late-stage laryngeal cancer patients is a significant cause of recurrence and treatment failure, making it essential to explore highly accurate and precise biomarkers or therapeutic targets [28].

The CCRI signaling pathway plays a critical role in innate and adaptive inflammatory host defense and cell death processes, and it is involved in the composition of the tumor microenvironment and mediates tumor immune responses [29]. Previous research has shown that the CCRI pathway is involved in tumor cell invasion, angiogenesis, and metastasis in various cancers such as cervical, colon, breast, and pancreatic cancers, and has become a sensitive marker for predicting the prognosis of some malignant tumors and monitoring disease progression [30]. This study analyzed laryngeal squamous cell carcinoma gene expression datasets GSE201777, GSE2379, GSE29330, and GSE39400 from the GEO database and validated with dataset GSE84957, identifying seven differentially expressed genes in the CCRI pathway through GO and KEGG pathway enrichment analysis, including upregulation of CXCL2, CXCL9, INHBA, IL6, CXCL11, and downregulation of CD27 and TNFRSF17. Utilizing public data and large-scale genome-wide association study data (GWAS), this study employed five Mendelian randomization methods, including inverse-variance weighted method, MR-Egger regression, weighted median method, simple mode method, and weighted mode method to investigate the causal relationship between cytokine receptors and laryngeal malignancy. To ensure the reliability and stability of the results, data were analyzed using Cochran Q test, Egger-intercept, leave-one-out test, and reverse MR analysis. The findings indicated a positive causal relationship between cytokine receptors and laryngeal malignancy, identifying it as a risk factor (IVW: OR = 1.629, 95 % CI: 1.060–2.504, P = 0.026). The reverse MR showed no significant correlation (IVW: P = 0.47 > 0.05).

This study uses bioinformatics combined with multiple related databases to analyze drug correlations of differentially expressed genes in the CCRI signaling pathway, construct mRNA-miRNA-lncRNA regulatory networks, and examine immune cell infiltration. Previous research has indicated that IL6 can regulate cancer-related markers and multiple signaling pathways, thereby affecting cellular apoptosis, proliferation, angiogenesis, invasiveness, and metastasis, promoting tumor development [31]. Elevated levels of IL6 correlate with the neutrophil-to-lymphocyte ratio (NLR) in tumor tissues, and in some malignancies, NLR has been identified as an independent prognostic factor for overall and progression-free survival [32]. In this study, IL6 is upregulated in laryngeal tumor tissues compared to adjacent normal tissues and correlates with 20 drugs, including RITUXIMAB, ETANERCEPT, and ELSILIMOMAB. IL6 expression positively regulates tumor-infiltrating macrophages M1 and M2.

Increasing evidence has revealed that chemokines play a crucial role in cancer development, progression, and prognosis [33]. The study discovered that CXCL2, by modulating tumor epithelial-stromal interactions, enhances tumor growth and infiltration, with high expression associated with poor prognosis in breast, prostate, bladder, and colorectal cancers [34]. CXCL9, induced by gamma-interferon, regulates the tumor microenvironment to exert anti-tumor effects and is expressed higher in various cancer tissues such as gastric, esophageal, urothelial, and head and neck squamous cell carcinoma, potentially related to longer survival periods. It is already considered a good prognostic marker and therapeutic target in breast cancer [35]. CXCL11, an interferon-induced T-cell alpha chemokine, under normal circumstances, participates in the transport of immune cells and induces cell proliferation and apoptosis. However, in tumor tissues, it binds to CXCR3 to regulate cancer cell proliferation, invasiveness, and metastasis, adversely affecting prognosis [36]. In this study, the chemokines CXCL2, CXCL9, and CXCL11 are upregulated in laryngeal tumor tissues. CXCL2 is associated with five drugs including DEFEROXAMINE and BATIMASTAT, and three lncRNAs (LINC00940, RP1-253P7.1, LINC00689) compete to bind with hsa-miR-877-3p, with CXCL2 also positively regulating monocytes and macrophages M1.No drugs interacting with CXCL9 were identified through the DGIdb database, but 22 LncRNAs competed with CXCL9 on the mRNA-miRNA-lncRNA regulatory network for binding to five miRNAs (hsa-miR-671-5p, hsa-miR-16-1-3p, has-miR-297, has-miR- 92a-1-5p, and hsa-miR-185-5p), and in immune cell infiltration, CXCL9 was able to positively regulate macrophage M1, macrophage M2, and resting mast cells, and negatively regulate monocytes; CXCL11 competed with LL22NC03-2H8.5 for binding to hsa-miR-1322, and with RP11-397O4.1 for binding to hsa-miR-9-5p, CXCL11 is involved in positive regulation of eosinophils, macrophage M1 and macrophage M2, and negative regulation of resting NK cells and regulatory T cells.INHBA, a member of the transforming growth factor-beta superfamily, is involved in the synthesis of inhibitors and activators and plays a significant role in the development and progression of cancer and other diseases [37]. Studies show that INHBA is overexpressed in cervical, gastric, esophageal, colorectal, and lung cancers, participating in tumor development and closely related to clinical prognosis [38]. In this study, INHBA is upregulated in laryngeal cancer tissues and correlates with the drug SOTATERCEPT, also competitively binding with seven lncRNAs to three miRNAs (hsa-miR-205-5p, hsa-miR-509-3-5p, hsa-miR-214-5p), positively regulating plasma cells, macrophages M1, M2, and negatively regulating resting dendritic cells.

CD27 and TNFRSF17 were downregulated among genes with differential CCRI signaling pathways. CD27 belongs to a transmembrane glycoprotein of the tumor necrosis factor receptor superfamily, and binding to its ligand CD70 is able to enhance immune cell function [39]. Some studies have reported that the low expression of CD27 in multiple myeloma may predict rapid disease progression, poor prognosis, and poor response rate to chemotherapeutic drugs [40]. In this study, CD27 was associated with the drug sarilumab, competing with RP3-470B24.5 for binding to hsa-miR-509-3-5p, and was able to positively regulate CD8 + T cells, plasma cells, helper T cells, regulatory T cells and activated NK cells, and negatively regulate memory CD4 + cells, resting NK cells, neutrophils and activated mast cells. TNFRSF17 Is a member of the tumor necrosis factor receptor superfamily 17, also known as CD269 and B cell mature antigen (BCMA), in the proliferation, differentiation, and apoptosis of participant cells in vivo. Song et al. [41]in a study on colon cancer found that TNFRSF17 showed low expression in colon cancer patients, and first confirmed that TNFRSF17 is a tumor suppressor gene in vitro experiments, and provided immune-related prognostic markers of colon cancer, immunotherapy sensitivity markers, and effective targets for the clinical treatment of colon cancer patients. In this study, TNFRSF17 showed interaction with the drug iivilxel, and TNFRSF17 expression was positively correlated with tumor-infiltrating plasma cells and negatively with monocytes and macrophages M1. Finally, the validation of CCRI signaling pathway differential genes using the validation set GSE84957 and the validation results were consistent with the training set.

5 Conclusion

In conclusion, this study, biological information technology, and multiple databases were used to analyze seven differentially expressed genes in the CCRI signaling pathway in LSCC. In the genetic aspect, MR analysis was used to illustrate the causal relationship between cytokine receptors and laryngeal malignancy, which laid the foundation for further exploring the causal relationship between LSCC and cytokines. Gene-drug interaction to explore the potential therapeutic drugs and therapeutic targets of LSCC can provide new ideas and clues for the therapeutic exploration of LSCC. To construct the mRNA-miRNA-lncRNA interaction in the ceRNA network to explore the potential regulation mechanism of LSCC and provide a new reference for future clinical diagnosis and treatment, prognosis, and rehabilitation strategies. Immune cell infiltration analysis of differential genes in cytokine-cytokine receptor interaction pathway in LSCC provides new ideas and biological markers for immune-related therapeutic targets. However, the results of this study are all based on data analysis, and a large number of clinical experimental data still need to be verified in the future.

Funding

There is no funding available.

Consent for publication

All authors approved the submitted version.

Data availability statement

The datasets analyzed in this study are publicly available summary statistics. The transcriptome datasets analyzed during the current study are available in the public NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/). Accession number: GSE201777,GSE2379,GSE29330,GSE84957 and GSE39400.The Mendelian randomization analysis dataset for this study was obtained from the GWAS database (https://gwas.mrcieu.ac.uk/).Drug-Gene Interaction database (DGIdb) at (https://dgidb.genome.wustl.edu/), freely available. ceRNA Regulatory Network, Immune Cell Infiltration Analysis, and Validation were analyzed by R. For further inquiries, please contact the corresponding author, data will be made available on request.

Ethics statement

Ethical oversight and endorsement for this investigation were exempted, given the utilization of data from a public repository. Affirmative consent had been secured from all participants whose data were included in the public database study.

CRediT authorship contribution statement

Qingyong Chen: Writing – original draft, Software, Resources, Methodology, Funding acquisition, Formal analysis, Data curation. Dongqing Wang: Formal analysis, Data curation. Zhipeng Chen: Software, Resources, Methodology. Liqiang Lin: Methodology, Formal analysis, Data curation. Qiang Shao: Visualization, Supervision, Software, Data curation. Han Zhang: Methodology, Funding acquisition, Formal analysis. Peng Li: Methodology, Data curation. Huaiqing Lv: Writing – review & editing, Writing – original draft, Software, Formal analysis, Data curation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors would like to acknowledge the used database that mentioned throughout the manuscript for providing the comprehensive and up-to-date data that played a pivotal role in our data processing and analysis and allowed us to conduct in-depth exploration and generate insightful results.
==== Refs
References

1 Barsouk A. Aluru J.S. Rawla P. Saginala K. Barsouk A. Epidemiology, risk factors, and prevention of head and neck squamous cell carcinoma Med. Sci. 11 2 2023 Jun 13 42 10.3390/medsci11020042 PMID: 37367741; PMCID: PMC10304137
2 Cavaliere M. Bisogno A. Scarpa A. D'Urso A. Marra P. Biomarkers of laryngeal squamous cell carcinoma: a review Ann. Diagn. Pathol. 54 2021 Oct 151787 10.1016/j.anndiagpath.2021.151787 Epub 2021 Jul 2. PMID: 34242969
3 Igissin N. Zatonskikh V. Telmanova Z. Tulebaev R. Moore M. Laryngeal cancer: epidemiology, etiology, and prevention: a narrative review Iran. J. Public Health 52 11 2023 Nov 2248 2259 10.18502/ijph.v52i11.14025 PMID: 38106821; PMCID: PMC10719707 38106821
4 Liu C. Yu Z. Huang S. Zhao Q. Sun Z. Fletcher C. Combined identification of three miRNAs in serum as effective diagnostic biomarkers for HNSCC EBioMedicine 50 2019 Dec 135 143 10.1016/j.ebiom.2019.11.016 Epub 2019 Nov 26. PMID: 31780396; PMCID: PMC6921333 31780396
5 Anderson G. Ebadi M. Vo K. Novak J. Govindarajan A. Amini A. An updated review on head and neck cancer treatment with radiation therapy Cancers 13 19 2021 Sep 30 4912 10.3390/cancers13194912 PMID: 34638398; PMCID: PMC8508236 34638398
6 Riechelmann H. Dejaco D. Steinbichler T.B. Lettenbichler-Haug A. Anegg M. Ganswindt U. Functional outcomes in head and neck cancer patients Cancers 14 9 2022 2135 10.3390/cancers14092135 35565265
7 Bikfalvi A. Billottet C. The CC and CXC chemokines: major regulators of tumor progression and the tumor microenvironment Am J Physiol Cell Physiol 318 3 2020 Mar 1 C542 C554 10.1152/ajpcell.00378.2019 Epub 2020 Jan 8. PMID: 31913695; PMCID: PMC7099520 31913695
8 Xue X. Falcon D.M. The role of immune cells and cytokines in intestinal wound healing Int. J. Mol. Sci. 20 23 2019 Dec 3 6097 10.3390/ijms20236097 PMID: 31816903; PMCID: PMC6929186 31816903
9 Morris R.M. Mortimer T.O. O'Neill K.L. Cytokines: can cancer get the message? Cancers 14 9 2022 Apr 27 2178 10.3390/cancers14092178 PMID: 35565306; PMCID: PMC9103018 35565306
10 Fu Y. Tang R. Zhao X. Engineering cytokines for cancer immunotherapy: a systematic review Front. Immunol. 14 2023 Jul 6 1218082 10.3389/fimmu.2023.1218082 PMID: 37483629; PMCID: PMC10357296
11 Ritchie M.E. Phipson B. Wu D. Hu Y. Law C.W. Shi W. Limma powers differential expression analyses for RNA-sequencing and microarray studies Nucleic Acids Res. 43 7 2015 Apr 20 e47 10.1093/nar/gkv007 Epub 2015 Jan 20. PMID: 25605792; PMCID: PMC4402510
12 Leek J.T. Johnson W.E. Parker H.S. Jaffe A.E. Storey J.D. The sva package for removing batch effects and other unwanted variation in high-throughput experiments Bioinformatics 28 6 2012 Mar 15 882 883 10.1093/bioinformatics/bts034 Epub 2012 Jan 17. PMID: 22257669; PMCID: PMC3307112 22257669
13 Greenacre M. Groenen P.J.F. Hastie T. D'Enza A.L. Markos A. Tuzhilina E. Principal component analysis Nat Rev Methods Primers 2 2022 100 10.1038/s43586-022-00184-w
14 Gene Ontology Consortium Aleksander S.A. Balhoff J. Carbon S. Cherry J.M. Drabkin H.J. The gene Ontology knowledgebase in 2023 Genetics 224 1 2023 May 4 iyad031 10.1093/genetics/iyad031 PMID: 36866529; PMCID: PMC10158837
15 Kanehisa M. Sato Y. Furumichi M. Morishima K. Tanabe M. New approach for understanding genome variations in KEGG Nucleic Acids Res. 47 D1 2019 Jan 8 D590 D595 10.1093/nar/gky962 PMID: 30321428; PMCID: PMC6324070 30321428
16 Yu G. Wang L.G. Han Y. He Q.Y. clusterProfiler: an R package for comparing biological themes among gene clusters OMICS 16 5 2012 May 284 287 10.1089/omi.2011.0118 Epub 2012 Mar 28. PMID: 22455463; PMCID: PMC3339379 22455463
17 Uffelmann E. Huang Q.Q. Munung N.S. Vries J.D. Okada Y. Martin A.R. Genome-wide association studies Nat Rev Methods Primers 1 2021 59 10.1038/s43586-021-00056-9
18 Sun B.B. Maranville J.C. Peters J.E. Stacey D. Staley J.R. Blackshaw J. Genomic atlas of the human plasma proteome Nature 558 7708 2018 Jun 73 79 10.1038/s41586-018-0175-2 Epub 2018 Jun 6. PMID: 29875488; PMCID: PMC6697541 29875488
19 Spiller W. Davies N.M. Palmer T.M. Software application profile: mrrobust—a tool for performing two-sample summary Mendelian randomization analyses 2019 10.1093/ije/dyy195
20 Mounier N. Kutalik Z. Bias correction for inverse variance weighting Mendelian randomization Genet. Epidemiol. 47 4 2023 Jun 314 331 10.1002/gepi.22522 Epub 2023 Apr 10. PMID: 37036286 37036286
21 Bowden J. Davey Smith G. Burgess S. Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression Int. J. Epidemiol. 44 2 2015 Apr 512 525 10.1093/ije/dyv080 Epub 2015 Jun 6. PMID: 26050253; PMCID: PMC4469799 26050253
22 Bowden J. Davey Smith G. Haycock P.C. Burgess S. Consistent estimation in mendelian randomization with some invalid instruments using a weighted median estimator Genet. Epidemiol. 40 4 2016 May 304 314 10.1002/gepi.21965 Epub 2016 Apr 7. PMID: 27061298; PMCID: PMC4849733 27061298
23 Hartwig F.P. Davey Smith G. Bowden J. Robust inference in summary data Mendelian randomization via the zero modal pleiotropy assumption Int. J. Epidemiol. 46 6 2017 Dec 1 1985 1998 10.1093/ije/dyx102 PMID: 29040600; PMCID: PMC5837715 29040600
24 Cannon M. Stevenson J. Stahl K. Basu R. Coffman A. Kiwala S. DGIdb 5.0: rebuilding the drug-gene interaction database for precision medicine and drug discovery platforms Nucleic Acids Res. 52 D1 2024 Jan 5 D1227 D1235 10.1093/nar/gkad1040 PMID: 37953380; PMCID: PMC10767982 37953380
25 Guan M. Jiao Y. Zhou L. Immune infiltration analysis with the CIBERSORT method in lung cancer Dis. Markers 2022 2022 Mar 18 3186427 10.1155/2022/3186427 PMID: 35340416; PMCID: PMC8956442
26 Ralli M. Grasso M. Gilardi A. Ceccanti M. Messina M.P. Tirassa P. The role of cytokines in head and neck squamous cell carcinoma: a review Clin. Ter. 171 3 2020 May-Jun e268 e274 10.7417/CT.2020.2225 PMID: 32323717 32323717
27 Ghiani L. Chiocca S. High risk-human papillomavirus in HNSCC: present and future challenges for epigenetic therapies Int. J. Mol. Sci. 23 7 2022 Mar 23 3483 10.3390/ijms23073483 PMID: 35408843; PMCID: PMC8998945 35408843
28 Zhang H. Chen Z. Huang Q. Guo Y. Wang M. Wu C. Preliminary study using a small plasma extracellular vesicle miRNA panel as a potential biomarker for early diagnosis and prognosis in laryngeal cancer Cell. Oncol. 46 4 2023 Aug 1015 1030 10.1007/s13402-023-00792-y Epub 2023 Mar 25. PMID: 36964893
29 Smeda M. Stojak M. Przyborowski K. Sternak M. Suraj-Prazmowska J. Kus K. Direct thrombin inhibitor dabigatran compromises pulmonary endothelial integrity in a murine model of breast cancer metastasis to the lungs; the role of platelets and inflammation-associated haemostasis Front. Pharmacol. 13 2022 Feb 28 834472 10.3389/fphar.2022.834472 PMID: 35295330; PMCID: PMC8918823
30 Khurshid Z. Zafar M.S. Khan R.S. Najeeb S. Slowey P.D. Rehman I.U. Role of salivary biomarkers in oral cancer detection Adv. Clin. Chem. 86 2018 23 70 10.1016/bs.acc.2018.05.002 Epub 2018 Jul 23. PMID: 30144841 30144841
31 Hirano, T. IL-6 in inflammation, autoimmunity and cancer. International immunology, 33(3), 127-148 10.1093/intimm/dxaa078.
32 Hou Y. Li X. Yang Y. Shi H. Wang S. Gao M. Serum cytokines and neutrophil-to-lymphocyte ratio as predictive biomarkers of benefit from PD-1 inhibitors in gastric cancer Front. Immunol. 14 2023 Oct 31 1274431 10.3389/fimmu.2023.1274431 PMID: 38022654; PMCID: PMC10643875
33 Ding J. Xu K. Zhang J. Lin B. Wang Y. Yin S. Xie H. Zhou L. Zheng S. Overexpression of CXCL2 inhibits cell proliferation and promotes apoptosis in hepatocellular carcinoma BMB Rep 51 12 2018 Dec 630 635 10.5483/BMBRep.2018.51.12.140 PMID: 30293547; PMCID: PMC6330937 30293547
34 Korbecki J. Bosiacki M. Barczak K. Łagocka R. Brodowska A. Chlubek D. Baranowska-Bosiacka I. Involvement in tumorigenesis and clinical significance of CXCL1 in reproductive cancers: breast cancer, cervical cancer, endometrial cancer, ovarian cancer and prostate cancer Int. J. Mol. Sci. 24 8 2023 Apr 14 7262 10.3390/ijms24087262 PMID: 37108425; PMCID: PMC10139049 37108425
35 Pein M. Insua-Rodríguez J. Hongu T. Riedel A. Meier J. Wiedmann L. Metastasis-initiating cells induce and exploit a fibroblast niche to fuel malignant colonization of the lungs Nat. Commun. 11 1 2020 Mar 20 1494 10.1038/s41467-020-15188-x PMID: 32198421; PMCID: PMC7083860 32198421
36 Yuan X. Zhou Q. Zhang F. Zheng W. Liu H. Chen A. Identification of immunity- and ferroptosis-related genes for predicting the prognosis of serous ovarian cancer Gene 838 2022 Sep 5 146701 10.1016/j.gene.2022.146701 Epub 2022 Jun 28. PMID: 35777713
37 Xiao Q. Xiao J. Liu J. Liu J. Shu G. Yin G. Metformin suppresses the growth of colorectal cancer by targeting INHBA to inhibit TGF-β/PI3K/AKT signaling transduction Cell Death Dis. 13 3 2022 Mar 2 202 10.1038/s41419-022-04649-4 PMID: 35236827; PMCID: PMC8891354 35236827
38 Zhang Y. Yan S. Li Y. Zhang J. Luo Y. Li P. Inhibin βA is an independent prognostic factor that promotes invasion via Hippo signaling in non-small cell lung cancer Mol. Med. Rep. 24 5 2021 Nov 789 10.3892/mmr.2021.12429 Epub 2021 Sep 10. PMID: 34505633; PMCID: PMC8441965 34505633
39 Wajant H. Therapeutic targeting of CD70 and CD27 Expert Opin. Ther. Targets 20 8 2016 Aug 959 973 10.1517/14728222.2016.1158812 Epub 2016 Mar 15. PMID: 26914723 26914723
40 Chu B. Bao L. Wang Y. Lu M. Shi L. Gao S. CD27 antigen negative expression indicates poor prognosis in newly diagnosed multiple myeloma Clin Immunol 213 2020 Apr 108363 10.1016/j.clim.2020.108363 Epub 2020 Feb 29. PMID: 32120013
41 Song Y. Zhang Z. Zhang B. Zhang W. CD8+ T cell-associated genes MS4A1 and TNFRSF17 are prognostic markers and inhibit the progression of colon cancer Front. Oncol. 12 2022 Sep 20 941208 10.3389/fonc.2022.941208 PMID: 36203424; PMCID: PMC9530608
