A likelihood ratio test for changes in homeolog expression bias

Background Gene duplications are a major source of raw material for evolution and a likely contributor to the diversity of life on earth. Duplicate genes (i.e., homeologs, in the case of a whole genome duplication) may retain their ancestral function, sub- or neofunctionalize, or be lost entirely. A primary way that duplicate genes evolve new functions is by altering their expression patterns. Comparing the expression patterns of duplicate genes gives clues as to whether any of these evolutionary processes have occurred. Results We develop a likelihood ratio test for the analysis of the expression ratios of duplicate genes across two conditions (e.g., tissues). We demonstrate an application of this test by comparing homeolog expression patterns of 1448 homeologous gene pairs using RNA-seq data generated from leaves and petals of an allotetraploid monkeyflower (Mimulus luteus). We assess the sensitivity of this test to different levels of homeolog expression bias and compare the method to several alternatives. Conclusions The likelihood ratio test derived here is a direct, transparent, and easily implemented method for detecting changes in homeolog expression bias that outperforms alternative approaches. While our method was derived with homeolog analysis in mind, this method can be used to analyze changes in the ratio of expression levels between any two genes in any two conditions.


Background
Gene duplications are a major source of raw material for evolution and a likely contributor to the diversity of life on earth [1][2][3][4][5][6][7][8][9]. Gene duplications are a special type of mutation resulting in the multiplication of intact functional components. These duplicate genes may either retain the ancestral function or individual portions of the gene's ancestral function may be partitioned (i.e., subfunctionalize) or evolve new functions entirely (i.e., neofunctionalize) [10][11][12]. Duplicate genes may evolve new functions either by changes in the primary coding sequence or altering where and when they are expressed. Previous work has indicated that changes to gene expression and their regulatory networks may be more important, rapid, or flexible than divergence of protein identities in the evolution of sub-and neofunctionlization [13][14][15][16][17][18][19]. There are multiple scenarios in which genes can be duplicated, ranging from small regional gene duplications to massive whole genome duplications (WGDs). The term polyploid refers to cells or organisms that have undergone a WGD event and contain more than two paired sets of chromosomes. Each complete set of chromosomes is referred to as a subgenome. Homologous genes located on separate subgenomes are referred to as homeologs.
WGDs are especially common in plants; indeed, all extant angiosperms (i.e., flowering plants) have at least two rounds of WGD in common [20], and up to 15% of speciation events in angiosperms may have been the product of WGDs [21]. Importantly, many major crops (corn, potato, wheat, etc.) are polyploid [22]. WGD events and the resulting polyploidy are not restricted to plants, but have occurred in both vertebrate and invertebrate lineages as well. For example, the African clawed frog, Xenopus, commonly used as an experimental model system and extensively studied in developmental biology, includes species ranging from diploid to dodecaploid [23].
Other examples of polyploids with ancient WGD events include the zebrafish Danio rerio [24], several salmonids [2], and some species of fungi [25]. Interestingly, there exists at least one polyploid mammal [26], a tetraploid rat from Argentina that mediates gene dosage by regulation of ribosomal RNA.
The biological consequences of gene duplications and subfunctionalization are significant and include examples such as the evolution of eyes [27], the evolution of hemoglobins [28], development of heat resistance in plants [29], and insecticide resistance [30]. Given the importance of duplicate genes in evolution, it is natural to ask how we might quantify differences in the activity or function of homeologous genes. One way to begin exploring this question is by analyzing gene expression levels.
Genome-wide gene expression levels are commonly quantified using high throughput RNA sequencing (RNAseq) [31]. In RNA-seq experiments, mRNA is extracted, purified, and reverse transcribed into cDNA. This cDNA is fragmented into smaller pieces and sequenced using next-generation technology. The resulting millions of sequence reads are then mapped to either a reference genome or reference transcriptome, and the number of sequences mapping to a particular gene is used as an indication of the expression level of that gene.
In differential expression analysis, high-throughput RNA-seq data is used to determine if gene expression levels vary under different experimental conditions, or in distinct tissues, etc. Several different approaches to this statistical analysis exist [32][33][34], some of which use methods based on maximum likelihood estimation and likelihood ratio tests.
Homeologous gene pairs frequently have distinguishing sequence differences. Therefore, sequencing reads derived from individual homeologs can be distinguished and expression levels can be determined for each homeolog. The term homeolog expression bias (HEB) refers to cases where homeologs are expressed at unequal levels in a single experimental condition [35]. The primary objective of this paper, development of a likelihood ratio test for statistical analysis of changes in homeolog expression bias (denoted HEB) is a non-trivial extension of the statistical analysis of differential expression.
The following sections begin with the derivation of a likelihood ratio test for HEB. This is our starting point for the development of a likelihood ratio test for HEB, i.e. changes in relative expression levels between homeologous genes in two conditions. We apply this method to RNA-seq data of homeologous gene expression in petals and leaves of the allotetraploid Mimulus luteus. Finally, using simulated data, we show that the likelihood ratio test for HEB derived here is the best choice among several alternative methods.

