A two-stage search strategy for detecting multiple loci associated with rheumatoid arthritis.

Gene × gene interactions play important roles in the etiology of complex multi-factorial diseases like rheumatoid arthritis (RA). In this paper, we describe our use of a two-stage search strategy consisting of information theoretic methods and logistic regression to detect gene × gene interactions associated with RA using the data in Problem 1 of Genetic Analysis Workshop 16. Our method detected interactions of several SNPs (single-SNP and SNP × SNP) that are located on chromosomal regions linked to RA and related diseases in previous studies.


Background
The risk of developing many common and complex diseases such as cancer and autoimmune disease involve complex interactions between multiple genes and several endogenous and exogenous environmental factors (or covariates). Rheumatoid arthritis (RA) is a complex genetic disease in which it is hypothesized that several loci contribute to disease susceptibility. Information theoretic methods are among the most promising approaches for genetic association studies and have been used for genetic analysis [1,2] and analysis of gene × gene interactions [3,4]. In this paper, we describe our use of a two-stage strategy consisting of an information theoretic search followed by logistic regression to detect gene × gene interactions associated with RA using selected genomic regions from the genome-wide scan data from the North American Rheumatoid Arthritis Consortium, which comprises 868 cases and 1194 controls. Data were provided as Problem 1 of Genetic Analysis Workshop 16.

Methods
Interaction information as measure of association Let X i denotes a genetic random variable representing the genotypes at locus L i . We assume L i is biallelic (with alleles A and a) with three possible genotypes (AA, Aa, and aa). The uncertainty of X i is given by Shannon's entropy [5] as

Open Access
Given a set of such genetic variables S = {X 1 ; X 2 ;...; X k }, the interaction information among the k variables (referred to as k-way interaction information or KWII) is defined as the amount of information (redundancy or synergy) present in the set of variables that is not present in any subset of these variables [4]. For the variables in set S, the KWII can be written succinctly as an alternating sum over entropies (H) of all possible subsets τ of S using the difference operator [6]: Let C be the random variable representing the disease status (phenotype variable) of RA. Then KWII(S;C) = KWII (X 1 ; X 2 ;...; X k ;C) is a measure of the association of the set of genetic variables in set S towards the disease phenotype variable C (i.e., how well the set explains the disease phenotype). The value of KWII(S;C) can be both positive and negative. We shall use only positive KWII values as the measure of association because larger positive values indicate stronger interaction (hence, higher association).

Redundancy between combinations of variables
Let S 1 = {X 1 ; ...; X m } and S 2 = {Y 1 ; ...; Y m } be two sets (or combinations) of variables. Then the redundancy between S 1 and S 2 is given by the maximized average of pairwise linkage disequilibrium (LD) (r 2 ) between variables from S 1 and S 2 : Such redundancies can arise because of LD between the variables across each set. For example, for a disease C that is caused by interactions between two untyped SNPs D 1 and D 2 , let four marker loci be designated X 1 , X 2 , X 3 , and X 4 such that X 1 and X 3 are in strong LD with D 1 , while X 2 and X 4 are in strong LD with D 2 . Then the KWII(X 1 ;X 2 ;C) and KWII(X 3 ;X 4 ;C) measure the association of the sets {X 1 ;X 2 } and {X 3 ;X 4 } for C, respectively. The redundancy between the combinations {X 1 ;X 2 } and {X 3 ;X 4 } is given by say, 0.5*(r 2 (X 1 , X 3 )+r 2 (X 2 , X 4 )) and existence of strong LD between X 1 and X 3 and between X 2 and X 4 will result in similar measures of KWII association for both sets, making one of the sets statistically redundant.
Stage I: Single-nucleotide polymorphism (SNP)-combination search strategy Let S be the set of all genetic (SNPs) and environmental (non-genetic) variables (e.g., sex) and C be the variable denoting the disease phenotype. The information theoretic metric KWII(X 1 ;...;X k ;C) is a measure of the association of the set of variables with the disease phenotype variable C (i.e., how well they explain the disease phenotype). Using this metric and a redundancy measure, we iteratively search for combinations of variables up to a fixed number (say τ) of iterations. Let the number of variables (except C) in a combination be defined as the "order" of the combination. In our method, we limit our search to up to second-order (or two-variable) combinations (i.e., we consider only {X i ;C} and {X i ;X j ;C} combinations). Let θ be the set of variables and ξ be the set of associated combinations output by our search method. Initially, both θ and ξ are empty. In iteration = 1, the variable X k having highest Also X k is removed from S. In a subsequent iteration = i (i > 1), a new variable X j S is considered for selection and its single variable and two-variable combinations are formed and KWII computed (using, Eq. (2)) with variables already selected in the previous iterations. At the same time each of the combinations formed are checked for redundancy with combinations already in ξ and of same order (using Eq. (3) and redundancy exceeding a threshold of 0.7). For example in iteration = 2, for X j S, the combinations {X j ;C} and {X k ;X j ;C} are formed and {X j ;C} is checked for redundancy with {X k ;C}. From all the new variables, the variable that has maximum KWII and all nonredundant combinations is selected. A variable with a redundant combination is dropped from consideration (i.e., removed from S) in subsequent iterations. Given the computational burden of determining redundancy with combinations of variables already selected, our selection procedure stops after a maximum of τ = 50 iterations. Thus, up to 50 variables with non-redundant combinations and highest KWII are selected. This stage yields a number of single and two-variable combinations and their KWII values, which are input to the second stage.

