Determining cell-type abundance and expression from bulk tissues with digital cytometry

Tissue dissection with scRNA-seq

CIBERSORTx was designed to enable large-scale tissue characterization using cell signatures derived from diverse sources, including single-cell reference profiles (Fig. 1). To achieve this goal, we developed analytical tools for deriving a signature matrix from single-cell or bulk sorted transcriptional data while minimizing batch effects as a source of confounding technical variation (Supplementary Figs. 1,2; Supplementary Note 1). We then investigated the utility of commonly used scRNA-seq technologies¹⁶ for enumerating cell proportions in RNA admixtures derived from bulk tissues.

We started by generating a single-cell RNA-seq library from peripheral blood mononuclear cells (PBMCs) obtained from a patient with non-small cell lung cancer (NSCLC) using 10x Genomics Chromium v2 (3’ assay). Unsupervised clustering and canonical marker gene assessment revealed six major leukocyte subsets (B cells, CD4 T cells, CD8 T cells, NKT cells, NK cells, and monocytes; Fig. 2a). To assess deconvolution performance, we built a signature matrix to distinguish these cell subsets and tested it on a validation cohort of bulk RNA-seq profiles of blood obtained from 12 healthy adults (Supplementary Tables 1 and 2).

Compared to ground truth cell proportions as determined by direct cytometry and fluorescence immunophenotyping, uncorrected deconvolution results showed clear estimation biases for some cell types in bulk admixtures (Fig. 2b, Supplementary Fig. 1d). We hypothesized that these biases could be driven by platform-specific variation between the signature matrix and bulk RNA-seq data, as might be introduced, for example, by the variable use of unique molecular identifiers during library preparation. Indeed, following application of a batch correction scheme that we developed (Supplementary Fig. 1a, Supplementary Note 1, Methods), deconvolution results substantially improved and compared favorably to ground truth cell proportions (Fig. 2b, Supplementary Fig. 1d). We observed similar gains in performance through batch correction when analyzing other datasets and signature matrices, including publicly available Chromium v2 PBMC data (3’ and 5’ kits; Supplementary Fig. 2a-c) and purified leukocyte subsets profiled using microarrays (Supplementary Fig. 2d). Given these systematic improvements, we therefore applied batch correction in all subsequent cross-platform analyses, unless stated otherwise (Supplementary Table 1).

We next extended our analysis to solid tumor biopsies where single cells were profiled by SMART-Seq2. Focusing on head and neck squamous cell carcinomas⁹ (HNSCC) and melanomas¹⁰ (n = 18 and 19 patients, respectively), we initially tested deconvolution performance on simulated tumors reconstructed from single cells. This allowed us to evaluate the utility of single-cell reference profiles in a manner that controlled for dissociation-related artifacts and heterogeneity in phenotypic definitions. In a dataset of 18 primary tumors and 5 lymph node metastases from patients with HNSCC⁹, we created a signature matrix from a training cohort consisting of 2 primary tumor specimens and 1 lymph node biopsy (Supplementary Table 2). This matrix distinguished malignant cells, CD4 and CD8 T cells, B cells, macrophages, dendritic cells, mast cells, endothelial cells, myocytes, and cancer-associated fibroblasts (Fig. 2c). When evaluated using reconstructed tumor samples, deconvolution results were highly concordant with ground truth cell proportions (Fig. 2d, Supplementary Fig. 3a). Strong performance was also maintained when considering deconvolution results across distinct tumor types and cell types, including within rare or difficult to isolate cell subpopulations, such as distinct CD8 T cell effector subsets infiltrating melanomas (Supplementary Figs. 3b-d, 4, Supplementary Table 2).

To explore the impact of key signature matrix-related parameters on single cell-guided deconvolution, we next examined CIBERSORTx performance as a function of the number of cells per phenotype, the number of donor samples, and the number of genes considered. Across a range of values for these factors, we observed a surprisingly modest effect on cell proportion estimates (Supplementary Fig. 5a-c). Moreover, regardless of their primary biological source, CIBERSORTx signature matrices exhibited strong generalizability across diverse expression profiling platforms, datasets, and tissues after batch correction was applied (Fig. 2e).

