Effect of false positive and false negative rates on inference of binding target conservation across different conditions and species from ChIP-chip data

Abstract

Background

ChIP-chip data are routinely used to identify transcription factor binding targets. However, the presence of false positives and false negatives in ChIP-chip data complicates and hinders analyses, especially when the binding targets for a specific transcription factor are compared across conditions or species.

Results

We propose an Expectation Maximization based approach to infer the underlying true counts of "positives" and "negatives" from the observed counts. Based on this approach, we study the effect of false positives and false negatives on inferences related to transcription regulation.

Conclusion

Our results indicate that if there is a significant degree of association among the binding targets across conditions/species (log odds ratio > 4), moderate values of false positive and false negative rates (0.005 and 0.4 respectively) would not change our inference qualitatively (i.e. the presence or absence of conservation) based on the observed experimental data despite a significant change in the observed counts. However, if the underlying association is marginal, with odds ratios close to 1, moderate to large values of false positive and false negative rates (0.01 and 0.2 respectively) could mask the underlying association.

Background

Transcription factors play an important role in gene regulation by binding to specific DNA sequences in the regulatory regions of their targets. Accurate identification of the binding targets of the transcription factors is paramount to the understanding of the regulatory mechanism. Chromatin immunoprecipitation (ChIP) experiments are commonly used to identify the regulatory targets in prokaryotes and eukaryotes. ChIP-chip experiments provide us with information about the binding targets of a particular regulator at the genome level [1–4].

The output of ChIP-chip experiments are often summarized in binary forms. Using replicate data, the statistical evidence for a gene being the binding target of a transcription factor is typically summarized as a p-value. A threshold for the p-value, e.g. 0.001 is then chosen, and genes with p-values less than the threshold are considered the binding targets for the transcription factor. Thus, for a transcription factor, we can enumerate a list of genes which are "positives", i.e. binding targets and a list of genes which are "negatives", i.e. non-binding targets. If the threshold is set at a very stringent level to control the number of false positives, this will be achieved at the expense of high false negatives. A more relaxed threshold will reduce the number of false negatives, but will end up with more false positive results. Over the past few years, ChIP-chip data has formed the basis of many transcription regulatory mechanism studies. Several groups have compared the binding of a regulator across multiple experimental conditions to determine condition dependence of binding [5–7]. Similarly, binding data of specific transcription factors across species has been used to investigate the presence of conserved binding targets [8]. Unfortunately, the presence of noise, in the form of false positives and false negatives as discussed above, may lead to inaccurate inference of the binding targets, and thus biased results and potentially incorrect conclusions on key aspects of transcription regulation, e.g. preservation of regulation targets across conditions and species. In this article, we develop a statistical approach to analyzing ChIP-chip data, appropriately incorporating false positives and false negatives. Based on our approach, we investigate the effect of false positives and false negatives on the inference of conservations of binding targets based on ChIP-chip data.

Methods

Summarizing Contingency Tables

As discussed above, the output of ChIP-chip experiments is typically summarized into binary forms and results across different experiments for the same transcription factor can be crosstabulated into a contingency table. A common question asked is whether a transcription factor has similar binding targets across conditions, and this is reflected as the dependency of outcome among the conditions. In the following, we give a brief discussion on two statistical measures that we will use to summarize the degree of dependency in a contingency table.

For the sake of clarity, we will focus on ChIP-chip experiments involving two different conditions or two species. The number of target genes in the two conditions/species can be cross tabulated into a 2 by 2 contingency table. We use two metrics to summarize such contingency tables – Odds Ratio and Positive specific agreement [9, 10].

Table 1 gives an illustration of a 2 by 2 contingency table. The goal is to identify whether a relationship, or association exists between the two categorical variables. In our scenario, it would correspond to whether the transcription factor exhibits condition dependent binding or condition independent binding. For such a contingency table, the odds ratio is a commonly used measure to quantify association among the categorical variables. An odds is defined as the ratio of the frequency of being in one category and the frequency of not being in that category. For example, from Table 1, the odds that a particular gene is a binding target in experimental condition 1 is equal to (c + d)/(a + b). This odds is called marginal odds, obtained from the total frequencies in one margin of the table, disregarding the effects of the other variable. Conditional odds are the chances of the transcription factor binding relative to not binding in one experimental condition, given a particular level (binding state) in the other experimental condition. The variables are deemed to be unassociated if the conditional odds are equal or close to each other, and hence equal to the marginal odds. To compare directly the two conditional odds, a single summary statistic, obtained by dividing the first conditional odds by the second is called odds ratio. Thus, for the data in Table 1, the odds ratio is defined as: Odds Ratio = (a/c)/(b/d) = ad/bc. Odds ratio takes only positive values and has no upper limit. An odds ratio of 1 indicates no relationship among the variables. In addition to the odds ratio, its logarithm is also commonly used. Logarithmic transformation of data has a number of advantages – the variation of log transformed data tends to be less dependent on the magnitude of values, while taking logs also reduces the skewness of the distributions. After log transformation, data tends to be spread out more evenly, also making it easier to examine visually.

