Highly efficient hypothesis testing methods for regression-type tests with correlated observations and heterogeneous variance structure

Background For many practical hypothesis testing (H-T) applications, the data are correlated and/or with heterogeneous variance structure. The regression t-test for weighted linear mixed-effects regression (LMER) is a legitimate choice because it accounts for complex covariance structure; however, high computational costs and occasional convergence issues make it impractical for analyzing high-throughput data. In this paper, we propose computationally efficient parametric and semiparametric tests based on a set of specialized matrix techniques dubbed as the PB-transformation. The PB-transformation has two advantages: 1. The PB-transformed data will have a scalar variance-covariance matrix. 2. The original H-T problem will be reduced to an equivalent one-sample H-T problem. The transformed problem can then be approached by either the one-sample Student’s t-test or Wilcoxon signed rank test. Results In simulation studies, the proposed methods outperform commonly used alternative methods under both normal and double exponential distributions. In particular, the PB-transformed t-test produces notably better results than the weighted LMER test, especially in the high correlation case, using only a small fraction of computational cost (3 versus 933 s). We apply these two methods to a set of RNA-seq gene expression data collected in a breast cancer study. Pathway analyses show that the PB-transformed t-test reveals more biologically relevant findings in relation to breast cancer than the weighted LMER test. Conclusions As fast and numerically stable replacements for the weighted LMER test, the PB-transformed tests are especially suitable for “messy” high-throughput data that include both independent and matched/repeated samples. By using our method, the practitioners no longer have to choose between using partial data (applying paired tests to only the matched samples) or ignoring the correlation in the data (applying two sample tests to data with some correlated samples). Our method is implemented as an R package ‘PBtest’ and is available at https://github.com/yunzhang813/PBtest-R-Package. Electronic supplementary material The online version of this article (10.1186/s12859-019-2783-8) contains supplementary material, which is available to authorized users.

