Functional network analysis of gene-phenotype connectivity based on pioglitazone.

Pioglitazone, a type of insulin sensitizer, serves as an effective anti-hyperglycemic drug. The mechanism of action of pioglitazone is through the activation of the peroxisome proliferator-activated receptor (PPAR), which results in enhanced insulin sensitivity of peripheral tissues and the liver, causing a reduction in the production and output of liver sugar. It has been reported that pioglitazone increases the risk of bladder cancer, but the underlying mechanisms have remained elusive. It was hypothesized that modulation of pioglitazone activity may be predicted by systematically analyzing data published on drugs. This hypothesis was tested by querying the Drug-Target Interactome (DTome), a web-based tool that provides open-source data from three databases (DrugBank, PharmGSK and Protein Interaction Network analysis). A total of 4 direct target proteins (DTPs) and further DTP-associated genes were identified for pioglitazone. Subsequently, an enrichment analysis was performed for all DTP-associated genes using Cytoscape software. A total of 12 Kyoto Encyclopedia of Genes and Genomes pathways were identified, including the 'PPAR signaling pathway' as well as 'pathways in cancer' as relevant pathways. Functional network analysis was able to identify direct and indirect target genes of pioglitazone, resulting in a list of possible biological functions based on published databases. Furthermore, Kaplan-Meier analysis indicated that pioglitazone may affect the survival rate of patients with bladder cancer through genetic alterations (missense mutation, truncating mutation, amplification, deep deletion and fusion) of target genes. Therefore, it should be used with caution. Copyright: © Wang et al.

Introduction

Pioglitazone is a type of thiazolidinedione anti-diabetic drug and a highly selective peroxisome proliferator-activated receptor (PPAR) agonist. PPAR, a nuclear hormone receptor, combines with the retinoic acid X receptor (RXR) to form a heterodimer, which then binds to the peroxisome proliferator response element of the target gene to activate the PPAR signaling pathway (1). It is able to significantly improve glycemic control in type 2 diabetes mellitus and was approved by the US Food and Drug Administration in 1999 (2,3). A prospective clinical trial has indicated that pioglitazone affects the incidence and mortality of macrovascular disease in type 2 diabetes mellitus (4). It has been suggested that pioglitazone not only improves insulin resistance and blood glucose, but also displays anti-oxidant activity. It exerts a protective effect on the body by affecting the formation of arteriosclerosis and instability of plaques. Its mechanism of action may include: i) Inhibition of the function of monocytes and/or macrophages, as well as inhibition of the production of inflammatory cytokines. ii) Inhibition of the expression of endothelial-cell adhesion molecules and reduction of the interaction between leukocytes and endothelial cells. iii) Direct action on vascular smooth muscle cells, inhibition of the growth and calcium intake of vascular smooth muscle cells and dilation of blood vessels. iv) Promotion of reverse transportation of cholesterol, increase of plaque stability and delay of the process of arteriosclerosis. v) Reduction of apoptosis of cardiomyocytes and improvement of myocardial function (5,6). Other studies have indicated that pioglitazone improves insulin resistance and hyperandrogenism in patients with polycystic ovary syndrome (PCOS), and this effect may also enhance the responsiveness of PCOS patients to clomiphene (7). In recent years, the correlation between pioglitazone and tumor incidence has been increasingly investigated. Certain observational and retrospective studies have indicated that pioglitazone increases the risk of bladder cancer in a dose-dependent manner, compared with that in a control group (8,9). The present bioinformatics study set out to identify the direct or indirect targets of pioglitazone, as well as to aim to calculate or predict the most significant or comprehensive pathways or mechanisms of action using data from published databases.