Table 1 Simple 2 by 2 contingency table.

Other measures of dependency are also often used in psychological and medical research. For example, the problem can be formulated as follows: Suppose two raters classify each subject in a sample from some target population according to the presence or absence of some characteristic of interest. The resulting data can then be summarized into a 2 by 2 table. The agreement between raters can be quantified by the metric simple agreement, which is defined as the proportion of cases for which both raters agree, or (a + d)/(a + b + c + d). However, if a is large, this would approach 1 regardless of the performance on positive cases. Positive specific agreement provides insight when the positive cases are rare. It estimates the conditional probability that one rater will agree that a case is positive given the other one rated it positive, where the role of the two raters is selected randomly. Positive specific agreement, p _pos is defined as: p _pos = 2d/(2d + b + c). Both (log) odds ratio and positive specific agreement will be considered in our following discussion.

Model Setup

Consider an experiment with a binary outcome. Let p₀ denote the proportion of true negatives, while p₁ be the proportion of true positives. We denote p = (p₀, p₁)^tas the vector of true proportions. Due to false positives and false negatives, the observed proportions likely differ from the true proportions. Let $\widehat{p}={({\widehat{p}}_0,{\widehat{p}}_1)}^t$ denote the vector of the observed proportions. The relationship between p and E($\widehat{p}$) can be written as:

$$\left(\begin{array}{c}E({\widehat{p}}_0)\\ E({\widehat{p}}_1)\end{array}\right)=\left(\begin{array}{cc}1-s& t\\ s& 1-t\end{array}\right)\left(\begin{array}{c}{p}_0\\ p_1\end{array}\right),$$

(1)

where s is the false positive rate and t is the false negative rate. Denoting the transformation matrix as M, Equation (1) can be written as:

$$E(\widehat{p})=Mp.$$

(2)

Thus, for different values of false positive and false negative rates, different observed proportions will be obtained based on Equation (2). If the false positive and false negative rates are known, the true proportions may be inferred based on the observed experimental proportions. Multiplying both sides of Equation (2) by M^-1 gives us:

$$M^{-1}E(\widehat{p})=p.$$

(3)

However, due to chance variations, p obtained through this approach based on the observed $\widehat{p}$ may have negative components, leading to uninterpretable results. Instead, we propose to estimate the true proportions using an Expectation Maximization (EM) based approach explained in detail in the following subsection.

Often, we are interested in the analysis of the binding of a particular transcription factor in multiple experimental conditions or across different species. In either case, we are interested in counts of similarity of binding across conditions or organisms. This would correspond to an extension of Equation (2) into a higher dimension. For simplicity, we present our analysis for a 2-dimensional case. For example, if we consider the binding targets of a transcription factor across two experimental conditions, the vector of true proportions can be represented as p = (p₀₀, p₀₁, p₁₀, p₁₁)^t. Here p₀₀ denotes the proportion of genes which are not targets of the regulator in either condition, p₀₁ denotes the proportion of genes which are targets of the regulator in the second condition but not in the first, p₁₀ denotes the proportion of genes which are targets of the regulator in the first condition but not in the second, while p₁₁ denotes the proportion of genes which are targets of the regulator in both conditions. Similarly, the vector for the observed proportions can be denoted as $\widehat{p}={({\widehat{p}}_{00},{\widehat{p}}_{01},{\widehat{p}}_{10},{\widehat{p}}_{11})}^t$. The relationship between the observed and true proportions can be then written as:

$$E(\widehat{p})=(M\otimes M)p.$$

(4)

If we consider Equation (2) to correspond to a 1-dimensional case, for the n-dimensional case, the new transformation matrix would simply be obtained by taking the tensor product of M with itself n times. Here we assume that the false positive and false negative rates to be the same across two conditions. In general that may not be the case. In such a scenario, for a 2-dimensional case, Equation 4 takes the general form:

