Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2862456

Formats

Article sections

- Abstract
- 1. Background
- 2. Methods
- 3. Simulation studies
- 4. Case study: ZNF217 ChIP-chip data
- 5. Discussion
- Supplementary Material
- References

Authors

Related links

Pac Symp Biocomput. Author manuscript; available in PMC 2010 May 3.

Published in final edited form as:

Pac Symp Biocomput. 2008 : 515–526.

PMCID: PMC2862456

NIHMSID: NIHMS193956

See other articles in PMC that cite the published article.

Whole genome tiling arrays at a user specified resolution are becoming a versatile tool in genomics. Chromatin immunoprecipitation on microarrays (ChIP-chip) is a powerful application of these arrays. Although there is an increasing number of methods for analyzing ChIP-chip data, perhaps the most simple and commonly used one, due to its computational efficiency, is testing with a moving average statistic. Current moving average methods assume exchangeability of the measurements within an array. They are not tailored to deal with the issues due to array designs such as overlapping probes that result in correlated measurements. We investigate the correlation structure of data from such arrays and propose an extension of the moving average testing via a robust and rapid method called
`CMARRT`. We illustrate the pitfalls of ignoring the correlation structure in simulations and a case study. Our approach is implemented as an R package called
`CMARRT` and can be used with any tiling array platform.

Whole genome tiling arrays utilize array-based hybridization to scan the entire genome of an organism at a user specified resolution. Among their applications are ChIP-chip experiments for studying protein-DNA interactions. These experiments produce massive amounts of data and require rapid and robust analysis methods. Some of the commonly used methods are ChIPOTle,^{1} Mpeak,^{2} TileMap,^{3} HMMTiling,^{4} MAT^{5} and TileHGMM.^{6} Although these algorithms have been shown to be useful, they don't address the issues due to array designs. The most obvious issue is the correlation of the measurements from probes mapping to consecutive genomic locations.^{15} The basis for such a correlation structure is due to both overlapping probe design and fragmentation of the DNA sample to be hybridized on the array. There are several hidden Markov model (HMM) approaches to address the dependence among probes but the current implementations are limited to first order Markov dependence.^{4} Generalizations to higher orders increase the computational complexity immensely. We investigate the correlation structure of data from complex tiling array designs and propose an extension of the moving average approaches^{1}^{,}^{7} that carefully addresses the correlation structure. Our approach is based on estimating the variance of the moving average statistic by a detailed examination of the correlation structure and is applicable with any array platform. We illustrate the pitfalls of ignoring the correlation structure and provide several simulations and a case study illustrating the power of our approach
`CMARRT` (Correlation, Moving Average, Robust and Rapid method on Tiling array).

Let *Y*_{1},…,*Y _{N}* denote measurements on the

A common test statistic for analyzing ChIP-chip data is a moving average of *Y _{i}*'s over a fixed number of probes or fixed genomic distance.

$${T}_{i}=\frac{1}{2{w}_{i}+1}\sum _{j=i-{w}_{i}}^{i+{w}_{i}}{Y}_{i}.$$

(1)

Then, standard variance calculation leads to

$$var({T}_{i})=\frac{1}{{(2{w}_{i}+1)}^{2}}\left((2{w}_{i}+1){\sigma}^{2}+\sum _{j=i-{w}_{i}}^{i+{w}_{i}}\sum _{k\ne j}cov({Y}_{j},{Y}_{k})\right).$$

(2)

The standardized moving average statistic is given by

$${S}_{i}=\frac{{T}_{i}}{\sqrt{var({T}_{i})}}.$$

(3)

Standard practice of using moving average statistics relies on (1) estimating *σ*^{2} based on the observations that represent lower half of the unbound distribution; (2) ignoring the covariance term in equation (2); (3) and obtaining a null distribution under the hypothesis of no binding at probe *i*. In particular, ChIPOTle considers a permutation scheme where the probes are shuffled and the empirical distribution of the test statistic over several shufflings is used as an estimate of the null distribution. As an alternative, a Gaussian approximation is utilized assuming that *Y _{i}*'s are independent and identically distributed as normal random variables under the null distribution. As discussed by the authors of ChIPOTle, both approaches assume the exchangeability of the probes under the null hypothesis. Exchangeability implies that the correlation within any subset of the probes is the same. However, empirical autocorrelation plots from tiling arrays often exhibit evidence against this (Fig. 1). In particular, in the case of overlapping designs, a correlation structure is expected by design. When the spacing among the probes is large, correlation diminishes as expected (the right panel of Fig. 1), and this was the case for the dataset on which ChIPOTle was developed.

