Genome-wide association studies and gene expression profiles of rheumatoid arthritis: An analysis.

OBJECTIVES: The molecular mechanism of rheumatoid arthritis (RA) remains elusive. We conducted a protein-protein interaction network-based integrative analysis of genome-wide association studies (GWAS) and gene expression profiles of RA.
METHODS: We first performed a dense search of RA-associated gene modules by integrating a large GWAS meta-analysis dataset (containing 5539 RA patients and 20 169 healthy controls), protein interaction network and gene expression profiles of RA synovium and peripheral blood mononuclear cells (PBMCs). Gene ontology (GO) enrichment analysis was conducted by DAVID. The protein association networks of gene modules were generated by STRING.
RESULTS: For RA synovium, the top-ranked gene module is HLA-A, containing TAP2, HLA-A, HLA-C, TAPBP and LILRB1 genes. For RA PBMCs, the top-ranked gene module is GRB7, consisting of HLA-DRB5, HLA-DRA, GRB7, CD63 and KIT genes. Functional enrichment analysis identified three significant GO terms for RA synovium, including antigen processing and presentation of peptide antigen via major histocompatibility complex class I (false discovery rate (FDR) = 4.86 × 10 - 4), antigen processing and presentation of peptide antigen (FDR = 2.33 × 10 - 3) and eukaryotic translation initiation factor 4F complex (FDR = 2.52 × 10 - 2).
CONCLUSION: This study reported several RA-associated gene modules and their functional association networks.Cite this article: X. Xiao, J. Hao, Y. Wen, W. Wang, X. Guo, F. Zhang. Genome-wide association studies and gene expression profiles of rheumatoid arthritis: an analysis. Bone Joint Res 2016;5:314-319. DOI: 10.1302/2046-3758.57.2000502.

The molecular mechanism of rheumatoid arthritis (RA) remains elusive. Integrating the information of genome-wide association studies (GWAS), microarray and protein interaction networks can provide a novel insight into the mechanism of RA.

dmGWAS analysis identified several RA-associated gene modules, mainly including HLA-A and GRB7.

Enrichment analysis detected three significant gene ontology terms for RA functionally involved in antigen processing and presentation.

Our study results may help to reveal the pathogenesis of RA. Further studies are needed to confirm our findings.

Introduction

Rheumatoid arthritis (RA) is one of the most common systemic autoimmune chronic inflammatory disorders. It is characterised by progressive joint destruction and multisystem inflammation, causing joint pain, swelling, stiffness and deformities. RA affects more than 1% of the population worldwide,[1] and results in a reduction of lifespan of RA patients of about ten years.[2] Due to its high morbidity and mortality, RA has a significant societal impact in terms of medical cost and lost productivity.

Benefiting from the rapid development of high-throughput omics technologies, researchers are able to monitor hundreds of thousands of genes and proteins involved in the development of RA. For instance, in the last three years more than ten genome-wide association studies (GWAS) have been conducted and have reported a set of RA-associated genes.[3] But due to the analysing strategy of focusing on the most significant loci, the susceptibility genes identified by GWAS are mostly limited and explain a small proportion of the genetic risk of RA. At the RNA level, extensive gene expression profile studies have been performed and have identified some important differentially expressed genes implicated in the development of RA.[4,5] However, the rich pathogenetic information contained in the abundant omics data was not fully used, limiting our efforts to clarify the elusive molecular mechanism and develop effective treatments for RA.

As the pathogenesis of RA is complicated, the traditional single-omic study approach is usually insufficient for fully unraveling the functional principles and dynamics of the global biological process of complex diseases. Multi-omics integrative studies have gained more and more attention in pathogenetic studies.[6-10] These studies are expected to provide more pathogenetic information and could potentially reduce the impact of data error and statistical bias on study results. Recently, the protein-protein interaction (PPI) network was introduced into omics data analysis, which could also enhance the identification of disease genes and help us to understand their relationships.[6,10,11]

In this study, we conducted a PPI network-based integrative analysis of GWAS and gene expression profiles of RA. Briefly, we first performed a dense search of RA-associated gene modules by integrative analysis of a large GWAS meta-analysis dataset (containing 5539 RA patients and 20 169 healthy controls) and gene expression profiles of RA synovium and peripheral blood mononuclear cells (PBMC). The biological function of the identified gene modules was then investigated by gene ontology (GO) enrichment analysis and protein-association network analysis. Our study identified several RA-associated gene modules and their functional association networks in the pathogenesis of RA. We also illustrated the performance of multi-omics data integration and mining for RA pathogenetic studies.

Materials and Methods

GWAS datasets

