System Level Meta-analysis of Microarray Datasets for Elucidation of Diabetes Mellitus Pathobiology.

BACKGROUND: Type 2 diabetes (T2D) is a common multi-factorial disease that is primarily ac-counted to ineffective insulin action in lowering blood glucose level and later escalates to impaired insu-lin secretion by pancreatic β cells. Deregulation in insulin signaling to its target organs is attributed to this disease phenotype. Various genome-wide microarray studies from multiple insulin responsive tis-sues have been conducted in past but due to inherent noise in microarray data and heterogeneity in dis-ease etiology; reproduction of prioritized pathways/genes is very low across various studies.
OBJECTIVE: In this study, we aim to identify consensus signaling and metabolic pathways through system level meta-analysis of multiple expression-sets to elucidate T2D pathobiology.
METHOD: We used 'R', an open source statistical environment, which is routinely used for Microarray data analysis particularly using special sets of packages available at Bioconductor. We primarily focused on gene-set analysis methods to elucidate various aspects of T2D. RESULT: Literature-based evidences have shown the success of our approach in exploring various known aspects of diabetes pathophysiology.
CONCLUSION: Our study stressed the need to develop novel bioinformatics workflows to advance our understanding further in insulin signaling.

INTRODUCTION

Diabetes has been considered a serious health problem; World Health Organization (WHO) has placed it among one of the four priority non-communicable diseases (NCD) and called for action on World Health Day 2016 [1]. There are an estimated 422 million adults with diabetes [1]. This complex metabolic condition depends on various factors such as high-calorie intake, sedentary lifestyle, obesity, and stress, as well as aging. Though global estimates of diabetes prevalence for type 1 and type 2 do not exist, the majority of people with diabetes are affected by type 2 diabetes which is generally attributed to insulin resistance in target organs and tissues; liver, skeletal, adipose and peripheral blood cells. T2D progresses through different disease phenotypes viz. Normal Glucose Tolerance (NGT) to Impaired Glucose Tolerance (IGT), Hyperinsulinemia (HI) and Insulin Resistance (IR) and eventually leads to overt diabetes. Different cultural settings clearly implicate the prevalence of T2D. There is definitely a pressing need to elucidate the mechanism underlying diabetes pathogenesis through the identification of disease-implicated signaling, and metabolic pathways at genome-wide level in order to design effective treatment.

Various genome-wide microarray studies pertaining to different disease phenotypes (IGT, IR, and T2D) have been conducted in past by various research groups and data is available at public repositories like Gene Expression Omnibus (GEO) [2] and Array Express [3] etc. However, due to inherent noise in microarray data as well as genetic differences in study population, disease heterogeneity, and non-uniform physiologic condition at the time of sampling and method of sampling, replication of most of the genes in different studies was very poor [4].

This situation clearly demands a need for robust bioinformatics meta-analysis of these microarray datasets, independent of geographical or other variables. In this study, we described a novel approach of the meta-analysis of microarray datasets to elucidate common signaling and metabolic pathways in an attempt to shed light on diabetes pathobiology.

MATERIALS AND METHODS

Microarray Data & Analysis Tools

