%\VignetteIndexEntry{Analysis of Bead-level Data using beadarray} %\VignetteKeywords{beadarray expression analysis} %\VignettePackage{beadarray} \documentclass{article} \textwidth=6.2in \textheight=8.5in %\parskip=.3cm \oddsidemargin=.1in \evensidemargin=.1in \headheight=-.3in \newcommand{\Rfunction}[1]{{\texttt{#1}}} \newcommand{\Robject}[1]{{\texttt{#1}}} \newcommand{\Rpackage}[1]{{\textsf{#1}}} \newcommand{\Rmethod}[1]{{\texttt{#1}}} \newcommand{\Rfunarg}[1]{{\texttt{#1}}} \newcommand{\Rclass}[1]{{\textit{#1}}} \begin{document} \title{Analysis of Bead-level Data using beadarray} \author{Mark Dunning} \maketitle \section*{Introduction} \Rpackage{beadarray} is a package for the pre-processing and analysis of Illumina BeadArray. The main advantage is being able to read raw data output by Illumina's scanning software. Data presented in this form are in the same format regardless of the assay (i.e expression, genotyping, methylation) being performed. Thus, beadarray is able to handle all these types of data. Many functions within beadarray have been written to cope with this flexibility. The BeadArray technology involves randomly arranged arrays of beads, with beads having the same probe sequence attached colloquially known as a bead-type. BeadArrays are combined in parallel on either a rectangular chip (BeadChip) or a matrix of 8 by 12 hexagonal arrays (Sentrix Array Matrix or SAM). The BeadChip is further divided into strips on the surface known as sections, with each section giving rise to a different image when scanned by BeadScan. These images, and associated text files, comprise the raw data for a beadarray analysis. However, for BeadChips, the number of sections assigned to each biological sample may vary from 1 on HumanHT12 chips, 2 on HumanWG6 chips or sometimes ten or more for SNP chips with large numbers of SNPs being investigated. This vignette demonstrates the processing of bead-level data using beadarray. The example dataset is taken from an early expression study using a BeadArray platform that is no longer commercially available. \section*{Citing beadarray} If you use \Rpackage{beadarray} for the analysis or pre-processing of BeadArray data please cite: Dunning MJ, Smith ML, Ritchie ME, Tavar\'e S, \textbf{beadarray: R classes and methods for Illumina bead-based data}, \emph{Bioinformatics}, \textbf{23}(16):2183-2184 \section{Asking for help on beadarray} Wherever possible, questions about \Rpackage{beadarray} should be sent to the Bioconductor mailing list\footnote{\texttt{http://www.bioconductor.org}}. This way, all problems and solutions will be kept in a searchable archive. When posting to this mailing list, please first consult the \emph{posting guide}. In particular, state the version of \Rpackage{beadarray} and {\tt R} that you are using\footnote{This can be done by pasting the output of running the function \texttt{sessionInfo()}.}, and try to provide a reproducible example of your problem. This will help us to diagnose the problem. <>= options(width = 80) @ <>= library("beadarray") data("BLData") @ \section{Reading bead-level data into beadarray} \subsection{File formats} The raw images and text files required to perform a bead-level analysis are produced by Illumina's BeadScan software. Usually, it will be necessary for you to modify BeadScan's default settings to obtain bead-level data, see \texttt{http://www.compbio.group.cam.ac.uk/Resources/illumina}. The command to read bead-level data from the current working directory is as follows. The \Rfunarg{dir} argument may be used to specifiy an alternative location. <>= BLData = readIllumina(useImages=FALSE, illuminaAnnotation = "Humanv3") @ The \Rfunarg{useImages} argument specifies whether beadarray will read foreground and background intensities from the TIFF images present in the directory, allowing users to experiment with strategies for image processing. In this example we set {\tt useImages=FALSE} (often a convenient choice), and locally background corrected intensities will simply be extracted from the \texttt{txt} files. The \emph{optical} background-correction that is refered to here is done by subtracting the \emph{background} pixel intensities surrounding each bead. It should not to be confused with another background correction further along the analysis pipeline, which may involve negative control beads to account for non-specific binding. Note that it is not compulsory to specify which type of Illumina assay was used to generate the data. However, for expression data it is convenient to specify the name of the platform using one of the strings Humanv4, Humanv3, Humanv2, Humanv1, Mousev2, Mousev1p1, Mousev1 or Ratv1. \Rpackage{beadarray} is able to use some of Illumina's files during analysis. These include {\tt .locs} files, which contain the locations of {\it all} beads on an array (not just those that were decoded), and {\tt .sdf} files, which contain information about the physical layout of the chip. In combination, using these files can result in significant time improvements to the detection of spatial artefacts and add additional information to some QA plots. These files are not read automatically, but if present, the path to these files is stored by beadarray for future use. If the \emph{metrics} file generated by BeadScan is present in the directory, it will be read unaltered and stored. \section{The beadLevelData class} Once imported, the bead-level data are stored in an object of class \Rclass{beadLevelData}. This class can handle raw data from both single channel and two-colour BeadArrays. Due to the random nature of the technology, each array generally has a variable number of rows of intensity data, and we use an R environment variable to store this information in a memory efficient way. The beadLevelData class contains a number of slots useful for describing Illumina data. The data that have been extracted from the text files are found in the beadData slot. This can be thought of as a list, which can be indexed by name or a numeric value representing a particular array-section. A data frame holds the data for that array-section, with the number of rows being the number of beads on the section. For convenience, the function \Rfunction{getBeadData} is used to access data held in the beadData slot. The function \Rfunction{insertBeadData} can be used to assign new data to this slot. Data types with one value per array-section can be stored in the sectionData slot. For instance, any metrics information present in the directory used by \Rfunction{readIllumina} will be stored here. This is also a convenient place to store any QC information derived during the pre-processing of the data, as we will see. The numeric identifiers for the bead-types in the \Rclass{beadLevelData} are known as ArrayAddress IDs in Illumina's annotation files. For downstream analysis it is convenient to convert these into the form ILMN\_... used in most annotation packages. Mapping objects to convert these IDs are supplied with beadarray in the {\tt extdata} directory, but this conversion may be performed automatically if the annotation of the \Rclass{beadLevelData} object is known. For this example dataset, beads that could not be decoded are assigned a special ArrayAddress ID of 0. For two-channel data, the intensities from the Red channel and associated coordinates are also stored in the object. <>= data("BLData") class(BLData) slotNames(BLData) ##Get the beadData for array-section 1 BLData[[1]][1:10,] ##Alternative using accessor function getBeadData(BLData, array=1, what="Grn")[1:10] ##Get unique ProbeIDs. These are the ArrayAddressIDs uIDs = unique(getBeadData(BLData, array=1, what="ProbeID")) uIDs[1:10] @ \section{Transformation Functions} A more flexible way to obtain per-bead data from a beadLevelData object is to define a transformation function that takes as arguments the beadLevelData object and an array index. The function then manipulates the data in the desired manner and returns a vector the same length as the number of beads on the array. The \Rfunction{logGreenChannelTransform} is the default transformation in many plotting / QA functions within beadarray. Users with two-channel data may also wish to experiment with the similarly defined \Rfunction{logRedChannelTransform} or \Rfunction{logRatioTransform} when plotting. <>= log2(BLData[[1]][1:10,2]) logGreenChannelTransform logGreenChannelTransform(BLData, array=1)[1:10] logRedChannelTransform @ In this example dataset, the local background-corrected intensities were not read from text files and separate foreground and background intensities were calculated for each bead (option {\tt useImages = TRUE}). The simple background correction that subtracts background from foreground is implemented in the \Rfunction{backgroundCorrectSingleSection} function, and creates a {\tt Grn.bc} column in the {\tt beadData} slot for each section. <>= for(i in seq_len(dim(BLData)["nArrays"])) { BLData = backgroundCorrectSingleSection(BLData, array=i) } head(BLData[[1]]) @ \section{Boxplots and imageplots} Two standard quality assessment plots supported by \Rpackage{beadarray} are the imageplot and boxplot. Boxplots can be used to compare foreground and background intensities between arrays. Image plots can be used to identify spatial artefacts on the array surface that can occur from mis-handling or scanning problems. With the raw bead-level data, we can plot false images of each array. This kind of visualisation is not possible when using the summarised BeadStudio output, as the summary values are averaged over spatial positions. Image plots in R are also more convenient than scrutinising the original tiffs, as multiple arrays can be visualised on the one page. By default, the array surface is plotted with the longest edge going horizontally. Both the \Rfunction{boxplot} and \Rfunction{imageplot} functions take a transformation function as an argument, with the default to do a $\log_2$ transformation on the green channel. In the code we show how to extract the background-corrected intensities that we have just calculated and display them on the boxplot. <>= getBackgroundCorrectionIntensities = function(BLData, array){ log2(getBeadData(BLData, array=array, what="Grn.bc")) } par(mfrow=c(1,2)) boxplot(BLData, las=2, outline=FALSE, ylim=c(4,12), main="Green Foreground\nwithout background correction", col="darkolivegreen") boxplot(BLData, las=2, transFun=getBackgroundCorrectionIntensities, outline=FALSE, ylim=c(4,12), main="Green Foreground\nafter background correction", col="darkseagreen") @ The imageplot can be configured in many ways (see manual page for more details). Sections from a BeadChip often have one edge that is much longer than the other, and it is important to recognise this when producing the plots. By default, \Rpackage{beadarray} makes imageplots with the longest edge on the x-axis (suitable for widescreen monitors). However, with {\tt horizontal = FALSE}, the imageplot will be displayed in the same orientation as the original TIFF image from the directory. With the \Rfunarg{squareSize} we can control how many pixels from the original image make up the pixels in the resulting imageplot. The following code produces imageplots for the first three array-sections in the example dataset. Note that we also change the colour scheme to represent low and high intensities by light and dark green respectively. If {\tt .locs} information is available to \Rpackage{beadarray}, it will be able to determine the optimal \Rfunarg{squareSize} parameter. If not (as with our example dataset), the user may have to experiment with different values for \Rfunarg{squareSize}. <>= par(mfrow=c(1,3), mar = c(2,2,2,2)) imageplot(BLData, array=1, low="lightgreen", high="darkgreen", horizontal = FALSE ,squareSize=25) imageplot(BLData, array=2, low="lightgreen", high="darkgreen", horizontal = FALSE ,squareSize=25) imageplot(BLData, array=3, low="lightgreen", high="darkgreen", horizontal = FALSE ,squareSize=25) @ \section{BASH} BASH is a method for managing the spatial artefacts that may be found on an array as described in Cairns et al (2008). BASH uses the methodology developed for the Harshlight package, but altered to exploit the availability of replicated observations on the same array. The algorithm first identifies Extended defects, where an array has gradual but significant shifts across the surface. BASH also seeks to find more localized artifacts on arrays by classifying features that have unusual intensities as outliers and then finding outliers close to each other on the array. Two separate algorithms then search for areas with a larger numbers of outliers than would be expected by chance (Diffuse Defects) and large connected clusters of outliers (Compact defects). The random nature (both in position and numbers of each feature type) of Illumina arrays mean that the Harshlight algorithm must proceed in a different way to the original Harshlight implementation. Whereas Affymetrix probes have replicates on other arrays, Illumina beads are replicated on the same array. We can therefore generate an error image based on how much each bead differs from the median of its replicates' intensities, instead of replicates on other arrays. Having performed manipulations to the error image, we can then find outliers on this image by bead type, determining which beads are more than 3 Median Absolute Devations, or MADs, from the median. Finally, since Illumina arrays are randomly arranged and use a hexagonal grid rather than rectangular, BASH has it's own method for creating networks of beads on the array. However, if {\tt .locs} files are available to \Rpackage{beadarray} the time taken for this step will be improved considerably. <>= bsh = BASH(BLData, array=1:10) @ The result of bash includes quality control scores; the number of beads masked in total and the extended score. <>= bsh$QC par(mfrow=c(1,2)) barplot(bsh$QC[,1], main="Number of beads masked") barplot(bsh$QC[,2], main="Extended Score") @ The weights themselves can be stored using setWeights. These will be taken into when summarizing the bead-level data. The QC tables themselves can be appended to the {\tt sectionData} slot of \Robject{BLData}. <>= for(i in 1:10){ BLData = setWeights(BLData, wts=bsh$wts[[i]], array=i) } BLData = insertSectionData(BLData, what="BASHQC", data = bsh$QC) @ \section{Using control information} Illumina have designed a number of control probes for each expression platform. For expression arrays, we store the ArrayAddressIDs of the control probes for in the ExpressionControlData object. Otherwise a data frame may be used to define these ids. As the example data described in this vignette were derived using an obsolete technology, we have stored the control information with the package in the \Robject{controlProfile} object. ArrayAddressIDs are listed in the first column, and the type of control in the second column. Objects of this form can be used in various quality assessment functions in beadarray. <>= data(controlProfile) head(controlProfile) table(controlProfile[,2]) @ Two particular controls on expression arrays are housekeeping and biotin controls. With the \Rfunction{poscontPlot} function, we can plot the intensities of any ArrayAddressIDs that are annotated as belonging to the Housekeeping or Biotin group in the \Robject{ExpressionControlData} object or \Robject{controlProfile}. The \Rfunction{poscontPlot} is flexible in allowing other "tags" in the \Robject{controlProfile}, in the example below we configure the plot to show both housekeeping and negative controls in the same plot. <>= par(mfrow=c(2,5), mar = c(2,2,2,2)) for(i in 1:10){ poscontPlot(BLData, array=i, controlProfile=controlProfile, ylim = c(10,16)) } @ <>= par(mfrow=c(2,5), mar = c(2,2,2,2)) for(i in 1:10){ poscontPlot(BLData, array=i, controlProfile=controlProfile,positiveControlTags = c("housekeeping", "negative"), ylim = c(10,16)) } @ With knowledge of which ArrayAddressIDs match control types, we can easily provide summaries of these control types on each array. In \Rfunction{quickSummary} the mean and standard deviation of all control types is taken for a specified array, using intensities of all beads that correspond to the control type. Note that these summaries may not correspond to similar quantities reported in Illumina's BeadStudio software, as the BeadStudio summaries are produced after removing outliers (see later). The \Rfunction{makeQCTable} function extends this functionality to produce a table of summaries for all sections in the \Rclass{beadLevelData} object. These data can be stored in the sectionData slot for future reference. It is also informative to compare the expression level of various control types to the background level of the array. This is done by the \Rfunction{controlProbeDetection} function that returns the percentage of each control type that are significantly expressed above background level. Obivously for positive controls we would prefer this percentage to be near 100 on a good quality array. <>= quickSummary(BLData, array=1, reporterIDs=controlProfile[,1], reporterTags=controlProfile[,2]) qcReport = makeQCTable(BLData, controlProfile=controlProfile) head(qcReport)[,1:5] BLData = insertSectionData(BLData, what="BeadLevelQC", data=qcReport) names(BLData@sectionData) for(i in 1:10){ print(controlProbeDetection(BLData, array = i, controlProfile=controlProfile, tagsToDetect=c("housekeeping", "biotin"), negativeTag="negative")) } @ The generation of QA plots for all sections in the beadLevelData object is provided by the expressionQCPipeline function. Results are generated in a directory of the users choosing. This report may be generated at any point of the analysis. If the overWrite paramater is set to FALSE, then any existing plots in the directory will not be re-generated. Futhermore, QC tables that have been stored in the beadLevelData object already can be used. <>= expressionQCPipeline(BLData, qcDir="QC") expressionQCPipeline(BLData, controlProfile=controlProfile, positiveControlTags=c("housekeeping", "biotin"), hybridisationTags = c("cy3_hyb", "high_stringency_hyb"), zlim=c(10,12), horizontal=F, negativeTag = "negative", tagsToDetect=list(Housekeeping = "housekeeping", Biotin = "biotin", hybridisation = "cy3_hyb")) @ \section{Outlier removal and plotting} Before combining the observations for each bead-type on an array, Illumina remove any observations with outlying intensity (more than 3 median absolute deviations from the median). This step can be repeated in \Rpackage{beadarray} and it is sometimes useful to see where these outliers are located on the array surface. Often, they will coincide with beads masked by BASH or with any spatial artefacts that may be seen. Users are able to define their own functions to identify outliers. Such functions must take a list of intensities and corresponding ArrayAddressIDs and return indices of which observations are found to be outliers. <>= par(mfrow=c(1,3), mar = c(2,2,2,2)) outlierplot(BLData, array=1, horizontal = FALSE) outlierplot(BLData, array=2, horizontal = FALSE) outlierplot(BLData, array=3, horizontal = FALSE) @ \section{Summarization} The summarization procedure takes the \Robject{BLData} object, where each bead-type is represented by differing numbers of observations on each array, and produces a summarized object to make comparisons between arrays. For each array section represented in the \Robject{BLData} object, all observations are extracted, transformed, and then grouped together according to their ArrayAddressID. Outliers are removed and the mean and standard deviation of the remaining beads are calculated. The \Rclass{illuminaChannel} class is used to define how summarization proceeds with specification of a transformation function, a function to remove outliers and function to calculate the means and standard deviation. The code below creates two different summarized objects; one which uses mean and standard deviations, and one which uses median and standard errors. <>= myMean = function(x) mean(x,na.rm=TRUE) mySd = function(x) sd(x,na.rm=TRUE) greenChannel = new("illuminaChannel", logGreenChannelTransform, illuminaOutlierMethod, myMean, mySd,"G") BSData <- summarize(BLData, list(greenChannel)) myMedian = function(x) median(x, na.rm=TRUE) mySe = function(x) sd(x, na.rm=TRUE)/sqrt(length(x)) greenChannel2 = new("illuminaChannel", logGreenChannelTransform, illuminaOutlierMethod, myMedian, mySe,"G") BSData2 <- summarize(BLData, list(greenChannel2)) @ <>= BSData head(exprs(BSData)[,1:4]) head(se.exprs(BSData)[,1:4]) BSData2 head(exprs(BSData2)[,1:4]) head(se.exprs(BSData2)[,1:4]) @ The \Robject{BSData} object is very similar to the \Rclass{ExpressionSet} class in Biobase. However, to accomodate the unique features of Illumina data we have added an {\tt nObservations} slot, which gives the number of beads that we used to create the summary values for each bead-type on each array after outlier removal. It is possible to have multiple channels, each of which is summarized in a different manner, in the same \Robject{ExpressionSetIllumina} object. This is achieved by passing a list of illuminaChannel objects to \Rfunction{summarize}. This would be especially useful for two-channel data, where the Red and Green channels, and some combination of the two would be of interest in the analysis. In the example code below we summarize both the non background-corrected and background-corrected intensities in the same object. The \Rfunction{channel} function is used to select the data for one of these channels, which returns an of object of type \Rclass{ExpressionSetIllumina} <>= greenBackgroundCorrected = new("illuminaChannel", getBackgroundCorrectionIntensities, illuminaOutlierMethod, myMean, mySd, "G.bc") BSData.multChannel = summarize(BLData, channelList = list(greenChannel, greenBackgroundCorrected)) channelNames(BSData.multChannel) G = channel(BSData.multChannel, "G") G.bc = channel(BSData.multChannel, "G.bc") @ The detection score is a standard measure for Illumina expression experiments, and can be viewed as an empirical estimate of the p-value for the null hypothesis that there is no expression. These can be calculated for summarized data provided that the identity of the negative controls on the array is known. For further analysis of the summarized object, see the separate beadsummary vignette {\tt beadsummary.pdf}. <>= status = rep("regular", as.numeric(dim(BSData.multChannel)[1])) negIDs = controlProfile[which(controlProfile[,2] == "negative"),1] status[match(negIDs, featureNames(BSData.multChannel))] = "negative" det = calculateDetection(G, status=status) head(det) Detection(G) = det @ <>= sessionInfo() @ \end{document}