t.test(). However, it is common to observe heterogeneous variabilities between high-throughput samples, which violates the assumption of the Student's t-test. For example, samples processed by a skillful technician usually have less variability than those processed by an inexperienced person. For two-group comparisons, a special case of the heterogeneity of variance, i.e., samples in different groups have different variances, is well studied and commonly referred to as the Behrens-Fisher problem. The best known (approximate) parametric solution for this problem is the Welch's t-test, which adjusts the degrees of freedom (hereinafter, DFs) associated with the t-distribution to compensate for the heteroscedasticity in the data. Unfortunately, the Welch's t-test is not appropriate when the data have even more complicated variance structure. As an example, it is well known that the quality and variation of the RNA-seq sample is largely affected by the total number of reads in the sequencing specimen [3,4]. This quantity is also known as sequencing depth or library size, which may vary widely from sample to sample. Fortunately, such information is available a priori to data analyses. Several weighted methods [5][6][7] are proposed to utilize this information and make reliable statistical inference.
As the technology advances and the unit cost drops, immense amount of data are produced with even more complex variance-covariance structures. In multi-site studies for big data consortium projects, investigators sometimes need to integrate omics-data from different platforms (e.g. microarray or RNA-seq for gene expression) and/or processed in different batches. Although many normalization [8][9][10] and batch-correction methods [11][12][13] can be used to remove spurious bias, the heterogeneity of variance remains to be an issue. Besides, the clustering nature of these data may induce correlation among observations within one center/batch. Correlation may arise due to other reasons such as paired samples. For example, we downloaded a set of data for a comprehensive breast cancer study [14], which contain 226 samples including 153 tumor samples and 73 paired normal samples. Simple choices such as Welch's t-test and paired t-test are not ideal for comparing the gene expression patterns between normal and cancerous samples, because they either ignore the correlations of the paired subjects or waste information contained in the unpaired subjects. To ignore the correlation and use a two-sample test imprudently is harmful because it may increase the type I error rate extensively [15]. On the other hand, a paired test can only be applied to the matched samples, which almost certainly reduces the detection power. In general, data that involves two or more matched samples are called repeated measurements, and it is very common in practice to have some unmatched samples, also known as unbalanced study design.
One of the most versatile tools in statistics, the linear mixed-effects regression (LMER), provides an alternative inferential framework that accounts both unequal variances and certain practical correlation structures. The standard LMER can model the correlation by means of random effects. By adding weights to the model, the weighted LMER is able to capture very complex covariance structures in real applications. Although LMER has many nice theoretical properties, fitting it is computationally intensive. Currently, the best implementation is the R package lme4 [16], which is based on an iterative EM algorithm. For philosophical reasons, lme4 does not provide p-values for the fitted models. The R package lmerTest [17] is the current practical standard to perform regression t-and F-tests for lme4 outputs with appropriate DFs. A fast implementation of LMER is available in the Rfast package, which is based on highly optimized code in C++ [2]; however, this implementation does not allow for weights.
Many classical parametric tests, such as two-sample and paired t-tests, have their corresponding rank-based counterparts, i.e. the Wilcoxon rank-sum test and the Wilcoxon signed rank test. A rank-based solution to the Behrens-Fisher problem can be derived based on the adaptive rank approach [18], but it was not designed for correlated observations. In recent years, researchers also extended rank-based tests to situations where both correlations and weights are presented. [19] derived the Wilcoxon ranksum statistic for correlated ranks, and [20] derived the weighted Mann-Withney U statistic for correlated data. These methods incorporate an interchangeable correlation in the whole dataset, and are less flexible for a combination of correlated and uncorrelated ranks. Lumley and Scott [21] proved the asymptotic properties for a class of weighted ranks under complex sampling, and pointed out that a reference t-distribution is more appropriate than the normal approximation for the Wilcoxon test when the design has low DFs. Their method is implemented in the svyranktest() function in R package survey. But most of the rank-based tests are designed for group comparisons; rank-based approaches for testing associations between two continuous variables with complex covariance structure are underdeveloped.
Based on a linear regression model, we propose two H-T procedures (one parametric and one semiparametric) that utilize a priori information of the variance (weights) and correlation structure of the data. In "Methods" section, we design a linear map, dubbed as the "PB-transformation", that a) transforms the original data with unequal variances and correlation into certain equivalent data that are independent and identically distributed; b) maps the original regression-like H-T problem into an equivalent one-group testing problem. After the PB-transformation, classical parametric and rank-based tests with adjusted DFs are directly applicable. We also provide a moment estimator for the correlation coefficient for repeated measurements, which can be used to obtain an estimated covariance structure if it is not provided a priori. In "Simulations" section, we investigate the performance of the proposed methods using extensive simulations based on normal and double exponential distributions. We show that our methods have tighter control of type I error and more statistical power than a number of competing methods. In "A real data application" section, we apply the PB-transformed t-test to an RNA-seq data for breast cancer. Utilizing the information of the paired samples and sequencing depths, our method selects more cancerspecific genes and fewer falsely significant genes (i.e. genes specific to other diseases) than the major competing method based on weighted LMER.
Lastly, computational efficiency is an important assessment of modern statistical methods. Depending on the number of hypotheses to be tested, our method can perform about 200 to 300 times faster than the weighted LMER approach in simulation studies and real data analyses. This efficiency makes our methods especially suitable for fast feature selection in high-throughput data analysis. We implement our methods in an R package called 'PBtest' , which is available at https://github.com/ yunzhang813/PBtest-R-Package.

Model framework
For clarity, we first present our main methodology development for a univariate regression problem. We will extend it to multiple regression problems in "Extension to multiple regressions" section.
Consider the following regression-type H-T problem: where μ, β ∈ R, y, x, , 1 = (1, · · · , 1) ∈ R n and ∼ N (0, ); Here, y is the response variable, x is the covariate, and is the error term that follows an n-dimensional multivariate normal distribution N with mean zero and a general variance-covariance matrix . By considering a random variable Y in the n-dimensional space, the above problem can also be stated as In this model, μ is the intercept or grand mean that is a nuisance parameter, and β is the parameter of interest that quantifies the effect size. We express the variancecovariance matrix of in the form where σ 2 is a nonzero scalar that quantifies the magnitude of the covariance structure, and S is a symmetric, positivedefinite matrix that captures the shape of the covariance structure. Additional constraints are needed to determine σ 2 and S; here, we choose a special form that can subsequently simplify our mathematical derivations. For any given , define From the above definition, we have the following nice property Hereinafter, we refer to S the standardized structure matrix satisfying Eq. 5.