We have downloaded twelve microarray datasets pertaining to various case-control diabetic studies from GEO. Those datasets containing samples from more than one disease phenotypes were further sub-divided and resulting twenty subsets were grouped tissue-wise: skeletal (11 sets), subcutaneous adipose (4 sets), peripheral (3 sets) and liver (2 sets). Selection of these datasets and their sub-division was based on a previous study of candidate gene identification in type 2 diabetes mellitus [5]. Some studies also included human subjects monitored for other diabetes-related clinical parameters i.e. disease family history, diet, and physical training regimen. Hence, these datasets can be considered an excellent collection of T2D gene-expression profiles due to the inclusion of disparate environmental conditions that seem vital for the elucidation of diabetes pathophysiology. Metadata of microarray datasets used in this study are shown in Table (. We used ‘R’ an open source Statistical software environment, used for Microarray data analyses particularly using a special set of packages available at bioconductor.

Table 1

Microarray samples used in meta-analysis.

S. No.	GEO Series	Tissue	Affymetrix chip Type		Number of Samples
				Place of Study	Normal Glucose Tolerant (NGT)/Healthy	Impaired Glucose Tolerant (IGT)/ Prediabetics/ Insulin Resistant (IR)	T2DM
1	GSE18732 (G1)	Skeletal	Affymetrix Gene Chip Human Genome U133 Plus 2.0 Array	Tissue Injury & Repair, Edinburgh University, Edinburgh, UK	20		20
2	GSE18732 (G2)					20	20
3	GSE19420(G1)		Affymetrix Human Genome U133 Plus 2.0 Array	Genetics and Cell biology, Maastricht University, Maastricht, Netherlands		12	10
4	GSE19420(G2)					12	8 (after 52 weeks Training)
5	GSE19420(G3)				12		10
6	GSE19420(G4)				12		8 (after 52 weeks Training)
7	GSE25462(G1)		Affymetrix Human Genome U133 Plus 2.0 Array	Research Division, Joslin Diabetes Center, Boston, MA, USA	15 (Family History-)		10
8	GSE25462(G2)					25 (Family History +)	10
9	GSE12643		Affymetrix Human Genome U95 Version 2 Array	Odense University Hospital, Odense, Denmark		10	10
10	GSE22309(G1)		Affymetrix Human Genome U95A Array	Biostatistics, University of Alabama at Birmingham, AL, USA	20		15
11	GSE22309(G2)					20	15
12	GSE26637(G1)	Adipose tissue	Affymetrix Human Genome U133 Plus 2.0 Array	Institute for Molecular Medicine Finland (FIMM), University of Helsinki, Helsinki, Finland	5 (Fasting)	5 (Fasting)
13	GSE26637(G2)				5 (Hyperinsulinemia)	5 (Hyperinsulinemia)
14	GSE15773 (G1)		Affymetrix Human Genome U133 Plus 2.0 Array	Department of Molecular Medicine, University of Massachusetts, MA, USA	5	4
15	GSE15773 (G2)				5	5
16	GSE13760	Peripheral Blood	Affymetrix Human Genome U133A 2.0 Array	Department of Hematology, Roskilde Hospital, Roskilde, Denmark	11
17	GSE9006		Affymetrix Human Genome U133A Array	UTSW Medical Center, Dallas, TX, USA	24		12
18	GSE15932		Affymetrix Human Genome U133 Plus 2.0 Array	Department of surgey, Second Affiliated Hospital, School of Medicine, Zhejiang University, Hangzhou, China	8
19	GSE15653	Liver	Affymetrix Human Genome U133A Array	Department of Research Division, Joslin Diabetes Center, Boston, MA, USA	5		5
20	GSE23343		Affymetrix Human Genome U133 Plus 2.0 Array	Kanazawa University, Kanazawa, Japan	7		10

R Analysis Workflow for Meta-analysis of Microarray Datasets

After acquisition of microarray datasets from GEO, we first downloaded Custom Chip Description Format (CDF) files from Brainarray [6] (www.brainarray.mbni.med.umich.edu) as due to rapid updates in genome and transcriptome databases, microarray probesets generally suffer from outdated/incorrect annotations that might affect downstream functional analysis for biological interpretation. Separate AffyBatch objects for each dataset were created. These were subsequently background corrected, normalized and summarized into ExpressionSet class objects using GeneChip Robust Multiarray Analysis (GC-RMA) method from package gcrma [7]. This method is capable of achieving precise gene expression values due to the use of sequence specific probe affinities.

Various algorithms for microarray statistical meta-analysis have been described in recent years [8]. In this study, the ‘direct merging of expression values’ approach was adapted as all studies used Affymetrix gene chip platform. Normalized values obtained in different studies were expected to have statistical coherence and appropriate for creating tissue-specific consolidated Expression Set objects.

In order to improve the analysis of differentially expressed genes (DEGs), a nonspecific filtering step was performed using package genefilter [9] that first estimates the overall variability across the array of each probe set and then uses it to estimate the shortest interval with half of the data. This value was chosen as variance cutoff to create the reduced ExpressionSet objects. As our Meta-analysis employed different Affymatrix chip types, to create tissue-specific consolidated Expression set, the probes must be converted to their official gene symbols and/or Entrezgene ids as later has been reported to be more effective in the computational analysis because gene symbols tend to inadvertently mutate into ‘date format’ during excel sheet operations [10]. We used package org.Hs.eg.db [11] to obtain annotations of probes based on Entrez Gene identifiers. Finally, five tissue-specific consolidated sets were obtained: Adipose, Liver, Peripheral, Skeletal-IGT, and Skeletal-NGT. Last two skeletal sets differ with each other based on T2D phenotype contrasted with either Impaired Glucose Tolerant or Normal Glucose Tolerant phenotypes. For Differential Gene Expression (DGE) analysis, an empirical Bayes moderation approach based on modified t-statistics was used that was implemented in package limma [12].

Downstream analysis of the obtained gene list for biological interpretation is a more challenging task than engineering a data analysis pipeline, particularly for a complex disease like type 2 diabetes. Genes in the gene list, obtained through DGE analysis can be grouped either on the basis of their numerical expression values using statistical clustering or classification. However these approaches have been criticized, due to the fact that weakly expressed genes may also contribute to disease phenotype but generally failed to qualify statistical thresholds [13]. Alternatively genes can also be grouped according to shared biological features i.e. involvement in metabolic/signaling pathway or Gene Ontology (GO) Consortium [14]. Therefore, gene-set based methods are biologically more meaningful than single gene methods due to the fact that even small expression changes in the members of a pathway dramatically alter the flux through the pathway. Various methods have been in use in microarray community for functional categorization of genes in the form of gene-sets that are statistically associated with expression data. We used two different types of gene set enrichment methods; over-representation analysis based on Hypergeometric or Fisher test; these methods create gene-sets using hard cut-off in the form of p-value or expression fold-change or a combination of both using package GOstats [15]. Another type of method for gene-set analysis is based on regression test to find gene-sets statistically associated with studied phenotypes; package globaltest [16] provided necessary software infrastructure to carry out this analysis. In addition to canonical Kyoto Encyclopedia of Genes and Genomes (KEGG) [17] and Gene Ontology (GO) [18] categories, globaltest also facilitates gene set analysis against Molecular Signature Database (MSigDB v5.1) maintained by Broad Institute [19]. MSigDB is an inventory of 13311 annotated gene sets divided into eight major collections. Though Broad Institute also offered a method for gene-set analysis termed, “Gene Set Enrichment Analysis (GSEA)”, [19] that is based on Kolmogorov-Smirnov-like statistics but regression model based Global test implementation has been reported to have higher power than GSEA [20]. We downloaded complete Molecular Signatures Database (MSigDB v5.1) in XML format and package GSEABase [21] was used to obtain complete gene set collection in R environment. This collection was sub-setted into the Hallmark gene set (H), Curetted gene set (C2) and Motif gene set (C3). Further splitting of C2 and C3 yielded reactome gene set collection and microRNA target collection, respectively.

RESULTS

Insulin Receptor (IR) being an Receptor Tyrosine Kinase (RTK) may bind to various downstream signaling molecules with SH2 domain (Src Homology 2) viz. members of Insulin Receptor Substrate family (IRS-1, IRS-2), Growth factor receptor-bound protein 2 (Grb2), Phospholipase C γ (PLC-γ), and Phosphatidylinositol-4, 5-bisphosphate 3-kinase (PI3-K).

Our microarray meta-analysis of various insulin-responsive tissues has enriched disparate disease-related

signaling and metabolic pathways and gene–sets. Tables , and enlisted top 10 pathways/gene-sets for each functional analysis. Furthermore, literature based exploration pointed their connection with T2D through SH2 domain-containing downstream molecules. Fig. () further illustrates the signaling cascades downstream of Insulin Receptor initiated by PLC-γ, Grb2, and PI3-K.

Results from GOstats Pathways Analysis

GOstats based KEGG analysis (P < 0.05) enriched ‘Focal adhesions’ among top-ranking pathways in all the tissue sets. The high glucose level in extra cellular fluid exerts a detrimental effect on Focal Adhesion Kinase (FAK) mediated wound healing process [22] which has been reported to associate with lower limb amputation (LLA) [23]; the prevalence of later was found to be eight times higher in diabetic individuals than in nondiabetic individuals [24]. Interestingly FAK can also bind with Grb2-SoS complex, PI 3-kinase, and phospholipase C-γ through its phosphotyrosine residues like that of insulin receptor and this might be a reason for clinical co-occurence of T2D and LLA. Another top enriched pathway was ‘Pathways in cancer’ and the link between metabolic syndromes and cancer has previously been reported [25]. Insulin-signaling induces two major signaling pathways: Mitogenic MAPK pathway and Anti-apoptotic PI3-K pathway. PI3K activates Protein Kinase B (PKB) or Akt resulting in the activation of mTOR/Raptor complex which mediate its effects on mitogenesis and cell growth. Impaired insulin- signaling may, therefore, leads to deregulation in mTOR-signaling which has been linked to numerous human cancers [26]. The analysis results in this study can be found in Online Supporting Information S1.

Results from Globaltest Pathways Analysis

Global test based KEGG analysis primarily enrich various immune system related pathways in all tissue types i.e., ‘Graft- versus- Host disease’, ‘Allograft reaction’, ‘Asthma’, ‘Intestinal Immune network for IgA’, ‘primary Immunodeficiency’ and ‘Type I Diabetes’. The analysis in this study, clearly established the role of immune deregulation in Type 2 diabetes as reported earliar by Jaffe & group in 2004 [27]. In addition, various metabolic pathways were also enriched including ‘Cysteine and methionine’, ‘D-Glutamine and D-glutamate’, ‘Nicotinate and nicotinamide’, ‘Phenylalanine’, ‘tyrosine and tryptophan’ and ‘Ubiquinone and other terpenoid- quinone metabolism’pathway. The implication of these metabolic pathways in T2D pathophysiology seems more challenging than the exposition of the role of signaling pathways in disease etiology and requires further literature support. The results of the analysis in this study can be found in Online Supporting Information S2.

Results from MSigDB’s Gene-sets Analysis

The MSigDB gene-sets analysis also contributed significantly in exposing various aspects of disease etiology. Reactome analysis revealed ‘T Cell receptor (TCR) signaling’ among top five pathways in all the studies on tissue. TCR-signaling involves the formation of a multi-molecular structure ‘Signalosomes’ that includes Grb2 and PLC- γ [28]; it is interesting to note that Grb2 is a positive regulator of insulin-signaling pathway. It binds to phosphorylated Tyr of IRS-1 and propagates signal downstream of insulin receptor. While another component of signalosome, PLC- γ is a negative regulator of insulin- signaling, PI3K phosphorylates Phosphatidylinositol 4, 5-bisphosphate (PIP2) into Phosphatidylinositol (3,4,5)-trisphosphate (PIP3), PIP3 intern recruits Akt to the plasma membrane and Akt facilitates the movement of glucose transporter GLUT4 from the internal membrane vesicle to the plasma membrane; hence, stimulating the uptake of glucose [29, 30]. PLC- γ hydrolyzes PIP2 into diacylglycerol (DAG) and Inositol triphosphate (IP3) and may probably lead to insulin resistance. Presumably, T cell Receptor signaling and impaired insulin signaling cross-talk with each other through these shared components and give rise to the diabetic or IR phenotype. Deregulated T cell receptor signaling may also lead to inappropriate Antibody Dependent Cell Cytotoxicity (ADCC) [31] which may result in various immune complications like those enriched previously in the KEGG analysis. Top enriched pathway reported in the liver study was ‘Glutamate synaptic neurotransmission’; this pathway includes vesicle-associated membrane protein (VAMP) which is also involved in insulin secretion from pancreatic β cells as well as insulin signaling to responsive tissues by SNARE-Mediated Vesicle Exocytosis [32]. Top enriched pathway in skeletal muscles was ‘Pyruvate dehydrogenase (PDH) complex-mediated glucose metabolism’ and ‘Hormone-sensitive lipase (HSL)-mediated triacylglycerol hydrolysis’; HSL has also been reported to be associated with abdominal obesity in patients with type 2 diabetes [33].

Top ranking Hallmark pathway included ‘IL6 signaling via JAK/STAT pathway’ in all the studies on tissue. IL6 is an inflammatory cytokine which is released from the expanded visceral adipose tissue and adversely affects the insulin signaling cascade. The key element in the JAK/STAT pathway is the Signal Transducer and Activator of Transcription 3 (STAT3) protein that acts as a transcription factor and regulates the Suppressor of Cytokine Signaling (SOCS3) mRNA synthesis. Interestingly SOCS3 is a negative regulator of both cytokine signaling as well as insulin signaling. It binds to phosphorylated tyrosine residues on IR/IRS2 and enhances their proteolytic degradation through ubiquitin-proteasomal pathway. STAT3 signaling has also been reported to contribute to skeletal muscle insulin resistance in T2D [34]. Another pathway that was observed as one of the top pathways was ‘TGF-β signaling pathway’. TGF-β signaling is involved in the regulation of insulin gene transcription, pancreatic β cell function and glucose tolerance and energy homeostasis [35-38]. Invitro studies have confirmed its role in the inhibition of insulin-stimulated tyrosine phosphorylation of IRS-1 [39]. Moreover, it also inhibits insulin/insulin-like growth factor-1 (IGF-1)-dependent adipose cell differentiation (adipogenesis) process. TGF-β exerts this effect through the induction of extracellular matrix (ECM) synthesis, which in turn inhibits adipogenesis and leads towards metabolic syndrome.

MicroRNA target analysis revealed miR-184 in all four studies on tissue as the top 10 microRNA targets. miR-184 expression has been reported to be negatively correlated with insulin secretion in the pancreas; [40] it also inhibits Protein Kinase C- Beta (PRKCB) expression, a positive regulator of insulin signaling. PRKCB functions downstream of IRS1 in the muscle cells and mediates insulin-dependent DNA synthesis through the RAF1-MAPK/ERK signaling cascade. miR-184 was also reported to participate in the regulation of glucose transport in adipocytes by negatively modulating the insulin-stimulated translocation of the glucose transporter SLC2A4/GLUT4 [40]. The aforementioned facts implicate its role in the pathogenesis of diabetes. IGT and NGT studies on skeletal tissue shared miR-346 among the top 10 targets and an miRNA microarray study on Diabetic nephropathy mice confirmed its role in negative regulation of SMAD 3/4 mediated signaling which has been reported in the pathogenesis of obesity and T2D [41]. Together, these findings proved miR-346 a novel therapeutic target against diabetes. The results of analysis can be found in Online Supporting Information S3.

CONCLUSION

Insulin-signaling and its deregulation in the form of diabetes phenotype is quite an intractable problem and our present understanding of insulin signaling is based on recent advances in biochemistry and medicine [42]. Though insulin is a peptide hormone, insulin-signaling also shares features with growth factor signaling, it is also mediated by Receptor tyrosine kinases (RTKs) like those of EGF, TGFα, PDGF, FGF, VEGFR and CSF-1 signaling. It is speculated that the main objective of insulin is to manage under-utilization of glucose in order to support anabolic reactions like DNA replication, protein synthesis, mitogenesis, cell growth and differentiation [42].

Experimental evidence suggests that insulin signaling actually cross-talks with the signaling pathways of growth factors. The present functional analysis revealed that this cross-talk might be mediated through SH-2 domain containing signaling molecules participating in different signaling pathways. It might also be a probable reason of clinical interlink of T2D with cancer, immune-system, and focal adhesion related pathways. In addition, the present study further highlights the need to develop bioinformatics tools and workflows to advance our understanding in insulin signaling in an effort to combat diabetes.

Table 2

Top 10 common KEGG pathways enriched in GOstat analysis.

S. No.	KEGG Pathway id	Pathway Name
1	5200	Pathways in cancer
2	4510	Focal adhesion
3	4012	ErbB signaling pathway
4	5220	Chronic myeloid leukemia
5	4722	Neurotrophin signaling
6	5100	Bacterial invasion of epithelial cells
7	5212	Pancreatic cancer
8	4520	Adherens junction
9	5211	Renal cell carcinoma
10	5223	Non-small cell lung cancer

Table 3

Top 10 common KEGG pathways enriched in global-test analysis.

S. No.	KEGG id	Pathway Name
1	4940	Type I diabetes mellitus
2	760	Nicotinate and nicotinamide metabolism
3	5330	Allograft rejection
4	5332	Graft-versus-host disease
5	5340	Primary immunodeficiency
6	5310	Asthma
7	4672	Intestinal immune network for IgA production
8	400	Phenylalanine, tyrosine and tryptophan biosynthesis
9	270	Cysteine and methionine metabolism
10	471	D-Glutamine & D-glutamate metabolism

Table 4

Top 10 MSigDB gene-sets enriched in global-test analysis.

S. No	Gene-set	Pathway/ Mi RNA
1	Reactome	Phosphorylation of CD3 and TCR ZETA Chains
2		Signaling by TGF beta Receptor Complex
3		Regulation of Pyruvate DehydrogenasePDG complex
4		Hormone Sensitive Lipase HSL Mediated Triacylglycerol Hydrolysis
5	Hallmark	IL6-JAK STAT3 Signaling
6		TGF BETA Signaling
7	MiRNA	TCCGTCC,MIR-184
8		CGGTGTG,MIR-220
9		ATACTGT,MIR-144
10		GGCAGAC,MIR-346