Leveraging the single cell-derived signature matrix from melanoma biopsies described above, we then applied CIBERSORTx to dissect melanoma RNA-seq profiles from The Cancer Genome Atlas (TCGA)¹⁷. We observed substantial differences in the fractional representation of B/T lymphocytes and macrophages when comparing predicted cell type proportions in bulk tumors and the original scRNA-seq results (Supplementary Fig. 3e). Since these cell subsets were unselected relative to one another in the scRNA-seq dataset¹⁰, such compositional distortions may have arisen either from technical artifacts owing to single-cell isolation and sequencing^18,19, or from the deconvolution approach itself. In support of the former, IHC estimates of TIL subsets in an independent melanoma cohort²⁰ were far more similar to TIL fractions estimated by CIBERSORTx than those determined by scRNA-seq (Supplementary Fig. 3e). Moreover, we observed the same distortion phenomenon in a dataset of human pancreatic islets profiled by scRNA-seq (SMART-Seq2), bulk RNA-seq, and IHC²¹ (Fig. 2f, Supplementary Fig. 3f–i). In a direct comparison of paired islet specimens, cell fractions determined by IHC in bulk tissues were significantly correlated with bulk islet deconvolution results, but not scRNA-seq (Fig. 2f, Supplementary Fig. 3i). These data further validate CIBERSORTx and highlight its value for mitigating dissociation-related distortions resulting from the physical isolation of intact single cells (for additional discussion, see Supplementary Note 2).

Cell type-specific gene expression without physical cell isolation

Cell-type-specific transcriptome profiles can provide valuable insights into cell identity and function. However, such profiles are generally derived from single cells or bulk sorted populations, which can be difficult to obtain for large cohorts and fixed clinical samples. Even when purified cell types are available, tissue dissociation and preservation conditions can cause non-biological alterations in gene expression that obscure downstream analyses^18,19. While mathematical separation of bulk tissue RNA profiles into cell type-specific transcriptomes can potentially overcome these problems^22–28, the accuracy of this technique on real tissue samples remains unclear. We therefore set out to evaluate whether a signature matrix, consisting of highly optimized marker genes, can be used to faithfully reconstruct cell type-specific transcriptome profiles from non-disaggregated tissue samples, including fresh/frozen and fixed tumors (Fig. 3a).

We began by profiling 302 fresh/frozen primary tumor biopsies from patients with untreated follicular lymphoma (FL) and tested a common approach²⁵ in which cell type proportions are used to infer a single representative GEP for each cell type from a group of mixture samples (Methods). Since B cells, CD8 T cells, and CD4 T cells comprise the vast majority of FL tumor cellularity and can be readily purified by FACS²⁹, we focused on these three subsets to assess the accuracy of the approach. To enumerate FL immune proportions, we applied LM22¹⁵, a microarray-derived signature matrix for distinguishing 22 human hematopoietic cell subsets in bulk tissues, including tumors^30,31. We started by examining B cells and CD8 T cells as examples of highly abundant and less abundant cell types in FL lymph nodes (>50% versus ~5–10%, respectively¹⁵). Although imputed and FACS-purified cell type transcriptomes were reasonably well correlated, considerable noise distorted expression estimates of many genes (Fig. 3b). We therefore developed an adaptive noise filter to eliminate unreliably estimated genes for each cell type (Fig. 3c). After implementing this novel strategy within CIBERSORTx, we observed consistently improved correlations between in silico purified transcriptomes and those from FACS-purified cells (Fig. 3d, Supplementary Fig. 6a).

We next investigated key factors that influence the accuracy of transcriptome purification (Fig. 3e, Supplementary Fig. 6). Using the set of 302 FL GEPs, we observed a predictable relationship between the number of tumors profiled and the accuracy of transcriptome imputation across B cells, CD4 T cells, and CD8 T cells (Fig. 3e, Supplementary Fig. 6b). The largest gains were achieved when analyzing at least 4–5 fold more mixture samples than cell types (Fig. 3e). Nevertheless, our filtration scheme uniformly improved performance over a previous approach²⁵ irrespective of cohort size, the number of cell types, and overall cell type abundance (Fig. 3e, Supplementary Fig. 6b-e,g). We also observed favorable performance for resolving previously identified markers of cell identity and additional cell type transcriptomes (Fig. 3f,g, Supplementary Fig. 6d-h). As expected, the fraction of recovered genes after adaptive filtration was proportional to both the number of evaluated samples and the proportion of each cell subset (Supplementary Fig. 6h). Moreover, the use of different single cell-derived signature matrices did not significantly impact results, provided that identical cell types were interrogated (Supplementary Fig. 5c).

