Gene expression profiles from the GSE142153 dataset (28), comprising data from peripheral blood samples obtained from 23 patients with DN and 10 healthy controls, were downloaded from the Gene Expression Omnibus database (ncbi.nlm.nih.gov/geo/), which was obtained from the GPL6480 platform (Table SI). This dataset included microarray data, which were presented in the form of raw signal intensities or gene expression levels. In addition, 904 lipid metabolism-related genes were acquired from the Molecular Signature Database (https://www.gsea-msigdb.org/gsea/msigdb).

The present study employed the R package limma (29) to analyze the differentially expressed genes (DEGs) between DN and healthy samples in the GSE142153 dataset, based on the criteria |log₂FC| >0.5 and P<0.05. Next, the ggplot2(30) and pheatmap (31) packages in R (version 4.0.2, 2020, R-project.org/) were used to generate the volcano plot and heatmap of DEGs, respectively.

To find the genes linked with DN, the WGCNA package in R (32) was performed to construct a co-expression network based on the 23 DN samples and 10 healthy samples in the GSE142153 dataset. First, to ensure the accuracy of the analysis, these samples were clustered to remove the outliers. Next, to ensure the interactions of genes accord with the scale-free distribution to the maximum extent, the soft threshold of all data was determined. Subsequently, the dissimilarity coefficients of genes were calculated and the systematic clustering tree was obtained. For each gene module, the minimum module size was set as 150 with the criteria of the dynamic tree cutting, and similar modules were merged. Ultimately, the relationships among the modules and clinical traits (DN and healthy) were calculated to identify the key modules (correlation coefficient |cor| ≥0.6, P≤0.05). Module membership (MM) correlation between key module genes and the modules, and the gene significance (GS) correlation between key module genes and clinical traits were further calculated to identify the genes in DN (|MM| >0.8 and |GS| >0.2) (33).

The present study intersected DEGs, key module genes and 904 lipid metabolism-related genes to identify lipid metabolism-related genes in DN using the Venn Diagram R package (34). The associations among these genes were determined. In addition, a PPI network was constructed using STRING (https://cn.string-db.org/) database (confidence=0.15) (35) to analyze the interaction between the proteins of lipid metabolism-related genes in DN.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to investigate the potential role of lipid metabolism-related genes in DN using the clusterProfiler R package (36). The results were visualized using GOplot (37) and the enrichplot package (38).

To identify the hub genes of lipid metabolism in DN, two classification models were established using 23 DN and 10 healthy samples: The Least Absolute Shrinkage and Selector Operation (LASSO) algorithm model in R package glmnet (39) and the Support Vector Machine-Recursive Feature Elimination (SVM-RFE) algorithm model in R package e1071(40). The efficiency of these two models was determined using a 10-fold cross-validation method. The intersections of the results obtained using these two models were employed as the lipid metabolism-related hub genes in DN for the subsequent research. In addition, the R package clusterProfiler (36) was used to perform GSEA of the hub genes. The correlations between hub genes and other genes was determined using the default gene sets in the org. Hs. eg. db package (version 3.12.0, bioconductor.org/packages/release/data/annotation/html/org.Hs.eg.db.html). The significance thresholds were |Normalized Enrichment Score (NES)| >1, P<0.05 and q<0.2.

To investigate discrepancies in immune cell infiltration between the DN and healthy samples, the present study evaluated 28 types of immune cell infiltrations in the 23 DN and 10 healthy samples in the GSE142153 dataset using the ‘ssGSEA’ method in the GSVA R package (41). A heatmap was plotted of the ssGSEA scores of immune cells for each sample. The various immune cells between DN and healthy samples were evaluated using Wilcoxon rank-sum test and the results were visualized through the Vioplot package (version 0.3.7; github.com/TomKellyGenetics/vioplot). Ultimately, the Pearson correlation coefficients between each immune cell, the immune cells and hub genes were calculated and visualized using corrplot package in R (42).

The present study collected peripheral blood mononuclear cell samples from 10 healthy subjects and 10 patients with DN from the First People's Hospital of Yunnan Province (Kunming, China). The patient samples were collected between 1 and 30 March, 2022. All subjects provided written informed consent before participating in the study. The present study was approved by the Ethics Committee of The First People's Hospital of Yunnan Province (approval no. 2022GJ227).

First, RNA was extracted using TRlzol^® (Invitrogen; Thermo Fisher Scientific, Inc.). RT was performed with SweScript First-Strand cDNA Synthesis Kit (Servicebio, Wuhan, China) according to the manufacturer's protocols. qPCR) was performed strictly according to the manufacturer's instructions for 2xUniversal Blue SYBR Green qPCR Master Mix (Servicebio, Wuhan, China). qPCR reaction mixture comprised 3 µl of cDNA, 5 µl of 2xUniversal Blue SYBR Green qPCR Master Mix (Servicebio, Wuhan, China) and 1 µl of forward and reverse primer (Tsingke, China) also at a concentration of 10 µM. qPCR was performed using the CFX96 Real-Time PCR Detection System (Bio-Rad, China) according to the following steps: pre-denaturation at 95˚C for 1 min, followed by 40 cycles of denaturation at 95˚C for 20 sec, annealing at 55˚C for 20 sec, and extension at 72˚C for 30 sec. GAPDH was used as the internal reference gene, and the relative expression levels of key genes were calculated using the 2^-ΔΔCq method (43). All systematic analyses were performed in triplicate. qPCR primer sequences are shown in Table I.

All analyses were carried out in R language (version 4.0.2). Differences between groups were analyzed by Wilcoxon rank-sum test. P<0.05. The Independent-samples t-test was employed. Assuming normal distribution and equal variances, the mean and SEM) were calculated for each group. All systematic analyses were performed in triplicate.

The present study identified 1,445 DEGs between the DN and healthy samples in the GSE142153 dataset. Among the DEGs present in DN samples, 707 genes were upregulated and 738 were downregulated. A volcano plot and a heatmap of these genes are in Fig. 1.

A sample clustering tree after excluding the outliers (GSM4221586) was plotted (Fig. 2A and B). The network accorded with the scale-free distribution to the maximum extent possible when the soft threshold was 7 (Fig. 2C). After determining the soft threshold, the gene dendrogram was constructed and 15 co-expression modules were found by setting the minimum module size at 150 (Fig. 2D). The correlations among the modules and DN and healthy samples showed that MEyellow was positively correlated with DN, whereas MEturquoise had a negative correlation with DN (Fig. 2E). Finally, 694 DN-related genes in MEyellow and MEturquoise with |MM| >0.8 and |GS| >0.2 were identified (Fig. 2F and G).

A Venn diagram was used to identify 17 genes associated with lipid metabolism in DN through the intersections of 1,445 DEGs, 694 DN-related genes in the MEyellow and MEturquoise modules, and 904 lipid metabolism-related genes (Fig. 3A). The correlation analysis among these 17 genes demonstrated a strong correlation between IDI1 and PTS (Fig. 3B). The aforementioned 17 lipid metabolism-related genes were then uploaded into the STRING database. However, the SAMD8 and TNFAIP8L2 genes were excluded from the STRING database as they did not show any interaction. Next, a PPI network of 15 genes with a confidence of 0.15 was constructed. The degree scores of these 15 genes are shown in a bar chart (Fig. 3C and D).

GO and KEGG functional enrichment analyses were conducted to investigate the potential functions of the 17 lipid metabolism-related genes implicated in DN. Notably, 5 KEGG pathways and 23 GO entries belonging to the biological process (BP) category were identified (P.adjust <0.05 and count >2). The BP results showed that these 17 lipid metabolism-related genes were closely associated with the ‘phospholipid biosynthetic process’, ‘alcohol biosynthetic process’ and ‘cholesterol biosynthetic process’ (Fig. 4A). The KEGG analysis demonstrated that the markedly enriched pathways associated with these genes included ‘glycerophospholipid metabolism’, ‘steroid biosynthesis’ and ‘fatty acid degradation’ (Fig. 4B).

To identify hub genes among the aforementioned 17 lipid metabolism-related genes, LASSO and SVM-RFE algorithm models were established based on the data obtained from 23 DN samples and 10 healthy samples. LASSO algorithm identified four genes as feature genes (Fig. 5A and B). Three genes were identified as feature genes by the SVM-RFE algorithm (Fig. 5C), which are the first three genes in Table SII. The intersections of these two algorithm models revealed two overlapping genes, SAMD8 and CYP51A1 (Fig. 5D). The top 10 most important GO and KEGG terms of the GSEA results of SAMD8 and CYP51A1 were screened. ‘mitochondrial matrix’ and ‘GTPase activity’ were identified as the common markedly enriched GO term in both SAMD8 and CYP51A1 (Fig. 5E and G). The following common markedly enriched KEGG pathways for SAMD8 and CYP51A1 were identified: ‘Alzheimer disease’, ‘IL-17 signaling pathway’, ‘pathways of neurodegeneration-multiple diseases’ and ‘TGF-beta signaling pathway’ (Fig. 5F and H).

The ssGSEA was employed to evaluate the infiltration levels of 28 immune cell types in both DN and healthy samples. A heatmap displaying the ssGSEA scores is shown in Fig. 6A. A total of nine immune cell types were significantly different between DN and healthy samples: Activated B cells, CD56bright natural killer cells, effector memory CD4 T cells, effector memory CD8 T cells, mast cells, memory B cells, plasmacytoid dendritic cells, T follicular helper cells and type 1 T helper cells. These cell types exhibited significant differences between the DN and healthy samples, as visualized using a violin chart (Fig. 6B). Furthermore, the Pearson correlation coefficients of the immune cell types showed that monocytes exhibited a markedly positive correlation with plasmacytoid and immature dendritic cells (Fig. 6C). Finally, the Pearson correlation coefficients between immune cells and hub genes (SAMD8 and CYP51A1) were calculated. SAMD8 and CYP51A1 were both correlated with the activated B cells, effector memory CD8 T cells, memory B cells and T follicular helper cells (|cor| >0.4; Fig. 6D and E).

The RT-qPCR analysis showed a higher expression of SAMD8 and CYP51A1 in DN samples than in those from normal subjects, which was consistent with the aforementioned results of the present study (Fig. 7A and B).