Discovery of Drug-Responsive Phenomic Alteration-Related Driver Genes in the Treatment of Coronary Heart Disease

Introduction

The Xueyu Zheng (XYZ) phenome, also called blood stasis syndrome in traditional Chinese Medicine (TCM), is characterized by a set of six common symptoms in coronary heart disease (CHD), including two main symptoms (chest pain and chest distress) and four secondary symptoms (palpitations, purple or dark lips, purple or dark tongue, and unsmooth pulse).^1–3 Previous studies have reported the genetic correlations and causal associations between different phenotypes or symptom clusters and CHD.^4–6 While disease describes a pathological phenotype in modern medicine, Zheng defines the special pathological phenome, which was used in the TCM clinical practice to differentiate the complex illness with similar presentations at different stages during the course of pathological progression.⁷ The XYZ phenome is considered to be a central symptom cluster in patients with CHD,^8,9 which is a complex pathophysiological state characterized by decreased or impeded blood flow ^10,11. Previous studies have made great efforts to reveal the underlying genetic associations of the XYZ phenome in CHD,¹² including the associations of inflammatory- and immune-related genes, such as the genes encoding interleukin (IL)-8,¹³ FcγRIIIA,¹⁴ tumor necrosis factor-α, and IL-1,¹⁵ as well as those encoding platelet activation-related factors,¹⁶ such as platelet-activating factor acetylhydrolase^17,18 and the platelet-activating factor receptor (PAFR) gene.¹⁹ However, most previous studies have focused on the differential changes in a single gene or the gene expression profiles of patients with CHD with and without the XYZ phenome, thus ignoring the genetic associations among the XYZ phenome, its six component symptoms, and CHD.

A phenomics strategy emphasizing the simultaneous study of multiple phenotypes across biological scales may be useful, particularly if the high heritability of CHD, the XYZ phenome, and its component symptoms are due to a large number of genetic variants with small effects. This strategy has been successfully applied to detect the genetic associations of the 46 CHD risk factors,⁵ depression and related phenotypes,²⁰ neuropsychiatric disorders and cognitive impairment,^21,22 and different syndromes in male fertility,²³ as well as the genetic associations among psychiatric disorders, education, socioeconomic status, and brain phenome.²⁴ It has also been used to assess the polygenic risk score for Alzheimer’s disease in electronic health records.²⁵ Nevertheless, studies on the genetic associations in drug-responsive phenome alterations are rare. Moreover, there is little information on which genes drive the dynamic genetic variations and the associations among the XYZ phenome and its component symptoms.

In our previous study, we found that treatment with Danhong injection (DHI) could ameliorate the angina-specific health status assessed with the Seattle Angina Questionnaire (SAQ), and improve the symptom score of the XYZ phenome for 90 days in patients with CHD and XYZ phenome, which could be related to anticoagulation and regulation of cholesterol metabolism.³ However, the genetic associations among the XYZ phenome, its six component symptoms, and CHD remain unclear. Therefore, in this study, we further applied the phenomics approach to study the quantitative relationship between the XYZ transcriptome and phenome in CHD using Spearman correlation and the canonical correlation analysis (CCA). The phenomic alteration-related genes (PARGs) in the XYZ phenome and its component symptoms were screened after treatment with DHI at day 14 and day 30. Finally, the driver genes of the PARGs were detected from protein-protein interaction (PPI) using a personalized network control (PNC) algorithm to reveal the dynamic genetic variations in the XYZ phenome and its component symptoms when treated with DHI from day 14 to day 30.

Methods

Ethics, Consent and Permissions

The data of this study were obtained from an adaptively designed, randomized, multicenter, double-blind, placebo-controlled trial including 920 patients with CHD with the XYZ phenome (NCT01681316).³ The trial protocol³ was approved by the Central Institutional Review Committee of the Chinese People’s Liberation Army General Hospital (IRB no. [2012] Yao (025)). This trial complied with the principles of the Declaration of Helsinki. Written informed consent was obtained from all patients.