Cell type-specific expression purification at high-resolution

Having shown that cell type-specific GEPs can be reliably estimated, we next turned to the problem of inferring cell type-specific differential expression from bulk specimens. Despite the utility of the above technique, it is limited to learning a single representative expression profile for each cell type given a group of mixture samples (i.e., group-mode GEPs, Fig. 1). While useful for comparing two groups of specimens, such profiles are not sample-specific and must be generated for each group of interest in order to study differentially expressed genes (DEGs) between them. Although approaches for sample-level deconvolution have been previously described, they only consider mixtures with two or three cellular components^32,33. We therefore developed a framework that generalizes to multiple (>3) components by modeling gene expression deconvolution as a unique non-negative matrix factorization problem with partial observations (Supplementary Fig. 7, Methods). Briefly, our approach attempts to separate a single matrix of mixture GEPs into a set of underlying cell type-specific expression matrices using imputed cell proportions (Fig. 4a,b). Once these expression profiles are obtained, they can be analyzed post hoc to gain insights into sample-level variation and patterns of gene expression for individual cell types of interest. To solve the matrix factorization problem, we implemented a novel divide and conquer algorithm that produces biologically realistic solutions (Supplementary Figs. 7 and 8, Methods).

To test the method’s capability for “high-resolution” cell purification (Fig. 1), we created a series of synthetic mixtures, each containing DEGs in one or more cell types. These DEGs were simulated to include overlapping block-like patterns, reminiscent of those seen in real tissues^34–36, and non-linear geometries, all of which would be difficult to ascertain by previous computational techniques. Remarkably, the method recovered expected DEG patterns in all tested cases, including an obscured target (“bullseye”) (Fig. 4b-d, Supplementary Fig. 9). Moreover, unlike group mode (Fig. 1), the resulting high-resolution profiles were amenable to standard methods for unsupervised analysis (e.g., Fig. 4c).

Using modeled tumor admixtures, we next evaluated the analytical performance of the method across several parameters, including cell type abundance and the magnitude of differential expression. First, simulated DEGs were ‘spiked’ into CD8 T cell transcriptomes to create two phenotypic classes. These CD8 GEPs were then randomly admixed in silico with three other immune subsets in modeled tumors, and a colon cancer cell line was included to simulate 50% unknown content (Fig. 4e, left). Following high-resolution purification, cell type-specific transcriptomes were grouped into defined DEG classes. Across a broad range of cell spike-in levels and expression fold changes, previously defined DEGs were recovered in CD8 T cells with high sensitivity and specificity (Fig. 4e, right). We observed similar performance in a second experiment, where we simulated melanoma tumors using pooled scRNA-seq data¹⁰ and assessed the above parameters in relation to cohort size (Supplementary Fig. 10).

High-resolution profiling of diverse tumor subpopulations

Diffuse large B cell lymphoma (DLBCL) can be classified into two major molecular subtypes based on differences in B cell differentiation states: germinal center-like (GCB) and activated B cell-like (ABC) DLBCL³⁴. Using GEPs of 150 DLBCL lymph node biopsies with previously annotated cell-of-origin subtypes³⁷, we asked whether high-resolution profiling of 10 major leukocyte subsets could correctly attribute known cell-of-origin differences to B cells. Although our approach was blinded to class labels, we identified DLBCL subtype-specific expression differences in malignant B cells that were (1) highly consistent with those of normal GC and activated B cells and (2) ~9-fold more significant, on average, than in bulk DLBCL tumors (Supplementary Fig. 11a-c). Notably, we obtained similar results when repeating this analysis using 10x Chromium-derived signature matrices derived from either peripheral blood or FL tumors, demonstrating the generalizability of our approach (Supplementary Fig. 11d,e).

We then compared our results with two alternative methods: (1) a common strategy for assigning bulk tissue expression patterns to individual cell types based on correlations with cell abundance³⁸, and (2) a previously described technique for imputing cell type-specific DEGs when phenotypic classes and cell type frequencies are known²⁵. In both cases, CIBERSORTx exhibited superior performance, both in relation to cell type specificity and the number of detectable DEGs at a given significance threshold (Supplementary Figs. 12, 13a-b).