Quantifying homeolog expression bias (HEB)
We will write A and B to denote a homeologous gene pair from which RNA-seq data is generated in n biological replicates. Typically, the mean expression levels of the homeologs (denotedā andb) are normalized by gene length and sequencing depth, as when reported in units of RPKM (reads per kilobase of coding sequence per million mapped reads). We define the homeolog expression bias (HEB) of the n replicates as a dimensionless quantity with HEB = 0 indicating no bias. If one uses the base 2 logarithm, HEB = −3 indicates 8fold bias towards homeolog A.

Likelihood ratio test for HEB
After accounting for the possibility of different gene lengths, the statistical test for HEB is essentially a likelihood ratio test for differential expression of a pair of homeologous genes. The goal is to determine whether there is sufficient evidence to reject the null hypothesis (H 0 ) that there is no bias (i.e., equal expression levels for homeologous genes) in favor of the alternative hypothesis (H 1 ) that bias is present, i.e., different expression levels for homeologous genes. In mathematical terms, the null hypothesis H 0 corresponds to the parameters (denoted by θ) of a probability model for generating the data being in a specified subset 0 of the parameter space , that is, Let θ = λ a , λ b denote the true but unknown expression levels (physical units of length −1 , e.g., RPKM). Assuming positive, i.e. non-zero, expression, the parameter space is = θ : λ a , λ b ∈ R + . The null (H 0 ) and alternative (H 1 ) hypotheses for the likelihood ratio test for homeolog expression bias are formalized as follows, Equivalently, let ω = λ b /λ a denote the ratio of expression levels and drop the superscript indicating the reference homeolog (λ = λ a ). In that case, λ b = ωλ and the hypotheses are written as follows, Once we specify a probability model for the data X , likelihood functions for each hypothesis, L 0 (θ|X ) and L 1 (θ|X ), can be derived (see next section). For composite hypotheses, the appropriate likelihood ratio test statistic is whereL 1 andL 0 are the maximized likelihoods, A critical value of the test statistic (W * ) is obtained from the Chi-squared distribution with significance level α = 0.05. The number of degrees of freedom δ is the difference in the number of free parameters in and 0 (here δ = 1) [36]. The null hypothesis H 0 is rejected in favor of the alternative H 1 when W (X ) > W * .