We used a large GWAS meta-analysis dataset involving 5539 autoantibody-positive RA patients and 20 169 healthy controls of European ancestry.[12] This large sample consisted of: 483 RA cases and 1449 healthy controls (genotyped by Affymetrix 6.0; Affymetrix, Santa Clara, California) from the Brigham Rheumatoid Arthritis Sequential Study; 589 RA cases and 1472 healthy controls (genotyped by Illumina 370K; Illumina, San Diego, California) from the CANADA RA study; 1173 RA cases and 1089 healthy controls (genotyped by Illumina 317K; Illumina) from the Epidemiological Investigation of Rheumatoid Arthritis; 1769 RA cases and 5551 healthy controls (genotyped by Illumina 317K and Illumina 370K; Illumina) from the North American Rheumatoid Arthritis Consortium I and III; 1525 RA cases and 10 608 healthy controls (genotyped by Affymetrix 500K;Affymetrix) from the Wellcome Trust Case Control Consortium (WTCCC). RA was diagnosed according to the 1987 American College of Rheumatology criteria for diagnosis of RA or by board-certified experts of rheumatology.

All study samples passed a genetic relatedness check based on pairwise identity by state testing.[12] Genetically related subjects were removed. The single nucleotide polymorphisms (SNPs) with call rates < 95.00% (< 99% for WTCCC GWAS data), minor allele frequencies < 0.01 or Hardy-Weinberg Equilibrium testing p-values < 10−6 (< 10−5 for WTCCC GWAS data) were excluded. EIGENSTRAT analysis was performed and did not detect significant impact of population stratification on the GWAS results (genomic control inflation factor λ = 1.04).[13] Ungenotyped SNPs across the genome were imputed by IMPUTE software (MathWorks, Natick, Massachusetts) against the reference panel of HapMap Phase 2 European Central European University founders.[14] After quality control, a total of 2 554 714 SNPs with high quality genotype data were used. Meta-analysis was performed by summing inverse variance-weighted β and z cores. The final genomic control corrected meta-analysis results were reported. A detailed description of study subjects, experimental procedures and analysis approaches can be found in the original study.[12]

Gene expression profile of RA synovium

The gene expression profile data of RA synovium were downloaded from the Gene Expression Omnibus (GEO, access number: GSE1919) database.[15] Briefly, synovial tissue samples were collected from eight RA patients undergoing synovectomy or arthroplasty, and from 15 normal subjects who died from fatal accidents, and were matched for age and gender. All donors were Caucasians living in Berlin, Germany. Total ribonucleic acid (RNA) was isolated from the synovium using RNeasy spin columns (Qiagen, Hilden, Germany). The isolated total RNA was labelled and hybridised to Affymetrix Human Genome U95A Array (containing 12 600 oligonucleotide probes) following standard protocol. Global normalisation of raw signal intensities was performed by GeneChip software (MAS 5.0, Affymetrix). Data analysis was conducted using DMT 3.0 software (Affymetrix). A detailed description of study subjects, experimental procedures and analysis approaches can be found in the previous study.[4]

Gene expression profile of RA PBMC

The gene expression profile of PBMC of RA was downloaded from the GEO database (access number: GSE15573). Briefly, peripheral blood specimens were collected from 18 RA patients and 15 healthy controls, matched for age and gender. All were Caucasians living in Paris, France. Total RNA was isolated from PBMC using the PAXgene RNA isolation kit (PreAnalytix, Hombrechtikon, Switzerland). The isolated total RNA was transcribed into cRNA and labelled biotin-16-UTP. Illumina human-6 v2.0 expression beadchip (Illumina, San Diego, California) (containing 12 600 oligonucleotide probes) was applied for microarray hybridisation following standard protocol. Hybridisation signals were scanned by Illumina BeadArray 500GX Reader and analysed by Illumina BeadScan software 2.3.0.13 (Illumina). BeadStudio software (version 1.5.0.34) was used for preliminary data analysis and quality control. The generated microarray data were expressed as log2 ratios of the fluorescence intensities of RA samples divided by the fluorescence intensities of control samples. Normalisation was conducted using the "normalise quantiles" function in the BeadStudio software. A detailed description of study subjects, experimental procedures and analysis approaches can be found in the original study.[5]

Dense gene module search

The integrative analysis of RA GWAS and gene expression profiles data was performed by dmGWAS 3.0 (Bioinformatics, Houston, Texas).[16,17] dmGWAS implemented a PPI network-based Dense Module Search (DMS) algorithm.[16,17] The detailed description of the dmGWAS algorithm can be found in previous studies[16,17] The brief analysis procedures of dmGWAS can be summarised as follows: the node weights of PPI network are derived from GWAS results. The edge weights of PPI network are derived from gene expression profiling results. A greedy algorithm was applied to the PPI network for searching dense modules enriched in the gene sets showing association evidence and abnormal expression for target diseases. dmGWAS has been successfully applied to the causal gene studies of multiple complex diseases.[6,10,18] The PPI network used by Goh et al[19] was applied here.[6] This PPI network consisted of 10 174 proteins and 61 070 protein-protein interactions.[20,21] In this study, dmGWAS was running under default parameters recommended by the developers. During the processing of RA GWAS data, a SNP was assigned to a gene if it was located within the gene or up to 20 kb immediately upstream or downstream of the gene. We assigned each gene a gene-level p-value, which was the smallest p-value of the SNPs assigned to the gene. The candidate genes residing in the top three identified gene modules were subjected to further function evaluation.

