Joint between-sample normalization and differential expression detection through ℓ0-regularized regression

Background A fundamental problem in RNA-seq data analysis is to identify genes or exons that are differentially expressed with varying experimental conditions based on the read counts. The relativeness of RNA-seq measurements makes the between-sample normalization of read counts an essential step in differential expression (DE) analysis. In most existing methods, the normalization step is performed prior to the DE analysis. Recently, Jiang and Zhan proposed a statistical method which introduces sample-specific normalization parameters into a joint model, which allows for simultaneous normalization and differential expression analysis from log-transformed RNA-seq data. Furthermore, an ℓ0 penalty is used to yield a sparse solution which selects a subset of DE genes. The experimental conditions are restricted to be categorical in their work. Results In this paper, we generalize Jiang and Zhan’s method to handle experimental conditions that are measured in continuous variables. As a result, genes with expression levels associated with a single or multiple covariates can be detected. As the problem being high-dimensional, non-differentiable and non-convex, we develop an efficient algorithm for model fitting. Conclusions Experiments on synthetic data demonstrate that the proposed method outperforms existing methods in terms of detection accuracy when a large fraction of genes are differentially expressed in an asymmetric manner, and the performance gain becomes more substantial for larger sample sizes. We also apply our method to a real prostate cancer RNA-seq dataset to identify genes associated with pre-operative prostate-specific antigen (PSA) levels in patients.


Introduction
A fundamental problem in RNA-seq data analysis is to identify genes or exons that are differentially expressed with varying experimental conditions based on the read counts. Some widely used methods for differential expression analysis in RNA-seq data are edgeR [1,2], DESeq2 [3] and limma-voom [4,5]. In edgeR and DESeq2, the read counts are assumed to follow negative binomial seq is to employ a sample-specific scaling factor, e.g., CPM/RPM (counts/reads per million) [7], upper-quartile normalization [8], trimmed mean of M values [7] and DESeq normalization [9]. A review of normalization methods in RNA-data data analysis is given in [6]. In most existing methods for DE analysis in RNA-seq, the normalization step is performed prior to the DE detection step, which is sub-optimal because ideally normalization should be based on non-DE genes for which the complete list is unknown until after the DE analysis.
In [10], a statistical method for robust DE analysis using log-transformed RNA-seq data is proposed, where sample-specific normalization factors are introduced as unknown parameters. This allows for more accurate and reliable detection of DE genes by simultaneously performing between-sample normalization and DE detection. An 0 penalty is introduced to enforce that a subset of genes are selected as being differentially expressed. The experimental conditions are restricted to be categorical (e.g., 0 and 1 for control and treatment groups, respectively), and a one-way analysis of variance (ANOVA) type model is employed to detect differentially expressed genes across two or more experimental conditions.
In [11], the model of [10] is generalized to continuous experimental conditions, and the sparsity-inducing 0 penalty is relaxed as the 1 penalty. An alternating direction method of multipliers (ADMM) algorithm is developed to solve the resultant convex problem. Due to the relaxation of the 0 regularization, the method in [11] may not be as robust against noise and efficient in inducing sparse solutions as that in [10]. In this paper, we again generalize the model in [10] from categorical to continuous experimental conditions. But different from [11], we retain the 0 penalty in our model to efficiently induce sparsity. We formulated two hypothesis tests suited to different applications: the first hypothesis test is that considered in [10] and answers the question of whether the expression of a gene is significantly affected by any covariate; and in addition, a second hypothesis is formulated to test whether the expression of a gene is significantly affected by a particular covariate, when all other covariates in the regression model are adjusted for.
Due to the use of the 0 penalty, the resulting problem is high-dimensional, non-differentiable and non-convex.
To fit the proposed model, we study the optimality conditions of the problem and develop an efficient algorithm for its solution. We also propose a simple rule for the selection of tuning parameters. Experiments on synthetic data demonstrate that the proposed method outperforms existing ones in terms of detection accuracy when a large fraction of genes are differentially expressed in an asymmetric manner, and the performance gain becomes more substantial for larger sample sizes. We also apply our method to a real prostate cancer RNA-seq dataset to identify genes associated with pre-operative prostatespecific antigen (PSA) levels in patients.

