scAEQN: A Batch Correction Joint Dimension Reduction Method on scRNA- seq Data

doi:10.21203/rs.3.rs-1759258/v1

Download PDF

Research Article

scAEQN: A Batch Correction Joint Dimension Reduction Method on scRNA- seq Data

https://doi.org/10.21203/rs.3.rs-1759258/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

With the continuous progress of single-cell RNA sequencing technology, the scale of single-cell RNA sequence(scRNA-seq) dataset is increasing, which produces batch effects and affects the downstream analysis of scRNA-seq datasets, such as cluster analysis, gene differential expression analysis, and so on.

Results

Based on the background, we propose a novel batch correction joint dimensionality reduction method named scAEQN, which utilizes the QuantNorm to obtain coefficient matrix and constructs autoencoder to approximate the coefficient matrix, and eventually obtains the low-dimensional representation and the reconstruct of the data. We compared different batch correction methods on six real single-cell RNA sequencing datasets, and scAEQN outperforms benchmark batch-effect correction methods in downstream analyses.

Conclusion

scAEQN can remove batch effects effectively and retain cell type cluster pattern greatly, which provides strong support for downstream analysis. scAEQN improves the performance of clustering and selects more representative differential expression genes in the analysis of gene differential expression.

batch correction

quantile normalization

autoencoder

dimensionality reduction

scRNA-seq data

Single-cell RNA sequencing (scRNA-seq) technology is the most widely used high-throughput sequencing technology at present. With the continuous progress of sequencing technology, the scale of cells is increasing. Plenty of samples from different batches will be sequenced together. These batches usually include different experimental platforms, library preparation, dates, experiment operators, and so on, which often produce technical deviation and batch effect, those differences will be confused with real biological differences, and it will affect downstream analysis [1]. For example, it may generate a cluster pattern that samples of the same batch with batches in cluster analysis, which will lead to wrong biological discovery in cell type identification, or in gene differential expression analysis, it will affect the identification of wrong differential genes, resulting in biological wrong discovery likewise. Therefore, it is very important to remove batch effects before downstream analysis.

At present, there are many methods to remove batch effects. Some tools developed for microarray data batch correction such as Combat [2] and Limma [3] have been employed on single-cell RNA-seq data. Combat-seq [4] is proposed based on the Combat method, which can directly remove batch effects on RNA sequence data instead of count matrix. Researchers also utilize the control probe technology in microarray [5] and spike-in gene technology [6] in RNA-seq to find and correct unknown batch effects in the microarray. However, these methods need to understand the knowledge of each kind of biological population, but it is almost impossible due to the high heterogeneity of scRNA-seq data [7].

Due to the above limitations of these methods, many strategies are motivated to remove batch effects for scRNA-seq data. For example, MNN method [8], assuming that there is an orthogonal relationship between batch information and biological information, mutual nearest neighbor pairs are performed on the data between batches to correct the batch of datasets. BBKNN [9], a graph-based batch correct algorithm, modifies the neighborhood construction step of k-nearest neighbors-based MNN method, like the Scanorama method [10]. Harmony [11] first employs PCA for dimensionality reduction, then iteratively removes batch effects present in the PCA space. QuantNorm [12] is proposed to correct the sample distance matrix through quantile normalization. Besides, based on deep learning technology, like the MMD-ResNet method [13], which utilizes the residual neural network to take the maximum mean discrepancy of two data sets as its loss function to remove batch effects, the data after dimensionality reduction is obtained for cluster analysis. The scVI method [14] with a variational autoencoder model to correct batch, reduce dimensionality and reconstruct data.

Although a large number of batch correction methods have been available for biological workers, we find that there is still room for improvement. For example, most batch correction methods such as MNN [8], Limma [3], Combat [2] and so on, obtained the corrected count expression matrix that retained the scale of raw data, it will not get satisfactory results in clustering by utilizing corrected matrices. Therefore, it is also necessary to perform dimensionality reduction before clustering. In contrast, there are some batch correction methods to obtain the distance matrix and low-dimensional matrix, such as BBKNN [9], QuantNorm [12], MMD-ResNet [13], scVI [14], Scanorama [10] and Harmony [11]. But they cannot support differential expression genes (DEGs) in downstream analysis. As we all know, differential expression gene analysis is a vital part of scRNA-seq data analysis. It is very useful to find the significant DEGs between distinct subpopulations or groups of cells. Although QuantNorm was improved and scBatch [15] was proposed to approximate the dissimilarity matrix obtained by the QuantNorm with a linear weighting, we found that scBatch will take about several hours on the data when the number of cells exceeds 1000, which is not the result what we expect.