$$E(\widehat{p})=(M_1\otimes M_2)p,$$

(5)

where M₁ and M₂ are the transformation matrices for the first and second conditions respectively.

EM Algorithm

Given a vector of observed proportions which we obtain from experimental output, for different values of false positive rates and false negative rates, we aim to infer the true proportions. This would give us an idea about how the observed and true proportions differ for different levels of noise in the form of false positives and false negatives. We infer the true proportions from the observed proportions using an EM based approach which we now discuss in detail.

Let us consider the binding patterns for a transcription factor in experimental conditions c₁ and c₂. We define the vector for the true binary binding pattern of a particular gene G as b = (b₁, b₂), where b₁ and b₂ take binary value 1 or 0 depending on whether the gene is a true binding target for the transcription factor in c₁ and c₂ respectively. Thus, for the experimental conditions c₁ and c₂, this binary binding pattern vector can take four possible values, {(0, 0), (0, 1), (1, 0), (1, 1)}. For example, a binary binding pattern vector equal to (1, 1) indicates that the gene is a binding target for the transcription factor in both c₁ and c₂. We aim to infer this true binary binding pattern for all the genes and thus obtain the true binary counts. Due to experimental errors, we have the observed counts as the experimental output.

We denote the observed binding pattern for a particular gene as g = (g₁, g₂), where each component is either 0 or 1 denoting whether the gene is observed to be the binding target of the particular transcription factor in c₁ and c₂ respectively, based on the experimental output. Thus, the vector g represents the observed data. The probability of the observed binding data is then given by

$$\begin{array}{lll}P(g)\hfill & =\hfill & P(b=(0,0))P(g|b=(0,0))+P(b=(0,1))P(g|b=(0,1))+\hfill \\ P(b=(1,0))P(g|b=(1,0))+P(b=(1,1))P(g|b=(1,1)).\hfill \end{array}$$

(6)

Thus, for N genes, the probability of the observed data is

$$P(g_1,g_2,\mathrm{...},g_N)={\displaystyle \prod _{i=1}^{N}P(g_i)}.$$

(7)

In this article, we propose to estimate the P(b) to maximize P(g₁, g₂, ..., g _N ) using the EM algorithm, by treating b as the missing data as follows.

E-Step: In the Expectation step, the conditional distribution of the missing data given the observed data is evaluated. We evaluate the posterior probabilities of each true binding state given the observed binding pattern. Thus, for every gene G with observed binding pattern g, we estimate:

$$P(b^{(m)}|g)=\frac{P(g|b^{(m)})\ast P(b^{(m)})}{{\displaystyle \sum P(g|b^{(m)})\ast P(b^{(m)})}}$$

(8)

where b^(m)is the estimate of the true binding state b probability at the m-th step. Since b^(m)can have four possible values, at each step, we estimate four probabilities. The probability P(g|b^(m)) can be expanded as:

$$\begin{array}{lll}P(g|b^{(m)})\hfill & =\hfill & P((g_1,g_2)|(b_1^{(m)},b_2^{(m)}))\hfill \\ =\hfill & P(g_1|b_1^{(m)})P(g_2|b_2^{(m)}),\hfill \end{array}$$

(9)

where $b_1^{(m)}$ and $b_2^{(m)}$ are the first and second components of the estimate b^(m)and take binary value 0 or 1.

The second equation in (9) results from the independence assumption for the data from two separate ChIP-chip experiments. Thus, the probability of observing g₁ would be independent of the estimate of binding state $b_2^{(m)}$, while the probability of observing g₂ would be independent of the estimate of binding state $b_1^{(m)}$. There are four possible cases for the expression $P(g_i|b_i^{(m)})$ in Equation (4). From Equation (1), they can be enumerated as:

$$\begin{array}{ll}P(g_i=0|b_i^{(m)}=0)=1-s,\hfill & P(g_i=0|b_i^{(m)}=1)=t,\hfill \\ P(g_i=1|b_i^{(m)}=0)=s,\hfill & P(g_i=1|b_i^{(m)}=1)=1-t.\hfill \end{array}$$

(10)

Thus, for each gene, we start with a set of estimates P(b^(m)) and obtain estimates of the posterior probabilities P(b^(m)|g) for each gene at the E-step.

