Exosome-related gene identification and diagnostic model construction in hepatic ischemia-reperfusion injury

Abstract

Hepatic ischemia-reperfusion injury (HIRI) may cause severe hepatic impairment, acute hepatic insufficiency, and multiorgan system collapse. Exosomes can alleviate HIRI. Therefore, this study explored the role of exosomal-related genes (ERGs) in HIRI using bioinformatics to determine the underlying molecular mechanisms and novel diagnostic markers for HIRI. We merged the GSE12720, GSE14951, and GSE15480 datasets obtained from the Gene Expression Omnibus (GEO) database into a combined gene dataset (CGD). CGD was used to identify differentially expressed genes (DEGs) based on a comparison of the HIRI and healthy control cohorts. The impact of these DEGs on HIRI was assessed through gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA). ERGs were retrieved from the GeneCards database and prior studies, and overlapped with the identified DEGs to yield the set of exosome-related differentially expressed genes (ERDEGs). Functional annotations and enrichment pathways of these genes were determined using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. Diagnostic models for HIRI were developed using least absolute shrinkage and selection operator (LASSO) regression and support vector machine (SVM) algorithms. Key genes with diagnostic value were identified from the overlap, and single-sample gene-set enrichment analysis (ssGSEA) was conducted to evaluate the immune infiltration characteristics. A molecular regulatory interaction network was established using Cytoscape software to elucidate the intricate regulatory mechanisms of key genes in HIRI. Finally, exosome score (Es) was obtained using ssGSEA and the HIRI group was divided into the Es_High and Es_Low groups based on the median Es. Gene expression was analyzed to understand the impact of all genes in the CGD on HIRI. Finally, the relative expression levels of the five key genes in the hypoxia-reoxygenation (H/R) model were determined using quantitative real-time PCR (qRT-PCR). A total of 3810 DEGs were identified through differential expression analysis of the CGD, and 61 of these ERDEGs were screened. Based on GO and KEGG enrichment analyses, the ERDEGs were mainly enriched in wound healing, MAPK, protein kinase B signaling, and other pathways. GSEA and GSVA revealed that these genes were mainly enriched in the TP53, MAPK, TGF[Formula: see text], JAK-STAT, MAPK, and NFKB pathways. Five key genes (ANXA1, HNRNPA2B1, ICAM1, PTEN, and THBS1) with diagnostic value were screened using the LASSO regression and SVM algorithms and their molecular interaction network was established using Cytoscape software. Based on ssGSEA, substantial variations were found in the expression of 18 immune cell types among the groups (p < 0.05). Finally, the Es of each HIRI patient was calculated. ERDEGs in the Es_High and Es_Low groups were enriched in the IL18, TP53, MAPK, TGF[Formula: see text], and JAK-STAT pathways. The differential expression of these five key genes in the H/R model was verified using qRT-PCR. Herein, five key genes were identified as potential diagnostic markers. Moreover, the potential impact of these genes on pathways and the regulatory mechanisms of their interaction network in HIRI were revealed. Altogether, our findings may serve as a theoretical foundation for enhancing clinical diagnosis and elucidating underlying pathogeneses.

Keywords: Diagnostic marker; Exosome; GEO; Hepatic ischemia-reperfusion injury; Immune microenvironment.

Fig. 1

Technology roadmap. PCA principal component analysis, GSEA gene set enrichment analysis, GSVA gene set variation analysis, DEGs differentially expressed genes, ERGs exosome-related genes, ERDEGs exosome-related differentially expressed genes, GO Gene Ontology, KEGG kyoto encyclopedia of genes and genomes, LASSO least absolute shrinkage and selection operator, SVM support vector machine, ROC receiver operating characteristic, RBP RNA-binding protein, TF transcription factor, Es Exosome scores.

Fig. 2

Debatching of the dataset. (A,B) Boxplot plot of CGD before (A) and after (B) normalization. (C,D) PCA plots of CGD before (C) and after (D) batch effect removal. Red represents the dataset GSE12720, purple represents the dataset GSE14951, and blue represents the dataset GSE15480. PCA principal component analysis, GEO gene expression omnibus, CGD combined GEO datasets.

Fig. 3

Differential gene expression analysis. (A) Volcano map of differentially expressed genes between HIRI and Normal group in CGD. (B) Venn diagram of DEGs and ERGs in CGD. (C) ERDEGs in CGD. DEGs differentially expressed genes, ERGs exosome-related genes, ERDEGs exosome-related differentially expressed genes, GEO gene expression omnibus, CGD combined GEO datasets.

