Model-based analysis of two-color arrays (MA2C)
- Jun S Song†1, 2,
- W Evan Johnson†1, 2,
- Xiaopeng Zhu†3,
- Xinmin Zhang4,
- Wei Li1, 2,
- Arjun K Manrai5,
- Jun S Liu2, 6,
- Runsheng Chen3 and
- X Shirley Liu1, 2Email author
© Song et al.; licensee BioMed Central Ltd. 2007
Received: 20 April 2007
Accepted: 29 August 2007
Published: 29 August 2007
A novel normalization method based on the GC content of probes is developed for two-color tiling arrays. The proposed method, together with robust estimates of the model parameters, is shown to perform superbly on published data sets. A robust algorithm for detecting peak regions is also formulated and shown to perform well compared to other approaches. The tools have been implemented as a stand-alone Java program called MA2C, which can display various plots of statistical analysis for quality control.
The technology is continuing to develop rapidly, but certainly not without difficulties that are imposed by the inherent complexity of biological systems and, as such, must be addressed by computational means for the foreseeable future. The main computational challenge lies in properly normalizing the data and distinguishing true peaks from the noisy background. Many problems that confound this type of microarray data actually arise from probe-specific biases, such as differential sequence copy numbers in the genome or variable melting temperature dependent upon the GC content. For Affymetrix tiling arrays, several good model-based methods already exist to account for probe biases and, thus, to adjust for probe-specific baseline signals. The recently introduced MAT , for instance, estimates probe affinity from probe sequence and copy number and provides a powerful tool for finding enriched regions in chromatin immunoprecipitation (ChIP) and other applications on Affymetrix tiling-array experiments. Incidentally, similar problems are also found in Affymetrix expression arrays, for which extensive effort has been previously exerted by various groups to develop robust methods for background correction and probe-level normalization (for example, [2–5]). It is relatively hard and expensive for Affymetrix to provide custom designed microarrays.
Commercial custom tiling arrays are relatively new in the field of microarray biotechnology and, just as expression arrays allow global assays of gene expression, provide an invaluable tool for investigating the locations and roles of DNA-binding proteins in the whole genome at high resolution. All currently available custom tiling arrays use the two-color technology. Considering the utility and power of high-resolution tiling arrays, it is thus imperative that reliable computational methods be developed now to facilitate the extraction of precise and accurate conclusions from such experiments.
For dual-channel cDNA arrays, several normalization methods have been proposed (for example, [2, 6]), but these procedures typically utilize methods that neglect probe sequence information and are also computationally expensive and, thus, unsuitable for currently available high-density tiling arrays. One common way of locally normalizing two-color arrays is the so-called M-A loess normalization. The fundamental assumption behind this procedure is that most probes should have similar values between the two-channels, an assumption violated in studies of chromatin structure such as nucleosome mapping described in [7, 8]. This method also does not account for sequence-specific effects, which may be significant in high-density tiling arrays, and also does not normalize the variance of M.
Single-channel normalization methods can be also applied to two-color arrays, such as those proposed by [3, 9], but they ignore the fact that the two channels are paired, and such approaches are thus likely to retain residual effects or correlation. Recently, Dabney and Storey  have introduced a normalization method that adjusts for intensity-dependent dye bias and array-to-array variations. However, their method, which was developed for expression arrays, does not model sequence-specific probe effects and is based on smoothing procedures that can be computationally demanding for tiling arrays; the approach also requires a dye swap and, thus, cannot be applied to single array experiments, which are often performed as test runs. In fact, as far as we are aware, there are, to date, only two published tools, MPeak [11, 12] and ChIPOTle [13, 14] for analyzing two-color high-density tiling arrays, but neither considers probe-specific normalization or is able to combine replicate experiments directly. This problem is rather serious since biological replicate experiments are perceived to be indispensable in any sound research utilizing microarrays.
In this paper, we address many of the issues discussed above and present robust algorithms for normalizing the raw data at probe-level and detecting peaks, implemented as a Java program called MA2C (model-based analysis of two-color arrays). Because our normalization method standardizes the probe intensities, our peak-detection algorithm naturally generalizes to combine replicate arrays.
Results and discussion
Comparison of normalization methods
We used the data (GEO GSE7523) from a recent spike-in experiment to test MA2C. The spike-in samples contained 96 clones in the ENCODE region of approximately 500 bp, at 8 different concentrations corresponding to (2n + 1)-fold enrichment compared to the human genomic DNA, for n = 1,...,8, and 12 different clones per concentration. The control sample contained sonicated genomic DNA without spike-ins. The spike-in and control samples were differentially labeled and hybridized to a NimbleGen ENCODE tiling array in triplicates, and the resulting data were used to assess the performance of MA2C against other currently available algorithms.
Comparison of MA2C with other algorithms using a spike-in experiment with a total of 96 regions and 47 unique non-overlapping regions
(C = 2 normalized)
The positive predictive value of MPeak was comparable to MA2C, but MA2C was more sensitive and also found more unique sites. MA2C again showed a better correlation with spike-in fold-changes than MPeak and, thus, provided better quantitative information about the enriched regions than both ChIPOTle and MPeak. We also tested the MA2C peak detection algorithm on the global median-scaled data without any GC-correction (the same data analyzed with MPeak and ChIPOTle) and still found MA2C to be more sensitive and to have a higher positive predictive value, indicating that MA2C can outperform other available algorithms even without its GC-specific normalization step (Table 1).
Furthermore, neither MPeak nor ChIPOTle can combine replicate data in a single test. As seen in Table 1, pooling data from replicate experiments can often increase the sensitivity and quantitativeness of analysis, and this option implemented in MA2C will prove to be useful. Since ChIP-chip experiments require biological replicates, which are much noisier than the technical triplicate spike-ins presented here, the ability to combine replicates at the probe-level will provide more sensitive and robust peak predictions than other methods of combining peaks. In addition, ChIP-chip experiments contain a PCR amplification step that often increases the GC bias of probes; in this regard, MA2C's GC-based probe normalization shows distinct advantages over ChIPOTle and MPeak on PCR amplified samples, as observed in a separate PCR amplified spike-in experiment (unpublished data).
ChIP-chip data in Caenorhabditis elegans
The protein DPY-27 functions as an essential dosage-compensator that suppresses the expression of genes on each X chromosome in hermaphrodite XX embryos of Caenorhabditis elegans, thereby reducing the expression level of the X-linked genes by half to the level in XO (male) counterparts. Chuang et al.  have shown that the basic suppression mechanism involves localization of DPY-27 to X chromosomes, likely leading to a subsequent modification of the chromatin structure of X chromosomes mediated by DPY-27. Davis and Meyer  later showed that SDC-3 also localizes to X chromosomes in XX hermaphrodites and associates with a dosage compensating complex involving DPY-27.
Numbers and annotation of SDC-3 binding sites detected by different methods
No. of peaks
Combined triplicate (p = 10-5)
Combined triplicate (p = 10-4)
Overlap of binding sites of SDC-3
MA2C (p = 10-5)
MA2C (p = 10-5)
ChIP-chip technology has quickly become popular among biologists, and high-density tiling microarrays are increasingly being used in novel genomic research. Some of the interesting applications involve finding novel transcripts in the genome, DNA methylation sites, nucleosome positions, DNA hypersensitivity regions, and alternative splicing events [7, 8, 18–21].
In all of these studies, which tend to combine experiments performed at various time points and under different conditions, the variability of array performance and sequence-specific effects must be addressed properly in order to remove any technical artifacts and to be able to formulate biologically sound conclusions. The problem of probe effects becomes more pronounced as the density of tiling increases, as one does not have the option of selecting probe sequences for similar melting temperature, or when the tiled regions predominantly cover promoter regions, which are known to be GC-rich. Our method of standardization explicitly accounts for such sequence-specific biases and inter-array variability. Together with the accompanying robust peak-detection algorithm, MA2C's standardization procedure is especially important for data sets with a significant noise level - for instance, stemming from PCR amplification, which tends to increase probe effects.
One issue we have not discussed so far is adjusting for the copy-number of probes or cross-hybridization of DNA with similar sequences. We chose not to model the sequence copy-number because both NimbleGen and Agilent use sufficiently long probes and also usually exclude repeat regions from their array design.
It is also instructive to note why our normalization method in equation 1 or equation 3 (See Materials and methods) gives a higher weight to the probes that are highly correlated between the two channels. Relying on the fact that the probes are long, NimbleGen tends to wash their arrays rather harshly after hybridization, minimizing cross-hybridization but also possibly leaving behind only random noise and causing a low correlation in low-GC probes between the two channels. Thus, as illustrated in Figures 2 and 3a, the low-GC probes are mostly measuring the background noise and also show a low inter-channel correlation; this relation between low intensity distribution and low inter-channel correlation in low-GC bins is the motivation behind MA2C's normalization method.
Materials and methods
where n k is the number of probes with GC = k. We further scale the t-values globally so that the rescaled t-values have variance 1.
The t-values thus yield log-ratios adjusted by the mean and normalized by the standard deviation within each GC bin.
Note that in equation 1, the covariance term ξ k has the effect of amplifying the difference between experiment and control probe intensities in GC bins that have a high baseline correlation between the two channels, while suppressing the difference in GC bins with low correlation. Therefore, the log-fold changes xi2 - xi1 are given more weight in GC bins with high correlation ξ k between the two channels than in low-correlation GC bins.
Robust estimation of parameters
where C is a fixed constant and M = median i |x i - μ*|, the median absolute distance. We then calculate the bi-weight for each data point as wi = (1 - )2 for -1 ≤ d i ≤ 1 and w i = 0 otherwise. Then, the mean is re-estimated as , and the process is repeated until a certain convergence criterion is satisfied.
and . Here, Z i is the projected vector previously described and Σ its variance matrix. In each iteration step, the mean is estimated as before, the variance as , and likewise for the covariance. Strictly speaking, the variance and covariance computed in this way are not consistent estimators, but as shown in the Results section, they do provide reasonable estimates of the parameters requisite for standardizing the data within and across arrays.
Detection of peak regions
To detect peak regions, we have implemented several adaptations of the powerful window-based approach proposed by Johnson et al.  for Affymetrix tiling arrays. More precisely, we consider a sliding window of some user-defined length (400 bp to 1,000 bp) centered at each probe. A MA2Cscore is assigned using a user-selected scoring function based on median, pseudo-median, median polish, or trimmed mean of the probes in the window. The median and trimmed mean options are implemented by calculating the median or trimmed mean of all the probes in the window; when replicates are available, the median t-value or trimmed mean of all pooled probes in identical windows across replicates is used. The pseudo-median of a distribution is the median of all pairwise arithmetic means, as discussed in . Median polish has been successfully applied in robust multi-chip analysis for Affymetrix gene expression arrays . We recommend using median polish for experiments with a large number of replicate samples, while trimmed mean is recommended for arrays with densely tiled probes. The pseudo-median and median provide robust alternatives that can be applied in experiments that are not densely tiled and have few available replicates.
To compare the performance of different scoring functions, the triplicate H3 acetylation data at 38 bp spacing from  were analyzed using a window size of 1,000 bp and a p value cutoff of 10-3. Median polish gave the most number of peaks while trimmed mean gave the least, the difference in number being around 3%. The agreement among median, pseudo-median and trimmed mean was around 97-99%, while median polish agreed with other methods by 93-97%. Comparable results were obtained, with 1-2% less agreement, when the data were re-analyzed at 76 bp spacing by skipping probes. The best agreement was found between trimmed mean and pseudo-median at 99-100% while the worst agreement was between median and median polish at 90-93%.
To increase reliability, windows containing less than k probes are discarded, where k is again defined by the user. Just like MATscores, MA2Cscores approximately follow a normal distribution, with the representative scores of peak regions corresponding to the right tail. This fact easily allows us to assign a p value to each MA2Cscore using the normal probability distribution. The lower-bound of MA2Cscores for determining peaks may be based on either false discovery rate (FDR) or p value computations. As in , we empirically estimate FDR as follows: for a given cutoff value M > 0 of MA2Cscore, we find all peaks with MA2Cscore greater than M and all peaks with MA2Cscore less than -M. Then, FDR is estimated as #(negative MA2Cscore peaks)/#(positive MA2Cscore peaks), and the number of true positive peaks as #(positive MA2Cscore peaks) - #(negative MA2Cscore peaks). The FDR table, along with other informative histograms, is generated by MA2C.
We have implemented our method as a user-friendly, stand-alone Java package called MA2C, which is fully automated and only requires the user to select the directory path and treatment channels.
The file structure of NimbleGen data consists of three main components, DesignFiles/, PairData/, and SampleKey.txt, which should all reside in the same parent directory. The text file SampleKey.txt contains the relevant design information about individual arrays; in particular, the file must contain DESIGN_ID, CHIP_ID and DYE for each array. The directory DesignFiles/contains the sequence and position files corresponding to each DESIGN_ID, while PairData/contains the single channel data for each CHIP_ID. Even though MA2C is primarily designed for NimbleGen arrays, we have also successfully tested the program on Agilent data by reformatting the necessary files and obtained excellent results.
When the user begins by selecting SampleKey.txt, MA2C reads the file and displays the content in a table. If DesignFiles/and PairData/are present in the parent path, MA2C also automatically lists the directory contents in two separate tables; otherwise, the user has to choose the corresponding folder locations. The user then selects the treatment channel for each experiment to be analyzed and clicks the Run button, which prompts MA2C to perform the normalization and peak detection steps as follows.
Step 1: DesignFiles/
For each DESIGN_ID, MA2C automatically reads the corresponding .ndf and .pos files and generates a .tpmap file containing the sequence, chromosome, position, and array coordinate information of probes.
Step 2: PairData/
For each chosen treatment channel with given CHIP_ID and DYE, MA2C searches for the correct two-channel data files. It is thus important that the pair data files contain a column corresponding to IMAGE_ID. For fast future access and also for compressed storage, the program combines each two-channel data into a single file named MA2C_CHIP_ID_raw.txt. Normalized data are similarly stored in files with the extension _normalized.txt.
Step 3: MA2C_output/
MA2C automatically creates this directory for writing files used in quality control of normalization and peak detection steps. The enriched regions are output in both .xls and .bed files that contain the chromosome, start, end, p value, MA2Cscore and peak-center information for each detected peak. MA2Cscore.bar and ratio.bar files are created for visualization using Affymetrix's Integrated Genome Browser .
MA2C is an open source Java package that can be downloaded from . MA2C runs on all platforms that support Java Runtime Environment 5.0 or higher and has been successfully tested on OS X, Linux and Windows operating systems. The program is written so as to economize the size of required files; once the .tpmap and _raw.txt files have been created, the subsequent runs of MA2C will use only those files and the user may remove the .ndf, .pos, and other pair data files. This approach can save hundreds of megabytes of disk space. In addition, our program is fast, the total execution time being usually less than a couple of minutes for multiple arrays. For example, on a laptop with a 2.13 GHz Intel M processor and 2 GB RAM, it takes 18 seconds to build a sequence file for 370,000 probes, 16 seconds to normalize the raw data, and 14 seconds to find peaks.
false discovery rate.
We thank Nathan Trinklein and Rick Myers' lab for generating the spike-in sample. We also thank the laboratories of Jason D Lieb and Bing Ren for their data and Xihong Lin and Nan Jiang for their insights on methodology. This project was partially funded by NIH grants 1R01 HG004069-01 and 1U01 HG004270-01.
- Johnson W, Li W, Meyer C, Gottardo R, Carroll J, Brown M, Liu X: Model-based analysis of tiling-arrays for ChIP-chip. Proc Natl Acad Sci USA. 2006, 103: 12457-12462. 10.1073/pnas.0601180103.PubMedPubMed CentralView ArticleGoogle Scholar
- Yang Y, Dudoit S, Luu P, Lin D, Peng V, Ngai J, Speed T: Normalization for cDNA microarray data: a robust composite method addressing single and multiple slide systematic variation. Nucleic Acids Res. 2002, 30: e15-10.1093/nar/30.4.e15.PubMedPubMed CentralView ArticleGoogle Scholar
- Bolstad B, Irizarry R, Astrand M, Speed T: A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003, 19: 185-193. 10.1093/bioinformatics/19.2.185.PubMedView ArticleGoogle Scholar
- Wu Z, Irizarry R, Gentleman R, Martinez Murillo F, Spencer F: A model based background adjustment for oligonucleotide expression arrays. JASA. 2004, 99: 909-917.View ArticleGoogle Scholar
- Hoffmann R, Seidl T, Dugas M: Profound effect of normalization on detection of differentially expressed genes in oligonucleotide microarray data analysis. Genome Biol. 2002, 3: RESEARCH0033-PubMedPubMed CentralGoogle Scholar
- Tseng G, Oh M, Rohlin L, Liao J, Wong W: Issues in cDNA microarray analysis: quality filtering, channel normalization, models of variations and assessment of gene effects. Nucleic Acids Res. 2001, 29: 2549-2557. 10.1093/nar/29.12.2549.PubMedPubMed CentralView ArticleGoogle Scholar
- Yuan G, Liu Y, Dion M, Slack M, Wu L, Altschuler S, Rando O: Genome-scale identification of nucleosome positions in S. cerevisiae. Science. 2005, 309: 626-630. 10.1126/science.1112178.PubMedView ArticleGoogle Scholar
- Ozsolak F, Song J, Liu X, Fisher D: High-throughput mapping of the chromatin structure of human promoters. Nat Biotech. 2007, 25: 244-248. 10.1038/nbt1279.View ArticleGoogle Scholar
- Schadt E, Li C, Ellis B, Wong W: Feature extraction and normalization algorithms for high-density oligonucleotide gene expression array data. J Cell Biochem Suppl. 2001, 120-125. 10.1002/jcb.10073. Suppl 37
- Dabney A, Storey J: A new approach to intensity-dependent normalization of two-channel microarrays. Biostatistics. 2007, 8: 128-139. 10.1093/biostatistics/kxj038.PubMedView ArticleGoogle Scholar
- Glynn E, Megee P, Yu H, Mistrot C, Unal E, Koshland D, DeRisi J, Gerton J: Genome-wide mapping of the cohesin complex in the yeast Saccharomyces cerevisiae. PLoS Biol. 2004, 2: E259-10.1371/journal.pbio.0020259.PubMedPubMed CentralView ArticleGoogle Scholar
- Kim T, Barrera L, Zheng M, Qu C, Singer M, Richmond T, Wu Y, Green R, Ren B: A high-resolution map of active promoters in the human genome. Nature. 2005, 436: 876-880. 10.1038/nature03877.PubMedPubMed CentralView ArticleGoogle Scholar
- Buck M, Nobel A, Lieb J: ChIPOTle: a user-friendly tool for the analysis of ChIP-chip data. Genome Biol. 2005, 6: R97-10.1186/gb-2005-6-11-r97.PubMedPubMed CentralView ArticleGoogle Scholar
- ChIPOTle. [https://sourceforge.net/projects/chipotle-perl/]
- Chuang P, Albertson D, Meyer B: DPY-27: A chromosome condensation protein homolog that regulates C. elegans dosage compensation through association with the X chromosome. Cell. 1994, 79: 459-474. 10.1016/0092-8674(94)90255-0.PubMedView ArticleGoogle Scholar
- Davis T, Meyer B: SDC-3 coordinates the assembly of a dosage compensation complex on the nematode X chromosome. Development. 1997, 124: 1019-1031.PubMedGoogle Scholar
- Ercan S, Giresi C, Whittle PG, Zhang X, Green R, Lieb J: X-chromosome repression by localization of the C. elegans dosage compensation machinery to sites of transcription initiation. Nat Genet. 2007, 39: 403-408. 10.1038/ng1983.PubMedPubMed CentralView ArticleGoogle Scholar
- Cheng J, Kapranov P, Drenkow J, Dike S, Brubaker S, Patel S, Long J, Stern D, Tammana H, Helt G, et al: Transcriptional maps of 10 human chromosomes at 5-nucleotide resolution. Science. 2005, 308: 1149-1154. 10.1126/science.1108625.PubMedView ArticleGoogle Scholar
- Schumacher A, Kapranov P, Kaminsky Z, Flanagan J, Assadzadeh A, Yau P, Virtanen C, Winegarden N, Cheng J, Gingeras T, et al: Microarray-based DNA methylation profiling: technology and applications. Nucleic Acids Res. 2006, 34: 528-542. 10.1093/nar/gkj461.PubMedPubMed CentralView ArticleGoogle Scholar
- Giresi P, Lieb J: How to find an opening (or lots of them). Nat Methods. 2006, 3: 501-502. 10.1038/nmeth0706-501.PubMedView ArticleGoogle Scholar
- Sabo P, Kuehn M, Thurman R, Johnson B, Johnson E, Cao H, Yu M, Rosenzweig E, Goldy J, Haydock A, et al: Genome-scale mapping of DNase I sensitivity in vivo using tiling DNA microarrays. Nat Methods. 2006, 3: 511-518. 10.1038/nmeth890.PubMedView ArticleGoogle Scholar
- Hubbell E, Liu W, Mei R: Robust estimators for expression analysis. Bioinformatics. 2002, 18: 1585-1592. 10.1093/bioinformatics/18.12.1585.PubMedView ArticleGoogle Scholar
- Emanuelsson O, Nagalakshmi U, Zheng D, Rozowsky J, Urban A, Du J, Lian Z, Stolc V, Weissman S, Snyder M, Gerstein M: Assessing the performance of different high-density tiling microarray strategies for mapping transcribed regions of the human genome. Genome Res. 2007, 17: 886-897. 10.1101/gr.5014606.PubMedPubMed CentralView ArticleGoogle Scholar
- Irizarry R, Hobbs B, Collin F, Beazer-Barclay Y, Antonellis K, Scherf U, Speed T: Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003, 4: 249-264. 10.1093/biostatistics/4.2.249.PubMedView ArticleGoogle Scholar
- Heintzman N, Stuart R, Hon G, Fu Y, Ching C, Hawkins R, Barrera L, Van Calcar S, Qu C, Ching K, et al: Distinct and predictive chromatin signatures of transcriptional promoters and enhancers in the human genome. Nat Genet. 2007, 39: 311-318. 10.1038/ng1966.PubMedView ArticleGoogle Scholar
- Affymetrix IGB. [http://www.affymetrix.com]
- MA2C. [http://liulab.dfci.harvard.edu/MA2C/MA2C.htm]
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/2.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.