Rapid Discovery of Potential Drugs for Osteonecrosis of Femoral Head Based on Gene Expression Omnibus Database and Connectivity Map.

OBJECTIVE: To use Gene Expression Omnibus (GEO) database coupled with Connectivity Map (CMap) databases to screen potential therapeutic drugs for osteonecrosis of femoral head (ONFH) rapidly.
METHODS: Raw genetic data with accession number GSE74089 that contained eight hip articular cartilage specimens from four ONFH patients and four healthy controls were obtained from the Gene Expression Omnibus (GEO) database and were then integrated using R to identify differentially expressed genes (DEGs). Subsequently, to identify several potential small molecular compounds that were most strongly negatively correlated with ONFH, a search query of DEGs was explored by using CMap.
RESULTS: Filtering revealed 1937 DEGs with log (fold-change) ≥1 and adjust P value < 0.001. Finally, a network of candidate targets for ONFH with 135 nodes and 660 edges was constructed through network topology analysis, including 96 up-regulated genes and 39 down-regulated genes. Several significant gene functions and signaling pathways associated with pathological processes of ONFH were identified via gene enrichment analysis. Based on the CMap database, some potential small molecular components that may be possible to counteract the effects of molecular signal imbalance for ONFH were identified. Neostigmine bromide with low CMap score and P value and specificity score was predicted to be the most candidate compound, involved in the "positive regulation of stem cell proliferation," "regulation of protein autophosphorylation," "VEGF signaling pathway," and "ECM-receptor interaction."
CONCLUSIONS: The GEO and CMap databases can be effectively used in understanding the molecular changes in ONFH and provide a systematic manner to identify potential drugs for ONFH prevention and treatment. However, additional clinical and experimental research of the candidate compound is warranted.

Introduction

Osteonecrosis of femoral head (ONFH) is one of the most refractory orthopedic diseases, which has become a common public health problem1. ONFH is a multifactorial progressive disease in which greater osteonecrosis and collapse lead to loss of joint function or disability2, 3. The mechanisms underlying ONFH are not particularly clear and there is no specific drug for ONFH nowadays. Total hip arthroplasty (THA) is recognized as the ultimate treatment. Although the surgical technique is mature, its high cost and complications are still unacceptable to many people. The most important point is that this treatment is not suitable for early patients. Therefore, for early ONFH patients, timely and accurate diagnosis and conservative treatment are very important. But unfortunately, according to data, the early diagnosis and remedy for ONFH remains limited4. Therefore, the identification of suitable targeted molecular drugs will significantly benefit ONFH patients, particularly young patients. Due to the huge investment and long clinical observation needed for new drug research and development, the progress of new drug research and development for ONFH has been slow and rarely reported.

The Gene Expression Omnibus (GEO) database, a National Center for Biotechnology Information (NCBI) database for gene expression and hybridization array data, contains a wide assortment of high‐throughput experimental data for various diseases5. The GEO database is available for querying, downloading, and analyzing gene expression data for ONFH. The Connectivity Map (CMap) database, a genome‐wide expression database of gene expression profiles of multiple human cells, can be used to compare disease gene expression patterns and analyze molecular mechanisms of different pathophysiological states6, 7, 8. In addition, researchers can screen a ranked list of small targeted molecules based on connectivity scores. If there exists a significant negative correlation between a disease‐related gene expression pattern and the functional profile of a compound, we can preliminarily speculate about the existence of a link between the small molecule in question and its therapeutic effects on the disease. Using the CMap database, we can quickly identify new indications for existing drugs, a process known as drug repositioning or drug rediscovery.

There have been several successful applications of CMap to drug rediscovery. Based on prior experiences, we designed and performed our research, which we present as a flowchart (Fig. 1). Data was collected from the hip articular cartilage specimens of four ONFH patients and four healthy controls. This data was analyzed using the GEO database; a series of differentially expressed genes (DEGs) was then identified using R. ClueGO was used for the Gene Ontology (GO) analysis including the three categories of cellular components, molecular functions and biological processes and KEGG pathways. By comparing characteristic DEGs with shared expressed genes in the CMap database, small molecular compounds or drugs with negative correlations with ONFH were obtained.

Figure 1

