A deconvolution framework that uses single-cell sequencing plus a small benchmark data set for accurate analysis of cell type ratios in complex tissue samples

Abstract

Bulk deconvolution with single-cell/nucleus RNA-seq data is critical for understanding heterogeneity in complex biological samples, yet the technological discrepancy across sequencing platforms limits deconvolution accuracy. To address this, we utilize an experimental design to match inter-platform biological signals, hence revealing the technological discrepancy, and then develop a deconvolution framework called DeMixSC using this well-matched, that is, benchmark, data. Built upon a novel weighted nonnegative least-squares framework, DeMixSC identifies and adjusts genes with high technological discrepancy and aligns the benchmark data with large patient cohorts of matched-tissue-type for large-scale deconvolution. Our results using two benchmark data sets of healthy retinas and ovarian cancer tissues suggest much-improved deconvolution accuracy. Leveraging tissue-specific benchmark data sets, we applied DeMixSC to a large cohort of 453 age-related macular degeneration patients and a cohort of 30 ovarian cancer patients with various responses to neoadjuvant chemotherapy. Only DeMixSC successfully unveiled biologically meaningful differences across patient groups, demonstrating its broad applicability in diverse real-world clinical scenarios. Our findings reveal the impact of technological discrepancy on deconvolution performance and underscore the importance of a well-matched data set to resolve this challenge. The developed DeMixSC framework is generally applicable for accurately deconvolving large cohorts of disease tissues, including cancers, when a well-matched benchmark data set is available.

Figure 1.

Assessing technological discrepancy between bulk and single-cell sequencing platforms using matched single-nucleus aliquots. (A) Workflow for generating a benchmark data set. We collect 24 healthy human retinal samples within 6 h postmortem. An illustration shows the layer and cell compositions of the human retina. Seven major cell types include photoreceptors (rod and cone cells), bipolar cells (BCs), retinal ganglion cells (RGCs), horizontal cells (HCs), amacrine cells (ACs), and Müller glia cells (MGs). Three minor cell types are not depicted in the illustration: astrocytes, microglia cells, and retinal pigment epithelial cells (RPEs). Samples are isolated into single-nucleus suspensions. The same aliquot of single nucleus is used for both bulk and snRNA-seq profiling. The matched pseudobulk mixtures are generated as conventionally done by summing UMI counts across cells from all cell types in each sample. This data generation pipeline guarantees the matched bulk and snRNA-seq data share the same cell type proportions, which enables us to evaluate the impact of technological discrepancy (i.e., the shot-gun sequencing procedure) on the bulk and snRNA-seq expression profiles. (B,C) The influence of technological discrepancy at the sample and gene level, respectively. (B) Spearman's correlation coefficient across genes between the matched real-bulk and pseudobulk RNA-seq data for one sample at a time for both batches. The correlations were calculated using quantile-normalized expression data (relative abundances). (C) MA-plots displaying the mean expression levels of all genes between matched real-bulk and pseudobulk data. Differentially expressed (DE) genes are identified using the paired t-test with Benjamini–Hochberg (BH) adjustment. Red represents genes expressed higher in the real bulk, and blue represents genes expressed higher in the pseudobulk. The horizontal dotted lines denote a twofold change between matched real-bulk and pseudobulk data. (adj.p) Adjusted P-values. (D) Venn diagrams showing genes consistently expressed higher in the bulk (top, overlap of red dots in panel C) or the snRNA-seq generated pseudobulk (bottom, overlap of blue dots in panel C) between the two batches, which were generated using different tissue samples and a different time.

Figure 2.