Methods
Given m genes and n samples, let y ij , i = 1, . . . , m, j = 1, . . . , n, be the log-transformed gene expression values of the i-th gene in the j-th sample. A small positive constant can be added prior to taking the logarithm to avoid taking logarithm of zeros. We formulate the following model: where α i is the intercept, is the regression coefficient vector of the linear model for gene i, and is a vector of p predictor variables for sample j representing its experimental conditions (drug dosage, blood pressure, age, BMI, etc.), d j represents the normalization factor for sample j, and ε ij ∼ N 0, σ 2 i is i.i.d. Gaussian noise. Our goal is for each gene to determine whether its expression level is significantly associated with the experimental conditions or not.

Remark 1
The α i and d j in (1) model gene-specific factors (e.g., gene length) and sample-specific factors (i.e., sequencing depth), respectively. Thus, model (1) can accommodate any gene expression levels summarized in the form of c ij /(l i · q j ), where c ij is the read count, l i is the gene-specific scaling factor (e.g., gene length) associated with gene i and q j is the sample-specific scaling factor (e.g., sequencing depth) associated with sample j. Special cases are read count (i.e., l i = q j = 1), CPM/RPM (i.e., l i = 1) [7], RPKM/FPKM [12,13] and TPM [14].
Since the random noise in gene expression measurements are independent across genes and samples, the likelihood is given by The negative log-likelihood is where C depends on σ 2 i but not on {α i }, {β i } and {d j }. In "Maximum likelihood estimation of noise variance" section, we will describe how to estimate σ 2 i , i = 1, . . . , m. Hereafter, we assume that σ 2 i 's are known and simply denote the negative log-likelihood as l α, In practice, typically only a subset of genes are differentially expressed. We introduce a sparse penalty to the negative log likelihood function: where λ i 's are tuning parameters that control the sparsity level of the solution, and p β i is a penalty function In this paper, we use the following two types of penalty functions.
i) Type I penalty: This penalty function applies to applications where all covariates are of interest and we want to identify genes for which at least one covariate is associated with its expression. ii) Type II penalty: This penalty applies to applications where only one (the p-th) covariate is of main interest (e.g., treatment) while we want to adjust for all other covariates (e.g., age, sex, etc).

Algorithm development
Note that without d j , model (1) would be decoupled as m independent linear regression models, one for each gene. The first step of our algorithm is to solve for d j and express it as a function of β i 's. Note that the optimization problem (6) is convex in (α, d). Therefore, the minimizer of (α, d) is one of its stationary points.
Taking partial derivatives of f α, {β i } m i=1 , d with respect to d j , j = 1, . . . , n, and setting them to zeros, we have The solution to model (1) is not unique because an arbitrary constant can be added to d j 's and subtracted from α i 's, while having the same model fit. To address this issue, we fix d 1 = 0. Therefore wherē y (w) Here the superscript (w) denotes "weighted mean". Calculating the partial derivatives of f α, {β i } m i=1 , d with respect to α i , i = 1, . . . , m, and setting them to zeros, we have wherē x := 1 n n j=1 x j .
From (10) it follows wherē Substituting (16) into (13) yields The sum of (10) and (18) yields Substituting (19) into (6), the problem becomes an 0regularized linear regression problem with {β i } m i=1 being the only variables to be optimized: wherẽ It is easy to see that In the next two sections, we will describe algorithms to solve Problem (20) with type I and type II penalties, respectively.

