Integrative Analysis of Acupuncture Targets and Immune Genes in Diabetes, Diabetic Peripheral Neuropathy, and Adjunct Therapy of Cancer

Introduction

Diabetes (DM) is a prevalent endocrine and metabolic disease characterized by chronic hyperglycemia and the systemic metabolic disorder syndrome.¹ In more than 60% of diabetic individuals, diabetic neuropathy is one of the most common chronic complications of diabetes.² About 75% of neuropathy is diabetic peripheral neuropathy (DPN), which is characterized by aggravated nerve damage, decreased muscle strength, muscle atrophy and weakened tendon reflexes.³ DPN greatly increases the risk of diabetic foot and other injuries in patients with diabetes, which causes the majority of diabetic amputees or handicapped people.⁴ Epidemiological studies have consistently demonstrated a higher cancer risk among diabetic patients, suggesting that hyperglycemia may play a pivotal role in tumorigenesis through mechanisms such as increased insulin levels, inflammation, and oxidative stress.⁵ Recent study has highlighted the efficacy of acupuncture in improving glycemic control and reducing neuropathic pain, suggesting its role as an adjunctive therapy in the management of DPN.⁶ Furthermore, acupuncture has been employed to manage chemotherapy-induced nausea and vomiting, fatigue, and pain in cancer patients.⁷

Acupuncture has gained significant attention as a complementary and alternative therapy for various medical conditions, particularly in managing chronic pain and metabolic disorders like DM and DPN.⁸ The World Health Organization has endorsed acupuncture for 107 indications, including diverse neurological disorders⁹, respiratory system,¹⁰ metabolic syndrome,¹¹ and as an adjunct therapy for cancers.¹² The growing body of research highlights the necessity of understanding the molecular mechanisms underlying acupuncture’s efficacy, particularly concerning its effects on gene expression and cellular signaling pathways that relate to immune responses and inflammation, which are pivotal in the context of DM and DPN. Acupuncture has been shown to reduce inflammatory responses by decreasing inflammatory cytokines, activating the cholinergic anti-inflammatory pathway, and modulating immunity.^13–15 In addition, it was recently shown that cancer patients who had cisplatin combined with electroacupuncture therapy had increased numbers of natural killer (NK) cells as well as decreased tumor sizes.¹⁶ While numerous studies have documented the beneficial effects of acupuncture, there remains a pressing need for systematic investigations that elucidate the gene expression changes associated with this therapeutic modality, particularly in the context of DM, DPN, and adjunct therapy of cancer patients.

Bioinformatics techniques have been employed to facilitate investigations and gain insights into the pathophysiological mechanisms of diseases.^17,18 One such technique, weighted gene coexpression network analysis (WGCNA), can identify connections between genes by creating a co-expression network and identifying modules related to phenotypes, making it a powerful tool for studying key pathways and genes involved in various human disorders. Recent studies have demonstrated the involvement of multiple human body systems and multi-target therapeutic effects of acupuncture.¹⁹ Consequently, it is essential to introduce innovative research ideas to investigate the potential mechanisms of acupuncture.

In this study, we first identified the molecular mechanisms of acupuncture’s anti-DM/DPN effects using bioinformatics, network topology, and network pharmacology. Then, we looked into their function in various cancer. We performed WGCNA in conjunction with immune infiltration studies for the first time to pinpoint potential immunological mechanisms of acupuncture in DM/DPN and other malignancies (Supplement Figure 1). This research might give a solid theoretical foundation for the future development of acupuncture therapy for DM, DPN, and cancer.

Materials and Methods

Data Acquisition and Differential Expression Analysis