Protein-association network

To understand the functional association network of the identified gene modules, the genes within the top three gene modules were input into STRING (version 9.1)[22] for protein-protein interaction analysis. STRING is a widely used protein-association network database[23,24] which can help to assemble and evaluate the relationships of disease-related proteins.

Functional enrichment analysis

The biological functional significance of identified gene modules was evaluated by functional enrichment analysis. The genes within the top three gene modules were input into the DAVID tool[25] for GO enrichment analysis. An enrichment p-value and a false discovery rate (FDR) were calculated by DAVID for each GO term. Significant GO terms were identified as FDR < 0.05.

Results

For RA synovium, 1149 gene modules were detected by dmGWAS analysis. Table I summarises the top three gene modules consisting of 16 genes. The top-ranked gene module is named HLA-A, which contains TAP2, HLA-A, HLA-C, TAPBP and LILRB1 genes. Within the HLA-A module, TAP2 achieved the most significant p-value (2.46 × 10−126) in the GWAS meta-analysis of RA.[12] The gene expression study also found that TAP2 was upregulated in RA synovium (ratio = 2.10 sd 1.66).[4] HLA-C (GWAS p-value = 1.57 × 10−40, ratio = 19.96 sd 19.96) is another notable gene in the HLA-A modules. The protein-association network analysis result of the top three gene modules of RA synovium is shown in Figure 1a.

Table I.

The top three gene modules identified for rheumatoid arthritis synovium

Module name	Gene	GWAS meta-analysis p-value
HLA-A	TAP2	2.46E−126
	HLA-C	1.57E−40
	HLA-A	1.37E−24
	TAPBP	7.63E−11
	LILRB1	2.00E−02
PTPN2	DAXX	9.19E−20
	PTPN2	4.33E−05
	EIF4E	2.12E−03
	EIF4G2	3.42E−03
	VCL	1.00E−02
	STAT3	9.00E−02
	EIF4B	2.00E−02
	PML	1.60E−01
STAT3	VCL	1.00E−02
	STAT3	9.00E−02
	DAXX	9.19E−20
	PML	1.60E−01
	EIF4B	2.00E−02
	EIF4E	2.12E−03
	EIF4G2	3.42E−03

GWAS, genome-wide association study

Images showing protein-association networks of gene modules identified for rheumatoid arthritis synovium (a) and peripheral blood mononuclear cells (b). The lines with different colours represent different types of association evidence. The red triangles denote the key genes of the top three significant gene modules.

In total, 1675 gene modules were detected for RA PBMC. Table II summarises the top three gene modules, which contain 26 genes. The top-ranked gene module is named GRB7, consisting of CD63, HLA-DRA, KIT, GRB7 and HLA-DRB5 genes. Within the GRB7 module, HLA-DRB5 achieved the most significant association signal (GWAS p-value = 3.51 × 10−235),[12] and was also upregulated (ratio = 4.53 sd 4.83) in RA PBMC.[5] Figure 1b presents the protein-association network analysis result of the top three gene modules of PBMC.

Table II.

The top three gene modules identified for rheumatoid arthritis peripheral blood mononuclear cells

Module name	Gene	GWAS meta-analysis p-value
GRB7	HLA-DRB5	3.51E−235
	HLA-DRA	1.00E−285
	KIT	2.65E−03
	GRB7	7.70E−01
	CD63	7.00E−02
MAP4K5	HLA-DQA2	4.08E−285
	HLA-DQB1	4.08E−285
	PIK3R1	1.72E−03
	CRKL	1.30E−01
	CD4	6.00E−02
	CD3E	5.00E−02
	MAP4K5	7.00E−02
	VAV3	1.00E−02
	IL1R1	5.05E−03
ZFP36L1	PTPN22	7.36E−27
	EGFR	1.94E−04
	NRG1	2.61E−03
	ITGA5	3.31E−03
	ZFP36L1	4.45E−03
	PTPRE	1.00E−02
	GRB2	2.00E−02
	MAPK14	3.00E−02

GWAS, genome-wide association study

