{'fraction_non_alphanumeric': 0.03301228364042434, 'fraction_numerical': 0.008724176437744277, 'mean_word_length': 4.634683444750295, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stratified Test Alleviates Batch Effects in Single-Cell Data

Analyzing single-cell sequencing data across batches is challenging. We find that the Van Elteren test, a stratified version of Wilcoxon rank-sum test, elegantly mitigates the problem. We also modified the common language effect size to supplement this test, further improving its utility. On both simulated and real patient data we show the ability of Van Elteren test to control for false positives and false negatives. The effect size also estimates the differences between cell types more accurately.


Introduction
Large-scale studies such as the Human Cell Atlas [15] involve hundreds of laboratories, thousands of patients, and millions of cells, bringing about opportunities and challenges in analyses. When comparing cell types or groups, discrepancies across experiments and differences among participants lead to omissions and false discoveries in differentially expressed genes. Even the trend (upregulated or downregulated) can be reversed in a phenomenon called Simpson's paradox [1]. Although multiple methods have been proposed to tackle such batch effects, no such option for the widely used Wilcoxon rank-sum test [11,18] has been applied to single-cell studies, to the best of our knowledge. Here, we show that the stratified rank-sum test (known as Van Elteren test [17]) and our modified common language effect size may benefit single-cell studies.
We briefly review and conceptually compare related works on correcting batch effect in Sect. 1.1. Then, in Sect. 2, we revisit Wilcoxon rank-sum test, and introduce the Van Elteren test supplemented by our direct extension of the common language effect size [5,12]. Simulation studies and applications to real data in Sect. 3 show that the test controls for the batch effects and leads to more accurate biological discovers, compared with Wilcoxon rank-sum test. More discussions and explanations, are shown in Sect. 4.

Related Works
Mainstream methods to mitigate batch effect fall into two categories, batch correction methods and batch-aware statistical tests. The former includes methods reducing batch effect in the data to facilitate downstream analysis, while the latter includes analyses that control for the batch effect.
Batch Correction Methods. Batch correction methods eliminate the discrepancy among batches to create an integrated dataset. The most conspicuous manifestation of batch effect is splitting one cell type into multiple clusters. To solve this problem, many methods match and combine clusters across samples based on similarities. A commonly adopted one, Mutual nearest neighbor (MNN) [3], uses similar cells across datasets as anchors, and based on them correct the gene expression of other cells. Scanorama [4] and Seurat [16] are both based on the MNN methodology. Another method, Harmony [7], iteratively corrects the data by clustering the cells and moving neighboring clusters toward each other. These methods typically produce a unified data matrix, which can be conveniently used in visualization and downstream analysis. However, these empirical corrections usually lack negative control and raise uncertainty in the discovery [14].
Batch-Aware Statistical Tests. Instead of manipulating the data directly, statistical methods may handle batch effect by considering it as a covariate in the model. This is possible in Student's t-test, Poisson test, negative binomial test, etc. Notably, all these tests are parametric, meaning that a distribution must be given in advance. However, the debate of the true distribution of singlecell gene expression has never ceased, which is a reason why the nonparametric Wilcoxon rank-sum test is widely used. To allow modeling covariates, one may use a generalized version of rank-sum test, the proportional odds model [2]. However, modeling batches by using a covariate also makes unnecessary assumptions upon them. Stratification, which only combines statistics from batches, is the "as simple as possible, but no simpler" way to handle batches. The Van Elteren test we use, is the stratified version of Wilcoxon rank-sum test.
It is worth noting that methods like scVI [9] have combined statistical modeling with batch effect correction. However, the effect of batches is modeled by a black-box neural network, making the interpretation elusive.

Wilcoxon Rank-Sum Test
We briefly revisit the Wilcoxon rank-sum test (also known as Mann-Whitney U test) [11,18]. The test statistics U is defined as where are the two samples to be compared (e.g., two cell types in one experiment), with sample sizes n 1 and n 2 , respectively. Function "1" takes value 1 when its condition holds true, and 0 otherwise. When n 1 and n 2 are both at least 10, which is common in single-cell studies, the distribution of U approximately follows a normal distribution N (μ, σ 2 ) where and in which t i stands for ties (corresponding to the second term in Eq. 1).

Van Elteren Test
The Van Elteren test [17] is the stratified version of Wilcoxon rank-sum test. For example, if there are m patients, they maybe treated as strata. In that case, a U statistic may be obtained from each patient g ∈ {1, · · · , m}, denoted as where w g is a weight for each sample to be discussed later. When m = 1, the formula degenerates to V = (U g − μ g ) 2 /σ 2 g ∼ χ 2 1 , which is consistent with the rank-sum test.

Weights.
As discussed by Van Elteren [17], the weights w g can be assigned in different ways. It should be noted that the U g for a batch g ranges from 0 to n g1 n g2 , the product of two sample sizes in the batch. Should the weights all be equal, a patient with more cells available will dominate the test results. It is proven in [17] that weight eliminates such effect, and a test utilizing such weight is thus named as "designfree test". However, given that a batch with more instances available (e.g., a patient with more cells sequenced) may be more convincing, another weight is introduced, which gives more power to larger batches. It also effectively assigns larger weights to batches whose samples are more balanced, when the batch sizes are the same. It is shown in [17] that this choice yields largest statistical power against randomized alternatives, and is thus named as "locally-best test". The comparison of two weights are shown in Fig. 1. Effect Sizes. For Wilcoxon rank-sum test, a simple definition of effect size is which is centered at 50%, meaning the probability P (a > b) when a and b are randomly drawn from sample A and sample B, respectively. An effect size greater than 50% generally means that A is higher, and vice versa. It may be easily extended for Van Elteren test by taking average using desired weights. For the design-free test, the effect size is as all batches are treated equally regardless of the sample sizes. It may be interpreted as the probability of P (a g > b g ) for a g and b g randomly drawn from A g and B g , after randomly choosing a batch g. For the locally-best test, the effect size becomes which changes the probability of choosing a group g to be in proportion to ng1ng2 ng1+ng2+1 , giving higher weights to batches with higher and balanced sample sizes (Fig. 1). Generally, any w g may be used to define f , as in the two previous options being its special cases.

Results
We implemented the Van Elteren test with the effect size in R, available at our GitHub repository (https://github.com/KChen-lab/stratified-tests-for-seurat), based on Seurat 3.0 by utilizing its differential expression analysis part (but irrelevant to the data integration) [16]. When there are two groups of cells, A and B, denoted as type, and patient identity, denoted as batch, the Van Elteren Test can be called as follows. The genre may be set to either locally-best or design-free, as introduced in Sect. 2, based on which p-values and effect sizes are calculated. Typical results are shown in Table 2. An effect size of larger than 0.5 indicates a higher expression in cell type A, and vice versa. The avg logFC, average logarithmic fold changes, are calculated automatically by Seurat, where a positive value indicates a higher expression. It may show different trends compared with the effect sizes. Generally, the effect sizes are more indicative after controlling for the batch effect.

Simulation Study
We simulated datasets to illustrate the key utilities of Van Elteren test. The parameters are specified in Table 1. Poisson distribution is used to model sequencing depth. Visualization is available in Fig. 2 for illustration. We assume that the library size of each sample is equalized by other genes beyond the simulated ones.  Suppressing False Negatives. Batch effect may introduce false negatives, where a significantly differentially expressed gene is overshadowed. For gene 1, on which cell type B always have higher expression on both patients, the Wilcoxon rank-sum test did not pass the threshold of 0.05 to reject the null hypothesis, while Van Elteren test yields a significant p-value. The effect size, smaller than 0.5, also correctly suggests that the expression of gene 1 in cell type B is higher than that in cell type A, compared with the average logarithmic fold change, which wrongfully indicates otherwise.  Suppressing False Positives. False discoveries are also possible outcome of batch effect. As is shown in gene 3, the distribution of both cell types are exactly the same in each patient. Nevertheless, Wilcoxon test yields a very significant p-value. The average logarithmic fold change also has a large magnitude. Van Elteren test returns a p-value of 0.7831, together with a effect size close to 0.5, suggesting that the difference is neither significant nor large.

Consistency.
As a negative control, when the three issues above are not present, p-value from Van Elteren test is consistent with Wilcoxon rank-sum test, as is shown by gene 4. The effect size and the log fold change also both show that the cell type B has higher expression in gene 4.

Retina Data
We have tested the Van Elteren test on real retina single-cell gene expression data gathered from three patients [8]. Two regions, macula (i.e., the center area) and peripheral, are labeled in the data. We question which genes differentially express for the same cell type between two regions. We run Wilcoxon rank-sum test and Van Elteren test on 2,295 rod cells and 203 cone cells. We compare the results in Fig. 3, where genes with large differences in p-values between two tests are labeled with gene names.

Rod Cells.
The results of two tests are largely comparable, showing a diagonal pattern. Meanwhile, some exceptions are present (see Table 3  and effect sizes), among which we observed that p-values for gene GNGT1 and SYNE2 change the most. For GNGT1, the reversed conclusion effect is also observed, as the Van Elteren f effect size suggests that the peripheral region has a higher expression, while the log fold change indicates otherwise. We further inspected the distributions to validate and interpret the differences. In Fig. 4A, the left panel does show generally higher expression of GNGT1 in each individual patient, while the aggregated distribution on the right shows a reversed effect, which is an instance of the aforementioned Simpson's paradox. For SYNE2 (Fig. 4B), conspicuous discrepancy among batches is also shown, which leads to a less precise rank-sum test result. Indeed, these two genes were found playing roles in macular degeneration diseases [6,10].
Cone Cells. For cone cells, we observed similar results. Some genes show changes while most genes are consistent across the tests. The p-values and effect sizes are shown in Table 4.
For gene PCBP4, Van Elteren test shows more significant p-value, and an effect size indicating smaller expression in macula, which is different from the log fold change. Decrease in PCBP4 has also been linked with age-related macular degeneration [13]. Figure 5 shows that batch effect in distribution of PCBP4 may have misled the rank-sum test and the logarithmic fold change.

Discussion
The results have clearly shown that Van Elteren test benefits biological studies in precisely identifying differentially expressed genes. Although the results we show do not include multiple comparison correction, Seurat 3.0 will automatically give corrected p-value based on the raw p-value using Bonferroni correction. Generally, any correction based on p-values will also apply. The result also indicates that stratified test is a neat way to handle batch effect. Although covariate has the ability to control for explanatory variables, it is generally more suitable for continuous variables. It also casts more assumptions when modeling covariate. Stratified test, on the other hand, does not infer the influence of the discrete batches. Rather, it directly aggregates the statistical power of multiple samples.
Admittedly, for the rod cells, although changes in p-values are observed, the significance threshold was well passed by both. However, it should be noted that the retina data are collected from relatively healthy tissues and are considered clean, while Van Elteren test is expected to make a more meaningful difference on noisy pathological and tumor data. In addition, rod cell is the most populous cell in retina. For rare cell types that take smaller proportions, like the cone cells, the difference Van Elteren test makes can be crucial.
The caveat of stratified test is that for it to work the strata shall not overlap with the variable of interest. For instance, it may not find the difference, meanwhile also control for the batch effect, between two patients. Nonetheless, neither is covariate applicable to such cases. As the batch effect and biological effect are convolved, more prior knowledge is generally needed to distinguish them.

Summary
We have adopted Van Elteren test, an underappreciated statistical test, and our weighted common language effect size to single-cell sequencing data. When batch effect is severe, the test control for false positives. Otherwise, it is consistent with Wilcoxon rank-sum test. The modified common language effect size also faithfully depicts the trends. This work may increase the precision of differential expression analysis to help identify genes of interests.