Classification patterns identification of immunogenic cell death-related genes in heart failure based on deep learning

ICD-related gene expression landscape

The overall flow chart of this paper is given in Figure 1. Firstly, the transcriptome data of the GSE141910 data set were analyzed for differential expression, and 8885 DEGs were obtained with adj.P.Val<0.05 as the standard. In addition, we obtained 20 ICDRGs from the previous literature. Fourteen DEICDRGs were obtained by crossing them (Figure 2A). Figure 2B, 2C show the differential expression thermogram and differential expression box diagram of 14 DEICDRGs in the sick and control groups, respectively. Figure 2D displays the chromosome position information of 14 DEICDRGs. Figure 2E shows the correlation thermogram among 14 DEICDRGs. It was obvious that most genes had a significant correlation. The details of DEGs, ICDRGs, and DEICDRGs can be found in Supplementary Table 1.

Figure 1. The technical roadmap of the article.

Figure 2. The expression landscape of DEICDRGs. (A) The intersection Wayne diagram of DEGs and ICDRGs (p=3.409506e-52). (B) The expression heat map of DEICDRGs obtained by differential expression analysis. (C) The box diagram of the differential expression of DEICDRGs in HF and its control group. (D) The chromosome map of DEICDRGs. (E) The correlation analysis heat map of DEICDRGs.

Identification of subtypes of HF based on DEICDRGs

The expression levels of 14 DEICDRGs in diseased samples were extracted in the study. We set the number of clusters to 2-5, respectively. Figure 3 displays the tsne dimensionality reduction graphs of AE, DAE, and K-means under different cluster numbers. Each point in the graph represented a sample. In order to evaluate the clustering performance of the three algorithms under different cluster numbers, this paper used three index columns: sample contour coefficient, Calinski-Harabasz Index, and Davies-Bouldin Index. We can observe the histograms of the three algorithms in four clustering situations in Figure 4A–4C. From the figure, we found that DAE had the best comprehensive performance of the three indicators when the number of clusters was 2. Therefore, this paper would interpret it as the result of subtype identification in the future. Figure 4D, 4E are box graphs indicating the differences in infiltration abundance and immune function of immune cells between the two types. We also obtained the difference in gene expression related to immune inspection sites between the two subtypes (Figure 4F). According to Figure 4F, there were significant differences in the infiltration abundance of most immune cells, immune function, and expression of genes related to immune inspection sites between the two typing samples, which confirmed the typing ability of the DAE algorithm. In addition, this paper also analyzed the GSVA results for two types (Figure 4G). We would analyze the biological significance of these pathways in detail in the discussion section.

Figure 3. Identification process of HF subtypes. (A–C) are the tsne dimensionality reduction scatter plot obtained by AE, DAE, and K-means algorithms when the number of clusters is set to two categories, respectively. (D–F) are the tsne dimension reduction scatter plot obtained by three algorithms when the number of clusters is set to three categories, respectively. (G–I) are the tsne dimension reduction scatter plot obtained by three algorithms when the number of clusters is set to four categories, respectively. (J–L) are the tsne dimension reduction scatter plot obtained by three algorithms when the number of clusters is set to five categories, respectively.

Figure 4. The performance of algorithm and the analysis of subtypes in the immune microenvironment and functional biological characteristics. (A–C) are the histogram of the sample contour coefficient, Calinski-Harabasz Index, and Davies-Bouldin Index of three algorithms under four clustering numbers, respectively. (D, E) are the difference box diagram of infiltration abundance and immune function of immune cells among different types, respectively. (F) The box diagram of the difference in immune examination sites between different subtypes. (G) The GSVA analysis between subtypes.

The expression levels of 14 DEICDRGs in diseased samples were extracted in the study. We set the number of clusters to range from 2 to 5. Figure 3 displays the t-distributed stochastic neighbor embedding (t-SNE) dimensionality reduction graphs of AE, DAE, and K-means clustering under different cluster numbers. Each point in the figure represented a sample. This paper used three index columns to evaluate the clustering performance of the three algorithms under other cluster numbers: sample silhouette coefficient, Calinski-Harabasz Index, and Davies-Bouldin Index. We can observe the histograms of the three algorithms in four clustering situations in Figure 4A–4C. From the figure, we found that DAE had the best overall performance of the three indicators when the number of clusters was 2. Therefore, this paper would interpret it as the outcome of subtype identification in the future. Figure 4D, 4E are box plots indicating the differences in infiltration abundance and immune function of immune cells between the two types. We also obtained the difference in gene expression related to immune inspection sites between the two subtypes (Figure 4F). According to Figure 4F, there were significant differences in the infiltration abundance of most immune cells, immune function, and expression of genes related to immune inspection sites between the two subtypes, which confirmed the classification ability of the DAE algorithm. This article introduces non-negative matrix factorization (NMF) and consensus clustering method (CC) to confirm the algorithm’s performance further. Three clustering performance indicators were calculated when the number of clusters was 2. The dimensionality reduction results of these two algorithms and the performance comparison results with the research algorithms included in this article are shown in Supplementary Figures 1, 2 in the Supplementary Table 2.