M-Step: In the Maximization step, the parameters P(b) are re-estimated to maximize the likelihood of the complete data. After obtaining P(b^(m)|g) for each gene, we cross-tabulate a two-way contingency table, with the "count" for each of the four values {(0, 0), (0, 1), (1, 0), (1, 1)} being the sum of the probabilities for that particular value across all the genes. These counts are then used to obtain updated P estimates for P(b). For example, $P(b=(0,0))=\frac{1}{N}{\displaystyle \sum P(b=0,0)|g_i)}$.

We iterate between the E-Step and the M-Step until convergence. The convergence criterion was set as: |P(b^(m)) - P(b^(m-1))| < 10^-12.

Results and discussion

In this section, we study the effect of false positives and false negatives on inferring regulatory target conservation across conditions/species through both simulations and real data analyses.

We consider an experiment involving the binding of a transcription factor in two different conditions with a total of 1000 genes. We consider the odds ratio as our metric of interest. For fixed true odds ratios, and different values of false positive and false negative rates, we plot the surface of the observed odds ratio in Figure 1. The observed odds ratio is obtained from Equation 4. It can be seen that the observed odds ratio is the largest for low values of false positive and false negative rates and its value decreases with increasing false positive and false negative rates. To visualize this phenomenon in two-dimensions, we fix the false negative rate, and plot the observed odds ratio as the false positive rate varies (Figure 2). We observe that with increasing false positive rates, the observed odds ratio decreases. This is expected, as with an increasing false positive rate, a larger number of true negatives are detected as positives. This reduces the count of genes which are observed negatives in both conditions, i.e. the cell in the contingency table corresponding to "00". Thus, there is a reduction in the observed odds ratio value. Similarly, we observe that for a fixed false positive rate with an increasing false negative rate, the observed odds ratio also decreases. This is also expected, as with increasing false negative rate, a larger number of true positives are detected as negatives. This reduces the number of genes which are observed positives in both conditions, i.e. the cell in the contingency table corresponding to "11", thereby causing a reduction in the observed odds ratio value. To study the effect of asymmetry between p₀₁ and p₁₀, we repeated this simulation for differing values of p₀₁ and p₁₀. We observed similar trends of decreasing observed odds ratios for increasing false positive rate for a fixed false negative rate and a fixed true odds ratio.

Figure 1

Surface of the Observed Odds Ratio for different values of false positive rate and false negative rate and different values of True Odds Ratio. We observe that the surface is highest for low values of the false positive rate and false negative rate, and falls for increasing values of the false positive rate and false negative rate.

Figure 2

Plot of the Observed Odds Ratio versus the false positive rate for different values of false negative rate and fixed True Odds Ratio. The Observed Odds Ratio decreases with increasing false positive rate for a fixed false negative rate.

In the following, we give an analytical proof for the reduction in the observed odds ratio for increasing false positive rates, with the false negative rate being fixed. Equation 4 can be expanded as:

$$\left(\begin{array}{c}E({\widehat{p}}_{00})\\ E({\widehat{p}}_{01})\\ E({\widehat{p}}_{10})\\ E({\widehat{p}}_{11})\end{array}\right)=\left(\begin{array}{cccc}{(1-s)}^2& (1-s)t& t(1-s)& t^2\\ (1-s)s& (1-s)(1-t)& ts& t(1-t)\\ s(1-s)& st& (1-t)(1-s)& (1-t)t\\ s^2& s(1-t)& (1-t)s& {(1-t)}^2\end{array}\right)\left(\begin{array}{c}{p}_{00}\\ p_{01}\\ p_{10}\\ p_{11}\end{array}\right).$$

(11)

The observed odds ratio is:

$$OOR=\frac{E({\widehat{p}}_{00})\ast E({\widehat{p}}_{11})}{E({\widehat{p}}_{01})\ast E({\widehat{p}}_{10})}$$

(12)

where from Equation 11 we get,

$$\begin{array}{l}E({\widehat{p}}_{00})=p_{00}{(1-s)}^2+(p_{01}+p_{10})(1-s)t+p_{11}{t}^2,\hfill \\ E({\widehat{p}}_{01})=p_{00}(1-s)s+p_{01}(1-s)(1-t)+p_{10}st+p_{11}(1-t)t,\hfill \\ E({\widehat{p}}_{10})=p_{00}(1-s)s+p_{01}st+p_{10}(1-s)(1-t)+p_{11}(1-t)t,\hfill \\ E({\widehat{p}}_{11})=p_{00}{s}^2+(p_{01}+p_{10})s(1-t)+p_{11}{(1-t)}^2.\hfill \end{array}$$

