Separation of phylogenetic and functional associations in biological sequences by using the parametric bootstrap

Abstract

Quantitative analyses of biological sequences generally proceed under the assumption that individual DNA or protein sequence elements vary independently. However, this assumption is not biologically realistic because sequence elements often vary in a concerted manner resulting from common ancestry and structural or functional constraints. We calculated intersite associations among aligned protein sequences by using mutual information. To discriminate associations resulting from common ancestry from those resulting from structural or functional constraints, we used a parametric bootstrap algorithm to construct replicate data sets. These data are expected to have intersite associations resulting solely from phylogeny. By comparing the distribution of our association statistic for the replicate data against that calculated for empirical data, we were able to assign a probability that two sites covaried resulting from structural or functional constraint rather than phylogeny. We tested our method by using an alignment of 237 basic helix-loop-helix (bHLH) protein domains. Comparison of our results against a solved three-dimensional structure confirmed the identification of several sites important to function and structure of the bHLH domain. This analytical procedure has broad utility as a first step in the identification of sites that are important to biological macromolecular structure and function when a solved structure is unavailable.

Figure 1

Inverse cumulative frequency distribution of MI values for the alignment of 237 bHLH protein sequences and 1,000 parametric bootstrap replicates using either the JTT substitution matrix or the rind substitution matrix. MI values were calculated by using Eq. 1 with n = 20 so that the maximum possible value is unity. Line A is the P < 0.01 threshold for the JTT replicates at MI = 0.188. Line B is the P < 0.001 threshold for the JTT replicates at MI = 0.250. Line A′ is the P < 0.01 threshold for the rind replicates at MI = 0.359. Line B′ is the P < 0.001 threshold for the JTT replicates at MI = 0.408. These were the MI values that were >99% (for P < 0.01) or 99.9% (for P < 0.001) of the MI values calculated in the parametric bootstrap replicates. Because MI is a pairwise measure, x(x − 1)/2 values were calculated in each replicate, where x is the number of nongapped sites in the alignment. For the alignment of 237 bHLH sequences, there were 32 sites without gaps, resulting in 496 MI values per replicate.