In addition, this paper also analyzed the GSVA results for two types (Figure 4G). Through the GSVA analysis of the two subtypes, we found that there were significant differences in several biological processes between the two subtypes (Figure 4G). Almost all of these pathways are associated with heart failure—for example, melanoma, apoptosis, and the JAK-STAT signal pathway. Early scholars reported a case where metastatic malignant melanoma could result in rapid occlusion of the right ventricle, thus leading to congestive HF [9]. Melanoma often involves the heart, resulting in cardiac issues such as HF and myocardial infarction [10, 11]. HF is also intricately linked to the complex pathophysiology of apoptosis [12, 13]. Researchers have confirmed that isoproterenol (ISO) can induce apoptosis, improve heart function, and relieve and treat ISO-induced HF models and cellular HF in rats [14]. The JAK-STAT signaling pathway has been proven to play an essential role in the pathophysiology of HF [15, 16].

Verification of subtypes of HF based on DEGs cluster

To further validate the rationality of the subtype identification mentioned above, we conducted differential expression analysis on the two subtypes, identifying 43 DEGs with p < 0.05. Utilizing these DEGs, this study employed three algorithms to generate t-SNE dimensionality reduction plots under four clustering numbers (Figure 5). The evaluation results of three metrics in Figure 6A–6C are presented. It can be observed that compared to Figure 5, the clustering effect has been significantly improved. This confirms the rationality of using DAE to classify patients into two subtypes. Similarly, we provided identification results of the immune microenvironment and functional biological characteristics of different gene clusters in Figure 6D–6G. There were significant differences between the two subtypes in terms of immune cell infiltration abundance, immune functions, immune checkpoint-related gene expression, GSVA analysis pathways, and others.

Figure 5. Re-identification of HF subtypes based on DEGs. (A–C) are the tsne dimensionality reduction scatter plot obtained by AE, DAE, and K-means algorithms when the number of clusters is set to two categories, respectively. (D–F) are the tsne dimension reduction scatter plot obtained by three algorithms when the number of clusters is set to three categories. (G–I) are the tsne dimension reduction scatter plot obtained by three algorithms when the number of clusters is set to four categories, respectively. (J–L) are tsne dimension reduction scatter plot obtained by three algorithms when the number of clusters is set to five categories, respectively.

Figure 6. The performance of algorithm and the analysis of subtypes based on DEGs identification in the immune microenvironment and biological function characteristics. (A–C) are the histogram of the sample contour coefficient, Calinski-Harabasz Index, and Davies-Bouldin Index of three algorithms under four clustering numbers, respectively. (D, E) are the difference box diagram of infiltration abundance and immune function of immune cells among different types, respectively. (F) The box diagram of the difference in immune examination sites between different subtypes. (G) The GSVA analysis between subtypes.

Construction and verification of lasso model and SVM model

This paper extracted diagnosis-related genes from 14 DEICDRGs using the LASSO and SVM-RFE algorithms to identify genes relevant to HF diagnosis. Figure 7A, 7B display the curves depicting the relationship and cross-validation results of the L1 norm and coefficients obtained by applying the LASSO algorithm. Figure 7C is the result of feature selection using SVM regression. When the number of DEICDRGs was 12, the AUC of 10-fold cross-validation was 0.978. Figure 7D is the Venn diagram of the gene intersection screened by LASSO and SVM. We identified 11 genes in the intersection. In the Supplementary Table 2, we provide detailed information on diagnosis-related genes selected by the LASSO algorithm, the SVM-RFE algorithm, and their intersection genes. Figure 7E, 7F are ROC curves of the diagnostic model constructed using a set of 12 genes in both the training and test sets. Among them, the AUC in the training set reached 0.995 (CI: 0.986-0.999). The AUC in the test set reached 0.95 (CI: 0.841-1), indicating a high level of diagnostic accuracy. In addition, we also evaluated the AUC for each diagnosis-related gene.

Figure 7. The diagnosis model and ROC analysis of the model based on LASSO and SVM. (A) The relationship curve between l1 norm and coefficient in the Lasso regression of DEICDRGs. (B) The cross-validation result of Lasso regression. (C) The result of using SVM regression to filter features. (D) The Wayne diagram of the intersection of characteristic genes screened by Lasso and SVM. (E, F) are ROC analysis of the diagnosis model in the training set and test set, respectively.