Based on the current challenges of batch correction, in this paper, we propose a novel batch correction joint dimensionality reduction method named scAEQN, which combines autoencoder and quantile normalization. We apply scAEQN to six real scRNA-seq datasets. Experimental results demonstrate that scAEQN can effectively correct the batch effects and retain cell type cluster pattern, and it also has superior performance in downstream analysis using corrected data. In addition, its applicability is also better than other methods.

scAEQN inherits the hypothesis of the QuantNorm method that is the distribution of different cell types among batches is roughly the same. This hypothesis is mentioned in Hicks et al. [16] and is widely accepted so that it helps to avoid deviations caused by confusion with batches. Under this assumption, although batch effects may destroy the overall data distribution, data types in each batch should have similar characteristics, including the similar quantile distribution in different batch blocks in the sample distance matrix.

2.1 Model building

2.1.1 Overview of scAEQN

Figure 1 shows the overall framework of scAEQN. Denotes gene expression count matrix as ${X}_{n\times p}$ containsbatches, whererepresents the number of cells,represents the number of genes. we first use QuantNorm method to calculate corrected sample Pearson correlation coefficient matrix${D}_{n\times n}$, then we construct autoencoder, and add a Pearson correlation coefficient layer to it, which is used to calculate the correlation coefficient matrix ${D}_{n\times n}$ of reconstructed data. The corrected correlation coefficient matrix ${C}_{n\times n}$obtained by the QuantNorm as loss function. Finally, we perform iterative learning on whole model, and take data of low dimensionality for cluster analysis, and identify the differential expression genes using reconstructed data of decoding layer.

2.1.2 Quantile normalization

Quantile normalization (QN) is a regularization technique widely used to remove technical differences because it can produce a good permutation distribution. QN was originally developed for gene expression microarray, but now it has been used in almost kinds of high-dimensional, high-throughput, and omics platforms, including RNA sequencing and proteomics [16]. However, previous studies have found that when using QN, the normalized samples will look very similar, even if the underlying classes are very different. In addition, QN can eliminate true signals and generate false signals during data analysis [17]. QuantNorm [11] is a nonparametric method that can correct the nonlinear batch effects. Because it is feasible that the sample distance matrix is limited. However, in the subsequent downstream analysis, the clustering results of the sample distance matrix are often unsatisfactory and it cannot detect differential expression genes.

It assumes that original data matrix ${X}_{n\times p}$ containsbatches, and $i$-th batch is a matrix of ${n}_{i}\times p$. Firstly, correlation coefficient matrix ${D}_{n\times n}$ is calculated by original gene expression matrix ${X}_{n\times p}$. Because there are negative values in the correlation coefficient matrix. QuantNorm subtracts correlation coefficient matrix by 1 to obtain a dissimilar matrix to ensure that the matrix is non negative. The size of the dissimilarity matrix is $n\times n$, including $m$ intra batch blocks and ${m}^{2}-m$ between batch blocks.

Secondly, the dimensionalities of each batch block are varied, so quantile normalization could be conducted for vectors with different length. QuantNorm adopt an interpolating quantile normalization to normalize the vector with respect to another vector with different lengths, as described below.

Assuming there are two positive real vectors, ${v_1}$ with length ${l_1}$ and ${v_2}$ with length ${l_2}$, the interpolating quantile normalization for ${v_2}$ with respect to ${v_1}$ by steps as follows:

Step 1. Calculate vector ${\tilde{v}}_{i}$,$i=\text{1,2}$ with length ${\tilde{l}}_{i}$,$i=1,2$so that ${\tilde{v}}_{i}$ contains all of non-zero entries in ${v}_{i}$ in ascending order.

Step 2. Calculate scaling vectors ${\tilde{z}}_{i}$,$i=1,2$with length ${\tilde{l}}_{i}$,$i=1,2$that the-th entry in ${\tilde{z}}_{i}$ equals $\frac{k}{{\tilde{l}}_{i}=1}$, where $k=\text{0,1},2,\cdots ,{\tilde{l}}_{i}$.

Step 3. Conduct linear interpolation for ${\tilde{v}}_{1}$ with respect to the scaling vector ${z}_{1}$, then use the predicted values ${\stackrel{-}{v}}_{2}$ at ${z}_{2}$ to replace ${\tilde{v}}_{2}$. Let ${\stackrel{-}{v}}_{i}\left(j\right),j=\text{1,2},\cdots ,{\tilde{l}}_{i}$ represents the $j$-th entry in vector ${\stackrel{-}{v}}_{i}$, then

${\stackrel{-}{v}}_{2}\left(j\right)={\tilde{v}}_{1}\left(j\right)+\frac{{\tilde{v}}_{1}\left(j+1\right)-{\tilde{v}}_{1}\left(j\right)}{{z}_{1}\left(j+1\right)-{z}_{1}\left(j\right)}\left({z}_{2}\left(j\right)-{z}_{1}\left(j\right)\right)$ ,$j=\text{1,2},\cdots ,{\tilde{l}}_{2}$. (1)

Step 4. The ${v_2}$ with all non-zero entries replaced by ${\tilde{v}}_{2}$ in original order of ${v_2}$ is the interpolation quantile normalized vector with respect to${v_1}$.

Thirdly, in order to facilitate optimization of information, QuantNorm selects data with the largest batch block as the reference batch and normalizes interpolation quantiles of other batch blocks. On the basis of interpolation quantile normalization, QuantNorm provides two calculation methods of normalizing other batch blocks with respect to reference batch, vectorize method and row/column iterative method (Table 1). Both they have their own advantages, and the utilization details will be described in Section 2.3.

Table 1

two calculation approaches of dissimilarity matrix in QuantNorm
Algorithm 1. Vectorization approach	Algorithm 2. Iterative approach by rows and columns
Input: original data matrix${X_{n \times p}}$. 1. Compute dissimilarity matrix; 2. Find reference block; 3. Based on the vectorized reference block, conducting the interpolating quantile normalization on the vectorized other blocks; 4. Restoring the vectorized blocks into their original dimensions. Output: dissimilar matrix ${D_{n \times n}}$.	Input: original data matrix${X_{n \times p}}$,tolerance$\varepsilon$,initial value tol. 1. Compute dissimilarity matrix ; 2. Find reference block; 3. While $tol>\varepsilon$do ${D_0}=D$; Compute the interpolating quantile normalization on each row/column of the other blocks ${D_{row}}$and ${D_{column}}$; $D=({D_{row}}+{D_{column}})/2$; $tol={\text{Eulidean}}_{{}}^{{}}{\text{distance}}({D_0},D)$; end Output: dissimilar matrix ${D_{n \times n}}$.

2.1.3 Autoencoder

Autoencoder is an unsupervised deep learning model and it has been widely used in the research of clustering methods, but it cannot be used alone to remove batch effect. Because of its special network structure, it can get the reconstructed data and low-dimensional embedding. Thus, we hope to utilize autoencoder to learn reconstructed data for differential expression genes analysis and low-dimensional embedding for clustering in scRNA-seq data.

We assume that original data is${X_{n \times p}}$as input of the autoencoder. The autoencoder includes three parts, encoder, decoder and correlation coefficient layer, where middle layer is latent low-dimension representation, and last layer is Pearson correlation coefficient of reconstruction layer in decoder. Correlation coefficient layer is defined by Pearson correlation coefficient as follows:

$$\rho \left({X}_{i},{X}_{j}\right)=\frac{cov({X}_{i},{X}_{j})}{\sigma {X}_{i}\sigma {X}_{j}}i,j=\text{1,2},\cdots ,n$$

Where${X_i}$and${X_j}$represent the-th and-th raw vectors of the gene expression matrix respectively. Thus, the Pearson correlation matrix ${C}_{n\times n}$ of reconstructed data is used as the output of total model, which is constructed with the dissimilar matrix ${D}_{n\times n}$ to optimize the whole model through back propagation. The loss function is MSE as follows:

$$MSE(C,D)=\frac{{\sum\limits_{{j=1}}^{n} {\sum\limits_{{i=1}}^{n} {{{({c_{ij}} - {d_{ij}})}^2}} } }}{{{n^2}}}$$

Whererepresents correlation coefficient matrix calculated from the reconstructed data, andrepresents correlation coefficient matrix obtained by the QuantNorm method.represents the number of cells,${c_{ij}}$and${d_{ij}}$represent $i$-th row and $j$-th column elements in matrixand matrixrespectively. In order to minimize the loss function, Adam [34] algorithm is conducted to optimize the total model.

2.2 Experimental design and parameter setting

In scAEQN, there are some hyperparameters, like activation function, threshold value, max iteration, learning rate, need to be set by ourselves, which affect the result of batch correct. For different datasets, the setting of hyperparameters is also different. Thus, we conducted experiments to adjust hyperparameters, utilized six scRNA-seq datasets with different number of samples. Through experiments, we find that the performance is better when threshold of loss function is set around 1e-3 to 1e-2 and step size is set within the range of 1e-4 ~ 1e-5. See the supplementary Table S1-S6 for specific details. We suggest that the loss function threshold is set 5e-3 and the learning rate is set 5e-4. The maximum number of iterations is set to 5000, and the activation function is set to Relu which is the most widely used and effective activation function in deep learning.

When utilizing the QuantNorm method to calculate interpolating quantile normalized similarity matrix, we find that when the number of cells exceeds 1000, the iterative method using “row/column” takes a long time but almost has not to difference from the “vectorize” method. Therefore, when the number of cells is less than 1000, the parameter of method is set to “row/column”, and when the number of cells is more than 1000, the parameter of method is set to “vectorize”.

2.3 Evaluation method

In order to evaluate the performance of scAEQN, we select six public scRNA-seq datasets that the number of cells varied from tens to thousands (Table 2), and different batch correction methods. Specific description and preprocessing work of datasets in the supplementary note S1-S6. Firstly, we evaluate the results of batch correction by low dimensional visualization method t-distributed stochastic neighbor embedding (t-SNE) [27]. We visualize the distribution of batches and cell types of the raw data, and corrected data obtained by MNN [8], Limma [3], Combat [2], QuantNorm [12], scBatch [15], Harmony [11], scVI [14], Scanorama [10] and our method scAEQN. To help to evaluate the effectiveness of scAEQN, we utilize some metrics that measure the level of batch integration quantitatively, such as adjusted rand index (ARI) [28], average silhouette width (ASW) [29], local inverse Simpson’s index (LISI) [30]. After obtaining the batch-corrected outputs, we compute the PCA vectors and used the top 10 PC as inputs to calculate the LISI, ASW, and ARI scores respectively. We first calculate the metrics for assessing cell type purity and batch mixing separately, to obtain stable results, we repeat 20 times to obtain 20 metric scores, normalize their median score to range between 0 to 1 and then combine the assessments into an F1 score, as described in the supplementary note S7-S9.

Cluster analysis is performed on the batch corrected data. We repeat k-means clustering on the corrected count matrix or similarity matrix 10 times. Hyperparameter k of k-means selects the real number of cell types in each dataset. Then, we use ARI [28] to evaluate the consistency between real cell types and k-means clustering labels. The closer the ARI value is to 1, the better the clustering result is. In this way, we can assess the performance of the different batch correction methods through cluster analysis. DEGs analysis plays an important role in downstream of scRNA-seq data analysis. To perform DEG analysis, we use the FindAllMarkers function from the Seurat package [31], with a log₂FC value greater than 0.05 as the threshold. We performed DEG analysis on the gene expression matrix of raw data, which contains batch effects, and on the batch-corrected expression matrix by scAEQN and then visualize the top 5 differential expression genes by heatmap function in ComplexHeatmap [32]. In addition, we also evaluate the feasibility of each batch correction method by comparing the running time on different numbers of cells.

Table 2

Summary of 6 single cell RNA-seq datasets
Datasets	cells	genes	batch	cell type
Mouse and human brain	62	15041	2	6
Mouse embryonic stem cells	469	38616	2	3
Mouse neuron	610	25334	8	4
Human pancreas (different donors)	651	6797	10	4
Mouse pancreas	1886	14878	2	13
Human pancreas (different protocols)	3289	12088	2	13

3.1 Batch correction and integration

In order to show the effectiveness of scAEQN, we have compared the results obtained by different batch correction methods using different datasets. Low dimensional data of latent layer in scAEQN was extracted for visualization and cluster.

We generate t-SNE plots to show correction and integration results in different methods from MESC dataset (Fig. 2B). As shown in Fig. 2B, raw data has obvious batches that it is separated two blocks with different batches in every cluster. As shown on cell type distribution plots, scAEQN and scVI has obvious three cluster blocks contrast with raw data, Combat, MNN, Harmony, Limma, QuantNorm, and scBatch method. They are more obvious for the separation of cell types. However, Scanorama overcorrected on the MESC dataset that cells of different cell types are mixed in every cluster block. In batch distribution plots, scAEQN, scVI, Combat, Harmony, Limma methods have similar levels of batch integration, MNN, QuantNorm, and scBatch methods still retain obvious batch blocks after batch corrected, Scanorama has superior effective batch integration. Likewise, we plot t-SNE plots in different methods on mouse pancreas and human pancreas2 datasets (Supplementary Figure S1 and Figure S2). For the mouse pancreas dataset, it does not have a more obvious batch effect from raw data. For such data, scAEQN still retains a superior clustering pattern of cell types after batch correction, compared with Harmony, scBatch, scVI, Scanorama methods, Combat, Limma, MNN, and QuantNorm methods show obvious batches after batch correction. For the human pancreas2 dataset, scAEQN largely preserves the clustering pattern of cell types although there are subtle batches within each cluster block, compared to Combat, Harmony, Limma, and MNN methods that cell type cluster blocks are mixed. Scanorama method has overcorrected likewise, and scVI method has superior performance on this dataset. In conclusion, by visualizing the distribution of batch and cell type on different datasets, it illustrates that scAEQN can effectively correct the batch effect and retain cell type cluster pattern in scRNA-seq data.

While a clear batch integration of small datasets with only a few cell types or batches is easy to evaluate by visualization method, the assessment can become subjective when the comparison is made between different competitive methods on large complicated datasets [33]. The difficulty increases when the batches are not clearly segregated, like brain, mouse neuron and human pancreas1 datasets. Therefore, we also use three metrics to help us assess the batch correction methods in terms of batch mixing and cell type purity (Fig. 2A). From the ARI F1 score (Fig. 2Aa), scAEQN get the highest ARI F1 score in brain, mouse neuron, and mouse pancreas datasets compared to other methods, the scores are 0.7742, 0.4997, and 0.8441 respectively. scAEQN is second only to scVI on MESC and human pancreas2 datasets, with scores of 0.9479 and 0.9886 on the MESC dataset and 0.7325 and 0.8664 on the human pancreas2 dataset. scAEQN does not perform well on the human pancreas1 dataset, only higher than raw data. Whereas, in the ASW F1 score (Fig. 2Ab), scAEQN is ahead of other batch correction methods with a score 0.6030 on the human pancreas1 dataset. Besides, scAEQN gets the highest ASW F1 scores in mouse pancreas and human pancreas2 datasets. scAEQN is comparable to other methods on brain, MESC, and mouse neuron datasets. For LISI F1score (Fig. 2Ac), scAEQN has superior scores on mouse pancreas and human pancreas2 datasets. scAEQN score of 0.2924 is lower than MNN (0.3790) and Scanorama (0.4610) methods on the brain dataset and its score of 0.2926 is lower than Combat (0.3250), Limma (0.4169), and Harmony (0.4601) methods on MESC dataset. On the mouse neuron dataset, the scAEQN score of 0.3558 is lower than Limma (0.3708), Scanorama (0.4399), and Harmony (0.3714) methods and its score of 0.4415 is lower than Limma (0.4935) and Harmony (0.4724) methods on human pancreas1 datasets.

In general, the metrics of scAEQN are improved compared with raw data, indicating that scAEQN has the ability of batch correction. From the quantitative comparison of three metrics, scAEQN has a better ability of batch correction compared with other methods.

3.2 Cluster analysis

The final purpose of batch correction is to conduct more accurate downstream analysis of data, so it is vital to conduct downstream analysis of batch corrected data. We select six common single-cell RNA-seq datasets (Table 2) and perform batch correction with different batch correction methods for cluster analysis. These methods are BBKNN, Combat, Harmony, Limma, MNN, QuantNorm, scBatch, scVI, and Scanorama. On six datasets, we utilized k-means to cluster the corrected expression data or similarity matrix, and calculate its ARI score. In addition, we also cluster the uncorrected raw matrix with the same clustering method as a control group to verify whether the batch correction is effective before clustering. The clustering results generated by k-means clustering are unstable, we adopt 10 times clustering to get the average value of ARI to get the final evaluation index (Table 3).

Table 3

ARI for different batch correction methods in six datasets
	Brain	MESC	Mouse neuron	Human pancreas1	Mouse pancreas	Human pancreas2
Uncorrected	0.2044	0.4705	0.3202	0.4320	0.1704	0.3255
BBKNN	0.1769	0.5037	0.3621	0.3814	0.1170	0.3401
Combat	0.2179	0.5477	0.3528	0.4308	0.1594	0.3739
Harmony	0.1769	0.5130	0.5127	0.4759	0.1227	0.3409
Limma	0.1635	0.5043	0.3730	0.4114	0.1518	0.3900
MNN	0.1291	0.5027	0.3434	0.2628	0.1129	0.2912
Quantnorm	0.3299	0.2488	0.1178	0.4621	0.5272	0.5207
scbatch	0.2022	0.5344	0.5451	0.5926	0.1434	/
scVI	0.2325	0.9842	0.4221	0.5002	0.3530	0.5867
Scanorama	0.0871	0.0769	0.1844	0.0187	0.1259	0.3418
scAEQN	0.4992	0.8469	0.6733	0.6612	0.7832	0.6077

From the comparison of ARI metric results, the performance of scAEQN is effective than the rest of nine batch corrected methods on brain, mouse neuron, human pancreas1, mouse pancreas and human pancreas2 datasets, which improves at least 2% higher than the best of the nine methods and improves 26% than the second on mouse pancreas dataset. On mouse embryonic stem cells (MESC) dataset, scAEQN is second to scVI but it still outperformed the other methods and improved by 37% over raw data. We did not get the clustering ARI metric of the human pancreas2 dataset after correction by the scBatch method, because there are too many cells in the human pancreas2 dataset, and scBatch cannot process datasets with too many cells, which is also mentioned by the author. From the bar plot of ARI (Fig. 3A), scAEQN has obvious advantages in the cluster. Based on this, we demonstrate that scAEQN improves clustering accuracy while removing batch effects, and is superior to most batch correction methods.

3.3 Differential expression genes analysis

To further assess scAEQN, DEGs analysis is performed from the corrected data by scAEQN and compared with the uncorrected original data. We utilize FindAllMarkers function in Seurat3 package to obtain DEGs on six datasets, where log₂FC threshold is 0.05 because of the corrected data are normalized and log transformed in scAEQN. The top 5 DEGs are selected by log2FC rank in every cell type on six datasets, and then plotted heatmap by ComplexHeatmap [28] method (Fig. 3B). We can see that the top 5 DEGs obtained by scAEQN have more distinct expression blocks in each cell type, compared with uncorrected raw data from Fig. 3B (a-d). scAEQN has distinct differential expression gene blocks in large cell types, like alpha, beta, ductal, delta, endothelial, and gamma cell types on the mouse pancreas (Fig. 3Be). And on the human pancreas2 dataset, scAEQN has also distinct differential expression gene blocks in large cell types, like alpha, beta, ductal, delta, acinar, and gamma cell types (Fig. 3Bf).

Based on the human pancreas1 dataset, we analyze the identification results of DEGs obtained by different batch correction methods using Seurat [15]. We compare our method scAEQN with the data of raw count matrix, MNN [10], Limma [3], Combat [2], and scBatch [13], eventually. We set the conditions of adjusted P-value < 1e-2 and log₂FC value > 0.05 to analysis the number of DEGs (Fig. 3D). More DEGs are detected by scAEQN under the same conditions, which avoid missing possible genes, and provide a solid foundation for gene function analysis. At the same time, the number of DEGs, detected jointly by scAEQN with uncorrected, MNN, Limma, Combat, scBatch, is 279, 444, 206, 224, 420. Nearly half of the DEGs detected by other methods are consistent with scAEQN. Finally, we conclude the common DEGs obtained by these six methods, and obtain 165 co-expressed DEGs. The specific detected DEGs and co-expressed genes are shown in the supplementary note S10. In conclusion, scAEQN is feasible in the identification and analysis of differential expression genes.

3.4 Running time

To verify the operability of scAEQN, we record running time in seconds of different batch correction methods under five datasets with different cell numbers varied from tens to thousands, and parameters setting are described in section 2.2.

As can be seen from Fig. 3C, the abscissa represents the number of cells arranged from small to large in each dataset, the ordinate represents the running time in seconds, and the lines of different colors represent the running time of different batch corrections methods. Figure 3C(b) is an enlarged plot of the line marked by the red circle in Fig. 3C(a). There is a dotted line between 1886 and 3289 in Fig. 3C(a), because a total of 12 hours were spent on the dataset of 3289 cells using the scBatch method and did not get results, and spent around 4 hours on the dataset of 1886 cells. The running time of the scBatch method increases exponentially with the increase of the number of cells. MNN method also takes more time with the increase in the number of cells similarly, but the running time is still acceptable. Combat, Limma, and BBKNN are stable in running time. The running time of scAEQN is related to the parameters of the selected QuantNorm. When the parameters of each dataset in scAEQN are set to “method = vectorization”, the running time of our method will increase with the increase of the number of cells, but the longest time is only a few minutes. For instance, on the dataset with 651 cells of the human pancreas from different donors, when the parameter is set to “method = row/column”, the running time for scAEQN is only six minutes. Although scAEQN is not the best in running time performance, it is still within our acceptable range.

Batch effect is a common problem in omics data analysis, but there is also batch effect in single cell data analysis, which affects the results of downstream analysis. Therefore, batch correction is very vital before downstream analysis. In this paper, we propose a batch correction joint dimension reduction method named scAEQN, which can obtain dimension reduction data and batch corrected data simultaneously, and provides great convenience for cluster analysis and differential expression genes analysis. We applied scAEQN to six datasets and compared it with the other nine batch correction methods, and found that scAEQN has been greatly improved in cluster analysis, which improves up to 26% than the best of the nine compared methods. At the same time, we analyzed differential expression genes on the corrected data with uncorrected raw data, and found that scAEQN can improve ability of differential expression analysis and detect more differential expression genes, and the rate of gene coincidence with other methods is at least 42%.

