The flowchart employed for this research is depicted in Figure 1 . The 2 NAFLD datasets (GSE48452 and GSE89632) and 2 IS datasets (GSE16561 and GSE58294) were downloaded from the GEO database. The GSE48452 dataset (platform: GPL11532) includes 41 controls and 32 NAFLD patients. The GSE89632 dataset (platform: GPL14951) includes 24 controls and 39 NAFLD patients. The GSE16561 dataset (platform: GPL6883) includes 24 controls and 39 IS patients. The GSE58294 dataset (platform: GPL570) includes 23 controls and 69 IS patients.

The 4 original datasets were first subjected to background correction and then normalized by using the “affy” package in the R software. Subsequently, probe-level data were converted to their respective gene expression values. When a gene corresponds to multiple probes, the gene expression value is replaced with the average expression value. Following the integration of the datasets, the batch effect was eliminated by using the combat function of the “SVA” package in R, as there was a notable batch difference across datasets from different platforms. Finally, in order to identify DEGs between the disease and control groups, the “Limma” package in R software was used with the false discovery rate (FDR) of <0.05 and |log₂ Fold change (FC)| >0.6 as the screening criterion.

The weighted co-expression network for the expression matrix of the NAFLD dataset was created using the “WGCNA” package in the R software. For further analysis, the genes with a median absolute deviation (MAD) >25% were selected. To ensure data integrity, the “goodSamplesGenes” function was applied. To determine and validate the optimal soft threshold (β) and construct a scale-free network, the “pickSoftThreshold” was used. Next, the adjacency matrix was turned into a topological overlap matrix (TOM), and the gene modules were constructed using hierarchical clustering and the dynamic pruning-tree algorithm. The gene significance (GS) and module significance (MS) were calculated by using phenotype and module data in order to examine the relationship between modules or genes and the clinical features.