Most of the diagnosis-related genes have been shown to play key roles in the development of HF. ATG5 is involved in the formation of autophagic vesicles, which may play an important role in the process of apoptosis. Autophagy is associated with HF, and autophagy activity has been detected in both patients with HF and animal models. The balance between myocardial apoptosis and autophagy in chronic HF can also be treated with drugs [17, 18]. The role of apoptosis in HF has also been determined [12–14]. CASP1 encodes a protein that is a member of the cysteine-aspartic acid protease (caspase) family. Caspases are involved in the signaling pathways of apoptosis, necrosis, and inflammation. IL1R1 is related to immune and inflammatory reactions induced by many cytokines. For a long time, we have been concerned that HF is related to systemic inflammation. Essentially, the progress of HF is attributed to the continuous signal transduction of pro-inflammatory cytokines, and the early stage of HF also shows an inflammatory state 345 [19–21]. The protein encoded by IL-10 gene is a cytokine that plays a pleiotropic role in immune regulation and inflammation and participates in the regulation of JAK-STAT signaling pathway. TNF encodes a multifunctional pro-inflammatory cytokine belonging to then TNF superfamily that is involved in regulating apoptosis. The relationship between HF and TNF was recognized as early as 1990 [22]. The mortality of patients with HF increases with the increase in TNF-α level [23]. TNF family members may represent a new target for HF treatment [24].

Figure 8A–8K show the ROC curves of the following genes (ATG5 (AUC:0.723,CI:0.67-0.775), CASP1(AUC:0.883,CI:0.848-0.917), CD8A(AUC:0.703,CI:0.644-0.755), ENTPD1(AUC:0.663,CI:0.610-0.713), IL1R1(AUC:0.748,CI:0.697-0.798), IL10(AUC:0.922,CI:0.893-0.949), IL17RA(AUC:0.842,CI:0.800-0.880), MYD88(AUC:0.690,CI:0.635-0.742), NT5E(AUC:0.913,CI:0.829-0.942), PRF1 AUC:0.869,CI:0.829-0.903, and TNF (AUC:0.677, CI:0.619-0.729)) in the training set. Figure 9A–9K display the ROC curves of the following genes (ATG5(AUC:0.723,CI:0.670-0.775), ASP1(AUC:0.883,CI:0.848-0.917), D8A(AUC:0.703,CI:-.644-0.755), NTPD1(AUC:0.663,CI:0.610-0.713), L1R1(AUC:0.748,CI:0.697-0.798), L10(AUC:0.922,CI:0.893-0.949), L17RA(AUC:0.842,CI:0.800-0.880), YD88(AUC:0.690,CI:0.635-0.742), T5E(AUC:0.913,CI:0.879-0.942), RF1(AUC:0.869,CI:0.829-0.903), and TNF(AUC:0.677,CI:0.619-0.729)) in the test set. All the genes had diagnostic significance for HF.

Figure 8. Diagnostic performance verification of diagnostic genes and diagnostic models in the training set. (A–K) are ROC curves of ATG5, CASP1, CD8A, ENTPD1, IL1R1, IL10, IL17RA, MYD88, NT5e, PRF1 and TNF in the training set, respectively.

Figure 9. Diagnostic performance verification of diagnostic genes and diagnostic models in the test set. (A–K) are ROC curves of ATG5, CASP1, CD8A, ENTPD1, IL1R1, IL10, IL17RA, MYD88, NT5e, PRF1, and TNF in the test set, respectively.

Construction of nomogram model

We developed a nomogram model by using diagnosis-associated genes (Figure 10A). The calibration curve in Figure 10B illustrated that the nomogram model had excellent diagnostic ability. The DCA of Figure 10C proved that the nomogram model has more excellent clinical utility than a single diagnosis-related gene. The clinical influence curve of Figure 10D demonstrated that the nomogram model had outstanding diagnostic ability.

Figure 10. Construction of column diagram model. (A) is a nomogram model constructed on selected diagnostically-relevant genes (ATG5, CASP1, CD8A, ENTPD1, IL1R1, IL10, IL17RA, MYD88, NT5E, PRF1, and TNF). (B) is the calibration curve used to evaluate the diagnostic capability of the nomogram model. (C) is based on DCA, and the Nomogram model had higher clinical utility than a single diagnostically-relevant gene. (D) is the clinical impact curve.

Correlation analysis of immune infiltration and construction of the regulatory network