Probability model for RNA-seq read counts
Denote the lengths of homeologous genes a and b as a and b (e.g., in kilobases) and let d i be the sequencing depth (e.g., in millions of mapped reads) of replicate i. The expected number of RNA-seq reads for gene a and replicate i is where in the second equality we have dropped the superscript for the reference homeolog (λ = λ a ). Similarly, the expected number of RNA-seq reads for gene b and replicate i is where ω = λ b /λ a = λ b /λ. In order to model the overdispersion commonly observed in RNA-seq data, the probability model assumes that the count data for each gene is drawn from a negative binomial distribution, where μ is the appropriate mean (μ a i or μ b i in Eqs. 2 and 3). That is, if X a i and X b i are random variables representing the count data for replicate i of homeologous genes A and B, In these expressions, the aggregation parameter r i is obtained from the observed mean-variance relation for all homeolog pairs of the ith experimental replicate (see Appendix 1).
Assuming independence of experimental replicates, the likelihood functions L 1 and L 0 are products of the likelihood functions for each observation, that is, and similarly for L 0 (X ), where X i = {a i , b i } indicates the observed read counts for replicate i and X = ∪ n i=1 X i . The likelihood function for the alternative hypothesis and the ith replicate is The likelihood function for the null hypothesis and the ith replicate, L i 0 (X ), is given by Eq. 4 with ω = 1.

Maximum likelihood estimation
Maximum likelihood estimation is performed using the the log-likelihood function corresponding to Eq. 4, namely, where and γ (·) = ln (·). The log-likelihood function for the null hypothesis (ln L 0 ) is given by Eq. 6 with ω = 1. The log-likelihood function ln L 1 (X ) is maximized by numerically solving forλ andω leading to zero partial derivatives, as described in Appendix 2. The log-likelihood function ln L 0 (X ) is maximized by solving forλ leading to The optimal parameter valuesλ andω are used to evaluate lnL 0 (X ;λ), lnL 1 (X ;λ,ω), and the test statistic W (see Eq. 1).

Quantifying changes in homeolog expression bias ( HEB)
Let A and B represent homeologous genes and RNAseq data is generated under conditions 1 and 2 in n biological replicates, leading to mean expression levels a 1 ,ā 2 ,b 1 ,b 2 . The change in homeolog expression bias ( HEB) is defined as where the last equality uses HEB 1 = logb 1 /ā 1 and HEB 2 = logb 2 /ā 2 .

Likelihood ratio test for HEB
The likelihood ratio test for HEB is designed to determine whether there is sufficient evidence to reject the null hypothesis (H 0 ) that homeolog expression bias is the same under two experimental conditions ( HEB = 0) in favor of the alternative hypothesis (H 1 ) that there is a difference in bias ( HEB = 0). Following notation similar to the previous section, our hypotheses are The difference in degrees of freedom of the alternative and null hypotheses is δ = 4 − 3 = 1.
The likelihood functions for the HEB test are similar to those for HEB, though the two different experimental conditions lead to twice as many terms (cf. Eq. 4). The likelihood function for H 1 is where L k,i 1 , the likelihood function for the ith replicate of the kth condition, has the form of Eq. 4 with parameters indexed by condition a k,i , b k,i , r a k,i , r b k,i , ω k . The loglikelihood function for H 1 is thus where ln L k,i and γ (·) = ln (·). The log-likelihood function for the null hypothesis (ln L 0 ) is given by the above expressions with ω 1 = ω 2 = ω. The aggregation parameters (r k,i ) are determined from the data with experimental conditions k = 1 and 2 considered separately (cf. Eqs. 17-19).
The log-likelihood function ln L 1 (X ) used in the analysis of HEB is maximized by numerically solving uncoupled systems of the form of Eqs. 7 and 8 for λ 1 ,ω 1 and λ 2 ,ω 2 . The log-likelihood function ln L 0 (X ) is maximized by solving forλ 1 ,λ 2 andω that lead to zero partial derivatives, The optimal parameter values are used to evaluate the likelihoods,L 0 (X ;λ 1 ,λ 2 ,ω) andL 1 (X ;λ 1 ,λ 2 ,ω 1 ,ω 2 ), and the test statistic W (see Eq. 1). The numerical solution of these equations was facilitated by transforming these equations in a manner that ensured both parameters are positive and symmetric with respect to the mean expression levels of homeolog A and B (see Appendix 2).