Overview of DeMixSC. The DeMixSC framework for deconvolution analysis of bulk RNA-seq data using sc/sn RNA-seq data as a reference. (A) The framework starts with a benchmark data set of matched bulk and sc/snRNA-seq data with the same cell type proportions. Pseudobulk mixtures are generated from the sc/sn data. DeMixSC identifies genes in G₁ and G₂ with the matched real-bulk and pseudobulk data. The non-DE genes are considered stably captured by both sequencing platforms (blue), whereas the DE genes are more impacted by the technological discrepancy (orange). (B) DeMixSC then employs a normalization procedure to perform the alignment between two bulk RNA-seq data sets (e.g., with ComBat). (C) DeMixSC estimates cell type proportions under a weighted nonnegative least square (wNNLS) framework with two improvements: (1) partitioning and adjusting genes with high technological discrepancy and (2) a new weight function. The final estimates are obtained when the algorithm either converges or reaches the prespecified maximum number of iterations. Here, G₁ is genes with low technological discrepancy, G₂ is genes with high technological discrepancy, a is a user-defined positive constant that serves as an adjustment factor, $\hat{\mathit{r}}$ is the reference matrix derived from the sc/snRNA-seq data, y is the observed expression in bulk RNA-seq data, $\hat{\mathit{p}}$ is the vector of estimated cell type proportions, and $\hat{\mathit{w}}$ is the estimated gene weights.

Figure 3.

Comparing the estimation accuracy of DeMixSC to existing deconvolution methods. (A) Workflow for the deconvolution benchmarking design. We use benchmark data from retinal samples. The cell count proportions for each cell type are used as ground truth for the corresponding tissue samples. We assess the deconvolution performance of DeMixSC and seven existing methods for both bulk and pseudobulk mixtures. In addition to the raw counts, we also test RPM, RPKM, and TPM. The deconvolution performance is assessed by RMSE and Spearman's correlation coefficient. Note the results by SQUID are discussed in the text only. (B,C) Boxplots showing the deconvolution performance of eight deconvolution methods for the bulk and pseudobulk data. RMSE and Spearman's correlation coefficient values are calculated across seven major cell types for each sample, with gray denoting pseudobulk and red denoting real bulk. Smaller RMSEs or larger Spearman's correlations indicate a higher accuracy in proportion estimation. (D) Boxplots showing the distributions of deconvolution estimates at the cell type level for all 24 retinal samples. Each color corresponds to a given deconvolution method, with black denoting the ground truth, and each panel corresponds to a given cell type. (E,F), An overview of deconvolution performance at the cell type level across the eight methods using RMSE and Spearman's correlation coefficient, respectively. Lighter colors correspond to lower RMSE or Spearman's correlation coefficient values. Gray indicates NA.

Figure 4.

Using DeMixSC to deconvolve a large cohort of human peripheral retinal samples. (A) PCA plots of both the retina cohort data and the benchmark data. Red denotes the bulk data to be deconvolved; blue denotes the benchmark bulk data; and green denotes the benchmark pseudobulk data. (B,C) Panels demonstrating the robustness of DeMixSC to different reference matrices at both the cell type and sample levels. Higher correlation coefficients indicate better performance. (D) Distributions of DeMixSC estimated cell type proportions of Ratnapriya et al. (2019) data using consensus references. Each panel corresponds to a given cell type. The P-values for Student's t-tests comparing the estimated cell type proportions between non-AMD (healthy) and AMD groups are denoted as follows: (ns) not significant, P-value > 0.05; (*) P-value ≤ 0.05; (**) P-value ≤ 0.01; and (***) P-value ≤ 0.001.

Figure 5.

Using DeMixSC to deconvolve HGSC samples. (A) Boxplots showing the deconvolution performance of eight deconvolution methods for the pseudobulk and three types of bulk data in the HGSC benchmark data set. RMSE values and Spearman's correlation coefficients are calculated across 13 cell types for each sample. Smaller RMSEs or larger Spearman's correlations indicate higher accuracy in proportion estimation. (B) Distributions of DeMixSC estimated cell type proportions of Lee et al. (2020) data using consensus references. Each panel corresponds to a given cell type. (NK cells) natural killer cells, (ILC) innate lymphoid cells, (DC) dendritic cells macrophages, and (pDC) plasmacytoid dendritic cells. The P-values for Student's t-tests comparing the estimated cell type proportions across R0, ER, and PR groups are denoted as follows: (ns) not significant, P-value > 0.05; (*) P-value ≤ 0.05; (**) P-value ≤ 0.01; and (***) P-value ≤ 0.001. (C) Scatter plot comparing DeMixSC estimates of macrophages with immunofluorescent measures (CD68/CD163) in 21 HGSC samples. The black dashed line represents the diagonal, and the gray solid line indicates the linear fit across the data points.