We show that for a true odds ratio greater than 1 and for s < 1/2 and s + t < 1, ∂ (OOR)/∂s is negative. These are reasonable assumptions for real data, where the false positives are low and false negatives are not very high. The denominator of ∂ (OOR)/∂s is always positive as it is a squared number. The numerator can be written as:

F = F 1 * F 2 * F 3

where,

F 1 = (p₀₁p₁₀ - p₀₀p₁₁)(-1 + s + t),

F 2 = (1 - s)(p₀₁ + p₁₀ + 2p₀₀s) + (-p₀₁ - p₁₀ + 2p₁₁ + 2(p₀₁ + p₁₀)s)t - 2p₁₁t²,

F 3 = (1 - p₀₀)(-1 + t)t + p₀₀(-1 + t - s(-2 + s + 2t)).

Let us consider each term separately.

If the true odds ratio is greater than 1 and s + t < 1, then p₀₁p₁₀ <p₀₀p₁₁ and (-1 + s + t) < 0. Thus, we have F 1 > 0. F 2 can be simplified as:

F 2 = (1 - s)(p₀₁ + p₁₀ + 2p₀₀s) + (p₀₁ + p₁₀)(-1 + 2s)t + 2p₁₁t(1 - t).

Thus, if (-1 + 2s) < 0, i.e. s < 1/2, then all the three product terms in F 2 are positive. Thus, for s < 1/2, F 2 > 0. F 3 can be simplified as:

F 3 = (-1 + t)(t(1 - p₀₀) + p₀₀) - p₀₀s(-2 + s + 2t) = -p₀₀s² + 2p₀₀s(1 - t) + (-1 + t)(t(1 - p₀₀) + p₀₀).

Thus, F 3 is a quadratic function of s. Since the coefficient of s² is negative, if the discriminant of the quadratic is negative, F 3 is always negative. The discriminant D is given by:

$$D=4p_{00}^2{(1-t)}^2+4p_{00}(-1+t)(t(1-p_{00})+p_{00})$$

Simplifying, we get, $D=-4p_{00}^2(1-t)t-4p_{00}(1-t)t(1-p_{00})$ which is clearly negative. Thus, F is negative for true odds ratio greater than 1, s < 1/2 and s + t < 1.

Simulation Results

To study the effect of false positives and false negatives on statistical inference of dependence between two conditions, we consider a similar setting in which the binding of a regulator in two conditions is studied for 1000 genes. We simulated data for a fixed true odds ratio, and fixed the false positive and false negative rates. We randomly added false positives and false negatives to the data based on the false positive rate and false negative rate. This manifests itself as the observed data, and we repeated this 1000 times. We performed a chi-squared test for independence between the two conditions and counted the number of times the null hypothesis was rejected at a significance level of 0.001. Further, for each observed dataset, we inferred back the true data using our EM algorithm. The inferred true counts are almost equal to the true counts before false positives and false negatives were randomly added. This is because the false positives and false negatives were randomly added based on the fixed false positive rate and false negative rate. These fixed rates are used in our EM algorithm to obtain the inferred true counts. For example, for a true odds ratio of 2, the vector of true counts was (800, 100, 80, 20)^t. We fixed the false positive rate to 0.01 and the false negative rate to 0.2. The vector of inferred true counts was determined to be (799.87, 101.92, 78.55, 20.66)^t. The EM algorithm was initialized by giving equal weights to each possible true binding pattern for each gene. We used the chi-squared test and counted the number of times the null hypothesis was rejected for the inferred true data at the same level of significance. We repeated this analysis for different values of the true odds ratio and different values of false positive and false negative rates. Figure 3 shows the plot of the number of null rejections versus the odds ratio for both the observed and inferred true data. Our results indicate that the number of null rejections for the inferred true data is consistently larger than that for the observed data. We also note that as the odds ratio increases, the difference between the number of null rejections for the inferred true data and observed data also increases. Instead of using the chi-squared test, we also used used thresholds for the odds ratio and positive specific agreement to ascertain the number of null rejections for observed and inferred true data. Thus, for each simulation, we rejected the null hypothesis if the odds ratio was greater than some threshold. We repeated this by applying a threshold to the positive specific agreement. The results are shown in Figures 4 and 5. Thus, in addition to the chi-squared, thresholds for the odds ratio and positive specific agreement also provide evidence that the number of null rejections for the inferred true data is consistently larger than that for the observed data.

