Inferring protein-protein interaction networks from inter-protein sequence co-evolution

Interaction between proteins is a fundamental mechanism that underlies virtually all biological processes. Many important interactions are conserved across a large variety of species. The need to maintain interaction leads to a high degree of co-evolution between residues in the interface between partner proteins. The inference of protein-protein interaction networks from the rapidly growing sequence databases is one of the most formidable tasks in systems biology today. We propose here a novel approach based on the Direct-Coupling Analysis of the co-evolution between inter-protein residue pairs. We use ribosomal and trp operon proteins as test cases: For the small resp. large ribosomal subunit our approach predicts protein-interaction partners at a true-positive rate of 70% resp. 90% within the first 10 predictions, with areas of 0.69 resp. 0.81 under the ROC curves for all predictions. In the trp operon, it assigns the two largest interaction scores to the only two interactions experimentally known. On the level of residue interactions we show that for both the small and the large ribosomal subunit our approach predicts interacting residues in the system with a true positive rate of 60% and 85% in the first 20 predictions. We use artificial data to show that the performance of our approach depends crucially on the size of the joint multiple sequence alignments and analyze how many sequences would be necessary for a perfect prediction if the sequences were sampled from the same model that we use for prediction. Given the performance of our approach on the test data we speculate that it can be used to detect new interactions, especially in the light of the rapid growth of available sequence data.


Introduction
: Contact map and protein-protein interaction network of small and large ribosomal subunits. The contact map and the protein-protein interaction network for A the small ribosomal subunit and B the large ribosomal subunit (proteins only), using a distance cutoff of 8Å between heavy atoms. The upper diagonal part shows the contact map, with red dots indicating intra-protein contacts, and blue dots inter-protein contacts. The lower triangular part shows the coarse graining into the corresponding protein-protein interaction networks, with the color levels indicating the number of intraresp. inter-protein contacts, cf. the scales. The sparse character of both the contact network and the interaction network is clearly visible.
ogy (residue contact map inference) to systems biology (PPI network inference). An obvious problem in this context is the sparsity of PPI networks, illustrated by the bacterial ribosomal subunits used in the following, cf. Fig. 1: The small subunit contains 20 proteins and 21 protein-protein interfaces (11% of all 190 possible pairs). In the large subunit, 29 proteins form 29 interfaces (7% of all 406 pairs). We see that while the number of potential PPI between N proteins is N 2 , the number of real PPI grows only linearly as O(N ). Furthermore, the number of potentially co-evolving residueresidue contacts across interfaces is much smaller than the number of intra-protein contacts. In the case of ribosomes, only 5.8% of all contacts in the small subunit are inter-protein contacts. In the large subunit this fraction drops down to 4.5%. So the larger the number of proteins, the more our problem resembles the famous search of a needle in a haystack. The noise present in the large number of non-interacting protein family pairs might exceed the co-evolutionary signal of interacting pairs. It should also be mentioned that the ribosomal structure relies on the existence of ribosomal RNA, which is not included in our analysis. We therefore expect many of the small PPI interfaces to be of little importance for the ribosomal stability and that only large interfaces constrain sequence evolution and thus become detectable by co-evolutionary studies.
Ribosomal proteins and their interactions are essential and thus conserved across all bacteria, and it appears reasonable to wonder whether this makes them a specialized example of a protein complex more amenable to co-evolutionary bias. As a second and smaller interaction network, we therefore considered the enzymes of the tryptophan biosynthesis pathway comprising a set of seven proteins in which only two pairs are known to interact (PDB-ID 1qdl for the TrpE-TrpG complex [27] and 1k7f for the TrpA-TrpB complex [40]). Also here the PPI network is very sparse; most pairs are not known to interact, but might show some degree of coordinated evolution due to the fact that in many organisms these genes show a common spatial co-localization in a single operon and also due to a number of gene fusion events, cf. the discussion below. While widespread, the tryptophan biosynthesis pathway is not essential for viability when environmental tryptophan is present.
In this paper we report the excellent performance of DCA in the prediction of protein-protein interaction partners in the systems tested. In a first step, we analyze the performance on data from an artificial model. This allows for a systematic analysis of the performance of different approaches and of the influence of the number of sequences in the alignment. With this artificial data set we are able to establish a lower-bound on the number of sequences that would make our predictions on the PPI scale completely accurate if the generating model was the same model we use for inference. Given the growth-rate of current protein sequence databases (notably UniProt [8]), we expect that such a lower bound could be met in few years. In a second step, we apply the method to the proteins of the bacterial ribosome and to the proteins of the trp operon, and show that the results obtained for simulated data translate well to the biological sequences of this test-set.

Materials and Methods
The goal of the present work is to analyze each of the N 2 possible pairs of multiple sequence alignments from a given set of N single-protein family alignments, and to extract a pairwise score that measures the co-evolution between the proteins in the alignments. A high co-evolutionary score is then taken as a proxy for interaction. In the spirit of [14] we describe in this section consecutively the data generation and matching, the model used for analyzing data and the inference and scoring mechanism.