The flowsheet for the DEGs query signature and the discovery of candidate molecular compounds for ONFH.

The purpose of this study is: (i) mining ONFH related chips through GEO database to obtain differentially expressed genes; (ii) the core differential genes of ONFH were obtained by network topology analysis. Afterwards, the related physiological and pathological processes and signaling pathways were obtained by enrichment analysis; and (iii) potential drugs for ONFH were obtained from CMap database. It will provide some reference for the follow‐up clinical, basic research and the development of new drugs.

Materials and Methods

Data Source

Only one Agilent microarray dataset related to ONFH on Homo sapiens, which was identified by the accession number GSE74089, and was available in the GEO database on 6 December 2018 (https://www.ncbi.nlm.nih.gov/geo/). The data of accession number GSE74089, which was based on the GPL13497 [Agilent‐026652 Whole Human Genome Microarray 4×44K v2], was contributed by Ruiyu L et al., which contained eight hip articular cartilage specimens that were collected from four ONFH patients and four healthy controls. Unfortunately, the descriptive text of the chip data does not mention whether the femoral head specimens included are from traumatic, steroid‐induced, or alcoholic ONFH patients, nor does it describe the stages.

Differentially Expressed Genes Analysis

Using the marray package in R to conduct and process probe level data in txt files, background correction was carried out by method of normexp, followed quantile normalization and probe summarization, and the gene matrix related data was extracted. DEGs between ONFH patients and healthy controls were identified by the limma package in R, which were displayed by the volcano plot. Data parameters were filtered via variance analysis with applicable thresholds of P value < 0.001 and fold‐change ≥ 2. Principle Component Analysis (PCA) were performed in R environment9 (prcomp function) using log2‐transformed probe set intensity of all genes and visualized using scatterplot3d package, after removing unwanted variants.

Protein–Protein Interaction Network Construction and Analysis

Based on the DEGs identified above, the Protein–Protein Interaction (PPI) network was formed by STRING database with combined score > 0. 4 (Version 10.5, ELIXIR, Europe, https://string-db.org/)10. And then, the network that only connected nodes were retained was visualized and analyzed by Cytoscape (Version 3.6.1, Cytoscape Consortium, U.S), and the degree measure of each node was calculated by Network Analyzer. The greater the degree value, the more important the node was in the network.

Referring to the previous study of Zhang and others,11, 12 we determined the extent of all the nodes of node more than twice the intermediate nodes as key targets. In order to explore more molecular complexes and specific regulatory relationship of proteins in above PPI network, we performed cluster analysis through MCODE (Version 1.5.1, The Bader Lab of University of Toronto, Canada)13 that calculates node information that included neighbors and density through each node in the network graph, and that builds the functional modules for clustering with K‐Core. Data parameters was set with applicable thresholds of K‐Core > 5.

GO and KEGG Pathway Enrichment Analysis

To investigate the main functional mechanisms of ONFH, the topological candidate DEGs were enriched by using ClueGO14. The results were divided into the following parts: biological processes, cellular components, molecular functions, and the KEGG pathway, which were displayed by the bubble plot using R.

To display and visualize the relationship between a list of candidate genes and terms, as well as the logFC of the genes, the GOChord plotting function was implemented by the GOChord package and GOCluster package in R, and only genes which are assigned to at least one process could be displayed.

CMap Analysis

To investigate potential drugs that may be possible to reverse the effects of molecular signal imbalance for ONFH, the topological candidate DEGs profiles were compared with the CMap database (Update 12 September 2017, https://portals.broadinstitute.org/cmap/index.jsp). Since the CMap database is compiled using Affymetrix U133A array, and application in our study was an Agilent‐026652 Whole Human Genome Microarray 4×44K v2, we needed to converted the chip by Affymetrix Batch Query (https://www.affymetrix.com/analysis/netaffx/batch_query.affx?netaffx=netaffx4_annot). We divided genes into files with “up” and “down” tags based on whether the genes were up‐ or down‐regulated, which were displayed by the volcano plot. And then we saved our results in “.grp” format as query terms. Detailed and permuted results were then exported to Excel by uploading a signature of “up” and “down” tags to the CMap Web Service and providing an appropriate description. In addition, the connectivity scores, P values, and specificities of relevant small molecular compounds or drugs were shown. Small molecules that were highly negatively correlated with the ONFH signature were screened out and might have the potential to treat ONFH.

RESULTS

Identification of DEGs in ONFH

Patients were clustered into two distinct clusters based on PCA: one cluster exclusively had samples from ONFH patients, and the other group had samples from healthy people (Fig. 2A). We identified 1937 DEGs between ONFH patients and health (fold‐change >2.0, q value <0.001), including 1223 up‐regulated and 714 down‐regulated genes (Fig. 2B).

Figure 2

PCA and Volcano plot shows the up‐ and down‐regulated genes. (A) PCA of four ONFH and four health samples showed clear distinction between ONFH and health patient. Red cycle represents the ONFH samples while green square represents the health samples. (B) The Volcano plot showed up‐regulated genes and down‐regulated genes between ONFH and health patient. The horizontal axis represents the value of Log2(Fold Change), the vertical axis represents the value of ‐Log10(P value). Red represents the up‐regulated gene while green represents the down‐regulated gene.

Protein–Protein Interaction

By integrating DEGs with combined score >0.4, PPI network was formed by STRING database, involving 645 nodes (DEGs) and 1529 edges, accounting for 33.30% of all DEGs. The PPI network was visualized in Cytoscape, in which several disconnected nodes were discarded in the STRING database, containing 407 nodes and 1157 edges. To analyze the complex network, we identified the extent of all 407 nodes of node more than twice the intermediate nodes referring to a previous study of Zhang and others11. Thus, we constructed a network of candidate targets for ONFH that had 135 nodes and 660 edges, including 96 up‐regulated and 39 down‐regulated genes, which were shown in heatmap from red to green (white in the middle; Figs 3, 4). In this network, DEGs which have high degree value and absolute value of fold change, included VEGFA (degree, 68; Log2(FC), 1.66), JUN (degree, 54; Log2(FC), 2.53), FGF2 (degree, 53; Log2(FC), 3.01), FN1 (degree, 34; Log2(FC), 2.47), CDH1 (degree, 33; Log2(FC), −2.31), PTGS2 (degree, 32; Log2(FC), 2.64), COL1A1 (degree, 29; Log2(FC), 5.57), ABL1 (degree, 29; Log2(FC), 2.42), COL1A2 (degree, 26; Log2(FC), 2.73), PRKCA (degree, 24; Log2(FC), 2.50), LUM (degree, 23; Log2(FC), 2.28), ITGAM (degree, 23; Log2(FC), −2.91), COL3A1 (degree, 21; Log2(FC), 3.66).

Figure 3

Identification of candidate DEGs for ONFH. The nodes representing candidate up‐regulated DEGs are shown as red circles and the down‐regulated DEGs are presented as green circles. The colors of the nodes are illustrated from red to green (white in the middle) in descending order of Log2(FC). The sizes of the nodes are illustrated from small to big in ascending order of degree values.

Figure 4

Heatmap shows the up‐ and down‐regulated genes. The horizontal axis represents the sample name. Red represents the up‐regulated gene while green represents the down‐regulated gene. The colors of the nodes are illustrated from red to green (black in the middle) in descending order of Log2(FC).

To analyze the subunits of the complex and their interactions, cluster analysis were performed by MCODE. Two modules were extracted from the PPI network through MCODE analysis with K‐Core > 5 (Fig. 5). One cluster performed included 13 nodes and 76 interactions (cluster rank 1; Score 12.667), the other cluster performed included 28 nodes and 82 interactions (cluster rank 2; Score 6.074).

Figure 5

The modules were extracted from the PPI network through MCODE analysis. The colors of the nodes are illustrated from red to green (white in the middle) in descending order of Log2(FC). The sizes of the nodes are illustrated from small to big in ascending order of degree values. K‐Core >5.

GO and KEGG

To further clarify the possible roles and relationship of the 135 DEGs for ONFH, we used ClueGO for enrichment analysis. Specifically, we obtained cellular components related with proteinaceous extracellular matrix, extracellular matrix, collagen trimer, extracellular matrix component, fibrillar collagen trimer, plasma membrane protein complex (Fig. 6A), and molecular functions related with platelet‐derived growth factor binding, growth factor binding, glycosaminoglycan binding, growth factor receptor binding, heparin binding, vascular endothelial growth factor receptor binding (Fig. 6B), and biological processes related with blood vessel development, extracellular matrix organization, circulatory system development, blood vessel morphogenesis, angiogenesis, and regulation of cell adhesion (Fig. 6C). The KEGG pathways were PI3K‐Akt signaling pathway, Focal adhesion, Pathways in cancer, AGE‐RAGE signaling pathway in diabetic complications, and ECM‐receptor interaction (Fig. 6D).

Figure 6

Enrichment analysis of 135 DEGs for ONFH in the bubble plot. (A) The bubble plot of enriched cellular component terms of DEGs. (B) The bubble plot of enriched molecular function terms of DEGs. (C) The bubble plot of enriched biological terms of DEGs. (D) The bubble plot of enriched KEGG pathway terms of DEGs. The colors of the nodes are illustrated from red to green in descending order of ‐Log10(P Value). The sizes of the nodes are illustrated from small to big in ascending order of gene counts. The horizontal axis represents the gene ratio, the vertical axis represents the GO or KEGG terms.

To visualize smaller subsets of high‐dimensional data, the DEGs were selected from top 10 enriched biological process terms and KEGG pathway terms in the chord and cluster plot, which was plotted with GOplot package in R (Figs 7, 8). Chord plot could display the relationship between a list of selected genes and terms, as well as the logFC of the genes. Hierarchical clustering is a common method of gene expression analysis. This study uses supervised clustering analysis. Genes are grouped together according to their expression patterns. Therefore, clusters may contain a group of genes that are related to common regulation or function and establish an organic relationship between differentially expressed genes and signaling pathways.

Figure 7

The relationship between top 10 enriched biological process terms and DEGs. (A) The relationship is represented by the chord plot. (B) The relationship is represented by the cluster plot. The colors of the nodes are illustrated from red to green in descending order of logFC. The genes were ordered according to their logFC values.

Figure 8

The relationship between top 10 enriched KEGG pathway terms and DEGs. (A) The relationship is represented by the chord plot. (B) The relationship is represented by the cluster plot. The colors of the nodes are illustrated from red to green in descending order of logFC. The genes were ordered according to their logFC values.

Potential ONFH Drugs Predicted by CMap

In order to identify compounds with molecular features that are capable of reversing ONFH defects, we upload above DEGs into the CMap database. Ranking based on negative connectivity scores (−0.875 to −0.7) was used to reveal the top 10 small molecular compounds that had complete P value and specificity score and that may counteract the observed gene expression pattern for ONFH and may therefore be the most promising novel candidates for ONFH treatment (Table 1).

Table 1

Top 10 small molecules with complete P value and specificity score ranking based on negative connectivity scores

Rank	Compound name	Cell line	Mean CMap Score	n	P value*	Specificity
4	neostigmine bromide	PC3	−0.369	4	0.00229	0.0115
6	ioxaglic acid	HL60	−0.582	3	0.00389	0.0053
11	dinoprostone	PC3	−0.278	4	0.0066	0.0067
14	tridihexethyl	PC3	−0.433	4	0.00869	0.0152
15	cefapirin	PC3	−0.411	4	0.00977	0.0055
17	rilmenidine	PC3	−0.199	4	0.01723	0.0274
20	harmol	PC3	−0.282	4	0.01884	0.0929
22	securinine	PC3	−0.501	4	0.02109	0.1266
23	spaglumic acid	PC3	−0.588	2	0.02284	0.0591
26	lomefloxacin	PC3	−0.431	6	0.02352	0.0359

Compound name: the name given to a perturbagen; cell: cell line; Mean CMap Score: the mean connectivity score, a combination of the up score and the down score. A high negative connectivity score indicates that the corresponding perturbagen reversed the expression of the query signature; n: number of repetitive samples; specificity score: the extent of connectivity found between the two groups is unexceptional and/or the candidate compounds have different biological effects.

The probability of the enrichment of a set of instances in the total set of instances by chance upon execution of a query.

Unexpectedly, we found that neostigmine bromide, an inhibitor of N‐choline receptor by cholinesterase on excitatory skeletal muscle, was highly ranked (rank 4) in the ranked list and had low P value and specificity score. The value of amplitudes(a) represented a measure of the extent of differential expression of a given probe set. Here, a = 0 indicates no differential expression, a > 0 indicates increased expression upon treatment, and a < 0 indicates decreased expression upon treatment. Based on amplitude for each probe set in a pair of tag lists in a given instance, the genes that could be reversed by neostigmine bromide are analyzed (Table 2). To indicate the target molecular mechanisms of the molecule, we conducted a functional analysis of these potentially reversed genes in ONFH (Fig. 9). The genes reversed by neostigmine bromide were mostly involved in the “positive regulation of stem cell proliferation,” “regulation of protein autophosphorylation,” “VEGF signaling pathway,” and “ECM‐receptor interaction”.

Table 2

DEGs reversed of neostigmine bromide in ONFH based on amplitude

Direction	Genes
Decrease	EEA1, TCF7L2, NOV, FZD1, ID2, PLOD2, PLA2G4A, VCAN, PTGFR, RHBDF2, CAV1, EGR2, NT5E, SMUG1, VEGFA, CD55, OMD, GLI3, COL5A1, SULF1, ADM, S100A4, LOX, ANK3, PDGFC, ADAM10, VEGFC, THY1, COL5A2, GJA1, ADCY7, COL6A1, ITGB8, FN1, P4HA2, KDELR3, CORO1C, PRKCA, PNP, RND3, CDK6, PLAT, UBE2D2, HNRNPU, COL6A3, STX1A
Increase	SPTB, CDH1, KALRN, KRT7, POLE, RARA, KRT8, HLA‐DQB1, SLC11A1, A2M

Figure 9

Enrichment analysis of DEGs reversed of neostigmine bromide in the bubble plot. (A) The bubble plot of enriched GO analysis terms of DEGs. (B) The bubble plot of enriched KEGG pathway terms of DEGs. The colors of the nodes are illustrated from red to green in descending order of ‐Log10(P Value). The sizes of the nodes are illustrated from small to big in ascending order of gene counts. The horizontal axis represents the gene ratio, the vertical axis represents the GO or KEGG terms.

Discussion

At present, ONFH is one of the most common chronic progressive diseases in the world, which leads to osteonecrosis and collapse, resulting in damage to hip function and permanent disability. There are many previous studies on etiological mechanism system, such as lipid metabolism disorder, coagulation circulatory disorder, intramedullary hypertension, cell dysfunction, and so on. According to the mechanism found at present, the corresponding drugs have been developed and have played a certain role in clinical practice. However, with the changes of life, environment, and diet, the treatment of above drugs for ONFH were gradually diluted. So, it is urgent to research or discover effective therapeutic strategies for ONFH. In addition, pharmacoeconomics should be comprehensively considered. Given the aforementioned situation, we identified and assessed candidate molecule using high‐throughput transcriptome technology coupled with Connectivity Map, a drug repositioning tool based on systematic analysis of transcriptomics data15. Many successful applications of drug repositioning and lead compound discovery have been reported using above strategy, such as cancer16, 17, muscle atrophy18, acute myelogenous leukemia19, Parkinson's disease20, and Alzheimer's disease8.

The whole‐genome transcriptome profiles of hip articular cartilage specimens from four ONFH patients and four healthy controls will provide an excellent opportunity to explore the molecular mechanisms of this complex disease and to identify candidate drugs. Although the femoral head contains bone cells, osteoblasts, osteoclasts, and bone marrow mesenchymal stem cells, but the drug is generally in the different cell types6. So, the CMap database that was built upon four types of cultured human cell lines could be utilized to identify potential drugs for ONFH.

There were 1937 DEGs between ONFH samples and controls (fold‐change > 2.0, q value < 0.001), including 1223 up‐regulated genes and 714 down‐regulated genes. In order to better identify the core targets, we further filtered them based on the topological features. Among these genes, COL1A1, FGF2, PTGS2, RUNX1, and MMP13 have been reported to be involved in the regulation of cell apoptosis, cell differentiation, cell proliferation, fibrillar forming collagen, osteogenic differentiation, and angiogenesis. COL1a1, a member of group I collagen, can activate fibrillar forming collagen; FGF2 acts as an integrin ligand which is required for FGF2 signaling, and can induce angiogenesis; PTGS2 is responsible for the production of inflammatory prostaglandins which is a key step in the production of prostaglandin E2 (PGE2), and plays important roles in modulating motility, proliferation and resistance. RUNX1 forms the heterodimeric complex core‐binding factor (CBF) with CBFB and activates the expression of IL2 and IFNG and down‐regulates the expression of TNFRSF18, IL2RA, and CTLA4, in conventional T‐cells; and MMP13 plays a role in the degradation of extracellular matrix proteins including fibrillar collagen, fibronectin, TNC, and ACAN, and functions by degrading fibrillar collagens of types IV, XIV, and X. The PPI network analysis and association analysis indicated that some DEGs between ONFH samples and controls exhibit a correlation with the progression or prognosis of ONFH. After identifying key targets through the above screening, we performed further functional analysis and annotation using ClueGO. This analysis highlighted that the most enriched GO categories and pathways of the DEGs were blood vessel development, circulatory system development, blood vessel morphogenesis, angiogenesis, PI3K‐Akt signaling pathway, focal adhesion, and ECM‐receptor interaction (Fig. 6). The aforementioned GO categories and KEGG pathways were mainly related to the pathogenesis of ONFH.

Moreover, based on the 135 identified DEGs, the top 10 potential therapeutic molecules that had complete P value and specificity score for ONFH treatment were screened using negative connectivity scores as determined via the CMap database. Using the strategy could enhance the reliability of the potential therapeutic drugs, which provided a good foundation for the in vitro and in vivo studies. According to the literature, neostigmine bromide – a medication used to treat myasthenia gravis, Ogilvie syndrome, and urinary retention without the presence of a blockage21, 22 – shows that the connectivity score of −0.849 is of particular interest among the candidate molecules identified in our study. The molecular mechanisms and biological functions of neostigmine bromide included regulation of endothelial cell proliferation, positive regulation of stem cell proliferation, endodermal cell differentiation, and VEGF signaling pathway. VEGF may be involved in many physiological factors, including embryonic development and early postpartum development. Bone growth, endochondral osteogenesis, and inflammation were up‐regulated; at the same time, oxidative stress, growth factor, oncogene, and other factors could also regulate the expression of VEGF. The activation of VEGF in human ONFH have been demonstrated23, 24. The VEGF signaling pathway is directly related to cellular proliferation, migration, differentiation, and angiogenesis. Neostigmine bromide which have emerged as promising candidates in treatment of ONFH, were effective in delaying and even preventing disease progression.

Up to now, many scholars have verified the reliability of data mining results from GEO and CMap databases through experiments. Wang et al. validated eight novel molecular biomarkers associated with the diagnosis, prognosis prediction, and therapeutic targets of low‐grade glioma by mining GEO and CMap databases and qpcr25; Zhang et al. discovered two potential compounds for the treatment of gastric cancer by using weighted gene co‐expression network analysis, GEO, and CMap database, and verified their efficacy by in vivo and in vitro experiments26; Ren et al. also used GEO and CMap databases and in vivo and in vitro experiments to verify that baclofen may be a radioprotective agent27. But nonetheless, there are a few limitations and weaknesses in the present study. First, the reliance on the reliability of the GEO database and variable data quality, which can substantially affect the results. Due to objective factors, the small number of femoral head cartilage samples collected in this study and the failure to classify them according to etiology and classification is also one of the disadvantageous factors. Second, we compared hip articular cartilage from ONFH patients and healthy controls, which might differ from subchondral bone that was the main manifestation of ONFH. Whereas ONFH was a systemic bone metabolic disease, the expression pattern of different tissues might be prejudiced to their homologous function. In addition, experimental and clinical research have not been designed to confirm the above hypothesis about potential drugs, which is the shortcoming in the current study.

As a result, we utilized the gene expression profiles of femoral head from ONFH patients coupled with the CMap database and bioinformatics methods. This study reveals the dysfunctional signaling pathways of ONFH and identified the possible compounds for onfh's treatment. Although this is only a simple step towards success, the above results are still very useful for better understanding the pathogenesis of ONFH and providing important information for further animal and clinical trials to prove the efficacy of neostigmine.