Figure 4. Schematic outlining the method for testing our hypothesis — correlation between expression divergence and DNA binding site motif similarity. For a set of human TF paralogs, and for a specific tissue, we estimate for each TF paralog, (i) tissue- specific expression (dE) and DNA binding site motif similarity (§B) (see text). Given these two values for each paralogous pair, we estimate their correlation using Kendall’s tau measure and estimate the significance of correlation based on expression randomization. doi:10.1371/journal.pone.0002345.g004 Shared functions between paralogs Next, we quantified the shared functions between paralogous TF pairs in order to determine if there is a correlation between DNA binding site motif similarity and biological function for a TF pair. For gene X, let F(X) represent the set of GO terms associated with X. We only consider terms at the functional hierarchy level of 3 or higher in order to exclude ubiquitous terms. For paralogs X and Y, we define functional coherence between the paralogs using the Jaccard’s coefficient as FC(X,Y) = |(F(X) > F(Y)| / |(F(X) < F(Y)|, which represents the number of shared GO terms normalized by the total number GO terms annotated for the two genes X and Y. For the list of paralogs, sorted by decreasing order of §B, Figure 6a shows the heat plot and Figure 6b shows the trend plot of FC. We found that FC values were higher for TF pairs with higher binding similarity and gradually decreased as the binding similarity decreased. The Kendall tau correlation coefficient is 0.165 with p-value = 0.02. This indicates that paralogs with diverged DNA binding site motifs, and presumably different target genes, are Diversification of TF Paralogs more likely to serve distinct functional roles. Our results were consistent when we substituted the Jaccard’s coefficient with a hypergeometric distribution based p-value to capture the overlap (Tau = 20.140, p-value = 0.04). This result further reinforces that paralogs with diverged DNA binding site motifs are more likely to serve distinct functional roles. Shared tissues between paralogs We quantified the shared tissues where paralogs are highly expressed and whether there is a correlation between the number of shared tissues and DNA binding site motif similarity. For paralogs X and Y, we defined tissue coherence between the paralogs as TC(X,Y) = |(T(X) > T(Y)|/(T(X) < T(Y)|. T(X) represents the set of tissues in which X is expressed above median expression value for that tissue. In the list of paralogs sorted by decreasing order of §B (Figure 7a), we found as expected from the above analysis, that TC values significantly decreased with increasing binding similarity. The Kendall tau correlation is 20.156 with p-value 0.03 (Figure 7b). In addition, the total number of tissues for which either of the TFs in a pair are significantly expressed also decreases as binding similarity decreases (Kendall tau correlation = 20.212, p-value = 3.7e-3). However, when we substituted the Jaccard’s coefficient with the hypergeometric p-value, the correlation between the tissue overall and motif similarity is no longer significant. This is likely because the overlaps themselves are not significant in most cases (only in 12 of the 95 cases the overlap p-value#0.05). Thus, while the tissue-overlap may not be significant, the degree of overlap is significantly correlated with the motif similarity. In other words, paralogs with similar DNA binding site motifs are expressed significantly less frequently within the same tissue at high levels, than paralogs with dissimilar DNA binding site motifs. Exclusive expression of one of the family members in a tissue The above results suggest that among the paralogous TF pairs with similar DNA binding site motifs, at most one TF is expressed at a high level in any given tissue. Hence, we next extend our analysis to TF families. We computed families using two methods, the first is a complete-linkage agglomerative clustering of the 98 pair-wise paralogy relationships, which resulted in 39 paralog families with 28 families of size 2, 7 of size 3 and 4 of size 4; the rest are singletons which could not be placed in a family. We have repeated the family analysis for BLAST alignable coverage of 50%, 60%, 70% and 80%. The results are consistent for all four choices of alignable coverage, including criteria identical to that of Makova and Li, but we only report the statistics for the 70% coverage threshold. The second method is based on KOG data and in terms of families, we obtained 24 families, 9 of size 2, 6 of size 3, and 9 of size greater than or equal to 4. The family-wise analysis was done for (i) the 11 families consisting of greater than two TFs constructed using clustering methods and (ii) the 15 families computed using KOG data. For a TF family we expect at most one member of the family among the members with similar DNA binding site motif to have high expression in any tissue. In order to test this hypothesis, for each family with k TFs (TF1, TF2, .., TFk) and each tissue T, we computed the ratio R = (Emax-Esim)/ (Emax-Emin), where Emax and Emin represent the maximum and the minimum expression levels, respectively, among the family members in tissue T. Thus, the R value is normalized for family-specific and tissue-specific gene expression. If TFi has the maximum expression Emax, then Esim is the expression level of the family member whose DNA binding site motif is most similar to that of TFi. This ratio effectively indicates how much the highest expressed gene in a family has diverged in expression relative to PLoS ONE | www.plosone.org 5 June 2008 | Volume 3 | Issue 6 | e2345