Figure 3

Simulation results showing the plot of the number of null rejections versus the odds ratio for both the observed and inferred true data. The number of null rejections for the inferred true data is consistently larger than that for the observed data.

Figure 4

Simulation results showing the plot of the number of null rejections versus the odds ratio for both the observed and inferred true data. The number of null rejections are obtained by applying a threshold of 1.5 to the Odds ratio. The number of null rejections for the inferred true data is consistently larger than that for the observed data.

Figure 5

Simulation results showing the plot of the number of null rejections versus the odds ratio for both the observed and inferred true data. The number of null rejections are obtained by applying a threshold of 0.15 to the Positive specific agreement. The number of null rejections for the inferred true data is consistently larger than that for the observed data.

Real Datasets

We considered two ChIP-chip datasets. Harbison et al. [5] described the binding profiles of 204 transcription factors for S. Cerevisiae in Rich medium, and 84 of these transcription factors were also profiled in at least one other experimental condition. In their study, transcription factors were selected for profiling in a particular environment if they were essential for growth in that environment, or if there was other evidence suggesting their role in gene regulation in that environment. Borneman et al. [8]. studied the divergence of binding sites of regulators Ste12 and Tec1 in the yeasts S. cerevisiae, S. mikatae and S. bayanus under pseudohyphal conditions. They listed genes which showed differing degrees of conservation across the three species, i.e. genes which were targets in only one species, the targets in two species, and the targets all three species.

For ChIP-chip data from Harbison et al. [5], we focussed on the binding data for the transcription factors Ste12 and Tec1 in three different experimental conditions – Rich medium, Filamentation inducing and Mating inducing. We used a p-value threshold of 0.001 to obtain the binding targets for these two regulators. For a pair of experimental conditions, we cross-tabulated the binding targets and created a 2 by 2 contingency table. The odds ratio was quite high, hence we used the log odds ratio and positive specific agreement as metrics to summarize the contingency tables. Thus, given the observed log odds ratio and observed positive specific agreement, for different values of the false positive rate and false negative rate, we inferred the underlying true log odds ratio and the true positive specific agreement using our EM based approach.

In the section describing the model setup, we stated that multiplication of the vector of the observed proportions with the inverse of the transformation matrix could lead to inferred true proportions with negative components. Here we illustrate the scenario. For the regulator Ste12, in Rich Medium and Mating inducing condition the vector of the observed proportions is p = (0.9761, 0.0144, 0.0040, 0.0056)^t. For a false positive rate of 0.001 and a false negative rate of 0.2, multiplying p by the inverse of the transformation matrix results in the vector of the inferred true proportions $\widehat{p}$ = (0.9743, 0.0150, 0.0020, 0.0087)^t. However, for a false positive rate of 0.002 and a false negative rate of 0.3, the vector of the inferred true proportions is $\widehat{p}$ = (0.9748, 0.0144, -0.0005, 0.0113)^t. Similarly, for a false positive rate of 0.004 and a false negative rate of 0.4, the vector of the inferred true proportions is $\widehat{p}$ = (0.9793, 0.0113, -0.0061, 0.0154)^t. Thus, the inferred true proportions obtained by simply multiplying the observed proportions with the inverse of the transformation matrix could contain negative components.

Table 2 shows how the inferred true proportions change for different values of the false positive rate and false negative rate. The vector of the observed proportions is (0.976, 0.014, 0.004, 0.006)^t. Table 3 gives the calculation of the inferred true odds ratios from the inferred true proportions. Figure 6 shows how the surface of the inferred true log odds ratio varies with different values of false positive rate and false negative rate. We notice that as the false positive rate and false negative rate increase, the inferred true log odds ratio differs quite significantly from the observed log odds ratio. We observe similar trends when we use positive specific agreement as the metric of interest (Table 4 and Figure 7). Harbison et al. reported that the false discovery rate in their data was likely to be approximately 4%, while the false negative rate was around 24% for a p-value threshold of 0.001. For binding data from Harbison et al., typically the number of "negatives" was close to 6000, while the number of "positives" was about 100 to 200 at a p-value threshold of 0.001. Thus, the false positive rate was close to 0.001. For our analysis, we studied the variation of observed and true outcomes by varying the false positive rate from 0.001 to 0.005, and the false negative rate from 0.2 to 0.4. Thus, our range of false positive rates would correspond to about 6 to 30 false positives, and about 20 (100 * 0.2) to 80 (200 * 0.4) false negatives, which appears to be quite reasonable. From Table 3, we see that for a false positive rate of 0.001 and false negative rate of 0.20, the inferred true log odds ratio is 5.60, while the observed log odds ratio is 4.56. Since both the log odds ratios are quite high, our inference of association among the two experimental conditions would not be affected by these values of the false positive rate and false negative rate.