Mimulus luteus
To demonstrate the application of the likelihood ratio test for HEB, five biological replicates of RNA-seq data were generated from petals of the tetraploid Mimulus luteus (monkeyflower), and another five replicates were generated from the leaves (see Appendix 3 for details). We have chosen M. luteus because it is a tetraploid with two distinct subgenomes, denoted A and B (mean synonymous divergence is ∼ 11.1%; for details on genome assembly, see [37]). In this section, we use the likelihood ratio test for HEB to find homeologous gene pairs where one homeolog is expressed at significantly different levels than the other, one tissue at a time. In the section on HEB we develop a likelihood ratio test to determine whether there is a significant difference in the bias between the two tissues. Figure 1 (top panel) shows the result of applying the likelihood ratio test for HEB to the petal data. There are 1853 homeologous gene pairs in M. luteus that can be identified as coming from separate subgenomes. Of these 1853 homoeologous pairs, 1560 were testable (measurable expression from each individual homeolog). Of testable pairs, a total of 676 gene pairs show significant bias (using a significance level of α = 0.05, and applying the Benjamini-Hochberg correction [38,39] to account for multiple testing error). In the 334 pairs biased towards the A homeolog the mean HEB is − 2.49 (5.6-fold change). In the 342 pairs biased towards the B homeolog, the mean HEB is 2.39 (5.2-fold change).

Homeolog expression bias in Mimulus luteus petals
These results may be indicative of a number of evolutionary processes. For example, one of the homeologs may have become sub-or neofunctionalized in this tissue, or one of the homeologs may simply be losing its function.

Homeolog expression bias in Mimulus luteus leaves
Next, the likelihood ratio test for HEB was applied to the leaf data (results shown in Fig 1, bottom panel). Of 1853 homoeologous pairs, 1498 were testable and a total of 399 gene pairs show significant bias. In the 199 pairs biased towards the A homeolog the mean HEB is − 2.83 (7.1-fold change). In the 200 pairs biased towards the B homeolog, the mean HEB is 2.80 (7.0-fold change).

Mimulus luteus
The likelihood ratio test for HEB requires each homeolog to have at least one read in each condition. Returning to the leaf and petal data from the previous sections on HEB, this gives 1448 testable pairs. Figure 2 shows the results of the likelihood ratio test for HEB. We find a total of 76 gene pairs show significant HEB. Of these, 35 are more biased towards the A homeolog in the leaf than they are in the petal. The remaining 41 gene pairs are more biased towards the B homeolog in the leaf than in the petal. Figure 3 shows a scatter plot of homeolog expression bias (HEB) in leaf and petal. Colored marks indicate gene pairs with statistically significant changes in homeolog expression bias ( HEB) (these points correspond to the colored bars in Fig. 2). Data points in the top-left and bottom-right quadrants of Fig. 3 represent homeologous pairs where one homeolog is more highly expressed in one tissue and its partner is more highly expressed in the other tissue. The top-right and bottom-left quadrants correspond to homeologous pairs where the difference in bias favors the same homeolog but has become more extreme. Finally, all of the marks that are colored blue or yellow show significant change in bias and are candidates for tissue specific sub-or neofunctionalization.
Although the change in homeolog expression bias is defined by Eq. 10 as the log-fold change in homeolog expression bias, the intercalation of significant (yellow and blue) and not significant (gray) HEB in Fig. 3 makes it clear that statistical evidence for HEB is not reducible to the difference between HEB leaf and HEB petal (the vertical or horizontal distance to the line of slope 1 where HEB leaf = HEB petal ).
Whether or not HEB is statistically significant also depends on differences in sequencing depths, mean expression levels (e.g., lowly expressed genes are more likely to be influenced by shot noise), and ratios of gene lengths. All of these factors are considered simultaneously in the likelihood ratio test presented here. Assessing the statistical significance of HEB using sequential HEB analysis would almost certainly result in a different set of genes being called significant.

