R/Bioconductor software for Illumina's Infinium whole-genome genotyping BeadChips.

UNLABELLED: Illumina produces a number of microarray-based technologies for human genotyping. An Infinium BeadChip is a two-color platform that types between 10(5) and 10(6) single nucleotide polymorphisms (SNPs) per sample. Despite being widely used, there is a shortage of open source software to process the raw intensities from this platform into genotype calls. To this end, we have developed the R/Bioconductor package crlmm for analyzing BeadChip data. After careful preprocessing, our software applies the CRLMM algorithm to produce genotype calls, confidence scores and other quality metrics at both the SNP and sample levels. We provide access to the raw summary-level intensity data, allowing users to develop their own methods for genotype calling or copy number analysis if they wish.
AVAILABILITY AND IMPLEMENTATION: The crlmm Bioconductor package is available from http://www.bioconductor.org. Data packages and documentation are available from http://rafalab.jhsph.edu/software.html.

1 INTRODUCTION

In recent years, large-scale genome-wide association studies have provided significant insight into the genetics underpinning many complex diseases (Grant and Hakonarson, 2008). High-density microarrays, which allow many single nucleotide polymorphisms (SNPs) to be genotyped simultaneously in a sample at low cost, have been the technology driving this research. Illumina Inc. (San Diego, CA, USA) is a major provider of such arrays.

Illumina BeadChips are composed of a number of rectangular strips, each containing many randomly arranged, replicated beads. For Infinium genotyping, beads are coupled with specific 50mer probes designed to be complementary to the sequence adjacent to the SNP site, and the two alleles (A, B) are discriminated using either a red or green dye (Steemers et al., 2006). Data are acquired by scanning each strip at different wave lengths using Illumina's scanning device followed by automatic image analysis (Galinsky, 2003). A robust summary of the intensity in each channel for each SNP assayed is reported in proprietary idat files. BeadChips of varying SNP density and sample format (single, duo, quad) are available for human genotyping. Some contain non-polymorphic probes for assessing copy number variation.

Many algorithms that take summarized alleles A and B signals as inputs to produce genotypes (AA, AB, BB) have been developed for Affymetrix SNP arrays (Carvalho et al., 2007; Hua et al., 2007; Rabbee and Speed, 2006; Xiao et al., 2007). A smaller number of Illumina-specific methods (Giannoulatou et al., 2008; Teo et al., 2007) including Illumina's GenCall algorithm in BeadStudio/GenomeStudio are also available. Software for the analysis of Illumina data such as beadarray (Dunning et al., 2007), beadarraySNP and lumi (Du et al., 2008) is available in R/Bioconductor (Gentleman et al., 2004); however, current packages do not deal specifically with Infinium BeadChip data.

In this article, we present the crlmm package for Illumina genotyping. Our software extracts summarized intensities, performs normalization and applies the CRLMM algorithm (Carvalho et al., 2007) to remove chip- and SNP-specific biases and call genotypes.

2 METHODS

To begin, summarized data are read from idat files (two per array, one for each channel) using the function readIdatFiles. Binary idat files are a convenient starting point, as they are routinely output by the scanning software, provide a compact representation of the data and have a consistent format (unlike output from Illumina's BeadStudio/GenomeStudio software, which is exported at the user's discretion, meaning the raw signals needed for the analysis are not always available). Access to the raw data allows for low-level plotting to help visualize trends and biases that may be present (Fig. 1A and B). It also allows alternative genotyping algorithms, which require data on the raw scale, to be applied.

Fig. 1.

(A) Plot of the log2allele B (green) intensity by strip (labelled by Row.Column position), which steadily increases from rows 1 to 10 down the BeadChip. This effect is less prominent in the allele A (red) channel (data not shown). The source of this trend is related to the way post-hybridization reagents are applied to the BeadChip and scan order, and its presence motivates strip-level normalization. (B) A smoothed scatter plot of M versus S for a typical array, where darker regions indicate a higher density of points. This plot shows intensity-dependent effects in M which vary for the AA and BB genotypes, and motivate the three-component mixture model in CRLMM. The curves represent the smoothing splines that model this effect. (C) SNR for 60 arrays, with the median (solid line) and median–median absolute deviation (dashed line) SNR values plotted. Lower scores correspond to poorer separation between the genotype clouds depicted in (B). This metric can be used to flag low-quality arrays to exclude from further analysis.