Fitting the model with type I penalty
Denote δ =β (w) , and let where β i 's are considered as functions of δ. The objective in Problem (20) can be written as f (β) = m i=1 g i β i . Assume that δ is fixed, f can be minimized by minimizing each g i β i separately.
Next we express the minimizing solution of g i β i as a function of δ.
The objective function value at The change in the objective value g i β i from β i = β Therefore, the solution is Now we only need to solve for δ. We havê where the second equality is due to the fact that is a constant independent of δ, and the third equality follows from (31). Problem (33) can be solved exactly using an exhausted grid search for p=1 or 2, and approximately using a general global optimization algorithm (e.g., the optim function in R) for larger p. A more efficient algorithm proposed in [15] can also be used.
After we obtain the estimate of δ, we substitute it into (32) to get the estimate of β i . Algorithm 1 describes the complete model fitting procedure.

Algorithm 1
Algorithm to fit the model with type I penalty Input: Log-transformed gene expression measurements: and design matrix: • Center columns of X to have mean zeros: x j and center y ij to have both zero row and column means: (14), (11) and (17), respectively.
• Select the tuning parameters λ i 's according to "Tuning parameter selection: regression with type I penalty" section. • Estimate the noise variance σ 2 i 's according to "Maximum likelihood estimation of noise variance" section.

Fitting the model with type II penalty
Denote δ =β (w) , and let where β i 's are considered as functions of δ. The objective is fixed, f can be optimized by minimizing each h i β i separately.
Next we find the solution for β i 's as a function of δ by minimizing h i β i . Denotẽ . .
Denote β i into (35) and after some matrix algebraic manipulation, we have When β ip = 0, The minimizing solution of h i β i is β shown in (27), and its p-th coordinate is where the second equality follows fromX = X − x p and the inverse formula for the partitioned matrix ofX TX . The The decrease in the objective value from where the second equality employs the following equality: which is obtained by partitioningX TX into a 2 × 2 block matrix and then substituting the formula for its inverse, and β (ols) ip is defined in Eq. (40). Therefore, the solution is Now we only need to solve for δ p . We havê where the second equality is due to the fact that h i β (ols) i is a constant independent of δ. Substituting (42) into (44) yieldŝ where the β Afterδ p is estimated, the estimate of β ip is obtained by substituting δ p =δ p into (43). Algorithm 2 describes the complete model fitting procedure.
Next, we introduce a simple method for the selection of the tuning parameters in our model, which is based on the property of the solution (32) or (43).

Tuning parameter selection: regression with type I penalty
Substituting (19) into (1) and assuming that δ =β (w) is fixed, we havẽ whereỹ ij + δ Tx j are the normalized data, which we use here as the response variables, and ε ij ∼ N 0, σ 2 i . The condition for β i = 0 in (32) can be rewritten as Under the null hypothesis, β i = 0; the left-hand side of (47) follows a chi-squared distribution with p degrees of freedom, i.e., χ 2 p . This suggests us is the cumulative distribution function of χ 2 p , and q is a pre-specified significance level.

Tuning parameter selection: regression with type II penalty
Letỹ ij + δ Tx j denote the normalized data: where ε ij ∼ N 0, σ 2 i . The condition for β ip = 0 in (43) can be rewritten as Algorithm 2 Algorithm to fit the model with type II penalty Input: Log-transformed gene expression measurements: and design matrix: • Center columns of X to have mean zeros: x j ← x j −x, withx := 1 n n j=1 x j and center y ij to have both zero row and column means: (14), (11) and (17), respectively.
• Select the tuning parameters λ i 's according to "Tuning parameter selection: regression with type II penalty" section. • Estimate the noise variance σ 2 i 's according to "Maximum likelihood estimation of noise variance" section.
whereX − is the submatrix ofX with its last column removed. Solvê For i = 1, . . . , m, estimate β ip : where β (ols) ip is defined in (40) and is the standard error of the estimate β (ols) ip . Under the null hypothesis, β ip = 0; the left-hand side of (49) follows the standard Gaussian distribution. This suggests us choose is the cumulative distribution function of the standard Gaussian distribution, and q is a pre-specified significance level.