Validation of the likelihood ratio tests using simulated data
A natural question to ask about HEB and HEB is, "How large does the change in expression levels between homeologs across conditions need to be before we can detect HEB most of the time?" Unsurprisingly, this depends largely on the number of biological replicates.
To explore this question, we generated simulated data with one expression level fixed at a constant value, μ a = 100, and varied the other expression level, μ b = 2 x μ a , with x ∈[ 0, 2] in steps of 0.1. For each value of x, we generated 10,000 sets of data from a negative binomial distribution for N = 3, 6, 12 and 24 replicates. We fixed the parameter r = 10 for simplicity; this is within the range of values typically observed in RNA-seq data. Figure 4 shows the results of the likelihood ratio test for HEB on this simulated data set. We find that a 4-fold change is almost always detectable, regardless of the number of replicates. However, detecting a 2-fold change at least 95% of the time requires at least 12 replicates.
To assess the sensitivity of HEB to different levels of bias shift, we created a similar data set. This time, we set 3 of the expression levels equal μ a 1 = μ b 1 = μ a 2 = 100 , and varied the fourth; μ b 2 = 2 x μ a 2 , with x ∈[ 0, 2] in steps of 0.1. The aggregation parameter was again fixed at r = 10. For each value of x, 10,000 sets of data were generated from a negative binomial distribution for N = 3, 6, 12 and 24 replicates. Fig. 4 Ability of the likelihood ratio test for HEB to detect different levels of bias. Simulation results show the fraction of times H 0 was rejected for 10,000 trials with the given values of x and n (parameters: α = 0.05, μ a = 100, r = 10). With n ≥ 3 replicates, a 4-fold change is detectable over 95% of the time. Detecting a 2-fold change greater than 95% of the time requires at least 6 replicates Figure 5 shows the results of the likelihood ratio test for HEB on this simulated data set. The results are similar to those for HEB, with the test for HEB being slightly less sensitive than the test for HEB. For HEB, a 4-fold change in bias is detected more than 95% of the time when N ≥ 6. As with the test for HEB, the ability to detect smaller changes increases significantly with the number of replicates.

Alternative methods
Our method is transparent, derived specifically for the analysis of HEB, and requires a minimal number of assumptions; nevertheless we wished to investigate whether other methods could achieve similar results. Because we found only one method analogous to ours in the literature (Homeoroq, [40]), we developed three additional ad hoc methods. To compare these methods we With n ≥ 6 replicates, a 4-fold change is detectable over 95% of the time. However, detecting a 2-fold change more than 95% of the time requires at least 12 replicates generated simulated data sets and analyzed ROC curves. Each data set contained 10,000 gene pairs, half of which had HEB fixed at a constant value (2, 8, and 16). Three replicates were generated from negative binomial distributions, and this was repeated 50 times for each value of HEB (150 simulations total). First, we took a naive approach and performed t-tests and z-tests on the ratio of log 2 -fold changes between conditions 1 and 2. Next, we ran DESeq2 and extracted the estimated shrunken log 2 -fold changes and their standard errors, and performed a z-test (we call this method 'DEZ'). Unsurprisingly, the naive methods (standard t-and ztests) underperformed the LRT, with area under the ROC curve (ROC area) typically less than the LRT by ≈ 0.05 to 0.36.
The LRT outperformed DEZ for HEB = 8 and 16 (Fig. 7, top region). For HEB = 2, both methods performed poorly with mean ROC area = 0.58080 for the LRT, while DEZ came out slightly ahead with mean ROC area = 0.58083 (not shown). In all cases, the Homeoroq method significantly underperformed both alternatives.
The test with largest ROC area is not necessarily the best choice, for example, when an ROC curve accumulates a small area for low FPR, and a large area for high FPR. To address this, we evaluated partial ROC area for false positive rates between 0 and 0.1, as researchers typically don't accept FPR> 0.1. An example truncated ROC curve is shown in Fig. 6. By this metric, the LRT outperforms DEZ for HEB = 8 and 16, while for HEB = 2 both methods performed poorly, with DEZ marginally better (Fig. 7, bottom region). Homeoroq significantly underperformed both alternatives at all levels of HEB.