FL is the most common indolent Non-Hodgkin lymphoma, and CREBBP mutations in FL tumors are associated with reduced antigen presentation in B cells^29,39. Since we originally discovered this association using FACS-purified FL B cells²⁹, we asked whether high-resolution purification could recapitulate this result starting from bulk tumor GEPs and paired tumor genotypes (Fig. 5a). Indeed, after stratifying tumors by CREBBP mutation status (Supplementary Table 3), previously described signatures were detectable in digitally sorted B cell GEPs, including loss of MHC II expression in CREBBP-mutant tumors (Fig. 5b-c). Notably, this result was reproducible across microarray and 10x Chromium-derived signature matrices covering leukocyte subsets derived from distinct biological sources (Fig. 5b-c, Supplementary Fig. 14a,b). Moreover, as observed for DLBCL (Supplementary Fig. 11), the majority of CREBBP mutation-associated genes did not correlate with B cell abundance, hindering their discovery in bulk tissues without deconvolution (Supplementary Fig. 14c).

To extend our assessment to solid tumors, we obtained surgically-resected primary NSCLC tumor biopsies (n = 26 patients), and for each tumor, generated RNA-seq libraries of four major subpopulations purified by FACS: epithelial/cancer, hematopoietic, endothelial, and fibroblast subsets (Fig. 5d, Supplementary Table 1). After deriving a signature matrix from a subset of four patients, we applied high-resolution profiling to digitally dissect these four populations in bulk tumors from the remaining 22 patients, all of which had RNA-seq profiles available. In a direct comparison against FACS-purified cell populations on matching patient samples, in silico profiles showed strong evidence of expression purification (Supplementary Fig. 13c-f). In addition, CIBERSORTx outperformed other methods^6,33 for purifying GEPs of epithelial cells from bulk tumors while also enabling the digital purification of more cell types (Supplementary Fig. 13c-e).

We applied the same signature matrix to resolve epithelial/cancer, hematopoietic, endothelial, and fibroblast GEPs from bulk RNA-seq profiles of 518 lung adenocarcinoma (LUAD) tumors, 504 lung squamous cell carcinoma (LUSC) tumors, and 110 adjacent normal tissues from TCGA^40,41 (Fig. 5d). Using t-SNE to visualize the resulting GEPs, we identified striking patterns of histology-specific gene expression for most cell types, including distinct phenotypic shifts in cancer-associated fibroblasts (CAFs) (Fig 5e). Histological differences were far less pronounced in tumor-associated endothelial cells, and adjacent normal tissues clustered together regardless of histology. We compared these results with bulk RNA-seq profiles of FACS-purified NSCLC cell subpopulations from 21 patients with LUAD or LUSC. We observed similar clustering tendencies at the whole transcriptome level and strong concordance in relation to patterns of cell type-specific differential expression between histological subtypes (Fig. 5e,f, Supplementary Table 4, Supplementary Note 1). Our results were further corroborated by histology-specific DEGs and tumor-specific DEGs identified from a recently published scRNA-seq atlas of NSCLC tumors and adjacent normal tissues⁴² (Supplementary Fig. 15, Supplementary Table 4).

Applications of CIBERSORTx to melanoma