We developed a nomogram model utilizing diagnosis-related genes (Figure 10A). The calibration curve in Figure 10B demonstrated that the nomogram model had excellent diagnostic ability. The DCA in Figure 10C indicated that the Nomogram model exhibited higher clinical utility than a single diagnosis-related gene. The clinical influence curve in Figure 10D showed that the nomogram model had remarkable diagnostic ability. We present a scatter plot of the correlation between immune cells/functions and diagnostic genes in Figure 11A–11V.

For the miRNA-mRNA interaction network (Figure 12), we confirmed that some miRNAs are related to the pathogenesis of HF by consulting the literature. MiR-423-5p was initially identified as a circulating biomarker of heart disease. Tijsen et al. proved that the circulating level of miR-423-5p in patients with clinical HF increased [25]. Deng et al. also studied and determined that miR-423-5p is a potential target for the diagnosis and treatment of HF [26]. The protein coding gene regulated by miR-107 and the gene regulated by miR-139-5p were identified as genes that play a role in HF [27]. For the miRNA-TF interaction network, we know that CREB3 can promote the expression of inflammatory genes. RELA, also known as NF-kappa-B, is a pleiotropic transcription factor, exists in almost all cell types and is the endpoint of a series of signal transduction events, which are triggered by a large number of stimuli related to many biological processes, such as inflammation, immunity, differentiation, cell growth, tumorigenesis, and apoptosis. ATF1 can regulate the expression of downstream target genes to affect cell physiological processes, which is related to soft tissue melanoma, and the relationship between melanoma and HF has been confirmed in previous literature [9–11].

Figure 11. (A–V) are scatter plots of immune cell/function and diagnostic gene correlations, respectively.

Figure 12. Exploration of the interaction between diagnosis-related genes and miRNA and TF. (A) The miRNA-mRNA interaction network of diagnosis-related genes. (B) The mRNA-TF interaction network of diagnosis-related genes.

Results qRT-PCR experimental verification

As shown in Figure 13, the expression levels of diagnosis-related genes MYD8, TNF, ATG5, CD8A, ENTPD1, IL17RA, NT5E, IL1R1, PRF1, IL-10 and CASP1 in the HF-cell model and control cell model. Among them, the p-value of MYD8, ATG5, ENTPD1, and IL1R1 expression in both groups was less than 0.05. MCAD, CD8A, IL17RA, NT5E, and CASP1 had p-value less than 0.01 in both groups. PRF1 and IL-10 had p-value less than 0.001 in both groups. The expression trends of these genes through qRT-PCR experiments were consistent with the results of differential expression analysis.

Figure 13. (A–K) are the mRNA expression levels of MYD8, TNF, ATG5, CD8A, ENTPD1, IL17RA, NT5E, IL1R1, PRF1, IL-10 and CASP1 by qRT-PCR. *p<0.05, **p<0.01, ***p<0.001.

Clustering algorithm

Autoencoders

Autoencoders (AE) is a deep neural network that consists of an encoder and a decoder. Both the encoder and decoder are composed of multilayer feedforward neural networks. They are connected by the bottleneck layer. The encoder and decoder are represented by Formula (2) and Formula (3), respectively.

$$z=f_{\text{encoder}}(x)$$ (1)

$$x^{\prime }=f_{\text{decoder}}(z)$$ (2)

Where z is the output of the encoder, which can be regarded as a reduced-dimension representation of data. x' is the output of the decoder. f_encoder and f_decoder are multilayer neural networks. In this paper, all encoder-based models are implemented using PyTorch. Among them, the loss function used by AE is MSE loss functions. For all encoder parts, the number of network layers was set to [10, 5, cluster_num]. Cluster_num represents the number of clusters. For all decoder parts, the network layer number was set to [10, 5, cluster_num]. For all models based on the self-encoder, epoch was set to 100 during training.

Denoising autoencoders

Different from AE, denoising autoencoders (DAE) construct partially damaged data by adding noise to the input data, and then restores it to the original input data by encoding and decoding. The newly generated $\tilde{x}$ can be expressed by the following formula.

$$\tilde{x}=q_D(\tilde{x}|x)$$ (3)

Where q_D represents random mapping and obeys the unit normal distribution $\text{N}(0,1)$. Its encoder and decoder can be represented by Formula (4) and Formula (5), respectively.

$$z=f_{\text{encoder}}(\tilde{x})$$ (4)

$$x^{\prime }=f_{\text{decoder}}(z)$$ (5)

K-means clustering algorithm

