How to Deconvolve Your Own Bulk RNA-Seq Into Cell-Type Fractions

By the end of this guide you will have a table with one row per sample and one column per cell type, holding estimated fractions for your own bulk RNA-seq, plus a set of diagnostics that tell you which of those numbers you should believe. You need a bulk RNA-seq run (whole blood, PBMC, or a tissue biopsy) with FASTQ files or a gene-level count matrix, roughly 30 million or more paired-end reads per sample, R 4.3+ with Bioconductor, and about 32 GB of RAM if you plan to run a single-cell reference through BayesPrism. If you have several timepoints, you will get considerably more out of this than from a single sample, because the useful signal in deconvolution is usually the change within one person rather than the absolute value.

Deconvolution is the inverse of a mixing operation. A bulk RNA-seq sample is a weighted sum: each cell type contributes its own expression profile in proportion to how many of those cells are present and how much RNA each one carries. Deconvolution takes the observed mixture and a reference describing the pure components, then solves for the weights ¹. The same idea appears throughout measurement science, from separating co-eluting compounds in GC-MS chromatograms ² to estimating cell-type-specific protein abundance in mass spectrometry proteomics ³. For transcriptomics the practical question is whether the reference you have matches the tissue you measured, and most failures trace back to that mismatch ⁴.

1. Quantify to gene-level TPM, and stay in linear space

Deconvolution methods assume the mixture is a linear combination of component profiles. Anything that breaks linearity before the solver sees the data will bias the result, often badly. That includes log transforms, variance-stabilizing transforms and z-scoring across samples. Quantify with Salmon or kallisto, then summarize to gene level with tximport.

salmon index -t gencode.v44.transcripts.fa.gz -i salmon_idx_v44 \
  --gencode -k 31 -p 16

salmon quant -i salmon_idx_v44 -l A \
  -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  --numBootstraps 30 -p 16 -o quant/sample01

library(tximport); library(AnnotationDbi); library(org.Hs.eg.db)

files <- file.path("quant", samples, "quant.sf"); names(files) <- samples
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "no", ignoreTxVersion = TRUE)

tpm <- txi$abundance                      # gene-level TPM, linear
sym <- mapIds(org.Hs.eg.db, sub("\\..*", "", rownames(tpm)),
              "SYMBOL", "ENSEMBL")
tpm <- tpm[!is.na(sym), ]; rownames(tpm) <- sym[!is.na(sym)]
# collapse duplicate symbols by max expression, not by sum
tpm <- rowsum(tpm, rownames(tpm)) / as.vector(table(rownames(tpm)))

Two details matter here. Signature matrices such as LM22 are keyed on HGNC gene symbols, so the mapping step is not optional, and duplicate symbols must be collapsed deterministically or the solver will silently use whichever row it hits first. Second, keep TPM rather than counts: TPM corrects for transcript length, and length bias between cell types is a real source of error in regression-based methods ⁵.

If your input is whole blood collected in PAXgene or Tempus tubes without globin depletion, hemoglobin transcripts will dominate. Check before you do anything else.

globin <- c("HBB","HBA1","HBA2","HBD","HBG1","HBG2")
colSums(tpm[intersect(globin, rownames(tpm)), , drop = FALSE]) / 1e4  # % of TPM

If that figure is above about 20%, remove the globin genes and renormalize the remaining TPM to sum to 1e6 per sample. Otherwise every other gene is compressed toward zero and the effective sequencing depth for the signature genes is a fraction of what you paid for.

2. Decide what you are estimating before you pick a method

There are three quantities people call “cell fractions,” and conflating them produces arguments that never resolve. The first is the fraction of RNA in the sample contributed by a cell type. The second is the fraction of cells. The third is an absolute score on some arbitrary scale. Methods differ in which one they return, because cell types differ several-fold in total RNA content: a neutrophil carries far less mRNA than a monocyte or a plasma cell, so RNA fractions systematically understate granulocytes relative to a blood count ¹.

