Brandon JewpersName, Marcus AlvarezpersName, Elior RahmanipersName, Zong MiaopersName, Arthur KopersName, Kristina M. GarskepersName, Jae Hoon SulpersName, Kirsi H. PietiläinenpersName, Päivi PajukantapersName, Eran HalperinpersName (2020date). “Accurate estimation of cell composition in bulk expression through robust integration of single-cell information.” Nature Communicationstitle·j.
Object of study — estimation of cell-type proportionsterm·key from bulk RNA-seqterm·key using single-cell references
Research question — how to decompose bulk expression when single-cell profiles break the proportionality assumption
Contribution — Bisqueterm·mth, a fast decompositionterm·key method with gene-specific bulk transformationsterm·mth
phenomenon · claim-type · cross-domain-transfer · mechanism · claim-relation · replication-status · limitation — grounded in the full text| phenomenon | gene-specific technical bias between single-nucleus and bulk RNA-seq — “we found that the majority of gene biases were preserved across individuals, tissues, and experiments ( R = 0.75 across experiments)” | |
| claim-type | method — “These transformations significantly improve decomposition performance compared to existing methods when there is significant technical variation in the generation of the reference profile and observed bulk expression.” | |
| method — evidence | — “Bisque implements a regression-based approach that utilizes single-cell RNA-seq (scRNA-seq) or single-nucleus RNA-seq (snRNA-seq) data to generate a reference expression profile and learn gene-specific bulk expression transformations to robustly decompose RNA-seq data.” | |
| object of study — evidence | — “We present Bisque, a tool for estimating cell type proportions in bulk expression.” | |
| research question — evidence | — “Thus, while single-cell RNA-seq technologies have provided unprecedented resolution in identifying expression profiles of cell types in heterogeneous tissues, these profiles generally may not follow the direct proportionality assumptions of regression-based methods, as we demonstrate here.” | |
| contribution — evidence | — “These transformations significantly improve decomposition performance compared to existing methods when there is significant technical variation in the generation of the reference profile and observed bulk expression.” | |
| mechanism | marker co-variation As a cell type's proportion varies, its marker genes' expression moves with it. — “The basis of most decomposition methods relies on the logic that as the proportion of a cell type varies across individuals, the expression of its marker genes will tend to correlate in the same direction as its cell-type proportion.” | mechanisms[0] |
| sense | decomposition - estimating the cell-type proportions that make up a bulk tissue sample — “We present Bisque, a tool for estimating cell type proportions in bulk expression.” | senses[0] |
| relation | scRNA-seq —same-as→ single-cell RNA-seq — ( ) — “Bisque implements a regression-based approach that utilizes single-cell RNA-seq (scRNA-seq) or single-nucleus RNA-seq (snRNA-seq) data to generate a reference expression profile and learn gene-specific bulk expression transformations to robustly decompose RNA-seq data.” — “utilizes single-cell RNA-seq (scRNA-seq)” | concept_edges[0] |
| weight | scRNA-seq —same-as→ single-cell RNA-seq — 0.95 — “Bisque implements a regression-based approach that utilizes single-cell RNA-seq (scRNA-seq) or single-nucleus RNA-seq (snRNA-seq) data to generate a reference expression profile and learn gene-specific bulk expression transformations to robustly decompose RNA-seq data.” | concept_edges[0] |
| relation | snRNA-seq —same-as→ single-nucleus RNA-seq — ( ) — “Bisque implements a regression-based approach that utilizes single-cell RNA-seq (scRNA-seq) or single-nucleus RNA-seq (snRNA-seq) data to generate a reference expression profile and learn gene-specific bulk expression transformations to robustly decompose RNA-seq data.” — “or single-nucleus RNA-seq (snRNA-seq) data” | concept_edges[1] |
| weight | snRNA-seq —same-as→ single-nucleus RNA-seq — 0.95 — “Bisque implements a regression-based approach that utilizes single-cell RNA-seq (scRNA-seq) or single-nucleus RNA-seq (snRNA-seq) data to generate a reference expression profile and learn gene-specific bulk expression transformations to robustly decompose RNA-seq data.” | concept_edges[1] |
| relation | NNLS —same-as→ non-negative least-squares regression — ( ) — “MuSiC 16 also leverages single-cell expression as a reference, instead using a weighted non-negative least-squares regression (NNLS) model for decomposition” — “non-negative least-squares regression (NNLS)” | concept_edges[2] |
| weight | NNLS —same-as→ non-negative least-squares regression — 0.95 — “MuSiC 16 also leverages single-cell expression as a reference, instead using a weighted non-negative least-squares regression (NNLS) model for decomposition” | concept_edges[2] |
| relation | CIBERSORT —is-a→ SVM-regression-based approach — is a — “CIBERSORT 14 is a SVM-regression-based approach, originally designed for microarray data that utilizes a reference generated from purified cell populations.” — “CIBERSORT 14 is a SVM-regression-based approach” | concept_edges[3] |
| weight | CIBERSORT —is-a→ SVM-regression-based approach — 0.9 — “CIBERSORT 14 is a SVM-regression-based approach, originally designed for microarray data that utilizes a reference generated from purified cell populations.” | concept_edges[3] |
| effect direction | transformations —causes→ decomposition performance — increases — “These transformations significantly improve decomposition performance compared to existing methods when there is significant technical variation in the generation of the reference profile and observed bulk expression.” — “significantly improve decomposition performance” | concept_edges[4] |
| note | transformations —causes→ decomposition performance — Conditional: the paper states the improvement holds when there is significant technical variation between the reference profile and the bulk data. — “These transformations significantly improve decomposition performance compared to existing methods when there is significant technical variation in the generation of the reference profile and observed bulk expression.” | concept_edges[4] |
| relation | transformations —causes→ decomposition performance — improve — “These transformations significantly improve decomposition performance compared to existing methods when there is significant technical variation in the generation of the reference profile and observed bulk expression.” — “These transformations significantly improve decomposition performance” | concept_edges[4] |
| weight | transformations —causes→ decomposition performance — 0.7 — “These transformations significantly improve decomposition performance compared to existing methods when there is significant technical variation in the generation of the reference profile and observed bulk expression.” | concept_edges[4] |
| relation | Bisque —contrasted-with→ CIBERSORT — compared to — “Bisque (4 s) and MuSiC (10 s) estimated cell proportions relatively quickly compared to BSEQ-sc (273 s), CIBERSORT (298 s), and CIBERSORTx (6566 s) (Fig. 5b ).” — “estimated cell proportions relatively quickly compared to BSEQ-sc (273 s), CIBERSORT (298 s)” | concept_edges[5] |
| weight | Bisque —contrasted-with→ CIBERSORT — 0.8 — “Bisque (4 s) and MuSiC (10 s) estimated cell proportions relatively quickly compared to BSEQ-sc (273 s), CIBERSORT (298 s), and CIBERSORTx (6566 s) (Fig. 5b ).” | concept_edges[5] |
| replication | original | |
| limitation | The method assumes cell populations are equally represented in single-cell and bulk data. — “Our reference-based decomposition method is based on the assumption that cell populations are equally represented in single-cell and bulk RNA sequencing of the same tissue samples.” | limitations[0] |
finding · quantity · variable · definition · key-term — grounded in the full text| finding | Gene-level biases between technologies are preserved across individuals, tissues and experiments. — “we found that the majority of gene biases were preserved across individuals, tissues, and experiments ( R = 0.75 across experiments)” | content_tags.findings[0] |
| finding | A core set of marker genes suffices; less important markers matter little. — “These results suggest that as long as a core set of marker genes are present, variations in less important marker genes will have little effect on downstream analyses.” | content_tags.findings[1] |
| finding | CIBERSORT depends on sorting cells to build its reference. — “A major limitation of this approach is the reliance on sorting cells to estimate a reference gene expression panel.” | content_tags.findings[2] |
| quantity | 0.30 — — “Specifically, we observed best fit lines relating these expression levels between technologies with a mean slope of roughly 0.30 and a variance in slope of 5.67.” | content_tags.quantities[0] |
| quantity | 4 s — — “Bisque (4 s) and MuSiC (10 s) estimated cell proportions relatively quickly compared to BSEQ-sc (273 s), CIBERSORT (298 s), and CIBERSORTx (6566 s) (Fig. 5b ).” | content_tags.quantities[1] |
| quantity | R = 0.75 — — “we found that the majority of gene biases were preserved across individuals, tissues, and experiments ( R = 0.75 across experiments)” | content_tags.quantities[2] |
| variable | decomposition performance — “These transformations significantly improve decomposition performance compared to existing methods when there is significant technical variation in the generation of the reference profile and observed bulk expression.” | content_tags.variables[0] |
| definition | Bisque — “We present Bisque, a tool for estimating cell type proportions in bulk expression.” | content_tags.definitions[0] |
| key-terms | Bisque, decomposition, cell-type proportions, bulk RNA-seq, snRNA-seq, reference profile, marker genes, CIBERSORT, MuSiC, regression |
| concept | scRNA-seq — single-cell RNA-seq — “Bisque implements a regression-based approach that utilizes single-cell RNA-seq (scRNA-seq) or single-nucleus RNA-seq (snRNA-seq) data to generate a reference expression profile and learn gene-specific bulk expression transformations to robustly decompose RNA-seq data.” | concepts[0] |
| concept | single-cell RNA-seq — not defined in this paper's own words; sequencing the transcriptomes of individual cells — “Bisque implements a regression-based approach that utilizes single-cell RNA-seq (scRNA-seq) or single-nucleus RNA-seq (snRNA-seq) data to generate a reference expression profile and learn gene-specific bulk expression transformations to robustly decompose RNA-seq data.” | concepts[1] |
| concept | snRNA-seq — single-nucleus RNA-seq — “Bisque implements a regression-based approach that utilizes single-cell RNA-seq (scRNA-seq) or single-nucleus RNA-seq (snRNA-seq) data to generate a reference expression profile and learn gene-specific bulk expression transformations to robustly decompose RNA-seq data.” | concepts[2] |
| concept | single-nucleus RNA-seq — not defined in this paper's own words; sequencing nuclear RNA from individual nuclei — “Bisque implements a regression-based approach that utilizes single-cell RNA-seq (scRNA-seq) or single-nucleus RNA-seq (snRNA-seq) data to generate a reference expression profile and learn gene-specific bulk expression transformations to robustly decompose RNA-seq data.” | concepts[3] |
| concept | NNLS — non-negative least-squares regression — “MuSiC 16 also leverages single-cell expression as a reference, instead using a weighted non-negative least-squares regression (NNLS) model for decomposition” | concepts[4] |
| concept | non-negative least-squares regression — not defined in this paper's own words; regression constrained to non-negative coefficients — “MuSiC 16 also leverages single-cell expression as a reference, instead using a weighted non-negative least-squares regression (NNLS) model for decomposition” | concepts[5] |
| concept | CIBERSORT — a SVM-regression-based approach, originally designed for microarray data that utilizes a reference generated from purified cell populations — “CIBERSORT 14 is a SVM-regression-based approach, originally designed for microarray data that utilizes a reference generated from purified cell populations.” | concepts[6] |
| concept | SVM-regression-based approach — not defined in this paper's own words; decomposition by support-vector regression — “CIBERSORT 14 is a SVM-regression-based approach, originally designed for microarray data that utilizes a reference generated from purified cell populations.” | concepts[7] |
| concept | Bisque — a tool for estimating cell type proportions in bulk expression — “We present Bisque, a tool for estimating cell type proportions in bulk expression.” | concepts[8] |
| concept | transformations — not defined in this paper's own words; gene-specific rescaling of bulk expression to match the single-cell reference — “These transformations significantly improve decomposition performance compared to existing methods when there is significant technical variation in the generation of the reference profile and observed bulk expression.” | concepts[9] |
| concept | decomposition performance — not defined in this paper's own words; accuracy of estimated cell-type proportions — “These transformations significantly improve decomposition performance compared to existing methods when there is significant technical variation in the generation of the reference profile and observed bulk expression.” | concepts[10] |
Full text — tagged
The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.
Accurate estimation of cell composition in bulk expression through robust integration of single-cell information Abstract We present Bisqueterm·mth, a tool for estimating cell type proportionsterm·key in bulk expression. Bisqueterm·mth implements a regression-based approach that utilizes single-cell RNA-seqterm·mth (scRNA-seqterm·mth) or single-nucleus RNA-seqterm·mth (snRNA-seqterm·mth) data to generate a reference expression profile and learn gene-specific bulk expression transformationsterm·mth to robustly decompose RNA-seq data. These transformationsterm·mth significantly improve decomposition performanceterm·mth compared to existing methods when there is significant technical variation in the generation of the reference profileterm·key and observed bulk expression. Importantly, compared to existing methods, our approach is extremely efficient, making it suitable for the analysis of large genomic datasets that are becoming ubiquitous. When applied to subcutaneous adipose and dorsolateral prefrontal cortex expression datasets with both bulk RNA-seqterm·key and snRNA-seqterm·mth data, Bisqueterm·mth replicates previously reported associations between cell type proportionsterm·key and measured phenotypes across abundant and rare cell types. We further propose an additional mode of operation that merely requires a set of known marker genesterm·key.
Introduction Bulk RNA-seqterm·key experiments typically measure total gene expression from heterogeneous tissues, such as tumor and blood samples 1 , 2 . Variability in cell-type composition can significantly confound analyses of these data, such as in identification of expression quantitative trait loci (eQTLs) or differentially expressed genes 3 . Cell-type heterogeneity may also be of interest in profiling changes in tissue composition associated with disease, such as cancer 4 or diabetes 5 . In addition, measures of cell composition can be leveraged to identify cell-specific eQTLs 6 , 7 or differential expression 6 from bulk data.
Traditional methods for determining cell-type composition, such as immunohistochemistry or flow cytometry, rely on a limited set of molecular markers and lack in scalability relative to the current rate of data generation 8 . Single-cell technologies provide a high-resolution view into cellular heterogeneity and cell-type-specific expression 9 – 11 . However, these experiments remain costly and noisy compared to bulk RNA-seqterm·key 12 . Collection of bulk expression data remains an attractive approach for identifying population-level associations, such as differential expression regardless of cell-type specificity. Moreover, many bulk RNA-seqterm·key studies that have been performed in recent years resulted in a large body of data that is available public databases such as dbGAP and GEO. Given the wide availability of these bulk data, the estimation of cell-type proportionsterm·key, often termed decompositionterm·key, can be used to extract large-scale cell-type-specific information.
There exist a number of methods for decomposing bulk expression, many of which are regression-based and leverage cell-type-specific expression data as a reference profileterm·key 13 . CIBERSORTterm·mth 14 is a SVM-regression-based approachterm·mth, originally designed for microarray data that utilizes a reference generated from purified cell populations. A major limitation of this approach is the reliance on sorting cells to estimate a reference gene expression panel. BSEQ-sc 15 instead generates a reference profileterm·key from single-cell expression data that is used in the CIBERSORTterm·mth model. MuSiCterm·key 16 also leverages single-cell expression as a reference, instead using a weighted non-negative least-squares regressionterm·mth (NNLSterm·mth) model for decompositionterm·key, with improved performance over BSEQ-sc in several datasets.
The distinct nature of the technologies used to generate bulk and single-cell sequencing data may present an issue for decompositionterm·key models that assume a direct proportional relationship between the single-cell-based reference and observed bulk mixture. For example, the capture of mRNA and chemistry of library preparation can differ significantly between bulk tissue and single-cell RNA-seqterm·mth methods, as well as between different single-cell technologies 17 , 18 . Moreover, some technologies may be measuring different parts of the transcriptome, such as nuclear pre-mRNA in single-nucleus RNA-seqterm·mth (snRNA-seqterm·mth) experiments as opposed to cellular and extra-cellular mRNA observed in traditional bulk RNA-seqterm·key experiments. As we show later, these differences may introduce gene-specific biases that break down the correlation between cell-type-specific and bulk tissue measurements. Thus, while single-cell RNA-seqterm·mth technologies have provided unprecedented resolution in identifying expression profiles of cell types in heterogeneous tissues, these profiles generally may not follow the direct proportionality assumptions of regression-based methods, as we demonstrate here.
We present Bisqueterm·mth, a highly efficient tool to measure cellular heterogeneity in bulk expression through robust integration of single-cell information, accounting for biases introduced in the single-cell sequencing protocols. The goal of Bisqueterm·mth is to integrate the different chemistries/technologies of single-cell and bulk tissue RNA-seq to estimate cell-type proportionsterm·key from tissue-level gene expression measurements across a larger set of samples. Our reference-based model decomposes bulk samples using a single-cell-based reference profileterm·key and, while not required, can leverage single-cell and bulk measurements for the same samples for further improved decompositionterm·key accuracy. This approach employs gene-specific transformationsterm·mth of bulk expression to account for biases in sequencing technologies as described above. When a reference profileterm·key is not available, we propose BisqueMarker, a semi-supervised model that extracts trends in cellular composition from normalized bulk expression samples using only cell-specific marker genesterm·key that could be obtained using single-cell data. We demonstrate using simulated and real datasets from brain and adipose tissue that our method is significantly more accurate than existing methods. Furthermore, it is extremely efficient, requiring seconds in cases where other methods require hours; thus, it is scalable to large genomic datasets that are now becoming available.
Results Method overview (Bisqueterm·mth) A graphical overview of Bisqueterm·mth is presented in Fig. 1 . Our reference-based decompositionterm·key model requires bulk RNA-seqterm·key counts data and a reference dataset with read counts from single-cell RNA-seqterm·mth. In addition, the single-cell data should be labeled with cell types to be quantified. A reference profileterm·key is generated by averaging read count abundances within each cell type in the single-cell data. Given the reference profileterm·key and cell proportions observed in the single-cell data, our method learns gene-specific transformationsterm·mth of the bulk data to account for technical biases between the sequencing technologies. Bisqueterm·mth can then estimate cell proportions from the bulk RNA-seqterm·key data using the reference and the transformed bulk expression data using non-negative least-squares (NNLSterm·mth) regressionterm·key. Fig. 1 Graphical overview of the Bisqueterm·mth decompositionterm·key method. We integrate single-cell and bulk expression by learning gene-specific bulk transformationsterm·mth (pictured on right ) that align the two datasets for accurate decompositionterm·key.
Evaluation of decomposition performanceterm·mth in adipose tissue We applied our method to 106 bulk RNA-seqterm·key subcutaneous adipose tissue samples collected from both lean and obese individuals, where 6 samples have both bulk RNA-seqterm·key and snRNA-seqterm·mth data available (Table 1 ). Each of the participants gave a written informed consent. The study protocol was approved by the Ethics Committee at the Helsinki University Hospital, Helsinki, Finland. Adipose tissue consists of several cell types, including adipocytes that are expected to be the most abundant population. Adipose tissue also contains structural cell types (i.e. fibroblasts and endothelial cells) and immune cells (i.e. macrophages and T cells) 19 . These 5 cell-type populations were identified from the snRNA-seqterm·mth data (Supplementary Fig. 1a ). Table 1 Summary of snRNA-seqterm·mth and bulk expression datasets used for benchmarking Bisqueterm·mth and existing methods. Tissue Number of samples Bulk RNA-seqterm·key platform snRNA-seqterm·mth platform snRNA-seqterm·mth samples Total nuclei Average nuclei per individual Number of cell types Subcutaneous adipose 106 Illumina NovaSeq 10x Genomics Chromium 6 10,947 1824 5 Dorsolateral prefrontal cortex 636 Illumina HiSeq 10x Genomics Chromium 8 68,028 8503 11 We observed significant biases between the snRNA-seqterm·mth and bulk RNA-seqterm·key data in samples that had both data available. We found that the linear relationship between the pseudo-bulk (summed snRNA-seqterm·mth reads across cells) and the true bulk expression varied significantly by each gene (Fig. 2a ). Specifically, we observed best fit lines relating these expression levels between technologies with a mean slope of roughly 0.30 and a variance in slope of 5.67. In our model, a slope of 1 would indicate no bias between technologies. We further investigated whether gene expression differences between the bulk and snRNA-seqterm·mth were the same across individuals and experiments. Comparing log-ratios of RNA-seq to snRNA-seqterm·mth expression levels, we found that the majority of gene biases were preserved across individuals, tissues, and experiments ( R = 0.75 across experiments) (Supplementary Fig. 3 ), providing evidence that technological differences drive consistent gene expression differences across bulk and snRNA-seqterm·mth methods. Fig. 2 The effect of discrepancies between a single-cell-based reference and bulk expression on decompositionterm·key. a Observed discrepancies in real data between single-nucleus and bulk expression for selected marker genesterm·key ( left ) for six individuals. Each color corresponds to a gene. On the left, observed bulk expression on the x axis is plotted against the pseudo-bulk expression on the y axis, where pseudo-bulk expression is calculated by summing the single-cell-based reference with cell proportions as weights. On the right, the Bisqueterm·mth transformationterm·mth of bulk expression is on the x axis. Bisqueterm·mth recovers a one-to-one relationship by transforming the bulk expression for improved decompositionterm·key accuracy ( right ). b Simulation of bulk expression for six individuals based on true adipose snRNA-seqterm·mth data with increasing gene-specific differences. These differences are modeled as a linear transformationterm·mth of the summed snRNA-seqterm·mth counts with coefficient and intercept sampled from half-normal distributions with parameter as indicated on the x axis. At σ = 0, the simulated bulk is simply the sum of the observed single-cell read counts. Performance on y axis measured in global Pearson correlation ( R ) ( left ) and root-mean squared deviation (RMSD) ( right ). Shaded regions indicate 95% confidence intervals based on bootstrapping with central lines indicating the mean observed value. Bisqueterm·mth remains robust to increasing gene-specific variation between single-cell and bulk expression levels. Source data are provided as a Source Data file.
We performed simulations based on the adipose snRNA-seqterm·mth data to demonstrate the effect of technology-based biases between the reference profileterm·key and bulk expression on decomposition performanceterm·mth. In these analyses, we benchmarked Bisqueterm·mth and three existing methods (MuSiCterm·key, BSEQ-sc, and CIBERSORTterm·mth). Briefly, we simulated bulk expression for 6 individuals by summing the observed snRNA-seqterm·mth read counts. To model discordance between the reference and bulk, we applied gene-specific linear transformationsterm·mth of the simulated bulk expression. For each gene, the coefficient and intercept of the linear transformationterm·mth were sampled from half-normal distributions with increasing variance. In this model, a higher variance corresponds to a larger bias between sequencing experiments. Although these transformationsterm·mth closely mirrored the Bisqueterm·mth decompositionterm·key model, they utilized the true snRNA-seqterm·mth counts for each individual whereas Bisqueterm·mth learned these transformationsterm·mth using the reference profileterm·key generated from averaging these counts across all cells. Hence, this simulation framework introduced additional noise that Bisqueterm·mth does not entirely model. We evaluated decomposition performanceterm·mth by comparing proportion estimates to the proportions observed in the snRNA-seqterm·mth data in terms of global Pearson correlation ( R ) and root-mean squared deviation (RMSD). Owing to the small number of samples, we applied leave-one-out cross-validation to predict the cell composition of each individual using the remaining snRNA-seqterm·mth samples as training data for each method. In these simulations, Bisqueterm·mth remained robust ( R ≈ 0.85, RMSD ≈ 0.07) at higher levels of simulated bias between the bulk and snRNA-seq-based reference (Fig. 2b ).
Next, we performed this cross-validation benchmark on the observed bulk RNA-seqterm·key data for these 6 individuals and found that Bisqueterm·mth ( R = 0.923, RMSD = 0.074) provided significantly improved global accuracy in detecting each cell type over existing methods (Table 2 , Supplementary Fig. 1b ). MuSiCterm·key ( R = −0.111, RMSD = 0.427), BSEQ-sc ( R = −0.113, RMSD = 0.432), and CIBERSORTterm·mth ( R = −0.131, RMSD = 0.416) severely underestimated the proportion of adipocytes (the most abundant population in adipose tissue) while overestimating the endothelial cell fraction. We also benchmarked CIBERSORTx 20 , which employs a batch correction mode to account for biases in sequencing technologies. Although CIBERSORTx ( R = 0.687, RMSD = 0.099) outperformed existing methods, Bisqueterm·mth provided improved accuracy. It should be noted that cell-specific accuracy is more informative than global R and RMSD; however, these small sample sizes did not provide robust measures of within-cell-type performance in this cross-validation framework (Supplementary Fig. 1c ). We were able to slightly improve the number of detected cell populations by MuSiCterm·key, BSEQ-sc, and CIBERSORTterm·mth when we considered only snRNA-seqterm·mth reads aligning to exonic regions of the transcriptome, indicating that intronic reads introduced increasing discrepancy between snRNA-seqterm·mth and bulk RNA-seqterm·key in the context of decompositionterm·key. However, given that a significant portion of the nuclear transcriptome consists of pre-mRNA, this filtering process removed over 40% of cells detected in the snRNA-seqterm·mth data. Moreover, Bisqueterm·mth provided improved accuracy over existing methods using this exonic subset of the snRNA-seqterm·mth data (Supplementary Fig. 1d ). Table 2 Leave-one-out cross-validation in subcutaneous adipose using 6 samples with snRNA-seqterm·mth and bulk RNA-seqterm·key data available. Method R RMSD Bisqueterm·mth 0.923 ± 0.064 0.074 ± 0.034 CIBERSORTx 0.687 ± 0.450 0.099 ± 0.046 MuSiCterm·key −0.111 ± 0.182 0.427 ± 0.058 BSEQ-sc −0.113 ± 0.180 0.432 ± 0.058 CIBERSORTterm·mth −0.131 ± 0.176 0.416 ± 0.059 Proportions based on snRNA-seqterm·mth were used as a proxy for the true proportions. Performance measured in Pearson correlation ( R ) and root-mean-square deviation (RMSD) across all 5 identified cell types in each sample. Reported values were averaged across the 6 samples with standard deviation indicated. Bold values indicate the highest performing method with respect to each metric. Source data are provided as a Source Data file.
We then applied these decompositionterm·key methods to the remaining 100 bulk samples and found that the distribution of cell-proportion estimates produced by Bisqueterm·mth were most concordant with the expected distribution inferred from the limited number of snRNA-seqterm·mth samples and previously reported proportions 21 , 22 (Fig. 3a ). Although these benchmarks provided a measure of calibration (i.e. the ability to detect cell populations in expected ranges), they did not provide measurements of cell-specific proportion accuracy across individuals. In order to evaluate cell-specific accuracy, we replicated previously reported associations between cell proportions and measured phenotypes. Specifically, we compared cell-proportion estimates from each method to body mass index (BMI) and Matsuda index, a measure of insulin resistance. We measured the significance of these association based on t -values estimated in a linear-mixed model accounting for age, age-squared, and sex as fixed effects and relatedness as a random effect. Fig. 3 Decompositionterm·key benchmark in human subcutaneous adipose tissue. a Comparison of decompositionterm·key estimates from 100 individuals with estimates from 6 individuals with snRNA-seqterm·mth data available. Each color represents a benchmarked method. Boxes indicate the quartiles of the estimated proportions with whiskers extending 1.5 times the interquartile range. Points are individual samples that are represented by the boxplot. b , c Scatterplots comparing decompositionterm·key estimates with measured phenotypes in 100 individuals. Reported ‘rho’ corresponds to Spearman correlation and p -values indicate the significance of these correlations, with an asterisk denoting significance after correction for covariates in a linear-mixed model. CIBERSORTterm·mth and BSEQ-sc are not shown as they did not detect these cell populations. These examples include the most abundant (adipocytes) and rarest (T cells) cell types identified in the snRNA-seqterm·mth data. Significance of associations reported in Supplementary Table 1 . b Adipocyte proportion has been observed to negatively correlate with BMI so we expected a negative correlation. Bisqueterm·mth ( p = 0.030) and CIBERSORTx ( p = 0.001) produced significant negative associations after correcting for sex, age, age-squared, and relatedness in a linear-mixed model. c T-cell proportion has previously been reported to positively correlate with insulin resistance. Matsuda index decreases with higher insulin resistance so we expected a negative correlation. Bisqueterm·mth ( p = 0.002) and CIBERSORTx ( p = 0.046) produced significant negative associations after correcting for diabetes status, sex, age, age-squared, and relatedness in a linear-mixed model. Source data are provided as a Source Data file.
Obesity is associated with adipocyte hypertrophy, the expansion of the volume of fat cells 23 ; thus, we expected a negative association between adipocyte proportion and BMI. Bisqueterm·mth, MuSiCterm·key, and CIBERSORTx produced adipocyte proportion estimates that replicate this behavior, while BSEQ-sc and CIBERSORTterm·mth were unable to detect this cell population (Fig. 3b ). The adipocyte proportion estimates produced by Bisqueterm·mth ( p = 0.030) and CIBERSORTx ( p = 0.001) had a significant negative association with BMI (Supplementary Table 1a ). In addition, macrophage abundance has been shown to increase in adipose tissue with higher levels of obesity, concomitant with a state of low grade inflammation 24 . Each method detected macrophage populations that positively associated with BMI; however, only Bisqueterm·mth ( p < 0.001), BSEQ-sc ( p = 0.004), and CIBERSORTx ( p = 0.034) reached significance (Supplementary Table 1b ).
T cells were the least abundant cell-type population identified from the snRNA-seqterm·mth data, constituting around 4% of all sequenced nuclei. The abundance of T cells has been observed to positively correlate with insulin resistance 25 . Thus, we compared decompositionterm·key estimates for T-cell proportions to Matsuda index. As a lower Matsuda index indicates higher insulin resistance, we expect a negative association between T-cell proportion and Matsuda index. Proportion estimates produced by Bisqueterm·mth and CIBERSORTx followed this trend while the remaining existing methods did not identify T cells in the bulk samples (Fig. 3c ). We found this association significant for Bisqueterm·mth ( p = 0.002) and CIBERSORTx ( p = 0.046) (Supplementary Table 1c ) after correcting for diabetes status, as Matsuda index may not be informative in these individuals 26 .
Evaluation of decomposition performanceterm·mth in cortex tissue We also benchmarked these decompositionterm·key methods using expression data collected from the dorsolateral prefrontal cortex (DLPFC). This dataset was generated by the Rush Alzheimer’s Disease (AD) Center 27 and includes 636 postmortem bulk RNA-seqterm·key samples. The Religious Orders Study and Rush Memory and Aging Project were approved by an IRB of Rush University Medical Center. Both bulk RNA-seqterm·key and snRNA-seqterm·mth data were collected from 8 of the individuals (Table 1 ). Using the same pipeline we used to process the adipose dataset, we identified 11 clusters: 3 neuronal subtypes, 2 interneuronal subtypes, 2 astrocyte subtypes, oligodendrocytes, oligodendrocyte progenitor cells, and microglia (Supplementary Fig. 2a ). We observed a higher overlap in marker genesterm·key for these clusters than in those identified in the adipose dataset (average of 10% of marker genesterm·key shared between clusters in DLPFC compared to 3% in adipose) (Supplementary Fig. 4a, b ).
We again applied leave-one-out cross-validation on the 8 individuals with both RNA-seq and snRNA-seqterm·mth data available. In this example, we randomly sampled 25% of the nuclei in the snRNA-seqterm·mth data to accommodate CIBERSORTx (which is currently web-based and restricts the size of files that can be processed). Bisqueterm·mth was able to detect each cell population identified from the snRNA-seqterm·mth data with high global accuracy ( R = 0.924, RMSD = 0.029) while MuSiCterm·key ( R = −0.192, RMSD = 0.173), BSEQ-sc ( R = 0.098, RMSD = 0.120), and CIBERSORTterm·mth ( R = −0.281, RMSD = 0.197) did not detect a number of cell populations (Table 3 , Supplementary Fig. 2b , c ). Bisqueterm·mth also provided higher accuracy than CIBERSORTx ( R = 0.671, RMSD = 0.070). However, we found that the performance of the existing methods improved when estimates with subtypes were summed together (Supplementary Fig. 2d ). Although each method was able to quantify major cell populations after merging subtypes, Bisqueterm·mth was able to distinguish between these closely related cell populations. Interestingly, we found that in both adipose and DLPFC, endothelial cell proportions were overestimated by each of the existing methods. Table 3 Leave-one-out cross-validation in dorsolateral prefrontal cortex using 8 samples with snRNA-seqterm·mth and bulk RNA-seqterm·key data available. Method R RMSD Bisqueterm·mth 0.924 ± 0.062 0.029 ± 0.010 CIBERSORTx 0.671 ± 0.153 0.070 ± 0.019 MuSiCterm·key −0.192 ± 0.107 0.173 ± 0.013 BSEQ-sc 0.098 ± 0.216 0.120 ± 0.023 CIBERSORTterm·mth −0.281 ± 0.049 0.197 ± 0.012 Proportions based on snRNA-seqterm·mth were used as a proxy for the true proportions. Performance measured in Pearson correlation ( R ) and root-mean-square deviation across all 11 identified cell types in each sample. Reported values were averaged across the 8 samples with standard deviation indicated. We performed these experiments with 25% of the snRNA-seqterm·mth data in order to accommodate the file size limit of the current web-based implementation of CIBERSORTx. Bold values indicate the highest performing method with respect to each metric. Source data are provided as a Source Data file.
We applied these decompositionterm·key methods to the remaining 628 individuals and compared the distribution of estimates to the proportions observed in the 8 snRNA-seqterm·mth samples. We found that Bisqueterm·mth was able to detect each cell population and produced estimates that were closest in mean to the snRNA-seqterm·mth observations (Fig. 4a ). The increased accuracy of Bisqueterm·mth over existing methods persisted when we merged closely related subtypes (Supplementary Fig. 2e ). Moreover, immunohistochemistry (IHC) analyses on a 70 of these samples found similar proportions of major cell populations 28 , confirming the relative accuracy of snRNA-seq-based estimates of cell proportions. Fig. 4 Decompositionterm·key benchmark in human dorsolateral prefrontal cortex tissue. We randomly sampled 25% of the nuclei in the snRNA-seqterm·mth data to accommodate the file size limit of the web-based implementation of CIBERSORTx at the time of writing. a Comparison of decompositionterm·key estimates from 628 individuals with estimates from 8 individuals with snRNA-seqterm·mth data available. Each color represents a benchmarked method. Boxes indicate the quartiles of the estimated proportions with whiskers extending 1.5 times the interquartile range. Points are individual samples that are represented by the boxplot. b , c Violin plots depicting association of decompositionterm·key estimates aggregated into major cell types with measured phenotypes in 628 individuals. Reported ‘rho’ corresponds to Spearman correlation and p -values indicate the significance of these correlations, with an asterisk denoting both an expected effect direction and significance after correction for covariates. Examples shown are for the most abundant (neurons) and least abundant (microglia) populations detected in the snRNA-seqterm·mth data. Significance of associations reported in Supplementary Table 2 . b Neuronal degeneration has been observed in patients diagnosed with Alzheimer’s disease (AD). Cognitive diagnostic category measures a physician’s diagnosis of cognitive impairment (CI), with 0 indicating no CI and 4 indicating a confident AD diagnosis. We expected a negative correlation between neuron proportion and cognitive diagnostic category. c Microglia proportion has been observed to positively correlate with increased severity of AD symptoms, such as neurofibrillary tangles. Braak stage provides a semiquantitative measure of tangle severity, so we expected an overall positive correlation between microglia proportion and Braak stage. In addition, a decrease in microglia abundance has been previously reported at Braak stages 5 through 6 in AD patients. Only Bisqueterm·mth produced estimates with a significant positive association ( p = 0.001) after correcting for sex, age, and age-squared in a linear regressionterm·key model. Source data are provided as a Source Data file.
Again, to determine cell-specific decompositionterm·key accuracy, we replicated known associations between cell-type proportionsterm·key and measured phenotypes in the 628 individuals. For these analyses, we compared cell-proportion estimates to each individual’s Braak stage and physician cognitive diagnostic category at time of death. Braak stage is a semiquantitative measure of neurofibrillary tangles, ranging in value from 0 to 5 with increasing severity. The cognitive diagnostic category provides a semiquantitative measure of dementia severity, where a code of 1 indicates no cognitive impairment and 5 indicates a confident diagnosis of AD by physicians. We determined the significance of these associations based on t -values estimated by a linear regressionterm·key model that accounted for age, age-squared, and sex.
Neuronal death is a hallmark symptom of AD 29 . Therefore, we expected to find a negative association between cognitive diagnosis and neuron proportion. We found that each decompositionterm·key method provides estimates of total neuron proportion that tend to decrease with cognitive diagnostic category (Fig. 4b ). Each method generates proportions with negative association with cognitive diagnosis. Each method, excluding BSEQ-sc, reached significance in this model ( p ≤ 0.001 for each method) (Supplementary Table 2a ). As another example, we compared each individual’s Braak stage to their estimated proportion of microglia, a relatively small cell population that constituted roughly 5% of the sequenced nuclei. Microglia activation has been observed to increase with AD severity 30 . We used Braak stage as a proxy for AD severity and expected a positive association between microglia proportion and Braak stage. Bisqueterm·mth and MuSiCterm·key provided estimates that follow this expected trend (Fig. 4c ). Only Bisqueterm·mth produced estimates with a significant positive association ( p = 0.001) (Supplementary Table 2b ). Interestingly, we observed a decrease in microglia proportions estimated by Bisqueterm·mth in Braak stage 6 individuals, which has been previously observed in AD patients 31 .
Runtime comparison of reference-based decompositionterm·key methods Given the large amounts of transcriptomic data that are becoming available, we also benchmarked these decompositionterm·key methods in terms of runtime. In the subcutaneous adipose dataset, which included 100 bulk RNA-seqterm·key samples and 6 snRNA-seqterm·mth samples with about 1800 nuclei sequenced per individual, Bisqueterm·mth was able to estimate cell proportions efficiently compared to existing methods. Bisqueterm·mth (1 s) and MuSiCterm·key (1 s) provided decompositionterm·key estimates faster than BSEQ-sc (26 s), CIBERSORTterm·mth (27 s), and CIBERSORTx (389 s) (Fig. 5a ). Bisqueterm·mth also provided improved efficiency in processing the reduced DLPFC dataset, which included 628 bulk RNA-seqterm·key samples and 8 snRNA-seqterm·mth samples with around 2125 nuclei per individual. Bisqueterm·mth (4 s) and MuSiCterm·key (10 s) estimated cell proportions relatively quickly compared to BSEQ-sc (273 s), CIBERSORTterm·mth (298 s), and CIBERSORTx (6566 s) (Fig. 5b ). Fig. 5 Runtime comparisons in log-transformed seconds for benchmarked reference-based decompositionterm·key methods. a Runtime for subcutaneous adipose dataset, which included 100 RNA-seq samples and 6 snRNA-seqterm·mth samples with around 1800 nuclei per individual. b Runtime for dorsolateral prefrontal cortex dataset, which included 628 RNA-seq samples and 8 snRNA-seqterm·mth samples. We benchmarked each method using around 2125 nuclei per snRNA-seqterm·mth sample. Source data are provided as a Source Data file.
Robustness of the reference-based decompositionterm·key model Our reference-based decompositionterm·key method is based on the assumption that cell populations are equally represented in single-cell and bulk RNA sequencing of the same tissue samples. As this assumption may be violated 32 , we explored the performance of our model as we relaxed this assumption in simulations. First, we simulated snRNA-seqterm·mth data where cell proportions were increasingly biased. Using the DLPFC snRNA-seqterm·mth data, we downsampled or upsampled the cells identified as microglia at varying levels and performed decompositionterm·key. Indeed, the absolute estimates produced by Bisqueterm·mth propagated these shifts in snRNA-seqterm·mth proportions. However, we found that our estimated microglia proportions, regardless of these shifts, maintained an expected positive association with Braak stage. This positive association served as evidence for the correlation between these estimates and the true microglia proportions (Supplementary Fig. 5a ). Given these results, we suggest that users take note of this behavior if both the mean abundances are important for downstream analysis and the single-cell reference data is known to be significantly biased against specific cell populations of interest.
Next, we simulated a situation where an unknown cell population contributes to bulk expression but is not represented in the snRNA-seqterm·mth reference data. For situations where this unknown contribution varies across the bulk dataset, we simulated bulk expression by mixing the observed bulk expression for the DLPFC dataset with increasing amounts of expression observed in the adipose dataset. To determine the effect of unknown cell populations on our model, we analyzed the distribution of residual norms produced by the method. These residual norms provide a measure of the difference between the vector of observed bulk and expression reference weighted by the estimated proportions across all genes for each individual. As we increased the contribution from unknown cell types, the residual norm values tend to increase (Supplementary Fig. 5b ). In our simulation framework, this variability in unknown cell-type contribution could be qualitatively identified by the presence of a multimodal residual norm distribution.
Given that single-cell datasets still remain relatively small compared to bulk datasets, we also explored the impact of sample size in the reference single-cell data on the performance of Bisqueterm·mth. In the DLPFC dataset, we saw a drop in performance when using less than four randomly selected snRNA-seqterm·mth samples (Supplementary Fig. 5c ). This threshold is likely to differ between experiments, though we recommend at least three single-cell samples to generate reference data.
Finally, as marker geneterm·key selection can vary between studies, we were interested in the performance of Bisqueterm·mth as we varied the number of marker genesterm·key. Again, we measured cell-type proportionterm·key estimation performance for microglia in the DLPFC dataset by correlating the estimates with Braak stage, which is known to have a positive association. We recalculated this correlation as we removed marker genesterm·key for this cell type. We removed marker genesterm·key in order of both decreasing and increasing log-fold change, which provides a measure of the importance of marker genesterm·key for identifying this cell type. In both procedures, we observe that as we remove an increasing percentage of the 102 identified marker genesterm·key, performance remains stable until a shared drop off point around 75% (Supplementary Fig. 5d ). As we observed this trend in both marker geneterm·key removal schemes, we assume that a relatively few number of marker genesterm·key, regardless of their log-fold change magnitude, can be used to accurately estimate cell-type proportionsterm·key. These results suggest that as long as a core set of marker genesterm·key are present, variations in less important marker genesterm·key will have little effect on downstream analyses.
Marker-based decompositionterm·key using cell-type marker genesterm·key Although a reference profileterm·key from snRNA-seqterm·mth can help to decompose bulk-level gene expression, it may not be available for the same dataset. The majority of bulk RNA-seqterm·key datasets do not have corresponding snRNA-seqterm·mth data in the same set of individuals. However, marker geneterm·key information from prior experiments can still be applied to distinct expression datasets of the same tissue. The basis of most decompositionterm·key methods relies on the logic that as the proportion of a cell type varies across individuals, the expression of its marker genesterm·key will tend to correlate in the same direction as its cell-type proportionterm·key. This linear co-variation can be captured in a principal components analysis (PCA). Under the same argument, the more cell-type-specific a marker geneterm·key is, the more its expression will reflect its cell-type proportionterm·key. These observations form the basis for BisqueMarker, a weighted PCA-based (wPCA) decompositionterm·key approach. Genes that are more specifically expressed within a cell type will provide more information than genes with shared expression across cell types. To estimate cell-type proportionsterm·key without the use of cell-type-specific gene expression information, we applied wPCA to bulk-level adipose tissue expression.
For each cell type, we extracted the first PC from a wPCA of the expression matrix of its markers. The expression matrix was corrected for the first global expression PC as a covariate so that wPCA estimates would not reflect technical variation. We first confirmed that these genes were distinct across cell types. If 2 cell types share a high proportion of marker genesterm·key, the wPCA estimates from bulk RNA-seqterm·key will correlate highly. We then investigated whether the second or third PC could have represented cell-type proportionsterm·key. The percent of variance explained by the first PC was typically 30–60% across adipose cell types, and additionally, over 90% of the markers correlated in the same direction as the first PC. In contrast, roughly 50–70% of markers correlated in the same direction as the second or third PC. As performed for reference-based decompositionterm·key, we correlated phenotypes with cell-type proportionsterm·key estimated by BisqueMarker. We identified the same associations as with reference-based decompositionterm·key, demonstrating its validity when a reference is not available (Supplementary Table 1 ). Similarly, we observed the same trends between estimated cell-type abundances and phenotypes as we did using our reference-based method in the DLPFC cohort (Supplementary Table 2 ).
Discussion Bisqueterm·mth effectively leverages single-cell information to decompose bulk expression samples, outperforming existing methods in datasets with snRNA-seqterm·mth data available. In simulations, we demonstrated that the decompositionterm·key accuracy of Bisqueterm·mth is robust to increasing variation between the generation of the reference profileterm·key and bulk expression, which is a significant issue when comparing snRNA-seqterm·mth and bulk RNA-seqterm·key data. In observed bulk expression, our reference-based method accurately estimates cell proportions that are consistent with previously reported distributions and reliably detects rare cell types. We found that these estimates consistently follow expected trends with measured phenotypes, suggesting that cell-specific estimates of proportion are sufficiently accurate to extract relevant biological signals. In addition, differences in tissue structure can lead to significant differences in the quality of single-cell expression data 33 . We demonstrated the improved performance of our method in adipose and DLPFC, two distinct tissues, suggesting that Bisqueterm·mth is robust across different tissue types.
The cell-type proportionterm·key estimates determined by Bisqueterm·mth may be utilized to effectively identify cell-type-specific interactions, such as expression quantitative trait loci (eQTLs), and adjust for confounding effects from variability in cell populations. With this reference-based approach, single-cell sequencing of a subset of samples from large-scale bulk expression cohorts can provide high power to detect cell-specific associations in complex phenotypes and diseases.
However, we note that there are limitations to this reference-based method that users should consider. First, if the number of individuals with single-cell data available is small, the reference profileterm·key and gene-specific transformationsterm·mth may become unreliable. In addition, a key assumption of our transformationterm·mth framework is that single-cell-based estimates of cell proportions accurately reflect the true proportions we wish to estimate. As a result of this assumption, Bisqueterm·mth provides estimates of cell proportions reported by the single-cell technology used to generate the reference data. Given that snRNA-seqterm·mth can provide less bias in isolating specific cell types compared to scRNA-seqterm·mth 34 , 35 , we expect these estimates to be useful for downstream analyses such as those previously discussed. Nevertheless, the accuracy of Bisqueterm·mth may decrease if the proportion of cell types captured by single-cell experiments differs significantly from the true physiological distributions. Therefore, we advise users to take caution if there is a known significant bias in the single-cell measurements of a tissue, such as severe underrepresentation of a cell type of interest 32 , 35 , that can affect downstream analysis. Our results demonstrate that even with these limitations, Bisqueterm·mth can be used to provide cell-type specific biological insight in relevant datasets.
In cases where these described issues may be significant, BisqueMarker provides cell-type abundance estimations using only known marker genesterm·key. Although this reference-free method may be less accurate than reference-based methods, it does not depend on single-cell based estimates of cell proportions or expression profiles, but rather on the fact that the expression in certain genes differs across different cell types; moreover, this method also does not model explicitly the expression level, and it is thus robust to biases in the single-cell sequencing protocol. We found that BisqueMarker estimates followed expected trends with measured phenotypes; however, it should be noted that this method estimates relative differences in abundances that cannot be compared across cell types. Also, given the semi-supervised nature of this method, these cell-type abundance estimates may include signals from technical or other biological variation in the data. Therefore, we highly suggest applying this method to data that is properly normalized with sources of undesired variation removed.
Methods Processing bulk expression data Paired-end reads were aligned with STAR v2.5.1 using default options. Gene counts were quantified using featureCounts v1.6.3. For featureCounts, fragments were counted at the gene-name level. Alignment and gene counts were generated against the GRCh38.p12 genome assembly. STAR v2.5.1 and GRCh38.p12 were included with CellRanger 3.0.2, which was used to process the single-nucleus data.
Processing single-nucleus expression data Reads from single nuclei sequenced on the 10x Genomics Chromium platform were aligned and quantified using the CellRanger 3.0.2 count function against the GRCh38.p12 genome assembly. To account for reads aligning to both exonic and intronic regions, each gene transcript in this reference assembly was relabeled as an exon as CellRanger counts exonic reads only. We perform this additional step since snRNA-seqterm·mth captures both mature mRNA and pre-mRNA, the latter of which includes intronic regions.
After aggregating each single-nucleus sample with the CellRanger aggr function, the full dataset was processed using Seurat v3.0.0 36 . The data were initially filtered for genes expressed in at least 3 cells and filtered for cells with reads quantified for between 200 and 2500 genes. We further filtered for cells that had percentage of counts coming from mitochondrial genes less than or equal to 5%. The data were normalized, scaled, and corrected for mitochondrial read percentages with sctransform v0.2.0 37 using default options.
To identify clusters, Seurat employs a shared nearest neighbor approach. We identified clusters using the top 10 principal components of the processed expression data with resolution set at 0.2. The resolution parameter controls the number of clusters that will be identified, and suggested values vary depending on the size and quality of the dataset. We chose a value that produced 6 clusters in the adipose dataset and 13 clusters in the DLPFC dataset and visualized the clustering results with UMAP 38 .
Marker genesterm·key were identified by determining the average log-fold change of expression of each cluster compared to the rest of the cells. We identified marker genesterm·key as those with an average log-fold change above 0.25. The significance of the differential expression of these genes was determined using a Wilcoxon rank sum test. Only genes that were detected in at least 25% of cells were considered. Clusters with many mitochondrial genes as markers (nine genes detected in both datasets) were removed from both datasets. In addition, a cluster with only three marker genesterm·key was removed from the DLPFC datasets. Finally, we remove mitochondrial genes from the list of marker genesterm·key for decompositionterm·key as we assume reads aligning to the mitochondrial genome originate from extra-nuclear RNA in the snRNA-seqterm·mth dataset (targeting nuclear RNA).
Clusters were labeled by considering cell types associated with the identified marker genesterm·key. Marker genesterm·key were downloaded from PanglaoDB 39 and filtered for entries validated in human cells. For each gene, we count the possible cell-type labels. Each cluster was labeled as the most frequent cell type across all of its marker genesterm·key, with each label associated with a gene weighted by the average log-fold change. If multiple clusters shared a cell-type label, we consider each cluster a subtype of this label.
Exon-aligned reads were processed in the same exact procedure but snRNA-seqterm·mth data was aligned to just exonic regions. Cluster names were manually changed for both datasets when aligned to exons to match the clusters from intronic reads as well. Specifically, for clusters identified in the exonic data not found in the full data, we relabeled as the label with the highest score found in the full data. These relabeled clusters were similar in proportion to the corresponding cluster in the full dataset.
Bisqueterm·mth reference-based decompositionterm·key model We assume that only a subset of genes are relevant for estimating cell-type composition. For the adipose and DLPFC datasets, we selected the marker genesterm·key identified by Seurat as described previously. Moreover, we filter out genes with zero variance in the single-cell data, unexpressed genes in the bulk expression, and mitochondrial genes. We convert the remaining gene counts to counts-per-million to account for variable sequencing depth. For m genes and k cell types, a reference profileterm·key Z ∊ R m × k is generated by averaging relative abundances within each cell type across the entire single-cell dataset.
Although there is a strong positive correlation between bulk and single-cell-based pseudo-bulk (summed single-cell counts) expression data, we observe that the relationship is not one-to-one and varies between genes. This behavior indicates that the distribution of observed bulk expression may significantly differ from the distribution of the single-cell profile weighted by cell proportions. We propose transforming the bulk data to maximize the global linear relationship across all genes for improved decompositionterm·key. Our goal is to recover a one-to-one relationship between the transformed bulk and expected convolutions of the reference profileterm·key based on single-cell based estimates of cell proportions. This transformed bulk expression better satisfies the assumptions of regression-based approaches under sum-to-one constraints.
Cell-type proportionsterm·key p ∊ R k × n ′ are determined by counting the cells with each label in the single-cell data for n ′ individuals. Given these proportions and the reference profileterm·key Z , we calculate the pseudo-bulk for the single-cell samples as Y = Zp , where Y ∊ R m × n ′ . For each gene j , our goal is to transform the observed bulk expression across all n bulk samples X j ∊ R n to match the mean and variance of Y j ∊ R n ′ ; hence, the transformationterm·mth of X j will be a linear transformationterm·mth.
If individuals with both single-cell and bulk expression are available, we fit a linear regressionterm·key model to learn this transformationterm·mth. Let X′ j ∊ R n ′ denote the expression values for these n ′ overlapping individuals. We fit the following model (with an intercept) and apply the model to the remaining bulk samples as our transformationterm·mth: 1 Y j = β j X j ′ + ϵ j If there are no single-cell samples that have bulk expression available, we assume that the observed mean of Y j is the true mean of our goal distribution for the transformed X j . We further assume that the sample variance observed in Y j is larger than the true variance of the goal distribution, as the number of single-cell samples is typically small. We use a shrinkage estimator of the sample variance of Y j that minimizes the mean squared error and results in a smaller variance than the unbiased estimator: 2 σ ^ j 2 = 1 n ′ + 1 ∑ i = 1 n ′ ( Y i , j − Ȳ j ) 2 We transform the remaining bulk as follows: 3 X j , transformed = X j − X ¯ j σ X j σ ^ j + Ȳ j where a bar indicates the mean value of the observed data and σ X j is the unbiased sample variance of X j .
To estimate cell-type proportionsterm·key, we apply non-negative least-squares regressionterm·mth with an additional sum-to-one constraint to the transformed bulk data. For individual i , we minimize the following with respect to the cell-proportion estimate p i : 4 ∣ ∣ Z p i − X i , transformed ∣ ∣ 2 s . t . p i ≥ 0 , ∑ p i = 1 Simulating bulk expression based on single-nucleus counts We simulate the base bulk expression as the sum of all counts across cells/nuclei sequenced from an individual. To introduce gene-specific variation between the bulk and single-cell data, we sample a coefficient β j and an intercept α j from a half-normal (HN) distributions: 5 β j ~ HN ( σ ) + 1 6 α J ~ HN ( σ ) At σ = 0, the base simulated bulk expression remains unchanged. We used a HN distribution to ensure coefficients and intercepts are positive. Although our method can handle negative coefficients, this simulation model assumes expression levels have a positive correlation across technologies. We performed 10 replicates of this data-generating process at each σ in 0, 5, 10, 20. Decomposition performanceterm·mth on these data were measured in terms of global R and RMSD and plotted with 95% confidence intervals based on bootstrapping.
Measuring significance of cell proportion-trait association Reported associations were measured in terms of Spearman correlation. To determine the statistical significance of these associations while accounting for possible confounding factors, we applied two approaches. For the adipose dataset, which consisted entirely of twin pairs, we applied a linear-mixed-effects model (R nlme package) with random effects accounting for family. For the DLPFC dataset, we assumed individuals were unrelated and fit a simple linear model (R base package). In each model, we include cell-type proportionterm·key, age, age-squared, and sex as covariates. We introduced an additional covariate for diabetes status when regressing Matsuda index due to a known significant association between these two variables. We test whether the cell proportion-effect estimates deviate significantly from 0 using a t -test. Each R method implements the described model fitting and significance testing.
Bisqueterm·mth marker-based decompositionterm·key model In order to estimate cell-type proportionsterm·key across individuals without the use of a cell-type-specific gene expression panel as reference, we use a weighted PCA approach. BisqueMarker requires a set of marker genesterm·key for each cell type as well as the specificity of each marker determined by the fold-change from a differential expression analysis. Typical single-cell RNA-seqterm·mth workflows calculate marker genesterm·key and provide both p -values and fold-changes, as in Seurat 36 . For each cell type, we take statistically significant marker genesterm·key (FDR < 0.05) ranked by p -value. A weighted PCA is calculated on the expression matrix using a subset of the marker genesterm·key by first scaling the expression matrix and multiplying each gene column by its weight (the log-fold change) XW , where X is the sample by gene expression matrix and W is a diagonal matrix with entries equal to log-fold change of the corresponding gene. The bulk expression X should be corrected for global covariates so that the proportion estimates do not reflect this global variation. The first PC calculated from XW is used as the estimate of the cell-type proportionterm·key. This allows cell-type-specific genes to be prioritized over more broadly expressed genes. Alternatively, if weights are not available, PCA can be run on the matrix X and the first PC can be used.
In order to select marker genesterm·key, we iteratively run the above PCA procedure on a specified range of markers (from 25 to 200) and calculate the ratio of the first eigenvalue to the second. We then select the number of marker genesterm·key to use that maximizes this ratio. This procedure is similar to other methods which select the number of markers to use by maximizing the condition number of the reference matrix 13 .
Software used Single-nucleus RNA-seqterm·mth data were aligned using CellRanger 3.0.2 against the GRCh38.p12 genome assembly. Bulk RNA-seqterm·key data were aligned with STAR 2.5.1 and quantified using featureCounts 1.6.3, both against the GRCh38.p12 genome assembly. R 3.5.1 was used for further processing and decompositionterm·key experiments. The Seurat v3.0.0 R package was used to filter, cluster, and identify cell-type marker genesterm·key from the single-nucleus data. The sctransform 0.2.0 R package was used to normalize and scale the single-nucleus data. Bisqueterm·mth 1.0, xbioc 0.1.7, Biobase 2.4.2, MuSiCterm·key 0.1.1, bseqsc 1.0, CIBERSORTterm·mth v1.06, and CIBERTSORTx were all used for decompositionterm·key using the processed bulk and single-nucleus RNA-seqterm·mth data. The R nlme 3.1-127 package was used for linear-mixed-model association. All visualizations and were generated with Python 3.7.2 using Seaborn 0.9.0, Matplotlib 3.0.3, Pandas 0.24.2, and Numpy 1.16.2, sklearn 0.20.3, and scipy 1.2.1.
Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Supplementary information Supplementary Information Peer Review File Reporting Summary