Simultaneous genome-wide association studies of anti-cyclic citrullinated peptide in rheumatoid arthritis using penalized orthogonal-components regression.

Genome-wide associations between single-nucleotide polymorphisms and clinical traits were simultaneously conducted using penalized orthogonal-components regression. This method was developed to identify the genetic variants controlling phenotypes from a massive number of candidate variants. By investigating the association between all single-nucleotide polymorphisms to the phenotype of antibodies against cyclic citrullinated peptide using the rheumatoid arthritis data provided by Genetic Analysis Workshop 16, we identified genetic regions which may contribute to the pathogenesis of rheumatoid arthritis. Bioinformatic analysis of these genomic regions showed most of them harbor protein-coding gene(s).


Background
Most of the available statistical approaches for genomewide association study (GWAS) have focused on studying one single-nucleotide polymorphism (SNP) at a time [1], thereby ignoring the multigenic nature of complex diseases [2,3] and the strong correlation between some SNPs due to linkage disquilibrium (i.e., some SNPs are inherited together in blocks of DNA). As pointed out by Waldron et al. [4], association studies using multiple SNPs have substantial advantages over those based on SNP associations. To capture the correlation between SNPs in regions of low recombination, haplotype-based methods that recognize the existence of linkage disequilibrium between genetic markers have been employed for multi-SNP analyses. However, such methods introduce additional problems, including the need to infer haplotypes, the impact of including (or excluding) rare haplotypes, and the need to define haplotype block boundaries. Even with haplotype-based association methods, high false-positive rates and a lack of reproducibility remain major concerns.

Open Access
Because hopes remain high for the value of GWAS, other approaches need to be pursued that account for the correlation structure among SNPs. Here we have developed a method for GWAS that incorporates regression of a phenotype simultaneously on all available SNPs, i.e., considering a multiple linear regression, where Y is the phenotypic value, X j counts one of the two alleles at the j th SNP, and b j is the additive effect of that allele, j = 1,..., p. To account for the issues raised by large p, such as lack of independence between SNPs due to linkage disequilibrium, we will conduct the GWAS using the penalized orthogonal-components regression (POCRE) [5]. POCRE offers a fast and efficient way to identify significant SNPs simultaneously from a large number of candidates.

Methods
GWAS usually entails the collection of a massive amount of SNPs (i.e., large p) for only a small number of biological individuals (i.e., small n). Therefore, identifying the few genetic variants underlying disease risk is equivalent to the task of finding "a very few needles in a haystack", and poses a challenging statistical issue. Zhang et al. [5] recently described the POCRE approach, which sequentially constructs sparsely loaded orthogonal components with proper regularization. They demonstrate this approach works well when fitting regression models with n <<p data.
Let Y be an n-dimensional column phenotype vector, and X be an n × p genotype matrix, where n and p are the number of individuals and number of SNPs, respectively. Further assume both Y and X are centralized, and accordingly assume μ = 0 in Eq. (1). POCRE sequentially constructs orthogonal components X X where g l (g) is a penalty function defined by a proper regularization on g with tuning parameter l.
When g l (g) ≡ 0, the optimal g solving Eq. (2) is proportional to the leading eigenvector of X YY X k T T k . Zhang et al. [5] employed empirical Bayes thresholding methods proposed by Johnstone and Silverman [6] to introduce proper penalty g l (g). Such penalty benefits estimating covariance between phenotype and genotypes, and provides adaptively sparse loadings of orthogonal components. The empirical Bayes implementation is also computationally efficient. The tuning parameter l also accounts for possible dependence structure among different SNPs, and 10-fold crossvalidation was employed to elicit its optimal values ranging from 0. 8  The sequential construction of the orthogonal components stops when the optimal g solving Eq. (2) is zero, which implies X is hardly correlated to Y Then, regressing Y on the orthogonal components, i.e., X X , , provides an estimate of b 1 ,..., b p in Eq. (1). Because nonzero loadings in ω j , j = 1,2,..., are sparse, most of estimated b 1 ,..., b p are therefore zero, reflecting the fact that most SNPs are insignificantly associated or are even completely uncorrelated to the phenotype of interest.

Results
Using our novel method (POCRE), the rheumatoid arthritis data from Genetic Analysis Workshop 16 were investigated for associations between SNPs and a serum biomarker for rheumatoid arthritis, i.e., antibodies anticyclic citrullinated peptide (anti-CCP). In this dataset, only 867 samples were positive for anti-CCP. The data set was preprocessed with the computer program PLINK [7] to control data quality, and with the computer program EIGENSTRAT [8] to control potential population structures. After preprocessing there were 490,613 SNPs remaining for the GWAS. POCRE was applied individually to each chromosome for the simultaneous association of the SNPs in that chromosome with the anti-CCP phenotype. The effects of the 10 principal components constructed by EIGENSTRAT were considered as covariates for the POCRE. Nonzero effects of SNPs were reported on seven chromosomes where the positiveβ indicates that the minor allele will increase the level of anti-CCP (Table 1).
Of the 12 SNPs identified to be associated with anti-CCP level, 5 SNPs are from chromosome 6, 2 SNPs from chromosome 12, and 1 SNP each is from chromosomes 7, 11, 13, 17, and 21, respectively. The location of each significant SNP was mapped to the human genome using the Ensembl database http://www.ensembl.org. Based on this analysis, eight of the SNPs were found to reside in seven genes (Table 1). For the other four SNPs, the nearest neighboring genes are listed in the table. None of these genes we identified have previously been linked to rheumatoid arthritis, but several of them encode proteins whose biological roles may be involved in the BMC Proceedings 2009, 3(Suppl 7):S20 http://www.biomedcentral.com/1753-6561/3/S7/S20 pathogenesis of this disease. For example, two SNPs (rs922898, rs6929401) that reside in the uronyl-2sulfotransferase (UST) gene were identified by our method. UST is involved in the chondroitin 3 sulfate and glycan structure biosynthesis pathways and could contribute to optimal cartilage development or repair.
PPIL4 is a member of cyclophilin family, a group of proteins crucial for protein folding and immunosuppression by cyclosporin A (CsA) [9]. Polymorphisms in this gene may influence development of inflammation during rheumatoid arthritis or the response of individuals with rheumatoid arthritis to treatment. DGKB encodes diacylglycerol kinase and participates in intracellular signalling processes via several pathways including the protein kinase C pathway [10]. Alteration in intracellular signalling could negatively influence inflammatory processes. Finally, CHODL encodes a type I transmembrane protein including a single carbohydrate recognition domain for the C-type lectins that can worsen inflammation [11].

Discussion
In addition to its diagnostic value for the general risk of rheumatoid arthritis, a high anti-CCP level is associated with a high risk of developing joint damage during the disease. As a result, an anti-CCP-test has been employed to monitor the progression of the disease. With our novel approach for GWAS, we identified several candidate SNPs associated with the level of anti-CCP in rheumatoid arthritis patients. Among the candidate genes identified, only PPIL4 is functionally related to the immune system. Further investigation will be necessary to define the potential roles that the other gene products play in rheumatoid arthritis. Due to the small sample size available for this study (867) and the large total number of SNPs measured (490,613), we expect some of the associations we identified to constitute false positives. Determining false detection rates will require development of additional procedures.

Conclusion
Our analyses using the newly developed method POCRE indicate that the genomic region 6q25.1 may harbor genes associated with anti-CCP level in rheumatoid arthritis patients. Further investigation is necessary to confirm this observation.