Methods that correct for RNA content using an explicit per-cell-type scaling factor (quanTIseq, EPIC, MuSiC) return something closer to cell fractions. CIBERSORT in relative mode returns RNA fractions normalized to sum to one, which means the numbers are compositional and moving one cell type moves all the others. Choose deliberately, and record the choice next to the output table, because you will not remember it in six months.

3. Run marker-and-signature methods first through immunedeconv

For blood and PBMC, start with the established signature-matrix methods. The immunedeconv R package wraps several of them behind one interface, which makes it cheap to run four methods and compare.

library(immunedeconv)

res_quantiseq <- deconvolute(tpm, "quantiseq", tumor = FALSE)
res_epic      <- deconvolute(tpm, "epic",      tumor = FALSE)
res_abis      <- deconvolute(tpm, "abis")
res_xcell     <- deconvolute(tpm, "xcell")

quanTIseq resolves ten immune populations and returns fractions that include an “uncharacterized” remainder, which is the single most useful diagnostic in the whole output. If that remainder is 60% in a PBMC sample, the reference does not describe your sample and you should stop and investigate rather than reporting the other numbers. EPIC works similarly and also exposes an “otherCells” term.

xCell is different in kind. It produces enrichment scores for 64 cell types. Those scores are comparable across samples for the same cell type. They are not fractions and do not sum to anything meaningful. No comparison across cell types within a sample is valid. That distinction is easy to lose once the numbers are in a data frame.

CIBERSORTx sits outside immunedeconv behind a licensed Docker image. Its LM22 signature matrix covers 22 leukocyte subsets and is still the most widely used reference for blood. If you run it, use 100 or more permutations so you get a p-value per sample, and treat samples with p above 0.05 as failures of fit rather than as low-confidence estimates.

docker run -v $PWD:/src/data -v $PWD/out:/src/outdir cibersortx/fractions \
  --username you@example.com --token $CSX_TOKEN \
  --mixture mixture_tpm.txt --sigmatrix LM22.txt \
  --perm 200 --QN FALSE --absolute TRUE

Set --QN FALSE for RNA-seq input. Quantile normalization was introduced for microarray mixtures and degrades RNA-seq results.

4. Add a single-cell reference for anything that is not standard blood

Signature matrices are fixed vocabularies. If your sample contains a cell state that LM22 does not name, the solver will distribute that expression across whatever columns fit best. A single-cell reference lets you define the vocabulary yourself, which matters for tissue, for bone marrow, and for blood when you care about subsets finer than LM22 provides.

The reference should come from the same tissue and ideally the same technology. For PBMC, the Hao et al. CITE-seq atlas and the Tabula Sapiens blood compartment are both reasonable starting points. MuSiC is the lightweight option: it weights genes by their cross-subject consistency in the reference, which downweights genes whose expression varies more between donors than between cell types.

library(MuSiC); library(SingleCellExperiment)

sce <- SingleCellExperiment(assays = list(counts = sc_counts),
                            colData = sc_meta)   # needs cellType, sampleID

est <- music_prop(bulk.mtx = as.matrix(round(tpm)),
                  sc.sce   = sce,
                  clusters = "cellType",
                  samples  = "sampleID",
                  select.ct = c("CD4 T","CD8 T","NK","B","CD14 Mono",
                                "CD16 Mono","DC","Platelet"))
fractions <- est$Est.prop.weighted

BayesPrism is slower and usually better behaved on tissue. It treats the problem as Bayesian inference over a joint distribution of cell-type fractions and cell-type-specific expression, so it can update the reference profiles toward the observed mixture instead of holding them fixed. That property is what makes it more tolerant of a reference drawn from a different individual or batch.

library(BayesPrism)

prism <- new.prism(reference = sc_ref_counts,   # cells x genes, raw counts
                   mixture   = t(bulk_counts),  # samples x genes
                   input.type = "count.matrix",
                   cell.type.labels = ct_labels,
                   cell.state.labels = cs_labels,
                   key = NULL, outlier.cut = 0.01, outlier.fraction = 0.1)

