A Bayesian latent class analysis for whole-genome association analyses: an illustration using the GAW15 simulated rheumatoid arthritis dense scan data

Rheumatoid arthritis (RA) is the most common cause of inflammatory polyarthritis in adults [1]. This chronic, inflammatory disease has a complex genetic component involving the human leukocyte antigen (HLA) region. For nearly two decades the association between the HLA region and RA has been known and confirmed in numerous population studies [2]. The HLA associations are extremely complex and their exact biologic role in RA is unknown. Furthermore, family and twin studies have indicated non-HLA genes may play an important role in RA.

Advancements in genotyping technology have facilitated the ability to generate large amounts of genetic data. The large numbers of single-nucleotide polymorphisms (SNPs) genotyped by whole-genome scans may overwhelm conventional statistical approaches such as maximum-likelihood estimation (MLE). A few novel statistical approaches reducing the dimensionality of large data sets and detecting the structural relationship between variables have been described elsewhere [3–5]. Simple association tests may produce a large number of false positives; therefore newer statistical approaches are needed to incorporate known information on disease etiology, thus reducing the potential for these false-positive associations. Bayesian analysis might be better able to detect true associations while minimizing false positives, because it can "borrow strength" from multiple loci to estimate association parameters and it can incorporate differences across loci in the prior probability of association. Appropriate Bayesian analysis should also reduce parameter estimate variability, similar to penalized regression methods such as ridge regression and the lasso (which are themselves special cases of Bayesian analysis [6]). As a methodological exercise, we contrast a novel Bayesian latent class analysis with MLE in a simulated data set of RA for chromosome 6 markers from the Genetic Analysis Workshop 15 (GAW 15 Problem 3). Although our analysis is restricted to a subset of the data that would be available from a genome-wide scan, in principle, our method could be applied to a whole-genome scan.

We randomly selected one RA case from the affected-sibling pair (ASP) in the first replicate of the GAW15 simulated data for Problem 3. After selecting all of the controls, our final data set included 1500 cases and 2000 controls. In order to ensure a significant finding we reviewed the answers to the simulated data prior to our analyses. Since the strongest signal for RA was simulated to the HLA region on chromosome 6, we limited our analyses to the dense genotyping for chromosome six. In total 17,820 SNPs were simulated on chromosome 6, yielding an average inter-marking spacing of 9586 base pairs. This corresponds to the density one would expect from a genome-wide 300,000 K SNP set.

We developed and implemented a Bayesian latent class analysis because we hypothesized that by "borrowing strength" across multiple loci and incorporating prior probabilities of association, such an analysis might be more sensitive and specific than p-values from maximum-likelihood based tests. We applied this latent class analysis to the GAW15 simulated chromosome 6 data, but found that the latent class models provided no obvious improvement over MLE. However, our overall results may not be generalizable due to the large simulated effect in the HLA-DR locus.

As can be seen in Figures 1 and 2, the Bayesian posterior means for β is less variable than the MLE. This suggests that if researchers are choosing markers to follow up on the basis of estimated effect size, some sort of smoothing procedure could be useful. MLEs for rare SNPs may be very unstable. Smoothing MLEs to a group mean will account for differences in information across SNPs and hence could reduce the false-positive report probability [11].

On the other hand, the posterior probability of association Pr(X > 0) from the latent class analyses was not as good as simple p-values distinguishing a signal from a non-signal. Negative consequences, although minor, occur when up-weighting regions that are not true candidates (Table 1). In the GAW15 simulated data set, when the HLA signal was extremely large, the weighting of candidate regions did not appear helpful.

There are clearly drawbacks to the latent class approach as we have implemented it. The posterior probability of belonging to a non-null class, Pr(X ≠ 0), is much too large for a majority of the loci. Setting f(β) equal to point mass at 0 may be too stringent. One potential solution is to allow for some noise to distinguish signals near, but not directly at, zero from true effects. Additionally, it appears the priors are overwhelmed by the large amount of data (loci). Because a majority of the loci are likely null, the prior sample size (e.g., κ₀ or the absolute magnitude of π _j ) could be increased. Other drawbacks to our implementation include the facts that the ratio of priors comparing candidate regions to non-candidate regions is fixed, and only loci in regions sharing the same prior contribute to estimating the posterior odds of association in those regions. A hierarchical Bayesian approach may be more appropriate because it estimates both baseline and relative odds of association from the data, rather than fixing them [12]. By not fixing the prior odds in non-candidate regions, the hierarchical Bayesian approach may also be more sensitive to true associations in unexpected regions (e.g., regions with little biologic annotation).