Table 2 Inferred true proportions of target genes of Ste12 in the Rich Medium and Mating Inducing conditions.

Table 3 Inferred true log odds ratios for target genes of Ste12 in the Rich Medium and Mating Inducing conditions.

Table 4 Inferred positive specific agreement for target genes of Ste12 in the Rich Medium and Mating Inducing conditions.

Figure 6

Surface of the inferred true log odds ratio for different values of false positive rate and false negative rate for real data. Observed log odds ratio is obtained from cross-tabulation of the binding targets of Ste12 in Rich Medium and Mating Inducing condition.

Figure 7

Surface of the inferred positive specific agreement for different values of false positive rate and false negative rate for real data. Observed positive specific agreement value is obtained from cross-tabulation of the binding targets of Ste12 in Rich Medium and Mating Inducing condition.

We also analyzed the results of ChIP-chip experiments performed by Borneman et al. [8]. We obtained the counts of genes which were the binding targets of the regulators Ste12 and Tec1 in one, two and all three species. We repeated our analysis as described in the previous paragraph for a pair of species (Tables 5, 6 and 7; Figures 8 and 9). Here too, we notice a considerable difference between the observed and inferred outcomes as the false positive rate and false negative rate increases. For example, for a false positive rate of 0.001 and a false negative rate of 0.2, compared to the observed log odds ratio of 4.50, the inferred true log odds ratio is 5.48; however, for a false positive rate of 0.005 and a false negative rate of 0.4, the inferred true log odds ratio is as high as 17.36. Further, we attempted to test the notion that genes falling under similar functional categories tend to be the conserved binding targets across the three species. We listed all the orthologous genes in Yeast. For all these genes, we used SGD GO Slim finder http://db.yeastgenome.org/cgi-bin/GO/goSlimMapper.pl to categorize the genes into broad functional categories. For the top categories which contained the largest number of genes, we cross-tabulated the genes and created a 2 by 2 contingency table based on counts of the genes which are binding targets (Tables 8, 9, 10, 11). For the genes falling in major categories, we notice that the log odds ratios are considerably high, indicating considerable degree of binding conservation. For example, for the 565 genes found to be enriched for Hydrolase activity, 551 were not the binding targets of Ste12 in either S. cerevisiae or S. mikatae. Of the 14 genes which were the binding targets in at least one of the two species, 5 (YNL053W, YDR452W, YGL163C, YIL118W, YMR305C) were the binding targets in both species. Of the remaining 9 genes, 5 (YIR027C, YER133W, YNL180C, YOR049C, YDL047C) were targets in only S. cerevisiae, while 4 (YNL141W, YHR005C, YOR126C, YOL011W) were targets in only S. mikatae.

Table 5 Inferred true proportions of target genes of Ste12 under the pseudohyphal condition in S. cerevisiae and S. mikatae.

Table 6 Inferred true log odds ratios for target genes of Ste12 under the pseudohyphal condition in S. cerevisiae and S. mikatae.

Table 7 Inferred positive specific agreement values for target genes of Ste12 under the pseudohyphal condition in S. cerevisiae and S. mikatae.

Table 8 Cross-tabulation of orthologous genes in functional category Hydrolase activity.

Table 9 Cross-tabulation of orthologous genes in functional category Tranferase activity.

Table 10 Cross-tabulation of orthologous genes in functional category Protein binding.

Table 11 Cross-tabulation of orthologous genes in functional category Transporter activity.

Figure 8

Surface of the inferred true log odds ratio for different values of false positive rate and false negative rate for real data. Observed log odds ratio is obtained from cross-tabulation of the binding targets of Ste12 under the pseudohyphal condition in S. cerevisiae and S. mikatae respectively.

Figure 9