Data extraction and matching for the ribosomal and trp operon proteins
The input data is given by N multiple sequence alignments D p consisting of M p sequences of length L p for every protein family p. These alignments are extracted from UniProt [8] using standard bioinformatics tools, in particular Mafft [25] and HM-Mer [16] (cf. Section 1 in S1 Text for details on the extraction pipeline and Tables S1,S2 in S1 Text listing the values of N , M p , L p for ribosomal and trp-operon proteins). For the analysis, it is necessary to concatenate the MSAs of two putative co-evolving protein families. This means to create, for each pair of protein families (p, p ), a new alignment D p,p of sequence length L p + L p . Each line contains the concatenation of two potentially interacting proteins. More precisely, in the case where families p and p actually interact, each line should contain a pair of interacting pro- Sketch of the matching procedure that allows us to concatenate two different MSAs, here MSA 1 ,MSA 2 . π represents the optimal permutation of the sequences on the second MSA computed using a standard linear programming routine.
teins. The general problem of producing a concatenated alignment out of single MSAs of two protein families is straightforward in two cases only: (i) we have prior knowledge which pairs of sequences represent interaction partners; (ii) no paralogs are present in the considered species (i.e. all species have at most a single homolog of each of the sequences to be matched). Often, as displayed schematically in Fig. 2, MSAs contain multiple protein sequences within a given species and no prior knowledge can be used to know who is (potentially) interacting with whom. In prokaryotes, interacting proteins are frequently found to be coded in joint operons. This suggests to use genomic co-localization as a matching criterion. To do so, as explained in Section 2 in S1 Text, we approximated the genomic distance between sequences using UniProt accession numbers. A better distance between sequences could be defined in terms of their genomic location. Unfortunately, genomic locations are available only in the context of whole genome sequencing projects. The majority of sequences in Uniprot originate from fragments or from incomplete genome sequencing projects. These difficulties lead us to content ourselves with the proxy of accession numbers. Having defined distances between each protein pair in the MSA, we calculate the matching which minimizes the average distance between matched sequences by linear programming. Additionally, we introduce a distance threshold used to discard matched distal protein sequence pairs. The numeric value for this threshold was determined using the small ribosomal subunit as a test case.
The average number of paralogs per species varies from system to system: For both ribosomal subunits the proteins have between 1.5 and 3 paralogous sequences per genome. The trp proteins on the other hand have considerably more paralogous sequences and the number of such sequences per genome varies between 4 and 24. This means that especially in the trp operon the matching procedure has the potential to generate much larger alignments than the competing approach of excluding species with paralogous sequences. In fact, using this last approach (which corresponds to setting our threshold parameter to 0) reduces the number of sequences in the alignments on the average by about 10% for the ribosomal proteins and by about 85% for the proteins of the trp operon (see Tables S3-S7 in S1 Text ).
Note that using paralogs may be dangerous since after duplication different paralogs often evolve different functions, and thus loose part of their interactions or gain others. However, our matching strategy based on genomic vicinity excludes proteins coming from isolated genes; it identifies mostly protein pairs coded in gene pairs colocalized inside an operon. It is therefore more likely that the two maintained interaction, when also the ancestral protein pair before duplication was interacting. We will show evidence that, in the interacting protein systems investigated here, this strategy leads to a reinforced coevolutionary signal. However, an independent and direct test whether protein pairs included in the alignment actually interact would constitute a big step forward.
Let us recall that the problem of finding a good matching between sequences has already been studied in the past using different strategies [6,35]. Unfortunately, both methods are computationally too demanding to be used in a case, where hundreds or thousands of protein family pairs have to be matched.

Statistical sequence model
Within DCA, the probability distribution over amino acid sequences x = (x 1 , ..., x L ) of (aligned) length L is modeled by a so-called Potts model, or pairwise Markov Random Field, which includes statistical couplings J ij (x i , x j ) between residue pairs and positionspecific biases h i (x i ) of amino-acid usage [38]. The number Z is the normalization constant of P(x), which is a probability distribution over all amino-acid sequences of length L. The variable x i represents the amino acid found at position i in the sequence and can take as values any of the q = 21 different possible letters in an MSA (gaps are treated as a 21st amino acid). The model parameters are inferred using MSAs of homologous proteins. In the case of two concatenated protein sequence (x, x ) = (x 1 , ..., x L , x 1 , ..., x L ), the joint probability takes the form The functions H(x) and H (x ) are the terms in the exponential in Eq. (1) referring to each single protein. The function describes the co-evolutionary coupling between the two protein families. In the last expression, x i is the ith amino acid in sequence x, and x j the jth amino acid in sequence x . The sum runs over all inter-protein pairs of residue positions. The q×q matrices J ij in this term quantify how strongly sites between the two proteins co-evolve in order to maintain their physicochemical compatibility. The matrix contains a real number for each possible amino acid combination at sites i and j and contributes to the probability in Equation 2 depending on whether an amino acid combination is favorable or not. The strongest inter-protein couplings are enriched for inter-protein contacts [32,38].
The same kind of model can be used to predict the interaction between more than two proteins, with a corresponding number of interaction terms. However, the number of parameters in the model is proportional to (L 1 + L 2 + .. + L N ) 2 for N proteins while the number of samples in the concatenated MSA D p1,....p N becomes smaller because one has to find matching sequences for N proteins simultaneously. This leads us to consider the case N > 2 only for artificial proteins where the total length and sample size are controllable.