Stage II
We conduct logistic regression analysis on the one-and two-variable combinations obtained by our information theoretic search using the methods outlined by Cordell [7] and the significance of each combination is determined. The full single-locus model is where r is the probability of each individual being a case, μ corresponds to the mean effect, the terms a and d correspond to the additive and dominance coefficient effects of the tested SNP variable, x and z are dummy variables with x = 1, z = -0.5 for one homozygote genotype (AA), x = 0, z = 0.5 for the heterozygote genotypes (Aa), and x = -1, z = -0.5 for the other homozygote type (aa). The chi-square is used to compare the full single-locus model with the null model given by 0 values for both a and d.
For SNP × SNP interactions, fully saturated model following Cordell's notation [7] is where r and μ are same as in Eq. (3), the terms a 1 , d 1 , a 2 , and d 2 are the dominance and additive effect coefficients of the two SNPs, i aa , i ad , i da , and i dd represent their interaction coefficients, and x i and z i are dummy variables with x i = 1, z i = -0.5 for one homozygous genotype (AA or BB), x i = 0, z i = 0.5 for the heterozygous genotypes (Aa or Bb), and x i = -1, z i = -0.5 for the other homozygous genotype (aa or bb). An interaction is tested by the deviance of the full two-locus model from the model minus the interaction terms with chi-square test.

Data
We have followed a candidate-gene-based approach and selected SNPs belonging to the candidate genes/regions in Table 1 for exploring both gene × RA and gene × gene × RA interactions using our two-stage approach. The start and end base-pair positions of each gene are obtained from http://www.pharmgkb.org/. Using the genes/regions from Table 1, we created the following three data sets for analysis: • 7087 SNPs selected for analysis using all genes/ regions (Data Set 1) • 5385 SNPs selected using all genes/regions except those that belong to only 6p21.3 and not to any other gene (Data Set 2) • 3263 SNPs selected using genes not on chromosome 6 (Data Set 3) Additionally, sex of the subjects and RA status were present in each data set as the environmental variable and the phenotype variable (C).

Results
We have obtained many single-variable and two-variable interactions with the disease phenotype, only the combinations with high values of KWII are presented in Tables 2, 3, 4. The SNPs shown to be in genomic regions 6p21.3 and 6q23 do not overlap with any other gene. We found no interaction between the covariates sex and RA. Table 2 shows the single-variable combinations with KWII values greater than or equal to 95 th percentile of all the single-variable KWII obtained using our method for the respective data sets. Tables 3 and 4 show the two-variable combinations with KWII values greater than or equal to 95 th percentile of all the two variable KWII obtained using our method for the respective data sets.      6q23 and HLA-C (with KWII95 1 overall <KWII<95 th percentile for respective data sets, not shown in tables). The strongest of the two-variable interactions are from the SNPs in region 6p21.3 (Table 3) obtained using Data Set 1. We have created the other two data sets because it was felt that several relatively weaker interactions are difficult to detect in the presence of the strongest interactions in 6p21.3. Using Data Sets 2 and 3, we found several two-variable KWII in genes MICA and RUNX1 (Table 3). Also several two-variable interactions are detected among SNPs in HLA-C, HLA-DR, MICA, and 6p21.3 (with KWII95 2 overall <KWII<95 th percentile for respective data sets, not shown in tables). Also we observed an interaction between rs11811771 (PTPN22) with rs2828104 (RUNX1) with p-value 2.7 × 10 -5 and KWII = 0.095. Separately we also calculated KWII for two-SNP combinations for Data Set 1 wherein the SNPs belong to different genes/genomic regions ( Table 4).

Discussion
We have used a two-stage strategy to search for single SNPs and SNP × SNP interactions associated with RA. Using our analysis on the candidate genes, we have found several strong interactions on 6p21.3 and interactions among SNPs on genes previously reported to be related with RA and other autoimmune diseases. For example, RUNX1 has been reported to be associated with systemic lupus erythematosus and psoriasis (two autoimmune diseases) [8,9] while associations of region 6q23 and MICA with RA has been reported by Thomson et al. [10] and Martinez et al. [11], respectively. Detecting genes and environmental factors interacting to increase the susceptibility to disease risk is a very challenging task for many reasons, particularly for the large size of the data and presence of confounding factors such as LD, presence of phenocopies, locus heterogeneity, and population stratification. Information theoretic methods have high power in detecting gene × gene interactions and have the advantage of being simpler and computationally faster; KWII-based interaction analysis has been employed in [3,4]. Also, our method can be used when the genetic and environmental variables have different numbers of classes or when the phenotype has more than two classes. Although we initially planned for a genome-wide analysis, given the large size of the data, we were able to execute only a few iterations using our computational resources. Therefore, we decided to follow a candidate-gene-based approach. We believe that with the help of additional hardware, it is possible to implement our search strategy in a distributed computing environment employing multiple processors and to explore many more interactions with moderate to low magnitudes that are potentially associated with RA. the development of the computational analysis. LS participated in the statistical genetics analysis and interpretations. MR conceived the study and was involved all aspects of design and coordination.