Modeling Expression Evolution after Gene Duplication as an OU Process

To design a model of expression evolution after gene duplication, we consider a pair of related species, Species 1 and Species 2, whose lineages diverged from that of a common ancestor at time $$ T_{\text{PCA}}$$ (fig. 2A–C). Suppose that the common ancestor had a single-copy gene that underwent duplication, giving rise to a pair of duplicate genes at time $$ T_{\text{PC}}$$ in the lineage of Species 1 after its divergence from the lineage of Species 2. Of the pair of duplicate genes in Species 1, we designate the copy corresponding to the original single-copy gene in the ancestor and in Species 2 as the parent, and the new copy that is absent in both the ancestor and Species 2 as the child. Further, suppose that optimal expression states for the parent, child, and ancestral genes are given by $$ {\theta }_{\mathrm{P}},\hspace{0.17em}{\theta }_{\mathrm{C}}$$, and $$ {\theta }_{\mathrm{A}}$$, respectively. Likewise, the optimal expression state for the single-copy gene in the ancestor prior to the divergence of Species 1 and Species 2 is given as $$ {\theta }_{\text{PCA}}$$, and for the single-copy gene in the lineage of Species 1 before duplication occurred as $$ {\theta }_{\text{PC}}$$. Additionally, assume that $$ {\theta }_{\text{PC}}={\theta }_{\text{PCA}}={\theta }_{\mathrm{A}}$$. We then model expression along the tree relating the parent, child, and ancestral genes as changing randomly through phenotypic drift with strength $$ {\sigma }^2$$, and toward the optimal expression state through selection with strength α, according to an OU process.

Fig. 2

Modeling expression evolution after gene duplication as an OU process. (A) Relationships between two species (black phylogeny) and their genes (green phylogeny). After the two species diverged, a blue gene in Species 1 (parent) underwent a duplication event to create a yellow copy (child). (B) Relationships among the parent gene copy (P) in Species 1, child gene copy (C) in Species 1, and ancestral single-copy gene (A) in Species 2. The duplication event occurred at time $$ T_{\text{PC}}$$, and both copies split from the ancestral gene at time $$ T_{\text{PCA}}$$. Optimal expression states for the parent, child, and ancestral genes are given by $$ {\theta }_{\mathrm{P}},\hspace{0.17em}{\theta }_{\mathrm{C}}$$, and $$ {\theta }_{\mathrm{A}}$$, respectively. The internal branch and the branch above the root have optimal expression states $$ {\theta }_{\text{PC}}$$ and $$ {\theta }_{\text{PCA}}$$, respectively. (C) Cartoon depicting expression profile changes (red lines) along the gene tree. Expression profiles change randomly through phenotypic drift with strength $$ {\sigma }^2$$, and toward the optimal expression state through selection with strength α. (D) Illustration of how we simulate multi-tissue expression vectors for parent, child, and ancestral genes. (E) Schematic of our feed-forward neural network architecture, which takes in p input units with values $$ x_1,x_2,\dots ,x_p$$, has K output units with values $$ y_1,y_2,\dots ,y_K$$, and has L hidden layers, where the number of units in layer $$ \ell $$ is $$ p[\ell ]$$ and the value of unit k in layer $$ \ell $$ is the activation $$ a_k^{[\ell ]}$$.

In each tissue, gene expression $$ \mathbf{e}=(e_{\mathrm{P}},e_{\mathrm{C}},e_{\mathrm{A}})\in {ℝ}^3$$ is therefore distributed as a multivariate normal (MVN) distribution with mean

Here, we assume that expression is independent across tissues. However, this approach can also be extended to account for the intertissue expression covariance structure using established approaches (Revell and Harmon 2008; Revell and Collar 2009; Clavel et al. 2015).

Neural Network Architecture for the CLOUD Classifier and Predictor

We denote the set of all genes with two copies in one species and one copy in the other as duplicate genes $$ \mathcal{D}$$, and the set of all genes with one copy in both species as single-copy genes $$ \mathcal{G}$$. Let

We transform and compare the expression vector $$ {\mathbf{e}}^{(d)}$$ of each duplicate gene $$ d\in \mathcal{D}$$ to the expression vector $$ {\mathbf{s}}^{(g)}$$ of each single-copy gene $$ g\in \mathcal{G}$$ to obtain the feature vector

Table 1

Set of $$ p=4m+84$$ derived features used as input to CLOUD.

