The GSE79962 gene expression profile was downloaded from the GEO database (http://www.ncbi.nih.gov/geo) using the GEOquery package of R software (version 4.2.1). The chip platform for GSE79962 was GPL6244 (Affymetrix Human Gene 1.0 ST Array), which consists of 20 SIC human heart tissue samples and 11 healthy human heart tissue samples. The ferroptosis-related genes were obtained from the GeneCards database (https://www.genecards.org/). All data are publicly available.

Raw data were downloaded as MINiML files from the Gene Expression Omnibus (GEO) database. Probes were converted to gene symbols according to the platform annotation information of the normalized data. Probes with more than one gene were eliminated, and the average of genes corresponding to more than one probe was calculated. The limma package in R software (version 4.2.1) was used to study differentially expressed genes. The adjusted p-value was determined to correct the false positive results in the GEO datasets. “Adjusted p < 0.05 and log(fold change) > 1 or log(fold change) < −1” were defined as the threshold for the differential expression of genes. The data for the listed DEGs were processed, and heatmaps and volcano plots were drawn using ComplexHeatmap and ggplot2 R packages.

To further confirm the underlying function of potential targets, the data were analyzed by functional enrichment. Gene Ontology (GO) is a widely used tool for annotating genes with functions, especially molecular function (MF), biological pathways (BP), and cellular components (CC). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis is a practical resource for studying gene functions and associated high-level genome functional information. To better understand mRNA involvement in pathogenesis, ClusterProfiler package (version: 3.18.0) in R was employed to analyze the GO function of potential targets and enrichment in the KEGG pathway. The R software package, pheatmap, was used to draw heatmaps. Then, gene set enrichment analysis (GSEA) was applied for further analysis of pathway regulation. We used KEGG rest API (https://www.kegg.jp/kegg/rest/keggapi.html) to obtain the latest gene annotation. Enrichment analysis was performed using the R package, clusterProfiler (version 3.14.3). In this analysis, the minimum gene set was 5 and the maximum gene set was 5,000. A p-value <0.05, false discovery rate (FDR) <0.25 and |NES| >1 indicated a significantly enriched term.

Human SIC heart and healthy heart tissue samples from the GSE79962 dataset were analyzed by the WGCNA R 1.70-3 package. First, Pearson’s correlation matrices and the average linkage method were both performed for all pair-wise genes. Then, a weighted adjacency matrix was constructed using a power function, β was a soft-thresholding parameter that could emphasize strong correlations between genes and penalize weak correlations. Meanwhile, a soft threshold was reasonably selected as the degree of scale independence reached 0.8. After choosing the power of β = 5, the adjacency was transformed into a topological overlap matrix (TOM), which could measure the network connectivity of a gene defined as the sum of its adjacency with all other genes for network Gene ration, and the corresponding dissimilarity (1-TOM) was calculated. To classify genes with similar expression profiles into gene modules, average linkage hierarchical clustering was conducted according to the TOM-based dissimilarity measure with a minimum size (gene group) of 30 for the genes dendrogram. The modules that correlated the most with the clinical traits were identified as SIC-related modules. All functions of hub genes with gene significance (GS) >0.2 and module membership (MM) >0.8 were analyzed by GO enrichment.

First, we obtained ferroptosis-related genes were from the FerrDb V2 database (http://www.zhounan.org/ferrdb/current/). Then, cuproptosis-related genes were obtained from the literature. The obtained cuproptosis-related and ferroptosis-related genes were inputed into the STRING website (https://string-db.org/), and the minimum required interaction score was set to 0.9 to obtain iron genes related to copper genes. Thus, we constructed a novel signature (CRF) by combing cuproptosis-related genes with ferroptosis-related genes.

A Venn diagram drawing tool (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to generate Venn diagrams of DEGs, WGCNA-based hub genes and CRF genes. Intersection genes were included in subsequent analyses.

The diagnostic values of intersection genes for SIC were detected by ROC curve and area under the ROC curve (AUC) analysis using the pROC R 1.17.0.1 package.

microRNAs (miRNAs) and TFs related to intersection genes were screened for based on the miRNet2/0 online database (https://www.mirnet.ca/). TFs and miRNAs related to intersection genes were identified and added to the network using Cytoscape software (version 3.8.2).

DGIdb (https://www.dgidb.org/) was used as a drug–gene interaction database to screen for drug–gene interactions and information from papers, databases, and web resources. Therapeutic drugs for intersection genes were identified based on the DGIdb.

The gene expression dataset, GSE79962, contained data from 20 SIC samples, and 11 normal myocardial tissue samples. As shown in Figure 1A, data normalization and cross comparability were assessed. Using the limma package in R software (version 4.2.1) for differential expression analysis, with adjusted p < 0.05 and |log2 FC| >1 as filtering conditions, we found that 173 genes were differentially expressed in the myocardial tissue of patients with SIC cardiomyopathy compared with normal myocardial tissue. Sixty-seven DEGs were upregulated and 106 were downregulated. The ferroptosis-related genes (NOX4, HMOX1, POR, SAT1, etc.) were highly expressed in sepsis-induced cardiomyopathy model. Recent studies showed that NOX4 was characterized in the cardiovascular system, HMOX1 can be induced by sepsis, and POR was selected as the central gene of SIC. The above results are consistent with our research. Clustering analysis of these DEGs was performed, as shown in a volcano plot (Figure 1B). The heatmap for the dataset indicated better clustering of samples and higher confidence (Figure 1C).

We performed KEGG and GO enrichment analyses on the up and downregulated DEGs in the GSE79962 dataset. The results showed that there were significant differences among the functions of DEGs. The upregulated DEGs were mainly enriched in pathways of muscle tissue development, regulation of small molecule metabolic process, regulation of actin cytoskeleton, leukocyte transendothelial migration, and AMPK signaling pathway (Figures 2A, B). The downregulated DEGs were mainly enriched in pathways of response to interleukin 1, regulation of peptidase activity, positive regulation of cytokine production, neutrophil activation involved in immune response, MAPK signaling pathway, JAK-STAT signaling pathway, and cytokine-cytokine receptor interaction (Figures 2C, D). As we all know, AMPK has a strong regulatory effect on cellular energy balance, metabolic homeostasis, inflammatory response, oxidative stress and myocardial cell survival, and is closely related to the pathogenesis of septic cardiomyopathy. MAPK signaling pathway and JAK/STAT signaling is an important pathway for the signal transduction of several key cytokines in the pathogenesis of sepsis, which can transcribe and modulate the host immune response. p38-MAPK in the MAPK family is involved in SIC signaling and apoptosis mechanism. Application of clinically used JAK/STAT inhibitors, tofacitinib and baricitinib, fully prevented IFNγ-induced cardiomyopathy, confirming the critical roles of this signaling pathway in inflammatory cardiac disease. Our results are consistent with previous studies.

GSEA showed that compared with control samples, the identified KEGG pathways were Huntington’s disease, Alzheimer’s disease, oxidative phosphorylation, citrate cycle TCA cycle, Parkinson’s disease, cardiac muscle contraction, valine leucine and isoleucine degradation, fatty acid metabolism and peroxisome (Figure 3A). The identified hallmark gene sets were oxidative phosphorylation, fatty acid metabolism, adipogenesis, UV response up, estrogen response early, apoptosis, androgen response, hypoxia, inflammatory response, estrogen response late, IL6-JAK-STAT3 signaling, P53 pathway, IL2-STAT5 signaling, unfolded protein response, TNFα signaling via NF-κB, and TGFβ signaling (Figure 3B). The human phenotype ontologies identified were arrhythmia and mitochondrion (Figures 3C, D).

A total of 22,828 genes were derived from the 31 samples of the GSE79962 dataset. These genes were used to construct a co-expression network. The cluster analysis results of the samples are shown in Figure 4. Clustering trees for each dataset were established and no outliers were found (Figure 4A). Soft threshold was reasonably selected as the degree of scale independence reached 0.8. The scale-free fit index and mean connectivity were calculated and the power of β = 5 (scale free R2 = 0.87) was selected (Figure 4B). The minimum number of genes per module was set to 30 according to the criteria of the dynamic tree-cutting algorithm. The final 37 transcriptional modules represented by different colors were identified (Figure 4C). The adjacencies of modules in the network are shown in Figure 4D. To correlate the modules with sample information, we analyzed the data according to the heatmap of module-clinical trait correlations, thereby correlating data for the clinical traits (Figure 4E). The black, floral white, magenta, and pale violet red 3 modules, which were identified as the hub modules associated with clinical traits, were used to explore the correlation between module membership (MM) and gene significance (GS) to identify the hub genes in SIC (Figures 4F–I). Furthermore, we demonstrated 411 hub genes were respectively identified from the above four modules with MM >0.8 and GS >0.2.

Based on the FerrDb V2 database, we obtained 612 ferroptosis-related genes. We also identified 16 cuproptosis-related genes from the literature (Li et al., 2022; Tsvetkov et al., 2022). A protein-protein interaction network was created using the STRING database to further explore relationships among these genes. We identified 128 ferroptosis-related genes to be closely associated with cuproptosis-related genes. Therefore, we constructed a novel signature (CRF) by combing cuproptosis-related genes with ferroptosis-related genes. The protein-protein interaction network was constructed based on the STRING online database and visualized using Cytoscape software (Figure 5A). By taking the intersection of the DEGs, WGCNA-based hub genes and CRF genes, three overlapping genes (POR, SLC7A5, and STAT3) were identified for SIC (Figure 5B). The diagnostic values of the three genes were confirmed by ROC curve and AUC analysis. As shown in Figure 5C, the AUC values of POR, SLC7A5 and STAT3 produced diagnosis powers for SIC of 0.922727, 0.990909, and 0.963636, respectively. Subsequent GSEA showed that the hallmark gene sets of the three genes were for fatty acid metabolism (Figures 6A–C).

We further investigated the regulatory mechanism of these three genes in SIC. The target miRNAs and TFs of the three genes were identified and then the TF–miRNA-hub gene network was constructed based on miRnet. Finally, a TF–miRNA-hub gene network, which included the three genes, 19 TFs, and 21 miRNAs, was constructed with 45 edges (Figure 7).

Potential therapeutic compounds for SIC associated with the three hub genes were screened for using DGIdb (Table 1). We identified 15 drugs as potential therapeutic compounds for SIC.