Inference and Scoring
Following [12], the parameters of the model were inferred by maximizing pseudolikelihood functions. This is an alternative to directly maximizing the likelihood and considerably faster (see Section 3 in S1 Text Text for details). Given that the model is mathematically equivalent to the one used in [12] we can use the output of the algorithm (plmDCA) with default parameters as presented there directly for our purposes. This output consists of scores F ij (the average-product corrected Frobenius norm of the matrices J ij ) that quantify the amount of co-evolution between sites i and j in the alignments. In order to quantify co-evolution between proteins, we took the F ij corresponding to inter-protein site pairs (i.e. i in x and j in x ) and calculated the mean of the 4 largest. These quantities, a real number for every protein pair, are used to rank protein-protein interaction partners. The number 4 was chosen because it performed well in the small ribosomal subunit, which we used as a test case when designing the algorithm. Subsequent tests on larger systems showed that any number between 1 and 6 performs almost equally well (see Section A Global View in Results).

Simulated data
As the basis for the simulated data we used a fictitious protein complex consisting of 5 proteins. Each protein has a length of 53 residues. The individual contact map of each one is given by the bovine pancreatic trypsin inhibitor (PDB ID 5pti [42]), which is a small protein performing well for the prediction of internal contacts by DCA.
Each P i has 551 internal contacts. Moreover, each protein interacts with two others in a circular way. The inter-protein contact matrices between P i and P i+1 (as well as between P 1 and P 5 ) are random binary matrices with a density of 10% of the internal contacts. This models the sparsity of the inter-protein contacts as compared to the intra-protein contacts. A contact map for the artificial complex can be found in Figure  S5 in S1 Text, There are no contacts between other pairs of proteins.
In order to define as realistically as possible the coupling parameters of the Potts model used for generating the artificial sequences, we used the Pfam protein family PF00014 of the pancreatic trypsin inhibitor [15]. Note that a member of this family was also used to define the structure. The couplings describing the co-evolution within the single proteins were directly extracted from the Pfam MSA using DCA. For the couplings corresponding to the co-evolution between the proteins, we used a random subset of the internal parameters and used them to couple sites that are in contact according the contact map as defined above. Non-contacting pairs of sites remain uncoupled between artificial proteins. Using this model, a joint MSA D 12345 of sequences of length 265 = 5 × 53 was generated using standard MC simulations.
The process of defining the contact map, choosing the parameters and generating the sequences is described in Section 5 S1 Text.

Results and Discussion
Testing the approach using simulated data As a first test of our approach, we use simulated data generated by Monte Carlo (MC) sampling of a Potts model of the form of Eq. (2), cf. Materials and Methods.
The main simplifying assumptions in this context are: (i) We assume intra-and inter-protein co-evolution strengths to be the same. (ii) We assume the distribution of inter-protein residues contacts within the possible contacts to be random. (iii) We assume the sequences to be identically and independently distributed according to our model. This model includes the assumption that non-contacting sites have zero couplings. The number of artificial sequences needed for a good performance of our method should therefore be taken at most as a lower bound for the number of biological sequences needed for a comparable performance.
In panel A of Fig. 3 we show the architecture of our artificial protein complex. It is composed of five fictitious, structurally identical proteins P 1 , ..., P 5 , each one consisting of 53 residues. In order to simulate co-evolution between the proteins, we generate a joint MSA D 12345 for all 5 proteins with a model that contains couplings between inter-protein site pairs. These couplings are modeled in a way to resemble couplings inferred from real proteins (see Materials and Methods).
To assess our capability to infer the PPI network of panel A from such data, we adopted two different strategies which we called combined and paired in panel B of Fig. 3. The combined strategy uses plmDCA on the full-length alignments of length 265 and models the interaction between all proteins pairs simultaneously. Given that in this artificial setting we use the same model to generate the data as to analyze it, the approach is guaranteed to infer the model correctly for a large number of analyzed sequences and therefore to assign a higher interaction score to any interacting protein pair than to any non-interacting pair.
To assess the coupling strength between two proteins, we average the four strongest residue coupling strengths between them. This leads to a score oriented toward the strongest signal while also reducing noise by averaging. In panel B of Fig , 000 is what we expect to be available in a few years from now, seen the explosive growth of sequence databases. The thickness of each link in Fig. 3 is proportional to the inferred inter-protein interaction score. The five strongest links are colored in green when they correspond to actual PPI according to panel A, and in red when they correspond to non-interacting pairs. For increasing sample size the predictions become more consistent and for M = 24, 000 any interacting protein-pair has a higher interaction score than any non-interacting pair.
Due to the running time of plmDCA only alignments for sequences of total length L 1000 can be analyzed. This is exceeded already by the sum of the lengths of the proteins of the small ribosomal subunit. Additionally, creating a combined multiple sequence alignment for more than two proteins would lead to very low sequence numbers due to the necessary matching (see Materials and Methods). Therefore, using the combined strategy is not generally applicable. In the paired strategy we therefore analyze each pair of proteins separately. This means that plmDCA is applied to all Fig. 3 we find that the paired strategy is also able to detect the correct PPI network for large enough M . We observe, however, that the performance of the paired strategy is slightly worse. Couplings between non-interacting proteins are estimated significantly larger than using the combined strategy for large M . Even in the limit M → ∞ we do not expect these links to disappear: Correlations between, e.g., P 1 and P 3 are generated via the paths 1 − 2 − 3 and 1 − 5 − 4 − 3, but in the paired strategy these correlations have to be modeled by direct couplings between P 1 and P 3 since the real direct coupling paths are not contained in the data.
After having answered the 'who-with-whom' question for the artificial protein network, we address the 'how' question of finding inter-protein contact pairs. Fig. 4 panel A displays individual residue contact pairs within and between proteins in the artificial complex. Panel B shows the 10 strongest intra-protein couplings for each protein and the 10 strongest inter-protein couplings inferred by plmDCA (M = 4000, combined strategy). Green links correspond to contact pairs and red links to non-contact pairs. We see that the intra-protein prediction is perfect, whereas a few errors appear for inter-protein predictions in agreement with the results of Fig. 3.