We implemented the set of CIBERSORTx techniques into a comprehensive toolkit (http://cibersortx.stanford.edu). We then explored three potential applications of CIBERSORTx for characterizing cellular heterogeneity in resected tumor biopsies from patients with melanoma. The following techniques were applied in turn: high-resolution expression purification (Fig. 6a-b), group-mode expression purification (Fig. 6c), and enumeration of cell composition across diverse platforms using single-cell reference profiles (Fig. 6d-f).

Oncogenic BRAF mutations occur in over half of melanomas and can be inhibited by approved targeted therapies^43,44 whereas NRAS mutations occur in approximately half of non-BRAF mutant melanoma tumors but lack such therapies⁴⁵. Understanding how key mutations influence cellular states could potentially lead to novel treatment strategies. Using single-cell reference profiles from melanomas to build a signature matrix, we applied high-resolution expression purification to dissect 8 major cell types from the transcriptomes of 342 bulk melanoma tumors profiled by TCGA⁴⁴ (Fig. 6a). Within digitally purified cell subsets, we discovered many significant DEGs within malignant cells and CAFs that distinguish melanomas according to BRAF or NRAS mutation status (Fig. 6b, Supplementary Table 4). We verified these findings using scRNA-seq data from primary melanomas where mutation data were available, allowing us to confirm GEPs associated with BRAF and NRAS genotypes within individual malignant cells and/or CAFs (Fig. 6b, Supplementary Fig. 16).

Tumor-infiltrating CD8 T cells are driven to a state of “exhaustion” by chronic antigen stimulation or by overexposure to inflammatory signals⁴⁶. Given the importance of these cells for current and emerging cancer immunotherapies^47,48, we next used CIBERSORTx to examine expression changes that characterize the exhaustion phenotype. Using LM22, which is derived from healthy peripheral blood leukocytes, we enumerated immune composition in fresh/frozen melanomas profiled by TCGA⁴⁴. We then performed group-mode expression purification to impute a representative CD8 TIL GEP. By rank-ordering the estimated CD8 TIL GEP against a baseline reference profile of normal peripheral blood CD8 T cells, we confirmed the expression of key exhaustion markers in the inferred CD8 TIL GEP (Fig. 6c). In addition, CD8 TIL-specific genes were consistent with those observed for CD8 TILs isolated from melanomas by single-cell RNA-seq¹⁰ and by FACS⁴⁹ (Fig. 6c). Similar results were obtained when repeating the analysis on FFPE tumors⁵⁰ (Fig. 6c, Supplementary Fig. 17).

The most effective regimens for metastatic melanoma currently employ checkpoint blockade targeting PD-1 and/or CTLA4 expression on exhausted T cells⁴⁷. Although a subset of patients achieve durable anti-tumor T cell responses, clinical outcomes remain heterogeneous and effective predictive biomarkers are lacking⁵¹. CD8 TILs expressing high levels of PDCD1 (encoding PD-1) or CTLA4 are key targets of these therapies^52,53, suggesting that CD8 TILs expressing both markers might correlate with response, as was recently shown in melanoma patients receiving PD1 blockade⁵⁴. To test this hypothesis, we used single-cell reference profiles of melanoma tumors to build a signature matrix containing PDCD1⁺/CTLA4⁺ CD8 T cells along with eight other major tumor cell types (Fig. 6d, Supplementary Figs. 3b, 18a,b). We then applied the signature matrix to interrogate three publicly available melanoma expression datasets (Fig. 6d). These included bulk expression data of FFPE and fresh/frozen melanoma tumors that were profiled by RNA-seq or NanoString^50,55. In support of our hypothesis, imputed levels of PDCD1⁺/CTLA4⁺ CD8 T cells were significantly associated with response in all three studies (P < 0.05; Fig. 6e). Moreover, the detection of these cells stratified overall survival in this meta-analysis, separated survival curves in individual datasets, and more robustly associated with survival and response than key marker genes and other cell types (Fig. 6f, Supplementary Fig. 18c–f).

Human subjects.

All patient samples in this study were collected with informed consent for research use and were approved by the Stanford Institutional Review Board in accordance with the Declaration of Helsinki. For a patient with metastatic NSCLC treated with an immune checkpoint inhibitor (Pembrolizumab, Merck), peripheral blood was obtained on the first day of treatment prior to infusion (Fig. 2a; NCT00349830 and NCT02955758). Fresh tumor biopsies from patients with early stage NSCLC were obtained during routine primary surgical resection (Figs. 3g and 5, Supplementary Fig. 13c-f). Fresh or frozen surgical biopsies of follicular lymphoma tumors were obtained from previously untreated FL patients enrolled in a phase III clinical trial (NCT00017290⁶³), as well as from patients seen as part of the Stanford University Lymphoma Program Project (NCT00398177; Figs. 3b-f, 5a-c, Supplementary Figs. 6a-c, 14). Whole blood samples from 12 healthy adult donors were obtained from the Stanford Blood Center (Fig. 2b,e and Supplementary Figs. 1d,k,l, 2a,d).

External datasets.

Full details of each dataset^{9,10,17,20,21,37,42,50,55,64–69}, including data type, sample type, source, and normalization approach, are available in Supplementary Table 1. Briefly, next generation sequencing datasets were downloaded and analyzed using the authors’ normalization settings unless otherwise specified; these consisted of transcripts per million (TPM), reads per kilobase of transcript per million (RPKM), or fragments per kilobase of transcript per million (FPKM) space. For analyses in log₂ space, we added 1 to expression values prior to log₂ adjustment. Affymetrix microarray datasets were summarized and normalized as described in ‘Gene expression profiling – Microarrays’ (Supplementary Note 1), using RMA in cases where bulk tissues and ground truth cell subsets were profiled on the same Affymetrix platform, and otherwise using MAS5 normalization. NanoString nCounter data were downloaded from the supplement of Chen et al.²⁰ and analyzed with batch correction in non-log linear space, but without any additional preprocessing.

Two publicly available PBMC datasets from healthy donors profiled by Chromium v2 (5’ and 3’ kits) were downloaded (Supplementary Table 1) and preprocessed as described in ‘Gene expression profiling – Single-cell RNA-seq’ (Supplementary Note 1), with the following minor modifications. During quality control, we excluded cells with >5000 expressed genes for 5’ PBMCs, >4000 expressed genes for 3’ PBMCs, and <200 expressed genes for both datasets. Seurat “FindClusters” was applied on the first 20 principal components, with the resolution parameter set to 0.6. Cell labels were assigned as described above. In addition, myeloid cells were defined by high CD68 expression, megakaryocytes by high PPBP expression, and dendritic cells by high FCER1A expression.

For the 3’ FL signature matrix in Supplementary Figs. 11d, and 14a-b, publicly available 10x Chromium v2 scRNA-seq data (3’ kit)⁷⁰ were downloaded (Supplementary Table 1) and preprocessed as described for the 10x PBMC signature matrices above, but with the following differences. Seurat “FindClusters” was applied on the first 10 principal components, with the resolution parameter set to 0.6. Cell labels were assigned based on the following canonical marker genes (MS4A1 = B cells; CD3E, CD8A and CD8B = CD8 T cells; CD3E and CD4 = CD4 T cells).

Single-cell signature matrix construction.

Expression data from input datasets (Supplementary Table 1) were summarized as described above. Given the variability in cell type representation and the inherent noise in scRNA-seq data, several techniques have been developed to address stochastic dropout, impute cell-level scaling factors, and model technical and biological noise components in single-cell differential expression analysis^71–73. Although useful for defining single cell phenotypes, we did not observe any significant gains in deconvolution performance when applying such techniques after cell labels were assigned (data not shown). Since the discovery of cell subpopulations in scRNA-seq data was not the focus of this work, we limited our pre-processing steps to the following approach, which performed well in our experiments. Single-cell expression values were first normalized to TPM and divided by 10 to better approximate the number of transcripts per cell¹⁰. For each cell phenotype, genes with low average expression (<0.75 transcript per cell) in log₂ space were then set to 0 as a quality control filter. Although an expression threshold of <0.75 was used in this work, downstream results remained comparable when using modestly different thresholds (data not shown). Owing to sparser gene coverage relative to SMART-Seq2, this filter was not applied to data generated by 10x Chromium. For each cell type represented by at least 3 single cells, we selected 50% of all available single-cell GEPs using random sampling without replacement (fractional sample sizes were rounded up such that 2 cells were sampled if only 3 were available). We then aggregated the profiles by summation in non-log linear space and normalized each population-level GEP into TPM. This process was repeated in order to generate 5 aggregated transcriptome replicates per cell type. Unless stated otherwise, scRNA-seq and bulk RNA-seq signature matrices were generated as described previously¹⁵ using the following parameters: minimum number of genes per cell type = 300, maximum number of genes per cell type = 500, q-value <0.01 for differential gene expression, and no quantile normalization.

All single cell-derived signature matrices are available in Supplementary Table 2.

Overview of CIBERSORTx analytical framework.

Cross-platform deconvolution.

While several prior studies have applied CIBERSORT to RNA-seq, including to reference phenotypes derived from single-cell transcriptome profiling^11,14, our original description of CIBERSORT did not explicitly handle technical variation between the signature matrix and bulk mixture profiles. To the best of our knowledge, there is no previously described approach that can be applied to eliminate technical variation between mixture sample expression profiles (denoted M) and a signature matrix comprised of cell type reference profiles (denoted B) while preserving biological signal. This is because, unlike technical batches of the same sample type, both biological differences and technical batches are inherently conflated in the setting of deconvolution. Since previous techniques for addressing technical dropout and cross-platform normalization in bulk GEP and scRNA-seq data are not directly applicable to this problem^57,58,78,79, we developed a novel approach for cross-platform deconvolution (Supplementary Fig. 1). By exploiting the fact that M can be modeled as a linear combination of B, and by estimating M from B (denoted M*), batch correction can be directly applied to M and M*. Our approach comprises two distinct strategies (B-mode and S-mode) to apply gene expression deconvolution across distinct platforms and tissue storage types (e.g., fresh/frozen versus FFPE). A decision tree to guide users in selecting the most appropriate strategy is provided in Supplementary Figure 1a. A comprehensive description of our approach is provided in Supplementary Note 1.

Software implementation and website.

Similar to its predecessor, CIBERSORTx was developed within a web framework with its back-end based on R and PHP and hosted at http://cibersortx.stanford.edu. This web framework minimizes inherent dependencies on specific hardware, software packages and libraries, and file-system attributes. Users are presented with a detailed guide employing several step-by-step Tutorials, and allowing the recreation of key figures in this work, including for each step depicted in Figure 1. Through this interface, CIBERSORTx allows users to process gene expression data representing a bulk admixture of different cell types, along with (1) a signature gene file that enumerates the genes defining the expression profile for each cell type of interest. For the latter, users can either use existing/curated signature matrices for reference cell types, or can create custom signature gene files by providing the reference gene expression profiles of pure cell populations. Specifically, to create a custom signature gene matrix, users can provide single-cell RNA sequencing data or data from bulk sorted samples, along with the phenotypic identities of single cell types or cell populations of interest.

Given these input files, CIBERSORTx allows (2) imputation of the fractional representations of each cell type present in the mixture, similar to its predecessor. However, unlike CIBERSORT, CIBERSORTx now supports deconvolution from bulk RNA-seq data by implementing the critical batch correction methods described above. CIBERSORTx also allows imputation of GEPs for individual in silico purified cell-types in two distinct modes as described above (i.e., (3) group-mode and (4) high-resolution). The resulting imputed cell fractions and imputed cell type-specific GEPs are then rendered as heat maps, tables, stacked bar plots for visualization and downloading. In addition, customizable t-SNE plots are automatically generated for high-resolution purification results.

The interactive CIBERSORTx user interface is powered by the jQuery JavaScript library and various open source libraries (including phpMailer, idiorm, blueimp jQuery-File-Upload, DataTables, phpExcel and mPDF), with the graphical user interface of the website powered by Twitter Bootstrap 2.3.2 and R Shiny. The site runs on an Apache server on a virtual machine and stores user and job data in a MySQL database. However, users have complete control over their data and can delete them at will. Each user’s environment includes example datasets used for benchmarking, tutorials for the use of CIBERSORTx and preparation of input data, and other example files.

Statistical analysis.

Linear concordance between known and predicted cell-type features (i.e., proportions or GEPs) was determined by Pearson correlation (r), Spearman correlation (rho), or linear regression (R²), as indicated, and a two-sided t test was used to assess whether the result was significantly nonzero. Lin’s concordance correlation coefficient (CCC) was determined by comparing predicted and expected log₂-adjusted expression profiles using the CCC function from the R package, DescTools. When data were normally distributed, group comparisons were determined using a two-sided t test with unequal variance or a paired t test, as appropriate; otherwise, a two-sided Wilcoxon test was applied. Multiple hypothesis testing was performed using the Benjamini and Hochberg method unless stated otherwise. Results with P < 0.05 were considered significant. Statistical analyses were performed with R and Prism v7 (GraphPad Software, Inc.). The investigators were not blinded to allocation during experiments and outcome assessment. No sample-size estimates were performed to ensure adequate power to detect a pre-specified effect size.

Code availability.

CIBERSORTx v1.0 was used to generate the results in this work and is freely available for academic research use at http://cibersortx.stanford.edu.

Data availability.

All expression datasets analyzed in this work, including accession codes, file names, and web links (if available), are listed in Supplementary Table 1. Expression data generated in this study are available at http://cibersortx.stanford.edu and through the Gene Expression Omnibus with accession code GSE127472.