Fig. 4

GO and KEGG enrichment analyses of the ERDEGs. (A,B) Bubble Diagram of GO and KEGG Enrichment Analysis for ERDEGs: BP, CC, MF and KEGG. GO terms and KEGG terms are shown on the abscissa. In the bubble plot, the size of the bubble corresponds to the number of genes, while the color represents the magnitude of the p value. A deeper red color indicates a smaller p value, and a bluer color signifies a larger p value. (C,D) Network Diagram of GO and KEGG Enrichment Analysis for ERDEGs, this diagram display BP, CC, MF (C), and KEGG (D). Red nodes represent terms, blue nodes represent molecules, and the lines represent the relationship between terms and molecules. (E,F) Bar graph of GO and KEGG enrichment analysis for ERDEGs. The criteria for GO KEGG) enrichment analysis were set at p.value < 0.05 and FDR value (q value) < 0.05. ERDEGs exosome-related differentially expressed genes, GO Gene Ontology, KEGG kyoto encyclopedia of genes and genomes, BP biological process, CC cellular component, MF molecular function.

Fig. 5

KEGG enrichment analysis of the ERDEGs. (A,B) The pathway map of KEGG enrichment analysis for ERDEGs (Endocytosis (A), Proteoglycans in cancer (B)) was shown. ERDEGs exosome-related differentially expressed genes, KEGG kyoto encyclopedia of genes and genomes. These pathway maps utilized information from the KEGG Pathway Database^–.

Fig. 6

GSEA for the combined datasets. (A) Presentation of five biological function mountain maps by GSEA for geneset in CGD. (B–F) GSEA showed that genesets in CGD were significantly enriched in TP53 pathway (B), MAPK pathway (C), TGFbeta pathway (D), JAK-STAT pathway (E). NFKB pathway (F). (G) Heat map of GSVA results between different groups of CGD. In the heat map, blue represents down-regulation and red represents up-regulation. The threshold value of GSVA were p value < 0.05 and FDR value (q value) < 0.25. Red represents the HIRI group and blue represents the Normal group. The threshold value of GSEA were p.value < 0.05 and FDR value (q value) < 0.25. GSEA gene set enrichment analysis, GSVA gene set variation analysis, HIRI hepatic ischemia reperfusion injury, CGD combined GEO datasets.

Fig. 7

Construction of the diagnostic model for hepatic ischemia-reperfusion injury. (A) Diagram illustrating the LASSO regression diagnostic model for ERDEGs inCGD. (B) Plot of variable trajectories of the LASSO diagnostic model. (C) Identification of the number of genes with the lowest error rate obtained by SVM algorithm. (D) Identification of the number of genes with the highest accuracy using the SVM algorithm. (E) Venn Diagram showing the intersection of genes obtained by both the LASSO and SVM algorithms. LASSO least absolute shrinkage and selection operator, SVM support vector machine, ERDEGs exosomes-associated differentially expressed genes, CGD combined GEO datasets.

Fig. 8

Diagnostic and validation analysis of HIRI. (A) Nomogram of model genes in CGD for the diagnostic model of HIRI. (B,C) Calibration curve plot (B) and DCA plot (C) of model genes in CGD for the Diagnostic Model of HIRI. (D) ROC analysis of linear predictors of Logistic regression models. (E) Heat map of GSVA results between High and Low groups of linear predictors of Logistic regression model. Blue represents down-regulation and red represents up-regulation in the heat map. The screening criteria of GSVA were p value < 0.05 and FDR value (q value) < 0.25. Red represents the High group of linear predictors of the Logistic regression model, and blue represents the Low group of linear predictors of the Logistic regression model. The ordinate of the Calibration Curve plot represents the net benefit, while the abscissa corresponds to the Threshold Probability. A value of AUC greater than 0.9 indicates a high level of accuracy. AUC area under the curve, DCA decision curve analysis, ROC receiver operating characteristic, GEO gene expression omnibus, GSVA gene set variation analysis, HIRI hepatic ischemia-reperfusion injury, CGD combined GEO datasets.

Fig. 9