The proposed method
As a special case of Model (3), if S is proportional to I, the identity matrix, it is well-known that regression t-test is a valid solution to this H-T problem. If S = I, e.g. the observed data are correlated and/or have heterogeneous variance structure, the assumptions of the standard t-test are violated. In this paper, we propose a linear transformation, namely PB : Y →Ỹ, which transforms the original data to a new set of data that are independent and identically distributed. Furthermore, we prove that the transformed H-T problem related to the new data is equivalent to the original problem, so that we can approach the original hypotheses using standard parametric (or later rank-based) tests with the new data.
To shed more lights on the proposed method, we first provide a graphical illustration in Fig. 1. The proposed procedure consists of three steps.
1 Estimateμ(Y) (i.e. the weighted mean of the original data), and subtractμ from all data. This process is an oblique (i.e. non-orthogonal) projection from R n to an (n − 1)-dimensional subspace of R n . The intermediate data from this step is Y (1) (i.e. the centered data). It's clear that EY (1) is the origin of the reduced space if and only if H 0 is true. 2 Use the eigen-decomposition of the covariance matrix of Y (1) to reshape its "elliptical" distribution to a "spherical" distribution. The intermediate data from this step is Y (2) . Step 1: Estimateμ(Y) (i.e. the weighted mean of the original data), and subtractμ from all data. This process is an oblique (i.e. non-orthogonal) projection from R n to an (n − 1)-dimensional subspace of R n . The intermediate data from this step is Y (1) , also called the centered data. If H 0 is true, Y (1) centers at the origin of the reduce space; otherwise, the data cloud Y (1) deviates from the origin.
Step 2: Use eigen-decomposition to reshape the "elliptical" distribution to an "spherical" distribution. The intermediate data from this step is Y (2) .
Step 3: Use QR-decomposition to find a unique rotation that transforms the original H-T problem to an equivalent problem. The equivalent problem tests for a constant deviation along the unit vector in the reduced space, thus it can be approached by existing parametric and rank-based methods. The final data from this step isỸ 3 Use the QR-decomposition technique to find a unique rotation that transforms the original H-T problem to an equivalent problem of testing for a constant deviation along the unit vector. The equivalent data generated from this step isỸ, and the H-T problem associated withỸ can be approached by existing parametric and rank-based methods.
In the proposed PB-transformation, B-map performs both transformations in Step 1 and 2; P-map from Step 3 is desi gned to improve the power of the proposed semiparametric test to be described in "A semiparametric generalization" section.