Surface of the inferred positive specific agreement for different values of false positive rate and false negative rate for real data. Observed positive specific agreement value is obtained from cross-tabulation of the binding targets of Ste12 under the pseudohyphal condition in S. cerevisiae and S. mikatae respectively.

Conclusion

In this article, we have studied the effect of false positives and false negatives in the analysis and interpretation of ChIP-chip data. We have derived a relationship between the observed and the underlying true binary outcomes. Given the observed binary outcome of an experiment, we have developed an EM based approach to infer the underlying true binary outcome for given values of false positive and false negative rates.

A common limitation with finding binding targets from ChIP-chip data is that typically an arbitrary threshold, e.g. 0.001, is applied to the data, and all genes with p-values less than this threshold are considered binding targets. The false positive rate and false negative rate for the binding data change with the threshold applied [5]. Datta and Zhao [11] proposed a statistical procedure to determine the binding targets without imposing a simple threshold to the ChIP-chip data. However, their approach relies on accurate inference of false discovery rate [12], which is a non-trivial task.

To summarize data in contingency tables we utilized two commonly used metrics – Odds Ratio and Positive specific agreement. Both these metrics are widely used to study dependency among categorical variables. Since we are interested in quantifying association among two categorial variables, i.e. whether there is association across two different conditions/species, the Odds Ratio and Positive specific agreement are appropriate metrics of interest. In our simulation, we used the chi-squared test of independence to test the null hypothesis that the binding targets are independent. Instead of using the chi-squared test, we could also use the Fisher's exact test to test the independence assumption on the two-way contingency table. The resulting p-values from both tests indicate the statistical evidence against the independence assumption. However, they do not provide a meaningful summary of the degree of dependence as they are also dependent on the sample size.

In general, for independently performed real world experiments, such as two separate ChIP-chip experiments, the independence assumption of equation (9) should hold. This is because we can assume that data points in a particular experiment are independent identically distributed random variables. However, for experiments with closely associated results, it is possible that the false positive and false negative data points for the experiments are not entirely independent. This could result in an under-estimation of the underlying association after the EM procedure.

Due to the limited degrees of freedom of the data, our EM algorithm cannot be used to estimate the false positive rate and false negative rate in experimental data. At each step in the EM algorithm, we estimate three parameters, and we have three equations to solve for them. If we also wish to estimate the false positive and false negative rates, we would have two additional parameters, but the number of equations would still be three. This would lead to an identifiability problem.

We initialized the EM algorithm by different initial estimates of the parameters. For each initial estimate of the parameters, the algorithm converged. The convergence criteria for the EM algorithm require that the log likelihood of the parameters l(b|g) be continuous and differentiable in the parameter space. Unfortunately, the M-step of our algorithm does not have a closed form. Hence, it is difficult to evaluate the gradient of the log likelihood function.

Harbison et al. performed their ChIP-chip experiments using microarrays consisting of spotted polymerase chain reaction (PCR) products representing all the intergenic regions of Saccharomyces cerevisiae. To obtain the binding targets a p-value threshold was applied to the binding intensities associated with the probes. One of the drawbacks of PCR based arrays is the low resolution of the DNA elements in the microarray chip. For PCR arrays designed for Yeast, the typical resolution achieved is less than 1 kb. In recent years, high density oligonucleotide arrays, comprising of large numbers (40, 000 to more than 6, 000, 000) of short oligonucleotides have been utilized for ChIP-chip studies [13–16]. A number of statistical algorithms have also been developed to determine the binding targets from such large scale tiling arrays [17–20]. Borneman et al. used high density oligonucleotide arrays to perform their experiment. The binding targets were obtained using Tilescope [21]. Since they report the target genes in each organism, we simply used their results to obtain the counts of target genes in each of the three organisms.

Our analysis can be applied to any experimental setting with binary outcomes. However, for the sake of simplicity, we have illustrated its application for ChIP-chip experiments. By applying our algorithm to ChIP-chip data from Harbison et al. and Borneman et al., we observe that for different values of the false positive and false negative rate, the observed and true metrics for the binary data can differ quite dramatically. However, we notice that when the true log odds ratio is greater than 4, i.e. there is a significant degree of association among the binding targets across conditions/species, such differences in the observed and true metrics would not change our inference. On the other hand, our simulation results indicate that when the true odds ratio is close to 1, i.e. for cases when the underlying association is marginal, moderate values of false positive and false negative rates (0.01 and 0.2 respectively) may not be able to provide conclusive evidence of any underlying association or independence.