Digital gene expression profiling (DGE) is a high-throughput sequencing technique used for detecting gene expression in specific tissues of a species. Through a bioinformatics search, analysis and comparison of gene expression profiles, specific information may be obtained, including that on gene regulation, gene transcription, signal transduction pathways and protein function (10). At present, the technology for gene expression profiling is mainly based on gene chip and high-throughput sequencing. Compared with traditional gene expression profiling chips, DGE technology has unique advantages: i) High-throughput-testing of almost all expressed genes in cells; ii) digitalization-the digital signals obtained by sequencing are direct and their resolution may reach 1 base difference between discriminable sequences; iii) good real-time performance-the sequence to be tested does not require to be designed beforehand; iv) ability to detect novel transcripts due to being independent on known genes. Therefore, it has been widely applied in the fields of life science and medicine (11–13). Mining these rich databases may become an effective novel method for studying the links and mechanisms between drugs and diseases (14,15). DrugBank is the only bioinformatics and chemical informatics database that combines detailed drug data with comprehensive drug target information (16,17). The DrugBank database contains information in numerous fields, including drug type, drug profile, chemical structure, drug ingredients, clinical trials, drug targets, enzymes, transporters, carriers, drug images, approval status, approved prescription drugs, foreign market names, drug interactions, manufacturers and packers. DrugBank has been widely used for computer-assisted retrieval of drug-associated data, drug recovery, retrieval of drug structural data, drug-target docking or screening, drug metabolism prediction and drug target prediction.

The purpose of the present study was to perform a functional network analysis of gene-phenotype connectivity based on pioglitazone. Direct and indirect target proteins of pioglitazone were identified. Enrichment analysis was then performed to identify the connective mechanism between pioglitazone and bladder cancer.

Materials and methods

Identification of drug target genes

DrugBank is the only bioinformatics and chemical informatics database that combines detailed drug data with comprehensive drug target information (17). Network interaction associations may be classified into four types: Drug-drug, drug-target, drug-gene and target-gene-protein associations. In the present study, drug-target associations were selected to build a pioglitazone-targeted network and the information collected is provided in Table I.

Table I.

Characterization of direct targets of pioglitazone using DrugBank.

Accession no.	Name	Target symbol	Uniprot ID	Uniprot name
DB01132	Pioglitazone	PPARG	P37231	Peroxisome proliferator-activated receptor gamma
DB01132	Pioglitazone	PPARD	Q03181	Peroxisome proliferator-activated receptor delta
DB01132	Pioglitazone	PPARA	Q07869	Peroxisome proliferator-activated receptor alpha
DB01132	Pioglitazone	MAOB	P27338	Monoamine oxidase B

Network visualization and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis

The Search Tool for the Retrieval of Interacting Genes and proteins (STRING; version 10.5) is a database that searches for interactions between proteins (18). It includes not only the direct physical interactions between proteins, but also the indirect functional correlations between proteins. In addition to experimental data, it also includes results mined from PubMed abstracts and data from other databases. It also uses bioinformatics methods to predict results, including interacting proteins. In the present study, STRING was used to construct the first- and second-level protein networks that are either directly or indirectly associated with pioglitazone. All protein data were imported into Cytoscape software (version 2.8.0) (19) in order to construct the protein network, while the Cluego plug-in was used for KEGG enrichment analysis in order to predict associated pathways and functions. The top 15 pathways with a significant P-value of <0.05 were selected for further analysis.

cBio cancer genomics portal cancer gene data visualization tool

The cBio Cancer Genomics Portal may aid in the analysis of molecular data obtained from cancer tissues and cytology, in order to recognize and understand heredity, epigenetics, gene expression and proteomics (20). The cBio portal provides graphical results that make complex cancer genomics data more understandable and acceptable, without the requirement of any special bioinformatics knowledge. Biological pathways, survival rates and numerous other pathways may also be retrieved through this portal.

The present study aimed to link pioglitazone-associated genes with all data on bladder cancer available in the cBio Portal in order to explore the association between target genes of pioglitazone and cancer. All pioglitazone-associated genes in all bladder cancer samples were explored and classified, while an integrative analysis of complex cancer genomics and clinical profiles was also performed (21).

Results

Pioglitazone-linked target gene search using Drugbank and STRING/visualization of the pioglitazone-target network using Cytoscape