Study Design and Data Collection

Danhong injection (DHI), made by Shandong Danhong Pharmaceutical Co. Ltd., with quality control via high-performance liquid chromatography (HPLC) fingerprinting for over 90% similarity among the batches used in this trial) is a polycomponent drug, which contains five main compounds including Danshensu sodium, protocatechualdehyde, p-coumaric acid, rosmarinic acid, salvianolic acid B, exacted from two kind of herbal medicine, Danshen [Salvia miltiorrhiza] and Honghua [safflower].³

Symptomatic patients aged 18–70 years who were diagnosed with CHD and the XYZ phenome, with a score of at least 15 on the Chinese Medicine Symptom Scale of Xueyu Zheng phenome for angina, and who were classified as Canadian Cardiovascular Society (CCS) class II or III, were enrolled in this study. All patients voluntarily participated. All patients underwent 14 days of treatment with DHI, or placebo (0.9% saline). The follow-up period was 76 days. The trial details have been reported previously.³

Sixty-two of 920 patients agreed to provide serum samples for RNA sequencing (28 from Chinese People’s Liberation Army General Hospital (301 Hospital) and 34 from Xuanwu Hospital) at baseline (day 0), after treatment (day 14), and at the 2-week follow-up (day 30), including 41 from the DHI group and 21 from the control group.

The status of the XYZ phenome in each patient was assessed according to the XYZ scale, which is a quantitative doctor-reported outcome scale on each symptom in the XYZ phenome, at days 0, 14, and 30. The scale includes the following items: (a) chest pain (0–10); (b) chest distress (0–10); (c) palpitations (0–5); (d) purple or dark lips (0–5); (e) purple or dark tongue (0–5); and (f) unsmooth pulse (0 or 5).³

Thirty-two patients in the DHI group and 19 patients in the control group with no missing data on the XYZ scale and mRNA sequencing profile were selected for further phenomic analysis.

Next-Generation Sequencing and Bioinformatics Analysis

Total RNA for the transcriptome analysis was isolated from the serum samples of the patients and sequenced using the Illumina HiSeq2000 platform according to the manufacturer’s instructions. We evaluated the raw data from two aspects: 1) filtering contaminants and 2) initial judgments on the data. Sequences were aligned to the human genome reference (hg19) sequence, Rfam, and the NCBI database using Bowtie 2²⁶ (one mismatch allowed) with default parameters. The raw Illumina reads were first preprocessed by cutting the 3’adapter sequence using the program fastx_clipper from the FASTX-Toolkit. Reads shorter than 18 nts after clipping were removed. The remaining reads are reduced to unique reads and their frequency per sample to make the mapping steps more time efficient.^27,28

Spearman Correlation Coefficient and Canonical Correlation Analysis

First, the changes in gene expression (δE_n), as well as the variations in total XYZ scale (δZ) and the score of each symptom in XYZ (δS_n, n = 6), were respectively obtained on day 14 compared with baseline and on day 30 compared with day 14. Then, Spearman correlation coefficient (ρ) (Equation 1) was calculated between the δE_n of each gene and the δZ or δS_n at day 14 and day 30 using RStudio (psych package). Spearman’s rank correlation coefficient is denoted as ρs for a population parameter and as rs for a sample statistic.²⁹ According to relevant published literature, 0.4 represents a moderate association.³⁰ An absolute correlation coefficient of ≥0.4 and a p value of <0.05 were used to determine the PARGs. PARGs identified in DHI group but not in the control group were considered as drug-responsive PARGs. The x and y are two random variables (can also be regarded as two sets), the number of elements of which is represented by N, and the i (1 ≤ i ≤ n) values of the two random variables are represented by x_i and y_i, respectively.

(1)

