Weighted composite likelihood ratio test using kinship coefficients

Here, we briefly introduce the weighted composite likelihood ratio (WCLR) test developed by Browning et al. [3]. Denote the individual i's genotype at a single marker as $g_i=(g_{i_1},g_{i_2})$ for i = 1, 2,...,n, where $g_{i_1}$ and $g_{i_2}$ are one of the alleles A₁,...,A _l with the corresponding frequency $p_{A_1},\mathrm{...},p_{A_l}$. The weight of each individual is denoted as w _i . A weighted composite likelihood of a single marker is then $L(p)={\displaystyle {\prod }_{i=1}^n{(p_{g_{i_1}}{p}_{g_{i_2}})}^{w_i}}$, where p = ($p_{A_1},\mathrm{...},p_{A_l}$). Based on this composite likelihood, the allele frequency can be estimated for only cases ($\widehat{p}$), only controls ($\widehat{q}$), and all cases and controls ($\widehat{r}$). A likelihood ratio test can be constructed as $2\ln (L(\widehat{p},\widehat{q})/L(\widehat{r}))$ and it asymptotically follows a chi-square distribution with l - 1 degrees of freedom. Similarly, the likelihood ratio test based on haplotypes and a given set of weights can also be constructed. In general, haplotypes for samples are unknown. In this situation, the haplotype frequencies can be estimated via the expectation-maximization (EM) algorithm and incorporated into the test.

The weight w _i of each individual is calculated based on their kinship coefficients. Specifically, w = (w₁,...,w _n ) satisfies 2(w₁,...,w _n )K = (1,...,1), where K is the n × n kinship matrix. The sum of weights, w₁ + ... + w _n , can be considered as the effective sample size, which is greater than the number of families but less than the total number of samples. Thus, we expect the power of this method to be greater than that of the method using one randomly selected case per family. Finally, it is worth noting that w _i = 1 for unrelated individuals. In this situation, the weighted composite likelihood becomes the ordinary likelihood.

Candidate gene analysis

We analyzed two candidate gene data sets from NARAC using CCREL software that implements the WCLR test of Browning et al. [3]. The first data set contains genotypes of 1256 cases from 665 nuclear families and 1519 unrelated controls at 14 SNP markers in the candidate gene PTPN22 [5]. The 14 SNPs are labeled as SNP1 to SNP14 corresponding to their positional order along chromosome 1. One family was excluded from the analysis because of genotyping errors. Most of the parents' genotypes were not available. The second data set contains genotypes of 816 cases from 461 nuclear families and 855 unrelated controls at 20 SNP markers in 14 candidate genes, including PTPN22, CTLA4, TNFRSF1B, PADI4, HAVCR1, IBD5, SLC22A4, IL3, IL4, SUMO4, ILG5, CARD15, RUNX1, and MFL [6]. We performed the case-control association analysis using two methods. In the first method we used all cases and controls and applied the WCLR test. In the second method, we randomly selected one case from each family and included all unrelated controls and applied the allelic chi-square test. For the first data set, we performed the single marker analysis, the multiple marker analysis (the stepwise logistic regression), and the haplotype analysis on two or three markers. We also obtained the linkage disequilibrium (LD) measure (r²) between SNPs from this data set. For the second data set, we performed the single-marker analysis because this data set only includes a few SNPs for each candidate gene.

Power comparison

The power of the two methods based on all related cases and all controls, the weighted composite likelihood ratio test [3] and the method suggested by Slager and Schaid [1], and the power of the method using only one randomly selected case per family and all controls were evaluated and compared based on simulated data. One thousand data sets were generated and each of them consisted of 200 affected sib pairs (400 cases) and 200 unrelated controls. Only genotypes at the disease locus were simulated and analyzed. The minor allele frequencies at the disease locus for cases and controls were set as 0.141 and 0.095 to calculate the power, and set as 0.10 and 0.10 to calculate the type I error rate. The significance level was set as 0.05.

Association of RA status with PTPN22

Using the data of one randomly selected case per family, we performed a stepwise logistic regression to see whether the association with SNP9 in PTPN22 accounted for the associations were observed with other SNPs. Six SNPs, SNP2, SNP4, SNP8, SNP9, SNP11, and SNP14 were included in the final model according to the Akaike Information Criterion (AIC). Among them, SNP2, SNP4, SNP8, and SNP11 did not show the associations with RA status in the single-marker analysis, suggesting the existence of the epistatic effect between them. SNP3, SNP5, SNP7, and SNP13 were in moderate LD with SNP9 and were not included in the final model, indicating that the associations observed with these SNPs were accounted for by SNP9.

Association of RA status with other candidate genes

Power comparison

Our simulation results showed that the weighted composite likelihood ratio test [3] had a power of 89.4%, which was very similar to the power of the method developed by Slager and Schaid [1] (89.2%). Both of them had higher power than the method using one randomly selected case per affected sib pair and all unrelated controls (79.0%). All three methods showed appropriate type I error rates close to the nominal level of 0.05 (0.047, 0.048, and 0.053, respectively).