Next, the allele A (X Raw) and allele B (Y Raw) signals are normalized between channels and samples simultaneously using strip-level quantile normalization. The between-channel aspect of the normalization [also recommended in Oosting et al., (2007) and Staaf et al., (2008)] aims to remove any dye-bias effects, while the strip-level component corrects for intensity gradients which can occur within BeadChips (Fig. 1A). Normalization at the strip-level has also proven useful for data from Illumina's gene expression BeadChips (Wei Shi, personal communication). By default, the strip-level quantiles are standardized against a reference distribution obtained from HapMap samples (International HapMap Consortium, 2007) run on the same platform to correct for lab and batch effects.

After normalization, the CRLMM genotyping algorithm (Carvalho et al., 2007; Lin et al., 2008) is applied. For each array, SNP-specific log-ratios (M = log2alleleA − log2alleleB) and average intensities [S = (log2alleleA + log2alleleB)/2] are calculated. As noted for Affymetrix data (Carvalho et al., 2007), S appears to have an effect on M. The effect appears to be a smooth function of S, but only applies to the AA and BB intensities (Fig. 1B). To remove this effect, we fit a three-component mixture model with a spline used to model the smooth function. This model is fitted per array via the expectation-maximization (EM) algorithm using a random sample of data-points. Due to the different chemistry used for Illumina genotyping, the fragment length covariate described in Carvalho et al. (2007) can be ignored.

Next, a two-level hierarchical model is applied. SNP-specific means and standard deviations (SDs) are obtained for each genotype via supervised learning using HapMap data. Independent genotype calls (available from http://www.hapmap.org/) provide the true states for samples that have been genotyped using the respective BeadChip platform. Normalized signals from these arrays are then used to estimate robustly the genotype means and SDs. The intensity-dependent splines from the EM (which explain the between-SNP variation) and the SNP-specific genotype means and SDs (obtained from training data) are combined in the model. New genotype calls are assigned by choosing the class that minimizes the negative log likelihood.

CRLMM produces a number of metrics for quality assessment (Lin et al., 2008). Confidence scores for each call are provided using the log-likelihood ratio test from the hierarchical model. The sample-specific SNR (signal-to-noise ratio) assesses the separation of the three genotypes within an array. Lower SNR values indicate poorer quality, and this metric can be used to exclude samples from further analysis (Fig. 1C). SNP-specific quality is measured as the minimum distance between the heterozygote centre and either of the two homozygous centres.

The preprocessing and genotyping steps above are performed by the crlmmIllumina function. All code is written in R (R Development Core Team, 2009) and existing Biobase classes are used to store the data.

The software requires chip-specific data packages (available at http://rafalab.jhsph.edu/software.html) that store basic SNP annotation information and various parameters used by CRLMM. We also provide the hapmap370k data package, which contains idats from 40 HapMap samples hybridized to HumanHap 370 K Duo BeadChips, and a user guide that provides example R code to analyse these samples (see Supplementary Material). A 64-bit Linux system takes ∼90 s and uses up to 1.2 GB of RAM to read these data, while normalization and genotyping takes a further 470 s and uses up to 3.3 GB of RAM. This equates to processing around 600 SNPs per second.

3 DISCUSSION

The crlmm package provides bioinformaticians with an additional tool outside of Illumina's proprietary software for analysing Infinium BeadChip data. Our software also facilitates the analysis of Affymetrix SNP chips and the use of a consistent algorithm and framework to process both the platforms allows data from different studies to be combined more easily.

Implementation in R/Bioconductor gives users the opportunity to exploit other tools that have been adapted for Illumina data. For example, if raw bead-level data were available, the BASH spatial artefact detection method (Cairns et al., 2008) in the beadarray package could be applied. Once summarized, the data could be further processed using crlmm.

The CRLMM algorithm can be applied to new versions of Illumina BeadChips for humans and other species, provided that the necessary training data and prior information on genotype calls are available. Future work will benchmark the performance of our method with other genotyping algorithms tailored to suit Illumina data, such as Illuminus (Teo et al., 2007) and Illumina's own algorithms in BeadStudio/GenomeStudio. Furthermore, tools for copy number analysis are being developed in the crlmm package.