Results and discussion
We demonstrate the performance of our proposed method (named rSeqRobust) by comparing it with other existing methods for DE gene detection from RNAseq data: edgeR-robust [1,2], DESeq2 [3], limma-voom [4,5], and ELMSeq (which fits a similar model but with 1 rather than 0 penalty) [11]. We consider the simple regression model (p = 1), in which case Algorithms 1 and 2 coincide. For ELMSeq, the tuning parameter is set as the 5th percentile of m val-  [11]. The tuning parameters λ i is set based on the significant level of q = 0.01.

Simulations on synthetic data
We simulate both log-normally distributed and negativebinomially distributed read counts, with m = 20, 000 genes and n = 7, 20 or 200 samples. The RNA-seq read counts are generated as c ij ∼ e N (log μ ij ,σ 2 i ) under the log-normal (LN) distribution assumption, and as c ij ∼ N B(μ ij , φ i ) [2] under the NB distribution assumption, where μ ij is the mean read counts of gene i in sample j. The generation of μ ij is described in Table 1. The variance of the normal distribution is set as σ 2 i = 0.01, and the dispersion parameter of the NB distribution is set as φ i = 0.25.
After estimating the sample-specific normalization factors d j 's using Algorithm 1, we substituted j 's into model (1) to obtain m decoupled gene-wise linear regression models. For each gene i, we test the null hypothesis that β i = 0, and calculate the p-value. We decide there is a significant linear association between the experimental variable x j and the gene expression y ij if the p-value is less than a predefined threshold (e.g., 0.05). With the p-value for each gene, we rank the genes and vary the p-value threshold from 0 to 1 to determine significant DE genes and calculate the area under the ROC curve (AUC). Table 2 shows the AUCs of the five methods on lognormally distributed data with sample size n = 20. We observe the followings: i) In relatively easy scenarios where a small percent of genes are differentially expressed (i.e., DE%=1% or 10%) or the up-and down-regulated Table 6 The AUCs of edgeR-robust, DESeq2, limma-voom, ELMSeq and rSeqRobust based on log-normally distributed data DE  genes are equal in portions (i.e., Up%=50%), all five methods perform equally well (within one standard error of AUC difference); ii) In challenging scenarios where a large percent of genes are differentially expressed (i.e., DE% ≥30%) and the up-and down-regulated genes are different in portions (i.e., Up%=75% or 100%), rSeqRobust outperforms ELMSeq, which in turn outperforms the rest; iii) In the most challenging scenarios with 70% or 90% DE genes that are all overexpressed (i.e., Up%=100%), only rSeqRobust achieves good results (AUC=0.9638 or 0.8276). Table 3 shows the AUCs of different methods for negative-binomially distributed data. The same observations are obtained as in Table 2: In relatively easy settings with small percent of DE genes or symmetric over-and under-expression pattern, rSeqRobust performs as well as other methods; In challenging settings with large percent of DE genes (DE% ≥10%) and asymmetric over-and under-expression pattern (Up%=75% or 100%), rSeqRobust consistently performs the best, and ELMSeq ranks second in all except extreme cases: (Up%, DE%)=(70%, 100%), (90%, 75%) or (90%, 100%) where most methods suffer from severe performance degradation or complete failure.
In Tables 4 and 5, the sample size is reduced to n = 7. Again, we observe similar patterns: when a small subset of genes are differentially expressed (i.e., DE%=1% or 10%), or the up-and down-regulated DE genes are imbalanced in numbers, rSeqRobust and other methods perform equally well; when most genes are differentially expressed (i.e., DE% = 50% or 70%) in an asymmetric manner (i.e., Up%=75% or 100%), rSeqRobust outperforms all other methods. Note that in the presence of 70% DE genes that are all up-regulated, only rSeqRobust achieves good results (AUC=0.9059 for LN data and AUC=0.8623 for NB data).
In Tables 6 and 7, the sample size is increased to n = 200. As the sample size increases from n = 20 to n = 200, the AUCs of edgeR-robust, DESeq2 and limma-voom increase for easy cases (small percent of DE genes or symmetric over-and under-expression patterns). However, for challenging cases (i.e., DE%=50%, 70% or 90%, Up%=75% or 100%), the AUCs decrease. On the contrary, the AUC of rSeqRobust increases consistently in all cases. The performance gain of rSeqRobust over other methods is more significant for more challenging cases. Note that rSeqRobust performs nearly as well in the most challenging cases (Up%, DE%)=(50%, 100%), (70%, 100%), (90%, 75%) or (90%, 100%) as in easy cases. In contrast, ELMSeq only works for (Up%, DE%)=(50%, 100%) and edgeR-robust, DESeq2, limma-voom completely fail in all these cases. This indicates that rSeqRobust is more robust than ELM-Seq, which in turn is more robust than edgeR-robust, DESeq2 and limma-voom. Table 8 shows the average running times (in seconds) of the five methods on an Intel Core i3 processor with 8GB of memory and a clock frequency of 3.9GHz. We can see that rSeqRobust is slower than limma-voom; however, it scales well for large sample sizes and is much faster than ELMSeq.

