# Genepleio

Authors: Wenhai Chen; Dandan Chen; Ming Zhao; Yangyun Zou; Yanwu Zeng; Xun Gu
Journal: BioMed Research International
 (2015)

---

## Abstract

Though pleiotropy, which refers to the phenomenon of a gene affecting multiple traits, has long played a central role in genetics, development, and evolution, estimation of the number of pleiotropy components remains a hard mission to accomplish. In this paper, we report a newly developed software package,Genepleio, to estimate the effective gene pleiotropy from phylogenetic analysis of protein sequences. Since this estimate can be interpreted as the minimum pleiotropy of a gene, it is used to play a role of reference for many empirical pleiotropy measures. This work would facilitate our understanding of how gene pleiotropy affects the pattern of genotype-phenotype map and the consequence of organismal evolution.


---

## Body

## 1. Introduction

Understanding the role of gene pleiotropy in the map from genotypes to phenotypes has been one of the central topics for biologists, which refers to the phenomenon of a gene affecting multiple traits As a major measure for the functional importance of a gene, this concept has played a fundamental role in genetics, development, and evolution (see [1–3] for recent reviews and comments). However, the degree of gene pleiotropy remains largely unknown. Historically, proposed the concept of universal pleiotropy; that is, a single mutation can potentially affect all phenotypic traits. Though Fisher’s model has been widely used as a theoretical basis for exploring the evolutionary interplay between the genotype and phenotype, the notion of universal pleiotropy has not been well tested.Indeed, compared with the wide availability of genomics data, the whole-range phenotype recourses, or “phenomics,” are highly limited. Nevertheless, recent advances have been able to bring high throughput data to bear on the nature and extent of pleiotropy [4–6]. These experiments showed that the number of phenotypic traits that may be affected by a gene may be actually limited implying the role of modularity in shaping the degree of gene pleiotropy. Many controversial issues such as the problem of arbitrary number of correlated traits may directly affect the count of phenotypes that are predicted to be affected by a gene. On the other hand, a new approach has emerged in the past decade, to estimate the gene pleiotropy from genetics or sequence data, rather than from the affected phenotypes [7–13] (Chen et al., 2013). In particular, Gu [8] developed a practically feasible approach. It proposed that molecular evolution of a gene occurs in a multidimensional space corresponding to the same canonical number of molecular phenotypes. Random mutations of the gene could affect these molecular phenotypes constrained by the stabilizing selection. Moreover, Gu [8] developed a statistical method to estimate the “effective pleiotropy” (Ke) of a gene from the multiple sequence alignment of protein sequences. Most genes have Ke in the range between 1 and 20 [11], with the medium of Ke=6.5 of these estimates that is comparable to some empirical pleiotropy measures [1, 3]. Yet the relationship between these two approaches remains complex. As the degree of gene pleiotropy is a basic parameter for understanding the complexity of genotype-phenotype map, to facilitate this line of research we develop a software packageGenepleio (freely available at http://www.xungulab.com) to estimate the effective gene pleiotropy from the protein sequences.


## 2. Material and Methods

### 2.1. Sequences

Three groups of datasets include eight vertebrates, twelve fruit flies, and seven yeast species. Each dataset contains 300 random selected orthologous sets. Eight vertebrates are fugu, zebrafish, xenopus, chicken, dog, cow, mouse, and human. TwelveDrosophila species areD. melanogaster, D. pseudoobscura, D. sechellia,D.simulans,D. yakuba,D. erecta,D.ananassae,D.persimilis,D.willistoni,D. mojavensis, D.virilis, andD.grimshawi.Seven yeast species areCandida glabrata, Debaryomyces hansenii, Kluyveromyces lactis, Saccharomyces bayanus, Yarrowia lipolytica, Saccharomyces cerevisiae,andSchizosaccharomyces pombe.Vertebrate CDS and protein sequences were extracted from Ensmart, Drosophila CDS and protein sequences were extracted from FlyBase, and yeast CDS was extracted from ORNA Man’s dataset (corresponding protein sequences were translated by Bioperl).


### 2.2. Sequences Alignment

Multiple protein sequence alignment for each orthologous group was obtained by ClustalW at default settings.


### 2.3. Estimation ofdN/dS

The number of synonymous substitutions per synonymous site (dS) and the number of nonsynonymous substitutions per nonsynonymous site (dN) between human and mouse orthologs were calculated by CODEML of PAML package [14]. When calculating the variance of dN and dS, we changed “getSE = 1” in CODEML control file; otherwise, we used the default CODEML parameters. We used the estimates of dN/dS between human and mouse for vertebrates, those betweenD.melanogaster (dmel) and D. sechellia (dsec)forDrosophila, and those betweenS. bayanus and S. cerevisiae for yeasts, respectively.


## 2.1. Sequences

Three groups of datasets include eight vertebrates, twelve fruit flies, and seven yeast species. Each dataset contains 300 random selected orthologous sets. Eight vertebrates are fugu, zebrafish, xenopus, chicken, dog, cow, mouse, and human. TwelveDrosophila species areD. melanogaster, D. pseudoobscura, D. sechellia,D.simulans,D. yakuba,D. erecta,D.ananassae,D.persimilis,D.willistoni,D. mojavensis, D.virilis, andD.grimshawi.Seven yeast species areCandida glabrata, Debaryomyces hansenii, Kluyveromyces lactis, Saccharomyces bayanus, Yarrowia lipolytica, Saccharomyces cerevisiae,andSchizosaccharomyces pombe.Vertebrate CDS and protein sequences were extracted from Ensmart, Drosophila CDS and protein sequences were extracted from FlyBase, and yeast CDS was extracted from ORNA Man’s dataset (corresponding protein sequences were translated by Bioperl).


## 2.2. Sequences Alignment

Multiple protein sequence alignment for each orthologous group was obtained by ClustalW at default settings.


## 2.3. Estimation ofdN/dS

The number of synonymous substitutions per synonymous site (dS) and the number of nonsynonymous substitutions per nonsynonymous site (dN) between human and mouse orthologs were calculated by CODEML of PAML package [14]. When calculating the variance of dN and dS, we changed “getSE = 1” in CODEML control file; otherwise, we used the default CODEML parameters. We used the estimates of dN/dS between human and mouse for vertebrates, those betweenD.melanogaster (dmel) and D. sechellia (dsec)forDrosophila, and those betweenS. bayanus and S. cerevisiae for yeasts, respectively.


## 3. Results and Discussion

### 3.1. Estimation of Effective Gene Pleiotropy

Gu [8] analyzed the pleiotropy model of molecular evolution under the following assumptions. (i) K-dimensional molecular phenotypes (y) of the gene are under Gaussian-like stabilizing selection, indicating a single fitness optima for multiple functions. Any deviation from the optima is under the purifying selection. (ii) The fitness optima of y may shift randomly during the course of evolution, according to a multivariate normal distribution. It generates the process ofmicroadaptation that could be caused by the external (environmental) or internal (physiological) perturbations or the functional compensation for the previously fixed slightly deleterious mutation. (iii) And the distribution of mutational effects, p(y), follows a multivariate normal distribution.The estimation procedure implemented in the softwareGenepleio is summarized in Figure 1. We address several key issues concisely to help in understanding how the software works. One may see Gu [8], Su et al. [11], and Gu (2014) for technical details.Figure 1
A flow chart to outline the computational pipeline implemented in softwareGenepleio.


#### 3.1.1. Calculation ofH-Measure

Calculation ofH is the key step to estimate Ke. Biologically, H measures the strength of rate variation among sites: H=0 when var(λ)=0, and H=1 when var(λ)=∞. After the gene phylogeny is given or inferred by the NJ option, the software implemented the methods of Gu and Zhang [15] to infer the number of amino acid changes along the phylogeny at each site, after correcting the multiple hits. The next step is to calculate the mean (M) and variance (V) of the estimated number of changes over sites. Under the Poisson-based evolutionary model, H can be estimated by H=(V-M)/[V+M(M-1)].


#### 3.1.2. Estimation of Effective Gene Pleiotropy (Ke)

Genepleio has implemented the following procedure to estimate K. (i) Calculate the dN/dS ratio (the ratio of nonsynonymous to synonymous rates) used as an empirical measure for the mean sequence conservation. (ii) Calculate the g-function g=dN/dS/(1-H). (iii) And the effective gene pleiotropy (Ke) can be estimated by numerically solving the following equation:
(1)dN/dS1-H=2-Ke/21+ϕKe,
where ϕ(Ke)=0.0208Ke(Ke+2)/(1+0.289Ke).


#### 3.1.3. Estimation of Selection Intensities

There are two types of selection intensity measures. The first one is the (overall) selection intensity of the gene under study,S, for the overall strength of purifying selection imposed on the protein sequence; the negative sign indicates the negative (purifying) selection. The second one is the baseline selection intensity, B0, which is a scaled measure for the contribution of a single pleiotropy component. The relationship between B0 and S is S=-K×B0. When Ke is obtained, the software estimates the effective selection intensity Se for S and the effective baseline selection intensity (Be) for the baseline selection intensity B0.


#### 3.1.4. Calculation of Sampling Variance ofKe

The sampling variance ofKe can be approximately calculated by the delta method. Numerical analysis of (1) found that the following formula is accurate enough in practice:
(2)Var⁡Ke≈11.037Var⁡dN/dSdN/dS2+Var⁡H1-H2≈11.037Var⁡dNdN2+Var⁡dSdS2+Var⁡H1-H2.In (2), Var⁡(dN/dS) can be estimated by the delta method so that Var⁡dN/dS≈Var⁡dN/dN2+dN2*Var⁡dS/dS4, where Var⁡(dN), Var⁡(dS) are the variances of dN and dS, respectively. The sampling variance H is difficult to compute analytically.Genepleio has implemented a bootstrapping approach to calculate the sampling variance of H, Var⁡(H), whereas sampling variances of dN and dS, Var⁡(dN) and Var⁡(dS), depend on the users’ input; their default values are set to be zero. Using this method, we can bootstrap 100 times within 1~2 minutes.We use triosephosphate isomerase gene (TPI1, SWISSPROT P60174) for illustration. (i) Infer the phylogenetic tree from the multiple alignment of vertebrate homologous TPI1 protein sequences (human, mouse, dog, cow, chicken, xenopus, fugu, and zebrafish), which is consistent with the known vertebrate phylogeny. (ii) EstimatedN/dS=0.045 between the human and mouse genes by the likelihood method using PAML. (iii) Estimate H-index for the rate variation among sites; H=0.614 for TPI1 gene. (iv) And we estimated Ke=7.29 and the mean selection intensity S=-11.65. Then, the baseline selection intensity is given by B0=11.65/7.29≈1.60.The estimated effective gene pleiotropy varies among different treatments but the scale of variation is small. On the other hand, we found that when the number of changes at each site is estimated by the parsimony method without any correction, gene pleiotropy tends to be overestimated. At any rate, we conclude that these 5–10% estimation differences should not affect the general pattern about the degree of gene pleiotropy.


### 3.2. Biological Interpretation ofKe

The key question is how one can acquire the number of pleiotropy components of a gene without biologically knowing each component (Su et al., 2010) [12, 16]. Gu (2014) [17] addressed this issue, showing that the method of Gu [8] actually aims to estimate the rank (K) of genotype-phenotype map. The main result can be concisely represented by the following simple formula: K=min⁡(r,Pmin⁡), where Pmin⁡ is the minimum pleiotropy among all legitimate pleiotropy measures and r is the rank of mutational effects. In short, the meaning of “effective gene pleiotropy” (Ke) estimated by Gu-2007 method is as follows. (i) Ke is an estimate of K=min⁡(r,Pmin⁡), the rank of genotype-phenotype map. (ii) Ke is an estimate for the minimum pleiotropy Pmin⁡ only if Pmin⁡<r. (iii) Gu-2007 method attempted to estimate the pleiotropy of amino acid sites, a conserved proxy to the true minimum pleiotropy. (iv) With a sufficiently large phylogeny such that the rank of mutational effect at an amino acid site is r→19 (the number of amino acid types minus one), one can estimate Pmin⁡ in the range from 1 to 19 by this method. And (v) Ke is a conserved estimate of K because those pleiotropy components that have small effects on fitness would be effectively removed by the estimation procedure.


### 3.3. Software Overview

We have developedGenepleio, a GUI-based software package that estimates the effective gene pleiotropy from the phylogenetic sequence analysis of amino acids.Genepleio has three inputs. (i) Input the file of multisequence alignment (MSA) of protein sequences:Genepleio supports the alignment format of CLUSTALW. As required by the method, the multiple alignment file should contain at least four sequences with reasonably large sequence divergences between them. (ii) Input two values dN (nonsynonymous distance) and dS (synonymous distance). Several methods such as PAML can be used to obtain these estimates. We suggest choosing closely related sequences, say, dS<1, to avoid large sampling variance when calculating the dN/dS ratio. Note that the dN/dS ratio should be less than 1; otherwise, the gene may not be suitable to do this type of analysis. (iii) Input the tree file in the Phylip format. Alternatively, one may use the neighbor-joining (NJ) method implemented in the software to infer the gene phylogeny. As illustrated in Figure 2, the interface ofGenepleio includes three main tab pages: the first page is for the MSA input and dN/dS values, the second page is for the input or the inference of the phylogenetic tree, and the third page will output the results of estimation.Figure 2
Screen illustration of the softwareGenepleio.We have conducted a preliminary analysis and found that the 75% quantile of estimatedKe is typically within Ke±2, suggesting that Ke estimation as a measure of gene pleiotropy is statistically reliable. Besides, we notice that the contributions from Var⁡(dN) and Var⁡(dS) are nontrivial. In other words, the sampling variance of Ke would be severely underestimated if the user has no input for the sampling variances of dN and dS.There are some notices about usage ofGenepleio. First, the multiple alignment file should contain more than four sequences; second, the dN/dS value should be within (0,1); third, the sequences similarity >90% should be cautious because of the lack of statistical power; fourth, in order to shorten the time consumed, we do not give the mean of Ke value through bootstrapping inGenepleio. Nevertheless, according to much simulation, the mean value is close to the estimated Ke value, so we only give the estimated Ke value.


### 3.4. Computer Simulations

We have carried computer simulations to evaluate the software performance. We setK to vary from 1 to 20, with the fixed baseline selection intensity B0=0.5, 1.0, or 2.0, respectively. In particular, we consider three simulation models. (a) Independent-equal model: pleiotropy components are identical and independent of each other. (b) Independent-unequal model: pleiotropy components are independent of each other but have different strengths. (c) Random-matrix model: the strengths and correlations between pleiotropy components are randomly drawn from a specified random matrix model. Our main results are summarized in Table 1. In general, the estimation bias of Ke decreases with the increasing of the baseline selection intensity B0. For instance, in model (a), underestimation of Ke is considerable only when B0 is very small, say, 0.5 or less. The estimation bias becomes intermediate when B0>1 and becomes negligible when B0>3 (not shown). Moreover, the estimation bias of Ke may increase when the simulation model becomes more complex. Indeed, in model (c), Ke only describes the canonical number of pleiotropy components that could be much less than the number of pleiotropy components used in the simulation model.Table 1
Simulation results ofKe-estimation by Genepleio.
K
B
0
K
e (model (a))
K
e (model (b))
K
e (model (c))
2
0.5
0.98 ± 0.03
0.94 ± 0.02
0.95 ± 0.02
4
0.5
1.98 ± 0.03
1.96 ± 0.03
1.58  ± 0.03
8
0.5
4.05 ± 0.05
3.67 ± 0.05
2.49 ± 0.03
12
0.5
6.16 ± 0.06
4.76 ± 0.06
3.89 ± 0.05
16
0.5
8.29 ± 0.13
6.65 ± 0.10
4.36 ± 0.06
2
1.0
1.34 ± 0.04
1.31  ± 0.04
1.31 ± 0.04
4
1.0
2.71 ± 0.06
2.13 ±  0.05
2.13 ± 0.05
8
1.0
5.44 ± 0.13
4.99  ± 0.10
3.26 ± 0.07
12
1.0
8.17 ± 0.29
6.56  ± 0.16
5.22 ± 0.10
16
1.0
10.9 ± 0.84
9.10  ± 0.37
5.85 ± 0.17
2
2.0
1.64 ± 0.06
1.60 ±  0.03
1.61 ± 0.06
4
2.0
3.28 ± 0.12
3.25  ± 0.05
2.61 ± 0.08
8
2.0
6.50 ± 0.40
6.12  ± 0.07
4.04 ± 0.15
12
2.0
9.67 ± 0.65
8.27  ± 0.07
6.59 ± 0.32
16
2.0
13.02 ± 0.67
11.30 ± 0.22
7.45 ± 0.46


### 3.5. Case Studies

To validate the performance of the newly developed software for the estimation ofKe, we analyzed three datasets, each of which includes eight vertebrates, twelve fruit flies, or seven yeast species, respectively (Figure 3). Each dataset contains 300 random selected (one-to-one) orthologous sets. We calculated Ke, Se (selection intensity), and B0 (the baseline selection intensity). Our analysis shows that all Ke estimates are in a reasonable range with only a few outliners. Interestingly, the distribution of Ke estimates is similar across these distantly related species. The underlying reason to explain this similarity remains unclear. Ke is an important parameter for evolutionary analysis. Indeed, the square of coefficient correlation (r2) between Ke and Se is 0.64 in vertebrates, 0.94 in yeasts, and 0.55 in fruit flies, suggesting that gene pleiotropy may be an important evolutionary constraint in molecular evolution. In short, in this paper, we have reported a new software packageGenepleio and demonstrated the steps of gene pleiotropy (K) estimation. We also examined the extent to which the statistical properties of dN/dS and H affect the estimation efficiency of K and S. Comparison among three different groups of species validates the stability of K estimation procedure.Estimation ofKe from three datasets: eight vertebrates, twelve fruit flies, or seven yeast species, respectively. Each dataset contains 300 random selected (one-to-one) orthologous sets.
(a)
Vertebrates
(b)
Fruit flies
(c)
Yeasts


## 3.1. Estimation of Effective Gene Pleiotropy

Gu [8] analyzed the pleiotropy model of molecular evolution under the following assumptions. (i) K-dimensional molecular phenotypes (y) of the gene are under Gaussian-like stabilizing selection, indicating a single fitness optima for multiple functions. Any deviation from the optima is under the purifying selection. (ii) The fitness optima of y may shift randomly during the course of evolution, according to a multivariate normal distribution. It generates the process ofmicroadaptation that could be caused by the external (environmental) or internal (physiological) perturbations or the functional compensation for the previously fixed slightly deleterious mutation. (iii) And the distribution of mutational effects, p(y), follows a multivariate normal distribution.The estimation procedure implemented in the softwareGenepleio is summarized in Figure 1. We address several key issues concisely to help in understanding how the software works. One may see Gu [8], Su et al. [11], and Gu (2014) for technical details.Figure 1
A flow chart to outline the computational pipeline implemented in softwareGenepleio.


### 3.1.1. Calculation ofH-Measure

Calculation ofH is the key step to estimate Ke. Biologically, H measures the strength of rate variation among sites: H=0 when var(λ)=0, and H=1 when var(λ)=∞. After the gene phylogeny is given or inferred by the NJ option, the software implemented the methods of Gu and Zhang [15] to infer the number of amino acid changes along the phylogeny at each site, after correcting the multiple hits. The next step is to calculate the mean (M) and variance (V) of the estimated number of changes over sites. Under the Poisson-based evolutionary model, H can be estimated by H=(V-M)/[V+M(M-1)].


### 3.1.2. Estimation of Effective Gene Pleiotropy (Ke)

Genepleio has implemented the following procedure to estimate K. (i) Calculate the dN/dS ratio (the ratio of nonsynonymous to synonymous rates) used as an empirical measure for the mean sequence conservation. (ii) Calculate the g-function g=dN/dS/(1-H). (iii) And the effective gene pleiotropy (Ke) can be estimated by numerically solving the following equation:
(1)dN/dS1-H=2-Ke/21+ϕKe,
where ϕ(Ke)=0.0208Ke(Ke+2)/(1+0.289Ke).


### 3.1.3. Estimation of Selection Intensities

There are two types of selection intensity measures. The first one is the (overall) selection intensity of the gene under study,S, for the overall strength of purifying selection imposed on the protein sequence; the negative sign indicates the negative (purifying) selection. The second one is the baseline selection intensity, B0, which is a scaled measure for the contribution of a single pleiotropy component. The relationship between B0 and S is S=-K×B0. When Ke is obtained, the software estimates the effective selection intensity Se for S and the effective baseline selection intensity (Be) for the baseline selection intensity B0.


### 3.1.4. Calculation of Sampling Variance ofKe

The sampling variance ofKe can be approximately calculated by the delta method. Numerical analysis of (1) found that the following formula is accurate enough in practice:
(2)Var⁡Ke≈11.037Var⁡dN/dSdN/dS2+Var⁡H1-H2≈11.037Var⁡dNdN2+Var⁡dSdS2+Var⁡H1-H2.In (2), Var⁡(dN/dS) can be estimated by the delta method so that Var⁡dN/dS≈Var⁡dN/dN2+dN2*Var⁡dS/dS4, where Var⁡(dN), Var⁡(dS) are the variances of dN and dS, respectively. The sampling variance H is difficult to compute analytically.Genepleio has implemented a bootstrapping approach to calculate the sampling variance of H, Var⁡(H), whereas sampling variances of dN and dS, Var⁡(dN) and Var⁡(dS), depend on the users’ input; their default values are set to be zero. Using this method, we can bootstrap 100 times within 1~2 minutes.We use triosephosphate isomerase gene (TPI1, SWISSPROT P60174) for illustration. (i) Infer the phylogenetic tree from the multiple alignment of vertebrate homologous TPI1 protein sequences (human, mouse, dog, cow, chicken, xenopus, fugu, and zebrafish), which is consistent with the known vertebrate phylogeny. (ii) EstimatedN/dS=0.045 between the human and mouse genes by the likelihood method using PAML. (iii) Estimate H-index for the rate variation among sites; H=0.614 for TPI1 gene. (iv) And we estimated Ke=7.29 and the mean selection intensity S=-11.65. Then, the baseline selection intensity is given by B0=11.65/7.29≈1.60.The estimated effective gene pleiotropy varies among different treatments but the scale of variation is small. On the other hand, we found that when the number of changes at each site is estimated by the parsimony method without any correction, gene pleiotropy tends to be overestimated. At any rate, we conclude that these 5–10% estimation differences should not affect the general pattern about the degree of gene pleiotropy.


## 3.1.1. Calculation ofH-Measure

Calculation ofH is the key step to estimate Ke. Biologically, H measures the strength of rate variation among sites: H=0 when var(λ)=0, and H=1 when var(λ)=∞. After the gene phylogeny is given or inferred by the NJ option, the software implemented the methods of Gu and Zhang [15] to infer the number of amino acid changes along the phylogeny at each site, after correcting the multiple hits. The next step is to calculate the mean (M) and variance (V) of the estimated number of changes over sites. Under the Poisson-based evolutionary model, H can be estimated by H=(V-M)/[V+M(M-1)].


## 3.1.2. Estimation of Effective Gene Pleiotropy (Ke)

Genepleio has implemented the following procedure to estimate K. (i) Calculate the dN/dS ratio (the ratio of nonsynonymous to synonymous rates) used as an empirical measure for the mean sequence conservation. (ii) Calculate the g-function g=dN/dS/(1-H). (iii) And the effective gene pleiotropy (Ke) can be estimated by numerically solving the following equation:
(1)dN/dS1-H=2-Ke/21+ϕKe,
where ϕ(Ke)=0.0208Ke(Ke+2)/(1+0.289Ke).


## 3.1.3. Estimation of Selection Intensities

There are two types of selection intensity measures. The first one is the (overall) selection intensity of the gene under study,S, for the overall strength of purifying selection imposed on the protein sequence; the negative sign indicates the negative (purifying) selection. The second one is the baseline selection intensity, B0, which is a scaled measure for the contribution of a single pleiotropy component. The relationship between B0 and S is S=-K×B0. When Ke is obtained, the software estimates the effective selection intensity Se for S and the effective baseline selection intensity (Be) for the baseline selection intensity B0.


## 3.1.4. Calculation of Sampling Variance ofKe

The sampling variance ofKe can be approximately calculated by the delta method. Numerical analysis of (1) found that the following formula is accurate enough in practice:
(2)Var⁡Ke≈11.037Var⁡dN/dSdN/dS2+Var⁡H1-H2≈11.037Var⁡dNdN2+Var⁡dSdS2+Var⁡H1-H2.In (2), Var⁡(dN/dS) can be estimated by the delta method so that Var⁡dN/dS≈Var⁡dN/dN2+dN2*Var⁡dS/dS4, where Var⁡(dN), Var⁡(dS) are the variances of dN and dS, respectively. The sampling variance H is difficult to compute analytically.Genepleio has implemented a bootstrapping approach to calculate the sampling variance of H, Var⁡(H), whereas sampling variances of dN and dS, Var⁡(dN) and Var⁡(dS), depend on the users’ input; their default values are set to be zero. Using this method, we can bootstrap 100 times within 1~2 minutes.We use triosephosphate isomerase gene (TPI1, SWISSPROT P60174) for illustration. (i) Infer the phylogenetic tree from the multiple alignment of vertebrate homologous TPI1 protein sequences (human, mouse, dog, cow, chicken, xenopus, fugu, and zebrafish), which is consistent with the known vertebrate phylogeny. (ii) EstimatedN/dS=0.045 between the human and mouse genes by the likelihood method using PAML. (iii) Estimate H-index for the rate variation among sites; H=0.614 for TPI1 gene. (iv) And we estimated Ke=7.29 and the mean selection intensity S=-11.65. Then, the baseline selection intensity is given by B0=11.65/7.29≈1.60.The estimated effective gene pleiotropy varies among different treatments but the scale of variation is small. On the other hand, we found that when the number of changes at each site is estimated by the parsimony method without any correction, gene pleiotropy tends to be overestimated. At any rate, we conclude that these 5–10% estimation differences should not affect the general pattern about the degree of gene pleiotropy.


## 3.2. Biological Interpretation ofKe

The key question is how one can acquire the number of pleiotropy components of a gene without biologically knowing each component (Su et al., 2010) [12, 16]. Gu (2014) [17] addressed this issue, showing that the method of Gu [8] actually aims to estimate the rank (K) of genotype-phenotype map. The main result can be concisely represented by the following simple formula: K=min⁡(r,Pmin⁡), where Pmin⁡ is the minimum pleiotropy among all legitimate pleiotropy measures and r is the rank of mutational effects. In short, the meaning of “effective gene pleiotropy” (Ke) estimated by Gu-2007 method is as follows. (i) Ke is an estimate of K=min⁡(r,Pmin⁡), the rank of genotype-phenotype map. (ii) Ke is an estimate for the minimum pleiotropy Pmin⁡ only if Pmin⁡<r. (iii) Gu-2007 method attempted to estimate the pleiotropy of amino acid sites, a conserved proxy to the true minimum pleiotropy. (iv) With a sufficiently large phylogeny such that the rank of mutational effect at an amino acid site is r→19 (the number of amino acid types minus one), one can estimate Pmin⁡ in the range from 1 to 19 by this method. And (v) Ke is a conserved estimate of K because those pleiotropy components that have small effects on fitness would be effectively removed by the estimation procedure.


## 3.3. Software Overview

We have developedGenepleio, a GUI-based software package that estimates the effective gene pleiotropy from the phylogenetic sequence analysis of amino acids.Genepleio has three inputs. (i) Input the file of multisequence alignment (MSA) of protein sequences:Genepleio supports the alignment format of CLUSTALW. As required by the method, the multiple alignment file should contain at least four sequences with reasonably large sequence divergences between them. (ii) Input two values dN (nonsynonymous distance) and dS (synonymous distance). Several methods such as PAML can be used to obtain these estimates. We suggest choosing closely related sequences, say, dS<1, to avoid large sampling variance when calculating the dN/dS ratio. Note that the dN/dS ratio should be less than 1; otherwise, the gene may not be suitable to do this type of analysis. (iii) Input the tree file in the Phylip format. Alternatively, one may use the neighbor-joining (NJ) method implemented in the software to infer the gene phylogeny. As illustrated in Figure 2, the interface ofGenepleio includes three main tab pages: the first page is for the MSA input and dN/dS values, the second page is for the input or the inference of the phylogenetic tree, and the third page will output the results of estimation.Figure 2
Screen illustration of the softwareGenepleio.We have conducted a preliminary analysis and found that the 75% quantile of estimatedKe is typically within Ke±2, suggesting that Ke estimation as a measure of gene pleiotropy is statistically reliable. Besides, we notice that the contributions from Var⁡(dN) and Var⁡(dS) are nontrivial. In other words, the sampling variance of Ke would be severely underestimated if the user has no input for the sampling variances of dN and dS.There are some notices about usage ofGenepleio. First, the multiple alignment file should contain more than four sequences; second, the dN/dS value should be within (0,1); third, the sequences similarity >90% should be cautious because of the lack of statistical power; fourth, in order to shorten the time consumed, we do not give the mean of Ke value through bootstrapping inGenepleio. Nevertheless, according to much simulation, the mean value is close to the estimated Ke value, so we only give the estimated Ke value.


## 3.4. Computer Simulations

We have carried computer simulations to evaluate the software performance. We setK to vary from 1 to 20, with the fixed baseline selection intensity B0=0.5, 1.0, or 2.0, respectively. In particular, we consider three simulation models. (a) Independent-equal model: pleiotropy components are identical and independent of each other. (b) Independent-unequal model: pleiotropy components are independent of each other but have different strengths. (c) Random-matrix model: the strengths and correlations between pleiotropy components are randomly drawn from a specified random matrix model. Our main results are summarized in Table 1. In general, the estimation bias of Ke decreases with the increasing of the baseline selection intensity B0. For instance, in model (a), underestimation of Ke is considerable only when B0 is very small, say, 0.5 or less. The estimation bias becomes intermediate when B0>1 and becomes negligible when B0>3 (not shown). Moreover, the estimation bias of Ke may increase when the simulation model becomes more complex. Indeed, in model (c), Ke only describes the canonical number of pleiotropy components that could be much less than the number of pleiotropy components used in the simulation model.Table 1
Simulation results ofKe-estimation by Genepleio.
K
B
0
K
e (model (a))
K
e (model (b))
K
e (model (c))
2
0.5
0.98 ± 0.03
0.94 ± 0.02
0.95 ± 0.02
4
0.5
1.98 ± 0.03
1.96 ± 0.03
1.58  ± 0.03
8
0.5
4.05 ± 0.05
3.67 ± 0.05
2.49 ± 0.03
12
0.5
6.16 ± 0.06
4.76 ± 0.06
3.89 ± 0.05
16
0.5
8.29 ± 0.13
6.65 ± 0.10
4.36 ± 0.06
2
1.0
1.34 ± 0.04
1.31  ± 0.04
1.31 ± 0.04
4
1.0
2.71 ± 0.06
2.13 ±  0.05
2.13 ± 0.05
8
1.0
5.44 ± 0.13
4.99  ± 0.10
3.26 ± 0.07
12
1.0
8.17 ± 0.29
6.56  ± 0.16
5.22 ± 0.10
16
1.0
10.9 ± 0.84
9.10  ± 0.37
5.85 ± 0.17
2
2.0
1.64 ± 0.06
1.60 ±  0.03
1.61 ± 0.06
4
2.0
3.28 ± 0.12
3.25  ± 0.05
2.61 ± 0.08
8
2.0
6.50 ± 0.40
6.12  ± 0.07
4.04 ± 0.15
12
2.0
9.67 ± 0.65
8.27  ± 0.07
6.59 ± 0.32
16
2.0
13.02 ± 0.67
11.30 ± 0.22
7.45 ± 0.46


## 3.5. Case Studies

To validate the performance of the newly developed software for the estimation ofKe, we analyzed three datasets, each of which includes eight vertebrates, twelve fruit flies, or seven yeast species, respectively (Figure 3). Each dataset contains 300 random selected (one-to-one) orthologous sets. We calculated Ke, Se (selection intensity), and B0 (the baseline selection intensity). Our analysis shows that all Ke estimates are in a reasonable range with only a few outliners. Interestingly, the distribution of Ke estimates is similar across these distantly related species. The underlying reason to explain this similarity remains unclear. Ke is an important parameter for evolutionary analysis. Indeed, the square of coefficient correlation (r2) between Ke and Se is 0.64 in vertebrates, 0.94 in yeasts, and 0.55 in fruit flies, suggesting that gene pleiotropy may be an important evolutionary constraint in molecular evolution. In short, in this paper, we have reported a new software packageGenepleio and demonstrated the steps of gene pleiotropy (K) estimation. We also examined the extent to which the statistical properties of dN/dS and H affect the estimation efficiency of K and S. Comparison among three different groups of species validates the stability of K estimation procedure.Estimation ofKe from three datasets: eight vertebrates, twelve fruit flies, or seven yeast species, respectively. Each dataset contains 300 random selected (one-to-one) orthologous sets.
(a)
Vertebrates
(b)
Fruit flies
(c)
Yeasts



---

*Source: 269150-2015-01-05.xml*