G Model BONSOI-3981; No. of Pages 6
ARTICLE IN PRESS Joint Bone Spine xxx (2014) xxx–xxx
Available online at www.sciencedirect.com
Original article
Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets Tae-Hwan Kim a , Sung Jae Choi b , Young Ho Lee b , Gwan Gyu Song b , Jong Dae Ji b,∗ a b
Hanyang University Hospital for Rheumatic Diseases, Seoul, South Korea Division of Rheumatology, College of Medicine, Korea University, 126-1, Anam-Dong 5-Ga, Sungbuk-Ku, Seoul 136-705, South Korea
a r t i c l e
i n f o
Article history: Accepted 15 January 2014 Available online xxx Keywords: Rheumatoid arthritis Anti-TNF Response Microarray
a b s t r a c t Objectives: Anti-tumor necrosis factor (TNF) therapy is the treatment of choice for rheumatoid arthritis (RA) patients in whom standard disease-modifying anti-rheumatic drugs are ineffective. However, a substantial proportion of RA patients treated with anti-TNF agents do not show a significant clinical response. Therefore, biomarkers predicting response to anti-TNF agents are needed. Recently, gene expression profiling has been applied in research for developing such biomarkers. Methods: We compared gene expression profiles reported by previous studies dealing with the responsiveness of anti-TNF therapy in RA patients and attempted to identify differentially expressed genes (DEGs) that discriminated between responders and non-responders to anti-TNF therapy. We used microarray datasets available at the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO). Results: This analysis included 6 studies and 5 sets of microarray data that used peripheral blood samples for identification of DEGs predicting response to anti-TNF therapy. We found little overlap in the DEGs that were highly ranked in each study. Three DEGs including IL2RB, SH2D2A and G0S2 appeared in more than 1 study. In addition, a meta-analysis designed to increase statistical power found one DEG, G0S2 by the Fisher’s method. Conclusion: Our finding suggests the possibility that G0S2 plays as a biomarker to predict response to anti-TNF therapy in patients with rheumatoid arthritis. Further investigations based on larger studies are therefore needed to confirm the significance of G0S2 in predicting response to anti-TNF therapy. © 2014 Société franc¸aise de rhumatologie. Published by Elsevier Masson SAS. All rights reserved.
1. Introduction Rheumatoid arthritis (RA) is a systemic autoimmune disease characterized by chronic inflammation of joints followed by destruction of cartilage and bone. The precise etiology of RA is not completely understood, but a large body of studies indicates that pro-inflammatory cytokines such as TNF␣ and IL-1 are significantly involved in the pathophysiology of joint inflammation and bone destruction in RA [1]. Biologic agents, especially anti-TNF agents, have produced significant and clinically relevant improvements in patients with active RA and appear to be very effective therapeutic molecules for RA patients. Although the use of anti-TNF agents has revolutionized the treatment of RA, a substantial proportion (up to one-third) of RA patients treated with anti-TNF agents fail to achieve a satisfactory clinical response [2]. Therefore, a strategy for predicting which patients will and will not respond prior to initiation of therapy is needed to prevent disease progression,
∗ Corresponding author. Tel.: +82 2 920 5489; fax: +82 2 920 5974. E-mail address: [email protected] (J.D. Ji).
unnecessary side effects, and high costs for non-responders. Although baseline demographics and clinical characteristics, including age, sex, disease duration, Disease Activity Score (DAS28), Health Assessment Questionnaire score, and concurrent use of disease-modifying anti-rheumatic drugs (DMARDs), have been extensively studied, potential prediction markers of therapeutic responsiveness have not yet been discovered [3]. Pharmacogenetic studies have also been performed to discover determinants of efficacy of anti-TNF agents. Some studies have shown that polymorphisms in genes encoding TNF␣, TNF receptors, the major histocompatibility complex region, and the Fc␥ receptor IIIA (FC␥RIIIA) are associated with a favorable response to anti-TNF therapy; however, these results were not consistent with those of other studies [3,4]. Gene expression profiling with microarrays has become a standard method for identifying the genes and biological pathways that are associated with various complex diseases [5,6]. Biomarkers identified by gene expression profiling can be used to make a diagnosis, monitor disease activity, and predict response to therapy. This technology has already been shown to identify gene expression profiles predictive of treatment response and clinical outcomes
1297-319X/$ – see front matter © 2014 Société franc¸aise de rhumatologie. Published by Elsevier Masson SAS. All rights reserved. doi:10.1016/j.jbspin.2014.01.013
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
G Model BONSOI-3981; No. of Pages 6
ARTICLE IN PRESS T.-H. Kim et al. / Joint Bone Spine xxx (2014) xxx–xxx
2
expression profiles predicting response to anti-TNF therapy were included in the analysis. 2.2. Individual analysis
Fig. 1. Flow chart of selection process for meta-analysis. A. Available studies were identified via Medline/PubMed searches. B. Available GEO datasets were identified through GEO datasets searching.
for some diseases [7–11]. Recently, several studies have attempted to identify gene expression profiles predicting the response to anti-TNF therapy in RA patients by using multiple gene expression microarrays [12–17]. Although these studies have successfully identified gene expression signatures predicting the response to anti-TNF therapy, the expression profiles identified in these studies were not consistent with each other, and the genes reported by the authors of each study had little overlap. Herein, we attempt to identify differentially expressed genes (DEGs) that differ between responders and non-responders to anti-TNF therapy by: • comparing gene signatures previously reported to predict response to anti-TNF therapy in RA patients; • analyzing microarray datasets of these studies, which are available at the National Center for Biotechnology Information’s Gene Expression Omnibus (GEO), using GEO2R; • performing a meta-analysis of these datasets. 2. Methods 2.1. Data collection We performed literature survey via Medline/PubMed searches to identify available studies that examined gene expression profiles of peripheral blood samples to predict treatment outcome of anti-TNF therapy in RA patients. We searched for a combination of KEY words such as “Infliximab”, “Adalimumab”, “Etanercept”, “Microarray”, “Gene expression profile”, “Anti-TNF” and “Rheumatoid arthritis” as both medical subject heading terms and text words. We also tried to identify additional studies by manually searching references of the original articles or review articles on this topic. Twenty-two relevant studies that investigated gene expression profiling in RA patients treated with anti-TNF agents were identified using PubMed. On the basis of their titles and abstracts, 9 studies were selected for full-text review (Fig. 1A). Three studies was excluded from analysis because they used synovial tissues for examining gene expression signatures [18,19] or did not focus on a predictor of the clinical response to anti-TNF therapy [20]. The remaining 6 studies were included in our analysis [12–17]. In addition, 10 sets of microarray data were identified by searching GEO DataSets. Two datasets, GSE15602 and GSE21537, were excluded from analysis because they used synovial tissues, and 3 datasets, GSE7524, GSE12897, and GSE19821, were excluded because they did not focus on a predictor of clinical response to anti-TNF therapy (Fig. 1B). Among 6 including studies, the dataset of one study was not deposited in the GEO database [15]. Finally, 6 studies and 5 sets of microarray data (GSE3592, GSE8350, GSE12051, GSE20690, and GSE33377) that used peripheral blood samples to identify gene
We analyzed each dataset individually by using the same analytical tool, GEO2R (GEO’s online tool for analyzing GEO data, available at http://www.nci.nlm.nih.gov/geo/geo2r/), to maintain consistency during individual analysis. GEO2R is an interactive, online tool for R-based analysis of GEO data to identify genes that are differentially expressed across experimental conditions [21]. GEO2R analyzes original, submitter-supplied data with the GEOquery and limma (Linear Models for Microarray Analysis) R packages from the Bioconductor project. The basic statistic used in limma for significance analysis is the moderated t-statistic. The Benjamini & Hochberg false discovery rate method is used to apply P-value adjustment for multiple-testing correction. The P-value and adjusted P-value were used to identify a list of DEGs from each dataset. Each list of DEGs was used for pathway enrichment analysis. This analysis was done using the Database for Annotation, Visualization and Integrated Discovery (DAVID) (http://david.abcc.ncifcrf.gov/home.jsp) [22]. Several annotation categories, including gene ontology (GO) terms, BIOCARTA pathways, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, were subjected to enrichment analysis. 2.3. Meta-analysis To combine the results of individual studies and obtain a list of more robust DEGs between responders and non-responders to antiTNF therapy, we followed the guideline provided by Ramasamy et al. for meta-analysis of gene expression microarray datasets [23]. We used 2 R packages, MetaQC and MetaDE, for quality control (QC) and identification of DEGs, respectively [24]. MetaQC implements 6 quantitative quality control measures, including internal quality control (IQC, internal homogeneity of co-expression structure among studies), external quality control (EQC, external consistency of co-expression structure correlating with pathway database), accuracy quality control (AQCg and AQCp, accuracy of DEG detection or pathway identification), and consistency quality control (CQCg and CQCp, consistency of differential expression ranking for genes or pathways) [25]. In addition, the mean rank of all QC measures in each dataset was computed as a quantitative summary score by calculating the ranks of each QC measure among all included datasets. The moderated t-statistic was used to calculate P-values in each dataset, and a meta-analysis was performed with the MetaDE package to identify DEGs by the Fisher method and the maximum P-value method. 3. Results 3.1. Study characteristics The analyses included 6 studies that used peripheral blood samples to identify gene expression profiles predicting response to anti-TNF therapy. The dataset of 1 of these studies was not deposited in the GEO database [15]. Therefore, the microarray dataset of this study was not included for individual analysis or meta-analysis. Detailed information of each study is described in Table 1. In the study of Lequerré et al., patients with very active disease (Disease Activity Score 28 > 5.1) were included in this study and mononuclear cell RNAs were collected at baseline. In the study of Sekiguchi et al., patients having inadequate control of RA symptoms were included and whole blood RNA was extracted before the infusion of infliximab. In the study of Julia et al., patients with active disease (Disease Activity Score > 3.2) were included and
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
G Model
ARTICLE IN PRESS
BONSOI-3981; No. of Pages 6
T.-H. Kim et al. / Joint Bone Spine xxx (2014) xxx–xxx
3
Table 1 Summary information of studies included in analysis. Study
Toonen et al. [17] Stuhlmüller et al. [15] Tanino et al. [16] Julia et al. [12] Sekiguchi et al. [14] Lequerré et al. [13]
Year
2012 2010 2009 2009 2008 2006
Platform
HuEx-1 0-st HG-U133A Agilent-014850 Illumina Human-6 v1 Customized Customized
Number of samples
Responder
Non-responder
18 5 43 37 8 6
24 2 25 7 10 7
Response determination
Response criteria
GEO dataset
Tissue
Anti-TNF
14 weeks – 14 weeks 14 weeks 22 weeks 3 months
EULAR ACR CRP EULAR ACR DAS28
GSE33377 – GSE20690 GSE12051 GSE8350 GSE3592
Whole blood Monocyte Whole blood Whole blood Whole blood PBMC
I, A A I I I I
I: infliximab; A: adalimumab.
whole blood samples were extracted before treatment. In the study of Tanino et al., patients were resistant to standard methotrexate treatment and whole blood samples were collected before treatment. In the study of Toonen et al., patients with active disease (Disease Activity Score > 3.2) and previous failure on at least two disease-modifying anti-rheumatic drugs (DMARDs) were included and whole blood was sampled before treatment. Gene expression signatures predictive of anti-TNF treatment outcome were reported in 5 studies, and Toonen et al. tried to validate previously reported gene expression signatures using microarray data on whole blood from RA patients treated with anti-TNF agents. We observed several clinical and experimental variations between these studies, including differences in microarray platform, types of anti-TNF therapy,
response criteria, and tissues used for analysis (whole blood versus monocytes). The gene expression profiles used in different studies were not consistent with each other, and only 2 genes, IL2RB and AP1S2, were used in more than 1 study (Table 2). 3.2. Individual analysis The included datasets were individually analyzed with GEO2R. When the Benjamini & Hochberg false discovery rate method was used to adjust statistics for multiple-testing (adjusted P-value < 0.05 was the cutoff for statistical significance), DEGs for different response groups were not identified in any of the datasets. Because no genes passed the cutoff for statistical significance after
Table 2 Transcript sets of studies included in analysis. Lequerré et al. (20 genes)
Sekiguchi et al. (18 genes)
Julia et al. (8 genes)
Tanino et al. (8 genes)
Stuhlmüller et al. (82 genes)a
AKAP9 COX7A2L ELMOD2 EPS15 FBXO5 HLA-DPB1 LAMR1 MCP MRPL22 MTCBP-1 PFKFB4 PSMB9 PTPN12 QIL1 RASGRP3 RPL35 RPS16 RPS28 SCAM-1 TBL2
AP1S2 CX3CR1 EIF2AK2 HLA-DQA1 IFI6 IFIT1 IFIT3 IFITM1 IFITM2 IGHM IL2RB MX1 MX2 OAS1 OAS2 OASL SRP9 ZFP36L2
CAMP GNLY HLA-DRB3 IL2RB MXD4 SH2D1B SLC2A3 TLR5
ANKRD55 ATP5I C21orf58 CLGN LOC643981 PSPH TBC1D8 TMEM141
ABCC3 ACTG1 ACTN1 AP1S2 APP AQP3 ATXN1 C14orf156 C20orf111 CALM1 CCDC91 CCL3L3 CD86 CD97 CDC2L6 CLC CLMN CNOT6 COX7B CPA3 CRYZ CTSZ CUGBP2 CYB5A DOCK5 DSC2 EDARADD EGR1 ENPP4 ERMAP ESD F2RL1 FAM3C FCER1A GATA2 GCLC GNA12 GNB1 GNE GNS
a
HLA-DQA2 HLA-DQB2 HMBOX1 IDH1 IL8 IL8RB ITGAX KIAA0746 KLF11 LBH LEPR MAL NADK NARF NDUFB1 NGRN NLRP3 PAIP1 PECAM1 PF4V1 PIM1 PLCL2 PTMAP7 RAPGEF2 RPS24 RPS26 SC4MOL SF3B1 SOD1 SUB1 TCF4 TM2D1 TMED2 TMEM176A TMEM176B TMEM45A TPD52 UCRC UQCRH WIPF1
Two Affymetrix ID (214807 at and AFFX-M27830 5 at) gene symbols not available.
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
G Model
ARTICLE IN PRESS
BONSOI-3981; No. of Pages 6
T.-H. Kim et al. / Joint Bone Spine xxx (2014) xxx–xxx
4
Table 3 The 50 top ranked genes that show a differential expression between responders and non-responders. GSE3592
GSE8350
GSE12051
GSE20690
GSE33377
LOC728392 SYT7 PIGA PSMA5 REPS2 SEPP1 MAPK11 MAP3K8 MUPCDH LOC647622 MAPKAP1 ACOT12 ATP6V1D RBM25 FKBP2 SRFBP1 MRPL22 NEDD4 ELP3 UCP2 RABEPK HBG1 SIAH2 UAP1 NRBP2 CARHSP1 HNRNPU DPY19L2P2 PANX1 C4BPA FZD1 HBBP1 PIGT CBARA1 AURKB LRRC32 FLJ12529 RRM2 SERP1 IFIT1 MFN2 RALA ATN1 YWHAE ADAMTS1 WDR13 COBL VAMP2 ASS1 SLC2A10
IL2RB CDKN1A TGFBR3 LTBP1 CD151 PTHR1 CKS2 EPS15 TGFB2 CYP2B6 FOLR2 HLA-DQA1 C1S PLSCR2 TBX21 CX3CR1 CP GZMA RPL41 IKBKB GZMB CD3E IFI30 TNFRSF11A FOXK2 CCL4 GPX2 HMGA1 CDKN2A GATA3 PTH SOAT1 COMT TNFSF11 HIST2H4A SELP ITGAL CCR3 HLA-DQB1 CD8B AOX1 PI4KA ITGAM YY1 PAK4 NEUROG1 PURA G0S2 NME1 CD2
SH2D1B HCK S100A12 RGS13 SPON2 NFKBIZ ACSL1 ENPP6 IL2RB CD247 CDR2L C1orf92 KRTAP6-2 GFOD1 SH2D2A ZNF24 OSMR KIR2DL5A MXD4 MTSS1 PROK2 LRCH3 CD160 UBL4B ZNF660 SLITRK2 GNLY KCNJ15 CD81 AQP9 HK3 EFHC1 CYTH4 RUNX1T1 DDX25 FOXP2 BCL6 FRAT2 SLC16A8 YWHAQ TRIM23 BASP1 PTGDR GINS3 CREB5 MKNK1 SEPHS1 HLA-DRB3 PLIN3 FGFBP2
BSN TMPRSS11A OBSL1 YOD1 G6PC GPER SSTR1 LOC100287076 TMEM98 CCDC105 RPSA C9orf57 VN1R4 PPBP KLK13 EPB41L5 VSTM2A EDNRB CTNNAL1 VN1R1 SLC1A3 LOC389740 LOC100288276 ASB14 RNF182 ARG2 SCRG1 GAFA3 UGT1A6 LOC728218 CDY2A LOC84931 PSPH TPST1 ZIC1 ADAMTS19 LOC197350 FGF2 C5orf20 SCML2 PLEKHA5 ADIPOQ NEUROD4 KCNN3 MECOM RXFP1 ABHD10 KLRG2 CA6 CLDN8
LOC100129033 MOP-1 G0S2 CLIC3 HIRIP3 PTGS2 PPP1R15A PCGF5 PMAIP1 TPM1 PIGV EGR3 RAPGEF1 GRAMD1B EGR2 SESTD1 CSRNP1 ZBTB6 DUSP1 C7orf41 CMIP GTPBP2 TMOD2 SH2D2A ZNF2 PDE3A KHDRBS1 LOC100128035 RAMP3 GBP2 c15orf40 HSPC159 ATP6V0A2 NR4A2 MED12L YPEL5 UQCRFS1 HERPUD2 GPD2 VEPH1 TAP2 OR10A7 PSTPIP2 SLK PPCDC HEBP1 FRMD3 RNF11 FOS RANBP17
multiple-testing correction, we sorted the genes by their raw P-values and manually ranked them. Lists of the top 50 genes in each dataset are shown in Table 3 and Table S1 (see the supplementary material associated with this article online); three genes, IL2RB, SH2D2A and G0S2, appeared in more than 1 study. We used DAVID to perform enrichment analysis to identify common biological themes in genes with raw P-values < 0.01 in each dataset. Lists of the top 20 enriched annotation terms associated with each dataset are presented in Table S2. The enriched annotation terms identified showed little overlap between datasets. In addition, we performed enriched pathway analysis using Kyoto Encyclopedia of Genes and Genomes (KEGG) in 247 genes, all the top 50 genes of each dataset. Lists of significant enriched pathways in these 247 genes are presented in Table S3. 3.3. Meta-analysis Among 5 datasets selected for individual analysis, we included 4 in the meta-analysis: GSE3592, GSE12051, GSE20690, and GSE33377. One dataset, GSE8350, was excluded from analysis
because of a low-density cDNA microarray covering a total of only 747 genes. Meta-analysis was performed using 2 R packages, MetaQC and MetaDE, as described previously [24]. After gene matching based on official gene symbols, preprocessing, and filtering, we assessed the quality and consistency of microarray datasets by using the MetaQC package. QC measures and mean rank scores were obtained and are summarized in Table 4. All included studies were very heterogeneous, and the findings indicated a limited potential for meta-analysis. The GSE3592 dataset was excluded from meta-analysis because of a poor IQC score, which indicates that this dataset has a heterogeneous co-expression structure with other datasets. The remaining datasets were included in meta-analysis: GSE12051, GSE20690, and GSE33377. When a raw P-value of < 0.01 was used as a cutoff, meta-analysis using the Fisher method and the maximum P-value method identified 57 and 47 genes, respectively. The lists of identified genes are shown in Table S4. We performed enriched pathway analysis using Kyoto Encyclopedia of Genes and Genomes (KEGG) in these lists of identified genes by meta-analysis. Arachidonic acid metabolism including AKR1C3, PTGS2, PTGDS was only one significant enriched
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
G Model
ARTICLE IN PRESS
BONSOI-3981; No. of Pages 6
T.-H. Kim et al. / Joint Bone Spine xxx (2014) xxx–xxx
5
Table 4 Results of QC measures and mean rank scores. Dataset
IQC
EQC
CQCg
CQCp
AQCg
AQCp
Mean score
GSE3592 GSE12051 GSE20690 GSE33377
1.3a 3.92 4.74 5.62
0.07a 0.2a 0.19a 1.35a
0.25a 0.0a 0.0a 0.01a
1.7a 13.75 7.73 3.17
0.52a 0.3a 0.48a 0.04a
0.12a 0.0a 0.0a 0.03a
2.50 2.50 2.83 2.17
a
Non-statistical significance and candidate of heterogeneous studies.
Table 5 Number of detected DEGs by two meta-analysis methods (Fisher and MaxP).
Raw P-value P < 0.01 P < 0.05 Adjusted P-value P < 0.05 P < 0.1
Fisher
MaxP
57 264
47 311
2 2
0 0
pathway in list of identified genes by meta-analysis using the Fisher method. In contrast, no significant enriched pathway was found in list of identified genes by meta-analysis using the maximum P-value method. After applying P-value adjustment for multiple-testing correction, two DEGs, S100A12 and G0S2, were identified by the Fisher method using MetaDE package (Table 5 and Table S5). However, S100A12 was down-regulated in responders of GSE20690 but up-regulated in responders of GSE12051. Although GSE8350 is a data of low-density microarray, it is possible that candidate genes were pre-selected with the aim to focus on relevant candidates. Therefore, we performed the meta-analysis including GSE8350. Only one DEG, G0S2, was identified by meta-analysis including GSE8350. 4. Discussion We performed a comprehensive literature search to identify microarray-based studies of gene expression profiles predicting response to anti-TNF therapy in RA patients and then attempted to find a DEG profile to discriminate between different response groups by combining the datasets of these studies. To our knowledge, this is the first study in which publicly available microarray datasets for anti-TNF response were combined and re-examined by means of the same analytical tool such as GEO2R and metaanalysis tools. We identified and analyzed 6 studies and 5 microarray datasets. In most included studies, RNA was isolated from whole blood. Because cell isolation procedure is a welldocumented source of variability in gene expression analysis, most transcriptional studies are now carried out using whole blood. Cell-type-specific gene expression profiling analyses are very useful approaches to identify genes specifically expressed in certain cell types that play an important role in the pathogenesis of RA [5]. However, gene expressions in peripheral blood cells are highly sensitive to the stresses, and therefore the expression patterns of genes might be altered during cell isolation and purification. Also, cell purification might be a source of bias by altering the steady state transcriptome in positive selection method or contamination of other cell types in enrichment purification method [26]. The expression profiles used in these studies for predicting response to anti-TNF therapy were not consistent with each other, and 2 genes appeared in more than 1 study. There are several reasons for the discordance of these studies: differences in microarray platforms, quality of microarray results, methods of analysis, tissues used for analysis, types of anti-TNF therapy, and response criteria for defining treatment outcome. To compensate for the influence of different analytical methods on gene expression
profiles, we analyzed each dataset individually by using the same tool, GEO2R. We ranked the genes for each study, and 2 genes, IL2RB and G0S2, appeared in the top 50 of more than 1 study. To increase statistical power, we combined the results of individual studies by using a meta-analysis approach [23,27–29]. Two DEGs, S100A12 and G0S2, were identified by the Fisher method, but not by the maximum P-value method using MetaDE package. However, S100A12 was down-regulated in responders of GSE20690 but up-regulated in responders of GSE12051, indicating that S100A12 is not a true predictive marker of anti-TNF response. Finally, we identified G0S2 as DEG predicting response to anti-TNF therapy. G0S2 is involved in cell proliferation, apoptosis, inflammation, metabolism, and carcinogenesis [30]. The roles and functions of G0S2 in RA have not been yet studied. Previous only one study reported that gene expression of G0S2 is increased in both the bone marrow-derived mononuclear cells (BMMCs) and peripheral blood mononuclear cells (PBMCs) of RA patients, compared to OA patients [31]. Meta-analysis is a statistical approach that combines results from independent but related studies and increases both the statistical power and generalizability of single-study analysis. We followed a stepwise approach suggested by Ramasamy et al. for conducting a meta-analysis of microarray datasets. One of the most important steps in a successful meta-analysis is to identify and remove any microarray datasets that are of poor quality [23,25]. To assess the quality of microarray datasets included in this meta-analysis, we used the MetaQC package [25]. The included studies were very heterogeneous with respect to the calculated QC measures and different microarray platforms. CQC and AQC quantify the reproducibility of DEGs or pathways detected in an individual study, as compared to those detected by meta-analysis of all other studies. Most CQCg, AQCg, and AQCp scores of the included studies were not statistically significant, indicating that the DEGs and pathways detected in individual studies were relatively less reproducible and this is major limitation of our study. Compared to other datasets, the GSE3592 dataset had a worse IQC score. Therefore, we excluded this dataset from meta-analysis. Lequerre et al. used a custom-made cDNA microarray platform and their sample size was small, which might explain the low IQC score. Following quality control assessment by MetaQC, metaanalysis was performed in three datasets, GSE12051, GSE20690 and GSE33377, which show the internal homogeneity of co-expression structure among studies. Several types of meta-analysis methods were used to combine information for DEG detection in microarray datasets, including combining P-values, combining ranks, combining effect sizes, and vote counting [23,28]. We used 2 popular P-value combination methods, the Fisher method and the maximum P-value method. The Fisher method detects genes that are differentially expressed in 1 or more studies, while the maximum P-value method detects genes that are differentially expressed in all studies [28,29]. Inflamed synovial tissues reflect pathogenesis at an affected joint more directly than peripheral blood cells; however, obtaining such samples is difficult. Thus, the number of studies using synovial tissues was limited. Therefore, our analysis only included studies that analyzed gene expression in blood cells. In two studies
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
G Model BONSOI-3981; No. of Pages 6
ARTICLE IN PRESS T.-H. Kim et al. / Joint Bone Spine xxx (2014) xxx–xxx
6
that used synovial tissue as samples, the patients with transcription profiles associated with high inflammatory activity were more likely to respond to anti-TNF therapy [32,33]. In contrast, Badot et al. showed that genes involved in cell division and regulation of immune responses were over-expressed at baseline in poor responders [18]. In the latter study, they showed that genes over-expressed in poor responders are induced in fibroblast-like synoviocytes (FLS) by several cytokines, including TNF-␣, IL-1, and IL-17. The studies that used synovial tissues as samples reported inconsistent gene expression profiles. Furthermore, a significant gene expression signature predictive of treatment outcome was not identified in a recent, large study that used synovial tissue as sample from 62 RA patients [19]. Gene expression profiling has been extensively investigated as a potential source of biomarkers for diagnosis, prognosis, and prediction of treatment outcome, and it is already being used in clinical settings. Array-based expression profiling has been applied to find RA biomarkers to aid in diagnostic and prognostic classification, quantify disease activity, and predict treatment outcome [6]. In this study, we performed a systematic literature search to identify array-based expression profiling studies of genes predicting response to anti-TNF therapy and attempted to combine and analyze these datasets by using a meta-analysis approach. To our knowledge, this is the first study to use meta-analysis to integrate publicly available microarray datasets related to the response of anti-TNF therapy in RA patients. In this study, we could identify one DEG, G0S2 by analyzing microarray datasets, which are publicly available, using GEO2R and performing a meta-analysis of these datasets. However, microarray-based studies included in our analysis were performed in small samples. Therefore, further investigations based on larger studies are needed to confirm the significance of G0S2 in predicting response to anti-TNF therapy. Disclosure of interest The authors declare that they have no conflicts of interest concerning this article. Acknowledgement This research was supported by a Korea University Grant. Appendix A. Supplementary data Supplementary data (tables S1–S5) associated with this article can be found, in the online version, at doi:10.1016/j.jbspin. 2014.01.013. References [1] McInnes IB, Schett G. Cytokines in the pathogenesis of rheumatoid arthritis. Nat Rev Immunol 2007;7:429–42. [2] Kievit W, Adang EM, Fransen J, et al. The effectiveness and medication costs of three anti-tumour necrosis factor alpha agents in the treatment of rheumatoid arthritis from prospective clinical practice data. Ann Rheum Dis 2008;67:1229–34. [3] Emery P, Dorner T. Optimising treatment in rheumatoid arthritis: a review of potential biological markers of response. Ann Rheum Dis 2011;70:2063–70. [4] Kooloos WM, de Jong DJ, Huizinga TW, et al. Potential role of pharmacogenetics in anti-TNF treatment of rheumatoid arthritis and Crohn’s disease. Drug Discov Today 2007;12:125–31. [5] Bauer JW, Bilgic H, Baechler EC. Gene-expression profiling in rheumatic disease: tools and therapeutic potential. Nat Rev Rheumatol 2009;5:257–65.
[6] Haupl T, Stuhlmuller B, Grutzkau A, et al. Does gene expression analysis inform us in rheumatoid arthritis? Ann Rheum Dis 2010;69:i37–42. [7] Bild AH, Yao G, Chang JT, et al. Oncogenic pathway signatures in human cancers as a guide to targeted therapies. Nature 2006;439:353–7. [8] Ebert BL, Galili N, Tamayo P, et al. An erythroid differentiation signature predicts response to lenalidomide in myelodysplastic syndrome. PLoS Med 2008;5:e35. [9] Holleman A, Cheok MH, den Boer ML, et al. Gene-expression patterns in drugresistant acute lymphoblastic leukemia cells and response to treatment. N Engl J Med 2004;351:533–42. [10] Pomeroy SL, Tamayo P, Gaasenbeek M, et al. Prediction of central nervous system embryonal tumour outcome based on gene expression. Nature 2002;415:436–42. [11] van’t Veer LJ, Dai H, van de Vijver MJ, et al. Gene expression profiling predicts clinical outcome of breast cancer. Nature 2002;415:530–6. [12] Julia A, Erra A, Palacio C, et al. An eight-gene blood expression profile predicts the response to infliximab in rheumatoid arthritis. PLoS One 2009;4:e7556. [13] Lequerre T, Gauthier-Jauneau AC, Bansard C, et al. Gene profiling in white blood cells predicts infliximab responsiveness in rheumatoid arthritis. Arthritis Res Ther 2006;8:R105. [14] Sekiguchi N, Kawauchi S, Furuya T, et al. Messenger ribonucleic acid expression profile in peripheral blood cells from RA patients following treatment with an anti-TNF-alpha monoclonal antibody, infliximab. Rheumatology (Oxford) 2008;47:780–8. [15] Stuhlmuller B, Haupl T, Hernandez MM, et al. CD11c as a transcriptional biomarker to predict response to anti-TNF monotherapy with adalimumab in patients with rheumatoid arthritis. Clin Pharmacol Ther 2010;87:311–21. [16] Tanino M, Matoba R, Nakamura S, et al. Prediction of efficacy of anti-TNF biologic agent, infliximab, for rheumatoid arthritis patients using a comprehensive transcriptome analysis of white blood cells. Biochem Biophys Res Commun 2009;387:261–5. [17] Toonen EJ, Gilissen C, Franke B, et al. Validation study of existing gene expression signatures for anti-TNF treatment in patients with rheumatoid arthritis. PLoS One 2012;7:e33199. [18] Badot V, Galant C, Nzeusseu Toukap A, et al. Gene expression profiling in the synovium identifies a predictive signature of absence of response to adalimumab therapy in rheumatoid arthritis. Arthritis Res Ther 2009;11:R57. [19] Lindberg J, Wijbrandts CA, van Baarsen LG, et al. The gene expression profile in the synovium as a predictor of the clinical response to infliximab treatment in rheumatoid arthritis. PLoS One 2010;5:e11310. [20] van Baarsen LG, Wijbrandts CA, Rustenburg F, et al. Regulation of IFN response gene activity during infliximab treatment in rheumatoid arthritis is associated with clinical response to treatment. Arthritis Res Ther 2010;12:R11. [21] Barrett T, Wilhite SE, Ledoux P, et al. NCBI GEO: archive for functional genomics data sets – update. Nucleic Acids Res 2012;41:D991–5. [22] Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009;4:44–57. [23] Ramasamy A, Mondry A, Holmes CC, et al. Key issues in conducting a metaanalysis of gene expression microarray datasets. PLoS Med 2008;5:e184. [24] Wang X, Kang DD, Shen K, et al. An R package suite for microarray metaanalysis in quality control, differentially expressed gene analysis and pathway enrichment detection. Bioinformatics 2012;28:2534–6. [25] Kang DD, Sibille E, Kaminski N, et al. MetaQC: objective quality control and inclusion/exclusion criteria for genomic meta-analysis. Nucleic Acids Res 2012;40:e15. [26] Pascual V, Chaussabel D, Banchereau J. A genomic approach to human autoimmune diseases. Annu Rev Immunol 2010;28:535–71. [27] Cahan P, Rovegno F, Mooney D, et al. Meta-analysis of microarray results: challenges, opportunities, and recommendations for standardization. Gene 2007;401:12–8. [28] Tseng GC, Ghosh D, Feingold E. Comprehensive literature review and statistical considerations for microarray meta-analysis. Nucleic Acids Res 2012;40:3785–99. [29] Wang X, Lin Y, Song C, et al. Detecting disease-associated genes with confounding variable adjustment and the impact on genomic meta-analysis: with application to major depressive disorder. BMC Bioinformatics 2012;13:52. [30] Heckmann BL, Zhang X, Xie X, et al. The G0/G1 switch gene 2 (G0S2): regulating metabolism and beyond. Biochim Biophys Acta 2013;1831:276–81. [31] Nakamura N, Shimaoka Y, Tougan T, et al. Isolation and expression profiling of genes upregulated in bone marrow-derived mononuclear cells of rheumatoid arthritis patients. DNA Res 2006;13:169–83. [32] Lindberg J, af Klint E, Catrina AI, et al. Effect of infliximab on mRNA expression profiles in synovial tissue of rheumatoid arthritis patients. Arthritis Res Ther 2006;8:R179. [33] van der Pouw Kraan TC, Wijbrandts CA, van Baarsen LG, et al. Responsiveness to anti-tumour necrosis factor alpha therapy is related to pre-treatment tissue inflammation levels in rheumatoid arthritis patients. Ann Rheum Dis 2008;67:563–6.
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
ARTICLE IN PRESS Joint Bone Spine xxx (2014) xxx–xxx
Available online at www.sciencedirect.com
Original article
Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets Tae-Hwan Kim a , Sung Jae Choi b , Young Ho Lee b , Gwan Gyu Song b , Jong Dae Ji b,∗ a b
Hanyang University Hospital for Rheumatic Diseases, Seoul, South Korea Division of Rheumatology, College of Medicine, Korea University, 126-1, Anam-Dong 5-Ga, Sungbuk-Ku, Seoul 136-705, South Korea
a r t i c l e
i n f o
Article history: Accepted 15 January 2014 Available online xxx Keywords: Rheumatoid arthritis Anti-TNF Response Microarray
a b s t r a c t Objectives: Anti-tumor necrosis factor (TNF) therapy is the treatment of choice for rheumatoid arthritis (RA) patients in whom standard disease-modifying anti-rheumatic drugs are ineffective. However, a substantial proportion of RA patients treated with anti-TNF agents do not show a significant clinical response. Therefore, biomarkers predicting response to anti-TNF agents are needed. Recently, gene expression profiling has been applied in research for developing such biomarkers. Methods: We compared gene expression profiles reported by previous studies dealing with the responsiveness of anti-TNF therapy in RA patients and attempted to identify differentially expressed genes (DEGs) that discriminated between responders and non-responders to anti-TNF therapy. We used microarray datasets available at the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO). Results: This analysis included 6 studies and 5 sets of microarray data that used peripheral blood samples for identification of DEGs predicting response to anti-TNF therapy. We found little overlap in the DEGs that were highly ranked in each study. Three DEGs including IL2RB, SH2D2A and G0S2 appeared in more than 1 study. In addition, a meta-analysis designed to increase statistical power found one DEG, G0S2 by the Fisher’s method. Conclusion: Our finding suggests the possibility that G0S2 plays as a biomarker to predict response to anti-TNF therapy in patients with rheumatoid arthritis. Further investigations based on larger studies are therefore needed to confirm the significance of G0S2 in predicting response to anti-TNF therapy. © 2014 Société franc¸aise de rhumatologie. Published by Elsevier Masson SAS. All rights reserved.
1. Introduction Rheumatoid arthritis (RA) is a systemic autoimmune disease characterized by chronic inflammation of joints followed by destruction of cartilage and bone. The precise etiology of RA is not completely understood, but a large body of studies indicates that pro-inflammatory cytokines such as TNF␣ and IL-1 are significantly involved in the pathophysiology of joint inflammation and bone destruction in RA [1]. Biologic agents, especially anti-TNF agents, have produced significant and clinically relevant improvements in patients with active RA and appear to be very effective therapeutic molecules for RA patients. Although the use of anti-TNF agents has revolutionized the treatment of RA, a substantial proportion (up to one-third) of RA patients treated with anti-TNF agents fail to achieve a satisfactory clinical response [2]. Therefore, a strategy for predicting which patients will and will not respond prior to initiation of therapy is needed to prevent disease progression,
∗ Corresponding author. Tel.: +82 2 920 5489; fax: +82 2 920 5974. E-mail address: [email protected] (J.D. Ji).
unnecessary side effects, and high costs for non-responders. Although baseline demographics and clinical characteristics, including age, sex, disease duration, Disease Activity Score (DAS28), Health Assessment Questionnaire score, and concurrent use of disease-modifying anti-rheumatic drugs (DMARDs), have been extensively studied, potential prediction markers of therapeutic responsiveness have not yet been discovered [3]. Pharmacogenetic studies have also been performed to discover determinants of efficacy of anti-TNF agents. Some studies have shown that polymorphisms in genes encoding TNF␣, TNF receptors, the major histocompatibility complex region, and the Fc␥ receptor IIIA (FC␥RIIIA) are associated with a favorable response to anti-TNF therapy; however, these results were not consistent with those of other studies [3,4]. Gene expression profiling with microarrays has become a standard method for identifying the genes and biological pathways that are associated with various complex diseases [5,6]. Biomarkers identified by gene expression profiling can be used to make a diagnosis, monitor disease activity, and predict response to therapy. This technology has already been shown to identify gene expression profiles predictive of treatment response and clinical outcomes
1297-319X/$ – see front matter © 2014 Société franc¸aise de rhumatologie. Published by Elsevier Masson SAS. All rights reserved. doi:10.1016/j.jbspin.2014.01.013
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
G Model BONSOI-3981; No. of Pages 6
ARTICLE IN PRESS T.-H. Kim et al. / Joint Bone Spine xxx (2014) xxx–xxx
2
expression profiles predicting response to anti-TNF therapy were included in the analysis. 2.2. Individual analysis
Fig. 1. Flow chart of selection process for meta-analysis. A. Available studies were identified via Medline/PubMed searches. B. Available GEO datasets were identified through GEO datasets searching.
for some diseases [7–11]. Recently, several studies have attempted to identify gene expression profiles predicting the response to anti-TNF therapy in RA patients by using multiple gene expression microarrays [12–17]. Although these studies have successfully identified gene expression signatures predicting the response to anti-TNF therapy, the expression profiles identified in these studies were not consistent with each other, and the genes reported by the authors of each study had little overlap. Herein, we attempt to identify differentially expressed genes (DEGs) that differ between responders and non-responders to anti-TNF therapy by: • comparing gene signatures previously reported to predict response to anti-TNF therapy in RA patients; • analyzing microarray datasets of these studies, which are available at the National Center for Biotechnology Information’s Gene Expression Omnibus (GEO), using GEO2R; • performing a meta-analysis of these datasets. 2. Methods 2.1. Data collection We performed literature survey via Medline/PubMed searches to identify available studies that examined gene expression profiles of peripheral blood samples to predict treatment outcome of anti-TNF therapy in RA patients. We searched for a combination of KEY words such as “Infliximab”, “Adalimumab”, “Etanercept”, “Microarray”, “Gene expression profile”, “Anti-TNF” and “Rheumatoid arthritis” as both medical subject heading terms and text words. We also tried to identify additional studies by manually searching references of the original articles or review articles on this topic. Twenty-two relevant studies that investigated gene expression profiling in RA patients treated with anti-TNF agents were identified using PubMed. On the basis of their titles and abstracts, 9 studies were selected for full-text review (Fig. 1A). Three studies was excluded from analysis because they used synovial tissues for examining gene expression signatures [18,19] or did not focus on a predictor of the clinical response to anti-TNF therapy [20]. The remaining 6 studies were included in our analysis [12–17]. In addition, 10 sets of microarray data were identified by searching GEO DataSets. Two datasets, GSE15602 and GSE21537, were excluded from analysis because they used synovial tissues, and 3 datasets, GSE7524, GSE12897, and GSE19821, were excluded because they did not focus on a predictor of clinical response to anti-TNF therapy (Fig. 1B). Among 6 including studies, the dataset of one study was not deposited in the GEO database [15]. Finally, 6 studies and 5 sets of microarray data (GSE3592, GSE8350, GSE12051, GSE20690, and GSE33377) that used peripheral blood samples to identify gene
We analyzed each dataset individually by using the same analytical tool, GEO2R (GEO’s online tool for analyzing GEO data, available at http://www.nci.nlm.nih.gov/geo/geo2r/), to maintain consistency during individual analysis. GEO2R is an interactive, online tool for R-based analysis of GEO data to identify genes that are differentially expressed across experimental conditions [21]. GEO2R analyzes original, submitter-supplied data with the GEOquery and limma (Linear Models for Microarray Analysis) R packages from the Bioconductor project. The basic statistic used in limma for significance analysis is the moderated t-statistic. The Benjamini & Hochberg false discovery rate method is used to apply P-value adjustment for multiple-testing correction. The P-value and adjusted P-value were used to identify a list of DEGs from each dataset. Each list of DEGs was used for pathway enrichment analysis. This analysis was done using the Database for Annotation, Visualization and Integrated Discovery (DAVID) (http://david.abcc.ncifcrf.gov/home.jsp) [22]. Several annotation categories, including gene ontology (GO) terms, BIOCARTA pathways, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, were subjected to enrichment analysis. 2.3. Meta-analysis To combine the results of individual studies and obtain a list of more robust DEGs between responders and non-responders to antiTNF therapy, we followed the guideline provided by Ramasamy et al. for meta-analysis of gene expression microarray datasets [23]. We used 2 R packages, MetaQC and MetaDE, for quality control (QC) and identification of DEGs, respectively [24]. MetaQC implements 6 quantitative quality control measures, including internal quality control (IQC, internal homogeneity of co-expression structure among studies), external quality control (EQC, external consistency of co-expression structure correlating with pathway database), accuracy quality control (AQCg and AQCp, accuracy of DEG detection or pathway identification), and consistency quality control (CQCg and CQCp, consistency of differential expression ranking for genes or pathways) [25]. In addition, the mean rank of all QC measures in each dataset was computed as a quantitative summary score by calculating the ranks of each QC measure among all included datasets. The moderated t-statistic was used to calculate P-values in each dataset, and a meta-analysis was performed with the MetaDE package to identify DEGs by the Fisher method and the maximum P-value method. 3. Results 3.1. Study characteristics The analyses included 6 studies that used peripheral blood samples to identify gene expression profiles predicting response to anti-TNF therapy. The dataset of 1 of these studies was not deposited in the GEO database [15]. Therefore, the microarray dataset of this study was not included for individual analysis or meta-analysis. Detailed information of each study is described in Table 1. In the study of Lequerré et al., patients with very active disease (Disease Activity Score 28 > 5.1) were included in this study and mononuclear cell RNAs were collected at baseline. In the study of Sekiguchi et al., patients having inadequate control of RA symptoms were included and whole blood RNA was extracted before the infusion of infliximab. In the study of Julia et al., patients with active disease (Disease Activity Score > 3.2) were included and
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
G Model
ARTICLE IN PRESS
BONSOI-3981; No. of Pages 6
T.-H. Kim et al. / Joint Bone Spine xxx (2014) xxx–xxx
3
Table 1 Summary information of studies included in analysis. Study
Toonen et al. [17] Stuhlmüller et al. [15] Tanino et al. [16] Julia et al. [12] Sekiguchi et al. [14] Lequerré et al. [13]
Year
2012 2010 2009 2009 2008 2006
Platform
HuEx-1 0-st HG-U133A Agilent-014850 Illumina Human-6 v1 Customized Customized
Number of samples
Responder
Non-responder
18 5 43 37 8 6
24 2 25 7 10 7
Response determination
Response criteria
GEO dataset
Tissue
Anti-TNF
14 weeks – 14 weeks 14 weeks 22 weeks 3 months
EULAR ACR CRP EULAR ACR DAS28
GSE33377 – GSE20690 GSE12051 GSE8350 GSE3592
Whole blood Monocyte Whole blood Whole blood Whole blood PBMC
I, A A I I I I
I: infliximab; A: adalimumab.
whole blood samples were extracted before treatment. In the study of Tanino et al., patients were resistant to standard methotrexate treatment and whole blood samples were collected before treatment. In the study of Toonen et al., patients with active disease (Disease Activity Score > 3.2) and previous failure on at least two disease-modifying anti-rheumatic drugs (DMARDs) were included and whole blood was sampled before treatment. Gene expression signatures predictive of anti-TNF treatment outcome were reported in 5 studies, and Toonen et al. tried to validate previously reported gene expression signatures using microarray data on whole blood from RA patients treated with anti-TNF agents. We observed several clinical and experimental variations between these studies, including differences in microarray platform, types of anti-TNF therapy,
response criteria, and tissues used for analysis (whole blood versus monocytes). The gene expression profiles used in different studies were not consistent with each other, and only 2 genes, IL2RB and AP1S2, were used in more than 1 study (Table 2). 3.2. Individual analysis The included datasets were individually analyzed with GEO2R. When the Benjamini & Hochberg false discovery rate method was used to adjust statistics for multiple-testing (adjusted P-value < 0.05 was the cutoff for statistical significance), DEGs for different response groups were not identified in any of the datasets. Because no genes passed the cutoff for statistical significance after
Table 2 Transcript sets of studies included in analysis. Lequerré et al. (20 genes)
Sekiguchi et al. (18 genes)
Julia et al. (8 genes)
Tanino et al. (8 genes)
Stuhlmüller et al. (82 genes)a
AKAP9 COX7A2L ELMOD2 EPS15 FBXO5 HLA-DPB1 LAMR1 MCP MRPL22 MTCBP-1 PFKFB4 PSMB9 PTPN12 QIL1 RASGRP3 RPL35 RPS16 RPS28 SCAM-1 TBL2
AP1S2 CX3CR1 EIF2AK2 HLA-DQA1 IFI6 IFIT1 IFIT3 IFITM1 IFITM2 IGHM IL2RB MX1 MX2 OAS1 OAS2 OASL SRP9 ZFP36L2
CAMP GNLY HLA-DRB3 IL2RB MXD4 SH2D1B SLC2A3 TLR5
ANKRD55 ATP5I C21orf58 CLGN LOC643981 PSPH TBC1D8 TMEM141
ABCC3 ACTG1 ACTN1 AP1S2 APP AQP3 ATXN1 C14orf156 C20orf111 CALM1 CCDC91 CCL3L3 CD86 CD97 CDC2L6 CLC CLMN CNOT6 COX7B CPA3 CRYZ CTSZ CUGBP2 CYB5A DOCK5 DSC2 EDARADD EGR1 ENPP4 ERMAP ESD F2RL1 FAM3C FCER1A GATA2 GCLC GNA12 GNB1 GNE GNS
a
HLA-DQA2 HLA-DQB2 HMBOX1 IDH1 IL8 IL8RB ITGAX KIAA0746 KLF11 LBH LEPR MAL NADK NARF NDUFB1 NGRN NLRP3 PAIP1 PECAM1 PF4V1 PIM1 PLCL2 PTMAP7 RAPGEF2 RPS24 RPS26 SC4MOL SF3B1 SOD1 SUB1 TCF4 TM2D1 TMED2 TMEM176A TMEM176B TMEM45A TPD52 UCRC UQCRH WIPF1
Two Affymetrix ID (214807 at and AFFX-M27830 5 at) gene symbols not available.
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
G Model
ARTICLE IN PRESS
BONSOI-3981; No. of Pages 6
T.-H. Kim et al. / Joint Bone Spine xxx (2014) xxx–xxx
4
Table 3 The 50 top ranked genes that show a differential expression between responders and non-responders. GSE3592
GSE8350
GSE12051
GSE20690
GSE33377
LOC728392 SYT7 PIGA PSMA5 REPS2 SEPP1 MAPK11 MAP3K8 MUPCDH LOC647622 MAPKAP1 ACOT12 ATP6V1D RBM25 FKBP2 SRFBP1 MRPL22 NEDD4 ELP3 UCP2 RABEPK HBG1 SIAH2 UAP1 NRBP2 CARHSP1 HNRNPU DPY19L2P2 PANX1 C4BPA FZD1 HBBP1 PIGT CBARA1 AURKB LRRC32 FLJ12529 RRM2 SERP1 IFIT1 MFN2 RALA ATN1 YWHAE ADAMTS1 WDR13 COBL VAMP2 ASS1 SLC2A10
IL2RB CDKN1A TGFBR3 LTBP1 CD151 PTHR1 CKS2 EPS15 TGFB2 CYP2B6 FOLR2 HLA-DQA1 C1S PLSCR2 TBX21 CX3CR1 CP GZMA RPL41 IKBKB GZMB CD3E IFI30 TNFRSF11A FOXK2 CCL4 GPX2 HMGA1 CDKN2A GATA3 PTH SOAT1 COMT TNFSF11 HIST2H4A SELP ITGAL CCR3 HLA-DQB1 CD8B AOX1 PI4KA ITGAM YY1 PAK4 NEUROG1 PURA G0S2 NME1 CD2
SH2D1B HCK S100A12 RGS13 SPON2 NFKBIZ ACSL1 ENPP6 IL2RB CD247 CDR2L C1orf92 KRTAP6-2 GFOD1 SH2D2A ZNF24 OSMR KIR2DL5A MXD4 MTSS1 PROK2 LRCH3 CD160 UBL4B ZNF660 SLITRK2 GNLY KCNJ15 CD81 AQP9 HK3 EFHC1 CYTH4 RUNX1T1 DDX25 FOXP2 BCL6 FRAT2 SLC16A8 YWHAQ TRIM23 BASP1 PTGDR GINS3 CREB5 MKNK1 SEPHS1 HLA-DRB3 PLIN3 FGFBP2
BSN TMPRSS11A OBSL1 YOD1 G6PC GPER SSTR1 LOC100287076 TMEM98 CCDC105 RPSA C9orf57 VN1R4 PPBP KLK13 EPB41L5 VSTM2A EDNRB CTNNAL1 VN1R1 SLC1A3 LOC389740 LOC100288276 ASB14 RNF182 ARG2 SCRG1 GAFA3 UGT1A6 LOC728218 CDY2A LOC84931 PSPH TPST1 ZIC1 ADAMTS19 LOC197350 FGF2 C5orf20 SCML2 PLEKHA5 ADIPOQ NEUROD4 KCNN3 MECOM RXFP1 ABHD10 KLRG2 CA6 CLDN8
LOC100129033 MOP-1 G0S2 CLIC3 HIRIP3 PTGS2 PPP1R15A PCGF5 PMAIP1 TPM1 PIGV EGR3 RAPGEF1 GRAMD1B EGR2 SESTD1 CSRNP1 ZBTB6 DUSP1 C7orf41 CMIP GTPBP2 TMOD2 SH2D2A ZNF2 PDE3A KHDRBS1 LOC100128035 RAMP3 GBP2 c15orf40 HSPC159 ATP6V0A2 NR4A2 MED12L YPEL5 UQCRFS1 HERPUD2 GPD2 VEPH1 TAP2 OR10A7 PSTPIP2 SLK PPCDC HEBP1 FRMD3 RNF11 FOS RANBP17
multiple-testing correction, we sorted the genes by their raw P-values and manually ranked them. Lists of the top 50 genes in each dataset are shown in Table 3 and Table S1 (see the supplementary material associated with this article online); three genes, IL2RB, SH2D2A and G0S2, appeared in more than 1 study. We used DAVID to perform enrichment analysis to identify common biological themes in genes with raw P-values < 0.01 in each dataset. Lists of the top 20 enriched annotation terms associated with each dataset are presented in Table S2. The enriched annotation terms identified showed little overlap between datasets. In addition, we performed enriched pathway analysis using Kyoto Encyclopedia of Genes and Genomes (KEGG) in 247 genes, all the top 50 genes of each dataset. Lists of significant enriched pathways in these 247 genes are presented in Table S3. 3.3. Meta-analysis Among 5 datasets selected for individual analysis, we included 4 in the meta-analysis: GSE3592, GSE12051, GSE20690, and GSE33377. One dataset, GSE8350, was excluded from analysis
because of a low-density cDNA microarray covering a total of only 747 genes. Meta-analysis was performed using 2 R packages, MetaQC and MetaDE, as described previously [24]. After gene matching based on official gene symbols, preprocessing, and filtering, we assessed the quality and consistency of microarray datasets by using the MetaQC package. QC measures and mean rank scores were obtained and are summarized in Table 4. All included studies were very heterogeneous, and the findings indicated a limited potential for meta-analysis. The GSE3592 dataset was excluded from meta-analysis because of a poor IQC score, which indicates that this dataset has a heterogeneous co-expression structure with other datasets. The remaining datasets were included in meta-analysis: GSE12051, GSE20690, and GSE33377. When a raw P-value of < 0.01 was used as a cutoff, meta-analysis using the Fisher method and the maximum P-value method identified 57 and 47 genes, respectively. The lists of identified genes are shown in Table S4. We performed enriched pathway analysis using Kyoto Encyclopedia of Genes and Genomes (KEGG) in these lists of identified genes by meta-analysis. Arachidonic acid metabolism including AKR1C3, PTGS2, PTGDS was only one significant enriched
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
G Model
ARTICLE IN PRESS
BONSOI-3981; No. of Pages 6
T.-H. Kim et al. / Joint Bone Spine xxx (2014) xxx–xxx
5
Table 4 Results of QC measures and mean rank scores. Dataset
IQC
EQC
CQCg
CQCp
AQCg
AQCp
Mean score
GSE3592 GSE12051 GSE20690 GSE33377
1.3a 3.92 4.74 5.62
0.07a 0.2a 0.19a 1.35a
0.25a 0.0a 0.0a 0.01a
1.7a 13.75 7.73 3.17
0.52a 0.3a 0.48a 0.04a
0.12a 0.0a 0.0a 0.03a
2.50 2.50 2.83 2.17
a
Non-statistical significance and candidate of heterogeneous studies.
Table 5 Number of detected DEGs by two meta-analysis methods (Fisher and MaxP).
Raw P-value P < 0.01 P < 0.05 Adjusted P-value P < 0.05 P < 0.1
Fisher
MaxP
57 264
47 311
2 2
0 0
pathway in list of identified genes by meta-analysis using the Fisher method. In contrast, no significant enriched pathway was found in list of identified genes by meta-analysis using the maximum P-value method. After applying P-value adjustment for multiple-testing correction, two DEGs, S100A12 and G0S2, were identified by the Fisher method using MetaDE package (Table 5 and Table S5). However, S100A12 was down-regulated in responders of GSE20690 but up-regulated in responders of GSE12051. Although GSE8350 is a data of low-density microarray, it is possible that candidate genes were pre-selected with the aim to focus on relevant candidates. Therefore, we performed the meta-analysis including GSE8350. Only one DEG, G0S2, was identified by meta-analysis including GSE8350. 4. Discussion We performed a comprehensive literature search to identify microarray-based studies of gene expression profiles predicting response to anti-TNF therapy in RA patients and then attempted to find a DEG profile to discriminate between different response groups by combining the datasets of these studies. To our knowledge, this is the first study in which publicly available microarray datasets for anti-TNF response were combined and re-examined by means of the same analytical tool such as GEO2R and metaanalysis tools. We identified and analyzed 6 studies and 5 microarray datasets. In most included studies, RNA was isolated from whole blood. Because cell isolation procedure is a welldocumented source of variability in gene expression analysis, most transcriptional studies are now carried out using whole blood. Cell-type-specific gene expression profiling analyses are very useful approaches to identify genes specifically expressed in certain cell types that play an important role in the pathogenesis of RA [5]. However, gene expressions in peripheral blood cells are highly sensitive to the stresses, and therefore the expression patterns of genes might be altered during cell isolation and purification. Also, cell purification might be a source of bias by altering the steady state transcriptome in positive selection method or contamination of other cell types in enrichment purification method [26]. The expression profiles used in these studies for predicting response to anti-TNF therapy were not consistent with each other, and 2 genes appeared in more than 1 study. There are several reasons for the discordance of these studies: differences in microarray platforms, quality of microarray results, methods of analysis, tissues used for analysis, types of anti-TNF therapy, and response criteria for defining treatment outcome. To compensate for the influence of different analytical methods on gene expression
profiles, we analyzed each dataset individually by using the same tool, GEO2R. We ranked the genes for each study, and 2 genes, IL2RB and G0S2, appeared in the top 50 of more than 1 study. To increase statistical power, we combined the results of individual studies by using a meta-analysis approach [23,27–29]. Two DEGs, S100A12 and G0S2, were identified by the Fisher method, but not by the maximum P-value method using MetaDE package. However, S100A12 was down-regulated in responders of GSE20690 but up-regulated in responders of GSE12051, indicating that S100A12 is not a true predictive marker of anti-TNF response. Finally, we identified G0S2 as DEG predicting response to anti-TNF therapy. G0S2 is involved in cell proliferation, apoptosis, inflammation, metabolism, and carcinogenesis [30]. The roles and functions of G0S2 in RA have not been yet studied. Previous only one study reported that gene expression of G0S2 is increased in both the bone marrow-derived mononuclear cells (BMMCs) and peripheral blood mononuclear cells (PBMCs) of RA patients, compared to OA patients [31]. Meta-analysis is a statistical approach that combines results from independent but related studies and increases both the statistical power and generalizability of single-study analysis. We followed a stepwise approach suggested by Ramasamy et al. for conducting a meta-analysis of microarray datasets. One of the most important steps in a successful meta-analysis is to identify and remove any microarray datasets that are of poor quality [23,25]. To assess the quality of microarray datasets included in this meta-analysis, we used the MetaQC package [25]. The included studies were very heterogeneous with respect to the calculated QC measures and different microarray platforms. CQC and AQC quantify the reproducibility of DEGs or pathways detected in an individual study, as compared to those detected by meta-analysis of all other studies. Most CQCg, AQCg, and AQCp scores of the included studies were not statistically significant, indicating that the DEGs and pathways detected in individual studies were relatively less reproducible and this is major limitation of our study. Compared to other datasets, the GSE3592 dataset had a worse IQC score. Therefore, we excluded this dataset from meta-analysis. Lequerre et al. used a custom-made cDNA microarray platform and their sample size was small, which might explain the low IQC score. Following quality control assessment by MetaQC, metaanalysis was performed in three datasets, GSE12051, GSE20690 and GSE33377, which show the internal homogeneity of co-expression structure among studies. Several types of meta-analysis methods were used to combine information for DEG detection in microarray datasets, including combining P-values, combining ranks, combining effect sizes, and vote counting [23,28]. We used 2 popular P-value combination methods, the Fisher method and the maximum P-value method. The Fisher method detects genes that are differentially expressed in 1 or more studies, while the maximum P-value method detects genes that are differentially expressed in all studies [28,29]. Inflamed synovial tissues reflect pathogenesis at an affected joint more directly than peripheral blood cells; however, obtaining such samples is difficult. Thus, the number of studies using synovial tissues was limited. Therefore, our analysis only included studies that analyzed gene expression in blood cells. In two studies
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
G Model BONSOI-3981; No. of Pages 6
ARTICLE IN PRESS T.-H. Kim et al. / Joint Bone Spine xxx (2014) xxx–xxx
6
that used synovial tissue as samples, the patients with transcription profiles associated with high inflammatory activity were more likely to respond to anti-TNF therapy [32,33]. In contrast, Badot et al. showed that genes involved in cell division and regulation of immune responses were over-expressed at baseline in poor responders [18]. In the latter study, they showed that genes over-expressed in poor responders are induced in fibroblast-like synoviocytes (FLS) by several cytokines, including TNF-␣, IL-1, and IL-17. The studies that used synovial tissues as samples reported inconsistent gene expression profiles. Furthermore, a significant gene expression signature predictive of treatment outcome was not identified in a recent, large study that used synovial tissue as sample from 62 RA patients [19]. Gene expression profiling has been extensively investigated as a potential source of biomarkers for diagnosis, prognosis, and prediction of treatment outcome, and it is already being used in clinical settings. Array-based expression profiling has been applied to find RA biomarkers to aid in diagnostic and prognostic classification, quantify disease activity, and predict treatment outcome [6]. In this study, we performed a systematic literature search to identify array-based expression profiling studies of genes predicting response to anti-TNF therapy and attempted to combine and analyze these datasets by using a meta-analysis approach. To our knowledge, this is the first study to use meta-analysis to integrate publicly available microarray datasets related to the response of anti-TNF therapy in RA patients. In this study, we could identify one DEG, G0S2 by analyzing microarray datasets, which are publicly available, using GEO2R and performing a meta-analysis of these datasets. However, microarray-based studies included in our analysis were performed in small samples. Therefore, further investigations based on larger studies are needed to confirm the significance of G0S2 in predicting response to anti-TNF therapy. Disclosure of interest The authors declare that they have no conflicts of interest concerning this article. Acknowledgement This research was supported by a Korea University Grant. Appendix A. Supplementary data Supplementary data (tables S1–S5) associated with this article can be found, in the online version, at doi:10.1016/j.jbspin. 2014.01.013. References [1] McInnes IB, Schett G. Cytokines in the pathogenesis of rheumatoid arthritis. Nat Rev Immunol 2007;7:429–42. [2] Kievit W, Adang EM, Fransen J, et al. The effectiveness and medication costs of three anti-tumour necrosis factor alpha agents in the treatment of rheumatoid arthritis from prospective clinical practice data. Ann Rheum Dis 2008;67:1229–34. [3] Emery P, Dorner T. Optimising treatment in rheumatoid arthritis: a review of potential biological markers of response. Ann Rheum Dis 2011;70:2063–70. [4] Kooloos WM, de Jong DJ, Huizinga TW, et al. Potential role of pharmacogenetics in anti-TNF treatment of rheumatoid arthritis and Crohn’s disease. Drug Discov Today 2007;12:125–31. [5] Bauer JW, Bilgic H, Baechler EC. Gene-expression profiling in rheumatic disease: tools and therapeutic potential. Nat Rev Rheumatol 2009;5:257–65.
[6] Haupl T, Stuhlmuller B, Grutzkau A, et al. Does gene expression analysis inform us in rheumatoid arthritis? Ann Rheum Dis 2010;69:i37–42. [7] Bild AH, Yao G, Chang JT, et al. Oncogenic pathway signatures in human cancers as a guide to targeted therapies. Nature 2006;439:353–7. [8] Ebert BL, Galili N, Tamayo P, et al. An erythroid differentiation signature predicts response to lenalidomide in myelodysplastic syndrome. PLoS Med 2008;5:e35. [9] Holleman A, Cheok MH, den Boer ML, et al. Gene-expression patterns in drugresistant acute lymphoblastic leukemia cells and response to treatment. N Engl J Med 2004;351:533–42. [10] Pomeroy SL, Tamayo P, Gaasenbeek M, et al. Prediction of central nervous system embryonal tumour outcome based on gene expression. Nature 2002;415:436–42. [11] van’t Veer LJ, Dai H, van de Vijver MJ, et al. Gene expression profiling predicts clinical outcome of breast cancer. Nature 2002;415:530–6. [12] Julia A, Erra A, Palacio C, et al. An eight-gene blood expression profile predicts the response to infliximab in rheumatoid arthritis. PLoS One 2009;4:e7556. [13] Lequerre T, Gauthier-Jauneau AC, Bansard C, et al. Gene profiling in white blood cells predicts infliximab responsiveness in rheumatoid arthritis. Arthritis Res Ther 2006;8:R105. [14] Sekiguchi N, Kawauchi S, Furuya T, et al. Messenger ribonucleic acid expression profile in peripheral blood cells from RA patients following treatment with an anti-TNF-alpha monoclonal antibody, infliximab. Rheumatology (Oxford) 2008;47:780–8. [15] Stuhlmuller B, Haupl T, Hernandez MM, et al. CD11c as a transcriptional biomarker to predict response to anti-TNF monotherapy with adalimumab in patients with rheumatoid arthritis. Clin Pharmacol Ther 2010;87:311–21. [16] Tanino M, Matoba R, Nakamura S, et al. Prediction of efficacy of anti-TNF biologic agent, infliximab, for rheumatoid arthritis patients using a comprehensive transcriptome analysis of white blood cells. Biochem Biophys Res Commun 2009;387:261–5. [17] Toonen EJ, Gilissen C, Franke B, et al. Validation study of existing gene expression signatures for anti-TNF treatment in patients with rheumatoid arthritis. PLoS One 2012;7:e33199. [18] Badot V, Galant C, Nzeusseu Toukap A, et al. Gene expression profiling in the synovium identifies a predictive signature of absence of response to adalimumab therapy in rheumatoid arthritis. Arthritis Res Ther 2009;11:R57. [19] Lindberg J, Wijbrandts CA, van Baarsen LG, et al. The gene expression profile in the synovium as a predictor of the clinical response to infliximab treatment in rheumatoid arthritis. PLoS One 2010;5:e11310. [20] van Baarsen LG, Wijbrandts CA, Rustenburg F, et al. Regulation of IFN response gene activity during infliximab treatment in rheumatoid arthritis is associated with clinical response to treatment. Arthritis Res Ther 2010;12:R11. [21] Barrett T, Wilhite SE, Ledoux P, et al. NCBI GEO: archive for functional genomics data sets – update. Nucleic Acids Res 2012;41:D991–5. [22] Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009;4:44–57. [23] Ramasamy A, Mondry A, Holmes CC, et al. Key issues in conducting a metaanalysis of gene expression microarray datasets. PLoS Med 2008;5:e184. [24] Wang X, Kang DD, Shen K, et al. An R package suite for microarray metaanalysis in quality control, differentially expressed gene analysis and pathway enrichment detection. Bioinformatics 2012;28:2534–6. [25] Kang DD, Sibille E, Kaminski N, et al. MetaQC: objective quality control and inclusion/exclusion criteria for genomic meta-analysis. Nucleic Acids Res 2012;40:e15. [26] Pascual V, Chaussabel D, Banchereau J. A genomic approach to human autoimmune diseases. Annu Rev Immunol 2010;28:535–71. [27] Cahan P, Rovegno F, Mooney D, et al. Meta-analysis of microarray results: challenges, opportunities, and recommendations for standardization. Gene 2007;401:12–8. [28] Tseng GC, Ghosh D, Feingold E. Comprehensive literature review and statistical considerations for microarray meta-analysis. Nucleic Acids Res 2012;40:3785–99. [29] Wang X, Lin Y, Song C, et al. Detecting disease-associated genes with confounding variable adjustment and the impact on genomic meta-analysis: with application to major depressive disorder. BMC Bioinformatics 2012;13:52. [30] Heckmann BL, Zhang X, Xie X, et al. The G0/G1 switch gene 2 (G0S2): regulating metabolism and beyond. Biochim Biophys Acta 2013;1831:276–81. [31] Nakamura N, Shimaoka Y, Tougan T, et al. Isolation and expression profiling of genes upregulated in bone marrow-derived mononuclear cells of rheumatoid arthritis patients. DNA Res 2006;13:169–83. [32] Lindberg J, af Klint E, Catrina AI, et al. Effect of infliximab on mRNA expression profiles in synovial tissue of rheumatoid arthritis patients. Arthritis Res Ther 2006;8:R179. [33] van der Pouw Kraan TC, Wijbrandts CA, van Baarsen LG, et al. Responsiveness to anti-tumour necrosis factor alpha therapy is related to pre-treatment tissue inflammation levels in rheumatoid arthritis patients. Ann Rheum Dis 2008;67:563–6.
Please cite this article in press as: Kim T-H, et al. Gene expression profile predicting the response to anti-TNF treatment in patients with rheumatoid arthritis; analysis of GEO datasets. Joint Bone Spine (2014), doi:10.1016/j.jbspin.2014.01.013