We illustrate the problem with ignoring the correlation structure on a ChIP-chip dataset from an E-coli RNA Polymerase II experiment utilizing a Nimblegen isothermal array (Landick Lab, Department of Bacteriology, UW-Madison). The probe lengths vary between 45 and 71 bp, tiled at a 22 bp resolution. Approximately half of the probes are of length 45 bp. We compute the standardized moving average statistic *S _{i}* (assuming cov(

Although it is desirable to develop a structured statistical model that captures the correlations, developing such a model is both theoretically and computationally challenging due to the complex, heterogeneous data generated by tiling array experiments. We propose a fast empirical method that estimates the correlation structure based on sample autocorrelation function. The covariance cov(*Y _{j}, Y_{j+k}*) can be estimated from the sample autocorrelation (

$$\widehat{\rho}(k)=\frac{{\sum}_{t=1}^{T-k}({Y}_{t}-\overline{Y})({Y}_{t+k}-\overline{Y})}{{\sum}_{t=1}^{T}{({Y}_{t}-\overline{Y})}^{2}},\widehat{cov}({Y}_{j},{Y}_{j+k})=\widehat{\rho}(k){\widehat{\sigma}}^{2}.$$

(4)

The following strategy is used in
`CMARRT` for estimating the correlation structure. The top *M%* of outlying probes which roughly correspond to bound probes are excluded in the estimation of (*k*). For the remaining probes, the sample autocorrelation at lag *k _{j}*(

In this section, we investigate the performance of
`CMARRT`, the conventional normal approximation approach under the independence assumption (Indep) and the HMM option in TileMap under various scenarios where we know the true bound regions in terms of sensitivity and specificity while controlling FDR at various levels used in practice.

We consider the following model

$${Y}_{i}={N}_{i}+{R}_{i},\phantom{\rule{1em}{0ex}}{N}_{i}=\sum _{k=1}^{p}{\alpha}_{i-k}{N}_{i-k}+{\epsilon}_{i},$$

(5)

where *N _{i}* is the autoregressive background component and

In this scenario, the data is simulated from hidden markov models (HMMs)^{12} with explicit state duration distribution to introduce direct dependencies at the probe level observations. Let the duration HMM densities be *p _{Si}* (

Each simulation scenario is repeated 50 times. A probe is declared as bound if its adjusted p-value^{11} is smaller than a pre-specified FDR level *α* when analyzing with
`CMARRT` and Indep. For TileMap, we use the direct posterior probability approach^{13} to control the FDR.

In Fig. 3, we summarize the sensitivity at the peak level and the specificity at various FDR thresholds from Simulation I for
`CMARRT`, Indep, and TileMap.
`CMARRT` is able to identify most of the bound regions at FDR of 0.05 and above while TileMap tends to be more conservative in declaring bound regions as shown in the sensitivity plots. Although Indep has the highest sensitivity, it also has a high proportion of false positives. The specificity of Indep is significantly lower compared to
`CMARRT`, even under the case of low correlation among the probes. Similar results are obtained in Simulation II under the duration HMM (Fig. 4). The left panels show the sensitivity and specificity for the case of smaller peaks with an average peak size of 10 probes while the right panels are for the case of larger peaks of size 20 probes on average. These results illustrate the superior performance of
`CMARRT` in terms of both sensitivity and specificity even when the data is generated from a complex model. The heuristic way of estimating the correlation structure in
`CMARRT` is able to reduce the false positives (specificity) significantly, but not at the expense of increasing false negatives (sensitivity). On the other hand, ignoring the correlation structure results in a higher proportion of false positives. Additionally, the HMM option in TileMap is more conservative than the moving average approach when the FDR is controlled at the same level.

We provide an illustration of
`CMARRT` with a ZNF217 ChIP-chip data tiling the ENCODE regions (available from Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/)^{14} with accession number GSE6624). The ENCODE regions were tiled at a density of one 50-mer every 38 bp, leading to ~ 380,000 50-mer probes on the array. We analyze two different replicates of this dataset separately and compare the analysis on these single replicates. In Krig et al.,^{14} the bound regions were identified with the Tamalpais Peaks program,^{9} which requires a bound region to have at least 6 consecutive probes in the top 2% of the log base 2 ratios. This criteria tends to be too stringent and fails to identify bound regions which contain a few outlier probes with log base 2 ratios below the top 2% threshold and may result in a higher level of false negatives. In the top right panel of Fig. 5, we show one potential peak missed by the Tamalpais Peaks program. In such cases, the sliding window approach is more powerful for finding peaks. Moreover, this method also assumes the observations are independent. As evident in the left panel of Fig. 1, observations from nearby probes in this tiling array are correlated. As shown in Fig 5, the histograms of p-values for the unbound probes under the independence assumption deviates from the expected distribution in both replicates. Similar problem is present in the normal quantile-quantile plots (online supp. mat.) when the correlation structure is ignored.

As in Krig et al.,^{14} we require the number of consecutive probes in each bound region to be at least 6. A set of peaks is obtained for each replicate at a given FDR control. We assess the extent of overlaps between the set of peaks in these two replicates. The results are summarized in Table 1. All the methods identified more peaks in replicate 1 than replicate 2. Therefore, using the peaks from rep 1 as reference, the common peaks are defined as the percentage of overlapping peaks in replicate 2. For all FDR thresholds (except 0.01),
`CMARRT` has the highest value of common peaks, followed by Indep and TileMap, which illustrates the consistency of the peaks identified by
`CMARRT`.

As an independent validation, we determine the location of bound regions relative to the transcription start site (TSS) of the nearest gene using GENECODE genes from UCSC Genome Browser as in Krig et al.^{14} (Table 1). For a given FDR control, the percentage of peaks located within ±2*kb*, ±10*kb* and ±100*kb* of the TSS is the highest in
`CMARRT`, followed by Indep and TileMap. As expected, these numbers decrease as we increase the FDR threshold for all the three methods. These results illustrate the power of
`CMARRT` in detecting biologically more plausible bound regions of ZNF217.

We have investigated and illustrated the pitfalls of ignoring the correlation structure due to tiling array design in ChIP-chip data analysis. We proposed an extension of the moving average approaches in
`CMARRT` to address this issue.
`CMARRT` is a robust and fast algorithm that can be used with any tiling platform and any number of replicates. Both the simulation results and the case study illustrate that
`CMARRT` is able to reduce false positives significantly but not at the expense of increasing false negatives, thereby giving a more confident set of peaks. We have recently became aware of the work of Bourgon^{15} who carefully studies the correlation structure in ChIP-chip arrays and proposes a fixed order autoregressive moving average model (ARMA(1, 1)) and we are in the process of comparing
`CMARRT` with this approach.

`CMARRT` is developed using the Gaussian approximation approach and the diagnostic plots illustrated can be utilized to detect whether a given dataset violates this assumption. One possible relaxation of this assumption is a constrained permutation approach that aims to conserve the correlation structure among the probes under the null distribution. Implementation of such an approach efficiently is a challenging future research direction.

We thank Professor Robert Landick for providing the E-coli ChIP-chip data for our analysis. Supplementary materials are available at http://www.stat.wisc.edu/~keles/CMARRT.sm.pdf. This research has been supported in part by a PhARMA Foundation Research Starter Grant (P.K. and S.K.) and NIH grants 1-R01-HG03747-01 (S.K.) and 4-R37-GM038660-20 (H.C.).

1. Buck MJ, Nobel AB, Lieb JD. ChIPOTle: a user-friendly tool for the analysis of ChIP-chip data. Genome Biol. 2005;6(11) [PMC free article] [PubMed]

2. Kim TH, Barrera LO, Zheng M, Qu C, Singer MA, Richmand TA, Wu Y, Green RD, Ren B. A high-resolution map of active promoters in the human genome. Nature. 2005;436:876–880. [PMC free article] [PubMed]

3. Ji H, Wong WH. TileMap: create chromosomal map of tiling array hybridizations. Bioinformatics. 2005;21(18):3629–3636. [PubMed]

4. Li W, Meyer CA, Liu XS. A hidden Markov model for analyzing ChIP-chip experiments on genome tiling arrays and its application to p53 binding sequences. Bioinformatics. 2005;21(1):i274–i282. [PubMed]

5. Johnson WE, Li W, Meyer CA, Gottardo R, Carroll JS, Brown M, Liu XS. MAT: Model-based Analysis of Tiling-arrays for ChIP-chip. Proc Natl Acad Sci USA. 2006;103:12457–12462. [PubMed]

6. Keles S. Mixture modeling for genome-wide localization of transcription factors. Biometrics. 2006;63(1):10–21. [PubMed]

7. Keles S, van der Laan MJ, Dudoit S, Cawley SE. Multiple Testing Methods for ChIP-Chip High Density Oligonucleotide Array Data. J of Comp Bio. 2006;13(3):579–613. [PubMed]

8. Royce TE, Rozowsky JS, Gerstein MB. Assessing the need for sequence-based normalization in tiling microarray experiments. Bioinformatics 2007 [PubMed]

9. Bieda M, Xu X, Singer MA, Green R, Farnham PJ. Unbiased location analysis of E2F1-binding sites suggests a widespread role for E2F1 in the human genome. Genome 2007 [PubMed]

10. Box GP, Jenkins GM. Time series analysis forecasting and control. Holden-Day; 1976.

11. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. JRSS-B. 1995;57:289–300.

12. Rabiner LR. A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition. Proceedings of the IEEE. 1989;77(2):257–286.

13. Newton MA, Noueiry A, Sarkar D, Ahlquist P. Detecting differential gene expression with a semiparametric hierarchical mixture method. Biostatistics. 2004;5:155–176. [PubMed]

14. Krig SR, Jin VX, Bieda MC, O'Geen H, Yaswen P, Green R, Farnham PJ. Identification of genes directly regulated by the oncogene ZNF217 using ChIP-chip assays. J Biol Chem. 2007;282(13):9703–9712. [PMC free article] [PubMed]

15. Bourgon RW. Ph D Thesis. UC Berkeley; 2006. Chromatin immunoprecipitation and high-density tiling microarrays: a generative model, methods for analysis and methodology assessment in the absence of a “gold standard”

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |