A Comprehensive Multi-Center Cross-platform Benchmarking Study of Single-cell RNA Sequencing Using Reference Samples

Microarray batch effect (BE) has been the primary bottleneck for large-scale integration of data from multiple experiments. Current BE correction methods either need known batch identities (ComBat) or have the potential to overcorrect, by removing true but unknown biological differences (Surrogate Variable Analysis SVA). It is well known that experimental conditions such as array or reagent batches, PCR amplification or ozone levels can affect the measured expression levels; often the direction of perturbation of the measured expression is the same in different datasets. However, there are no BE correction algorithms that attempt to estimate the individual effects of technical differences and use them to correct expression data. In this manuscript, we show that a set of signatures, each of which is a vector the length of the number of probes, calculated on a reference set of microarray samples can predict much of the batch effect in other validation sets. We present a rationale of selecting a reference set of samples designed to estimate technical differences without removing biological differences. Putting both together, we introduce the Batch Effect Signature Correction (BESC) algorithm that uses the BES calculated on the reference set to efficiently predict and remove BE. Using two independent validation sets, we show that BESC is capable of removing batch effect without removing unknown but true biological differences. Much of the variations due to batch effect is shared between different microarray datasets. That shared information can be used to predict signatures (i.e. directions of perturbation) due to batch effect in new datasets. The correction can be precomputed without using the samples to be corrected (blind), done on each sample individually (single sample) and corrects only known technical effects without removing known or unknown biological differences (conservative). Those three characteristics make it ideal for high-throughput correction of samples for a microarray data repository. We also compare the performance of BESC to three other batch correction methods: SVA, Removing Unwanted Variation (RUV) and Hidden Covariates with Prior (HCP). An R Package besc implementing the algorithm is available from http://explainbio.com.

Based on benchmarking various methods and analyzing multiple real scRNA-seq datasets, a computational platform/workflow and a set of tips for best practices are developed for analyzing causal interactions or relationships in single cells.

Mass spectrometry is a powerful tool for identifying and analyzing biomolecules such as metabolites and lipids in complex biological samples. Liquid chromatography and gas chromatography mass spectrometry studies quite commonly involve large numbers of samples, which can require significant time for sample preparation and analyses. To accommodate such studies, the samples are commonly split into batches. Inevitably, variations in sample handling, temperature fluctuation, imprecise timing, column degradation, and other factors result in systematic errors or biases of the measured abundances between the batches. Numerous methods are available via R packages to assist with batch correction for omics data; however, since these methods were developed by different research teams, the algorithms are available in separate R packages, each with different data input and output formats. We introduce the malbacR package, which consolidates 11 common batch effect correction methods for omics data into one place so users can easily implement and compare the following: pareto scaling, power scaling, range scaling, ComBat, EigenMS, NOMIS, RUV-random, QC-RLSC, WaveICA2.0, TIGER, and SERRF. The malbacR package standardizes data input and output formats across these batch correction methods. The package works in conjunction with the pmartR package, allowing users to seamlessly include the batch effect correction in a pmartR workflow without needing any additional data manipulation.

In a single-cell RNA-seq (scRNA-seq) data set, a high proportion of missing values (or an excessive number of zeroes) are frequently observed. For the related follow-up tasks, such as clustering analysis and differential expression analysis, a data set without missing values is generally required. Many imputation approaches have been proposed for this purpose. Multiple imputation (MI) is a well-established approach to address possible biases in a follow-up analysis result based on one-time imputed data. There is a lack of investigation on this in the analysis of scRNA-seq data. In this study, we have investigated how to efficiently apply the MI approach to the clustering analysis and the differential expression analysis of scRNA-seq data. We proposed an MI procedure for clustering analysis and an MI procedure for differential expression analysis. To demonstrate the improvements achieved by MI in clustering analysis and differential expression analysis of scRNA-seq data, we analyzed three well-known scRNA-seq data sets. scIGANs, an scRNA-seq imputation method based on the generative adversarial networks (GANs), has been recently proposed for scRNA-seq data imputation. Multiple randomly imputed data sets can be conveniently generated by this method. We implemented our MI procedures based on scIGANs. We demonstrated that MI yielded improved performances on the clustering analysis and differential expression analysis results. Our applications to experimental scRNA-seq data illustrated the advantages of MI over one-time imputation of missing values in scRNA-seq data.

Introduction: Tissue microarrays (TMAs) have been used in thousands of cancer biomarker studies. It is unknown to what extent TMAs are affected by batch effects, i.e., measurement error in biomarker levels between batches (slides from TMAs), what impact batch effects have on scientific inference from TMAs, and how to correct for batch effects in the analysis phase. Methods: We leveraged 14 prostate cancer TMAs from the tumors of 1,464 men, prospectively diagnosed 1982-2009 and followed for survival through 2017. Twenty protein biomarkers were measured, using eye-scoring or computational scoring of a median of 3 cores per tumor. We quantified the maximal extent of potential batch effects (between-TMA variance) using intraclass correlation coefficients. We developed a free R package (batchtma) and implemented several methods for batch effect correction, focusing on batch means (simple means, marginal standardization, inverse probability weighting), mean and variance (ComBat), background and dynamic range (quantile regression), ranks (quantile normalization), and stratification-based approaches. Some methods attempted to retain biomarker differences from differing clinical and pathological characteristics of TMAs. We evaluated the performance of correction methods using plasmode simulation and quantified the impact of batch effect correction on common longitudinal biomarker analyses. Results: In 10 of the 20 biomarkers, more than 10% of biomarker variance was attributable to between-TMA differences (range, 1 to 47%). We observed differences in clinical and pathological characteristics between TMAs with secular changes; for example, between 11% and 32% of tumors per TMA were high-grade (Gleason 8-10). In parallel, we confirmed batch effects between multiple cores from the same tumors included on >1 TMA and between tumors with similar characteristics across TMAs. Batch effect-corrected biomarker levels were more similar between different correction methods compared to uncorrected values. For biomarkers with little between-TMA differences, associations of biomarker levels and survival outcomes were similar with and without batch effect correction. For biomarkers with higher between-TMA differences, Cox regression results after batch effect correction were substantially different from results based on uncorrected values, with up to twofold differences in hazard ratios. In plasmode simulation, all correction methods were better than no batch effect correction. Conclusions: Some extent of batch effects should be expected in TMA-based cancer biomarker studies. Batch effects should be quantified and, if necessary, corrected in any biomarker study using more than one TMA. Anticipating batch effects during study design and, for biomarkers found to have batch effects, utilizing straightforward analytical correction approaches increase the validity of TMA-based studies. Citation Format: Konrad H. Stopsack, Molin Wang, Svitlana Tyekucheva, Travis A. Gerke, J. Bailey Vaselkiv, Kathryn L. Penney, Philip W. Kantoff, Stephen P. Finn, Michelangelo Fiorentino, Massimo Loda, Giovanni Parmigiani, Lorelei A. Mucci. Batch effects in tumor biomarker studies using tissue microarrays: Extent, impact, and remediation [abstract]. In: Proceedings of the American Association for Cancer Research Annual Meeting 2021; 2021 Apr 10-15 and May 17-21. Philadelphia (PA): AACR; Cancer Res 2021;81(13_Suppl):Abstract nr 893.

Single-cell RNA sequencing (scRNA-seq) data allows to create cell type specific transcriptome profiles. Such profiles can be aligned with genome-wide association studies (GWASs) to implicate cell type specificity of the traits. Current methods typically rely only on a small subset of available scRNA-seq datasets, and integrating multiple datasets is hampered by complex batch effects. Here we collated 43 publicly available scRNA-seq datasets. We propose a 3-step workflow with conditional analyses within and between datasets, circumventing batch effects, to uncover associations of traits with cell types. Applying this method to 26 traits, we identify independent associations of multiple cell types. These results lead to starting points for follow-up functional studies aimed at gaining a mechanistic understanding of these traits. The proposed framework as well as the curated scRNA-seq datasets are made available via an online platform, FUMA, to facilitate rapid evaluation of cell type specificity by other researchers.

The incidence and mortality rates of lung cancer are exceptionally high. Many patients are diagnosed with early stage lung cancer but experience rapid recurrence post-surgery. Many research studies have shown that the unfavorable prognosis of patients may be associated with micro-metastasis in the lymph nodes. Our research aimed to develop a nomogram to predict the prognosis of lung adenocarcinoma (LUAD). Single-cell RNA sequencing (scRNA-seq) data were analyzed to identify 11 cell clusters. Patterns of incoming and outgoing signals were identified across the entire cell population. A weighted gene co-expression network analysis (WGCNA) was conducted to uncover critical genes in LUAD. The intersecting marker genes were used to construct the prognostic model. scRNA-seq data were analyzed to identify 19 cell clusters. We identified 3,464 marker genes from the scRNA-seq dataset, 1,994 differentially expressed genes from the bulk RNA sequencing (RNA-seq) dataset, and 1,863 genes associated with a key module identified by the WGCNA. After performing the intersection, univariate Cox, and least absolute shrinkage and selection operator analyses, a prognostic model was established based on the expression levels of 13 signature genes. Subsequent functional experiments confirmed the role of selected regulated genes. Through the integration of scRNA-seq data and bulk RNA-seq data, we developed an innovative model to predict the prognosis of patients. The risk score was found to be a significant independent predictor and clinical-pathological features of LUAD.

Mapbatch: Conservative Batch Normalization for Single Cell RNA-Sequencing Data Enables Discovery of Rare Cell Populations in a Multiple Myeloma Cohort

In this paper, we introduce Gene Knockout Inference (GenKI), a virtual knockout (KO) tool for gene function prediction using single-cell RNA sequencing (scRNA-seq) data in the absence of KO samples when only wild-type (WT) samples are available. Without using any information from real KO samples, GenKI is designed to capture shifting patterns in gene regulation caused by the KO perturbation in an unsupervised manner and provide a robust and scalable framework for gene function studies. To achieve this goal, GenKI adapts a variational graph autoencoder (VGAE) model to learn latent representations of genes and interactions between genes from the input WT scRNA-seq data and a derived single-cell gene regulatory network (scGRN). The virtual KO data is then generated by computationally removing all edges of the KO gene—the gene to be knocked out for functional study—from the scGRN. The differences between WT and virtual KO data are discerned by using their corresponding latent parameters derived from the trained VGAE model. Our simulations show that GenKI accurately approximates the perturbation profiles upon gene KO and outperforms the state-of-the-art under a series of evaluation conditions. Using publicly available scRNA-seq data sets, we demonstrate that GenKI recapitulates discoveries of real-animal KO experiments and accurately predicts cell type-specific functions of KO genes. Thus, GenKI provides an in-silico alternative to KO experiments that may partially replace the need for genetically modified animals or other genetically perturbed systems.

BackgroundWe developed a novel approach to minimize batch effects when assigning samples to batches. Our algorithm selects a batch allocation, among all possible ways of assigning samples to batches, that minimizes differences in average propensity score between batches. This strategy was compared to randomization and stratified randomization in a case–control study (30 per group) with a covariate (case vs control, represented as β1, set to be null) and two biologically relevant confounding variables (age, represented as β2, and hemoglobin A1c (HbA1c), represented as β3). Gene expression values were obtained from a publicly available dataset of expression data obtained from pancreas islet cells. Batch effects were simulated as twice the median biological variation across the gene expression dataset and were added to the publicly available dataset to simulate a batch effect condition. Bias was calculated as the absolute difference between observed betas under the batch allocation strategies and the true beta (no batch effects). Bias was also evaluated after adjustment for batch effects using ComBat as well as a linear regression model. In order to understand performance of our optimal allocation strategy under the alternative hypothesis, we also evaluated bias at a single gene associated with both age and HbA1c levels in the ‘true’ dataset (CAPN13 gene).ResultsPre-batch correction, under the null hypothesis (β1), maximum absolute bias and root mean square (RMS) of maximum absolute bias, were minimized using the optimal allocation strategy. Under the alternative hypothesis (β2 and β3 for the CAPN13 gene), maximum absolute bias and RMS of maximum absolute bias were also consistently lower using the optimal allocation strategy. ComBat and the regression batch adjustment methods performed well as the bias estimates moved towards the true values in all conditions under both the null and alternative hypotheses. Although the differences between methods were less pronounced following batch correction, estimates of bias (average and RMS) were consistently lower using the optimal allocation strategy under both the null and alternative hypotheses.ConclusionsOur algorithm provides an extremely flexible and effective method for assigning samples to batches by exploiting knowledge of covariates prior to sample allocation.

SummarySingle-cell RNA-sequencing (scRNA-seq) technology enables studying gene expression programs from individual cells. However, these data are subject to diverse sources of variation, including ‘unwanted’ variation that needs to be removed in downstream analyses (e.g. batch effects) and ‘wanted’ or biological sources of variation (e.g. variation associated with a cell type) that needs to be precisely described. Surrogate variable analysis (SVA)-based algorithms, are commonly used for batch correction and more recently for studying ‘wanted’ variation in scRNA-seq data. However, interpreting whether these variables are biologically meaningful or stemming from technical reasons remains a challenge. To facilitate the interpretation of surrogate variables detected by algorithms including IA-SVA, SVA or ZINB-WaVE, we developed an R Shiny application [Visual Surrogate Variable Analysis (V-SVA)] that provides a web-browser interface for the identification and annotation of hidden sources of variation in scRNA-seq data. This interactive framework includes tools for discovery of genes associated with detected sources of variation, gene annotation using publicly available databases and gene sets, and data visualization using dimension reduction methods.Availability and implementationThe V-SVA Shiny application is publicly hosted at https://vsva.jax.org/ and the source code is freely available at https://github.com/nlawlor/V-SVA.Contactleed13@miamioh.edu or duygu.ucar@jax.orgSupplementary informationSupplementary data are available at Bioinformatics online.

Bulk transcriptomes are an essential data resource for understanding basic and disease biology. However, integrating information from different experiments remains challenging because of the batch effect generated by various technological and biological variations in the transcriptome. Numerous batch-correction methods to deal with this batch effect have been developed in the past. However, a user-friendly workflow to select the most appropriate batch-correction method for the given set of experiments is still missing. We present the SelectBCM tool that prioritizes the most appropriate batch-correction method for a given set of bulk transcriptomic experiments, improving biological clustering and gene differential expression analysis. We demonstrate the applicability of the SelectBCM tool on analyses of real data for two common diseases, rheumatoid arthritis and osteoarthritis, and one example to characterize a biological state, where we performed a meta-analysis of the macrophage activation state. The R package is available at https://github.com/ebi-gene-expression-group/selectBCM.

The notable progress in single-cell RNA sequencing (ScRNA-seq) technology is beneficial to accurately discover the heterogeneity and diversity of cells. Clustering is an extremely important step during the ScRNA-seq data analysis. However, it cannot achieve satisfactory performances by directly clustering ScRNA-seq data due to its high dimensionality and noise. To address these issues, we propose a novel ScRNA-seq data representation model, termed Robust Graph regularized Non-Negative Matrix Factorization with Dissimilarity and Similarity constraints (RGNMF-DS), for ScRNA-seq data clustering. To accurately characterize the structure information of the labeled samples and the unlabeled samples, respectively, the proposed RGNMF-DS model adopts a couple of complementary regularizers (i.e., similarity and dissimilar regularizers) to guide matrix decomposition. In addition, we construct a graph regularizer to discover the local geometric structure hidden in ScRNA-seq data. Moreover, we adopt the l2,1-norm to measure the reconstruction error and thereby effectively improve the robustness of the proposed RGNMF-DS model to the noises. Experimental results on several ScRNA-seq datasets have demonstrated that our proposed RGNMF-DS model outperforms other state-of-the-art competitors in clustering.

Joint learning of multiple gene networks from single-cell gene expression data

BackgroundNeuroblastoma (NB) is a heterogeneous pediatric malignancy with highly variable outcomes. Traditional clinical factors, such as stage, MYCN status, and patient age, often fail to fully capture disease complexity. Recent advances in single-cell sequencing and integrative transcriptomic analyses provide an opportunity to identify more precise prognostic biomarkers and guide individualized therapies. This study aimed to develop and validate a robust prognostic model for NB by integrating single-cell and bulk transcriptomic data.MethodsWe integrated seven publicly available single-cell RNA sequencing (scRNA-seq) datasets to form the Neuroblastoma Atlas and stratified tumor cells into high-risk and intermediate/low-risk groups. Differentially expressed genes (DEGs) were identified using defined fold-change and expression thresholds. Candidate genes were further validated using bulk RNA-sequencing data (HRA002064) and intersected with essential genes defined by DepMap Computational estimation of gene dependency using Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR) screening (CERES) scores. From the resulting set, we constructed a multivariate Cox regression-based prognostic model using the GSE49710 dataset. The model’s performance was evaluated via Kaplan-Meier curves, time-dependent receiver operating characteristic analysis, and decision curve analysis. External validation was performed in the E-MTAB-8248 dataset, and comparisons with standard clinical indicators (International Neuroblastoma Staging System, MYCN status, age) were conducted.ResultsIntegrating scRNA-seq and bulk RNA-seq data identified 123 overlapping DEGs, of which seven genes (BIRC5, CDC2, GINS2, MAD2L1, ORC6L, RRM2, TOP2A) were further prioritized based on CERES dependency scores. Multivariate Cox regression and collinearity screening yielded a four-gene prognostic model (RiskScore) that significantly discriminated high- and low-risk patients in the GSE49710 cohort. Compared to traditional indicators, the four-gene model demonstrated superior time-dependent area under the curve (AUC) values and clinical decision-making benefits. External validation in E-MTAB-8248 confirmed the model’s robust predictive performance, with high AUC values at 1-, 3-, and 5-year time points and consistent superiority over clinical parameters.ConclusionsThis study presents a novel four-gene prognostic signature derived from integrative scRNA-seq and bulk RNA-seq analyses. The model outperforms established clinical factors in predicting NB outcomes and provides added clinical decision-making value. Further prospective validation and mechanistic investigations may facilitate the translation of this prognostic signature into routine clinical practice, enabling more refined risk stratification and personalized treatment strategies for children with NB.