bp <- run.prism(prism = prism, n.cores = 16)
theta <- get.fraction(bp, which.theta = "final", state.or.type = "type")

Before running either, filter the reference. Drop mitochondrial and ribosomal protein genes, and drop genes detected in fewer than about 1% of cells. Also drop cell types represented by fewer than 50 cells. Rare types with thin support produce noisy profiles that the solver will happily assign real weight to.

There is a newer class of methods that avoids supplying a reference at all. UniCell Deconvolve is trained across a large corpus of single-cell transcriptomes and returns cell-type probabilities without a user-specified signature matrix, which is useful when you genuinely do not know what is in the sample ⁶. Autoencoder approaches extend the same idea across modalities, jointly dissociating bulk transcriptome and other omics layers into cell-type components ⁷. We treat these as a second opinion rather than a primary estimate, because when they disagree with a reference-based run it is rarely obvious which is right.

5. Check the answer against something you measured independently

Run at least four methods and compare. Agreement across methods on the major populations is weak evidence of correctness. Disagreement is strong evidence that something is wrong, and it is common: benchmarking studies consistently find that method choice and reference choice together account for more variation in the output than the underlying biological differences between samples ^{5 4}.

The decisive check for blood is a complete blood count with differential drawn from the same blood draw. Correlate your estimated fractions against the measured neutrophil, lymphocyte, monocyte, eosinophil, and basophil percentages.

cbc <- data.frame(sample = samples,
                  neut = c(...), lymph = c(...), mono = c(...))
cor(fractions[,"Neutrophil"], cbc$neut/100, method = "spearman")

Expect a good correlation for monocytes and total lymphocytes, and a weaker one for neutrophils (RNA content again). Expect essentially nothing for basophils, which sit near the detection floor of every method. Interpreting the clinical meaning of a CBC differential, or any discrepancy between it and your estimates, is a matter for a physician. That interpretation does not belong to a deconvolution pipeline.

A second check uses canonical marker genes directly. Readers often ask which genes mark which immune cells, and the short answer for blood is the following.

• CD3D, CD3E, and IL7R for T cells, with CD8A and CD8B separating the cytotoxic compartment
• MS4A1 and CD79A for B cells
• GNLY, NKG7, and KLRD1 for NK cells
• LYZ, CD14, and S100A8/S100A9 for classical monocytes, with FCGR3A marking the non-classical subset
• FCER1A and CD1C for conventional dendritic cells
• PPBP and PF4 for platelets

If your estimated CD14 monocyte fraction rises across timepoints, LYZ and S100A8 TPM should rise with it. When a fraction moves and its markers do not, the solver is redistributing weight. The change does not come from detecting cells. The five broad classes most people mean by “main immune cell types” are lymphocytes, monocytes and macrophages, neutrophils, eosinophils, and basophils, with dendritic cells often counted separately.

6. Analyze longitudinal fractions as compositions

Cell fractions sum to one, which means they live on a simplex and ordinary statistics mislead. An apparent 15% drop in T cells may reflect a genuine expansion of monocytes and nothing else. Transform before you model.

library(compositions)
clr_frac <- as.matrix(clr(acomp(fractions + 1e-4)))

The centered log-ratio transform maps the composition to unconstrained space so that a change in one component no longer mechanically forces the others. For a single person across timepoints, fit each cell type’s CLR value against time and look at the residual spread. You also need a noise floor: rerun the same library through the pipeline with a different random seed, or better, sequence a technical replicate of one sample and compute the absolute difference in each fraction. In our experience with whole blood, run-to-run variation of one to two percentage points on the major populations is normal, so a change smaller than that is not a change.

Finally, keep the input matrix, the reference, the method version, and the exact parameters alongside the output. Deconvolution results are not portable across pipeline versions, and a fraction with no provenance is not interpretable six months later ⁸.

Common problems

The estimated fractions are nearly identical across all your samples. This usually means the signature genes are dominated by a few high-expression housekeeping-adjacent transcripts, or the mixture was z-scored before input. Verify that the matrix you passed is linear TPM with row-wise dynamic range intact, and check whether one gene accounts for an outsized share of the signature-gene sum.