Centering data
Using weighted least squares, the mean estimation based on the original data isμ(Y) = 1 S −1 Y (for details please see Additional file 1: Section S1.1). We subtractμ from all data points and define the centered data as where J = 1 · 1 (i.e. a matrix of all 1's). With some mathematical derivations (see Additional file 1: Section S1.1), we have

The B-map
Now, we focus on S−J, which is the structure matrix of the centered data. Let T T denote the eigen-decomposition of S − J. Since the data are centered, there are only n − 1 nonzero eigenvalues. We express the decomposition as follows where T n−1 ∈ M n×(n−1) isasemi-orthogonalmatrixcontaining the first n − 1 eigenvectors and n−1 ∈ M (n−1)×(n−1) is a diagonal matrix of nonzero eigenvalues. Based on Eq. 6, we define (see Additional file 1: Section S1.2) so that Y (2) := BY ∈ R n−1 have the following mean and covariance We call the linear transformation represented by matrix B the "B-map". So far, we have centered the response variable, and standardized the general structure matrix S into the identity matrix I. However, the covariate and the alternative hypothesis in the original problem are also transformed by the B-map. For normally distributed Y, the transformed H-T problem in Eq. 7 is approachable by the regression t-test; however, there's no appropriate rankbased counterpart. In order to conduct a rank-based test for Y with broader types of distribution, we propose the next transformation.

The P-map
From Eq. 7, define the transformed covariate We aim to find an orthogonal transformation that aligns z to 1 n−1 in the reduced space. We construct such a transformation through the QR decomposition of the following object where A ∈ M (n−1)×2 is a column-wise concatenation of vector z and the target vector 1 n−1 , Q ∈ M (n−1)×2 is a semi-orthogonal matrix, and R ∈ M 2×2 is an upper triangular matrix. We also define the following rotation matrix Geometrically speaking, ξ = cos θ, where θ is the angle between z and 1 n−1 .
With the above preparations, we have the following result.
Q is the unique orthogonal transformation that satisfies the following properties: Proof See Additional file 1: Section 1.3.
We call the linear transformation P defined by Theorem 1 the "P-map". Equation 9 ensures that this map is an orthogonal transformation. Equation 10 shows that the vector z is mapped to 1 n−1 scaled by a factor ζ . Equation 11 is an invariant property in the linear subspace L ⊥ z , which is the orthogonal complement of the linear subspace spanned by 1 n−1 and z, i.e. L z = span (1 n−1 , z). This property defines a unique minimum map that only transforms the components of data in L z and leaves the components in L ⊥ z invariant. A similar idea of constructing rotation matrices has been used in [22].
With both B and P, we define the final transformed data asỸ := PY (2) = PBY, which has the following joint distributioñ The normality assumption implies that eachỸ i follows an i.i.d. normal distribution, for i = 1, · · · , n − 1. The location parameter of the common marginal distribution is to be tested with unknown σ 2 . Therefore, we can approach this equivalent H-T problem with the classical one-sample t-test and Wilcoxon signed rank test (more in "A semiparametric generalization" section).

Correlation estimation for repeated measurements
If is unknown, we can decompose in the following way where W is a diagonal weight matrix and Cor is the corresponding correlation matrix. By definition, the weights are inversely proportional to the variance of the observations. In many real world applications including RNA-seq analysis, those weights can be assigned a priori based on the quality of samples; but the correlation matrix Cor needs to be estimated from the data. In this section, we provide a moment-based estimator of Cor for a class of correlation structure that is commonly used for repeated measurements. This estimator does not require computationally intensive iterative algorithms. Let Y be a collection of repeated measures from L subjects such that the observations from different subjects are independent. With an appropriate data rearrangement, the correlation matrix of Y can be written as a block-diagonal matrix We assume that the magnitude of correlation is the same across all blocks, and denote it by ρ. Each block can be expressed as Cor l (ρ) = (1−ρ)I n l ×n l +ρJ n l ×n l , for l = 1, · · · , L, where n l is the size of the lth block and n = L l=1 n l . We estimate the correlation based on the weighted regression residualsˆ defined by Eq. (S3) in Additional file 1: Section S2.1. Define two forms of residual sum of squares whereˆ l is the corresponding weighted residuals for the lth block. With these notations, we have the following Proposition.
An estimator of ρ based on the first moments of SS 1 and SS 2 iŝ Moreover, ifˆ ∼ N (0, ) and n 1 = · · · = n L = n/L (i.e. balanced design), the above estimator coincides with the maximum likelihood estimator of ρ, which has the form Proof See Additional file 1: Section S2.1.
Standard correlation estimates are known to have downward bias [23], which can be corrected by the Olkin and Pratt's method [24]. With this correction, our final correlation estimator iŝ

Kenward-roger approximation to the degrees of freedom
The degree of freedom (DF) can have nontrivial impact on hypothesis testing when sample size is relatively small. Intuitively, a correlated observation carries "less information" than that of an independent observation. In such case, the effective DF is smaller than the apparent sample size. Simple examples include the two-sample t-test and the paired t-test. Suppose there are n observations in each group, the former test has DF = 2n − 2 for i.i.d. observations, and the latter only has DF = n − 1 because the observations are perfectly paired. These trivial examples indicate that we need to adjust the DF according to the correlation structure in our testing procedures. We adopt the degrees of freedom approximation proposed by [25] (K-R approximation henceforth) for the proposed tests. The K-R approximation is a fast momentmatching method, which is efficiently implemented in R package pbkrtest [26]. In broad terms, we use the DF approximation as a tool to adjust the effective sample size when partially paired data are observed.

Alternative approach using mixed-effects model
As we mentioned in "Background" section, the H-T problem stated in Model (3) for repeated measurements can also be approached by the linear mixed-effects regression (LMER) model. Suppose the ith observation is from the lth subject, we may fit the data with a random intercept model such that The LMER model is typically fitted by a likelihood approach based on the EM algorithm. Weights can be incorporated in the likelihood function. The lmer() function in R package lme4 [16] provides a reference implementation for fitting the LMER model. The algorithm is an iterative procedure until convergence. Due to relatively high computational cost, the mixed-effects model has limited application in high-throughput data.
The R package lmerTest [17] performs hypothesis tests for lmer() outputs. By default, it adjusts the DF using the Satterthwaite's approximation [27], and can optionally use the K-R approximation.

A semiparametric generalization
In the above sections, we develop the PB-transformed t-test using linear algebra techniques. These techniques can be applied to non-normal distributions to transform their mean vectors and covariance matrices as well. With the following proposition, we may extend the proposed method to an appropriate semiparametric distribution family. By considering the uncorrelated observations with equal variance as a second order approximation of the data that we are approaching, we can apply a rank-based test on the transformed data to test the original hypotheses. We call this procedure the PB-transformed Wilcoxon test.

Proposition 2
LetY := Y 1 , . . . ,Y n−1 be a collection of i.i.d. random variables with a common symmetric density function g(y), g(−y) = g(y). Assume that EY 1 = 0, var(Y 1 ) = σ 2 . Let Y * be a random number that is independent ofY and has zero mean and variance σ 2 . For every symmetric semi-definite S ∈ M n×n , x ∈ R n and μ, β ∈ R, there exists a linear transformation D : R n−1 → R n and constants u, v, such that is an n-dimensional random vector with Furthermore, if we apply the PB-transformation to Y, the result is a sequence of (n − 1) equal variance and uncorrelated random variables with zero mean if and only if β = 0.
The essence of this Proposition is that, starting with an i.i.d. sequence of random variables with a symmetric common p.d.f., we can use linear transformations to generate a family of distributions that is expressive enough to include a non-normal distribution with an arbitrary covariance matrix and a mean vector specified by the effect to be tested. This distribution family is semiparametric because: a) the "shape" of the density function, g(y), has infinite degrees of freedom; b) the "transformation" (D, u, and v) has only finite parameters.
As mentioned before, applying both the B-and P-maps enables us to use the Wilcoxon signed rank test for the hypotheses with this semiparametric distribution family. This approach has better power than the test with only the B-map as shown in "Simulations" section . Once the PBtransformed data are obtained, we calculate the Wilcoxon signed rank statistic and follow the testing approach in [21], which is to approximate the asymptotic distribution of the test statistic by a t-distribution with an adjusted DF. Note that Wilcoxon signed rank test is only valid when the underlying distribution is symmetric; therefore, the symmetry assumption in Proposition 2 is necessary. In summary, this PB-transformed Wilcoxon test provides an approximate test (up to the second order moment) for data that follow a flexible semiparametric distributional model.