K-means clustering is a classical clustering algorithm, and its implementation steps are as follows: Firstly, k clustering centers are randomly selected. Then the distance from each sample point is calculated, and the cluster center will divide it into the nearest cluster, forming k clusters. Next, the centroid (mean) of each cluster is recalculated. Repeat the above process until the position of the center of mass no longer changes or reaches the set number of iterations. In this paper, the algorithm was implemented with default parameters based on the scikit-learn package of Python.

Algorithm evaluation index

In this paper, three evaluation indexes of clustering performance were introduced, including the sample contour coefficient (the value was between -1 and 1, the closer to 1, the better the clustering effect), the Calinski-Harabasz Index (the value was greater than 0, the better the clustering effect), and the Davies-Bouldin Index (the value was greater than 0, the closer to 0, the better the clustering effect). All of them were implemented by Python’s scikit-learn package.

Data acquisition

All the data in this paper came from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Specifically, we used GSE141910 data set (126 diseased samples and 240 control samples) as the training set, and the GSE116250 data set (50 diseased samples and 14 control samples) as the test set. The GSE141910 data set is derived from left ventricular free wall tissue harvested during cardiac surgery from HF subjects undergoing transplantation and from unused donor hearts with apparently normal function. Cold myocardial paralysis was perfused before cardiac resection to block contraction and prevent ischemic damage, and tissue specimens were frozen in liquid nitrogen. The GSE141910 data set comes from 64 samples of human left ventricular tissue.

The expression of ICDRGs before and after renal ischemia-reperfusion

The differential expression analysis of the GSE141910 data set was carried out using the “limma” package, and 8885 differentially expressed genes were obtained. We set Adj.P.Val<0.05 as the threshold for screening differentially expressed genes. Then, this paper collected 20 ICDRGs from the previous work and intersected them with differentially expressed genes to obtain intersection genes (DEICDRGs). The expression of ICDRGs in the diseased group and the control group in the GSE141910 data set was displayed as a box graph. In order to evaluate the correlation between ICDRGs, the Pearson correlation coefficient of DEICDRGs in the sample is calculated and visualized by “corrplot” in R software.

Enrichment analysis of different clusters

In this paper, the enrichment analysis of gene ontology (GO) and genome encyclopedia (KEGG) was realized by using the R package “clusterProfile.” And the analysis of gene set variation (GSVA) was realized by using R-packet “GSVA.” Among them, c2.cp.kegg.v7.4.symbols.gmt was downloaded from the MSigDB database and used to study the changes in biological signal pathways. The R package “ggplot2” was used to visualize the enrichment results.

Immunoassay

The ssGSEA algorithm was used to estimate the infiltration abundance of immune cells and the score of immune function in the sick group and the control group. Then, we also collected the genes related to immune inspection sites and explored the difference in expression levels of genes related to immune inspection sites between the two groups by box chart.

Construction and verification of HF-related diagnosis model

For the purpose of screening the diagnosis-related genes of HF, we adopt a small absolute contraction and selection operator (LASSO) and support vector machine recursive feature elimination (SVM-RFE). The R-packet “glmnet” was used to implement the LASSO algorithm in the case of ten times cross-validation. SVM-RFE was also implemented by R script under the condition of ten times cross-validation. Finally, we draw ROC curves of diagnostic genes and diagnostic models by using the R-package “pROC.”

Construction of nomogram

The study constructed the nomogram model based on diagnosis-related genes using the R package “rms”. Then the validity of the nomogram model was evaluated by the calibration curve. The clinical practicability of the nomogram model was evaluated by a decision curve. Finally, the high-risk probability stratification was predicted by the clinical impact curve when the population size was 1000.

Construction of interaction network

We conducted online analysis of hub genes using the NetworkAnalyst database (https://www.networkanalyst.ca/NetworkAnalyst/) to construct a transcription factor (TF) hub gene network and a miRNA hub gene interaction network. In this research, miRNA and TF interacting with diagnosis-related genes were queried by using the miRTarBase database (http://mirtarbase.cuhk.edu.cn/php/index.php) and TargetScan database (http://www.targetscan.org/vert_72/) respectively. Then, an interaction network between diagnostic-related genes and miRNA was constructed, the same between diagnostic-related genes and TF.

Experimental validation of diagnosis-related genes

In this paper, the mRNA expression of diagnosis-related genes was detected by qRT-PCR. Total RNA was isolated from the cell culture of H9C2 cells for HF and HL-1 cells for normal using the TRI Reagent. Five samples were in the HF group and five in the normal group. The mRNA expression level was calculated with the 2-ΔΔCt. All data were expressed by means ± SD, and the statistical differences between groups were tested by T-Test, and p < 0.05 indicated a significant difference. The primer sequences are shown in Table 1.

Data availability

The data used in the paper was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/).