quanTIseq or EPIC reports a large uncharacterized fraction. Your tissue contains cell types the reference does not cover. For a biopsy this is expected and you need a tissue-matched single-cell reference. For blood it usually means globin contamination, heavy rRNA carryover, or a degraded library.

Neutrophils come out far below the CBC value. Two causes. If the sample was PBMC rather than whole blood, neutrophils were removed by density gradient and should be near zero. If it was whole blood, this is the RNA-content effect, and you should switch to a method that applies per-cell-type RNA scaling.

Two methods give opposite directions for the same cell type over time. Trust neither until the marker genes adjudicate. Plot raw TPM for three or four canonical markers of that cell type across timepoints and see which method’s trajectory they follow.

BayesPrism runs for hours or exhausts memory. Downsample the reference to a few hundred cells per type, filter to about 5,000 to 8,000 informative genes, and drop cell states with fewer than 20 cells. Runtime scales with cells times genes times samples, and the fraction estimates are stable well below the full reference size.

Results shift after you add new samples to the batch. Any method that normalizes across the mixture matrix couples your samples together. For longitudinal tracking, deconvolute each sample independently, or fix the batch and reprocess everything whenever it changes. The underlying difficulty is that the reference, the normalization, and the mixture all interact, and no current method is insensitive to all three ^{9 4}.

Footnotes

1. Yebin Im, Yongsoo Kim. A Comprehensive Overview of RNA Deconvolution Methods and Their Application. Molecules and Cells, 2023. https://doi.org/10.14348/molcells.2023.2178 ↩ ↩²

2. Xiuxia Du, Steven H Zeisel. SPECTRAL DECONVOLUTION FOR GAS CHROMATOGRAPHY MASS SPECTROMETRY-BASED METABOLOMICS: CURRENT STATUS AND FUTURE PERSPECTIVES. Computational and Structural Biotechnology Journal, 2013. https://doi.org/10.5936/csbj.201301013 ↩

3. Yingnan Song, Qi Zhou, Chen Huang. Evaluating and Optimizing Mass Spectrometry Proteomics Data to Deconvolve Cell-Type-Specific Protein Expression in Tumors. Journal of Proteome Research, 2025. https://doi.org/10.1021/acs.jproteome.5c00980 ↩

4. Lana X. Garmire, Yijun Li, Qianhui Huang, et al. Challenges and perspectives in computational deconvolution of genomics data. Nature Methods, 2024. https://doi.org/10.1038/s41592-023-02166-6 ↩ ↩² ↩³

5. Francisco Avila Cobos, Jo Vandesompele, Pieter Mestdagh, et al. Computational deconvolution of transcriptomics data from mixed cell populations. Bioinformatics, 2018. https://doi.org/10.1093/bioinformatics/bty019 ↩ ↩²

6. Daniel Charytonowicz, Rachel Brody, Robert Sebra. Interpretable and context-free deconvolution of multi-scale whole transcriptomic data with UniCell deconvolve. Nature Communications, 2023. https://doi.org/10.1038/s41467-023-36961-8 ↩

7. Jiao Sun, Ayesha A. Malik, Tong Lin, et al. MOADE: a multimodal autoencoder for dissociating bulk multi-omics data. Genome Biology, 2025. https://doi.org/10.1186/s13059-025-03805-1 ↩

8. Maik Wolfram-Schauerte, Thomas Vogel, Hanati Tuoken, et al. Approaching the holistic transcriptome—convolution and deconvolution in transcriptomics. Briefings in Bioinformatics, 2025. https://doi.org/10.1093/bib/bbaf388 ↩

9. Kavoos Momeni, Saeid Ghorbian, Ehsan Ahmadpour, et al. Unraveling the complexity: understanding the deconvolutions of RNA-seq data. Translational Medicine Communications, 2023. https://doi.org/10.1186/s41231-023-00154-8 ↩