Extension to multiple regressions
In this section, we present an extension of the proposed methods for the following multiple regression Here the error term is assumed to have zero mean but does not need to have scalar covariance matrix. For example, can be the summation of random effects and measurement errors in a typical LMER model with a form specified in Eq. 4.
To test the significance of β k , k = 1, . . . , p, we need to specify two regression models, the null and alternative models. Here the alternative model is just the full Model (16), and the null model is a regression model for which the covariate matrix is X −k , which is constructed by removing the kth covariate (X k ) from X Compared with the original univariate problem, we see that the nuisance covariates in the multiple regression case are X −k β −k instead of 1μ in Eq. 1. Consequently, we need to replace the centering step by regressing out the linear effects of X −k The new B-transformation is defined as the eigendecomposition of cov (E) = σ 2 S − X −k X −k . The Ptransformation is derived the same as before, but with the new B matrix.

Simulations
We design two simulation scenarios for this study: SIM1 for completely paired group comparison, and SIM2 for regression-type test with a continuous covariate. For both scenarios we consider three underlying distributions (normal, double exponential, and logistic) and four correlation levels (ρ = 0.2, ρ = 0.4, ρ = 0.6, and ρ = 0.8). We compare the parametric and rank-based PB-transformed test with oracle and estimated correlation to an incomplete survey of alternative methods. Each scenario was repeated 20 times and the results of ρ = 0.2 and 0.8 for normal and double exponential distributions are summarized in Figs. 2 and 3, and Tables 1 and 2. See Additional file 1, Section S3 for more details about the simulation design, additional results of ρ = 0.4 and 0.6, and results for logistic distribution. Figures 2 and 3 are ROC curves for SIM1 and SIM2, respectively. In all simulations, the proposed PBtransformed tests outperform the competing methods.
The PB-transformed t-test has almost identical performance with oracle or estimated ρ. Using the estimated ρ slightly lowers the ROC curve of the PB-transformed Wilcoxon test compared with the oracle curve, but it still has a large advantage over other tests. Within the parametric framework, the weighted LMER has the best performance among the competing methods. It achieves similar performance as our proposed parametric test when the correlation coefficient is small; however, its performance deteriorates when the correlation is large. Judging from the ROC curves, among the competing methods, the svyranktest() is the best rank-based test for the A B C D group comparison problem, primarily because it is capable of incorporating the correlation information. However, it fails to control the type-I error, as shown in Table 1. Tables 1 and 2 summarize the type-I error rate and power at the 5% significance level for SIM1 and SIM2, respectively. Overall, the PB-transformed tests achieve the highest power in all simulations. In most cases, the proposed tests tend to be conservative in the control of type-I error; and replacing the oracle ρ by the estimated ρ does not have significant impact on the performance of PB-transformed tests. The only caveat is the rank-based test for the regression-like problem. Currently, there's no appropriate method designed for this type of problem. When the oracle correlation coefficient is provided to the PB-transformed Wilcoxon test, it has tight control of type I error. With uncertainty in the estimated correlation coefficient, our PB-transformed Wilcoxon test may suffer from slightly inflated type I errors; but it is still more conservative than its competitors. Of note, other solutions, such as the naive t-test and rank-based tests, may have little or no power for correlated data, though they may not have the lowest ROC curve.