scAEQN is based on the assumption that the distribution between batches and cell types are roughly balanced. This assumption is reasonable, and many algorithms based on this assumption. In human pancreas dataset from different donors, the distribution of batch and sample type is slightly unbalanced, but the method under this assumption is still effective. However, it is worth mentioning that datasets often have complex batch effects, and real datasets with a balanced distribution of batch and cell types are often difficult to find. Moreover, due to the limitations of quantile normalization, it is easy to fall into extremes, so it may also produce non-existent cluster blocks [12], and cause the lack of some biological information. Besides, there is still room for improvement in scAEQN. In the design of the autoencoder, the use of the random gradient descent method (Adam) makes our results less robust [34], and we need to run many times to obtain better results. In addition, scAEQN only uses the most common Pearson correlation coefficient matrix to measure the correlation of count data. There are also many correlation coefficients, such as Spearman correlation coefficient, Kendall correlation coefficient, and Cosine correlation coefficient that are often utilized to measure the similarity between cells in single-cell RNA-seq data. We believe that there is still room for improvement in the measurement of similarity data. Despite such shortcomings, scAEQN still has very important application significance. It provides convenient and beneficial help for downstream analysis in single-cell RNA-seq data analysis.

Supplementary information

Supplementary material is available online.

Acknowledgements

This work was supported by the open fund of Information Materials and Intelligent Sensing Laboratory of Anhui Province (Grant No. IMIS202105); the Xinjiang Autonomous Region University Research Program (No. XJEDU2019Y002); and the National Natural Science Foundation of China (No. U19A2064, 61873001).

Author contributions

SSN contributions are processed a batch correction joint dimensionality reduction method scAEQN and the numerical simulation results of scAEQN model. ZJP, SYS and LDM contribution lies in the embellishment of the article. WHY contribution lies in the discussion and revision of the article. ZCH contribution lies in the thought guidance of the method. All authors read and approved the final manuscript.

Funding

Not applicable

Availability of data and materials

The datasets supporting this study are available from public databases: human and mouse brain scRNA-seq dataset (brain): GSE73721 in the NCBI Gene Expression Omnibus (GEO); mouse embryonic stem cell dataset (MESC): E-MTAB-2600 in the Array Express database; mouse neurons dataset: GSE59739 in the NCBI GEO; human pancreas datasets from different donors (human pancreas1): GSE81608 in the NCBI GEO; human pancreas dataset from different protocols (human pancreas2): GSE81076 and GSE85241 in the NCBI GEO; mouse pancreas dataset: GSE84133 in the NCBI GEO. MESC, mouse neurons, and human pancreas1 datasets that have been preprocessed available in the public database of the Hemberg group in https://hemberglab.github.io/scRNA.seq.datasets. All preprocessed datasets can be downloaded on GitHub: https://github.com/SiningSong/scAEQN/datasets. The source code of scAEQN is available at https://github.com/SiningSong/scAEQN.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interest.

Author details

¹College of Mathematics and System Sciences, Xinjiang University, Urumqi, China.²Institute of Mathematics and Physics, Xinjiang University, Urumqi, China. ³School of Artificial Intelligence, Anhui University, Hefei, China

