Identification and analysis of mitochondria-related differentially expressed genes in sepsis
The overall experimental procedure is depicted in Fig. 1. We downloaded three sepsis-related datasets, GSE26378, GSE57065, and GSE95233, from the GEO database and analyzed DEGs from sepsis samples and normal control samples using the "limma" package after quality control, normalization, and background correction. We used an adjusted P value < 0.05 and |log2-fold change|≥1 as thresholds to obtain differentially expressed genes. The difference analysis revealed 879 DEGs in the GSE26378 dataset, of which 653 genes were up-regulated and 226 genes were downregulated in the DCM samples compared with the normal samples; 1190 DEGs in the GSE57065 dataset, including 560 upregulated genes and 630 down-regulated genes; and 1621 DEGs in the GSE95233 dataset, including 860 up-regulated genes and 761 downregulated genes.
We retrieved mitochondria-related genes from the GeneCards database and selected genes overlapping with DEGs as MitoDEGs from the three datasets. There were a total of 26 MitoDEGs in the GSE26378 dataset (22 upregulated and 4 downregulated), 30 MitoDEGs in the GSE57065 dataset (14 upregulated and 16 downregulated), and 26 MitoDEGs (18 upregulated and 10 downregulated) in the GSE95233 dataset (Fig. 2A, B).
We visualized the MitoDEGs as volcano plots and heatmaps (Fig. 3A-F). Additionally, we visualized the proportion of each MitoDEG in the sample as a circular histogram (Fig. 3G-I). Merging the MitoDEGs for each dataset yielded 67 overlapping MitoDeGs, of which 35 genes were upregulated and 32 genes were downregulated in the sepsis samples compared to the normal samples.
Functional and pathway enrichment analysis of MitoDEGs
We performed further functional enrichment of MitoDEGs using GO and KEGG pathway analyses. GO analysis revealed that in biological processes (BP), regulation of inflammatory response, positive regulation of cytokine production, and in cellular component (CC), secretory granule lumen, cytoplasmic vesicle lumen were significantly enriched in several datasets. In molecular function (MF), protein tyrosine kinase activity, protein serine/threonine/tyrosine kinase activity, and protein phosphatase binding were significantly enriched in several datasets (Fig. 4A-I). The most enriched KEGG pathways of MitoDEGs were mainly involved in PD-L1 expression and the PD-1 checkpoint pathway in cancer, Th1 and Th2 cell differentiation, Th17 cell differentiation, and hepatitis B pathways (Fig. 5A-F).
PPI network analysis and hub MitoDEG identification
We merged the MitoDEGs of each dataset to obtain 49 overlapping MitoDEGs, in which 33 genes were upregulated and 16 genes were downregulated in the sepsis samples compared with the normal samples. We imported the MitoDEGs into the String database and constructed a network graph with 49 nodes, 320 edges, and an average node degree of 13.1 based on the default implementation criteria of the String database. We downloaded the network data from the String database and imported it into Cytoscape software for visualization (Fig. 6A). We used the plug-in MCODE of Cytoscape to identify the significant module and obtain a network graph containing 18 MitoDEGs (Fig. 6B). Using the MCC algorithm in the plug-in CytoHubba, we identified 10 candidate hub genes from the PPI network, including TP53, GAPDH, MYC, PTEN, SMAD3, MMP9, PPARG, TLR4, STAT1, and MAPK14 (Fig. 6C).
Hub MitoDEGs-TF network and ceRNA regulatory network
We used the NetworkAnalyst website to construct a gene regulatory network including a network of 44 transcription factors interacting with hub MitoDEGs (Fig. 7A). In R, we used eight databases, including ENCORI, miRDB, miRWalk, RNA22, RNAInter, TargetMiner, TargetScan, and miRTarBase to predict the miRNAs of Hub MitoDEGs. We found TP53, SMAD3, PTEN, and MMP9 in 8 databases, and a total of 10 miRNAs were predicted. Subsequently, we obtained 13 target lncRNAs by the starBase database. We constructed a ceRNA regulatory network consisting of 3 mRNAs, 10 miRNAs, and 13 lncRNAs and visualized it in Cytoscape software (Fig. 7B).
Immune cell infiltration in sepsis
We combined the sepsis datasets GSE26378 and GSE57065 after batch effect removal. After batch correction, we obtained and normalized the integrated sepsis dataset, including 164 sepsis samples and 46 normal samples. The differences between the two datasets were significantly reduced after removal of batch effects (Fig. 8A, B). To better understand the differences in immune function in sepsis, we applied the ssGSEA algorithm to analyze the immune infiltration of GSE26378 and GSE57065 samples. The results revealed significant differences between sepsis samples and normal samples in terms of infiltration of 26 immune cell types (p < 0.05). Among them, activated dendritic cells and neutrophils were more abundantly infiltrated in sepsis samples, and effector memory CD4 T cells and activated CD8 T cells were more abundantly infiltrated in normal samples (Fig. 8C-E). Analysis of multiple correlations between infiltrating immune cells in sepsis revealed the strongest synergistic effect of central memory CD4 T cells and effector memory CD4 T cells, followed by central memory CD4 T cells and activated CD8 T cell (Fig. 9A).
Hub MitoDEG correlation with immune cells and GSEA enrichment analysis
We conducted correlation analysis between hub MitoDEGs and ssGSEA immune infiltration, which showed that most genes were significantly correlated with immune cells (Fig. 9B). Among them, TP53 was positively correlated with effector memory CD4 T cells, MDSC4 and negatively correlated with type 17 T helper cells (Fig. 10A); SMAD3 was positively correlated with effector memory CD4 T cells and negatively correlated with mast cell and activated dendritic cells (Fig. 10B); PTEN was positively correlated with neutrophils and negatively correlated with activated CD4 T cells and effector memory CD4 T cells (Fig. 10C); MMP9 was positively correlated with mast cells, activated dendritic cells, activated CD8 T cells, and effector memory CD4 T cells (Fig. 10D). GASE enrichment analysis of the four hub MitoDEGs showed that the respirasome and respiratory chain complex were significantly enriched in TP53 and MMP9, and azurophil granules, specific granule, tertiary granules, and primary lysosomes were significantly enriched in SMAD3 and PTEN (Fig. 10E-H).