Computational cost and degrees of freedom
We record the system time for testing 2000 simulated hypotheses using our method and lmer(), since they are the most appropriate methods for the simulated data with the best statistical performance. Our method takes less than 0.3 s with given , and less than 0.9 s with the estimation step; lmer() takes 182 s. We use a MacBook Pro equipped with 2.3 GHz Intel Core i7 processor and 8GB RAM (R platform: x86_64-darwin15.6.0). Of note, lmer() may fail to converge occasionally, e.g. 0 -25 failures (out of 2,000) in each repetition of our simulations. We resort to a try/catch structure in the R script to prevent these convergence issues from terminating the main loop.
We also check the degrees of freedom in all applicable tests. In this section, we report the DFs used/adjusted in SIM1, i.e. the completely paired group comparison.

A real data application
We download a set of RNA-seq gene expression data from The Cancer Genome Atlas (TCGA) [14] (see Additional file 1: Section S4). The data are sequenced on the Illumina GA platform with tissues collected from breast cancer subjects. In particular, we select 28 samples from the tissue source site "BH", which are controlled for white female subjects with the HER2-positive (HER2+) [28] biomarkers. After data preprocessing based on nonspecific filtering (see Additional file 1: Section S4.1), a total number of 11,453 genes are kept for subsequent analyses. Among these data are 10 pairs of matched tumor and normal samples, 6 unmatched tumor samples, and 2 unmatched normal samples. Using Eq. 13, the estimated correlation between matched samples across all genes iŝ ρ = 0.10.
The sequencing depths of the selected samples range from 23.80 million reads to 76.08 million reads. As mentioned before, the more reads are sequenced, the better is the quality of RNA-seq data [4]; thus it is reasonable to weigh samples by their sequencing depths. Since this quantity is typically measured in million reads, we set the weights w i = sequencing depth of the ith sample × 10 −6 , (18) for i = 1, · · · , 28.
With the above correlation estimate and weights, we obtained the covariance structure using Eq. 12. For properly preprocessed sequencing data, a proximity of normality can be warranted [29]. We applied the PB-transformed t-test and the weighted LMER on the data.
Based on the simulations, we expect that if correlation is small, the PB-transformed t-test should have tighter control of false positives than alternative methods. At 5% false discovery rate (FDR) level combined with a fold-change (FC) criterion (FC < 0.5 or FC > 2), the PB-transformed t-test selected 3,340 DEGs and the weighted LMER selected 3,485 DEGs (for biological insights of the DEG lists, see Additional file 1: Section S4.4).
To make the comparison between these two methods more fair and meaningful, we focus on studying the biological annotations of the top 2,000 genes from each DEG list. Specifically, we apply the gene set analysis tool DAVID [30] to the 147 genes that uniquely belong to one list. Both Gene Ontology (GO) biological processes [31] and KEGG pathways [32] are used for functional annotations. Terms identified based on the 147 unique genes in each DEG list are recorded in Additional file 1: Table S6. We further pin down two gene lists, which consist of genes that participate in more than five annotation terms in the above table: there are 11 such genes (PIK3R2, AKT3,  MAPK13, PDGFRA, ADCY3, SHC2, CXCL12, CXCR4,  GAB2, GAS6, and MYL9) for the PB-transformed t-test, and six (COX6B1, HSPA5, COX4I2, COX5A, UQCR10, and ERN1) for the weighted LMER. Expression level of these genes are plotted in Fig. 4. These DEGs are biologically important because they are involved in multiple biological pathways/ontology terms.
Those 11 genes uniquely identified by the PBtransformed t-test are known to be involved in cell survival, proliferation and migration. The CXCR4-CXCL12 chemokine signaling pathway is one of the deregulated signaling pathway uniquely identified by PB-transformed t-test in HER2+ breast cancer cells. This pathway is known to play a crucial role in promoting breast cancer metastasis and has been reported to be associated with poor prognosis [33,34]. Compared with the state-of-the-art method (weighted LMER), the PBtransformed t-test identifies more genes whose protein products can be targeted by pharmaceutical inhibitors. CXCR4 inhibitors have already demonstrated promising anti-tumor activities against breast [35,36], prostrate [37] and lung [38] cancers. Additional downstream signaling molecules identified by our analysis to be significantly associated with HER2+ breast tumor such as PI3K, p38, adaptor molecule GAB2 and SHC2 can also be potential therapeutic targets for selectively eliminating cancer cells. Please refer to Additional file 1: Section S4.5 for full list of functional annotation terms.

Discussion
In this paper, we present a data transformation technique that can be used in conjunction with both the Student's t-type test and rank-based test. In the simulation studies, our proposed tests outperform the classical tests (e.g. twosample/regreesion t-test and Wilcoxon rank-sum test) by a large margin. In a sense, this superiority is expected, because the classical methods do not consider the correlation nor heteroscedasticity of the data.
In our opinion, the most practical comparison in this study is the one between the PB-transformed t-test and the weighted LMER. The fact that the PB-transformed t-test outperforms the weighted LMER, and this advantage is more pronounced for data with higher correlation (see e.g., Figs. 2 and 3), is the highlight of this study, which may have profound implications for applied statistical practice.  Table S6. These genes are not only differentially expressed, but also biologically meaningful