Application to a real RNA-seq dataset
We further assess the proposed method on a real RNA-seq dataset from The Cancer Genome Atlas (TCGA) project [17], which contains 20,531 genes from 187 prostate adenocarcinoma patient samples. The dataset was downloaded from the TCGA data portal (https://portal.gdc. cancer.gov). In this experiment, we aim at identifying genes associated with pre-operative prostate-specific antigen (PSA), which is an important biomarker for prostate cancer. The data are pre-processed using the procedures described in [11]. We use the Bonferroni correction and determined DE genes using a p-value threshold of 0.05/m. Figure 1 shows the Venn diagram based on the sets of differentially expressed genes discovered by five methods.
There are twelve genes that are detected by rSe-qRobust and ELMSeq, but not by edgeR, DESeq2 and limma-voom: EPHA5, RNF126P1, BCL11A, RIC3, AJAP1, CDH3, WIT1, PRSS16, CEACAM1, DCHS2, CRHR1 and SRD5A2. For the majority of these twelve genes, there are existing publications reporting their associations with prostate cancer. For instance, EPHA5 is known to be underexpressed in prostate cancer [18]. CEACAM1 is known to suppress prostate cancer cell proliferation when overexpressed [19]. Two of the twelve genes, CRHR1 and SRD5A2, are identified only by rSeqRobust, where SRD5A2 is associated with racial/ethnic disparity in prostate cancer risk [20].

Conclusion & discussion
In this paper, we present a unified statistical model for joint normalization and differential expression detection in RNA-seq. Different from existing methods, we explicitly model sample-specific normalization factors as unknown parameters, so that they can be estimated simultaneously together with detection of differentially expressed genes. Using an 0 -regularized regression approach, our method is robust against large proportion of DE genes and asymmetric DE pattern, and is shown in empirical studies to be more accurate in detecting differential gene expression patterns.
This model generalizes [10] from categorical experimental conditions using an ANOVA-type model to continuous covariates using a regression model. In addition, two hypothesis tests are formulated: i) Is the expression Fig. 1 Venn diagram based on the set of differentially expressed genes identified by edgeR, DESeq2, limma-voom, ELMSeq and rSeqRobust level of a gene associated with any covariates of the regression model? This is the test considered in [10]; ii) Is the expression level of a gene associated with a specific covariate of our interest, when all other variables in the regression model are adjusted for? Although the model is high-dimensional, non-differentiable and non-convex due to the 0 penalty, we manage to develop an efficient algorithms to find their its solution by making use of the optimality conditions of the 0 -regularized regression. It can be shown that for categorical experimental data, the developed algorithm for the first hypothesis test for the slopes in a regression model with p binary covariates reduces to that in [10] for the (p + 1)-group model.