Regulatory network of key genes. (A) mRNA-RBP regulatory network of key genes. (B) mRNA-TF regulatory network of key genes. (C) mRNA-drug regulatory network of key genes. (D) Key Genes predict the interaction network of genes with similar functions. Circular nodes represent genes, and the size is determined by the attributes and characteristics of the genes. Lines represent relationships, interactions or functional connections between genes, and line thicknesses represent strong associations or important interactions. Yellow is mRNA, purple is RBP, green is TF, and blue is Drug. TF transcription factor, RBP RNA-binding protein.

Fig. 10

Differential expression analysis of key genes between the HIRI and normal groups in the CGD. (A) The group comparison diagram illustrates the expression patterns of key genes in the Normal versus the HIRI group in CGD. (B) The correlation heat map of the key genes in CGD. (C) Functional similarity analysis of Key Genes. (D) Chromosomal localization map of Key Genes in human. (E–I) ROC curve analysis of Key Genes ICAM1 (E), ANXA1 (F), PTEN (G), THBS1 (H), and HNRNPA2B1 (I) in CGD. The symbol *** is equivalent to P < 0.001 and extremely statistically significant. The symbol ** represents P < 0.01, indicating highly statistical significance. The symbol * represents P < 0.05, indicating statistical significance. The closer the AUC is to 1 on the ROC curve, the better the diagnostic performance. An AUC above 0.9 indicates high accuracy. An AUC between 0.7 and 0.9 suggests moderate accuracy. ROC receiver operating characteristic, AUC area under curve. Red represents the HIRI group in the dataset and blue represents the Normal group in the dataset, GEO gene expression omnibus, HIRI hepatic ischemia reperfusion injury, CGD combined GEO datasets.

Fig. 11

Immune Infiltration analysis (ssGSEA and MCPCounter). (A) Group comparison graph for 28 types of immune cells in different groups in CGD by ssGSEA. (B) Heat map of correlation analysis between Key Genes and the infiltrating abundance of immune cells by ssGSEA. (C) Heat map of correlation analysis between Key Genes and immune cell infiltration abundance by MCPCounter. In the correlation heat map, the red circle represents the positive correlation between the genes and the infiltration abundance of immune cells. The larger the circle is, the stronger the correlation is. Blue circles represent the negative correlation between genes and the infiltrating abundance of immune cells, and the larger the circle, the stronger the correlation. The symbol ns is equivalent to P < 0.05,indicating no statistical significance. The symbol * is equivalent to P < 0.05, indicating statistical significance. The symbol ** is equivalent to P < 0.01, indicating a highly statistical significance. The symbol *** is equivalent to P < 0.001, indicating a extremely statistical significance. ssGSEA single-sample gene-set enrichment analysis, MCPCounter microenvironment cell populations-counter, HIRI hepatic ischemia reperfusion injury, CGD combined GEO datasets.

Fig. 12

GSEA for combined datasets. Five biological function mountain maps of GSEA in CGD were presented (B–F). GSEA showed that ERDEGs were significantly enriched in IL18 pathway (B), TP53 pathway (C), MAPK pathway (D), TGFbeta pathway (E). JAK-STAT pathway (F). The screening criteria of GSEA were p.value < 0.05 and FDR value (q value) < 0.25. GSEA Gene set enrichment analysis, ERDEGs exosomes-associated differentially expressed genes, GEO gene expression omnibus, CGD combined GEO datasets.

Fig. 13

Protein structure of the model genes. (A–E) Protein domains of Key Genes ANXA1 (A), HNRNPA2B1 (B), ICAM1 (C), PTEN (D) and THBS1 (E) are shown. When pLDDT < 50, the predicted structure has low confidence. When 50 < pLDDT < 70, the predicted structure has moderate confidence. When 70 < pLDDT < 90, the predicted structure has high confidence. When pLDDT > 90, the predicted structure has extremely high confidence. AlphaFoldDB AlphaFold protein structure database, pLDDT predicted local distance difference test.

Fig. 14

Verification of different expression levels of the five key genes in the H/R model. (A–E) Relative mRNA levels of ANXA1 (A), HNRNPA2B1 (B), ICAM1 (C), PTEN (D) and THBS1 (E) between H/R and control group were shown in a vivo experiment. All data were calculated as mean ± SD. The symbol * is equivalent to P < 0.05, indicating statistical significance. The symbol ** is equivalent to P < 0.01, indicating a highly statistical significance. The symbol *** is equivalent to P < 0.001, indicating a extremely statistical significance.