The PPI network of bacterial ribosomes
As a more realistic test we apply the method to the bacterial large and small ribosomal subunits (LRU, SRU). To define contacts and protein interaction partners we used high-resolution crystal structures with PDB-IDs 2z4k (SRU) and 2z4l (LRU) [4]. The contact network is summarized by the contact maps in Fig. 1. The ribosomal RNA is ignored in our analysis.
Panels C, E of  (1,439) inter-protein contacts, so the contacts relevant for our study comprise only 5.8% (4.5%) of all contacts. Fig. 3 panel D shows the inferred SRU PPI architecture. As expected, the biological case is harder than the artificial case where the data are independently and identically distributed according to the generating model. Even though the histograms of the inferred interaction scores for both cases are very similar (see Figure S2 in S1 Text, biological data are expected to show non-functional correlations due to the effect of phylogeny or sequencing efforts which are biased to model species and known pathogens. Nonetheless, among the top ten predicted interacting protein pairs the method makes only three errors (true-positive rate 70% as compared to 21/190 11% true PPI between all protein pairs, with an overall area under ROC curve (AUC) of 0.69 (see Fig. 7). The method spots correctly the pairs with larger interaction surfaces whereas the small ones are lost. Two of the false-positive (FP) predictions include protein RS21, which has the smallest paired alignments with other proteins (M between 1468 and 1931). Also the third FP, corresponding to the pair RS4-RS18, is probably due to a small MSA with M = 2064. At the same time, the interaction of RS21 with RS11, which is one of the largest interfaces (199 contacts), is still detected despite the low M = 1729. The same procedure for the LRU (406 protein pairs) performs even better: 9 out of the 10 first PPI predictions are correct (see Fig. 3 panel F), and the AUC is 0.81.
The results on the residue scale for both SRU and LRU are depicted in panels D and F of Fig. 4. Shown are the first 20 intra-protein residue contact predictions for each protein (excluding contacts with linear sequence separations below 5 to concentrate on non-trivial predictions) and the first 20 inter-protein residue contact predictions. In the SRU case of panel D for example, the results are qualitatively similar to the artificial case, albeit with a slightly reduced true-positive rate of 60% among the first 20 inter-protein residue contact predictions (compared to the ratio of 1401 actual inter-protein residue contacts and 2,403,992 possible inter-protein residue contacts, i.e., 0.058%). Again 3 out of the 8 false positives are related to RS21, which due to the smaller MSA size is also the only one having a considerable false-positive rate in the intra-protein residue contact prediction. About 95% of the displayed 400 highest intraprotein residue contacts are actually contacts (see Figure 1 in S1 Text). Analogous considerations with a somewhat larger accuracy (85%) hold for LRU as displayed in Fig. 4

panel F.
The PPI network of the tryptophan biosynthetic pathway As a distinct test case for our methodology we analyzed the 7 enzymes (TrpA, B,C,D,E,F,G) that comprise the well characterized tryptophan biosynthesis pathway. In contrast to the ribosomal proteins, these enzymes are only conditionally essential in the absence of environmental tryptophan and their genes are only expressed under deplete tryptophan conditions. In this particular system, only two protein-protein interactions are known and resolved structurally: TrpA-TrpB (PDB-ID 1k7f [40]) and TrpG-TrpE (PDB-ID 1qdl [27]). Whereas the TrpG-TrpE pair catalyzes a single step in the pathway and their interaction is thus essential for correct functioning, the TrpA-TrpB pair catalyzes the last two steps in tryptophan biosynthesis. Both enzymes function in isolation but their interactions are known to increase substrate affinity and reaction velocity by up to two orders of magnitude. All other proteins catalyze individual reactions, but one might speculate that the efficiency of the pathway could benefit from co-localization of enzymes involved in subsequent reactions. Interestingly, the Pfam database [15] reports that in many species pairs of genes in the operon appear to be fused, suggesting that some of the fused pairs are actually PPI candidates. An example is the TrpCF protein, which is fused in Escherichia coli and related species (but not in the majority of species).
After applying our method to all 21 protein pairs we find elevated interaction scores only for TrpA-TrpB and TrpE-TrpG, which are the only known interacting pairs (see Fig. 5 and Table S10 in S1 Text for the interaction scores of all pairs). Those two pairs have interaction scores of 0.375 and 0.295, while the other pairs are distributed between 0.071 and 0.167. Even though we do not define a significance threshold for prediction (see Section A global view ), these two pairs would be discernible as interesting candidates even if we did not have the 3D structures.
We speculate therefore that the fusions in many species do not imply strong interprotein co-evolution. To further investigate this aspect, we took a closer look at the protein pair TrpC-TrpF. For this protein pair, a high resolution structure of a fused version exists (PDB-ID 1pii [41]). We ran our algorithm on the complete multiple sequence alignment, the multiple sequence alignment with fused sequence pairs removed and only on the fused sequences. In none of these cases did we observe a statistically significant interaction score or a statistically significant prediction of inter-protein contacts present in the structure of the fused protein.
Our results are corroborated by the finding that all scores measuring the co-evolution between a ribosomal protein and an enzyme from the tryptophan synthesis pathway are small (see the following subsection). No indication for an interaction between the two systems is found, as to be expected from the disjoint functions of the two systems.

A global view
It is interesting to assemble a larger-scale system out of the three systems (SRU, LRU, Trp). To this end, we created all possible pairings between the proteins used in the present study (SRU vs. RU, SRU vs. Trp, LRU vs Trp, SRU vs SRU, LRU vs. LRU, and Trp vs. Trp). This leads to a total of 1540 pairs, out of which only 49 pairs are known to interact (which we defined as true positives). We present the findings in Fig. 7 and in Figs. S7-S9 in S1 Text. Fig. S7 in S1 Text shows the true-negative rate, which is the fraction of true negatives in the indicated number of predictions with the lowest interaction scores. As it can be seen our scoring produces a false negative just after 420 true negatives. Figs. 7 and S8 show true positive rates for the complete system and the individual systems. We also show true positive rates for alternative ways to calculate the interaction score between protein pairs, i.e. a different number of inter-protein residue-residue interaction scores to average. We notice that in the complete system, the performance is similar to the performance in individual systems. All of the 10 highest-scoring protein pairs are known to interact, and 75% of the first 20 protein-pairs. After these first 20 pairs, the true positive rate drops to around 45% in the first 40 predictions. This is analogous to the case of protein contact prediction, where methods based on the same model are able to extract a number of high confidence contacts but see a large drop in performance afterwards [31]. The area under the AUC for the whole system is 0.83 (see Fig. 7). This is stable when averaging different numbers of residue-contact scores to arrive at a protein-protein interaction score, but the performance seems to worsen when using more than 6. This is probably because only a few inter-protein residue contacts have a large score and averaging over too many only adds noise. It can also be seen that averaging over 4 performs very well in the small ribosomal subunit, which is why we have chosen this value for the large part of the analysis. On the larger-scale system, though, any number between 1 and 6 performs almost identically.
A further question is whether it is possible to define a threshold allowing to reliably discriminate between interacting and non interacting pairs in terms of the interaction score. Figure S9 in S1 Text shows two normalized histograms of the interaction scores. The rightmost tail of the interacting pairs distribution is well separated from the rightmost tail of non-interacting one, but the highest scores of non-interacting pairs are strongly overlapping with the lowest scores of the interacting ones. The situation is therefore analogous to what is observed in the case of the inference of contacts within single protein families [2,12,31,38], where the same technique is known to produce relatively few high confidence contacts in the topmost scoring residue pairs. To conclude, while high scores seem to reliably predict interacting pairs, and low scores non-interacting pairs, there is a large gray zone prohibiting a clear discrimination between interacting and non-interacting pairs.  Shown are means for all protein pairs that have at least 100 residue pairs in contact. The ribosomal and the trp proteins were tested independently. The red curves correspond to a matching including only protein sequences without paralogs inside the same species ("matching by uniqueness in genome"). The low performance of this approach on Trp proteins is due to a very low number of species without homologs, which leads to very small matched alignments. The blue curves show the results for our matching procedure as described in the text. The green curves correspond to alignments that have been obtained by first applying our matching procedure and then randomizing the matching within individual species. The definition of "contact" was the same as used above (a distance of less than 8.0Åbetween two heavy atom in the residues).

Conclusions
To conclude, we have shown that DCA performs excellently in the systems tested when used to predict protein-protein interaction partners. In the small and large ribosomal subunit our tests resulted in a true positive rate of 70% and 90% in the first 10 predictions (AUC of 0.69 and 0.81) while in the trp operon the two largest interaction scores corresponded to the only two interactions experimentally known (AUC 1). The performance is summarized in Fig. 7. The Figure shows both the high quality of the first predictions, but also a drop in performance after a fraction of all interacting pairs (about 40% in our test case). This is analogous to the case of protein contact prediction by DCA and related methods, where the performance drops after a limited number of high-confidence predictions [31]. In the same context and with the same caveat, an excellent performance in predicting inter-protein contacts on the residue level has been shown. The artificial data have shown that the performance of our approach depends crucially on the size of the alignments. Only for very large MSA (M = 24, 000 sequences in our data) a perfect inference of the artificial PPI network was achieved. MSA for real proteins pairs are typically much smaller. Even for pairs of ribosomal proteins, which exist in all bacterial genomes, only about 1500-3200 sequence pairs could be recovered. This places these data towards the lower detection threshold of PPI. We therefore expect the performance of the presented approach to improve in the near future thanks to the ongoing sequencing efforts (the number of sequence entries in Uniprot [8] has been growing from about 10 millions in 2010 to 90 millions in early 2015) and improved inference schemes. The strong performance of the same algorithm on different and dissimilar systems naturally prompts us to expect that the approach can be used to detect interactions experimentally unknown so far. In fact, if we trust our results on the trp operon we can already draw some speculative biological inferences. While there are many high-resolution structures of the ribosome available, one might have expected that in the trp operon there could be more transient previously unreported interactions in the tryptophan biosynthesis pathway beyond the two interactions that have been structurally characterized. As mentioned, various enzyme fusions can be observed in the databases, suggesting that there is an evolutionary benefit to co-localizing the enzymes of the pathway in the cell. An obvious benefit of such co-localization would be that the pathway intermediates do not have to diffuse throughout the cell from one enzyme component to the next. In the tryptophan biosynthesis pathway in particular, there are numerous phosphorylated intermediates that need to be protected from unspecific cellular phosphatase activity. Organizing the enzymes in the pathway in a multi-protein complex would seem like an efficient way to protect the intermediates from decay. However, our data indicate that the only statistically relevant co-evolutionary signals that can be observed are restricted to the known strong interactions between TrpA with TrpB and TrpE with TrpG. This could be interpreted in a number of ways: (i) The most obvious explanation is that there are no additional protein-protein interactions beyond those that are known and that no multi-enzyme complex exists for the tryptophan biosynthesis pathway. Alternatively (ii) it seems plausible that there are numerous structural solutions to form a tryptophan biosynthesis complex and that there is no dominant structure from which a co-evolutionary pattern can be observed in the sequence databases. Lastly (iii) it is not out of question that the enzymes of the pathway do not directly form a complex but that they are jointly interacting with an unidentified scaffold component. Of course we cannot exclude that our method is not able to capture other potentially present interactions. From a methodological point of view, one possible algorithmic improvement is creating better MSAs for protein pairs. The vast majority of protein families show genomic amplification within species. This raises the issue of which sequence in one MSA should be matched with which sequence in the other MSA when concatenating the two MSAs, as shown in Fig. 2. In the absence of prior knowledge and as long as only prokaryotes are concerned, we showed that it is possible to use the simple criterion of matching by genomic proximity. This criterion is based on the observation that two sequences are more likely to interact if they are genomically co-localized. Our results have shown that in the case of the ribosomal network better inference results can be obtained by using this matching criterion than by using a random matching or using a conservative matching taking only species with a single sequence in both MSAs into account, cf. Fig. 6. However, we found it beneficial for the predictive performance to introduce a threshold distance above which we simply discarded candidate sequences. This is not based on biological principles.
We believe that our naive matching strategy can be improved substantially. Even if closeness of sequence pairs on the genome is a good proxy for interaction in some cases, for example if they belong to the same operon, excluding all distal pairs is a very crude criterion. This criterion is known to be erroneous in many cases, for example in the bacterial two component signal transduction system [6,7,35]. It would therefore be interesting to include the matching into the inference procedure itself, e.g. to find a matching that maximizes the inter-protein sequence covariation, cf. [6] for a related idea. However, for highly amplified protein families this leads to a computationally hard optimization task. Simple implementations get stuck in local minima and do not lead to improvements over the simple and straight-forward scheme proposed here.
Supporting Informations S1 Text 1 Multiple Sequence Alignments

Multiple Sequence Alignments
The data we use are multiple sequence alignments (MSA). Each such MSA is a rectangular matrix, with entries coming from a 21-letter alphabet containing the 20 standard amino acids and a gap symbol "-". In the following we denote this alignment by a matrix X = (x a i ) , i = 1, ..., L, a = 1, ..., M with L being the number of residues of each MSA row, i.e., the number of residues in each considered protein, and M the number of MSA rows, i.e., the number of proteins collected in the alignment. For simplicity of notation we assume that the 21 amino acids are translated into consecutive numbers 1,...,21.

Alignment Generation
For all proteins of the small ribosomal subunit (SRU) and the large ribosomal subunit (LRU) the sequence names were extracted from the corresponding PFAM alignments [15]. Using these names, the following procedure was used to create the alignments for the single proteins: 1. Extract sequences corresponding to names from Uniprot [8] 2. Run MAFFT [26] on them using mafft --anysymbol --auto 3. Remove columns from the alignment that contain more than 80% gaps 4. Create an Hidden Markov Model (HMM) using hmmbuild from the hmmer suite [17] 5. Search Uniprot using hmmsearch [17] 6. Remove inserts 7. If there exist in one species two or more sequences that are more than 95% identical, remove all but one.
The number of sequences for the single files can be found in Table 1 The alignments for the proteins of the Trp Operon where constructed analogously with some modifications to ensure that only full-length sequences where extracted. Also, we chose the linsi program of the MAFFT package to create the initial MSAs. The number of sequences for the Trp alignments can be found in Table 2

Internal Sensitivity Plots
As an assessment of quality for the alignments, sensitivity plots using the pdb files 2Z4K and 2Z4L were made. Figure 8 shows results for contact predictions based on the GaussDCA [2] and plmDCA alghorithm [13].  2 Matching Procedure

Pipeline for Matching
The problem of generating a concatenated alignment from two MSAs of two different protein families (say MSA 1 and MSA 2 ) is to decide which sequence from the first alignment should be concatenated to which sequence from the other alignment. This means to find for any protein p 1 i in MSA 1 a matching partner p 2 j in MSA 2 belonging to the same species. The problem is trivially solved in the case when no paralogs are present and each species has one and only one sequence in each individual MSA. In this case we can simply concatenate these two sequences (we term this case matching by uniqueness). The problem is that species often have several paralogs. In this case, given that we would like to observe a co-evolutionary signal between protein interaction partners, one would like to match sequences of proteins that are (possibly) interacting.
As long as Prokaryotes are concerned, it turns out empirically that proteins are more likely to interact if their genes are co-localized on the DNA [6,38]. This suggests to try to match proteins that are close on the genome when creating a concatenated MSA.
As a proxy to the genomic distance we use a distance between Uniprot accession numbers (UAN). This UAN consists of a 6 digit alphanumeric sequence for every sequence and can be extracted from the sequence annotation, e.g. the "D8UHT6" part of the sequence annotation "D8UHT6_PANSA".
We define the distance between UANs as follows: Different positions in the UAN can take on different values, some only numeric (0-9) and some alphanumeric values (0-9,A-Z). We define for every position i ∈ 1 . . . We further map the possible single position values in the UAN to the natural numbers in ascending order, i.e. we assign to the numeric symbols 0−9 the natural numbers 0 − 9 and to the letters the natural numbers following 9 (so to A we assign 10, to B we assign 11 etc.). This leads for example for the the UAN L9XG27 to the numeric sequence A = (21,9,33,16,2,7). Now we can define a unique number N for any UAN that has been mapped to the sequence of natural numbers A i as The distance between two UANs that have been mapped to the numbers N 1 and N 2 can now be defined as This procedure induces a distance D ij for any sequence p i ∈ MSA 1 and p j ∈ MSA 2 , where both p i , p j belong to the same species. In this way we define a complete weighted bipartite graph, and the problem of finding the proper pairing can thus be translated into a minimum weighted bipartite matching problem. This problem can be readily solved using a standard linear programming techniques. Finally we discard from the optimal solution sequence pairs whose distance is above a given threshold of 100 (manually optimized on the small ribosomal subunit). In the cases we analyzed, such a threshold moderately increases the quality of the prediction of interaction partners.

Inference technique
As a simple but meaningful statistical model, we consider a pairwise generalized 21 states (to mimic the 20 amino acids + 1 insert symbol alphabet of MSAs) Potts model with the following Hamiltonian We can now assume to have a dataset D = {x 1 , . . . , x M }, where x represents one sequence, either artificially generated, or extracted using the bioinformatic pipeline discussed above. Notice that if the sequences x are concatenations of two sequences (x, x ), the sums in Equation 7 can be split into three parts: One in which appear only sites in x, one in which appear only sites in x and one interaction part with J ij for which i is in x and j in x . By labeling the first part H(x), the second H (x ) and the third H int (x, x ) one arrives at the representation referred to in the main text. Given that the representations are mathematically equivalent, we will here in supplemental information treat the sequence as one simple sequence x.
The inference proceeds by assuming as a working hypothesis that the dataset D is composed by configuration sampled uniformly from the equilibrium Boltzmann-Gibbs distribution P ( x) = exp(−H)/Z (as an inference process, we are free to consider T = β = 1). We are now ready to use D to infer the topology of the network. To do so -as discussed in the main text -in the last years different maximum-likelihood techniques have been proposed [1,12,23,29,31,39]. So far the most promising in terms of accuracy seems to be the pseudo-likelihood maximization introduced in [12] where from the previously defined Boltzmann-Gibbs measure we consider the following conditional probability distribution: Given a data set D we can thus maximize the conditional likelihood by maximizing as a function of J i,\i , h i . As customary in many maximum-likelihood inference techniques, we add to the maximization an L 2 regularization term, so that eventually the extremization procedure turns out to be: with J ij 2 = 21 a,b=1 J 2 ij (a, b), and h i 2 = In order to produce an interaction score for the two proteins, we run the PLM algorithm [12] on the concatenated alignments. This results in a list of residue pairs of the alignment ordered by their interaction strength. We filtered out the pairs that contain one residue of one protein and one of the other. This results in a list of possibly interacting inter-protein residue pairs ordered by the interaction score. In order to arrive at an interaction score for the two proteins we took the mean of the scores for the 4 highest scoring pairs (PPI-score). The number 4 was used because it performed best on the small ribosomal subunit, but the predictive performance on a larger-scale network is virtually identical for any value between 1 and 6 (see Figure  S8). The list of protein pairs ordered by this score was used for prediction. The first few predictions are shown in Table 8. For completeness, we show the same table but with the score calculated by the Gaussian approximation of [2] in Table 9. Finally in Table 11 we display for the LSU the number of intra/inter-protein contacts, while in Table 12 we do the same for the LRU. Table 10 shows the interaction scores for the protein pairs of the Trp Operon.

Artificial Data
An artificial large network consisting of 5 proteins was created in two steps: 1) First, a contact map was defined. This contact map contains the information which residues are in contact. This includes internal residue contacts (where both residues belong to one of the 5 proteins) and inter-protein residue contacts (where one residue belongs to one protein and the other to a different protein). The contact map is therefore a binary, symmetric matrix of size N all × N all with N all = N 1 + N 2 + N 3 + N 4 + N 5 where N i is the number of residues in the i th protein. We decided to use the Kunitz domain (PF00014) as a model for the proteins and set all N i = 53. The 53 × 53 submatrices that define the contacts within each protein were defined by extracting the contacts of the PDB structure 5pti of the Kunitz domain. This implies that the internal structure of every protein is the same.
We defined as contacting proteins the protein pairs 1 − 2, 2 − 3, 3 − 4, 4 − 5 and 1 − 5. For the 53 × 53 submatrices that define the contacts between contacting protein pairs we used random binary matrices with 10% of the number of internal contacts. This was done individually for each contacting protein pair such that no two contact matrices between two proteins were the same. For non-contacting protein pairs all entries of the contact matrices were set to 0.
The resulting contact map can be seen in Fig. 12.
2) Couplings for every contact in the contact map were defined. As a basis for this, couplings and fields inferred from the PF00014 PFAM alignment (Kunitz Domain) were used. This inference was done using a masking with the PDB structure, such that only couplings corresponding to PDB-contacts were allowed to differ from zero. Given that the same PDB-contacts were used to define the contacts within one protein in the artificial complex, we could use the couplings thus inferred without change for the couplings within the artificial proteins.
Then we defined the couplings for residue contacts between two proteins. For every such a resiue contact we chose randomly a coupling of an internal contact as inferred from the Kunitz domain alignment and assigned it to the residue contact.
Notice that the 'coupling' between two sites i and j is actually a 21 × 21 matrix J ij (a, b) where a and b can be any of the 21 amino acids. Given that the internal structure of these matrices might be important we decided to treat the matrices J ij as single entities and not change their internal structure.
The fields for every residue, a vector of length 21 for every of the 5 · 53 residues, were randomly chosen from the inferred fields.
From these couplings and fields, sequences were generated by MC (see section below) and inferred by plmDCA. Interestingly, a crude comparison between the histogram of the scores in the artificial model seem to be very close to that obtained for instance for the LRU case as shown in Fig. 9.
In Table 13 we compare the ranks of the strongest inter-protein residue interaction scores in the generating model and the inferred model. The first column represents the rank of the inter-protein residue interaction in the generating model, the second column the rank of the same residue interaction in the inferred model. The model was inferred with the combined strategy and with 4000 sequences. The numbering is  Inferred Rank  1  101  2  13806  3  10658  4  64  5  4  6  9575  7  1  8  15890  9  6712  10  1035  7  1  32  2  41  3  5  4  11  5  11473  6  22464  7  53  8  1877  9  26 10 Table 13: Original vs. inferred rank for the 10 largest original inter-protein residue interaction scores and the 10 largest inferred inter-protein residue interaction scores treating the complex as one large protein.

Monte Carlo Sequence Generation
Given the parameters of the artificial model, a simple MCMC algorithm was run to generate samples from the corresponding distribution. We used one million MC steps to equilibrate the chain and took a sample every one million steps.  Figure 13: Inferred protein network for different sample sizes; the line-thickness is proportional to the inferred interaction scores between the proteins (mean of the 4 highest residue interaction scores). The thickness has been normalized in the sense that the scores have been divided by the mean of the scores of the network. The color code is applied for the first 5 predictions and shows a green line if the prediction is a true positive and a red line if the prediction is a false positive. Predictions after the first 5 are grey. Combined Analysis: The complete sequences in their whole length were used for the inference and calculation of the scores Paired Anlysis: Every protein family was independently cut out of the generated sequences and thus a MSA for only this protein created. These single MSAs were then paired for all protein pairs and used for inference and calculation of the scores.  teins are considered and the protein-protein interaction score is defined as the average of the 4 largest interaction scores on the residue level (as in the main paper). The true negative rate is the fraction of true negatives in the N pairs with the lowest interaction score, where N is the value indicated by the x-axis.