Abstract Objective This study aimed to explore potential associations between polyunsaturated fatty acid (PUFA) metabolism and laryngeal cancer (LC) using a combination of Mendelian randomization (MR) and integrative bioinformatics, while investigating its potential immunogenic implications. Methods Transcriptomic data from the GEO database were integrated to identify differentially expressed genes (DEGs) in LC. Gene Set Enrichment Analysis (GSEA) was performed in addition to Gene Ontology and KEGG enrichment analyses to characterize biological functions. MR analysis was employed to assess the potential causal link between PUFA metabolism and LC risk. PUFA-associated genes were screened for druggability using the DGIdb database. Furthermore, an lncRNA–miRNA–mRNA (ceRNA) regulatory network was constructed, and immune cell infiltration in LC was profiled using CIBERSORT with LM22 signature matrix to uncover immunoregulatory associations with PUFA pathway genes. Results A total of 730 DEGs were identified from a combined cohort of 30 LC and 31 normal samples. GSEA and enrichment analysis revealed significant involvement in fatty acid metabolism and immune-related pathways. MR analysis demonstrated a statistically significant but modest association between PUFA metabolism and LC risk (OR = 0.9996, 95% CI = 0.9992–0.9999, p = 0.020). A network of key PUFA pathway genes (e.g., ADH1B, ADH1A, ALDH7A1, CYP2J2, PTGIS) was linked with multiple drug targets. A focused sub-network of regulatory interactions among lncRNAs, miRNAs, and mRNAs was elucidated. Immune infiltration analysis revealed distinct correlations between pathway gene expression and immune cell populations, particularly showing PTGIS positively correlated with activated mast cells and negatively with naive B cells. Conclusions Our findings suggest a potential mechanistic link between PUFA metabolism and LC susceptibility, generating hypotheses that PUFA-associated pathways may modulate immune microenvironment characteristics. This computational study provides preliminary insights into potential immunometabolic therapeutic strategies in laryngeal cancer that warrant experimental validation. Supplementary Information The online version contains supplementary material available at 10.1007/s12672-025-03578-w. Keywords: Laryngeal cancer, Polyunsaturated fatty acids, Mendelian randomization, Immune infiltration, Bioinformatics, ceRNA network Introduction Two kinds of laryngeal carcinoma (LC) exist: primary and secondary. Primary LC refers to the primary site of laryngeal tumors, and the most frequently occurring type is squamous cell carcinoma. Secondary LC refers to the metastasis of a tumor from another location to the larynx, which is relatively sparse [[38]1–[39]3]. LC is marked by a hoarse voice, cough, difficulty in breathing, lymph node metastasis, dysphagia, etc [[40]4]. Early detection and treatment are essential to reduce the risk of LC. Moreover, it can increase the patient’s chances of survival after surgery [[41]5]. But it is not known exactly what causes LC, and it could be due to a variety of factors, including drinking, smoking, sex hormone levels, genetics, HPV infection, and environmental changes [[42]6]. Numerous studies demonstrate the critical function gene regulation plays in LC. The body goes through several biochemical events during the fat metabolism process, which refers to the process of digesting, absorbing, synthesizing, and decomposing body fat with the help of different types of related enzymes, and processing it into substances needed by the body to ensure the normal functioning of physiological functions. These are of great significance for life activities [[43]7]. Lipid metabolism, in particular fatty acid synthesis, is an important process in which nutrition is transformed into metabolic intermediates for membrane biosynthesis, energy storage and signaling molecules. The changes of lipid metabolism influence the migration of tumor cells, induce angiogenesis, co-metabolism, monitoring drug resistance and the ability to evade immunity, and are the important metabolic phenotypes of cancer cells [[44]8, [45]9]. However, the precise nature of the relationship between fatty acid metabolism and LC risk remains to be fully elucidated. In this research, we gathered and examined multiple data sets from the GEO database to identify DEGs in LC. We performed both traditional enrichment analysis and GSEA to identify important molecular functions and metabolic pathways of LC. Using Mendelian randomization (MR), we explored the potential causal relationship between the metabolism of polyunsaturated fatty acids (PUFA) and the risk of LC. Furthermore, we identified genes and potential drugs that target PUFAs’ metabolic pathways and constructed an lncRNA -- miRNA -- mRNA regulatory network in LC using the DGIdb database. Finally, we investigated the immune cell infiltration in LC and the correlation between the pathway genes and the immune cells. Methods Data source and collation The GEO database is an extensive and comprehensive repository established by one of the world’s leading research centers in the field of biotechnology - the North American National Center for Biotechnology Information. It collects high-throughput gene expression data submitted by research institutions worldwide, covering multiple fields such as tumors, non tumors, NGS, chips [[46]10]. In the course of our research, we utilized specific datasets from the GEO database, namely [47]GSE201777, [48]GSE51985, and [49]GSE117005. These data sets provide us with three distinct transcriptome gene expression matrices. In order to make the data meaningful and useful, we first transformed the probe identifiers from the data set into gene symbols. We divided the samples into two distinct groups based on the associated clinical data: LC group and normal group. These two groups were used as control groups so that we could conduct comparative analysis and evaluate the differences between the two states. This grouping method ensures the accuracy and reliability of the research results, providing a solid foundation for subsequent experimental design. To further enhance our analytical capabilities, we adopted R language software tools and utilized limma and sva packages. Through these advanced data processing tools, we can delve deeper into the hidden patterns and patterns behind the data, providing more precise and powerful support for our research. Using these tools, the three gene expression matrices were successfully combined by us, producing a more sophisticated and cohesive gene expression matrix [[50]11]. Differentially expressed gene (DEG) analysis in LC Using the R packages limma and sva, we first applied batch correction to the combined gene expression matrix in order to address batch effects that were present across different GEO datasets [[51]12]. To identify differentially expressed genes (DEGs) between the LC and normal control groups, we used R software to analyze them in detail. This required using the Wilcoxon test function to compare gene expression levels across the two distinct groups. Based on this in-depth analysis, a comprehensive expression matrix has been generated for each distinct DEG, in conjunction with differential results parameters. This matrix is essential for detailed interpretation of how the individual components are interconnected and contributes to the overall biological function of the DEGs. The limma package, which was utilized in this process, stands as an integral tool in facilitating the generation of such matrices and providing critical insights into the complex molecular interactions that drive these genes. Then, we visualized the DEGs using the heat map and the volcano diagram to make the data more clearly visible [[52]13]. The assessment of statistical significance was conducted with stringent criteria, applying a modified p-value cutoff of less than 0.05, in addition to an actual log fold change (|logFC|) of at least 0.585. This robust analytical method ensures that we can identify DEGs with certainty and significance, which helps to fully understand the difference in gene expression between LC and normal groups. Function and pathway enrichment analysis of DEGs Using Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG), we have gained a deeper understanding of DEGs’ molecular functions and metabolic pathways [[53]14]. At first, we did a complete GO enrichment analysis of the DEGs using the R programming environment, namely the packages clusterProfiler and org.Hs.eg.db. This analysis included Biological Processes (BP), Cellular Components (CC), and Molecular Functions (MF), which provided a variety of perspectives on the roles and interactions of genes. After conducting GO analysis, we performed enrichment analysis using an enhanced KEGG function. This step is a crucial step in our research process, as it can reveal which pathways or networks are significantly activated or inhibited in a given dataset. Through this approach, we can not only gain a deeper understanding of these differences, but also identify potential biological functions and metabolic pathways, providing powerful clues for further experimental design and molecular mechanism exploration. The objective of this step was to identify important biological pathways associated with identified DEGs [[54]15, [55]16]. Additionally, we performed Gene Set Enrichment Analysis (GSEA) using the entire gene expression profile without applying arbitrary p-value cutoffs, which allowed us to detect more subtle but coordinated changes in gene expression patterns across pathways. The results of the enrichment were visualized by means of bubbles, which effectively showed the significance level, which was determined by P-values below 0.05. Not only did these visual representations improve our understanding of the underlying biology, but they also helped to identify critical pathways that could be targeted for future research. Mendelian randomization analysis In order to explore the possible causal relationship between biological pathways and laryngeal cancer susceptibility, a two-sample Mendelian randomisation (MR) method was used [[56]17]. In our recent research effort, we employed single nucleotide polymorphisms, or SNPs, as instrumental variables (IV) in the context of data analysis. These variant sites were identified from the whole genome using high-throughput sequencing technology, aiming to provide strong genetic evidence for genetic pathology research. In addition, to further explore the association between these genetic markers and diseases, we also obtained and utilized an open access large-scale genome-wide association study (GWAS) database ([57]https://gwas.mrcieu.ac.uk/) [[58]18]. The target pathway (ID: met-d-PUFA) and laryngeal cancer data set (ID: ieu-b-4913) were selected for further MR study. The MR study was performed with the “TwoSampleMR” package, and the association between the risk of laryngeal cancer and the metabolism of polyunsaturated fatty acids was assessed using the Inverse Variance Weighting (IVW) technique. In addition, the MR-Egger method was used to analyze the sensitivity to ensure that the results were robust [[59]19]. Analysis of interaction between target pathway-related genes and molecular drugs The Drug-Gene Interaction database (DGIdb) is a groundbreaking resource for researchers and scientists who are dedicated to the understanding of how drugs interact with their target genes. This database, accessible at [60]https://dgidb.genome.wustl.edu/, offers a wealth of genetic data that helps elucidate the molecular interactions between existing pharmaceutical agents and potential candidates [[61]20]. In this study, we first populated the search query with pathway genes. They then limited the medication classes to approved, anti-tumor, and immunotherapies using filtering parameters. The database source’s standard set up, the type of gene, and the type of interaction were kept. Finally, the interactions between the genes and medications were graphically shown using the Cytoscape program [[62]21, [63]22]. Development of a ceRNA (lncRNA–miRNA–mRNA) regulatory network At first, we took four bioinformatics tools to predict miRNAs related to the pathway: miRanda, miRDB, miRWalk, and TargetScan. In each database, a list of possible interactions between the miRNAs of interest is provided. Then, we discovered the miRNAs that were reliably predicted to be targets by comparing the data from all four databases. Only the miRNAs identified by the four tools were chosen as the ultimate candidate [[64]23]. In addition, the spongeScan database was used for miRNA binding prediction lncRNAs, and names were assigned to identified targeted lncRNAs. Finally, Cytoscape software was used to show the lncRNA-miRNA-mRNA regulatory network [[65]24]. Immune cell infiltration and its correlation analysis The R-package e1071 was used to examine the infiltration of 22 immune cell types in data from the LC and Normal groups using the CIBERSORT algorithm with the LM22 signature matrix as reference, with the corresponding degrees of infiltration computed. We utilized powerful software packages such as limma, dplyr, tidyverse, ggplot2, and linkET in R language to analyze in depth the differences in immune cell infiltration. Meanwhile, we also explored the association between these 22 immune cells and 15 pathway genes. Through scientific statistical methods and precise data processing, we can precisely assess the variations in the infiltration of various immune cell types and reveal their complex connections with specific pathway genes [[66]25, [67]26]. Then, the findings were represented visually and considered statistically significant at a P < 0.05. Results Gene expression matrix of LC Upon first consideration, the datasets [68]GSE201777, [69]GSE51985, and finally [70]GSE117005 were meticulously compiled to form a comprehensive collection of information for these three different samples. Each dataset consisted of distinct groups: [71]GSE201777 encompassed 15 liver cancer (LC) specimens while encompassing 16 control or normal tissues as well; [72]GSE51985 included 10 LC samples and 10 controls in their respective sections; and lastly, [73]GSE117005 featured 5 LC samples and 5 NC samples, each with its own set of associated features that required careful scrutiny and analysis. Then, an expression matrix was obtained for each data set containing a variable number of genes (including Matrix_GSE201777.xls, Matrix_GSE51985.xls, and Matrix_GSE117005.xls). The matrices are combined to form a consolidated gene expression matrix (consisting of 30 LC and 31 Normal samples) (see the additional file Matrix_merge.xls), which is a combination of the LC and the normal samples. Identification of DEGs between LC and normal groups Using differential expression analysis, 730 DEGs were found between the LC and Normal groups. Of those DEGs, 320 showed up-regulation and 410 showed downregulation, as shown in Fig. [74]1a. In addition, a heat map of the top 50 DEGs displaying both up- and down-regulation is displayed in Fig. [75]1b. Fig. 1. [76]Fig. 1 [77]Open in a new tab Identification and visualization of differentially expressed genes (DEGs) in laryngeal cancer. (a) Volcano plot showing the distribution of all genes based on their log2 fold change (x-axis) and -log10(adjusted P-value) (y-axis). Red dots represent significantly upregulated genes (n = 320), blue dots represent significantly downregulated genes (n = 410), and gray dots represent non-significant genes. The horizontal dashed line indicates the significance threshold (adjusted P-value = 0.05), and vertical dashed lines indicate the fold change threshold (|log2FC| = 0.585). (b) Hierarchical clustering heatmap of the top 50 DEGs (25 upregulated and 25 downregulated) across all samples. Each row represents a gene, and each column represents a sample. The color scale indicates normalized gene expression levels, with red indicating high expression and blue indicating low expression. The dendrogram shows the clustering of samples, clearly separating laryngeal cancer samples (LC, n = 30) from normal control samples (Normal, n = 31) GO and KEGG enrichment of 730 DEGs in LC The enrichment of “GO” and “KEGG” in bioinformatics has been thoroughly investigated and analyzed. We have gained insight into the nature of differentially expressed genes (DEGs) in biological processes, molecular functions, and metabolic pathways. According to the GO enrichment study, the BPs were linked to the development of the epidermis, nuclear division, extracellular matrix organization, hormone metabolism and prostate development, while CCs was associated with collagen-containing extracellular matrix, apical and basal membrane, kinetochore, cytosol and collagen trimer. The MFs are mainly concentrated in extracellular matrix structure, peptidase regulator and inhibitor activity, heme binding and tetrapyrrole, actin and microtubule activity. Through KEGG enrichment analysis, the pathways identified by DEGs mainly focus on coronavirus COVID-19, human papillomavirus infection, cell cycle, motor proteins, complement and coagulation cascade, relaxin signaling, fatty acid degradation, polyunsaturated fatty acid biosynthesis, and other pathways (Fig. [78]2). Within the fatty acid degradation pathway, we observed that genes such as ACADM and ACOX3 were downregulated, while ADH1B and ADH1A were upregulated, suggesting a complex metabolic reprogramming in LC that aligns with the protective effect suggested by our MR analysis. Fig. 2. [79]Fig. 2 [80]Open in a new tab Functional enrichment analysis of differentially expressed genes in laryngeal cancer. (a) Gene Ontology (GO) enrichment analysis bubble plot. The top 10 enriched terms for each GO category are shown: Biological Process (BP, red), Cellular Component (CC, green), and Molecular Function (MF, blue). The x-axis represents the gene ratio (number of DEGs in the pathway/total number of DEGs), and the y-axis shows the GO terms. Bubble size indicates the number of genes enriched in each term, and color intensity represents the adjusted P-value. (b) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis bubble plot. The top 20 significantly enriched pathways are displayed. The x-axis represents the gene ratio, and the y-axis shows the pathway names. Bubble size indicates the number of genes enriched in each pathway, and color gradient from blue to red represents increasing significance (-log10(adjusted P-value)). Notable pathways include fatty acid degradation, polyunsaturated fatty acid biosynthesis, and multiple immune-related pathways Gene set enrichment analysis (GSEA) results To complement our traditional enrichment analysis, we performed GSEA on the entire gene expression dataset. GSEA revealed significant enrichment of several metabolic pathways without requiring arbitrary cutoffs for differential expression. Notably, the PUFA metabolism pathway showed a normalized enrichment score (NES) of -1.82 (FDR q-value = 0.023), indicating coordinated downregulation of this pathway in LC samples. This finding supports our hypothesis that altered PUFA metabolism plays a role in LC pathogenesis. Additionally, immune-related pathways including “T cell receptor signaling” (NES = -1.75, FDR q-value = 0.031) and “Natural killer cell mediated cytotoxicity” (NES = -1.68, FDR q-value = 0.042) were also significantly downregulated, suggesting an immunosuppressive microenvironment in LC. PUFA metabolism was causally associated with the risk of LC A thorough summary of the SNP features pertaining to PUFA and LC metabolism was given by the supplementary SNP.xls file. No SNPs were found to be weak instruments. The causality of individual genetic variations on LC is shown in Fig. [81]3a and b. Utilizing the IVW method, a statistically significant association was observed between PUFA metabolism activity level and LC risk, with an odds ratio (OR) of 0.9996 (95% confidence interval (CI) = 0.9992 ~ 0.9999, p = 0.020) (Table [82]1). While this OR is very close to 1, suggesting a small effect size, the statistical significance indicates that even small changes in PUFA metabolism may have meaningful biological implications when considered at a population level. This modest effect size is consistent with the complex, multifactorial nature of cancer development. The funnel diagram’s symmetrical features demonstrate that no significant heterogeneity exists (see Fig. [83]3c). The MR Egger regression’s lack of horizontal pleiotropy (p = 0.8291) underscores the validity of our causal inference. To confirm that our findings are reliable, we made an MR analysis of the remaining SNPs after the removal of each individual SNP, as illustrated in Fig. [84]3d. The study’s consistent findings imply that each SNP contributes equally to the observed association, suggesting that there is no single dominant SNP that influences the metabolism of PUFA and LC. Fig. 3. [85]Fig. 3 [86]Open in a new tab Mendelian randomization analysis of the causal relationship between PUFA metabolism and laryngeal cancer risk. (a) Scatter plot of SNP effects on PUFA metabolism (x-axis) versus SNP effects on laryngeal cancer risk (y-axis). Each black dot represents an individual SNP (n = 24). The colored lines represent different MR methods: IVW (light blue), MR-Egger (dark blue), weighted median (green), weighted mode (red), and simple mode (pink). The slope of each line indicates the causal estimate. (b) Forest plot showing the individual causal effects of each SNP on laryngeal cancer risk. SNPs are listed on the y-axis with their corresponding odds ratios (OR) and 95% confidence intervals (CI) displayed on the x-axis. The diamond at the bottom represents the pooled estimate from the IVW analysis (OR = 0.9996, 95% CI: 0.9992–0.9999, P = 0.020). (c) Funnel plot for assessing heterogeneity and potential bias. The x-axis shows the SNP effect on laryngeal cancer risk, and the y-axis shows the inverse of the standard error (1/SE). The vertical line represents the IVW estimate. Symmetrical distribution of points around the line suggests absence of heterogeneity and directional pleiotropy. (d) Leave-one-out sensitivity analysis. Each row shows the MR estimate (with 95% CI) when the corresponding SNP is excluded from the analysis. The red line indicates the overall IVW estimate using all SNPs. Consistent results across all iterations confirm the robustness of the findings Table 1. Statistical outcomes of Mendelian randomization analysis for the association between PUFA metabolism and laryngeal cancer risk Method Number of SNPs OR (95% CI) Beta SE P-value Heterogeneity P-value Pleiotropy P-value Inverse variance weighted (IVW) 24 0.9996 (0.9992–0.9999) -0.0004 0.0002 0.020 0.432 - MR-Egger 24 0.9994 (0.9987–1.0001) -0.0006 0.0003 0.091 0.415 0.829 Weighted median 24 0.9995 (0.9990-1.0000) -0.0005 0.0003 0.048 - - Weighted mode 24 0.9994 (0.9988-1.0000) -0.0006 0.0003 0.074 - - Simple mode 24 0.9993 (0.9985–1.0001) -0.0007 0.0004 0.112 - - [87]Open in a new tab Abbreviations: OR, odds ratio; CI, confidence interval; SE, standard error; SNP, single nucleotide polymorphism; PUFA, polyunsaturated fatty acid. Note: The heterogeneity P-value was calculated using Cochran’s Q test. The pleiotropy P-value was obtained from the MR-Egger intercept test. Dashes indicate that the test is not applicable for the specific method Pathway genes and targeted drugs: a regulatory network Firstly, we selected 15 metabolic pathway genes (SCD5, ACOX3, FADS1, ACOT7, ACADM, ADH1B, ADH1A, ALDH7A1, ADH1C, CYP2J2, EPHX2, ALOXE3, HPGD, PTGDS, PTGIS) as candidate genes (Fig. [88]4a). A complete pathway gene-drug interaction table has been created from the DGIdb database, per the Supplementary File’s description (Drug-gene.xls). Following the construction of the regulatory network with Cytoscape, the targeted medications and the pathway genes were graphically shown. As shown in Fig. [89]4b, those targeted therapeutic drugs are identified by green ellipses, and a red rectangle is used to represent the pathway genes. The dotted lines connecting these nodes indicated the interactions between them. There are eight pathway genes in all (ADH1B, ADH1A, ALDH7A1, ADH1C, CYP2J2, EPHX2, HPGD, and PTGIS), 44 molecular drugs, and 51 interacting relationships among these genes in this drug regulatory network. Of particular interest, Fomepizole, an alcohol dehydrogenase inhibitor, targets both ADH1B and ADH1A. Given our findings that these genes are upregulated in LC and our MR analysis suggesting a protective role for PUFA metabolism, inhibition of these enzymes might restore normal PUFA metabolic flux and potentially provide therapeutic benefit. Fig. 4. [90]Fig. 4 [91]Open in a new tab Expression patterns of PUFA metabolism genes and drug-gene interaction network. (a) Volcano plot highlighting the 15 PUFA metabolism pathway genes among all DEGs. Red dots indicate significantly upregulated PUFA metabolism genes (ADH1B, ADH1A, ADH1C, EPHX2, HPGD), green dots indicate significantly downregulated PUFA metabolism genes (SCD5, ACOX3, FADS1, ACOT7, ACADM, CYP2J2, ALOXE3, PTGDS, PTGIS), and gray dots represent other DEGs. Gene symbols are labeled for the 15 pathway genes. (b) Drug-gene interaction network constructed using the Drug Gene Interaction Database (DGIdb). Red rectangles represent the 8 PUFA metabolism genes with known drug interactions, and green ellipses represent potential therapeutic compounds (n = 44). Edge thickness indicates the strength of evidence for drug-gene interactions. Key interactions include Fomepizole targeting ADH1B/ADH1A, and multiple compounds targeting CYP2J2 and PTGIS. The network visualization was created using Cytoscape v3.9.1 A ceRNA regulatory network based on the lncRNA-miRNA-mRNA pairs of LC A comprehensive dataset consisting of 127 pairs of mRNAs and miRNAs, along with 211 pairs of lncRNAs and miRNAs, was compiled from diverse databases (see Supplementary files mRNA-miRNA.xls and miRNA-lncRNA.xls). To improve interpretability, we focused on creating a sub-network centered on the key PUFA metabolism genes identified as DEGs. The ceRNA (lncRNA-miRNA-mRNA) regulatory network (depicted in Fig. [92]5) was utilized to illustrate the intricate interplay among these RNA molecules. The focused sub-network contains 5 key PUFA metabolism mRNAs (FADS1, ADH1B, PTGIS, CYP2J2, and ALDH7A1), 18 miRNAs, and 52 lncRNAs with 87 interactions. This refined network clearly shows that FADS1 is regulated by multiple miRNAs including miR-524-5p and miR-520d-5p, which are in turn sponged by several lncRNAs such as XIST and MALAT1, suggesting complex post-transcriptional regulation of this key PUFA metabolism enzyme. Fig. 5. [93]Fig. 5 [94]Open in a new tab Competing endogenous RNA (ceRNA) regulatory network centered on key PUFA metabolism genes in laryngeal cancer. This focused sub-network illustrates the post-transcriptional regulation of 5 key PUFA metabolism genes through lncRNA-miRNA-mRNA interactions. Red squares represent mRNAs (FADS1, ADH1B, PTGIS, CYP2J2, ALDH7A1), green ovals represent miRNAs (n = 18), and blue diamonds represent lncRNAs (n = 52). A total of 87 regulatory interactions are shown. The network reveals complex regulatory patterns, particularly for FADS1, which is targeted by miR-524-5p and miR-520d-5p. These miRNAs are in turn sponged by lncRNAs including XIST and MALAT1, suggesting multilevel regulation of PUFA metabolism. Edge width represents the prediction confidence score from integrated analysis of miRanda, miRDB, miRWalk, and TargetScan databases. Network visualization was performed using Cytoscape v3.9.1 with a force-directed layout algorithm Immune cell infiltration and its correlation in LC By examining the immune cells’ infiltration, 22 immune cells (including LC and Normal) were obtained (see Appendix CIBERSORT-Results.xls). As shown in Fig. [95]6, the correlation matrix between different immune cells was obtained. The association between the two types of immune cells is represented by each square in the correlation matrix. The red square shows positive correlation, the blue square shows the negative correlation, the darker the color, the bigger the square, the stronger the correlation. Specifically, we observed strong positive correlations between activated mast cells and memory B cells (r = 0.68, p < 0.001), while naive B cells showed negative correlations with multiple activated immune cell types. In addition, this study accurately described the correlation between pathway genes and immune cells. Visually, a unique color gradient has been established to represent the dynamic interaction between two entities. The green line represents the positive correlation between gene expression and immune cell infiltration, where an increase in pathway gene expression leads to an increase in immune cell abundance. Conversely, the orange line has a negative correlation, suggesting that the higher the expression of the pathway genes, the lower the immune cell infiltration. Of particular biological significance, PTGIS expression showed strong positive correlation with activated mast cells (r = 0.52, p = 0.003) and negative correlation with naive B cells (r = -0.48, p = 0.007). This suggests that PTGIS, which encodes prostacyclin synthase involved in anti-inflammatory prostacyclin production, may promote a pro-tumorigenic immune microenvironment characterized by mast cell activation and B cell dysfunction. In addition, the intensity of the color further indicates the degree of correlation; the deeper the shade, the stronger the correlation. This graphical representation not only provides insights into how gene expression affects immune cell infiltration, but also highlights the complex interdependence between these two important components. Fig. 6. [96]Fig. 6 [97]Open in a new tab Correlation analysis between PUFA metabolism pathway genes and immune cell infiltration in laryngeal cancer. The heatmap displays the Pearson correlation coefficients between the expression levels of 15 PUFA metabolism genes (columns) and the infiltration levels of 22 immune cell types (rows) as determined by CIBERSORT analysis with the LM22 signature matrix. Green colors indicate positive correlations, while orange colors indicate negative correlations. Color intensity reflects the strength of correlation (scale: -0.8 to 0.8). Asterisks denote statistical significance: * P < 0.05, ** P < 0.01, *** P < 0.001. Notable findings include strong positive correlation between PTGIS and activated mast cells (r = 0.52, P = 0.003) and negative correlation with naive B cells (r = -0.48, P = 0.007). The dendrogram shows hierarchical clustering of immune cells based on their correlation patterns with pathway genes. Analysis was performed using R packages including limma, dplyr, and linkET Discussion LC is common in the head and neck tumors, ranking second. The incidence rate is 2.2/100,000, which is 1% ~ 5%. Due to the impact of LC on mortality rates, it has become a serious global health hazard [[98]27–[99]29]. The genetic defects and abnormal expression of RNA in LC play a significant part in its incidence and growth. PUFA, Polyunsaturated fatty acids refer to long-chain fatty acids that contain more than two double bonds (including two of them). The location of the first unsaturated bond in these fatty acids allows for additional classification, mainly into four categories: ω-3, ω-6, ω-7, and ω-9. Among these four categories, fatty acids ω-3 and ω-6 are essential for controlling a number of bodily physiological processes. In addition, there is also a correlation between inflammatory reactions and ω-3 and ω-6 fatty acids, and therefore play an important role in preventing cardiovascular disease [[100]30, [101]31]. They play an important role in combating cancer, regulating immunity, delaying weight loss, aging, and beauty [[102]32, [103]33]. In this study, 730 DEGs were identified in LC and Normal groups. The biological processes, molecular functions, and metabolic pathways of the DEGs were elucidated through GO, KEGG and GSEA enrichment analysis. A potential causal relationship between the risk of LC and PUFA metabolism was suggested by the MR analysis. In addition, a regulatory network was constructed to explain the interaction of the pathway gene with the target drug, and the ceRNA regulatory network of LC. Finally, we examined the immune cells’ infiltration in LC and their relationship to the pathway genes. Our research demonstrated a potential association between PUFA metabolism and LC risk, and identified key genes including SCD5, ACOX3, FADS1, ACOT7, ACADM, ADH1B, ADH1A, ALDH7A1, ADH1C, CYP2J2, EPHX2, ALOXE3, HPGD, PTGDS, and PTGIS, which participate in PUFA metabolism and LC development. Additionally, our research has uncovered a diverse array of potential drug candidates that can be further developed into effective therapies. The comprehensive assessment we conducted seamlessly integrates the wealth of genetic data obtained through genome-wide association studies (GWAS) with detailed transcriptomic analyses. This cross-disciplinary approach allows us to identify novel molecular mechanisms underlying complex diseases and develop targeted treatments that are tailored to specific gene expression profiles. Our findings suggest that the metabolism of PUFA is inversely related to the risk of LC. In particular, there is a corresponding decrease in the risk of LC development as shown by the odds ratio (OR) below 1 and the beta coefficient (B) below 0 as the PUFA metabolic pathway is activated. These results imply a potentially protective causal relationship between PUFA metabolism and LC. Albino L et al. [[104]34] studied fatty acid metabolism in human LC cells Hep2 and found that the alpha-linolenic acid conversion to higher homologues in PUFA metabolism was superior to linoleic acid, but its efficiency was much lower compared to other non-highly undifferentiated human tumor cells. Therefore, understanding the potential regulatory role of specific PUFAs in cell proliferation regulation can further elucidate the molecular mechanisms of LC proliferation and cell death. Fatty acid desaturase 1 (FADS1), a crucial enzyme in the biosynthesis of PUFAs, is capable of catalyzing the dihomogamma-linolenic acid (DGLA) conversion to arachidonic acid (AA) [[105]35]. Zhao R’s research [[106]36] findings indicated that elevated FADS1 expression is observed in patients with LC, and this heightened expression is closely related to advanced clinical stages and unfavorable prognosis in individuals with recurrent LC following chemotherapy. Researchers have found that when FADS1 protein is inhibited, the proliferation rate of LC cells slows down, their invasion ability weakens, and their migration activity is also affected. This finding suggests that FADS1 is an essential regulator of LC cells’ physiological processes. Conversely, when FADS1 levels are elevated, i.e. over-expression, the reverse is true. This discrepancy may be due to different biological mechanisms and signaling pathways, providing important clues for further exploring the specific regulatory mechanisms of FADS1 in LC cells. Through these findings, researchers can better understand the role of FADS1 in LC cell function and could use it as a potential therapeutic target to improve the effectiveness of LC disease treatment. ADH1B (Alcohol Dehydrogenase 1B (Class I), Beta Polypeptide), ADH1A (Alcohol Dehydrogenase 1 (Class I), Alpha Polypeptide), and ADH1C (Alcohol Dehydrogenase 1 C, Class I, Gamma Polypeptide), have the capacity to metabolize a broad range of substrates, such as ethanol, retinol, aliphatic alcohols, hydroxysteroids and lipid peroxidation products. Comprising diverse homo- and heterodimers of alpha, beta, and gamma subunits, these proteins exhibit notable efficacy in oxidation of ethanol and play a major role in ethanol catabolism [[107]37–[108]39]. Huang HY et al. [[109]40] reviewed 68 cases of LC with 10 years of follow-up, and found that the ADH1B polymorphism of rs17033 had an independent relationship with high risk of LC recurrence (HR = 3.325, 95% CI = 1.684–6.566, p = 0.001). Visapaa JP et al. [[110]41] analyzed the ADH1C genotypes in 107 heavy drinkers with upper respiratory tract and digestive tract cancer and 103 age-matched cancer-free drinkers (both groups consumed similar amounts of alcohol) and discovered that the latter group’s ADH1C*1 allele frequency had significantly increased (61.7% v 49.0%, p = 0.011), which may be the molecular mechanism of the latter’s susceptibility to upper respiratory tract and digestive tract cancer. Our research is a unique combination of MR and advanced bioinformatics techniques for analyzing transcriptome data, which serves as a powerful instrument for understanding intricate biological mechanisms. This multifaceted approach is of great significance in the clinical field as it enables us to analyze molecular pathways that may indicate potential diseases or illnesses. However, it is important to acknowledge that this is primarily a hypothesis-generating computational study. Nonetheless, the study is constrained by its scope and the amount of work conducted, necessitating future validation through protein-level analysis and large-scale clinical trials. Furthermore, we did not examine how LC’s PUFA metabolism pathway genes and potential medications affected animal models, which will be a key focus of our forthcoming research endeavors. The title has been revised to more accurately reflect the exploratory nature of this work, and we emphasize that our findings represent preliminary insights that require experimental validation before clinical translation can be considered. Limitations * The study relies on transcriptomic data from only three GEO datasets ([111]GSE201777, [112]GSE51985, and [113]GSE117005) with a relatively small combined sample size (30 LC and 31 Normal samples), which may limit the generalizability of the findings. * Our Mendelian randomization (MR) analysis, while suggesting a causal relationship between PUFA metabolism and LC risk, is based on publicly available GWAS summary statistics that may not capture all genetic variants relevant to both PUFA metabolism and laryngeal cancer. * The study lacks experimental validation of the identified pathway genes and their proposed functional roles in laryngeal cancer development at the protein level. * The drug-gene interaction network is constructed based on in silico prediction tools without experimental validation of the efficacy of the identified drug candidates in laryngeal cancer models. * The ceRNA (lncRNA-miRNA-mRNA) regulatory network is based on computational predictions from multiple databases, but direct binding interactions among the identified RNA molecules have not been experimentally validated in laryngeal cancer. * The immune cell infiltration analysis does not account for the spatial distribution and functional states of immune cells within the tumor microenvironment, limiting the interpretation of immune-metabolic associations. * The study does not include animal models or in vitro experiments to validate the functional significance of the identified PUFA metabolism pathways in laryngeal cancer progression. * The cross-sectional nature of the data used prevents analysis of temporal changes in gene expression patterns during laryngeal cancer progression. * The identified associations between PUFA metabolism genes and immune cell infiltration may be influenced by confounding factors not accounted for in the analysis. * The clinical relevance of the findings, including their potential impact on patient prognosis or treatment response, requires further validation in prospective clinical studies with larger cohorts. Conclusions Our research suggests a potential causal relationship between LC risk and PUFA metabolism, identifying ADH1B, ADH1A, ALDH7A1, ADH1C, CYP2J2, EPHX2, HPGD, and PTGIS as key players in the molecular mechanisms underlying PUFA metabolism in LC development. We have also identified potential therapeutic targets. Furthermore, analysis of the ceRNA regulatory network (lncRNA-miRNA-mRNA) revealed complex interactions among RNA molecules involved in this pathway. Lastly, we investigated immune cell infiltration in LC tissue and its correlation with pathway genes. These computational findings provide a foundation for future experimental studies aimed at validating the role of PUFA metabolism in laryngeal cancer pathogenesis and exploring its potential as a therapeutic target. Supplementary Information Below is the link to the electronic supplementary material. [114]12672_2025_3578_MOESM1_ESM.xls^ (15.8MB, xls) Supplementary Material 1. Combined gene expression matrix from three GEO datasets (GSE201777, GSE51985, GSE117005) after batch correction, containing 30 laryngeal cancer samples and 31 normal control samples. [115]12672_2025_3578_MOESM2_ESM.xlsx^ (21.2KB, xlsx) Supplementary Material 2. Summary of single nucleotide polymorphisms (SNPs) used as instrumental variables in the Mendelian randomization analysis, including SNP identifiers, effect alleles, effect sizes, standard errors, and P-values for associations with both PUFA metabolism and laryngeal cancer risk. [116]12672_2025_3578_MOESM3_ESM.xls^ (13.8KB, xls) Supplementary Material 3. Complete list of drug-gene interactions identified from the DGIdb database for the 15 PUFA metabolism pathway genes, including drug names, interaction types, sources, and evidence scores. [117]12672_2025_3578_MOESM4_ESM.xls^ (3.9KB, xls) Supplementary Material 4. Predicted miRNA-mRNA interactions for PUFA metabolism genes based on consensus predictions from miRanda, miRDB, miRWalk, and TargetScan databases (127 interaction pairs). [118]12672_2025_3578_MOESM5_ESM.xls^ (5.6KB, xls) Supplementary Material 5. Predicted lncRNA-miRNA interactions from the spongeScan database, showing potential lncRNA sponges for miRNAs targeting PUFA metabolism genes (211 interaction pairs). [119]12672_2025_3578_MOESM6_ESM.xls^ (12.4KB, xls) Supplementary Material 6. Complete immune cell infiltration analysis results showing the relative proportions of 22 immune cell types in each sample as determined by CIBERSORT algorithm with LM22 signature matrix. [120]Supplementary Material 7.^ (8.3MB, xls) [121]Supplementary Material 8.^ (2.3MB, xls) [122]Supplementary Material 9.^ (9.4MB, xls) Acknowledgements