Feature number	Feature value
1	$$ t_{\text{PC}}=T_{\text{PC}}/T_{\text{PCA}}$$
2 to $$ 3m+1$$	$$ {\mathbf{e}}^{(d)}=(e_{\mathrm{P}1}^{(d)},e_{\mathrm{C}1}^{(d)},e_{\mathrm{A}1}^{(d)},\dots ,e_{\mathrm{P}m}^{(d)},e_{\mathrm{C}m}^{(d)},e_{\mathrm{A}m}^{(d)})$$
3 m + 2 to $$ 4m+1$$	$$ e_{\mathrm{P}1}^{(d)}+e_{\mathrm{C}1}^{(d)},\dots ,e_{\mathrm{P}m}^{(d)}+e_{\mathrm{C}m}^{(d)}$$
4 m + 2	$$ \text{dist}(\mathrm{P},\mathrm{C})=\sqrt{{(e_{\mathrm{P}1}^{(d)}-e_{\mathrm{C}1}^{(d)})}^2+\cdots +{(e_{\mathrm{P}m}^{(d)}-e_{\mathrm{C}m}^{(d)})}^2}$$
4 m + 3	$$ \text{dist}(\mathrm{P},\mathrm{A})=\sqrt{{(e_{\mathrm{P}1}^{(d)}-e_{\mathrm{A}1}^{(d)})}^2+\cdots +{(e_{\mathrm{P}m}^{(d)}-e_{\mathrm{A}m}^{(d)})}^2}$$
4 m + 4	$$ \text{dist}(\mathrm{C},\mathrm{A})=\sqrt{{(e_{\mathrm{C}1}^{(d)}-e_{\mathrm{A}1}^{(d)})}^2+\cdots +{(e_{\mathrm{C}m}^{(d)}-e_{\mathrm{A}m}^{(d)})}^2}$$
4 m + 5	$$ \text{dist}(\text{PC},\mathrm{A})=\sqrt{{(e_{\mathrm{P}1}^{(d)}+e_{\mathrm{C}1}^{(d)}-e_{\mathrm{A}1}^{(d)})}^2+\cdots +{(e_{\mathrm{P}m}^{(d)}+e_{\mathrm{C}m}^{(d)}-e_{\mathrm{A}m}^{(d)})}^2}$$
4 m + 6	$$ \text{branch}(\mathrm{P})=[\text{dist}(\mathrm{P},\mathrm{C})+\text{dist}(\mathrm{P},\mathrm{A})-\text{dist}(\mathrm{C},\mathrm{A})]/2$$
4 m + 7	$$ \text{branch}(\mathrm{C})=[\text{dist}(\mathrm{P},\mathrm{C})+\text{dist}(\mathrm{C},\mathrm{A})-\text{dist}(\mathrm{P},\mathrm{A})]/2$$
4 m + 8	$$ \text{branch}(\mathrm{A})=[\text{dist}(\mathrm{P},\mathrm{A})+\text{dist}(\mathrm{C},\mathrm{A})-\text{dist}(\mathrm{P},\mathrm{C})]/2$$
4 m + 9	rank of $$ \text{dist}(\mathrm{P},\mathrm{C})$$ among $$ \text{dist}(\mathcal{G})$$
4 m + 10	rank of $$ \text{dist}(\mathrm{P},\mathrm{A})$$ among $$ \text{dist}(\mathcal{G})$$
4 m + 11	rank of $$ \text{dist}(\mathrm{C},\mathrm{A})$$ among $$ \text{dist}(\mathcal{G})$$
4 m + 12	rank of $$ \text{dist}(\text{PC},\mathrm{A})$$ among $$ \text{dist}(\mathcal{G}$$
4 m + 13 to 4 m + 20	kth moment of $$ [\text{dist}(\mathrm{P},\mathrm{C})-\text{dist}(\mathcal{G})]/\max \{\text{dist}(\mathcal{G})\}$$ for $$ k=1,2,\dots ,8$$
4m + 21 to 4 m + 28	kth moment of $$ [\text{dist}(\mathrm{P},\mathrm{A})-\text{dist}(\mathcal{G})]/\max \{\text{dist}(\mathcal{G})\}$$ for $$ k=1,2,\dots ,8$$
4m + 29 to 4 m + 36	kth moment of $$ [\text{dist}(\mathrm{C},\mathrm{A})-\text{dist}(\mathcal{G})]/\max \{\text{dist}(\mathcal{G})\}$$ for $$ k=1,2,\dots ,8$$
4m + 37 to 4 m + 44	kth moment of $$ [\text{dist}(\text{PC},\mathrm{A})-\text{dist}(\mathcal{G})]/\max \{\text{dist}(\mathcal{G})\}$$ for $$ k=1,2,\dots ,8$$
4m + 45	$$ \text{cor}(\mathrm{P},\mathrm{C})=\text{PearsonCorrelation}(e_{\mathrm{P}1}^{(d)},\dots ,e_{\mathrm{P}m}^{(d)};e_{\mathrm{C}1}^{(d)},\dots ,e_{\mathrm{C}m}^{(d)})$$
4 m + 46	$$ \text{cor}(\mathrm{P},\mathrm{A})=\text{PearsonCorrelation}(e_{\mathrm{P}1}^{(d)},\dots ,e_{\mathrm{P}m}^{(d)};e_{\mathrm{A}1}^{(d)},\dots ,e_{\mathrm{A}m}^{(d)})$$
4 m + 47	$$ \text{cor}(\mathrm{C},\mathrm{A})=\text{PearsonCorrelation}(e_{\mathrm{C}1}^{(d)},\dots ,e_{\mathrm{C}m}^{(d)};e_{\mathrm{A}1}^{(d)},\dots ,e_{\mathrm{A}m}^{(d)})$$
4 m + 48	$$ \text{cor}(\text{PC},\mathrm{A})=\text{PearsonCorrelation}(e_{\mathrm{P}1}^{(d)}+e_{\mathrm{C}1}^{(d)},\dots ,e_{\mathrm{P}m}^{(d)}+e_{\mathrm{C}m}^{(d)};e_{\mathrm{A}1}^{(d)},\dots ,e_{\mathrm{A}m}^{(d)})$$
4 m + 49	rank of $$ \text{cor}(\mathrm{P},\mathrm{C})$$ among $$ \text{cor}(\mathcal{G})$$
4 m + 50	rank of $$ \text{cor}(\mathrm{P},\mathrm{A})$$ among $$ \text{cor}(\mathcal{G})$$
4 m + 51	rank of $$ \text{cor}(\mathrm{C},\mathrm{A})$$ among $$ \text{cor}(\mathcal{G})$$
4 m + 52	rank of $$ \text{cor}(\text{PC},\mathrm{A})$$ among $$ \text{cor}(\mathcal{G})$$
4 m + 53 to 4 m + 60	kth moment of $$ \text{cor}(\mathrm{P},\mathrm{C})-\text{cor}(\mathcal{G})$$ for $$ k=1,2,\dots ,8$$
4m + 61 to 4 m + 68	kth moment of $$ \text{cor}(\mathrm{P},\mathrm{A})-\text{cor}(\mathcal{G})$$ for $$ k=1,2,\dots ,8$$
4m + 69 to 4 m + 76	kth moment of $$ \text{cor}(\mathrm{C},\mathrm{A})-\text{cor}(\mathcal{G})$$ for $$ k=1,2,\dots ,8$$
4m + 77 to 4 m + 84	kth moment of $$ \text{cor}(\text{PC},\mathrm{A})-\text{cor}(\mathcal{G})$$ for $$ k=1,2,\dots ,8$$

Given the input feature vector $$ {\mathbf{x}}^{(d)}$$, we seek to predict the output vector

When performing classification of duplicate gene retention mechanisms, $$ {\mathbf{y}}^{(d)}$$ is the vector of K = 5 class probabilities, corresponding to class labels “Conserved” for conservation, “Neofunctionalized parent” for neofunctionalization in which the parent copy acquires a new function, “Neofunctionalized child” for neofunctionalization in which the child copy acquires a new function, “Subfunctionalized” for subfunctionalization, and “Specialized” for specialization. In contrast, when predicting evolutionary parameters of duplicate genes, $$ {\mathbf{y}}^{(d)}$$ is the vector of $$ K=5m$$ parameter predictions in each of the m tissues, where in each tissue we obtain parameter estimates for $$ {\theta }_{\mathrm{P}},\hspace{0.17em}{\theta }_{\mathrm{C}},\hspace{0.17em}{\theta }_{\mathrm{A}},\hspace{0.17em}{\sigma }^2$$, and α.

We consider a dense feed-forward neural network with $$ L\in \{0,1,2,3\}$$ hidden layers. The first hidden layer has $$ p[1]=256$$ hidden units, and hidden layer $$ \ell \in \{1,2,\dots ,L\}$$ has $$ p[\ell ]=256/2^{\ell -1}$$ hidden units, such that each hidden layer has half the number of hidden units as the previous hidden layer. For the purposes of condensing notation, we also consider the input layer as hidden layer zero, such that $$ p[0]=p=4m+84$$ is the number of input features, and we consider the output layer as hidden layer L + 1, such that $$ p[L+1]=K$$.

We define the values at unit $$ k\in \{1,2,\dots ,p[\ell ]\}$$ of hidden layer $$ \ell \in \{0,1,2,\dots ,L\}$$ for duplicate gene $$ d\in \mathcal{D}$$ by its activation $$ a_k^{(d)[\ell ]}$$. Because hidden layer zero is the input layer and hidden layer L + 1 is the output layer, then

For hidden layer $$ \ell \in \{1,2,\dots ,L\}$$, we define the activation for unit k as a linear combination of the activations from the previous hidden layers, followed by a non-linear transformation (Goodfellow et al. 2016). Here we choose the rectified linear unit (ReLU; Goodfellow et al. 2016) function defined as $$ \text{ReLU}(x)=\max (0,x)$$, such that the activation for unit k in hidden layer $$ \ell $$ of duplicate gene d is

For the prediction problem, we instead use the linear activation function (Goodfellow et al. 2016), such that the output for parameter prediction $$ k\in \{1,2,\dots ,5m\}$$ of duplicate gene d is

This neural network was implemented in R (R Core Team 2013), using Keras (Allaire and Chollet 2017) with a TensorFlow backend (Abadi et al. 2015). A schematic of the neural network architecture is provided in figure 2E. Note that when L = 0, the neural network simplifies to a multinomial regression model (Hastie et al. 2009) for the classification problem, and to a linear regression model (Hastie et al. 2009) for the prediction problem.

Classification Power and Accuracy of CLOUD Relative to CDROM

To evaluate the classification power and accuracy of our multi-layer neural network classifier CLOUD, we trained and tested it on independent data sets simulated under each class of duplicate gene retention mechanisms (see Materials and Methods section). We assumed two hidden layers when training and testing CLOUD, as this resulted in the best cross-validation performance (see Materials and Methods section). Our training set consisted of 50,000 observations, of which 10,000 were simulated under each class. We trained CLOUD on these data, and explored evolutionary parameters drawn on a logarithmic scale across many orders of magnitude. Specifically, we independently drew the five parameters $$ {\theta }_{\mathrm{P}},{\theta }_{\mathrm{C}},{\theta }_{\mathrm{A}}\in [{10}^{-4},{10}^4],\hspace{0.17em}\alpha \in [1,{10}^3]$$, and $$ {\sigma }^2\in [{10}^{-2},{10}^3]$$ for each of the six tissues, for a total of 30 random parameters per simulated replicate. We then tested CLOUD on a separate set of 5,000 observations, of which 1,000 were simulated under each class, with evolutionary parameters drawn from the same broad space as that of the training set. For comparison, we also applied the existing classifier CDROM (Assis and Bachtrog 2013; Perry and Assis 2016) to the same simulated test set (see Materials and Methods section).

Analysis of the resulting classifications reveals that CLOUD generally has substantially higher power (fig. 3A) and accuracy (fig. 3B) than CDROM. Specifically, across the wide range of test parameter values explored, CLOUD achieved an accuracy of 80.18%, whereas CDROM only reached 68.76% accuracy (fig. 3B). This represents a boost in over 16% classification accuracy of CLOUD relative to CDROM. In addition to increased overall accuracy, CLOUD yields similar accuracies across classes, illustrated by narrow ranges of correct classification rates (between 77.1 and 84.2%; diagonal cells of fig. 3B) and mis-classification rates (between 1.3 and 9.2%; non-diagonal cells of fig. 3B). In contrast, CDROM demonstrates a much higher correct classification rate for the “Specialized” class (96.9%) than for other classes (between 45.7 and 67.9%; fig. 3B), and a higher mis-classification rate toward the “Specialized” class (between 28.7 and 30.6%; fig. 3B). Moreover, CDROM experiences additional issues when classifying true “Subfunctionalized” observations, with an 18.3% mis-classification rate toward the “Conserved” class (fig. 3B).

Fig. 3

Classification results for CDROM and CLOUD with L = 2 hidden layers applied to data simulated under parameters $$ \alpha \in [1,{10}^3]$$ and $$ {\sigma }^2\in [{10}^{-2},{10}^3]$$. (A) Receiver operating characteristic curves across the full range of false positive rates (left) and truncated at a false positive rate of 5% (right). Because CDROM is a decision tree classifier, its true positive and false positive rates are plotted as points. (B) Confusion matrices depicting the classification rates of each of the five duplicate gene retention classes for CDROM (left) and CLOUD (right).

In addition, CLOUD is much more conservative than CDROM for pairs of α and $$ {\sigma }^2$$ values that are difficult to classify (supplementary figs. S1–S4, Supplementary Material online). For example, both methods typically have higher power (supplementary figs. S3 and S4, Supplementary Material online) and accuracy (supplementary figs. S1 and S2, Supplementary Material online) when either selection is strong (large α) or random phenotypic drift is weak (small $$ {\sigma }^2$$). In contrast, when selection is weak (small α) and phenotypic drift is strong (large $$ {\sigma }^2$$), then classification is more difficult for both methods. However, in these cases, CLOUD tends to choose classes at similar rates, whereas CDROM is overconfident and chooses the “Specialized” class regardless of the true class (compare supplementary figs. S1 and S2, Supplementary Material online). Therefore, CLOUD not only demonstrates uniformly higher power and accuracy than CDROM across a wide array of evolutionary settings but is also unbiased unlike CDROM.

Parameter Prediction Accuracy of CLOUD

In addition to its vastly improved classification performance relative to CDROM, a unique attribute of CLOUD is its ability to learn parameters underlying the expression evolution of duplicate genes. Thus, we next assessed the accuracy of the CLOUD predictor by training and testing it on the same independent simulated data sets that we employed for training and testing the CLOUD classifier. In particular, we trained CLOUD (again assuming two hidden layers) to make predictions for each of the five parameters ($$ {\theta }_{\mathrm{P}},\hspace{0.17em}{\theta }_{\mathrm{C}},\hspace{0.17em}{\theta }_{\mathrm{A}}$$, α, and $$ {\sigma }^2$$) in six tissues (total of 30 parameters) from the training set, and then applied it to make predictions for these parameters from the test set (see Materials and Methods section).

To investigate prediction accuracy, we examined the distributions of mean parameter prediction errors across the six tissues (fig. 4). In general, all parameter estimates appear unbiased, with mean prediction errors centered on zero. Moreover, estimates of optimal expression states ($$ {\theta }_{\mathrm{P}},\hspace{0.17em}{\theta }_{\mathrm{C}}$$, and $$ {\theta }_{\mathrm{A}}$$) are more precise than those of selection strength (α), which are more precise than those of phenotypic drift ($$ {\sigma }^2$$). Further, parameter predictions for the “Specialized” class are less precise than those for other classes, likely due to the additional degrees of freedom in estimating parameters for this class. In particular, for the “Specialized” class, all optimal expression values are unconstrained (table 2), whereas at least two of the three optimal expression states are constrained to be identical in the “Conserved” and “Neofunctionalized” classes, and $$ {\theta }_{\mathrm{P}}$$ and $$ {\theta }_{\mathrm{C}}$$ are constrained to sum to $$ {\theta }_{\mathrm{A}}$$ for the “Subfunctionalized” class (table 2).

Fig. 4

Prediction results for application of CLOUD with L = 2 hidden layers to data simulated under parameters $$ \alpha \in [1,{10}^3]$$ and $$ {\sigma }^2\in [{10}^{-2},{10}^3]$$ for each of the five classes of duplicate gene retention mechanisms. Violin plots display distributions of mean parameter prediction errors across the m = 6 tissues for each simulated test set.

Table 2

Optimal expression states under OU processes used to simulate the five classes of duplicate gene retention mechanisms.

Class	Optimal expression state at a given tissue
Conserved	$$ {\theta }_{\mathrm{P}}={\theta }_{\mathrm{C}}={\theta }_{\mathrm{A}}=\theta $$
Neofunctionalized parent	$$ {\theta }_{\mathrm{C}}={\theta }_{\mathrm{A}}=\theta $$ and $$ {\theta }_{\mathrm{P}}\ne \theta $$
Neofunctionalized child	$$ {\theta }_{\mathrm{P}}={\theta }_{\mathrm{A}}=\theta $$ and $$ {\theta }_{\mathrm{C}}\ne \theta $$
Subfunctionalized	$$ {\theta }_{\mathrm{A}}=\theta ,\hspace{0.17em}{\theta }_{\mathrm{P}}\ne \theta ,\hspace{0.17em}{\theta }_{\mathrm{C}}\ne \theta $$, and $$ {\theta }_{\mathrm{P}}+{\theta }_{\mathrm{C}}=\theta $$
Specialized	$$ {\theta }_{\mathrm{A}}=\theta ,\hspace{0.17em}{\theta }_{\mathrm{P}}\ne \theta ,\hspace{0.17em}{\theta }_{\mathrm{C}}\ne \theta $$, and $$ {\theta }_{\mathrm{P}}+{\theta }_{\mathrm{C}}\ne \theta $$

As with classification, confidence in parameter predictions made by CLOUD also vary with α and $$ {\sigma }^2$$ (supplementary fig. S5, Supplementary Material online). Though precision in estimation tends to be highest when selection is strong (large α) or phenotypic drift is weak (small $$ {\sigma }^2$$), it decreases as selection becomes weaker (smaller α) or phenotypic drift becomes stronger (larger $$ {\sigma }^2$$). Further, as with our general results across a wide parameter space (fig. 4), estimates of optimal expression states ($$ {\theta }_{\mathrm{P}},\hspace{0.17em}{\theta }_{\mathrm{C}}$$, and $$ {\theta }_{\mathrm{A}}$$) appear to be unbiased even in narrow regions of the space, with mean prediction errors centered on zero. In contrast, estimates of α and $$ {\sigma }^2$$ are biased for some pairs of values. Specifically, estimates of α and $$ {\sigma }^2$$ are both upwardly biased for weak selection (small α) with weak phenotypic drift (small $$ {\sigma }^2$$), and downwardly biased for strong selection (large α) with strong phenotypic drift (large $$ {\sigma }^2$$).

CLOUD Behavior under Non-uniform Retention Mechanisms across Tissues

We showed that under ideal settings, CLOUD is a superior classifier to CDROM, and is also adept at predicting underlying evolutionary parameters. Thus, we next explored the performance of the trained CLOUD classifier and predictor on test data generated under scenarios that violated model assumptions of the training data. In particular, we considered test data in which the duplicate gene retention mechanism was non-uniform across the simulated tissues. Specifically, we evaluated scenarios in which $$ k\in \{1,2,\dots ,m-1\}$$ tissues shared one retention mechanism (denoted Tissue Mechanism A) and the remaining m − k tissues shared a different mechanism (denoted Tissue Mechanism B). As for the trained CLOUD classifier and predictor, we assumed m = 6 tissues, and explored all possible distinct scenarios in which k tissues shared one mechanism and m − k tissues shared a different mechanism. For each setting, we evaluated 1,000 independent replicate test data sets, with tissue model parameters drawn from the same wide distributions as in the training data set.

Comparisons among these diverse scenarios illustrates that the classification performance of CLOUD is dependent on a combination of the difference between flexibilities in model parameters of the two retention mechanisms and the numbers of tissues sharing each retention mechanism (supplementary fig. S6, Supplementary Material online). In particular, the retention mechanism with greater flexibility in model parameters is most frequently chosen by CLOUD unless a majority of tissues share the more constrained retention mechanism. For example, when Tissue Mechanism A is “Conserved” (the most constrained retention mechanism) and shared by $$ k\in \{1,2,3,4,5\}$$ tissues, and Tissue Mechanism B is any other (more flexible) retention mechanism and shared by $$ 6-k$$ tissues, Tissue Mechanism B is most frequently chosen by CLOUD unless the “Conserved” retention mechanism is shared by the majority (either four or five) of the tissues. Moreover, an intriguing pattern emerges when both retention mechanisms have the same flexibility, but differ in directionality of expression divergence, that is, when one is “Neofunctionalized parent” and the other is “Neofunctionalized child.” In such cases, one of these mechanisms is still chosen most often when it is shared by five of the six tissues. However, for all other scenarios, “Specialized” is the most prominently inferred retention mechanism. This result is sensible, as “Neofunctionalized parent” is characterized by a change in the parent copy, “Neofunctionalized child” by a change in the child copy, and “Specialized” by unconstrained changes in both parent and child copies. In summary, if two different retention mechanisms are each shared by a subset of the tissues, then the more flexible retention mechanism will be predominantly chosen unless the more constrained retention mechanism is shared by the majority of the tissues.

Despite its difficulty in classifying retention mechanisms under these mixed scenarios, CLOUD is generally unbiased in its parameter predictions (supplementary figs. S7–S11, Supplementary Material online). Though mixtures of retention mechanisms occasionally lead to biased parameter predictions, such as the overestimation of $$ {\sigma }^2$$ for the setting of one “Conserved” tissue and five “Specialized” tissues (supplementary fig. S7, Supplementary Material online), or the underestimation of α for the setting of three “Conserved” and “Subfunctionalized” tissues (supplementary fig. S9, Supplementary Material online), estimates of expression states ($$ {\theta }_{\mathrm{P}},\hspace{0.17em}{\theta }_{\mathrm{C}}$$, and $$ {\theta }_{\mathrm{A}}$$) are always unbiased. Thus, the parameters supporting specific retention mechanisms are well-estimated on average (supplementary figs. S7–S11, Supplementary Material online). These results highlight that under complex settings, CLOUD should accurately predict underlying evolutionary parameters, even in situations for which no single classified retention mechanism is correct or sensible (supplementary fig. S6, Supplementary Material online).

Application of CLOUD to Empirical Data from Drosophila

Our simulation experiments highlight the exceptional classification performance of CLOUD relative to CDROM, as well as the unique ability of CLOUD to predict parameters underlying the evolution of duplicate genes. Hence, we next sought to use CLOUD to classify retention mechanisms and predict parameters of 208 duplicate genes in Drosophila (Assis and Bachtrog 2013) from their expression data in six tissues (Assis 2019b). Specifically, we first used PhyML (Guindon et al. 2010) to estimate a gene tree relating each parent, child, and ancestral gene in this data set of duplicate genes (Assis and Bachtrog 2013) (see Materials and Methods section). Next, as in our simulation studies, we trained CLOUD (assuming two hidden layers) on a large balanced simulated training set of 50,000 observations (10,000 from each of five classes), with evolutionary parameters $$ {\theta }_{\mathrm{P}},{\theta }_{\mathrm{C}},{\theta }_{\mathrm{A}}\in [-4,4],\hspace{0.17em}\hspace{0.17em}{\text{log}}_{10}(\alpha )\in [0,3]$$, and $$ \hspace{0.17em}{\text{log}}_{10}({\sigma }^2)\in [-2,3]$$ drawn independently across six tissues, for a total of 30 random parameters per simulated training observation (see Materials and Methods section). We tailored CLOUD to this data set of duplicate genes (Assis and Bachtrog 2013) by generating p = 108 input features (table 1) from comparisons to the empirical distribution of single-copy genes identified in this lineage (Assis and Bachtrog 2013) (see Materials and Methods section). Then, we used CLOUD to classify retention mechanisms and predict parameters of the 208 Drosophila duplicate genes (Assis and Bachtrog 2013) from these features.

Analysis of the resulting classifications reveals that the predominant mechanism of duplicate gene retention in Drosophila is neofunctionalization in which the child copy acquires a new function (61.43%; fig. 5), mirroring the findings of Assis and Bachtrog (2013). Moreover, classifications of CLOUD are generally concordant with those of CDROM (59.29%), with three key differences. In particular, of the 167 duplicates classified as “Neofunctionalized child” by CDROM, 16 are classified as “Conserved” by CLOUD. In addition, of the 53 duplicates classified as “Conserved” by CDROM, 18 are classified as “Neofunctionalized child” and 14 as “Specialized” by CLOUD. Finally, of the 41 duplicates classified as “Specialized” by CDROM, 18 are classified as “Neofunctionalzied child” by CLOUD. Based on our simulation results (supplementary figs. S1 and S2, Supplementary Material online), it is likely that these discrepancies reflect differences in the abilities of the CLOUD and CDROM classifiers to handle gene expression stochasticity.

Fig. 5

Classification and prediction results for application of CLOUD with L = 2 hidden layers to empirical data from Drosophila (Assis and Bachtrog 2013; Assis 2019b). Box plots overlaid onto strip plots show distributions of log-transformed parameter estimates for each of the five classes of duplicate gene retention mechanisms. Note that six estimates, corresponding to the six tissues in the empirical data set, are plotted for each parameter. The confusion matrix in the bottom right illustrates the high concordance in classifications of CLOUD and CDROM for these empirical data, with both methods classifying the majority of duplicate genes as retained by neofunctionalization of the child copy.

We next examined the parameter predictions of CLOUD. Here, our major question was whether these predictions match theoretical expectations of duplicate gene retention mechanisms. To answer this question, we examined distributions of empirical parameter estimates for each class obtained with the CLOUD classifier (fig. 5). We first considered optimal expression estimates $$ {\theta }_{\mathrm{P}},\hspace{0.17em}{\theta }_{\mathrm{C}}$$, and $$ {\theta }_{\mathrm{A}}$$. For the “Conserved” class, the distributions of estimated $$ {\theta }_{\mathrm{P}},\hspace{0.17em}{\theta }_{\mathrm{C}}$$, and $$ {\theta }_{\mathrm{A}}$$ are not significantly different from one another, consistent with expectations (table 2). For the “Neofunctionalized parent” class, the distribution of $$ {\theta }_{\mathrm{P}}$$ is different (though not significantly) from those of $$ {\theta }_{\mathrm{C}}$$ and $$ {\theta }_{\mathrm{A}}$$, whereas the distributions of $$ {\theta }_{\mathrm{C}}$$ and $$ {\theta }_{\mathrm{A}}$$ are not significantly different from one another. This qualitative pattern is as expected (table 2), with the lack of a significant difference of $$ {\theta }_{\mathrm{P}}$$ likely due to the small number of duplicate genes in this class. For the “Neofuntionalized child” class, the distribution of $$ {\theta }_{\mathrm{C}}$$ is significantly different from those of $$ {\theta }_{\mathrm{P}}$$ and $$ {\theta }_{\mathrm{A}}$$, whereas the distributions of $$ {\theta }_{\mathrm{P}}$$ and $$ {\theta }_{\mathrm{C}}$$ are not significantly different from one another, also consistent with theoretical expectations (table 2). It is also interesting that, relative to other values of θ, $$ {\theta }_{\mathrm{P}}$$ is increased in the “Neofunctionalized parent” class, whereas $$ {\theta }_{\mathrm{C}}$$ is decreased in the “Neofunctionalized child” class. Though the sample size for the “Neofunctionalized parent” class is small, decreased $$ {\theta }_{\mathrm{C}}$$ in the “Neofunctionalized child” class is consistent with lower gene expression levels of testis-specific genes (Brawand et al. 2011; Assis and Bachtrog 2013), which compose a majority of the child duplicate gene copies in our data set (Assis and Bachtrog 2013). For the “Subfuntionalized” class, distributions of $$ {\theta }_{\mathrm{P}}$$ and $$ {\theta }_{\mathrm{C}}$$ are different (though not significantly) from one another, with the center of the distributions of $$ {\theta }_{\mathrm{P}}$$ and $$ {\theta }_{\mathrm{C}}$$ located at approximately the center of the $$ {\theta }_{\mathrm{A}}$$ distribution. This qualitative pattern matches expectations, though formal tests of significance were again underpowered due to only a handful of duplicates classified as “Subfunctionalized.” Finally, for the “Specialized” class, $$ {\theta }_{\mathrm{P}},\hspace{0.17em}{\theta }_{\mathrm{C}}$$, and $$ {\theta }_{\mathrm{A}}$$ are all significantly different from one another, matching theoretical expectations (table 2). Analogously, we also observe a general concordance between estimates of α and $$ {\sigma }^2$$ and theoretical expectations of classified duplicate gene retention mechanisms. In particular, duplicate genes classified as neofunctionalized or specialized have significantly elevated estimated selection strengths (α) compared to those classified as conserved or subfunctionalized (fig. 5). These differences are consistent with theoretical expectations, as both neofunctionalization and specialization result in acquisitions of new functions that are hypothesized to be driven by strong selection, whereas both conservation and subfunctionalization result in preservations of ancestral functions that may occur in the absence of selection. Further, estimates of phenotypic drift ($$ {\sigma }^2$$) are also significantly larger for duplicate genes classified as neofunctionalized or specialized than as conserved or subfunctionalized. This result supports the hypothesis that traits require some minimum threshold of plasticity to effectively explore the space of novel states on which selection may act.

As a final analysis of the empirical data, we performed a case study of the child gene Dntf-2r and its parent Dntf-2. We chose this duplication event because it was well-characterized in earlier studies (Betrán and Long 2003; Assis and Bachtrog 2013; Jiang and Assis 2017), providing us with a baseline for comparing our findings. In particular, Dntf-2r arose in the D. melanogaster lineage after its divergence from the D. pseudoobscura lineage (Betrán and Long 2003; Assis and Bachtrog 2013). Several studies showed that Dntf-2r is specifically expressed in the testis and evolving under positive selection, whereas its parent Dntf-2 is expressed broadly across tissues and evolving under negative selection (Bhattacharya and Steward 2002; Betrán and Long 2003; Assis and Bachtrog 2013; Jiang and Assis 2017). Hence, it has been hypothesized that Dntf-2r underwent neofunctionalization and acquired a new male-specific function after duplication (Betrán and Long 2003; Assis and Bachtrog 2013). Consistent with this hypothesis, Dntf-2r and Dntf-2 are classified by both CDROM and CLOUD as retained by neofunctionalization of the child copy. Moreover, CLOUD estimates mean $$ \hspace{0.17em}{\text{log}}_{10}(\alpha )=4.254$$ and mean $$ \hspace{0.17em}{\text{log}}_{10}({\sigma }^2)=-3.429$$, supporting previous findings that strong selection and weak phenotypic drift underlie neofunctionalization of Dntf-2r (Betrán and Long 2003; Jiang and Assis 2017).

Training the Neural Network on Data Simulated from OU Processes

Consider a set of N_k training observations for class $$ k\in \{1,2,3,4,5\}$$, such that the total training sample size is $$ N=N_1+N_2+N_3+N_4+N_5$$. For observation $$ i\in \{1,2,\dots ,N\}$$ and output $$ k\in \{1,2,\dots ,K\}$$, let $$ y_k^{(i)}$$ denote the true value and $$ {\widehat{y}}_k^{(i)}$$ denote the estimated value. We wish to train a neural network model to minimize the overall discrepancy between $$ y_k^{(i)}$$ and $$ {\widehat{y}}_k^{(i)}$$, which we measure with the loss function $$ ℒ({\widehat{\mathbf{y}}}^{\left(i\right)},{\mathbf{y}}^{\left(i\right)})$$, across the N samples and K outputs. Let

To train the neural network, we wish to identify the set of parameters (weights and biases) $$ \mathcal{W}=\{\mathbf{w},{\mathbf{W}}^{[0]},\dots ,{\mathbf{W}}^{[L]}\}$$ that minimize the cost

To prevent overfitting, we include an elastic net-style regularization penalty term (Zou and Hastie 2005) on the weights of each layer with two tuning hyper parameters. Specifically, we reduce the complexity of the fitted model with the tuning parameter $$ \lambda \ge 0$$, which shrinks the magnitude of the weights to zero. We also perform simultaneous weight shrinkage and feature selection with the elastic net tuning parameter $$ \gamma \in [0,1]$$, such that we are performing L₂-norm regularization when γ  =  0, L₁-norm regularization that incorporates feature selection when γ = 1, and both types of regularization when $$ \gamma \in (0,1)$$. In particular, we seek to find the model parameters $$ \mathcal{W}$$ that minimize the penalized cost function

In the classification problem, for training observation $$ i\in \{1,2,\dots ,N\}$$, we define the indicator random variable $$ y_k^{(i)}=1$$ if observation i is from class k, and 0 otherwise. Hence, all output values are zero except for that corresponding to class k, which has a value of one. Based on this output, we employ the loss function used in $$ J(\mathcal{W},L,\lambda ,\gamma )$$ as the categorical cross entropy deviance (Goodfellow et al. 2016)

In the prediction problem, $$ y_k^{(i)}$$ is instead the kth parameter value from simulated replicate i. Based on this output, we employ the loss function used in $$ J(\mathcal{W},L,\lambda ,\gamma )$$ as the residual sum of squared error

Simulated data have been successfully used to train models for learning about evolution from genomic data in many recent studies (Lin et al. 2011; Schrider and Kern 2016; Sheehan and Song 2016; Kern and Schrider 2018; Schrider et al. 2018; Sugden et al. 2018; Flagel et al. 2019; Mughal et al. 2020; Mughal and DeGiorgio 2019; Adrion et al. 2020). Therefore, to train CLOUD for both classification and prediction, we generated a balanced simulated data set with 10⁴ observations from each of the five classes, for a total of N = 50,000 training observations. We assumed that tissues were independent, and that there were a total of m = 6 tissues as in an empirical data from Drosophila (Assis 2019b) that we later applied our method to, for a total of p = 108 input features.

To make the simulated data set more realistic, we drew model parameters $$ T_{\text{PC}}$$ and $$ T_{\text{PCA}}$$ from empirical gene tree estimates for the set of Drosophila duplicate genes used by Assis and Bachtrog (2013). The procedure for estimating these gene trees is detailed in subsection Application of CDROM and CLOUD to empirical data from Drosophila below. For all analyses, we scaled the root of the gene tree to have height one, and considered a new scaled time for the duplication event of $$ t_{\text{PC}}=T_{\text{PC}}/T_{\text{PCA}}$$, such that $$ t_{\text{PC}}$$ represented the time of the duplication relative to the height of the root of the gene tree. For a given class, we drew parameters $$ \Omega ={\{t_{\text{PC}},{\theta }_{\mathrm{P}j},{\theta }_{\mathrm{C}j},{\theta }_{\mathrm{A}j},{\alpha }_j,{\sigma }_j^2\}}_{j=1}^m$$ uniformly at random, assuming that $$ {\theta }_{\text{PC}j}={\theta }_{\text{PCA}j}={\theta }_{\mathrm{A}j}$$ for tissue j (schematic provided in fig. 2D). In particular, we drew $$ t_{\text{PC}}$$ from the distribution of empirical gene tree estimates, $$ {\theta }_j\in [-4,4]$$ for $$ j\in \{\mathrm{P},\mathrm{C},\mathrm{A}\}$$, α from $$ \hspace{0.17em}{\text{log}}_{10}(\alpha )\in [0,3]$$, and $$ {\sigma }^2$$ from $$ \hspace{0.17em}{\text{log}}_{10}({\sigma }^2)\in [-2,3]$$. We chose this specific range for $$ {\theta }_{\mathrm{P}},\hspace{0.17em}{\theta }_{\mathrm{C}}$$, and $$ {\theta }_{\mathrm{A}}$$ because we found that differences in $$ \hspace{0.17em}{\text{log}}_{10}$$-transformed empirical expression data were normally distributed, matching expectations under an OU model. For this reason, all empirical expression data were also $$ \hspace{0.17em}{\text{log}}_{10}$$-transformed prior to applying CLOUD. The class k is determined by $$ \{{\theta }_{\mathrm{P}},{\theta }_{\mathrm{C}},{\theta }_{\mathrm{A}}\}$$, which is summarized in table 2. Then, we simulated gene expression data $$ {\mathbf{e}}^{(i)}\in {ℝ}^{3m}$$ under model parameters for a given class k (table 2), assuming independence among tissues, and generated N_k simulated replicates of parameter values $$ {\Omega }_k^{(i)}$$ for $$ i=1,2,\dots ,N_k$$.

Given a number of hidden layers L and the pair of regularization tuning parameters λ and γ, we estimated the set of parameters $$ \mathcal{W}$$ using the Adam optimizer (Kingma and Ba 2014) with learning rate $$ {10}^{-3}$$ and exponential decay rates for the first and second moment estimates of $$ {\beta }_1=0.9$$ and $$ {\beta }_2=0.999$$ (Kingma and Ba 2014), respectively. This optimizer was used as it efficiently traverses the cost function surface $$ J(\mathcal{W},L,\lambda ,\gamma )$$ to rapidly identify the minimum (Kingma and Ba 2014). We also used mini-batch optimization with a batch size of 5,000 observations for 500 epochs. To estimate L, λ, and γ, we performed five-fold cross-validation (Hastie et al. 2009). Specifically, we used 80% (40,000) of observations for training, and held out the remaining 20% (10,000) for validation. We also ensured that the training and validation sets were balanced in class representation, such that there were equal numbers of observations from each class in the training (8,000) and validation (2,000) sets. To assess method performance for a given fold, we computed the validation loss

Previous studies have found that neural networks with enough hidden layers or units can approximate any function, and therefore lead to overfitting (Cybenko 1989; Goodfellow et al. 2016). Hence, based on simulations, we estimated that a neural network with $$ \widehat{L}=2$$ hidden layers provides the best cross-validation performance, with the validation loss for the classifier of approximately 0.918 with optimal tuning parameters $$ \widehat{\lambda }\approx 1.778\times {10}^{-4}$$ and $$ \widehat{\gamma }=1$$, and the validation loss for the predictor of approximately 0.899 with optimal tuning parameters $$ \widehat{\lambda }\approx 7.499\times {10}^{-8}$$ and $$ \widehat{\gamma }=0.8$$. Comparisons of classification and prediction performances across the four network architectures $$ L\in \{0,1,2,3\}$$ are given in supplementary figures S12 and S13, Supplementary Material online, highlighting the generally superior performance of the architecture with two hidden layers.

Application of CDROM and CLOUD to Simulated Test Data

To compare the relative classification powers and accuracies of CDROM and CLOUD and explore the prediction accuracy of CLOUD, we simulated training and test data sets for duplicate genes based on an OU process, which is described in subsection Training the neural network on data simulated from OU processes above. However, in that subsection, we assumed that expression vectors for single-copy genes $$ \mathcal{G}$$ were given. These would typically be extracted from the genome-wide distribution of single-copy genes for the pair of species being studied, such that trained models are based on the level of expression divergence typically observed in the study system.

For our simulated training and test sets, we generated a background set of 10,000 six-tissue expression vectors for single-copy genes that was inspired by those of the single-copy genes identified in Drosophila (Assis and Bachtrog 2013; Assis 2019b). Specifically, we applied the Brownian motion model (Felsenstein 1973) implemented in mvmorph (Clavel et al. 2015) to expression vectors of single-copy genes between Species 1 and 2, assuming that tissues were independent and that the root of the two-species phylogeny had height one, to estimate ancestral expression θ and variance $$ {\sigma }^2$$ parameters consistent with the empirical distribution in Drosophila at each tissue and single-copy gene. Given the set of parameters, we then sampled values of θ and $$ {\sigma }^2$$ uniformly at random from the estimated empirical distribution, and generated simulated single-copy expression vectors in Species 1 and Species 2 for m = 6 independent tissues, giving us the simulated set $$ \mathcal{G}$$.

To test either the classifier or predictor, we generated a balanced set of duplicate gene expression vectors, such that each of the five classes had 1,000 observations, for a total of N = 5,000 test observations. We assumed that tissues were independent, and that there was a total of m = 6 tissues as in the training set. For a given class, we drew parameters $$ \Omega ={\{t_{\text{PC}},{\theta }_{\mathrm{P}j},{\theta }_{\mathrm{C}j},{\theta }_{\mathrm{A}j},{\alpha }_j,{\sigma }_j^2\}}_{j=1}^m$$ uniformly at random. In particular, as with the training set, we drew $$ t_{\text{PC}}$$ from the distribution of empirical gene tree estimates, $$ {\theta }_j\in [-4,4]$$ for $$ j\in \{\mathrm{P},\mathrm{C},\mathrm{A}\}$$, α from $$ \hspace{0.17em}{\text{log}}_{10}(\alpha )\in [0,3]$$, and $$ {\sigma }^2$$ from $$ \hspace{0.17em}{\text{log}}_{10}({\sigma }^2)\in [-2,3]$$. The class k was determined by $$ \{{\theta }_{\mathrm{P}},{\theta }_{\mathrm{C}},{\theta }_{\mathrm{A}}\}$$ (table 2), and gene expression data were generated $$ {\mathbf{e}}^{(i)}\in {ℝ}^{3m}$$ under model parameters for a given class k (table 2), assuming independence among tissues and N_k simulated replicates of parameter values $$ {\Omega }_k^{(i)}$$ for $$ i=1,2,\dots ,N_k$$.

To assay how CLOUD performs in different portions of the parameter space, we also examined its accuracy on test sets drawn from restricted parameter values. Specifically we considered three distinct ranges for α of $$ [1,10],\hspace{0.17em}[10,100]$$, and[100, 1,000], and five distinct ranges for $$ {\sigma }^2$$ of [0.01, 0.1], [0.1, 1], [1, 10], [10, 100], and [100, 1,000]. For each combination of a range for α and a range for $$ {\sigma }^2$$, we sampled α and $$ {\sigma }^2$$ uniformly at random.

Application of CDROM and CLOUD to Empirical Data from Drosophila

We applied CDROM and CLOUD to empirical data consisting of Drosophila duplicate and single-copy genes (Assis and Bachtrog 2013) and their expression abundances in six tissues (Assis 2019b). In particular, duplicate and single-copy genes in D. melanogaster and D. pseudoobscura were obtained from Assis and Bachtrog (2013). In that study, pairs of duplicate genes in each species were identified via BLAST searches (Altschul et al. 1990) and supplemented with those from Chen et al. (2010). A table of orthologs, or genes that arose from the same common ancestor in 12 sequenced Drosophila species (Drosophila 12 Genomes Consortium 2007), was downloaded from FlyBase at https://www.flybase.org. Orthologs were used to determine presence or absence of duplicate gene copies in D. melanogaster and D. pseudoobscura. Gene duplication events that occurred after the divergence of the D. melanogaster and D. pseudoobscura lineages were defined as those for which one duplicate gene copy is present in both species (parent) and the second is only present in one species (child). Quantile-normalized gene expression abundances for carcass, female head, ovary, male head, testis, and accessory gland tissues in D. melanogaster and D. pseudoobscura were obtained from the Dryad data set associated with Assis (2019b) at https://doi.org/10.5061/dryad.742564m. In that study, paired-end RNA-sequencing reads were downloaded from modENCODE (Celniker et al. 2009) at https://www.modencode.com, and aligned to reference transcriptomes of each species with Bowtie 2 (Langmead et al. 2009). Abundances in fragments per kilobase of exon per million fragments mapped (FPKM, Trapnell et al. 2013) were calculated with eXpress (Roberts and Pachter 2013), quantile-normalized, and log-transformed. After examination of the distribution of these values, genes with little or no expression in all tissues were removed.

Because CLOUD requires estimates of $$ T_{\text{PC}}$$ and $$ T_{\text{PCA}}$$, we first generated multiple sequence alignments with MACSE (Ranwez et al. 2018), which accounts for underlying codon structure, and then inferred a gene tree with PhyML (Guindon et al. 2010) for each parent, child, and ancestral gene in the duplication data set (Assis and Bachtrog 2013). The empirical distributions of estimated $$ T_{\text{PC}}$$ and $$ T_{\text{PCA}}$$ across these gene trees were used as input to an OU process to generate a balanced training set with N = 50,000 observations as described in subsection Training the neural network on data simulated from OU processes above. Gene expression data from single-copy genes in Drosophila (Assis and Bachtrog 2013; Assis 2019b) were used as the set $$ \mathcal{G}$$ necessary for application of both CDROM and CLOUD. We trained a classifier and predictor for CLOUD assuming L = 2 hidden layers, and through five-fold cross-validation (Hastie et al. 2009), estimating the regularization tuning parameters as $$ \widehat{\lambda }\approx 1.778\times {10}^{-4}$$ and $$ \widehat{\gamma }=0.9$$ for the classifier and $$ \widehat{\lambda }={10}^{-6}$$ and $$ \widehat{\gamma }=0.5$$ for the predictor.