Conclusion
Gene duplication and polyploidy are extremely important factors in generating the diversity of life on earth. As Ohno stated in his seminal work on gene duplication [1], "Natural selection merely modified while redundancy created" the raw materials necessary for the diversification of life on earth.
In this paper we have developed a robust statistical framework specifically designed for the comparison of  Each trial contained 20,000 genes, half of which had HEB fixed at a constant value. Results for HEB = 2 are not shown as they were too low (mean area for LRT= 0.5616, for DEZ=0.5620). The bottom part of the plot shows ROC area constrained to false positive rates less than 0.1. In both regions, boxes indicate interquartile ranges, whiskers indicate 5th and 95th percentiles, and black lines indicate medians duplicate gene expression patterns. Importantly, this technique is consistent and reproducible. Through analysis of simulated data we have shown that these methods perform well, especially given the small sample sizes typical of RNA-seq experiments. We have shown that the ability to detect small differences in expression levels increases as a function of sample size, a fact that can be used to aid experimental design. Other authors have noted this in the context of traditional differential expression analysis and made similar recommendations [41][42][43]. Moreover, we demonstrate the usefulness of the likelihood ratio test for HEB using homeolog expression RNA-seq data derived from a polyploid plant. While we have developed this test for the purpose of analyzing changes in expression patterns of homeologous genes, we emphasize that the method is suitable for the expression analysis of any two genes (they need not be homeologs) across any two conditions.

Appendix 1: Estimation of aggregation parameters
Due to the typically small number of replicates in RNAseq experiments, accurate estimation of the aggregation parameter is not realistic on a gene-by-gene basis [34,44]. Instead, we use the mean-variance relation of a negative binomial distribution, namely, to compute an aggregation parameter r for each experimental replicate, after rescaling to account for each replicates sequencing depth. In brief, let x i j denote the count data for the jth pair of homeologous genes obtained for experimental replicate i ∈ {1, 2, . . . , n}. For each of the n replicates, we produce an auxiliary data set y i k,j by rescaling the count data for all replicates as though each were obtained in an experiment with the sequencing depth of replicate k, For each gene (j), we compute a scaled mean (μ k,j ) and variance σ 2 k,j of y i k,j over replicates (i). To obtain the aggregation parameter r k , we perform a nonlinear least squares fit of the observed mean-variance relation across all genes. That is, r k minimizes the sum of squares error,

Appendix 2: Numerical scheme for maximum likelihood estimation
For the analysis of both HEB and HEB, parameter values maximizing the likelihood functionsL 0 andL 1 were obtained using the built-in MATLAB command fsolve applied to Eqs. 7-9 and 14-16. In both cases, the numerical procedure was facilitated by changing variables from (λ, ω) to (v, y) through λ = e v−y ω = e 2y , that is, v = ln λ + y and y = (ln ω)/2. This ensures positivity of λ and ω and leads to a system of equations that is symmetric in λ a ↔ λ b . The new variable v is the logarithm of the geometric mean of the expression levels λ a = λ and that is, λ a = λ = e v−y and λ b = ωλ = e v+y . The transformed partial derivatives used to maximize the loglikelihood ln L 1 (Eqs. 7-8) are where The transformed partial derivative used to maximize ln L 0 are found by substituting y = 0 in Eq. 20, For the analysis of HEB, the partial derivatives used to maximize ln L 1 are two uncoupled systems of the form of Eq. 20-23, one for each experimental condition (k = 1 and 2), For the null hypothesis y 2 = y 1 = y we numerically solve a system of three equations, including for k = 1 and 2. These are coupled via