The “clusterProfiler” package in the R software, the functional enrichment analysis of the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) was performed, and P < 0.05 was considered to indicate statistical significance. GO and KEGG analyses were performed in accordance with the intersection genes of the most significant module genes and DEGs of NAFLD, and the intersection genes of the most significant module genes of NAFLD and DEGs of IS, respectively. In addition, we also utilize the String database (http://string-db.org/) for Protein–Protein Interaction (PPI) network construction. The intersection genes of the most significant module genes of NAFLD and the DEGs of IS were imported into the String database, and then the species was selected as “Homo sapiens”, and finally it was considered significant when it was greater than the minimum interaction score of 0.4.

To identify the important biomarkers, candidate genes for the diagnosis of NAFLD with IS were further screened by using 3 machine-learning algorithms. LASSO regression, which can be used for variable selection and regularization to increase the prediction accuracy, was performed using the “glmnet” package in the R software. RFE of the RF algorithm is a supervised machine-learning algorithm that can be applied to rank the intersection genes of NAFLD and IS as well as to determine the genes with relative importance >2 as feature genes. SVM-RFE is an SVM-based machine-learning algorithm that can determine the most suitable subset of genes while avoiding overfitting. Finally, the genes selected by the three machine learning algorithms were intersected, and the intersected genes were used for subsequent analysis.

The nomograms were constructed using the ‘rms’ package in R software based on the intersection genes screened by using three machine-learning algorithms. Then, the calibration curve was applied to evaluate the nomogram’s predictive capability. When establishing the ROC curve of the candidate genes, the area under the curve (AUC) was computed to assess the prediction accuracy of NAFLD with IS.

Using the ssGSEA algorithm of the “GSVA” package in R software, the immune infiltration of IS was quantified. Spearman’s correlations were computed to investigate the correlation between immune-infiltrating cells and the intersection of genes of machine-learning algorithms.

The “forestplot” package in R software was used to run univariate Cox regression analyses, after which the p-value, HR, and 95% CI values were computed. Subsequently, the data was visualized through a forest plot.

TIMER, xCELL, MCPCOUNTER, CIBERSORT, EPIC, and QUANTISEQ algorithms were used for further exploring the correlation between candidate gene expression and immune infiltration in all TCGA tumors. Then, we looked at the relationship between candidate genes and the expression of genes relevant to immunological checkpoints in diverse cancers. Moreover, we examined the relationship between MSI or TMB in different cancers and the expression of candidate genes.

R software (version 4.2.0; https://www.r-project.org/) was utilized for all statistical analyses and graph generation. The Student t-test was utilized to compare group differences. The predictive performance of the candidate genes utilized to build the predictive model was evaluated using the ROC curve. P<0.05 was considered to indicate statistical significance.

In the comprehensive dataset of NAFLD, 117 DEGs were identified using the “Limma” package in R. A total of 98 DEGs were identified in the combined dataset of IS. Volcano maps and heatmaps were used to visualize the DEGs of IS and NAFLD, as shown in Figures 2 , 3 , respectively. The detailed DEGs of NAFLD or IS are listed in Supplementary Table S1 and Supplementary Table S2 , respectively.

WGCNA was used to identify the most significantly related modules in the NAFLD combined dataset. The soft threshold was set at β = 6 (scale-free R² = 0.9) to fit the gene expression associated with a scale-free network ( Figure 4A ). After removing the abnormal samples, the clustering dendrogram of NAFLD and controls was obtained, as shown in Figure 4B . Next, using dynamic hybrid shearing, four gene coexpression modules were produced ( Figure 4C ). The correlation of the gene modules with NAFLD and controls is depicted in Figure 4D , and the turquoise module (containing 279 genes) demonstrated the most significant correlation with NAFLD (correlation coefficient = 0.6, p = 3e-09). Additionally, in the turquoise module, a strong association was noted between module membership and GS (correlation coefficient = 0.47, p = 9.7e-17) ( Figure 4E ). Therefore, this module was used for subsequent analyses.

First, the turquoise module genes of NAFLD with DEGs were overlapped, and 95 intersection genes were obtained using Venn diagrams ( Figure 5A ). Then, GO and KEGG functional enrichment analysis were employed to comprehend the biological functions of the intersection genes in NAFLD. GO enrichment analysis revealed that “response to extracellular stimulus,” “fat cell differentiation,” and “response to steroid hormone” were the mainly enriched biological process (BP) ( Figure 5B ). Cellular component (CC) analysis showed that the intersection genes were mainly enriched in “extrinsic component of membrane,” “endoplasmic reticulum lumen,” and “phosphatidylinositol 3-kinase complex” ( Figure 5B ). In terms of molecular function (MF), “DNA-binding transcription activator activity,” “receptor ligand activity,” and “signaling receptor activator activity” were the most important items ( Figure 5B ). According to KEGG enrichment analysis ( Figure 5C ), the intersecting genes were primarily enriched in the pathways of “IL-17 signaling pathway,” “TNF signaling pathway,” and “JAK-STAT signaling pathway.”

To further explore whether the key genes related to NAFLD were associated with the pathogenesis of IS, the turquoise module genes of NAFLD and the DEGs of IS were intersected to obtain nine intersection genes, which were displayed with Venn diagrams ( Figure 6A ). According to GO enrichment analysis, these intersection genes were primarily enriched in the BP of “epithelial cell migration,” “regulation of inflammatory response,” “neutrophil activation involved in immune response,” and “regulation of neuroinflammatory response”, CC of “clathrin-coated vesicle membrane,” “outer membrane,” and “secretory granule lumen,” and MF of “calcium-dependent protein binding,” “RAGE receptor binding,” and “IgG binding” ( Figure 6B ). According to KEGG enrichment analysis, the nine genes were primarily enriched in “transcriptional misregulation in cancer,” “IL-17 signaling pathway,” and “NF-kappa B signaling pathway” ( Figure 6C ). In summary, these nine common genes were highly related to immune regulation and inflammatory response and also associated with transcriptional dysregulation in cancer. The functional enrichment results of the common genes were similar to those of NAFLD, which implies that the common genes could be used for subsequent analysis.

In addition, after the identification of the nine intersection genes relevant to immunity and cancer, a PPI network was created to understand the interaction ( Figure 6D ).

In this study, candidate genes were screened using the three machine learning algorithms of LASSO regression, SVM-RFE, and RF. First, LASSO regression analysis was applied for the intersection genes. After removing the redundant variables, six characteristic genes were identified as potential biomarkers ( Figure 7A ). Second, for the results of SVM-RFE analysis, when the eigengene was nine, the error of the classifier was the smallest ( Figure 7B ). Next, the relative importance of the genes was ranked using the RF algorithm, and eight characteristic genes were identified ( Figures 7C, D ). Finally, the genes analyzed with the three machine algorithms were intersected, and six candidate characteristic genes (PTGS2, FCGR1A, MMP9, VNN3, S100A12, and RRS1) were recognized ( Figure 7E ).

To augment the clinical utility, a nomogram was constructed based on the six candidate genes ( Figure 8A ). The calibration curve revealed that the discrepancy between the actual and predicted values of IS risk was small, which demonstrated the high diagnostic value of the nomogram ( Figure 8B ). Moreover, ROC analysis was used to determine the AUC and 95% CI of each candidate gene. The findings were as follows: FCGR1A (AUC: 0.806, 95% CI: 0.725–0.880), MMP9 (AUC: 0.837, 95% CI: 0.775–0.897), PTGS2 (AUC: 0.768, 95% CI: 0.688–0.841), RRS1 (AUC: 0.904, 95% CI: 0.844–0.956), S100A12 (AUC: 0.865, 95% CI: 0.790–0.925), and VNN3 (AUC: 0.847, 95% CI: 0.782–0.907) ( Figure 8C ). The above results show that all candidate genes have a high diagnostic value.

The immune regulation mode in IS was further elucidated with immune cell infiltration analysis because the intersection genes of NAFLD and IS were primarily enriched in the inflammatory response and immune regulation. In the combined IS dataset, Figure 9A displays the percentage of immune cells in IS and controls. Compared with controls, lower levels of Effector.memory.CD8.T.cell, Type.2.T.helper.cell, Activated.B.cell, Activated.CD4.T.cell, Effector.memory.CD4.T.cell, Central.memory.CD8.T.cell, Activated.CD8.T.cell, CD56bright.natural.killer.cell, Memory.B.cell, and Central.memory.CD4.T.cell were observed in patients with IS. On the contrary, higher levels of Eosinophil, MDSC, Activated.dendritic.cell, Macrophage, Plasmacytoid.dendritic.cell, Gamma.delta.T.cell, Immature.B.cell, Mast.cell, Neutrophil, Regulatory.T.cell, Immature.dendritic.cell, Natural.killer.cell, and Type.17.T.helper.cell were noted ( Figure 9B ). Additionally, a correlation heatmap was used to demonstrate the correlation between the six candidate genes and immune cells ( Figure 9C ). The detailed results of correlation analysis between immune cells and candidate genes are listed in Supplementary Table S3 .

According to the KEGG enrichment assessment, the intersection genes of NAFLD and IS were associated with cancer transcriptional dysregulation. In addition, RRS1 showed the highest diagnostic value for NAFLD with IS (AUC: 0.904, 95% CI: 0.844–0.956). Therefore, the immune-related gene RRS1 was selected for further pan-cancer analysis. In the TCGA database, RRS1 was highly expressed in 18 cancer types, including CHOL, BRCA, and BLCA, with a statistical significance ( Figure 10A ). The data for normal tissue in the GTEx database were then collected, and RRS1 was found to be weakly expressed only in LAML but strongly expressed in 27 cancers, including ACC, BLCA, and BRCA ( Figure 10B ). The expression status of RRS1 in each tumor cell line is depicted in Figure 10C .

The relationship between the expression level of RRS1 and the patient prognosis was investigated by determining progression-free survival (PFS), disease-free survival (DFS), disease-specific survival (DSS), and overall survival (OS). RRS1 expression was substantially linked to LIHC and LUAD in OS analysis and acted as a risk factor in both LIHC and LUAD ( Figure 11A ). The DFS study found a significant relationship among three cancers—OV, LIHC, and PAAD—and RRS1 expression ( Figure 11B ). Furthermore, RRS1 was a protective factor in OV but a risk factor in LIHC and PAAD. In DSS analysis, the expression of RRS1 was significantly associated with four cancers, namely, KIRP, PAAD, LUAD, and UCS, and was a risk factor in all four cancers ( Figure 11C ). PFS analysis indicated that the expression of RRS1 was correlated with the PFS of seven malignancies, namely, ACC, KIRP, LIHC, OV, PAAD, PRAD, and UVM ( Figure 11D ). RRS1 served as a protective factor only in OV and was a risk factor in the remaining six malignancies.

The TIMER algorithm was used to assess the relationship between RRS1 expression and the degree of immune cell infiltration in various malignancies. RRS1 was related to B cells in 13 cancers, CD4+ T cells in 14 cancers, CD8+ T cells in 15 cancers, myeloid dendritic cells in 17 cancers, macrophages in 16 cancers, and neutrophils in 13 cancers according to TIMER analysis ( Figure 12A ). Moreover, RRS1 was significantly negatively correlated with immune cells in LGG, COAD, LUSC, BRCA, SKCM, LUAD, and STAD but significantly positively correlated with immune cells in LIHC, PCPG, and THCA. In addition, the xCELL algorithm ( Figure 12B ), MCPCOUNTER program ( Figure 12C ), CIBERSORT algorithm ( Figure 13A ), EPIC algorithm ( Figure 13B ), and QUANTISEQ algorithm ( Figure 13C ) demonstrated that diverse immune cells were highly associated with RRS1 expression. Furthermore, the expression levels in different malignancies were highly linked to immune checkpoint-related genes ( Figure 13D ). RRS1 showed a positive correlation in KICH, PCPG, and KIRP, whereas it showed a negative correlation in SKCM, LUSC, and COAD.

MSI and TMB are expected to guide tumor immunotherapy. According to the results of MSI analysis, RRS1 was significantly correlated with LUSC, MESO, DLBC, and PCBG ( Figure 14A ). In the TMB study, RRS1 was significantly associated with STAD, PRAD, LUAD, and COAD ( Figure 14B ).