Yeast two-hybrid screens are an important method for mapping pairwise physical interactions between proteins. The fraction of interactions detected in independent screens can be very small, and an outstanding challenge is to determine the reason for the low overlap. Low overlap can arise from either a high false-discovery rate (interaction sets have low overlap because each set is contaminated by a large number of stochastic false-positive interactions) or a high false-negative rate (interaction sets have low overlap because each misses many true interactions). We extend capture–recapture theory to provide the first unified model for false-positive and false-negative rates for two-hybrid screens. Analysis of yeast, worm, and fly data indicates that 25% to 45% of the reported interactions are likely false positives. Membrane proteins have higher false-discovery rates on average, and signal transduction proteins have lower rates. The overall false-negative rate ranges from 75% for worm to 90% for fly, which arises from a roughly 50% false-negative rate due to statistical undersampling and a 55% to 85% false-negative rate due to proteins that appear to be systematically lost from the assays. Finally, statistical model selection conclusively rejects the Erdös-Rényi network model in favor of the power law model for yeast and the truncated power law for worm and fly degree distributions. Much as genome sequencing coverage estimates were essential for planning the human genome sequencing project, the coverage estimates developed here will be valuable for guiding future proteomic screens. All software and datasets are available in Datasets S1 and S2, Figures S1–S5, and Tables S1−S6, and are also available from our Web site, http: //www. baderzone. org.
Maps of pairwise protein–protein interactions are being generated in increasing numbers by the two-hybrid method [1]. Genome-scale two-hybrid screens have now been conducted for Saccharomyces cerevisiae (yeast) [2,3], Caenorhabditis elegans (worm) [4], and Drosophila melanogaster (fly) [5]. More recently, screens have been reported for herpesviruses and human [6–8]. These datasets have stimulated large-scale analysis of the topology of protein interaction networks. Limitations in the data, both false positives (spurious interactions reported from high-throughput screens) and false negatives (true interactions missing from the screens), continue to make it difficult to infer network properties [9–11], including distinctions as basic as the difference between Erdös-Rényi (ER), power law [12–14], and other network degree distributions [15]. A recent review points out the challenges in estimating false-positive rates, false-negative rates, and completion to full coverage of protein interaction networks [16]. Virtually every published method falls back to an estimate based on intersections of datasets. For false-positive rates, these methods have large variance when assays have little overlap, and indeed could not be used to analyze the existing large-scale maps for worm and fly. Estimates for false-negative rates based on overlap of datasets may have even larger uncertainty. Finally, global estimates of false-positive and false-negative rates say little about protein-specific properties, including whether certain classes of proteins behave well or badly in two-hybrid screens. The goal of this work is to develop and apply a statistical model for two-hybrid pairwise interaction screens. Previous methods typically summarize the presence or absence of an interaction as a 1/0 binary variable, and possibly split off a high-confidence core dataset. The method we describe reaches back to the raw counts of observed bait–prey clones. This frees the statistical method from the need for an external gold standard of true-positive and true-negative interactions, or even a second dataset. It permits protein-specific predictions that for the first time permit tests of hypotheses that some classes of proteins are more or less likely to have nonspecific interactions. Finally, estimates of false-negative rates permit statistically grounded confidence intervals for the total number of pairwise interactions present in model organism proteomes. A flowchart of a two-hybrid screen orients the discussion by showing where true-positive interaction partners can be lost and where false-positive, spurious interactions may arise (Figure 1). In a two-hybrid assay, one protein is fused to the binding domain (bait construct) of a yeast transcription factor, and a second protein is fused to the activation domain (prey construct). Physical interactions between bait and prey proteins reconstitute transcription factor activity. Due to the expense of the assay, not every protein may be selected to be made into a bait or prey construct. Furthermore, some constructs may not be functional at all due to improper folding or incompatibility with the two-hybrid system. These missing interactions are important to consider when estimating the total number of interactions in a proteome. High-throughput two-hybrid screens have used multiplexed pairwise tests, either by testing a single bait versus a pool of preys [4,5], or by pooling both baits and preys [3]. Unnormalized prey pools can be generated from mRNA extracted from growing cells. With access to clone collections, pools can be normalized by designing baits and preys individually for each protein or protein domain, then mixing preys in equal proportion. The yeast screen considered here [3] tested 62 normalized bait pools versus 62 normalized prey pools, each pool having approximately 96 genes. The fly screen and worm screen each tested one bait in turn versus both normalized and unnormalized pools. The testing occurs by using mating or transformation to express both the bait and prey construct in a single yeast cell. True-positive interactions drive reporter genes that permit the yeast cell to grow in selective media. Yeast cells whose bait–prey constructs do not interact are expected to drop out during the population expansion. True positives may also be lost during the population expansion for at least two reasons. First, the mating or transformation may lack enough cells to ensure that every combination is tested. Second, a particular construct may have domain-specific misfolding, making it functional for some interactions but nonfunctional for others. True interactions that are not represented in the cells following the population expansion are systematic false negatives for a particular screen. False negatives due to insufficient mating/transformation and due to nonfunctional domains could in principle be discriminated by repeating the mating or transformation step and the selective population expansion. Without this additional step, however, losses during the population expansion combine to yield a systematic false-negative rate termed 1 − psyst, with psyst representing the true-positive rate for an interacting pair to survive the population expansion. Some cells expressing noninteracting proteins may also survive the population expansion, and the final population of cells will be a mixture of true positives and false positives. In Figure 1, the mass fraction of true-positive cells is 1 − α, and of false-positive cells is α. The ratio of false positives to the total number of true negatives is the false-positive rate. Usually, however, the ratio is with respect to the total number of observed interactions (Equation 31), defined as the false-discovery rate and synonymous with the parameter α. An ongoing point of contention in two-hybrid screens is the possibility that two proteins that never interact in vivo in the host organism might have a strong, reproducible interaction in vitro in the engineered two-hybrid system. Conversely, proteins with a strong two-hybrid interaction might nevertheless fail to interact in vivo. For the purposes of this work, we assume that such cases are rare and we classify any pair of proteins with a reproducible two-hybrid interaction as a true positive. While the total false-positive fraction α may be large, it represents a sum over many different false-positive pairs. Most models, including ours, assume that any particular false positive is rare, with vanishing probability of observing a specific false-positive interaction more than once. Interactions detected in pooled screens often require sequencing to identify the interacting partners, although advanced pooling designs may improve deconvolution efficiency [17]. Cost constraints limit the number of interactions that can be sampled for sequencing. If the number of clones selected for sequencing is smaller than the number of true interaction partners of a bait, some true partners will certainly be lost. Limited sampling depth also truncates the observed degree distribution for baits. The false-negative rate due to undersampling is termed 1 − psamp in Figure 1. False-discovery rates have typically been estimated by comparing datasets [18–20], suggesting up to 50% false positives, but these analyses can confound false-positive and false-negative error sources. Estimated error rates have large uncertainty because few interactions are observed in multiple datasets. For example, comparing the Uetz and Ito two-hybrid datasets for yeast reveals only 9. 1% of the total interactions in common [3], and comparing the two-hybrid interactions with mass spectrometry interactions reveals only 0. 6% in common [20]. Similarly, comparison of two fly screens reveals few interactions in common [5,21]. Cross-species comparisons have also revealed little overlap in the reported interactions [4], although protein and network evolution are additional confounding factors. Efforts to estimate the true number of interaction partners of a protein have used contingency tables for observing an interaction in multiple screens. These methods require that all the interactions be true positives, for example by excluding singleton observations [22], which can reduce the estimated interaction count. A notable exception is previous work in the context of mass spectrometry of protein complexes [23], which used a Bayesian model to infer global parameters for screen-specific false-positive and false-negative rates. These parameters then provided posterior estimates for the probability of a true interaction given results of one or more screens. This work is important in using the number of trials and successes, rather than a single summary yes/no observation, in its probability model; it serves as motivation for developing similar models for the more complicated two-hybrid sampling process involving strong protein-specific effects. Quantitative predictions of the amount of work required to identify some fraction of true interactions would be analogous to formulas for genome sequencing [24] and would be useful for planning new experiments [25]. The new work presented here uses the raw screening data to estimate the false-negative rate from undersampling, together with the false-positive rate. A schematic illustrates the sampling process (Figure 2). Interactions are sampled with replacement from two sets, one representing true positives and the other true negatives. The observations are the number of times that each interaction is sampled, which we summarize with three variables: n, the total number of samples drawn; w, the number of unique interactions within the n samples; and s, the number of interactions observed exactly once. From these observations we are to estimate the unknown (hidden) values of k, the total number of true interaction partners, and f, the number of false positives within the sample n. We also estimate the parameter α representing the fraction of false positives in the mixture (the false-discovery rate), as well as parameters representing the probability distribution for k. For simplicity, the illustration suggests sampling interactions in the entire network; in reality, this sampling process occurs separately for each bait, and the estimation of k and f is performed separately for each bait. This estimation problem is akin to estimating population sizes or species counts from capture–recapture experiments, estimating vocabulary size from word counts, estimating the number of distinct alleles at a particular locus, and estimating the number of facts in the scientific literature [26–33]. Classic capture–recapture theory permits heterogeneous capturability rates, here analogous to different probabilities of observing each true interaction partner of a bait. The canonical estimator has a simple form: k̂ = w + s2/2k2 [34–36], where k2 is the number of partners observed exactly twice. The classic estimator fails in the two-hybrid setting because it does not account for false positives. To our knowledge, false positives have never been discussed in the capture–recapture setting. False positives will vastly inflate the interaction count by adding to the number of singleton observations, s, and to the total observed count, w. The standard estimator has high variance when the number of observations is small, yielding a small value for the denominator k2. The estimator fails to converge when each partner is observed only once, yielding n = w = s, k2 = 0, and k̂ → ∞. We present a front-to-back statistical model for both false-positive and false-negative error rates in two-hybrid screens. A glossary of model terms is provided (Table 1). The overall approach is to start by estimating the parameters of a mixture model for true positives and false positives following the population expansion. This permits us to estimate bait-specific false-discovery rates and false-negative rates due to undersampling. We can then back-calculate the false-negative rate due to systematic effects. Putting the results together yields an overall estimate for the false-negative rate of a screen and a basis for comparing interaction lists produced by different efforts. Along the way we examine issues that our model is able to address quantitatively: selecting the best model for the protein degree distribution; correlating false-discovery rates with bait properties such as “sticky” or “promiscuous” domains or hydrophobic regions; and determining the relative performance of prey libraries generated from cDNA libraries or ORFeome collections.
We applied our methods to experimental data from two-hybrid screens conducted in the model organisms yeast, worm, and fly. The key parts of the datasets are the numbers of times that a specific bait identifies each prey, from which all other required values may be calculated. Yeast data was taken from ITO FULL, with clone counts from the IST HIT column [3]. Worm data was from WI5 with clone counts in the NumHitADcDNA and NumHitADORF columns [4]. The worm interactions were from the CORE_1, CORE_2, and NON_CORE sets; interactions annotated as SCAFFOLD (previous screens by the same group), LITERATURE (interactions reported in the scientific literature), and INTEROLOG (interactions inferred cross-species) were excluded. Fly data was from the CuraGen screen with clone counts in the baitprey and preybait columns [5]. A summary of the data sources is provided (Table 2), and a compendium of the data sources is available (Dataset S1). In collecting these datasets, we noted that many two-hybrid screening publications do not report the clone counts that are required for capture–recapture analysis. This includes one of the two major yeast high-throughput screens [2], a screen for Helicobacter pylori interactions [37], and important recent screens for human protein–protein interactions [6,7]. Part of the motivation of this work is to demonstrate the value of making this type of raw data available for analysis. The relevant variables describing a two-hybrid screen are listed in Table 1 and summarized here. Each of N baits is screened against a prey library. For bait i, ni clones from a two-hybrid screen are sampled and the preys are identified. The number of times that prey j occurs within bait i' s sample is termed nij. The number of unique preys within the ni clones is termed wi. The number of preys observed exactly once (singletons) is si. The ni clones comprise a mixture of false positives and true positives, but it is not known a priori which are the false positives, or even the total number fi of false positives. The goal of our analysis is to estimate the number of false positives, fi, and the number of true positives that were left unsampled for each bait. Our statistical model makes the following assumptions. Prey constructs are either functional (with probability psyst) or systematically lost (with probability 1 − psyst) with respect to a particular true interaction partner bait construct. Due to possible differences in binding sites, a prey may be functional for one bait and nonfunctional for a different bait. The total number of true positives for a particular bait i is termed κi, of which ki ≡ psystκi are functional. The parameter psyst is estimated from the observed probabilities of bidirectional interactions. Prey libraries are normalized, with each prey present at equal concentration. True-positive interaction partners are sampled with equal probability with replacement from the ki functional preys. False-positive preys occur stochastically, not systematically, with a low probability per prey and negligible probability that any single true negative is sampled twice for a given prey. Thus, clones observed once are a mixture of false positives and true positives; clones observed two or more times are assumed to be true positives. The cumulative probability that a particular clone is a false positive may be large because it sums over all the possible true negatives. This cumulative false-positive rate is the false-discovery rate per clone, termed αi for bait i, and may be different from bait to bait. These assumptions are justified in the Materials and Methods section. Even if restrictive, they still provide a necessary starting point for building more complicated models. Given these assumptions, we show in Materials and Methods how false-discovery rates and corrected counts of interaction partners can be determined for each bait. The posterior estimates for false-discovery rates and interaction counts depend on the functional forms selected for the bait-to-bait heterogeneity in the false-positive rate and the protein interaction degree distribution. We used a variety of model selection criteria, also described in Materials and Methods, that had perfect performance on simulated data. While false positives are a recognized byproduct of two-hybrid screens, there has been little work to investigate bait-to-bait variation in the false-discovery rate. We investigated three models for bait-specific false-discovery rates, described in words here and mathematically in Materials and Methods, Equation 2. The false-positive rate in the model is expressed per sampled clone, rather than per prey in the library (which would be a much smaller error rate) or per unique interaction (which would be a larger error rate). We selected three representative functional forms as possible models for the probability that a bait protein has k functional interaction partners in the prey library, described in text here and more formally in the section Theory. The summary results, Table 3, extrapolate the number of interaction partners from the estimated number of true positives within the preys screened to the total number in the proteome. The results suggest about 40 pairwise interaction partners per protein in yeast, and roughly 100 pairwise interaction partners per protein in worm and fly. These numbers, however, are based on the estimated means. For long-tailed degree distributions, the median values may provide greater intuition, and may in fact be more robust by discounting outliers with high interaction counts. Median numbers of interaction partners obtained from parametric degree distributions (see the section Total Interaction Counts), are provided at the bottom of Table 3. The final values obtained are roughly ten partners per yeast protein, 61 for worm, and 46 for fly. The 1. 5-fold difference between worm and fly might point to built-in biases in the screens (different baits and preys, different selection thresholds, etc.) rather than any fundamental biological differences. Using the median and mean estimates as brackets, our results suggest between 30,000 and 140,000 pairwise interactions in yeast; 600,000 to 1,200,000 in worm; and 300,000 to 600,000 in fly. Other work, using a contingency table approach similar to k∩, has suggested a 95% confidence interval of about 40,000 to 75,000 interactions in yeast, and 150,000 to 370,000 in human [16]. This previous work was unable to make predictions for worm or fly, however, due to the lack of multiple datasets for comparison.
The methods introduced here provide a new model for false-positive and false-negative rates for two-hybrid screens. To our knowledge, this is the first model that considers the number of observations of each prey, as opposed to a binary interaction / no-interaction summary statistic, to calculate these rates. We have validated the model thoroughly using simulated data and using published biological datasets. The applications to published data demonstrate the crucial ability to predict how many new interactions will be observed as more preys are collected, together with the true-positive and false-positive fractions. One of the major criticisms of the two-hybrid method has been a high false-positive rate. Unlike previous methods that produce average false-discovery rates over an entire screen, our method provides bait-by-bait estimates. False-discovery rates are heterogeneous: some baits perform better than others. As others have suggested [38], this permits the possibility of correlating false-positive rates with hydrophobicity and related protein properties. We find strong evidence for higher false-discovery rates for membrane proteins, but not for hydrophobic proteins in general. Two-hybrid screens such as the split-ubiquitin system [55] have been developed to detect interactions between membrane proteins. These assays could very well show a correlation of false-positive rates with other classes of proteins. Classification of proteins according to enzymatic function reveals that those in signaling pathways have lower false-discovery rates than those in metabolic pathways. This suggests greater evolutionary pressure to maintain specificity of information-carrying networks. One suggested mechanism of network evolution is that recent paralogs may continue to share interaction partners. This would imply that proteins within a single family should show cross-reactivity with each other' s binding partners, eventually leading to false positives due to weak remnants of ancestral interactions. We rejected this hypothesis by finding no significant correlation between false-discovery rate and family size. Analysis of false-positive rates also provides a quantitative estimate of the value of using constructs from a sequence-verified ORF collection rather than from cDNA libraries. When we classify bait constructs as “good” or “bad, ” we find that the “good” category is 90% for an ORF collection and 45% for a cDNA library. On the prey side, using an ORF library produces one-third fewer false positives than a cDNA library. This model yields estimates of false-negative rates from screening statistics, and to our knowledge is the first attempt to discriminate between false negatives due to undersampling and false negatives due to biological and systematic effects. We find that sampling and systematic factors are both important contributors to false negatives, with undersampling yielding a roughly 2× reduction in interactions, and systematic effects yielding an additional 2× to 6× reduction. False-negative rates estimated from the statistical model are in general agreement with those estimated from comparisons between datasets or to a gold standard. The statistical framework provides a convenient route to assessing the likelihood of different population-level functional forms for the protein degree distribution and the false-discovery rate. We provide conclusive evidence that, among one-parameter degree distribution models, a PL model is far superior to an ER model. We find evidence for exponential truncation of the degree distribution in worm and fly, but not in yeast. The number of interactions per protein is predicted to increase from about ten partners for yeast to about 50 partners for the metazoans worm and fly. These results suggest that more complex organisms have more interactions per protein component, as well as more components overall. This model will have value in application to ongoing pool-based assays for protein–protein interactions in model organisms and human. An immediate demonstration is the ability to predict the total number of pairwise protein–protein interactions based on two-hybrid data. We suggest that the total number of pairwise interactions observable by the two-hybrid system is roughly 140,000 in yeast, and 600,000 to 1,300,000 in worm and fly, with about 95% remaining to be discovered. An attractive extension of the model presented here is to include unequal capture probabilities for true interaction partners. The current model represents true-positive preys as a two-component mixture: a fraction 1 − psyst of true-positive preys are considered absent from the pool, with capture probability 0; the remaining k true-positive preys have uniform capture probability. It would be possible to include more components, or even a continuous variable representing an inhomogeneous capture probability of a prey. This is important for libraries generated directly from mRNAs with varying abundances, and could still be important for libraries generated from normalized clone collections due to varying effective nuclear concentrations and binding constants. Including heterogeneous capture rates for baits could be accomplished by extending the model to represent the true-positive rate psyst for each bait as an additional hidden variable to be optimized within the Expectation–Maximization (EM) framework. In this work, psyst is a global parameter calculated after the sampling-based parameters have been estimated. Both of the above extensions would involve a probability model that considers interactions in both directions, bait–prey and prey–bait, and would necessarily add complication to what is already a mathematically detailed model. While a more complicated model would seem unlikely to lead to different conclusions from those presented here, it could answer questions relating capture probabilities to protein physical properties, protein abundances in libraries, and transient versus stable protein interactions. Developing related statistical models for other types of protein interaction screens will also be important. A constant proviso attached to interaction screens is the suspicion that methods such as two-hybrid screens, affinity pull-downs [56,57], and protein binding chips [58] will identify different subsets of interactions. Quantitative comparisons are difficult, however, because systematic assay-specific differences are confounded with random loss of interactions due to incomplete sampling. Methods such as the one presented here will contribute to understanding what different screening technologies tell us about the proteome.
An overview of notation is provided (Table 1). Consider a particular protein j used as one of N baits in a two-hybrid screen against a pool of Γ species of preys of which κj are true interaction partners. We assume that κj = Γ so that Γ − κj ≈ Γ is a good approximation for the number of true negatives for each bait. We model the first stage of a two-hybrid screen as an all-or-none process reflecting whether a bait mates successfully with a prey and yields progeny that survive selection. For simplicity, and to reduce the number of free parameters, we assume an identical systematic true-positive rate psyst for each bait j with each of its κj true interaction partners. The parameter psyst includes systematic biological effects, such as generating functional fusion proteins in the two-hybrid system. The number of surviving true positives is kj with binomial distribution We further assume that the Γ true negatives continue to grow slowly in the selective media with a population expansion that is only p* times the population expansion of surviving true positives. We assume a stochastic, not systematic, model for false positives, with p* = 1 and identical for each prey. Also for simplicity, and in accord with prey libraries constructed from normalized ORF collections, we assume that each prey is present initially at equal concentrations. The final mass fraction of false positives is denoted αj = p*Γ[kj + p*Γ], yielding a scaled error model that depends on a single constant p*Γ ≡ a. More generally, a variety of error models are possible: The MIXTURE model introduces an index zj ∈ {1,2…, m} to one of m possible values of α and prior probabilities π (1) + π (2) + …+ π (m) = 1 for the m components. With m = 2, this permits “good” baits (zj = 1) and “bad” baits (zj = 2) with α (1) ≤ α (2). The second stage of screening bait j is to sequence nj clones from the mixture of true positives and false positives. We assume that each of the kj true positives is sampled with uniform probability (1 − αj) /kj and each of the Γ false positives is sampled stochastically with uniform probability αj/Γ. The number of times that prey species m, either a true positive or a false positive, is sampled within the nj clones is njm, with 0 ≤ njm ≤ nj and ∑mnjm = nj. The probability of the observed counts njm is a multinomial, [nj! /Πmnjm! ]Πm, with θm = (1 − αj) /kj or αj/Γ. As is typical in a capture–recapture setting, it is more convenient to work in the context of abundance classes. Let for i ≥ 1 represent the observed data as the number of preys observed exactly i times within the nj samples. For convenience, we introduce sj as a synonym for, the number of singleton preys observed only once. The total number of distinct preys observed is wj, where δ (arg) is 1 if its argument is true and 0 if false. The standard generalized multinomial distribution is obtained by summing over the {S} permutations that yield wj distinct species, identical to Equation 3 of [34]. A rough motivation for this formula is that nj/Πi≥1 is the number of distinct permutations of the nj clones, |S|/Πi≥1 is the number of distinct permutations of the observed species, and Πm is the probability of selecting the species in specified order. Our final approximation is that each true negative occurs at most once as a false positive, njm = 0 or 1 when m is a true negative. The expected number of false-positive clones within the nj clones is αnj. The probability that each of these, selected at random from the Γ total possibilities, is distinct is [1 − (m − 1) /Γ], or approximately exp[− (m − 1) /Γ] = exp[−αjnj (αjnj − 1) 2Γ], analogous to the Birthday Paradox (the probability that two people in a large random group share a birthday). An appropriate constraint ensuring distinct false positives is that nj ≤ /αj. With genome-size prey libraries, Γ ≥ 5,000, and we anticipate that αj ≤ 0. 5, making this approximation valid for nj ≤ 140. For yeast, ten baits (0. 67%) violate this constraint; for worm, 18 baits (2. 5%); for fly, ten baits (0. 27%). Denote fj as the number of false-positive observations within the sample nj. By the above assumption, the false positives must be within the sj singletons, and 0 ≤ fj ≤ sj. Using the uniform capture probabilities, Πm =. The number of permutations |S| can be calculated under the above assumption of singleton false positives as The first factor is the number of ways that false positives can be assigned to a subset of fj of the sj singleton species. We have used Γ! / (Γ − fj)! ≈, which is valid because Γ ≫ nj ≥ fj. The second factor is the number of permutations that select the wj − fj observed true positives out of kj. Combining results yields with an additional factor of π (zj) depending on the hidden variable zj that indicates the component for the MIXTURE error model, Equation 2. The probability distribution for the hidden variables fj and kj are obtained through the Bayesian relation For the MIXTURE model, the analogous equation includes the hidden variable zj. When kj is independent of nj and αj, Pr (kj | nj, αj) = Pr (kj) ≡ Pr (kj | Φ), where Φ comprises one or more global parameters describing the interaction degree distribution. The simplified Bayesian result is or Pr (yj | xj, Q) where the hidden variables yj = { kj, fj}, and possibly zj; the observed variables xj are the counts of singletons (sj), distinct preys (wj), and total samples (nj); and the parameters Q are the global parameters for the error model (a, α, or {α (1) … α (m); π (1) …π (m) }) and the protein degree distribution. The three summary statistics {sj, wj, nj} are sufficient statistics for the observed data due to the assumption of homogeneous probabilities for observing each true-positive and true-negative species. The sum over k formally starts at wj − f, which may equal 0 when each of the nj observations is a singleton. In the results, however, we restrict attention to probability distributions for which Pr (k = 0 | Φ) = 0 and start the summation at k = 1. Three distributions are considered: The normalization of the Poisson distribution by 1 − e−λ in Equation 11 reflects that the summation begins at k = 1 rather than at 0. In keeping with the definition of 1 − psyst as the systematic false-negative rate, it may be more appropriate to use parametric distributions for Pr (κ), than to obtain Pr (k) as the convolution Pr (k) = ∑κ≥k Pr (k | κ) Pr (κ). We are in a sense replacing Pr (k | κ) of Equation 1 by a delta function near the mean value κpsyst in order to retain the form of a simpler parametric distribution. Estimates for {kj, f} for each bait could in principle be obtained using Equations 10–11. This requires, however, estimates for the global parameters Q. Furthermore, the asymptotic form of the summand in Equation 10 is Pr (k | Φ). Writing the asymptotic form of Pr (k | Φ) as k−ɛ, existence of a maximum a posteriori estimator requires nj + ɛ − wj > 0; convergence of the sum requires nj + ɛ − wj > 1; and convergence of the mean of kj requires nj + ɛ − wj > 2. For the ER prior, λ > 0 guarantees convergence of all powers of kj; for the TPL prior, c > 0 guarantees convergence. Convergence could be achieved by normalization of k in Equation 11 to a finite cutoff Γ rather than to ∞. In practice, however, results for the PL model are sensitive to the cutoff value when ɛ < 2. The TPL and ER models are not sensitive to a cutoff, as both provide a natural cutoff as part of the model parameters. If a cutoff is appropriate, we anticipate that these models will provide improved descriptions of a degree distribution. To overcome both these difficulties, we use EM to obtain parameter estimates Q̂ that maximize the probability of the observed data [59,60], assuming a uniform prior Pr (Q). We introduce the notation for the mean of a generic function F (x, y) of the hidden and observed variables. The sum over the hidden variables expands to for the scaled α and single α error models, and to for the m-component mixture. As mentioned before, while some models permit kj = 0, the power law model does not; for consistency, we start all degree distributions at kj = 1 and the lower limit of the sum over kj is effectively max (1, wj − fj). The standard equations giving a new parameter estimate Q in terms of a previous estimate Q′ are where the simplification holds because Pr (xj, yj | Q) = Pr (xj | yj, Q) Pr (yj | Q) and Pr (xj | yj, Q) = Pr (xj | yj) is independent of Q. Update equations for the error models are as follows: Update equations for the degree distribution are as follows: An interesting and unfortunately common boundary case occurs when only a single clone is sequenced for a bait, nj = 1. In these cases, sj and wj must also be 1, and Pr (xj | Q) = 1 regardless of Q. Thus, baits with n = 1 do not affect the final model parameters because the partial derivatives of their contributions to the log-likelihood are always 0. The appearance of the expectation of log k rather than k in the EM equations for the power law parameter ɛ in the PL and TPL models suggests the use of the posterior mean of log k as a route to estimating the hidden variable decay. We define this estimator as k̂, The three error models and the three degree distribution models yield a total space of nine possible models, with varying degrees of freedom (df): 1 df for the scaled and single error models; 2m − 1 df for the m-component mixture error model; 1 df for the ER and PL degree distributions; 2 df for the TPL distribution. We used three separate criteria to assess which model provides the best fit: log-likelihood cross-validation (CV); full data Bayesian information criterion (BIC); and bootstrap BIC. The CV method with F-fold cross-validation divides the full data into F subsets. For subset f, model parameters Qf are estimated using the remaining F − 1 subsets, and the log-likelihood of subset f is calculated as loglikf = log Pr ({xf}|Qf). This procedure is repeated F times, once for each subset. Thus, each subset is used F − 1 times to obtain model parameters and 1 time to obtain an unbiased log-likelihood. The final log-likelihoods, ∑f loglikf, can be compared directly. The statistical significance of a difference in log-likelihoods for two models can be assessed using a paired test, such as the nonparametric Wilcoxon rank signed test, for the differences − for pairs of models M and M′. The BIC is an appropriate heuristic for performing model selection in the context of maximum likelihood parameter estimation for Q and a uniform prior over model classes M: where d is the number of df in the model and N is the number of sets of observations, here baits. A smaller BIC indicates a more likely class of models, and the term d log N penalizes more complex models. Overfitting is unlikely for our models: the typical number of sets of observations N ∼ 1,000, while the models have only two to four free parameters. Sometimes, the BIC heuristic may indicate a small preference for one model over another. Bootstrap replicates may be used to assess the stability of the BIC results. Bootstrap replicates are constructed by selecting N examples from the full data of N examples uniformly and with replacement. Thus the number of times n that an example occurs in a bootstrap replicate is approximately Poisson with Pr (n) = 1/ (n! e). The BIC heuristic for each model is then calculated for each bootstrap replicate, and the number of times that each model has the best BIC score is recorded. We calculated the cumulative number of clones sampled for a domain, ndom, and the cumulative posterior estimate for the number of false positives, f̂dom, by summing over the counts for each protein annotated as having that domain: p-Values for the upper and lower tail, p> and p<, were calculated assuming a binomial distribution with ndom trials and a success rate equal to the overall false-discovery rate α̂ for each organism (0. 093 for yeast, 0. 122 for worm, 0. 157 for fly). To ensure a conservative test, fractional values of f̂dom were rounded down for the upper-tail test and rounded up for the lower-tail test, Finally, the single-value p-values were adjusted for the number of domains observed among baits in each species (783 for yeast, 473 for worm, 1,310 for fly). When two domains refer to an identical subset of proteins, results for only a single domain are displayed. Parameter estimates. Simulations were performed separately for each of the three protein degree distributions (ER, PL, TPL) combined with each of the three error models (scaled, single, mixture), yielding nine total models. Simulations over a range of parameter values used N = 1,000 baits and n = 10 clones sampled per bait, and were usually repeated three times at each parameter value (Figures S1–S3). The agreement between known and estimated parameters is very good over parameter values that span the estimated values for the published datasets (Table S1). Agreement is quantified by the root-mean-square (RMS) difference between the known and estimated parameters and the R2 goodness of fit measure, defined as 1 − ∑t (θt − θ̂t) 2/∑t (θt − $theta; ̄) 2 with sum t over trials, θt the true parameter value for trial t, $theta; ̄ the mean of θ over the trials, and $theta; ̄t the estimated value for trial t. The RMS values are generally less than 0. 05 in absolute units. The R2 values for the TPL-mixture range from 0. 65 for the false-discovery rate parameter to 0. 98 for the power law parameter. The R2 values depicted for the TPL-mixture model are not typical but rather the worst results obtained over all combinations of degree distribution and error model (Figure S1–S3). Other models with fewer parameters are more accurately estimated. For example, simulations using the PL-mixture model yield an R2 of 0. 99 for ɛ, 0. 97 for π (1), and 0. 88 for α (1) and α (2). Simulations using the PL-single model and PL-scaled model yield R2 ≥ 0. 99 for all parameters (Table S1). Hidden variable estimates. Simulations also assessed the ability to predict hidden variables (Figure S4). We simulated a dataset using parameters obtained by fitting the fly experimental data with a TPL degree distribution and a two-component error model. The simulated dataset had N = 1,000 baits and n = 10 baits per prey. Next, the EM algorithm was used to estimate the model parameters ɛ̂, α̂ (1), α̂ (2), and π̂ (1). The estimated parameter values, rather than the true known values that generated the data, were used to avoid introducing a favorable bias in the hidden variable predictions. The converged model parameters were then used to obtain the posterior estimates 〈log k〉, 〈f〉, and 〈δ (z = 1) 〉 for each bait, with the notation 〈…〉 defined by Equation 13 (Figure S4). Using the log-transform of k is more natural than using k due to the long tail of the power law distribution. It is also motivated by the EM equations for the power law exponent ɛ, which depends on 〈log k〉 rather than 〈k〉 as shown in Equation 18. The hidden variable 〈log k〉 is estimated with R2 = 0. 84, indicating good correlation between true and estimated values. The RMS error in the estimate is 0. 42, indicating the ability to predict the true value of k within a factor of exp (±0. 42), or 1. 5-fold. The RMS for the number of false positives is 1. 1, which means that the estimate for the number of false positives for a bait is usually within 1 of the true count. Estimates of the error component of a bait are accurate for low values (low error rate) and high values (high error rate) of ẑ. Baits with intermediate estimated values, 1. 2 ≤ ẑ ≤ 1. 6, may come from either error component. For these baits, all but one or two of the preys are singletons, and it is difficult to determine whether this is due to a large bait degree or a large error rate. Model selection for simulated data. We next validated the BIC heuristic for model selection (Table S2). In this test, we used each of the nine possible models (three degree distributions × three error models) to generate datasets with N = 1,000 baits and n = 10 preys per bait, then calculated the log-likelihood for each of the nine models. A total of 81 fits were performed (nine generative models × nine fitting models). The parameters of the generative models were deliberately selected to yield similar data by using the values obtained by fitting the experimental worm data (Table 2). In each case, the BIC identified the model accurately (Table S2). The probability to obtain a perfect result by chance is approximately (1/9) 9, or 2. 6 × 10−9. The BIC results indicate the TPL model can provide a good fit for data generated by all three models: ER, PL, and TPL. The TPL model includes the PL model as a special case with the exponential decay constant c = 0. A large value of c, which truncates the degree distribution, permits the TPL model to mimic the ER model. The BIC score adds a penalty of log 1,000 = 6. 91 to the TPL fit to account for the extra parameter. In several of the entries of Table S2, this penalty is essential to select the true generative model over the TPL model. Properties of the experimental datasets for yeast, worm, and fly are summarized at the top of Table 2: total number of baits, N; number of preys sampled per bait, n; number of unique observed interaction partners per bait, w; and number of interaction partners observed a single time, s. Each dataset was fit using three separate degree distribution models (ER, PL, TPL) and three error models (SCALED, SINGLE, MIXTURE), a total of nine possible combinations of degree distribution and error model. Model selection. The BIC heuristic selects the PL-MIXTURE model for yeast and the TPL-MIXTURE model for worm and fly (Table 2). In general, the PL and TPL models are much better than the ER models. The MIXTURE error model is somewhat better than the SINGLE error model, which in turn is much superior to the SCALED error model. To explore the robustness of this conclusion, we used 10-fold cross-validation (CV) to compute the log-likelihood of the data under each of the models. The CV method identifies the TPL-MIXTURE model as the best for worm and fly, and finds no significant difference between PL-MIXTURE and TPL-MIXTURE for yeast (p-value = 0. 35). Finally, we generated 100 bootstrap replicates each of the yeast, worm, and fly datasets, calculated the BIC scores for each of the nine models, and tabulated the number of times that each model had the best score. The PL-MIXTURE model won 94/100 times for yeast, and the TPL-MIXTURE won 98/100 times for worm and 100/100 times for fly. Model parameters. The PL models yield power law parameters ɛ that are robust to the choice of error model: ɛ = 1. 67–1. 72 for yeast, 1. 48–1. 53 for worm, and 1. 50–1. 58 for fly. In contrast, we can mimic a conventional fit by estimating the PL ɛ using Equation 18 but with k (the corrected number of interaction partners) replaced with w (the observed number of unique partners, which may include false positives). The conventional fit introduces two sources of error: it inflates bait degrees by including false positives, and it deflates bait degrees by excluding false negatives. The conventional estimates for yeast, worm, and fly yield exponents of 2. 22,1. 66, and 1. 61, which are larger (have steeper decay of the degree distribution) than the model results. Thus, the error due to false negatives may dominate the error due to false positives when PL parameters are estimated. Parameter estimates based on prey degree rather than bait degree might be less sensitive to these sources of error. Work by others connects the inverse of c to the typical domain size in a network [15]. The best-fitting TPL models for worm and fly have c ≈ 0. 04 to 0. 07, suggesting a domain size of ten to 30 proteins in a subnetwork. These estimates may not be overly precise as the TPL parameters are sensitive to the error model. The worm data, for example, yields (ɛ, c) = (1. 25,0. 012) for scaled, (0. 46. 0. 035) for single, and (0. 95,0. 040) for mixture error models. The extra variability arises because a larger value of ɛ can compensate for a smaller value of c. The yeast network, which is best fit by a PL-MIXTURE model, yields a very small value of c when fit by any TPL model. This suggests that the yeast network may show less modular structure than either the worm or fly networks. Interaction networks from viruses, parasites, and simpler organisms have been shown to be less clustered than interaction networks from more complex organisms [52]. The fraction of false positives is estimated consistently by the SCALED, SINGLE, and MIXTURE error models regardless of the choice of degree distribution. The false-positive fraction is calculated as ∑if̂i/∑ini from the estimated false-positive count f̂i and number of preys ni for each bait i. This fraction is estimated as 0. 08–0. 09 for yeast, 0. 12 for worm, and 0. 16 for fly. The false-positive fraction is, of course, larger when defined relative to the number of unique interactions identified rather than the number of preys. Although the generative model, the model parameters, and the hidden variables are unknown for the experimental yeast, worm, and fly data, cross-validation is still possible. The cross-validation used half of the dataset to predict the number of new and single interactions in the remainder of the dataset. First, for each bait, we extracted a random half of the preys to serve as a training set. For baits with an odd number n of preys, we selected (n − 1) / 2 and (n + 1) / 2 for the training set with equal probability. Next, we used the training half of the data to estimate model parameters for each organism. For simplicity, we restricted attention to TPL–MIXTURE model. With the maximum likelihood parameter estimates, we then obtained posterior estimates for the false-positive rate and true interaction count of each bait protein. The false-positive rate was determined directly from the mixture model. Based on experience with simulated data, we used exp (〈log k〉) for the posterior estimate of the bait degree. As noted before, the EM equations for PL-like networks, Equation 18, suggest that the logarithm of the degree is a more natural variable than the degree itself (which is not even guaranteed to converge). Finally, we use the statistical model to predict how many unique interaction partners and singleton interaction partners are expected to be observed as the remaining test half of the data is added back. The predictions of the model based on the training half can then be compared with the actual results for the number of unique partners, w, and the number of singleton partners, s. Since the posterior estimate of the bait degree may be fractional, we replaced the factorial function with the Gamma function. Also, rather than starting with the values for s and w for the training set and calculating the marginal increase, we performed a more demanding comparison by using the model parameters to estimate s and w for the observed training set as well. Fitting error was observed only for the count of yeast singletons. Our model can distinguish between random or stochastic false negatives due to undersampling, which could be detected by sampling additional clones, and systematic false negatives that cannot be corrected by deeper sampling of clones from a two-hybrid screen. As described in the section Theory, the overall true-positive rate is the product of the random and systematic rates, psamp × psyst. The true-positive rate from sampling for bait i is (wi − fi) / ki, the number of true positives observed out of the ki true positives represented in the pool. The values of fi and ki are hidden, however, and psamp must be estimated. An appropriate estimator, weighting each true positive equally, is using Equation 19 to define k̂i. The false-negative rate due to systematic losses may be estimated from interactions observed in both directions. We first restrict attention to proteins that were found to have at least one interaction as a bait and as a prey. Using the notation that nij indicates the number of times that bait i recovers prey j, we extract the subset with nij ≥ 2, ensuring that the interaction between i and j is a true positive. These N2 cases are then subdivided into N0 not observed in the reverse orientation, with j as bait and i as prey, and N+ observed at least once in the reverse orientation: The indicator function δ (arg) is 1 for a true argument and 0 for a false argument. The expected value of δ (nji ≥ 1) is the true-positive rate for bait j, equal to the product of the true-positive rates accounting for systematic losses and for random undersampling, Similarly, the expectation of δ (nji ≥ 0) is 1 − psyst (wj − f̂j) /k̂j. The estimated true-positive rate from just the sampling step for the interactions contributing to N2, denoted, is In practice, we find that is somewhat larger than the value of (w − f̂) /k̂ averaged over all baits. For yeast, worm, and fly, the values for are (0. 81,0. 84,0. 81), and the values for averaged over all baits are (0. 58,0. 60,0. 71). Note that, the average of the ratios, is distinct from p̂samp, the ratio of the averages of w − f̂ and k̂ (Table 10). An estimator for N+ in terms of the unknown psyst and other quantities that are known is An obvious route to estimating psyst is to assume that N+ follows a binomial distribution for N2 trials with success rate. The corresponding maximum likelihood estimate p̂syst and its standard error are The estimated number of true interaction partners corrected for sampling losses and systematic losses is then κ̂i, with An alternative estimate for the false-negative rate of a high-throughput screen may be obtained by comparison to literature data. To accomplish this, we downloaded protein interaction datasets from the Database of Interacting Proteins (DIP) [54]. As is common in studies such as this, we required a “gold-standard” dataset that did not use information from the two-hybrid screens we are studying and also was unlikely to be contaminated by false positives. We therefore filtered the entire DIP database to include interactions from only small-scale experiments. Because there is no firm definition of small scale, we used cutoffs of 10 and 100 for the number of interactions reported. For each unique gold-standard interaction, denote the pair of proteins (i, j). We defined bi = 1/0 if protein i was / was not used as a bait in the high-throughput screen, and pi = 1/0 if protein i was / was not used as a prey. We defined similar terms for bj and pj. Also define I (i, j) as 1/0 if the high-throughput screen detected or missed the interaction between i as bait and j as prey. We then calculated the following values: These summary statistics consider the gold-standard interaction in both orientations. For gold-standard interactions between identical proteins, only one of the two identical terms was included in the sum. The bait fraction screened and the prey fraction screened indicate possible correlation in the choice of baits and preys in the small-scale experiments and the high-throughput screens. The summary statistic for the false-negative rate should correct for this bias by only considering interactions that were within the space considered by the high-throughput experiment. This false-negative rate is interpreted as including both the undersampling loss and the systematic loss. The mean number of true interaction partners per bait, corrected for false positives and for undersampling, is Using the definition of the false-negative rate from undersampling, psamp from Equation 25, the mean may also be written as Correcting the mean for systematic losses, psyst from Equation 30, yields This value requires a final correction for the actual search space relative to the entire genome size. The correction we use is the number of preys with at least one interaction, Nprey, relative to the number of proteins annotated in the entire genome, Nproteome. The final proteome-wide mean interaction count is (Nproteome/Nprey) 〈κ̂〉. The number of pairwise interactions in the entire proteome is then An alternative method for estimating interaction counts is to use the inferred probability distribution directly: Pr (k | λ) for Poisson; Pr (k | ɛ) for PL; and Pr (k | ɛ, c) for TPL, Equation 11. For brevity, denote each of these Pr (k | θ̂) for the appropriate parameter estimate θ̂. An overall value for k denoted k̂ may be obtained using an analog of the bait-specific estimate Equation 19, A more typical value for the interaction count corrected for undersampling may be obtained using the median defined by the parametric distribution. We use linear interpolation to estimate a non-integer median, where k+ is the smallest value of k such that CumPr (k) ≥ 0. 5 and CumPr (k) is the cumulative probability Pr (k′ | θ̂) The median interaction count corrected for systematic losses, κ̂med, is then and the median number of interaction partners per protein, corrected for the space screened, is (Nproteome/Nprey) κ̂med.
The genome sequence of an organism provides a parts list of proteins, but not an instruction manual for assembling the parts into a cell. Assembly instructions now come from experiments such as two-hybrid screens that detect physical interactions between pairs of proteins. Defining the resources required for generating a full interaction map requires accurate estimates of the false-negative and false-positive rates of genome-scale screens. Two-hybrid screens often select a query protein and sample its interaction partners. True partners may be missed, and false partners may be spuriously identified. This sampling process resembles a capture–recapture experiment, except that classical capture–recapture theory assumes no false positives. Novel extensions to capture–recapture theory permit its application to proteomic screens. This new theory provides statistically grounded answers to long-standing questions: false-discovery rates of high-throughput screens (possibly over 50% per unique interaction, but probably no more than 15% per clone); the quality of different screening libraries; protein properties leading to “sticky” or “promiscuous” interactions; the global network topology; and, most importantly, the coverage of existing two-hybrid maps. Models estimate roughly 30,000 total pairwise interactions in yeast and 500,000 to 1,000,000 in metazoans. The majority of these interactions remain to be discovered.
mathematics
caenorhabditis
computational biology
drosophila
saccharomyces
Where Have All the Interactions Gone? Estimating the Coverage of Two-Hybrid Protein Interaction Maps