Hicks SC, Townes FW, Teng M, et al. Missing data and technical variability in single-cell RNA-sequencing experiments. Biostatistics, 2018, 19: 562–578.
Johnson W E, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics, 2007,8: 118–127.
Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res.,2015, 43: e47.
Zhang Y, Parmigiani G, Johnson W E. Combat-seq: batch effect adjustment for RNA-seq count data. NAR Genom Bioinform. 2020, 2(3): lqaa078.
Yang I V. Use of external controls in microarray experiments. Methods Enzymol.,2006, 411: 50–63.
Jiang L, Schlesinger F, Davis C A, et al. Synthetic spike-in standards for RNA-seq experiments. Genome Res.,2011, 21: 1543–1551.
Luo X and Wei Y. Batch effects correction with unknown subtypes. J. Am. Stat. Assoc.,2019, 114: 581–594.
Haghverdi L, Lun ATL, Morgan MD, et al. Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. Nat Biotechnol.,2018, 36(5): 421–427.
Polański K, Young MD, Miao Z, et al. BBKNN: fast batch alignment of single cell transcriptomes. Bioinformatics, 2020, 36(3): 964–965.
Hie B, Bryson B, Berger B. Efficient integration of heterogeneous single-cell transcriptomes using Scanorama. Nat Biotechnol. 2019;37: 685–91
Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, Baglaenko Y, Brenner M, Loh P-r, Raychaudhuri S. Fast, sensitive and accurate integration of single-cell data with Harmony. Nature Methods; 2019.
Fei T, Zhang T, Shi W, et al. Mitigating the adverse impact of batch effects in sample pattern detection. Bioinformatics, 2018, 34: 2634–2641.
Shaham U, Stanton KP, Zhao J, et al. Removal of batch effects using distribution-matching residual networks. Bioinformatics, 2017, 33: 2539–2546.
Lopez R, Regier J, Cole M B, et al. Deep generative modeling for single-cell transcriptomics. Nature methods, 2018, 15(12): 1053–1058.
Fei T, Yu T. scBatch: batch-effect correction of RNA-seq data through sample distance matrix adjustment. Bioinformatics, 2020, 36(10): 3115–3123.
Zhao Y, Wong L, Goh WWB. How to do quantile normalization correctly for gene expression data analyses. Sci. Rep., 2020, 10(1): 15534.
Wang D, Cheng L, Wang M, et al. Extensive increase of microarray signals in cancers calls for novel normalization assumptions. Comput. Biol. Chem., 2011, 35:126–130.
Hicks SC, Townes FW, Teng M, et al. Missing data and technical variability in single-cell RNA-sequencing experiments. Biostatistics, 2018,19: 562–578.
Zhang Y, Sloan SA, Clarke LE, et al. Purification and characterization of progenitor and mature human astrocytes reveals transcriptional and functional differences with mouse. Neuron, 2016, 89: 37–53.
Kim JK, Kolodziejczyk AA, Ilicic T, et al. (2015) Characterizing noise structure in single-cell RNA-seq distinguishes genuine from technical stochastic allelic expression. Nat.Commun., 2015, 6: 8687.
Usoskin D, Furlan A, Islam S, et al. Unbiased classification of sensory neuron types by large-scale single-cell RNA sequencing. Nat. Neurosci., 2015, 18: 145–153.
Xin Y, Kim J, Okamoto H, et al. RNA sequencing of single human islet cells reveals type 2 diabetes genes. Cell Metab., 2016,24: 608–615.
Grün D, Muraro MJ, Boisset JC, et al. De Novo Prediction of Stem Cell Identity using Single-Cell Transcriptome Data. Cell Stem Cell, 2016, 19(2):266–277.
Muraro MJ, Dharmadhikari G, Grün D, et al. A Single-Cell Transcriptome Atlas of the Human Pancreas. Cell Syst, 2016, 3(4):385–394.
Wolf F A, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol.2018, 19(1):15.
Baron M, Veres A, Wolock S L, et al. A single-cell transcriptomic map of the human and mouse pancreas reveals inter-and intra-cell population structure. Cellsystems,2016,3(4): 346–360.
Maaten LVD, Hinton G. Visualizing data using t-SNE. Mach. Learn. Res.,2008,9: 2579–2605.
Hubert L and Arabie P. Comparing partitions. Classification,1985,2:193– 218.
Rousseeuw PJ. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. J Comput Appl Math. 1987; 20:53–65
Buttner M, Miao Z, Wolf FA, Teichmann SA, Theis FJ. A test metric for assessing single-cell RNA-seq batch correction. Nat Methods. 2019; 16:43–9.
Satija R, Farrell JA, D Gennert, et al. Spatial reconstruction of single-cell gene expression data. Nat. Biotechnol.,2015, 33: 495–502.
Gu Z, Eils R, and Schlesner M, et al. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics, 2016, 32(18): 2847–9.
Tran, H.T.N., Ang, K.S., Chevrier, M. et al. A benchmark of batch-effect correction methods for single-cell RNA sequencing data. Genome Biol., 2020, 21(1).
Kingma DP, Ba J. Adam: A method for stochastic optimization. arXiv preprint arXiv,2014,1412:6980.
Chen G, Ning B, Shi T. Single-Cell RNA-Seq Technologies and Related Computational Data Analysis. Front. Genet.,2019,10:317.

No competing interests reported.

Supplementarydocument.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

scAEQN: A Batch Correction Joint Dimension Reduction Method on scRNA- seq Data

Status:

Version 1

Abstract

Background

Results

Conclusion

Figures

1. Introduction

2. materials And Methods

2.1 Model building

2.1.1 Overview of scAEQN

2.1.2 Quantile normalization

2.1.3 Autoencoder

2.2 Experimental design and parameter setting

2.3 Evaluation method

3. Results

3.1 Batch correction and integration

3.2 Cluster analysis

3.3 Differential expression genes analysis

3.4 Running time

4. Conclusion And Discussion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1