The relationship between δS_n and δZ was calculated using the CCA (Equation 2) in R 3.6.0 (mvstats package). The CCA is a useful dimension-reduction technique to explore the relationship between two sets of variables.³¹ The u and v represent the linear combination between the respective features of x and y, and a and b are the weights that make u and v have the greatest correlation.

(2)

Gene Enrichment Analysis

The Gene Ontology (GO) and pathway enrichment analyses were conducted using Metascape (https://metascape.org/).³² The species was selected as “Homo sapiens (human)”. The settings used were as follow; Min Overlap: 3, P value Cutoff: 0.05, and Min Enrichment:1.5. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched from the drug-responsive PARGs in the XYZ phenome were mapped to each symptom to detect their associated pathways.

Detection of Driver Genes for Phenomic Changes

To identify driver genes that may contribute to phenomic alterations in the XYZ phenome in patients with CHD, a PPI network of PARGs involved in the enriched KEGG pathways was constructed using STRING 11.0 (https://cn.string-db.org/).³³ The minimum set of driver nodes, which direct the network of a specific condition to the desired control aim, can be identified from a network controllability perspective. Therefore, based on the genetic data of the patients, personalized driver genes (PNGs) associated with phenomic changes in each patient were first detected using PNC algorithm (https://github.com/NWPU-903PR/PNC) programmed in MATLAB software (The MathWorks Inc., Natick, MA, USA).³⁴ Those identified as PNGs in more than 15 patients treated with DHI were considered to be drug-responsive driver genes affecting phenomic changes in XYZ from day 14 to day 30. The distribution of the driver genes in the PPI network of the PARGs was visualized using Cytoscape 3.6.0 (https://cytoscape.org/).

Furthermore, to validate the controllability of the driver genes, the node importance of these genes in the PPI network of the PARGs, including the degree centrality, betweenness centrality, and closeness centrality, as well as the characteristic path length for the PPI network of the PARGs after sequentially deleting these driver genes, were calculated using the network analyzer application in Cytoscape 3.6.0. We also compared the Spearman correlation coefficients of these driver genes for phenomic changes in XYZ between the two groups to assess their specificity for the response to DHI. Finally, the relationships among the symptoms, driver genes and related pathways, and GO biological processes were further analyzed to reveal the core mechanism of temporal phenomic changes from day 14 to day 30.

Results

Improvement in the XYZ Phenome in Patients with CHD Treated with Danhong Injection (DHI)

The distribution of demographic and clinical characteristics was reasonably well balanced between the two groups (32 patients in DHI group and 19 patients in the control group) (Table 1). In these 51 cases, although there were no differences between the two groups in terms of the XYZ phenome changes or in the six symptoms, the scores of the XYZ phenome and the six symptoms improved significantly at day 14 and day 30 from baseline (p < 0.001) (Figure 1a and b). The correlation coefficient between the change in each symptom and the XYZ score was >0.6, except for “unsmooth pulse”. The changes in the main symptoms of chest pain and chest distress were strongly correlated (0.98 and 0.92, respectively) with the variation in the XYZ phenome (Figure 1c).

Table 1 Baseline Characteristics Between the Two Groups from the Selected Population

Figure 1 Changes of the XYZ phenome in CHD patients after Danhong injections. (a) On day 14 vs day 0, all XYZ phenome symptoms were improved except for unsmooth pulse (p < 0.01). (b) On day 30 vs day 0, all XYZ phenome symptoms were improved except for unsmooth pulse (p < 0.01). (c) Heatmap of correlation coefficients between symptoms and XYZ.

Phenomic Alterations Related Genes (PARGs) Associated with the Improvement in the XYZ Phenome in Patients with CHD

There were 544 and 318 genes correlated with XYZ phenome changes at day 14 and day 30 of Danhong injections, respectively (Figure S1). After comparing these results with the changes in the XYZ phenome in the control group, 525 and 309 genes were considered as drug-responsive PARGs from day 0 to day 14 and from day 14 to day 30 (Figure 2a, Table S1 and Table S2), respectively, seven of which were observed at both time points (UBE2F, ALKBH3, ARHGEF18, ILF2, CHMP5, RAB6A, LOC646938) (Figure 2b). The top 10 PARGs at day 14 and day 30, most of which had an absolute correlation coefficient (ρ) of >0.57, are listed in Table 2. The genes related to the changes in five symptoms are shown in Figure 2a and Table 2, with the tendency that the PARGs at day 14 from baseline were mostly related to the two main symptoms (chest pain and chest distress) and those at day 30 from day 14 were mostly related to the three secondary symptoms (palpitations, purple or dark lips, and purple or dark tongue).

Table 2 Top ten Genes That Significantly Related to Each XYZ Symptom and XYZ Phenome (P<0.05)

Figure 2 Molecular mechanism of the improvement in the XYZ phenome in patients with CHD after treatment with Danhong injections. (a) Number of PARGs in the XYZ phenome and its component symptoms. (b) Comparison of genes, pathways, and GO biological processes between the day 14 and day 30. (c) Relationship between the involved pathways and the alterations in the XYZ phenome at day 14 and day 30.

GO biological Processes and Pathways Related to Changes in the XYZ Phenome in CHD

The 525 drug-responsive PARGs identified at day 14 from baseline were associated with 48 pathways and 543 GO biological processes, while the 309 drug-responsive PARGs identified from day 14 to day 30 were enriched in 35 pathways and 312 GO biological processes (Figure 2b, Table S3–14). Nine pathways and 71 GO biological processes were common between the two time points (Figure 2b and c).

As shown in Figure 2c, the common mechanism for the variation in the XYZ phenome in CHD may be related to regulation of PI3K-Akt signaling pathway (Figure 3a), endocytosis (Figure 3b), and HIF-1 signaling pathway(Figure 3c). We also identified pathway variations from immune response (Th1 and Th2 cell differentiation, Th17 cell differentiation) and endothelial function (oxidative phosphorylation) at day 14 to atherosclerosis (Ras signaling pathway) and platelet activation (Rap1 signaling pathway) at day 30.

Figure 3 The overlapping KEGG pathways in patients treated with Danhong injections between the day 14 and day 30. (a) PI3K-Akt signaling pathway; (b) endocytosis; and (c) HIF-1 signaling pathway. The gene highlighted in red represents the day 14 and the gene highlighted in blue represents the day 30.

Biological Interactions Responsible for the Temporal Changes in XYZ Phenome Symptoms

As shown in Figure 4a and c, at day 14, 286 genes and 29 pathways were common between the two main symptoms (chest pain and chest distress). At day 30, only 20 genes were common between the two main symptoms and no common pathways were found (Figure 4b and d). The common genes and pathways related to the secondary symptoms (palpitations, purple or dark lips, and purple or dark tongue) showed the opposite tendency from day 14 to day 30.

Figure 4 The internal genetic relationship of the symptoms of the XYZ phenome. (a and c) The common PARGs and pathways among five component symptoms at day 14; (b and d) the common PARGs and pathways among five component symptoms at day 30; (e and g) comparison of PARGs and pathways between five component symptoms with XYZ phenome at day 14; (f and h) comparison of PARGs and pathways between five component symptoms with XYZ phenome at day 30.

To detect the genetic associations among the symptoms of the XYZ phenome, we found that the two main symptoms of chest pain and chest distress accounted for a large proportion of the XYZ phenome (genes: 58.67% and 48.57%, respectively; pathways: 52.08% and 81.25%, respectively) at day 14 (Figure 4e and g), but this proportion was dramatically decreased at day 30 (genes: 11.65% and 23.62%, respectively; pathways: 2.86% and 45.71%, respectively) (Figure 4f and h). The common pathways of the main symptoms included the HIF-1 signaling pathway, the Notch signaling pathway, and the chemokine signaling pathway (Figure 5).

Figure 5 The internal relationship among the XYZ phenome, internal symptoms, and involved pathways. Day 14 is marked in red and day 30 is marked in blue. The red font indicates the pathways associated with the XYZ phenome.

Compared with the main symptoms, the secondary symptoms showed the opposite phenomenon. At day 14, three secondary symptoms accounted for a small proportion of the XYZ phenome, including palpitations (genes: 1.71%; pathways: 4.17%) and purple or dark lips (genes: 8.57%; pathways: 0%). However, the proportion of genes and pathways associated with palpitations and purple or dark lips increased at day 30 (genes: 11% and 28.8%, respectively; pathways: 45.71% and 62.86%, respectively) (Figure 4e-h). The common pathways included the Ras signaling pathway and the Apelin signaling pathway (Figure 5).

Driver Genes Related to the PPI Network of PARGs

The PPI network of 189 PARGs involved in the enriched KEGG pathways (Figure 5) was established (Figure 6a). The sub-PPI network at day 14 (Figure 6b) was denser than at day 30 (Figure 6c). Five driver genes, including four genes related to the day 14 subnetwork (Figure 6b) and one gene related to the day 30 subnetwork (Figure 6c), were detected, namely EIF2S1, HSPA8, UBA52, AURKA, and VCP. From the node distribution of the network, it was found that all five driver genes were located at the center of the network or the temporal subnetwork (Figure 6).

Figure 6 The PPI network of PARGs. (a) The PPI network of 189 PARGs. (b) The sub-PPI network at day 14. (c) The sub-PPI network at day 30.

Furthermore, the four driver genes at day 14 were associated with four symptoms (chest pain, chest distress, purple or dark tongue, and purple or dark lips). AURKA was specific to day 30 and was associated with chest distress (Figure 7).

Figure 7 Symptoms, pathways, and GO biological processes related to driver genes. The red and blue bars in the middle of the graph represent the day 14 and day 30, respectively. The different colored lines indicate the different genes associated with symptoms, pathways, and GO biological processes. The black font indicates the pathway name, and the blue font indicates the GO biological process.

Compared with the control group, the correlation coefficients between the driver genes and the phenomic changes were significantly higher in DHI than in the control group (Figure 8a). The degree centrality, betweenness centrality, and closeness centrality of the five driver genes were higher than the average values of the PPI network (Figure 8b). Furthermore, after sequentially deleting these driver genes, the characteristic path length of the PPI network of PARGs increased and the number of shortest paths decreased (Figure 8c). This proved that these five driver genes were important nodes in the PPI network.

Figure 8 Validation of driver gene controllability. (a) The correlation coefficients between the driver genes and phenomic changes in the Danhong injection group and the control group, **p<0.01, *p<0.05. (b) The degree centrality, betweenness centrality, and closeness centrality of the driver genes. (c) The changes in the characteristic path length and the number of shortest paths for the PPI network of the PARGs with sequential deletion of driver genes.

Discussion

The XYZ phenome is an essential symptom cluster in patients with CHD. In our study, we applied a phenomics strategy to screen drug-responsive PARGs for the XYZ phenome and its component symptoms by correlation analyses. We detected five driver genes from the PPI network of PARGs using a structure-based PNC network control model, and we uncovered the factors driving dynamic genetic variations in the XYZ phenome and its component symptoms.

We found that the genes associated with changes in the XYZ phenome were mostly related to the regulation of endoplasmic reticulum (ER) stress-mediated apoptosis and autophagy, which could suppress the progression of atherosclerosis³⁵ (Figure 9). The driver gene EIF2S1, associated with four symptoms (chest pain, chest distress, purple or dark lips and purple or dark tongue), encoding eukaryotic initiation factor 2 alpha (eIF2α), contributes to reduced protein translation and thereby the overload of unfolded and misfolded proteins in the ER lumen, which is a key protein in PERK/eIF2α/ATF4 signaling pathway highly correlated with ER stress.³⁵ UBA52 (ubiquitin A-52 residue ribosomal protein fusion product 1), a major source of ubiquitin protein, associated with three symptoms (chest pain, chest distress, and purple or dark lips), is crucial in the ubiquitylation of the heat shock protein family proteins,³⁶ and the latter ones are unfolded protein response (UPR)-associated proteins to induce ER stress and and trigger the ER apoptotic pathway.³⁷ VCP (valosin containing protein), a newly identified calcium-associated ATPase, reduces ER stress with cardioprotective effects by inhibition of ATP consumption.³⁸ AURKA (aurora kinase A), which induces activation of the PI3K/Akt signaling pathway,³⁹ and the latter could inhibit ER stress.⁴⁰ PI3K/Akt signaling pathway was found to be associated with protective autophagy.⁴¹ It is reported that PI3K/Akt signaling, as a key component of RISK pathway, is involved in cardiac protection. The activation of endogenous Akt plays a critical role in controlling cell survival by resisting ER stress-induced cell death.⁴² There is growing evidence that ER stress plays an important role in various stages of atherosclerosis initiation and progression.⁴³ It can be both protective and pro-apoptotic. ER sensor proteins, such as IRE1 and PERK, are involved in the UPR pathway regarding adaptive autophagy and pro-apoptosis. PERK is activated by autophosphorylation after dissociation from GRP78. Activated PERK blocks cellular protein synthesis by phosphorylating EIF2S1, which is required for the induction of cap-dependent translation. IRE1, as an important ER stress sensor, is activated by dimerization and autophosphorylation after dissociation from GRP78.⁴⁴ HSPA8 (heat shock protein family A (Hsp70) member 8) protects against endogenous or exogenous production of ROS, which promotes the progression of coronary artery disease.⁴⁵ HSPA8 also induces ER stress and UPR, in which Perk-p-eIF2α-Atf4 likely functions as a major activated pathway.⁴⁶ These observations support the utility of autophagy induction as a therapeutic strategy to neutralize ER stress and apoptosis in CHD (Figure 9). While proper autophagy removes aggregated proteins to maintain ER homeostasis and delay apoptosis under ER stress, excessive autophagy triggered by sustained ER stress may ultimately lead to apoptosis. Notably, sustained ER stress can also disrupt autophagy, leading to inadequate cellular clearance, which can lead to inflammatory, cardiovascular, and metabolic diseases.⁴⁷

Figure 9 Molecular mechanism of the phenomic alterations in XYZ in patients with CHD treated with Danhong injections from the day 14 to day 30. Created with BioRender.com.

We first used the phenomics strategy to analyze the genetic associations between XYZ phenome in CHD and its component symptoms of altered drug responsiveness. Particularly, the molecular mechanisms of phenomic alterations in response to drug interventions were uncovered in our study. The phenomics strategy enables the transition from traditional genome-wide association studies to phenome-wide association (PheWAS) studies. For example, the identification of multiple disease genes or proteins at the genomic and proteomic levels provides a solid platform for the systematic integration (or aggregation) of new definitions of hypertension phenotypes across the genomic/proteomic/metabolomic spectrum.⁴⁸ Since 2010, PheWAS has been used to investigate whether one genotype that correlated to one phenotype was also associated with other phenotypes.⁴⁹

Our study has some limitations that should be noted. We only analyzed the intrinsic relationship in the XYZ phenome of CHD, and we did not exclude the influence of CHD on XYZ. In addition, our study was a pilot study of 32 patients. Thus, drug-responsive PARGs of the XYZ phenome should be further validated in further independent large cohort studies.

Conclusion

Our study identified the drug-responsive PARGs of the XYZ phenome in patients with CHD treated with Danhong injections, which were mostly related to the regulation of ER stress-mediated apoptosis and autophagy. Our results provide an exemplary strategy to investigate the genetic associations among common phenomes and their component symptoms.