ChIPmix: mixture model of regressions for two-color ChIP–chip analysis
1UMR AgroParisTech/INRA MIA 518, 16 rue Claude Bernard, 75231 Paris Cedex 05 and 2URGV UMR INRA/CNRS/UEVE, 2 rue Gaston Crémieux, CP5708, 91057, Evry Cedex, France
*To whom correspondence should be addressed.
 |
ABSTRACT |
|---|
 TOP ABSTRACT 1 INTRODUCTION2 Statistical framework3 RESULTS4 DISCUSSIONAppendixACKNOWLEDGEMENTSREFERENCES |
|---|
Motivation: Chromatin immunoprecipitation (ChIP) combined with DNA microarray is a high-throughput technology to investigate DNA–protein binding or chromatin/histone modifications. ChIP–chip data require adapted statistical method in order to identify enriched regions. All methods already proposed are based on the analysis of the log ratio (Ip/Input). Nevertheless, the assumption that the log ratio is a pertinent quantity to assess the probe status is not always verified and it leads to a poor data interpretation.
Results: Instead of working on the log ratio, we directly work with the Ip and Input signals of each probe by modeling the distribution of the Ip signal conditional to the Input signal. We propose a method named ChIPmix based on a linear regression mixture model to identify actual binding targets of the protein under study. Moreover, we are able to control the proportion of false positives. The efficiency of ChIPmix is illustrated on several datasets obtained from different organisms and hybridized either on tiling or promoter arrays. This validation shows that ChIPmix is convenient for any two-color array whatever its density and provides promising results.
Availability: The ChIPmix method is implemented in R and is available at http://www.agroparistech.fr/mia/outil_A.html
Contact: marie_laure.martin{at}agroparistech.fr
|
1 INTRODUCTION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 Statistical framework3 RESULTS4 DISCUSSIONAppendixACKNOWLEDGEMENTSREFERENCES |
|---|
Chromatin immunoprecipitation (ChIP) is a well-established procedure used to investigate proteins associated with DNA. ChIP on chip involves analysis of DNA recovered from ChIP experiments by hybridization to microarray. In a two-color ChIP–chip experiment, two samples are compared: DNA fragments crosslinked to a protein of interest (IP) and genomic DNA (Input). The two samples are differentially labeled and then co-hybridized on a single array. The goal is then to identify actual binding targets of the IP, i.e. probes whose IP signal is significantly larger than the Input signal.
Many authors have already pointed out the need for efficient statistical procedures to detect enriched probes (Buck and Lieb, 2004; Keles, 2007). Recently, two strategies have been widely applied for the detection of enriched DNA regions. The first strategy takes advantage of the spatial structure of the data. Since probes are positioned all along the genome, if one region is enriched we expect several adjacent probes to obtain high ratio measurements, resulting in a ‘peak’ of intensity. Spatial methods such as sliding windows (Cawley et al., 2004; Keles, 2007) or Hidden Markov Models (Ji and Wong, 2005; Li et al., 2005) have been proposed to detect these peaks. Alternatively, the second strategy is to consider that the whole population of probes can be divided into two components: the population of IP-enriched genomic fragments, and the population of genomic DNA that is not IP enriched. Different statistical methods have been proposed to distinguish between the two populations by considering the distribution of the ratios (or their associated rank). Assuming that a non-negligible proportion of the fragments are enriched, the log ratio distribution is bimodal, the highest mode corresponding to the enriched population. A probe is then declared enriched when its ratio exceeds a selected cutoff, which is fixed according to the data distribution (Buck and Lieb, 2004).
Importantly, both strategies assume that the log ratio measurement is a pertinent statistical quantity to assess the probe status (enriched or not). This assumption is correct if the distribution of the ratio mostly depends on the status (normal/enriched) of the probe. Figure 1A shows the ideal situation described in Buck and Lieb (2004), where the distribution is bimodal. In many applications, the distribution of the log ratios is closer to Figure 1B, and the performance of log ratio-based methods may be poor. At least two technical reasons may explain the difference between the ideal and real cases. First, there are some technical difficulties to obtain the IP sample: it requests the use of a very specific antibody and a careful experimental process to avoid a high level of contamination. The second reason comes from the possible cross-hybridization phenomena.
From observation of Figure 1, we argue that it is worth working directly with the two measurements of each probe (Input and IP) rather than with the log ratio. In Figure 1C, we observe that the relationship between the two measurements is almost linear. Working on log ratio amounts to stating that the slope of the linear relationship is the same whatever the status of the probe. In many cases the slopes are different: Figure 2 (synthetic data) shows that even a slight difference between the two slopes may turn the distribution of the log ratios into unimodal rather than bimodal, as observed for the NimbleGen slide in Figure 1.
In this work, we propose a new statistical method that we call ChIPmix, based on a mixture model of regressions. This framework allows us to well characterize the IP–Input relationship, and to provide a statistical procedure to control the proportion of probes wrongly classified as enriched. The article is organized as follows. The statistical model and the procedure for false positive control are described in Section 2. In Section 3, we consider several large datasets obtained from different organisms and hybridized on different array types (tiling or promoter). We show that the method outperforms competing methods in terms of sensitivity. The main conclusions and some possible extensions are discussed in Section 4.
|
|
|
2 Statistical framework |
|---|
|
TOPABSTRACT1 INTRODUCTION2 Statistical framework3 RESULTS4 DISCUSSIONAppendixACKNOWLEDGEMENTSREFERENCES |
|---|
2.1 Model and inference
Let (xi, Yi) be the log-Input and log-IP intensities of probe i, respectively. The (unknown) status of the probe is characterized through a label Zi which is 1 if the probe is enriched and 0 if it is normal (not enriched). We assume the Input–IP relationship to be linear whatever the population, but with different slope and intercept. More precisely, we have:
|
|
i is a Gaussian random variable with mean 0 and variance 
2. Such a model is named a mixture model of regressions.
The marginal distribution of Yi for a given level of Input xi is
|
| (1) |
is the proportion of enriched probes, and 
j(·|x) stands for the probability density function (PDF) of a Gaussian distribution with mean aj+bjx and variance 2.
The mixture model is used to classify probes as normal or enriched. To do this, we calculate the probability of a probe to be enriched given its Input and IP intensities. This probability is called the posterior probability and is defined from Equation (1) by
|
| (2) |
The mixture parameters (proportion, intercepts, slopes and variance) are estimated using the EM algorithm. The EM algorithm is dedicated to the class of incomplete data models where the status of the observations is unknown. In the E step, the posterior probability for each observation to belong to each class is calculated. In the M step, the parameters of each class are estimated using a weighted regression, in which the weights are given by the posterior probabilities. This algorithm is implemented in the mixreg function of the mixreg R package (Turner, 2000). Figure 4A shows the application of Chipmix on the NimbleGen high-density array data, presented in Section 3.3.
In the mixreg function, the initial values of the parameters must be given by the users otherwise they are chosen randomly. Nevertheless, the EM algorithm is well-known to be sensitive to the initial values (Bohning and Seidel, 2003; Karlis and Xekalaki, 2003) and to solve this difficulty, we propose initial values derived from the first axis of the Principal Component Analysis (PCA) of the whole dataset (see the ChIPmix R function for details).
The mixture model with two linear regressions is adapted if the protein under study has some targets. When the protein has no target, all probes belong to the normal class. In this case, a simple linear regression is sufficient to fit the data. For each dataset the two models (one or two classes) are fitted and the best model is selected according to the BIC criterion (Schwarz, 1978).
2.2 False discovery control
Posterior probabilities are used to classify probes into the normal or enriched class, using the following classification rule
|
|
In ChIP–chip experiments, where false positives are of concern, it is important to control the false positive proportion and to fix s accordingly. In the hypothesis test theory, the false discovery control is performed by controlling the probability to reject wrongly the null hypothesis. We propose an analogous concept in the mixture model framework. Our aim is to control the probability for a probe to be wrongly assigned to the enriched class. Therefore, we want Pr{
i > s|xi, Zi=0} to be equal to a predefined level 
. In practice, we fix and we find the threshold s depending on and xi (see Appendix).
|
3 RESULTS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 Statistical framework3 RESULTS4 DISCUSSIONAppendixACKNOWLEDGEMENTSREFERENCES |
|---|
We present three applications of ChIPmix to assess the performance of the method whatever the specificity and density of the array (tiling or promoter array). The first two applications validate the method on already published data. The third dataset is used to compare ChIPmix with existing methods.
3.1 Promoter DNA methylation in the human genome
Weber et al. (2007) measured DNA methylation using a NimbleGen microarray representing 15 609 promoter regions of the human genome. Each promoter region is covered by 15 probes and is classified into a category according to its CpG rate. We focus on the analysis of the class ICP (intermediate CpG promoter). Weber et al. based their classification on the mean log ratio value for the 15 probes per promoter region. If this value was larger than 0.4 (threshold based on bisulfite sequencing), the promoter region was declared hypermethylated. Among the 2056 promoter regions under study, 460 were declared hypermethylated.
We applied ChIPmix to these data without averaging the 15 values per promoter region. The estimated proportion of enriched probes was 0.794. This is in keeping with a large proportion of targets expected in such experiments. The estimated regression slopes were 
for the normal class and 
for the enriched one, which shows that the Input–Ip relations substantially differ between the two status. At the level =0.01, a total of 1706 promoter regions were found to have at least one probe enriched. Except for one region, all the promoter regions of the Weber's list have at least enriched probe, and 403 have 5 or more enriched probes. Besides, ChIPmix identified 38 promoter regions with 9 probes or more classified as enriched that were not detected in Weber et al. (2007).
3.2 Histone modification in Arabidopsis thaliana
Turck et al. (2007) studied several histone modifications of A.thaliana using a custom genomic tiling array of Chromosome 4. To declare a tile enriched, they developed a two-step method based on a Gaussian mixture model and a total of 2775 tiles were found to be marked by histone H3 trimethylated at lysine 27 (H3K27me3) according to their analysis.
We analyzed the same dataset using ChIPmix. The estimated proportion and slopes were 
, 
and 
. The tiles classified as enriched at risk =0.01 include all the tiles found by Turck et al. (2007) plus 2346 others: 1404 tiles extend the genomic region already found marked by H3K27me3 and 942 tiles form 62 new genomic regions. The difference between the two slopes enables us to better discriminate the two classes for high Input intensities. This may explain the higher number of enriched probes detected by ChIPmix.
3.3 NimbleGen high-density array (Histone modification H3K9me3)
In this last example, we considered Chip–chip data produced on a two-color NimbleGen array of 1 132 140 probes. Each chromosome of the model plant A.thaliana is covered by about 200 000 probes. Such very high-density arrays are more and more popular, so we need to assess the efficiency of ChIPmix on such a very large dataset.
|
From a biological point of view, the same IP and Input samples were already hybridized on a custom genomic tilling array covering the Chromosome 4 (Turck et al., 2007). Regions identified in Turck et al. (2007) were biologically validated and are used as true positives. In addition, the chloroplastic genome can be used as a negative control, since no histone modification target is expected in this region. We did not use the mitochondrial genome as a negative control since some regions have been duplicated in the nuclear genome. From a statistical point of view, since ChIPmix does not take the spatial structure into account, it is important to compare it with methods using this information. We compared our results with those provided by the NimbleGen software and ChIPOTle method (Buck et al., 2005). NimbleGen software uses a permutation-based algorithm to find statistically significant peaks, using scaled log ratio data, and ChIPOTle method uses a sliding window approach.
Two biological replicates were available, for which hybridizations were performed in dye-swap. We performed a normalization step to remove technical biases as well as dye bias. Since the Input and IP samples differed substantially, array-by-array normalization such as lowess could not be applied. We quantified biases by an ANOVA model (Kerr et al., 2002), and removed them from the raw data. The IP and Input signals for each biological replicate were averaged on the dye-swap to remove the gene-specific dye bias. Analyses per chromosome were performed on the normalized data.
For a risk =0.01, a total of 30 477 probes were detected in the first replicate and 27 553 in the second. The intersection contains more than two-thirds of the probes declared enriched in at least one replicate (23 546 probes). Although ChIPmix does not take the spatial structure of the genome into account, enriched probes are clustered in genomic regions (Fig. 3). These regions are rich in genes and corroborate the results of Turck et al. (2007), who have shown that H3K9me3 is actually a euchromatin mark. Moreover, more than 80% of the probes classified as enriched in this experiment cover genomic regions already found in Turck et al. (2007).
For the chloroplastic genome, the BIC criterion selected a two component regression model. For the first biological replicate, two 2-probe clusters and one 3-probe cluster were declared enriched. On the same replicate, ChIPOTle (window = 500 and step = 100) found five 2-probe clusters, two 3-probe cluster and one 6-probe cluster. With other parameters (window = 200 and step = 50), the number of detected peaks increased. Results are similar for the second biological replicate, and one cluster was declared enriched with ChIPmix in both biological replicates (two with ChIPOTle). NimbleGen did not provide the analysis of the chloroplastic genome.
We also compared ChIPmix to the results given by NimbleGen and ChIPOTle on Chromosome 4, studied in Turck et al. (2007). The probes declared enriched by ChIPmix include almost all of those found enriched by the NimbleGen software, but cover much larger genomic regions. Moreover ChIPmix identifies other genomic regions not found by the NimbleGen software (Fig. 3), that are validated by a comparison with results of Turck et al. (2007). ChIPmix detects 30 477 enriched probes, including 24 575 in common with Turck et al. (2007). ChIPOTle detects 24 357 probes [20 866 common with Turck et al. (2007)] and NimbleGen detects 19 837 probes [16 600 common with Turck et al. (2007)] (Fig. 4B). Among the three methods ChIPmix provides the closest results to the reference publication.
|
. Colors change every 20% (blue: 
, red: 
). The two lines are the two estimated linear regressions of the mixture. (B) Venn diagram summarizing the results of the three methods. |
4 DISCUSSION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 Statistical framework3 RESULTS4 DISCUSSIONAppendixACKNOWLEDGEMENTSREFERENCES |
|---|
We propose a statistical method based on mixture of regression to classify probes in ChIP–chip experiments. Our approach accounts for different relations between IP and Input intensity in the two classes of probes (enriched and normal). The ChIPmix method outperforms the standard approaches based on the log ratio.
We presented various applications each dedicated to one specific biological question (histone modification and DNA methylation on different organisms). ChIPmix can also be applied to the detection of transcription factor binding sites (TFBS, results not shown). The method is valid when the proportion of positive probes is expected to be large (e.g. histone modification), or small (e.g. TFBS). Through the examples we have shown that ChIPmix is convenient for any two-color chip whatever its density (array size from thousands to hundreds of thousands of probes) and the nature of the probe (tiling and promoter arrays).
ChIPmix does not account for the spatial structure of the data. While this could be seen as a drawback, we showed that enriched probes are clustered into genomic regions in the presented applications. Moreover, this may become perfectly relevant for specific experiments as well as RIP-chip, which investigates interactions between protein and RNA (Schmitz-Linneweber et al., 2005) or ChIP–chip experiments performed on array where promoter are represented by only one probe (see project SAP at www.psb.ugent.be/SAP/).
The only parameter of the ChIPmix method is the risk , which can be easily interpreted. In contrast, two parameters have to be tuned in the ChIPOTle method (window size and step). The tuning of this two parameters depends on both the experimental protocol and the array type. The results are very sensitive to this tuning.
The proposed strategy can be extended in different ways. The ChIPmix extension to the unequal variance case is straightforward. However, the equal variance case provides an efficient framework for the false discovery control. If the equality of variance is not assumed, the calculations given in appendix do not hold anymore, and the solving of the equation system becomes much more complex.
The proposed regression models allow us to correct the IP intensity with respect to the Input one. Other elements may influence the level of IP signal. Weber et al. (2007) show that the CpG rate has to be taken into account to classify probes. The specificity of the probes (number of hits) may also alter the IP intensity. All this information can be considered as covariates and added in the model. This will lead to a mixture of multiple regression for which the statistical framework is almost the same as the one we propose.
The proportion of false negative results can be controlled in the same way as the false discovery described in Section 2.2. This allows us to evaluate the sensitivity of the classification at each Input level. Moreover, the two criteria (false negative and false discovery) can be combined to derive a threshold s that optimizes some trade-off between them.
|
Appendix |
|---|
|
TOPABSTRACT1 INTRODUCTION2 Statistical framework3 RESULTS4 DISCUSSIONAppendixACKNOWLEDGEMENTSREFERENCES |
|---|
We propose to control the probability for a normal probe i to be wrongly assigned to the enriched class:
|
| (3) |
and we find the threshold s depending on and xi. Using definition 2, Pr{i>s|xi, Zi=0} can be rewritten as |
|
0(Yi|xi) and 1(Yi|xi) with their expression, we get Equation (3) equivalent to
|
| (4) |
|
|
2, and we deduce that Equation (4) is equivalent to solve |
|
is the (1–) quantile of Gaussian with mean 0 and variance 1. Using the definition of 
(s, xi), the expression of threshold s is given by 
where |
|
|
ACKNOWLEDGEMENTS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 Statistical framework3 RESULTS4 DISCUSSIONAppendixACKNOWLEDGEMENTSREFERENCES |
|---|
The authors want to thank Vincent Colot, Alain Lecharny and Michel Caboche from the URGV unit for helpful discussions and advice.
Funding: This work was supported by the TAG ANR/Genoplante project.
|
FOOTNOTES |
|---|
The authors wish it to be known that, in their opinion, the first two authors should be regarded as joint First Authors. 
|
REFERENCES |
|---|
|
TOPABSTRACT1 INTRODUCTION2 Statistical framework3 RESULTS4 DISCUSSIONAppendixACKNOWLEDGEMENTSREFERENCES |
|---|
Bohning D, Seidel W. Editorial: recent developments in mixture models. Comput. Stat. Data Anal (2003) 41:349–357.[CrossRef]
Buck MJ, Lieb JD. Chip-chip: considerations for the design, analysis, and application of genome-wide chromatin immunoprecipitation experiments. Genomics (2004) 83:349–360.[CrossRef][Web of Science][Medline]
Buck MJ, et al. ChIPOTle: a user-friendly tool for the analysis of ChIP-chip data. Genome Biol (2005) 6:R97.[CrossRef][Medline]
Cawley S, et al. Unbiased mapping of transcription factor binding sites along human chromosomes 21 and 22 points to widespread regulation of noncoding RNAs. Cell (2004) 116:499–509.[CrossRef][Web of Science][Medline]
Ji H, Wong WH. Tilemap: create chromosomal map of tiling array hybridizations. Bioinformatics (2005) 21:3629–3636.
Karlis D, Xekalaki E. Choosing initial values for the EM algorithm for finite mixtures. Comput. Stat. Data Anal (2003) 41:577–590.[CrossRef]
Keles S. Mixture modeling for genome-wide localization of transcription factors. Biometrics (2007) 63:10–21.[CrossRef][Web of Science][Medline]
Kerr MK, et al. Statistical analysis of a gene expression microarray experiment with replication. Stat. Sin (2002) 12:203–217.
Li W, et al. A hidden markov model for analyzing chip-chip experiments on genome tiling arrays and its application to p53 binding sequences. Bioinformatics (2005) 21:i274–i282.[Abstract]
Schmitz-Linneweber C, et al. RNA immunoprecipitation and microarray analysis show a chloroplast pentatricopeptide repeat protein to be associated with the 5' region of mRNAs whose translation it activates. Plant Cell (2005) 17:2791–2804.
Weber M, et al. Distribution, silencing potential and evolutionary impact of promoter dna methylation in the human genome. Nat. Genet (2007) 39:457–466.[CrossRef][Web of Science][Medline]
Schwarz G. Estimating the dimension of a model. (1978) 6:461–464.
Turck F, et al. Arabidopsis tfl2/lhp1 specifically associates with genes marked by trimethylation of histone h3 lysine 27. PLoS Genet (2007) 3:e86.[CrossRef][Medline]
Turner TR. Estimating the rate of spread of a viral infection of potato plants via mixtures of regressions. Appl. Stat (2000) 49:371–384.
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  Y. Kim, S. Bekiranov, J. K. Lee, and T. Park Double error shrinkage method for identifying protein binding sites observed by tiling arrays with limited replication Bioinformatics, October 1, 2009; 25(19): 2486 - 2491. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||