The functional enrichment analysis results of identified gene modules are presented in Table III. For RA synovium, we identified three significant GO terms, including antigen processing and presentation of peptide antigen via major histocompatibility complex (MHC) class I (FDR = 4.86 × 10−4), antigen processing and presentation of peptide antigen (FDR = 2.33 × 10−3) and eukaryotic translation initiation factor 4F complex (FDR = 2.52 × 10−2). For RA PBMC, seven significant GO terms were identified, including transmembrane receptor protein tyrosine kinase signalling pathway (FDR = 4.34 × 10−5), MHC class II protein complex (FDR = 7.87 × 10−3), plasma membrane (FDR = 1.95 × 10−2), protein kinase cascade (FDR = 2.25 × 10−2), antigen processing and presentation of peptide or polysaccharide antigen via MHC class II (FDR = 2.50 × 10−2), positive regulation of kinase activity (FDR = 3.31 × 10−2) and positive regulation of transferase activity (FDR = 3.98 × 10−2). Interestingly, we found that antigen processing and presentation-related GO terms were enriched in the top three gene modules of both RA synovium and PBMC.

Table III.

Gene ontology enrichment analysis results of identified gene modules

	Gene ontologyname	Gene ontology number	Enrichment score[*]	p-value	False discovery rate
Synovium	Antigen processing and presentation of peptide antigen via MHC class I	0002474	244.85	3.60E−07	4.86E-04
	Antigen processing and presentation of peptide antigen	0048002	148.66	1.73E−06	2.33E-03
	Eukaryotic translation initiation factor 4F complex	0016281	355.06	2.42E−05	2.52E-02
Peripheral blood	Transmembrane receptor protein tyrosine kinase signaling pathway	GO:0007169	21.96	2.97E−08	4.34E-05
	MHC class II protein complex	GO:0042613	97.95	6.99E−06	7.87E-03
	Plasma membrane part	GO:0044459	3.87	1.74E−05	1.95E-02
	Protein kinase cascade	GO:0007243	11.63	1.54E−05	2.25E-02
	Antigen processing and presentation of peptide or polysaccharide antigen via MHC class II	GO:0002504	74.53	1.71E−05	2.50E-02
	Positive regulation of kinase activity	GO:0033674	15.97	2.26E−05	3.31E-02
	Positive regulation of transferase activity	GO:0051347	15.37	2.72E−05	3.98E-02

Enrichment score calculated by the DAVID tool

MHC, major histocompatibility complex

Discussion

In this study, we performed a PPI network-based integrative analysis of GWAS and gene expression profiles of RA. We identified several RA-associated genes, gene modules and their association networks. Further functional enrichment analysis of identified gene modules highlighted the importance of antigen processing and presentation in the pathogenesis of RA.

It is well documented that genetic factors contribute greatly to the risk of RA.[26] The identified gene modules contain several reported candidate genes of RA, such as PTPN22,[27] CD4,[28] HLA-C,[29] MAPK14[30] and TAP2.[29] Integrating the GWAS and gene expression profile data confirmed the relevance of several reported susceptibility genes to RA. However, the molecular mechanism of multiple RA relevant genes interactively leading to RA remains elusive. The identified gene modules and their functional association networks may help to answer this question. For instance, the protein-association network analysis of RA-associated gene modules suggests that PTPN22, CD4 and MAPK14 indirectly interacted with each other to contribute to the development of RA (Fig. 1b). General functional interactions were also observed among HLA, TAP2 and LILRB1 genes (Fig. 1a). Furthermore, we provide functional association networks of RA-associated gene modules, which may help to clarify the complex pathogenesis of RA. Further studies are required to confirm our findings.

Another finding of this study is the antigen processing and presentation-related GO terms enriched in the identified gene modules of both RA synovium and PBMCs. Antigen processing and presentation is an immunological process, in which antigen-presenting cells prepare antigen and present generated peptide–MHC class II complexes to CD4 + T cells for activation of the immune system. The mechanism of antigen processing and presentation implicated in the imbalance of the immune system of RA remains elusive. Our study results confirm the importance of antigen processing and presentation in the pathogenesis of RA. We also provide some potential key gene modules involved in the imbalance of the immune system of RA.

This study also illustrates the application of multi-omics integrative analysis in the pathogenetic studies of complex diseases. Multi-omics integrative studies are expected to provide more pathogenetic information. Through the mutual complement and validation of various omics data, multi-omics integrative studies have the potential to reduce the impact of data error and statistical bias on study results. Compared with single-omic studies, multi-omics integrative studies have significant advantages for revealing the functional principles and dynamics of the global pathogenetic process of complex diseases. However, the availability of multi-omics integrative study approaches and tools is limited and they warrant further study.

We performed a PPI network-based integrative analysis of GWAS and gene expression profiles of RA. We identified several RA-associated gene modules and their functional association networks in the development of RA. This study provides novel clues for understanding the pathogenesis of RA. We also demonstrate the performance of multi-omics integrative analysis for mechanism studies of complex diseases.