Diversification of TF Paralogs Figure 3. Measuring tissue-specific expression divergence. Theoretical plot of the function dE = (1/EX+1/EY) where EX and EY are expression levels of the two paralogs The expression values for each TF gene in a paralogous pair are indicated on the x and y-axes, and the value of the function dE is depicted by the color of the point (darker color indicates high dE). (a) Using gene expression in Cerebellum (b) Using gene expression in Heart. The function dE roughly approximates the trend for the DNA binding motif similarity in Figure 2. We used a log scale for expression to accommodate extreme expression values. doi:10.1371/journal.pone.0002345.g003 of paralog genes. We have performed a number of additional analyses to ensure the robustness of our general conclusions. These include (1) investigating the effect of excluding paralogous pairs with low expressions, (2) using alternative measures of expression divergence, and (3) using several way of aggregating multiple probe data and multiple PWMs. All these analyses yielded consistent results and the details are provided in Data file S1. Tissue-specific expression diversification of paralogous transcription factor genes We next sorted the 95 pairs of paralogous TF genes in decreasing order of their DNA motif similarity §B. Figure 5a exhibits a ‘heat plot’ of these 95 pairs of paralogous TFs; the columns indicate the 23 non-redundant tissues and the rows comprise the TF pairs. The expression value for each TF was normalized in a tissue-specific fashion by subtracting the median expression and dividing the result by the Median Absolute Deviation (MAD) statistic. The median absolute deviation for a data sample [25] is defined as MAD = median(|Yi – median(Yi)|), where |Yi| is the absolute value of Yi. The MAD statistic is preferred over the traditional z-score normalization (mean of 0 and standard deviation of 1) if the data is not normally distributed, because in the latter case the mean and standard deviation are severely affected by outliers[26]. The color of each cell represents the normalized expression level in each tissue for each TF in the pair. Black cells indicate that the level of expression for this tissue was below our threshold for noise and hence, effectively zero. The general trend observed in the section above is underscored by Figure 5 in that for the paralogous pairs with similar DNA binding site motifs (high §B towards the top of the list), there are very few cases where both paralogs are expressed at high levels, while this pattern is a relatively common occurrence among paralogs with low DNA binding similarity. At a glance, there is more homogeneity of high expression as the binding similarity of the TF pair decreases, reinforcing that the amount of expression divergence increases as the binding similarity increases. Figure 5b quantifies this trend by showing a plot of the number of tissues wherein both TFs in a pair are highly expressed as a function of the binding site motif similarity, using a cubic smoothing spline with 5 degrees of freedom. A cubic spline with 5 degrees of freedom essentially fits 5 piece-wise cubic polynomials over the data range to provide a smooth interpolation of the data. The resulting Kendall tau correlation between the numbers of tissues for which both TFs in a pair are highly expressed versus the binding similarity score is 20.21 with p-value 0.008. This plot indicates a clear inverse trend between the number of tissues for which both TFs in a pair are highly expressed and the binding similarity score. In Figure 5a, the expression values are discretized into three levels. In Figure 2 in supplementary Data File S1 the expression values are discretized into five levels, thus providing a finer grained depiction of the analysis. Age of paralogs Because of the stringency of our criteria for paralogy, the TF pairs are likely to be recent paralogs. Therefore, we next investigated whether or not a greater DNA binding similarity is a simple reflection of shorter divergence time since duplication. In other words, it is expected that paralogs with similar DNA binding site motifs (high §B) are likely to be a result of recent duplication events and indeed if this were the case, the conservation in DNA binding site motif (and presumably in DNA binding domain of the gene) would merely be a reflection of divergence time. To test this hypothesis, for the entire set of TF paralogs we compared their BLAST-based percent identity score over the entire coding region (a high score is likely to correspond to recent duplication) with their §B values. The Kendal tau correlation was 0.063 with p-value = 0.36. We have repeated this analysis by substituting the BLAST-based percent identity with the synonymous substitution rate, Ks, using PAML [27] (see Methods), which yielded a consistent result with Tau = 20.067 and the p-value = 0.35. Thus, the observed relation- ships between hE and §B are independent of the age of the paralogs. Figure 3 in supplementary Data File S1 demonstrates this result. PLoS ONE | www.plosone.org 4 June 2008 | Volume 3 | Issue 6 | e2345