GSE95849 was chosen to acquire DM and DPN mRNA expression data from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/gds) databases. The collection yielded a total of 18 samples, including 6 healthy, 6 DM and 6 DPN samples. Then, the mRNA expression matrices of the GSE95849 dataset underwent log2 normalization through Sangerbox, and subsequently, volcano plots and heatmaps of the differentially expressed genes (DEG) were efficiently generated and conducted by Sangerbox (http://www.sangerbox.com/).²⁰ DEG with |fold change| ≥1.5 and FDR <0.05 were deemed significantly differentially expressed.

Weighted Gene Co-Expression Network Analysis (WGCNA)

We used WGCNA²¹ to build scale-free co-expression networks for the GSE95849 dataset in order to find groups of highly correlated genes and hub genes. To begin, Pearson correlation-based hierarchical grouping was used to group all of the genes and data. Then, in line with a scale-free co-expression network, a soft threshold power value (power = 0.7) of network building was chosen to create the co-expression network. Third, the adjacency matrices were converted into a topological overlap matrix for the identification of gene modules. Similar components were clustered and merged in the subsequent stage. Finally, the module-trait association maps were successfully created, leading to the identification of gene lists for modules that showed a significant correlation with DM and DPN, specifically those with a correlation coefficient (|r|) greater than 0.6.

Identification of Genes Targeted by Acupuncture in DM/DPN and Functional Enrichment

To identify acupuncture target, several databases were examined using “acupuncture” or “electroacupuncture” as the keyword. Among the datasets investigated were Coremine Medical (http://www.coremine.com/), GenCLiP 3 (http://cismu.net/genclip3/analysis.php) and GeneCards (https://www.genecards.org/). The target genes were then combined to identify all acupuncture-related genes. For WGCNA-related genes, the top two modules, which are most closely linked to DM/DPN, were found using the GSE95849 dataset’s module-trait association diagrams. The intersection targets of acupuncture and the effective targets of DM or DPN (DEGs and WGCNA-related genes) were loaded into Venny 2.1 (https://bioinfogp.cnb.csic.es/tools/venny/) to acquire and display the findings. The intersection genes were applied for functional enrichment analysis through Metascape software,²² which combines functional enrichment, interactome analysis, gene annotation, and membership search to harness over 40 distinct databases. GeneMANIA (http://genemania.org/search/),²³ a database encompassing protein–protein, protein–DNA, and genetic interactions, pathways, reactions, gene and protein expression information, protein domains, and phenotypic screening profiles, was applied to build a gene–gene interaction network for intersection genes to predict the gene functions.

Immune and Stromal Cell Infiltration Profiles

The immune and stromal cell composition of “bulk tissue” can be determined using the well-known xCell²⁴ using a microarray or RNA-seq data matrix that contains 64 different cell types. Based on GSE95849 microarray normalized matrix, xCell analysis was performed using the xCell website (https://xcell.ucsf.edu/). The cell enrich score of stromal and immune cells was extracted, and the difference was conducted using the Wilcoxon test to determine the differently expressed cell types in DM/DPN.

Common Immune Genes in DM and DPN

Since the immune response is important not only in DM/DPN but also in cancer, we used overlapping immune genes to see if there is any link among these diseases. Immune-related genes (IRGs) data were retrieved from the ImmPort database (https://immport.niaid.nih.gov) and GeneCards (using “immune system” as the keyword). Overlapping immune-related genes from the DEGs, WGCNA-related genes and IRGs were selected for further analysis.

Analysis of the Hub Immune Genes’ Signaling Information Networks and TF-miRNA Regulatory Networks

We then explore the interplay among TFs, miRNAs, and the designated hub genes. The previous hub immune genes’ signal transduction connection was investigated using SIGNOR 2.0 (https://signor.uniroma2.it/),²⁵ a public repository that holds manually annotated causal links between targets and other biologically relevant items. The hub genes were then subjected to TF-miRNA coregulatory network analysis using the RegNetwork repository²⁶ and NetworkAnalyst online tools²⁷, two curated regulatory interaction database.

Pan-Cancer Differential Gene Expression and Survival Analysis

We obtained the TCGA TARGET GTEx (PANCAN, N = 19131) pan-cancer data set from the UCSC (https://xenabrowser.net/) database. Furthermore, we retrieved the ESR1, CCR2, SLC1A1, CTSD, MAPK3, IL1RN, NPY and SOD2 gene expression data from each sample and screened the samples from Solid Tissue Normal, Primary Solid Tumor, Primary Tumor, Normal Tissue, Primary Blood Derived Cancer-Bone Marrow, and Primary Blood Derived Cancer-Peripheral Blood. In addition, each expression value was log2(x + 0.001) transformed. And we removed cancer types that had less than three samples in a single cancer type. Finally, 34 cancer types’ expression data were acquired. The expression difference between normal and tumor samples was calculated using Sangerbox, and the difference was analyzed using unpaired Wilcoxon Rank Sum and Signed Rank Tests.

For survival analysis, we obtained a high-quality TCGA prognosis dataset from a previously published TCGA prognosis study and supplemented it with the TARGET follow-up data from UCSC. Samples with a follow-up time of less than 30 days were eliminated. We removed cancer types that had less than ten samples in a single cancer. Finally, the expression data of 44 tumors were acquired, as well as the overall survival of the corresponding samples. To investigate the association between gene expression and prognosis in each tumor, we utilized the coxph function of the Sangerbox to create a Cox proportional hazards regression model, and the Logrank test statistical test was employed to determine prognostic significance.

Pan-Cancer Analysis of Single Nucleotide Variations

GSCALite (http://bioinfo.life.hust.edu.cn/web/GSCALite/),²⁸ an integrated genomic and immunogenomic method, was utilized to assess the single nucleotide variations of the hub genes across 33 TCGA tumor subtypes in 10,995 samples. This study included seven different types of mutations: Nonsense_Mutation, Frame_Shift_Ins, Splice_Site, Frame_Shift_Del, In_Frame_Del, In_Frame_Ins. These were referred to as detrimental mutations. The heatmap summarized the prevalence of harmful mutations in various cancer types. The waterfall plot depicted the mutation distribution of mutated genes from the inputted gene set in the sample set of malignancies. It also provided SNV type categorization, and only samples with at least one SNV were included. The variant classification method was used to calculate the number of each type of detrimental mutation in the inputted gene set in various cancer types.

Pan-Cancer Analysis of Copy Number Variation and Methylation

We also utilized GSCALite to assess the hub genes’ copy number variation and methylation. CNV data from 11495 samples were obtained from the TCGA database and analysed using GISTIC2.0, which looked for substantially changed areas of amplification or deletion across groups of patients. This study examined four kinds of CNV: homozygous deletion, heterozygous deletion, heterozygous amplification, and homozygous amplification. The CNV pie plot was used to display the Heterozygous/Homozygous CNV of each gene in each tumour. The association between gene mRNA expression and CNV was determined using Spearman correlation analysis. The P-value was adjusted by FDR, with the black outline border indicating an FDR of ≤0.05.

The TCGA database was used to get Illumina HumanMethylation 450k level 3 data. This research focused on types that had more than ten matched tumor and surrounding non-tumor samples. In general, numerous methylation sites can be found in the same gene area. As a result, each site had numerous tags that store the methylation level. We used correlation analysis to choose the location that was most adversely linked with gene expression for inclusion in this module. The association between gene mRNA expression and methylation levels was determined using Spearman correlation analysis. The P-value was adjusted by FDR, with the black outline border indicating an FDR of ≤0.05.

Developing and Testing a Nomogram Model

GSVA (Gene Set Enrichment Analysis) analysis was used to estimate the variation of the eight hub gene set activity (represents as GSVA score) over a specific cancer sample population in an unsupervised manner. In addition, the value of angiogenesis, epithelial to mesenchymal transition (EMT), cell cycle gene set was enumerated as angiogenesis GSVA-score, EMT GSVA-score, G1S GSVA-score, and G2M GSVA-score, respectively. The eight hub potential biomarkers were used to create a nomogram model using the “regplot” software. The calibration curve was also utilized to assess the model’s capacity for prediction. Last but not least, to calculate the clinical value of the model, decision curve analysis (DCA) and clinical impact curve (CIC) were used.

Results

Identification of DEGs in Diabetes Mellitus (DM) and Diabetic Peripheral Neuropathy (DPN)

The GSE95849 dataset contained the gene expression matrix, including 6 healthy, 6 DM and 6 DPN samples. After data preprocessing and differential expression analysis, we identified 3217 DEGs in DM compared with healthy samples (|Fold change| ≥ 1.5 and FDR value < 0.05), including 2982 up-regulated and 235 down-regulated genes (Figure 1A). We also identified 2191 DEGs in DPN compared with DM patients (|Fold change| ≥ 1.5 and FDR value < 0.05), including 916 up-regulated and 1275 down-regulated genes (Figure 1B). The top 100 DEGs of DM (Figure 1C) were visualized in heatmaps, including FAM69A, PREPL, SRPK2, MITD1, STT3A, OTUD7A, RCL1, ACBD7, ROS1 and B3GALT5. The top 100 DEGs of DPN (Figure 1D) were visualized in heatmaps, including LOC401397, ZNF248, NHLRC3, DCAF7, PRDX3, SLC17A6, PITX3, VENTXP7, CCDC85C and PNPLA2.

Figure 1 Identification of differentially expressed genes (DEG) in DM and DPN. Volcano plots of DEG in DM (A) and DPN (B). Heatmaps of the top 50 DEG in DM (C) and DPN (D). Red, up-regulated DEG; blue, down-regulated DEG.

WGCNA of the Whole Transcriptome Expression Matrix

Genes with similar expression patterns tended to exert similar biological functions. Therefore, we performed WGCNA according to the transcriptome expression matrix. We analyzed the soft threshold powers of the network topology and selected β = 10 as the optimal soft-thresholding parameters (Figure 2B and C). We identified 12 modules from the GSE95849 dataset (Figure 2A), and plotted module-trait diagrams to explore the relationships between gene modules and DM/DPN. The modules with the highest correlation with DM (|r| > 0.6) were MEpurple (r = 0.82, p = 3e-05) and MEbrown (r = 0.65, p = 0.004), which contained 587 hub genes. The modules with the significant correlation (|r| > 0.6) with DPN were MEred (r =−0.88, p = 1e-06), MEblue (r = 0.83, p = 2e-05) and MEpurple (r = −0.66, p = 0.003), which contained 1254 hub genes (Figure 2D).

Figure 2 Identification of modules correlated with the clinical phenotype in DM and DPN. Hierarchical clustering dendrograms of samples in the GSE958497 (A). Analysis of the scale Independence and mean connectivity for the optimal soft threshold powers in the dataset (B and C). Module-trait relationship diagrams in the GSE dataset (D). Each row corresponds to a color module and each column corresponds to a clinical trait (DM or DPN). Each cell contains the corresponding correlation and P-value.

Identification of Genes Targeted by Acupuncture in DM/DPN and Functional Enrichment

The target information of acupuncture was obtained from Coremine Medical, GenCLiP 3 and GeneCards. Once duplicate genes were deleted, a total of 1830 targets for acupuncture were obtained. Two methods were applied to achieve the key genes in DM or DPN. WGCNA was used to identify the related modules. Differential expression analysis was performed to obtain the most dysregulated genes. In total, we identified 21 concurrent genes in acupuncture genes, DM DEGs and WGCNA-related genes (Figure 3A). These genes were significantly found to participate in the formation of autolysosome, tertiary granule, promyelocytic leukemia protein (PML) body, fibrillar center, extracellular matrix; the biological processes included cellular chemical homeostasis, response to peptide, regulation of system process, macroautophagy, inflammatory response; the molecular functions included receptor tyrosine kinase binding, protein kinase C binding, G protein-coupled receptor binding, phosphotransferase activity, alcohol group as acceptor, peptide binding; signal pathways enriched in Cytokine Signaling in Immune system, CAMKK2 pathway, RANKL/RANK signaling pathway, Metabolism of Angiotensinogen to Angiotensins, Adrenergic signaling in cardiomyocytes (Figure 4). The hub genes (enriched in more than 4 terms) were MAPK3, CTSD, IL1RN, AGER, ADM, ATP1A3, AGT, NLRP6, SOD2 and HK2.

Figure 3 Pharmacological targets of acupuncture in DM and DPN. (A) twenty one intersection genes of acupuncture against DM using text mining, differential expression analysis and WGCNA analysis. (B) forty three intersection genes of acupuncture against DM using text mining, differential expression analysis and WGCNA analysis.

Figure 4 Chord plots (top 6) of the functional and pathway enrichment analyses of the overlapping genes between acupuncture and DM. Biological process (A), cellular component (B), molecular function (C) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (D) enrichment analyses.

In addition, a total of 43 concurrent genes were found in acupuncture genes, DPN DEGs and WGCNA-related genes (Figure 3B). They were mainly involved in biological processes, such as negative regulation of catalytic activity, regulation of inflammatory response, positive regulation of locomotion, cellular metal ion homeostasis, blood vessel morphogenesis; cellular components including axon terminus, cell body membrane, peroxisomal matrix, focal adhesion, euchromatin; molecular functions including G protein-coupled receptor binding, molecular function inhibitor activity, steroid binding, receptor ligand activity, RNA polymerase II-specific DNA-binding transcription factor binding; signal pathways including GPCR ligand binding, cAMP signaling pathway, Formation of Senescence-Associated Heterochromatin Foci (SAHF), FOXO-mediated transcription, Relationship between inflammation, COX-2 and EGFR (Figure 5). The hub genes (enriched in more than 4 terms) were NPY, PIK3CB, PTGER2, CCR2, AGT, ESR1, SLC1A1, KLF4 and RECK.

Figure 5 Chord plots (top 6) of the functional and pathway enrichment analyses of the overlapping genes between acupuncture and DPN. Biological process (A), cellular component (B), molecular function (C) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (D) enrichment analyses.

GeneMANIA Analysis

GeneMANIA constructed a DM network of 41 genes, including 20 related genes and 146 total links. A total of 21 DM hub genes displayed a complex PPI network with co-expression of 69.85%, predicted of 14.99%, physical interactions of 9.98%, pathways of 2.92%, and co-localization of 2.26%. The top ten proteins with the highest degrees were SQSTM1, MAPK3, SOD2, IL1RN, CTSD, MAP1LC3A, HCAR3, FTH1, GAB2 and PI3 (Figure 6A).

Figure 6 Network of the overlapping genes of acupuncture in DM (A) and DPN (B). Black protein nodes indicate target proteins, and different connecting colors represent different correlations. Genes in black circles were submitted as query terms in searches. Grey circles indicate genes associated with query genes.

In addition, we also built a DPN network of 63 genes, including 20 related genes and 288 total links. A total of 43 DPN hub genes displayed a complex PPI network with physical interactions of 57.11%, co-expression of 34.12%, shared protein domains of 3.14%, predicted of 2.61%, co-localization of 1.54%, and genetic interactions of 1.47%. The top ten proteins with the highest degrees were ESR1, MTNR1B, CAST, CAT, SLC1A1, SCP2, PROX1, RB1, CCR2 and NPY (Figure 6B).

Altered of Immune Cell Infiltration in DM and DPN

Compared to normal samples, DM showed two significant up-regulated (Macrophages M0 and Neutrophils) and one down-regulated cell types (Macrophages M2) (Figure 7). We observed 6 significant differences in the infiltration of cells between DPN and normal samples. In DPN, 2 cell types were up-regulated, including Monocytes and Neutrophils; the four down-regulated cell types were T cells follicular helper, T cells regulatory, Macrophages M2 and Eosinophils. However, T cells CD4 naïve, T cells follicular helper and Macrophages M0 were down-regulated in DPN compared to DM patients.

Figure 7 The results of immune and stromal cell infiltration analysis.

Common Immune Genes in DM and DPN

Since the immune response is important not only in DM/DPN but also in cancer, we used overlapping immune genes to see if there is any link among these diseases. 3045 IRGs were obtained from ImmPort and GeneCards database. The overlapping genes were further retrieved from the DM (MAPK3, IL1RN, SOD2, and CTSD) and DPN (ESR1, SLC1A1, NPY, and CCR2), respectively (Figure 8).

Figure 8 Venn plot displaying the identification of acupuncture-immune-related hub genes in DM (A) and DPN (B).

TF–miRNA Coregulatory Network

NetworkAnalyst online platform is used to generate TF–miRNA co-regulatory network. The interaction between TFs and miRNAs and the selected hub genes was delivered via the TF–miRNA co-regulatory network analysis. This interaction may be the underlying mechanism for the regulation of hub genes expression. The TF–miRNA co-regulatory network of DM was consisted of 83 nodes and 86 edges. A total of 13 miRNAs and 64 TF genes interacted with the validated hub genes, including MAX, MYC, NR3C1, hsa-miR-125b, hsa-miR-132, and hsa-miR-147. Figure 9A shows the TF–miRNA co-regulatory network. MAPK3 owned the largest number of neighbors, followed by IL1RN, SOD2, and CTSD.

Figure 9 TF–miRNA coregulatory and signaling information network of the hub genes. DM TF–miRNA coregulatory network (A), DPN TF–miRNA coregulatory network (B), DM signaling information network (C), DPN signaling information network (D). The round nodes are validated hub genes, the diamond nodes are TF genes, and the triangle nodes are miRNA genes.

The TF–miRNA co-regulatory network of DPN was consisted of 787 nodes and 790 edges. A total of 99 miRNAs and 128 TF genes interacted with the validated hub genes, including TBP, SP1, POU2F1, hsa-miR-26a, hsa-miR-26b, and hsa-miR-9. Figure 9 shows the TF–miRNA co-regulatory network. ESR1 owned the largest number of neighbors, followed by SLC1A1, NPY, and CCR2 (Figure 9B).

Construct the Signaling Information Network Based on SIGNOR

As shown in the Figure 9C, we constructed the DM signaling information network of the four hub genes based on SIGNOR2.0. MAPK3 showed the highest number of interaction, followed by CTSD, SOD2 and IL1RN. Moreover, STAT1 interacted with more complex than other hub genes, including GSK3B/Axin/APC, mTORC1 complex. In addition, MAPK3 showed the interactions with other two genes: MAPK3 can activate CTSD through TP53; increase IL1RN via STAT3.

In DPN signaling information network, ESR1 owned the highest number of interaction, followed by NPY, CCR2 and SLC1A1. Moreover, ESR1 interacted with more complex than other hub genes, including CyclinA2/CDK2, TFIIH, PI3K, AP1, MLL2 complex. In addition, CCR2 can increase ESR1 through ERK1/2. We also identified NPY repressed the effects of food intake (Figure 9D).

Pan-Cancer Hub Gene Expression

In order to study the changes of hub gene expression levels in different tumor tissues compared with normal tissues. We first analyzed TCGA, TARGET, and GTEx databases. The analysis showed that MAPK3 expression increased significantly in STAD, KIRC, LIHC, PAAD, PCPG, and CHOL; decreased significantly in GBM, UCEC, BRCA, CESC, LUAD, KIRP, PRAD, HNSC, LUSC, WT, SKCM, BLCA, THCA, READ, OV, TGCT, UCS, ALL, LAML, ACC, and KICH (Figure 10A). CTSD expression increased significantly in GBM, LGG, UCEC, BRCA, CESC, LUAD, ESCA, KIRP, COAD, PRAD, STAD, HNSC, KIRC, LIHC, SKCM, THCA, OV, PAAD, TGCT, UCS, KICH, and CHOL; decreased significantly in LUSC, WT, READ, ALL, and ACC (Figure 10B). SOD2 expression increased significantly in GBM, LGG, ESCA, KIRP, STAD, KIRC, PAAD, and TGCT; decreased significantly in BRCA, LUAD, PRAD, LUSC, WT, SKCM, BLCA, THCA, OV, UCS, ALL, LAML, PCPG, ACC, and KICH (Figure 10C). IL1RN expression increased significantly in GBM, LGG, UCEC, BRCA, CESC, ESCA, COAD, STAD, LUSC, THCA, OV, PAAD, and UCS; decreased significantly in LUAD, PRAD, HNSC, LIHC, WT, SKCM, TGCT, ALL, LAML, ACC, KICH, and CHOL (Figure 10D).

Figure 10 Expression status of the eight hub genes in the pancancer analysis. (A) MAPK3, (B) IL1RN, (C) SOD2, (D) CTSD, (E) ESR1, (F) SLC1A1, (G) NPY and (H) CCR2.

ESR1 expression increased significantly in BRCA, OV, PAAD, ALL, and LAML; decreased significantly in LGG, UCEC, CESC, LUAD, ESCA, KIRP, COAD, PRAD, KIRC, LUSC, LIHC, WT, SKCM, BLCA, THCA, READ, TGCT, UCS, PCPG, ACC, KICH, and CHOL (Figure 10E). NPY expression increased significantly in PRAD, WT, UCS, ALL, and PCPG; decreased significantly in GBM, BRCA, CESC, LUAD, ESCA, KIRP, COAD, STAD, KIRC, LUSC, LIHC, SKCM, BLCA, THCA, READ, TGCT, LAML, ACC, and KICH (Figure 10F). CCR2 expression increased significantly in GBM, LGG, BRCA, KIRP, KIRC, OV, PAAD, TGCT, and LAML; decreased significantly in UCEC, LUAD, ESCA, COAD, PRAD, LUSC, WT, BLCA, READ, UCS, ALL, ACC, and KICH (Figure 10G). SLC1A1 expression increased significantly in LGG, ESCA, KIRP, STAD, KIRC, OV, PAAD, LAML, and PCPG; decreased significantly in GBM, UCEC, LUAD, HNSC, LUSC, LIHC, WT, SKCM, BLCA, THCA, READ, TGCT, UCS, ALL, and KICH (Figure 10H).

Hub Gene’s Prognostic Value in Pan-Cancer

We looked into the relationship between hub gene expression and pan-cancer patient overall survival (OS) prognosis. Cox regression of 33 tumors revealed that MAPK3 expression was substantially associated with OS in eight cancers: LAML, LUSC, LIHC, and BLCA were risk factors, whereas LGG, KIRC, MESO, and PCPG were protective factors (Figure 11A). CTSD expression was substantially associated with OS in nine cancers, with risk factors in LGG, LAML, GBM, LUSC, THCA, TGCT, and protective factors in KIRP, KIRC, and DLBC (Figure 11B). SOD2 expression was substantially associated with OS in nine cancers, with risk factors in LGG, CESC, LAML, KIPAN, THYM, TGCT, PAAD, and UVM and protective factors in SARC and SKCM (Figure 11C). IL1RN expression was substantially associated with OS in nine cancers: LGG, LAML, KIPAN, KIRC, THYM, TGCT, PAAD, and UVM were risk factors, whereas LIHC was a protective factor (Figure 11D). ESR1 expression was substantially associated with OS in nine cancers: LGG, STES, and STAD were risk factors, whereas HNSC, KIRC, LIHC, and SKCM were protective factors (Figure 12A). NPY expression was substantially associated with OS in nine cancers: HNSC and OV were risk factors, whereas LUSC and PAAD were protective factors (Figure 12B). CCR2 expression was substantially associated with OS in nine cancers: KIPAN, LGG, and UVM were risk factors, whereas LUAD, SARC, HNSC, THCA, MESO, and SKCM were protective factors (Figure 12C). SLC1A1 expression was substantially associated with OS in nine cancers: UVM and UCS were risk factors, whereas LGG, KIRP and KIRC were protective factors (Figure 12D).

Figure 11 Survival analysis of the DM hub genes in the pancancer analysis. (A) MAPK3, (B) IL1RN, (C) SOD2, (D) CTSD.

Figure 12 Survival analysis of the DPN hub genes in the pancancer analysis. (A) ESR1, (B) SLC1A1, (C) NPY and (D) CCR2.

SNVs of the Hub Genes are Co-Occurred with Other Oncogenic Mutations

Our analysis of the single nucleotide variation in the gene signature revealed that ESR1, CCR2, SLC1A1, CTSD, MAPK3, IL1RN, NPY and SOD2 have mutation frequencies of 32%, 20%, 19%, 16%, 11%, 9%, 8% and 7%, respectively, across the TCGA cancer types (Figure 13A). According to the variant classification, a missense mutation is the most prevalent SNV of these genes in TCGA cancer cohorts (Figure 13B). Specifically, the top three mutations transitions were the C > T, C > A and T > C (Figure 13C). Among the genes, ESR1, CCR2 and SLC1A1 were the most frequently mutated, while IL1RN, NPY and SOD2 were less frequently mutated across the cancer types. However, according to the cancer types, SNV occurred most frequently in the order of UCEC, SKCM, COAD, READ, STAD, LUAD, BLCA, ACC, LUSC and HNSC (Figure 13D).

Figure 13 The eight hub genes’ SNVs frequently coexist with other oncogenic alterations. (A) The waterfall plot displaying the distribution of mutations and the SNV types’ SNV categorization. The summary plot of the eight hub genes’ SNV variant categorization (B) and class (C) in the entire TCGA database. (D) A heat map displaying the eight hub genes’ SNV prevalence across various cancer types.

Hub Genes’ Expression Exhibited a Tumor-Context-Dependent Association with Copy Number Variation, and Gene Methylation

Somatic copy number alterations (SCNAs) are widespread in human cancers and have been suggested to drive tumorigenesis. We found that CNVs of the 8 hub immune-related genes showed a high frequency in READ, KIRP, COAD, ACC, DLBC, UCS, SARC, TGCT, LUAD, MESO, CESC, KIRC, LUSC and PAAD, while less frequently in THCA, PRAD, LAML, THYM, LGG and PCPG. And gene heterozygous amplification and heterozygous deletion were the most frequently occurring CNVs (Figure 14A). Among the genes, MAPK3, CTSD and SOD2 were the highest correlations of CNV with mRNA expression across the cancer types. However, according to the cancer types, significant CNV occurred most frequently in the order of LUAD, OV, LUSC, SKCM, BRCA, HNSC, LIHC and STAD (Figure 14B). Furthermore, the 8 hub immune-related genes were differential methylated (Figure 14C) and related to the mRNA expression in different cancers. We found that methylation of CTSD, IL1RN and ESR1 were the most frequently methylated genes, while NPY, CCR2 and SOD2 were less frequently methylated across the cancer types (Figure 14D).

Figure 14 CNV and gene methylation are associated with the expression of the eight hub genes. (A) A pie chart depicting the distribution of the eight hub genes’ CNVs across TCGA cancer types. Hete Del stands for heterozygous deletion; Hete Amp stands for heterozygous amplification; Homo Del stands for homozygous deletion; Homo Amp stands for homozygous amplification. (B) Correlation between CNVs and mRNA expression of the eight hub genes in various cancer types. (C) A bubble plot illustrating the methylation changes between tumor and control samples. (D) The relationship between methylation and mRNA expression in various malignancies. The negative and positive correlations are shown by the blue and red bubbles, respectively. The stronger the association, the deeper the hue. The size of the point indicates statistical significance; the larger the size, the greater the statistical significance.

Establishment and Evaluation of the Nomogram Model

Then, the nomogram model was built to evaluate the 1-year, 3-year, 5-year, 7-year, and 10-year OS based on the gsva score, age, stage, T, M, N, sex, and cancer type (Supplement Figure 2). The C-index was 0.783 and 0.79 in TCGA training and test datasets, respectively. The calibration curve of the model for the possibility of OS at 1-year, 3-year, 5-year, 7-year, and 10-year suggested brilliant agreement among the apparent curve, bias-corrected curve, and ideal curve (Figure 15A and B). In addition, the time-dependent ROC curves showed good prediction efficiency of the nomograms compared for 1-year, 3-year, 5-year, 7-year, and 10-year OS in TCGA training and test datasets (all AUC > 0.8, Figure 15C and D). DCA showed that the nomogram model showed the best net benefit for 1-year, 3-year, 5-year, 7-year, and 10-year OS as well (Figure 15E and F). These results indicated that the prediction efficiency of the model was good.

Figure 15 Evaluation of the nomogram model. (A) Calibration curve in the training dataset. (B) Calibration curve in the test dataset. (C) The time-dependent ROC analysis in the training dataset. (D) The time-dependent ROC analysis in the training dataset. (E) Decision curve analysis in the training dataset. (F) Decision curve analysis in the test dataset.

The AMNC Signature and Malignant Features of the Tumor

During the transformation of cells from a benign to a malignant phenotype, they acquire characteristics such as accelerated growth, enhanced EMT, and the ability to stimulate the formation of new blood vessels, all of which are defining traits of cancer. To explore the relationship between AMNC signature and these cancerous properties, we utilized a GSVA-score algorithm to evaluate the influence of AMNC signature on tumor promotion, angiogenesis, EMT, and cell cycle progression. Our analysis revealed significant positive correlations between the AMNC signature score and the scores for angiogenesis (R = 0.2, p = 1.7e-84), G1S (R = −0.33, p = 3.5e-244), G2M (R = −0.29, p = 4.1e-194), and EMT (R = 0.06, p = 1.4e-8) across the entire TCGA pan-cancer dataset (Supplement Figure 3), as well as in the majority of cancer types analyzed (Figure 16). This suggests that tumors with a high capacity to induce AMNC are typically associated with more vigorous angiogenic activity in the tumor microenvironment and more aggressive tumor behavior.

Figure 16 The AMNC score was highly correlated with many malignant features of the tumor. (A) angiogenesis, (B) G1S, (C) G2M, and (D) EMT features.

Discussion

Although acupuncture may be useful in treating cancer pain, diabetes, and DPN,^6,11,29 no research has rigorously examined how it could work. Here, we sought to identify the acupuncture’s integrative mechanism for treating DM, DPN, and cancer pain by employing a powerful method of bioinformatics and network topology. This research proved that acupuncture employed multiple targets to produce a broad regulatory effect. Our findings implied that several immunological, inflammatory, and nervous system-related gene activities and signaling pathways were simultaneously implicated, which forms the basis of acupuncture treatment for cancer, DM, and DPN.

The implications of our findings are profound for clinical practice and policy-making, particularly concerning the integration of acupuncture into standard treatment protocols for DM and DPN. Given the pressing need for effective, multifaceted therapeutic strategies in managing these conditions, our results suggest that acupuncture could serve not only as a symptomatic relief measure but as a potential modulator of underlying pathophysiological processes. This aligns with the growing body of evidence supporting the inclusion of complementary therapies in conventional medicine, as seen in recent meta-analyses demonstrating improved patient outcomes when acupuncture is utilized alongside pharmacological interventions.³⁰ As health systems increasingly recognize the importance of holistic approaches, our research provides a foundation for advocating for acupuncture’s wider acceptance and implementation in diabetes care protocols. Furthermore, the study highlights immune mechanisms in acupuncture’s effects on DM and DPN, suggesting shared inflammatory pathways with cancer that could benefit multiple diseases. Acupuncture’s modulation of immune signaling may alleviate DM/DPN symptoms and aid in treating autoimmune disorders and cancer. Understanding these relationships can lead to targeted acupuncture interventions, enhancing immune modulation and clinical outcomes, marking a shift in acupuncture’s role in modern therapy for comprehensive treatment strategies.³¹

In addition, the prognostic model aids oncologists in treatment decisions by predicting patient outcomes through hub gene expression, allowing for tailored therapies that improve treatment effectiveness and survival rates. This is crucial given the historical challenges in cancer management due to varied treatment responses. The model’s strong predictive performance, with C-index values over 0.78, highlights its clinical relevance and the need for its integration into standard oncology. Additionally, the results emphasize the role of gene function and immune dynamics in DPN, particularly the upregulation of immune genes like FCER1G and SYK, which may influence neuropathic changes, offering insights into sensory deficits and potential targeted interventions for restoring immune and neuronal health.

KEGG analysis revealed that most enriched pathways were associated with inflammatory response, immune, CAMKK2 pathway, RANKL/RANK signaling pathway, GPCR ligand binding, cAMP signaling pathway. By mediating the downstream effects of hormones, metabolites, inflammatory agents, and neuroendocrine signals, CAMKK2, a calcium-responsive serine-threonine kinase, played an important role in cellular metabolism, and was linked to an increase in cancer cell proliferation via activating various signaling modules.³² Ca2+/cAMP signaling relationship is currently being extensively addressed as a factor in both cancer and diabetes development,³³ and plays a crucial role in coordinating the activation and functions of neutrophils.³⁴ An RNA-sequencing investigation discovered that electroacupuncture helped reduce experimental chronic inflammatory pain by suppressing calcium voltage-gated channel-mediated inflammation.³⁵ Moreover, acupuncture can treat depression through activating the cAMP signaling pathway and elevating the BDNF level.³⁶ Therefore, these pathways might be important pathways underlying the immunological, inflammatory, and nervous system-related role of acupuncture.

Mitogen-activated protein kinase 3 (MAPK3) belonged to the serine/threonine protein kinase family,³⁷ which was connected to the insulin signaling system.³⁸ The MAPK3 pathway was typically activated in cholangiocarcinoma, resulting in cell proliferation, differentiation, and survival.³⁹ MAPK3 was found to be a contributor to malfunction in numerous early neuronal development pathways, with altered soma and electrophysiological characteristics in adult neuronal cells.⁴⁰ CTSD was an aspartyl protease that was prevalent in the brain.⁴¹ CTSD has previously been linked to insulin resistance in newly diagnosed type 2 diabetic patients.⁴² CTSD gene mutations caused the congenital form of neuronal ceroid lipofuscinosis-10, providing further evidence for CTSD’s critical function in the brain.⁴³ SOD2 encoded the manganese superoxide dismutase, and regulated oxidative stress in the cell by dismutating superoxide radicals (O2) to H2O2 and molecular oxygen, which increased mutagenesis and carcinogenesis.^44,45 SOD2 was down-regulated in newly diagnosed DM patients,⁴⁶ and protected neurons from damage in diabetic neuropathy.⁴⁷ IL1RN encoded a protein that belonged to the interleukin 1 cytokine family. IL1RN functioned as an inhibitor of interleukin 1 alpha (IL1A) and interleukin 1 beta (IL1B) activities, and regulated various immune and inflammatory responses related to interleukin 1.⁴⁸ ESR1 belonged to the family of ligand-activated transcription factors that exerted pivotal roles in metabolism, growth, sexual development, and other reproductive functions.⁴⁹ It was also implicated in the pathogenesis of various cancers and metabolic disorders.^50–53 NPY played a pivotal role in various physiological systems such as the nervous, immune, endocrine systems, and cancer progression,^54,55 demonstrating significant effects on the proliferation, apoptosis, differentiation, migration, mobilization, and cytokine secretion of different types of cells.⁵⁶ CCR2 was a part of the chemokine receptor family that played a crucial role in promoting the spread of cancer cells and linked to the development of metabolic diseases, including type-1 diabetes, diabetic nephropathy, multiple sclerosis, atherosclerosis.^57,58 SLC1A1 was responsible for transporting glutamate and aspartate amino acids, and its overexpression can lead to impaired cytotoxic T-cell function and was linked to disease progression and clinical outcomes in solid malignancies.⁵⁹ Reduced SLC1A1 expression in mice was linked to neuroinflammation, reduced sensorimotor gating, and poor cognitive performance.⁶⁰ Therefore, these genes may be important targets for acupuncture therapy of DM, DPN and cancer.

This study, while providing valuable insights into the genetic and immune landscape associated with acupuncture in DM and DPN, is not without its limitations. One notable constraint is the reliance on bioinformatic data derived from publicly available databases, which may introduce variability due to differences in sample collection methods, data processing, and inherent biases within the datasets. Furthermore, the intersection of identified genes and their functional roles may not fully capture the complexity of the biological systems involved, as gene expression does not always correlate with protein activity or clinical outcomes. Additionally, the absence of experimental validation for the identified hub genes and immune gene networks limits the translational potential of these findings. Lastly, the generalizability of our results may be constrained by the heterogeneity of the cancer types analyzed and the relatively small sample sizes in certain categories, highlighting the need for further investigations in larger, more diverse cohorts.

Conclusions

In conclusion, our exhaustive analysis has illuminated the intricate interplay between acupuncture-targeted genes and immune responses in DM, particularly DPN, with MAPK3, IL1RN, SOD2, CTSD, ESR1, SLC1A1, NPY, and CCR2 emerging as key biomarkers. These genes, found among 3,217 DEGs in DM/control and 2,191 in DPN/DM comparisons, are not only potential therapeutic targets for diabetes-related complications but also signify their importance in cancer biology. Our study revealed their links to inflammatory responses, immune processes, and signaling pathways, supported by distinct immune cell infiltration patterns and regulatory interactions with microRNAs and transcription factors. Dysregulation of these hub genes correlated with cancer prognosis, emphasizing their role in cancer development and progression. We developed a pan-cancer nomogram model integrating clinical variables to predict patient survival, enhancing our understanding of diabetes-related complications and cancer biology. This work underscores the crucial roles of these genes in mediating the immune response to acupuncture in DM, DPN, and cancer therapy, paving the way for personalized therapeutic strategies leveraging traditional and modern medicine to reduce disease incidence and improve outcomes.