The core of Cytoscape is the identification of networks. A simple network graph includes nodes and edges, with each node being a gene, microRNA or other entity. The edge between nodes represents the interaction between these nodes, including protein-protein interaction and DNA-protein interaction. In the present analysis, 4 directly associated target proteins [PPAR gamma (PPARG), PPAR delta (PPARD), PPAR alpha (PPARA) and monoamine oxidase B (MAOB) were identified using DrugBank and are listed in Table I. STRING software was then used to predict the interacting proteins of the 4 direct target proteins. The confidence value was set to 0.5 and the number of interacting proteins was set to 50. All proteins identified were either directly or indirectly associated with pioglitazone, and Cytoscape was used to create the relevant network presented in Fig. 1. All of these procedures were performed to obtain the pioglitazone network and visualize pioglitazone-associated target protein interactions.

Figure 1.

Drug-target interaction network of pioglitazone. Direct target proteins: PPARA, PPARG, PPARD and MAOB (indicated in yellow) and indirect target proteins (indicated in red). PPARA, peroxisome proliferator-activated receptor alpha; MAOB, monoamine oxidase B.

KEGG enrichment analysis of pioglitazone-associated genes using Cluego

KEGG is a database for the systematic analysis of gene function and genome information. It helps researchers to study genes and their expression information as a whole network. Pathway databases store information on gene functions and graphically represent biological processes within cells, including metabolism, membrane transport, signal transduction and cell growth cycle. Therefore, a KEGG pathway analysis of pioglitazone-associated genes was performed using Cluego. The top 15 pathways were selected, which included the PPAR signaling pathway (19 genes), Wnt signaling pathway (8 genes), thyroid hormone signaling pathway (15 genes), glucagon signaling pathway (8 genes), non-alcoholic fatty liver disease (11 genes), AGE/RAGE signaling pathway in diabetic complications (9 genes), fat digestion and absorption (5 genes), cholesterol metabolism (8 genes), cocaine addiction (5 genes), mitophagy (6 genes), JAK/STAT signaling pathway (11 genes), pathways in cancer (36 genes), transcriptional misregulation in cancer (13 genes) and adenosine monophosphate kinase signaling pathway (13 genes). The results are provided in Table II and Fig. 2. All of the pathways obtained through the enrichment analysis may be linked to the mechanism of action of pioglitazone. Of these pioglitazone-associated genes, certain candidates may also be involved in the pathogenesis of tumors, which warrants further research. PPARG is a member of the nuclear hormone receptor family and has been indicated to be expressed in numerous tumor types (22). At present, it is thought that it mainly induces cell differentiation, and that cell terminal transformation may stop cell growth, but it may also induce apoptosis of cancer cells (23,24). The Wnt signaling pathway is a conserved signaling pathway that has an important role in regulating normal embryonic development, cell proliferation and differentiation. Abnormal activation or imbalanced regulation of this gene may lead to tumor formation (25).

Table II.

List of KEGG pathways enriched by pioglitazone-associated genes determined using Cluego.

KEGG ID	Associated Pathway name	Genes genes(%)	number	Associated genes
KEGG:03320	PPAR signaling pathway	25.68	19	ACOX1, ADIPOQ, APOA1, APOA2, APOA5, APOC3, CD36, CPT1A, CYP7A1, EHHADH, FABP1, FABP4, LPL, PCK1, PLIN1, PPARA, PPARD, PPARG, RXRA
KEGG:04310	Wnt signaling pathway	5.48	8	CCND1, CREBBP, CSNK1E, EP300, JUN, NFATC1, PPARD, TP53
KEGG:04919	Thyroid hormone signaling pathway	12.93	15	CCND1, CREBBP, EP300, ESR1, FOXO1, HDAC3, MAPK1, MED1, MED24, NCOA1, NCOA2, NCOA3, NCOR1, RXRA, TP53
KEGG:04922	Glucagon signaling pathway	7.77	8	CPT1A, CREBBP, EP300, FOXO1, PCK1, PPARA, PPARGC1A, SIRT1
KEGG:04932	Non-alcoholic fatty liver disease	7.38	11	ADIPOQ, CEBPA, IKBKB, JUN, LEP, LEPR, NFKB1, PPARA, RELA, RXRA, TGFB1
KEGG:04933	AGE/RAGE signaling pathway in diabetic complications	9.09	9	CCND1, CDK4, FOXO1, JUN, MAPK1, NFATC1, NFKB1, RELA, TGFB1
KEGG:04975	Fat digestion and absorption	12.2	5	ABCA1, APOA1, APOB, CD36, FABP1
KEGG:04979	Cholesterol metabolism	16	8	ABCA1, APOA1, APOA2, APOB, APOC3, CD36, CYP7A1, LPL
KEGG:05030	Cocaine addiction	10.2	5	DDC, JUN, MAOB, NFKB1, RELA
KEGG:04137	Mitophagy	9.23	6	CITED2, E2F1, FOXO3, JUN, RELA, TP53
KEGG:04630	JAK/STAT signaling pathway	6.79	11	CCND1, CREBBP, EP300, IL2, IL2RA, IL2RB, IL2RG, IL4, IL4R, LEP, LEPR
KEGG:05200	Pathways in cancer	6.84	36	AR, CCND1, CDK2, CDK4, CEBPA, CHUK, CREBBP, E2F1, E2F2, E2F3, EP300, ESR1, FOXO1, HSP90AA1, HSP90AB1, IKBKB, IL2, IL2RA, IL2RB, IL2RG, IL4, IL4R, JUN, MAPK1, NCOA1, NCOA3, NFKB1, NFKBIA, PPARD, PPARG, RARA, RB1, RELA, RXRA, TGFB1, TP53
KEGG:05202	Transcriptional misregulation in cancer	6.99	13	CEBPA, CEBPB, FOXO1, IL2RB, NCOR1, NFKB1, PER2, PPARG, RARA, REL, RELA, RXRA, TP53
KEGG:04152	AMPK signaling pathway	10.83	13	ADIPOQ, CCND1, CD36, CPT1A, FOXO1, FOXO3, LEP, LEPR, PCK1, PPARG, PPARGC1A, SIRT1, SLC2A4
KEGG:05216	Thyroid cancer	13.51	5	CCND1, MAPK1, PPARG, RXRA, TP53

KEGG, Kyoto Encyclopedia of Genes and Genomes; PPAR, peroxisome proliferator-activated receptor.

Figure 2.

Kyoto Encyclopedia of Genes and Genomes enrichment analysis of pioglitazone-associated genes using Cluego. Red points represent pioglitazone-associated genes; Green points represent enrichment pathways and the larger the point, the more genes are enriched.

Mining of genes (PPARD, PPARG and RXRA) associated with pioglitazone-associated genes in bladder cancer using the cBio portal

KEGG enrichment analysis was performed on pioglitazone-associated genes and the relevant pathways obtained are presented in Table II. In order to further explore the association between pioglitazone-associated genes and tumor pathways, a web-based data mining database (cBio portal) was used to search for genes associated with bladder cancer. Since the PPAR signaling pathway is the major pathway among the pioglitazone-associated genes and the aim of the present study was to identify a common link with gene expression in bladder cancer, the genes in the PPAR signaling pathway and pathways in cancer were intersected, which led to the discovery of three overlapping genes (PPARD, PPARG and RXRA). These three overlapping genes were selected and queried for bladder cancer using the cBio portal. All of the 11 bladder cancer studies published prior to 2018 (26–35) were analyzed; however, certain studies are not listed on PubMed. Alterations, including gene mutation and amplification, ranged from 0–27%, and a total of 2,887 samples were queried (Fig. 3). For PPARG, most alterations were amplifications, but they also included deep deletions, fusions and missense mutations (Fig. 4). Alterations of the three selected genes in bladder cancer are provided in Fig. 4. Analysis of The Cancer Genome Atlas data published in 2014 (31) indicated that the most obvious alteration among all of the studies was in the one by Oncoprint (27). The query contained 1 gene pair (PPARD and PPARG) with mutually exclusive alterations (none significant) and 2 gene pairs with co-occurrent alterations (1 significant). Co-occurrence of PPARD and RXRA across samples was indicated to be significant (P=0.020), as presented in Table III.

Figure 3.

Summary of alterations for pioglitazone-associated target genes (PPARD, PPARG, RXRA) in bladder cancer by cBioportal. TCGA, The Cancer Genome Atlas; PPARD, peroxisome proliferator-activated receptor delta; RXRA, retinoid X receptor alpha; BLCA, bladder cancer; DFCI, dana-farber cancer institute; MSKCC, memorial sloan kettering cancer center; BGI, beijing genomics institution; NMIBC, nonmuscle invasive bladder cancer; PV, provisional.

Figure 4.

Characteristics of alterations of three pioglitazone-associated target genes (PPARD, PPARG, RXRA) in bladder cancer determined by oncoprint. PPARD, peroxisome proliferator-activated receptor delta; RXRA, retinoid X receptor alpha.

Table III.

Mutual exclusivity of gene sets.

Gene A	Gene B	Neither	A not B	B not A	Both	Log odds ratio	Adjusted P-value	Characteristic
PPARD	RXRA	2711	42	127	7	1.269	0.020	Co-occurrence
PPARD	PPARG	2597	48	241	1	−1.494	0.220	Mutual exclusivity
PPARG	RXRA	2527	226	118	16	0.416	0.273	Co-occurrence

Neither: Numbers of samples with alterations in neither gene A nor gene B; A Not B: Number of samples with alterations in gene A but not in gene B; B not A: Number of samples with alterations in gene B but not in gene A; Both: Number of samples with alterations in both gene A and gene B; Log odds ratio: Quantifies how strongly the presence or absence of alterations in gene A are associated with the presence or absence of alterations in gene B in the selected samples. PPARD, peroxisome proliferator-activated receptor delta; RXRA, retinoid X receptor alpha.

Creating a survival curve for bladder cancer using the cBio portal

Kaplan-Meier curves provide prognostic maps regarding overall or disease-free survival. They are plots of the percentage of surviving patients in different groups against the time, with stratification by the presence or absence of least one variant in the tumor sample, and the P-value regarding the difference between the curves may be determined to assess the significance of the feature in question as a prognostic factor. Using the cBio portal, the survival curves for mutation and non-mutation of three target genes in the bladder cancer samples were obtained, as provided in Fig. 5. Mutations in either of the three target genes (PPARD, PPARG and RXRA) resulted in significant differences in survival rates, compared with those of bladder cancer patients without the mutation (P<0.05). This suggests that pioglitazone may affect the survival rate of bladder cancer patients by targeting these three genes and associated pathways.

Figure 5.

Survival curves of mutation and non-mutation of three target genes in bladder cancer samples.

Discussion

At present, the association between pioglitazone and bladder cancer remains controversial. Meta-analysis also proves that high doses of pioglitazone exposure may be linked to a higher risk of developing bladder cancer (36,37). However, the relevant mechanisms have remained to be elucidated. Pioglitazone may produce a hypoglycemic effect by activating PPARG (2), and it also acts on bladder tissue (38). It was also reported that activation of PPARG may affect the proliferation of tumor cells (39). Bishop-Bailey and Warner (40) first proposed the role of the epidermal growth factor receptor signaling pathway in bladder epithelial cell differentiation. PPARγ has been proven to affect cell differentiation by regulating epithelial differentiation-associated antigens and keratin. However, this study did not determine when PPARγ has a role in the cell cycle. Sato et al (41) indicated that pioglitazone induces solid crystallization in urine of male rats, while an acidic diet inhibits pioglitazone-induced bladder tumors in male rats. Microcrystallization of advanced proliferative lesions, particularly calcium crystallization, may induce cytotoxicity and stimulate the proliferation of epithelial cells, which may be the mechanism of development of pioglitazone-induced bladder cancer (42).

New methods and technologies are required to explore intermolecular interactions and mechanisms in order to develop precise medical treatment methods. Bioinformatics is a discipline that studies the collection, processing, storage, dissemination, analysis and interpretation of biological information. In the present study, a functional network analysis of gene-phenotype connectivity based on pioglitazone was created with the aim to explore the molecular mechanisms of pioglitazone and its association with bladder cancer using web-based tools and software (DrugBank, Cytoscape and the cBio portal). First, searches with DrugBank and STRING were performed to identify direct and indirect target proteins of pioglitazone. The alterations of the target genes were then verified in bladder cancer samples using the cBio portal. Using this process, the pioglitazone-associated genes were linked with bladder cancer in particular common pathways. This is a novel type of bioinformatics application.

The 4 direct target proteins of pioglitazone identified are PPARA, PPARD, PPARG and MAOB. Genes associated with those target genes were also predicted in order to then perform pathway analysis. Analysis of all associated genes yielded 15 significant pathways. The PPAR signaling pathway and Wnt signaling pathway are the first two most significant pathways that were found to be enriched by the pioglitazone-associated genes. It has been reported that activation of PPARγ may affect the proliferation of tumor cells (39). Therefore, it is indicated that the PPAR signaling pathway is associated with cancer. The Wnt signaling pathway is a conserved signaling pathway that has an important role in regulating normal embryonic development, cell proliferation and differentiation. Abnormal activation or imbalanced regulation of genes may lead to tumor formation (25). It is an interesting phenomenon that most enriched pathways shared certain mechanisms with cancer. The PPARD, PPARG and RXRA genes were indicated to participate in the PPAR signaling pathway and also in pathways in cancer. PPARG and PPARG are involved in regulating the cell cycle, which may have an effect on bladder cancer. Mutations of the RXRA gene are also associated with bladder cancer, as a result of activating peroxisome proliferator-activated receptors, and this mechanism has been previously observed in urothelial cell proliferation (43). Mutation of RXRA has been reported to be accountable for 20–25% of bladder cancers. The results of the current study suggested that mutations of PPARD include amplification, missense mutations and deep deletions, with the major change being amplifications. Mutations of PPARG include amplifications, fusions, missense mutations and deep deletions, while amplifications were indicated to be the most significant change. Mutations in RXRA comprise deep deletions, missense mutations, amplifications and truncating mutations. This indicates that multidirectional mixed mutations may be involved in the mechanisms of the development of bladder cancer.

Biological processes or pathways in cancer are frequently regulated by a variety of genes or mechanisms. Mutual exclusivity may be discovered using the cBio portal in order to reveal previously unknown mechanisms of cancer, which may have an important role in tumorigenesis and cancer progression. According to this concept, genes associated with a particular tumor type tend to be mutually exclusive. Genetic exclusion is when a specific tumor only involves one specific gene. On the contrary, if there are multiple genes involved in tumors, these genes may coexist and have a role in the occurrence and development of tumors, while these tumors may not be the result of a single gene problem. The present results demonstrate that the PPARD and RXRA genes have a co-occurrence association, which indicates multiple gene alterations in bladder cancer. The results of the survival curve analysis reveal that alterations of the three genes selected (PPARD, PPARG and RXRA) lead to significant differences in overall survival (P<0.01). The use of pioglitazone may induce bladder cancer through causing variations in PPARD, PPARG and RXRA. These alterations result from amplifications, fusions, missense mutations or deep deletions.

Bioinformatics has become an important part and a frontier of biomedical research. Bioinformatics aims to clarify and interpret the biological significance from massive quantities of biomedical data, reveal and understand the complexity of genomic information structure and the fundamental laws of genetic language, and it features a combination of genomes, information structure and complexity. In fact, along with the rapid development of high-throughput sequencing technology, the large amount of biomedical data available has changed the focus of bioinformatics from data generation to data analysis. Methods of searching for and identifying trends using the data available in order to further interpret mechanisms associated with human health and diseases has become one of the major problems that require to be urgently solved at present. Functional analysis of drugs using relevant databases, including DrugBank and cBio portal, may help us create a more in-depth understanding of the mechanisms of action of drugs without experimental error that is included in traditional research (44,45). Disease-specific gene expression analysis may be used as a novel method of bioinformatics (46). Similarly, functional network analysis may serve as a strategy for further basic research in the traditional sense, as a method that may be more conducive for the discovery of mechanisms of disease development and drug action, and to also further develop biomedicine.

Pioglitazone, as a commonly used hypoglycemic agent, may act through PPARG, but at the same time, its association with side effects requires to be considered (47). A previous functional network analysis has indicated that pioglitazone may increase the risk of prostate and pancreatic cancer (48). The effect of PPARG agonists on bladder cancer remains controversial. Yang et al (49) indicated that pioglitazone did not promote malignant alterations of normal urothelial transitional epithelium cells or stimulate the proliferation of bladder cancer cells. Lv et al (50) demonstrated that PPARγ activation by pioglitazone markedly induced cell cycle arrest in G2 phase and apoptosis in bladder cancer cells, which resulted in suppression of tumor growth. Survival analysis was used in the present study to propose that pioglitazone could alter the survival rate of bladder cancer through associated the presence of target genes. Therefore, it should be used with caution.