Immune Deconvolution

Immune Deconvolution

Estimate immune and stromal cell composition from bulk RNA-seq using multiple algorithms. Wraps CIBERSORT, quanTIseq, EPIC, xCell, MCP-counter, TIMER, and ESTIMATE through the immunedeconv unified interface.

When to Use This Skill

Activate when the user requests:

• Immune cell type estimation from bulk RNA-seq or microarray
• Tumor microenvironment characterization
• Immune subtype classification across a cancer cohort
• Tumor purity estimation from expression data
• Comparison of immune infiltration between conditions or subtypes
• Multi-method deconvolution benchmarking

Inputs

All methods except TIMER and ESTIMATE need TPM with HGNC gene symbols as rownames. Raw counts and Ensembl IDs will produce wrong results silently.

Preparing Input from TCGA

library(TCGAbiolinks)
library(SummarizedExperiment)

query <- GDCquery(
  project = "TCGA-BRCA",
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification",
  workflow.type = "STAR - Counts"
)
GDCdownload(query, directory = "GDCdata")
se <- GDCprepare(query, directory = "GDCdata")

# Extract TPM — immunedeconv needs TPM, not counts
tpm <- assay(se, "tpm_unstrand")

# Convert Ensembl IDs to gene symbols
library(org.Hs.eg.db)
symbols <- mapIds(org.Hs.eg.db,
  keys = sub("\\..*", "", rownames(tpm)),
  keytype = "ENSEMBL", column = "SYMBOL", multiVals = "first")

# Drop unmapped genes, resolve duplicates by keeping highest mean expression
tpm <- tpm[!is.na(symbols), ]
rownames(tpm) <- symbols[!is.na(symbols)]

# Deduplicate: quanTIseq crashes on duplicate rownames
dups <- duplicated(rownames(tpm))
if (any(dups)) {
  means <- rowMeans(tpm)
  keep <- !duplicated(rownames(tpm)) |
    (duplicated(rownames(tpm)) & means == ave(means, rownames(tpm), FUN = max))
  tpm <- tpm[keep, ]
  tpm <- tpm[!duplicated(rownames(tpm)), ]  # safety net
}

# Keep only tumor samples for deconvolution
tumor_barcodes <- colData(se)$sample_type == "Primary Tumor"
tpm_tumor <- tpm[, tumor_barcodes]

Method Selection

What do you need from the deconvolution?

Absolute cell fractions (values sum to 1, interpretable as proportions)?
  -> quanTIseq: 10 immune types + "Other". Validated against flow cytometry and IHC.
  -> EPIC: 6 immune types + cancer cells. Built for solid tumors.
  -> CIBERSORT absolute mode: 22 immune subtypes. Requires registration.

Relative proportions (which immune cells dominate, not how much)?
  -> CIBERSORT relative: 22 types. Gold standard, widely published.
  -> quanTIseq also works (ignore "Other", compare immune fractions).

Enrichment scores (rank samples by infiltration, not quantify)?
  -> xCell: 64 cell types. Broadest coverage, good for discovery.
  -> MCP-counter: 8 immune + 2 stromal. Fewest assumptions, works well with noisy data.

Immune vs stromal vs tumor purity?
  -> ESTIMATE: immune score, stromal score, purity estimate.
  -> Or use "Other" fraction from quanTIseq/EPIC as purity proxy.

Cancer-type-specific correction?
  -> TIMER: 6 immune types, pre-built models per TCGA cancer type.
  -> ConsensusTME: cancer-specific gene sets, uses ssGSEA.

Don't know which to pick?
  -> Run quanTIseq + MCP-counter + EPIC. If all three agree on a trend,
     the signal is real. Disagreement means the effect is method-dependent.

Score Interpretation

This matters more than method choice. Getting it wrong leads to wrong biological conclusions.

xCell and MCP-counter scores are arbitrary units. A CD8 score of 0.4 from xCell and 3.7 from MCP-counter cannot be compared to each other, and neither means "40% CD8 T cells."

immunedeconv Setup

# Install: wraps 9 methods in one API
# remotes::install_github("omnideconv/immunedeconv")
library(immunedeconv)  # v2.1.0+

# For CIBERSORT (optional, requires registration at cibersortx.stanford.edu):
# Download CIBERSORT.R and LM22.txt, then:
# set_cibersort_binary("/path/to/CIBERSORT.R")
# set_cibersort_mat("/path/to/LM22.txt")

quanTIseq

qt <- deconvolute(tpm_tumor, method = "quantiseq")
# Returns tibble: cell_type column + one column per sample
# Cell types: B cell, Macrophage M1, Macrophage M2, Monocyte, Neutrophil,
#   NK cell, T cell CD4+ (non-regulatory), T cell CD8+,
#   T cell regulatory (Tregs), Myeloid dendritic cell,
#   uncharacterized cell

# Fractions sum to 1 per sample. "uncharacterized cell" = tumor + stroma.
# In solid tumors, "uncharacterized cell" should be 0.5-0.9.
# If < 0.3, check if sample is actually a lymphoma or immune-rich tissue.

EPIC

epic <- deconvolute(tpm_tumor, method = "epic")
# Cell types: B cell, CAFs, CD4+ T cell, CD8+ T cell,
#   Endothelial, Macrophage, NK cell, otherCells

# otherCells = uncharacterized (mostly tumor in solid cancers)
# EPIC was specifically built for solid tumors — preferred over
# quanTIseq when you need a cancer cell fraction estimate

CIBERSORT

# Requires setup (see above). LM22 has 22 immune subtypes:
# naive/memory B cells, plasma, CD8/CD4 subtypes, gamma-delta T,
# NK resting/activated, monocytes, M0/M1/M2 macrophages,
# DC resting/activated, mast resting/activated, eosinophils, neutrophils

cib <- deconvolute(tpm_tumor, method = "cibersort")
# Relative fractions (sum to 1 across 22 types)

cib_abs <- deconvolute(tpm_tumor, method = "cibersort_abs")
# Absolute mode: scales by total immune content

# CIBERSORT reports a p-value per sample from permutation testing.
# p > 0.05 means the deconvolution fit is unreliable — flag or exclude.
# Access via running CIBERSORT directly (not through immunedeconv wrapper).

xCell

xc <- deconvolute(tpm_tumor, method = "xcell")
# 64 cell types including fine subtypes:
#   T cell subtypes (Th1, Th2, Th17, Treg, CD8 naive, CD8 effector, etc.)
#   B cell subtypes (naive, memory, plasma)
#   Myeloid (M1, M2, DC subsets)
#   Stromal (fibroblasts, endothelial, pericytes)
# Plus composite: ImmuneScore, StromalScore, MicroenvironmentScore

# xCell uses ssGSEA followed by spillover correction.
# Scores are enrichment-based — only compare within a cell type across samples.
# All-zero scores for a cell type = not enough gene signature overlap.
# Check that your gene symbols are standard HGNC.

MCP-counter

mcp <- deconvolute(tpm_tumor, method = "mcp_counter")
# 10 populations: T cells, CD8+ T cells, Cytotoxic lymphocytes,
#   B lineage, NK cells, Monocytic lineage, Myeloid dendritic cells,
#   Neutrophils, Endothelial cells, Fibroblasts

# Scores are arbitrary units based on marker gene expression.
# Uses manually curated transcriptomic markers with minimal cell-type overlap
# with minimal overlap between cell types.
# Works well even with noisy or poorly normalized data.

TIMER

# TIMER uses tumor-type-specific regression — requires specifying cancer type
timer <- deconvolute(tpm_tumor, method = "timer",
  indications = rep("brca", ncol(tpm_tumor)))
# Valid indications: lowercase TCGA codes (brca, luad, coad, gbm, skcm, etc.)

# Returns 6 types: B cell, CD4+ T cell, CD8+ T cell,
#   Neutrophil, Macrophage, Dendritic cell

# If your cancer type is not in TIMER's pre-built models, use another method.
# TIMER needs at least 2 samples — single-sample deconvolution will fail.

ESTIMATE

# ESTIMATE: stromal score, immune score, tumor purity
# Oldest method (2013) but still the standard for purity estimation

est <- deconvolute_estimate(tpm_tumor)
# Returns: StromalScore, ImmuneScore, ESTIMATEScore, TumorPurity

# TumorPurity is derived from: cos(0.6049872018 + 0.0001467884 * ESTIMATEScore)
# This conversion was calibrated on Agilent microarrays.
# For RNA-seq, treat purity as approximate — the rank ordering is reliable,
# the absolute values less so.

Tumor Purity Correction

Why purity matters:
  Low purity sample = fewer tumor cells = immune cells are a larger fraction
  High purity sample = more tumor cells = immune signal is compressed

  If group A has lower purity than group B on average, deconvolution
  will show group A as "more immune" even if the actual immune
  infiltration per unit tissue is identical.

When to correct:
  Comparing immune infiltration between groups
    -> Include ESTIMATE purity as covariate in linear model
  Correlating deconvolution scores with mutations or expression
    -> Partial correlation controlling for purity
  Filtering bad samples
    -> Remove samples with ESTIMATE purity < 0.2 (likely failed or lymphomas)

When NOT to correct:
  Reporting absolute immune fractions from quanTIseq/EPIC
    -> These already account for non-immune content via "Other"/"otherCells"
  Comparing within the same tumor type with similar purity distributions
    -> Correction adds noise without removing bias

# Purity as covariate
est_scores <- deconvolute_estimate(tpm_tumor)
purity <- as.numeric(est_scores["TumorPurity", ])
sample_info$purity <- purity

# Linear model with purity adjustment
cd8_scores <- as.numeric(mcp[mcp$cell_type == "T cell CD8+", -1])
fit <- lm(cd8_scores ~ subtype + purity, data = sample_info)
summary(fit)
# Subtype coefficients are now purity-adjusted

# Partial correlation
library(ppcor)
pcor_result <- pcor.test(
  x = cd8_scores,
  y = sample_info$mutation_load,
  z = purity,
  method = "spearman"
)

Cross-Method Comparison

# Run multiple methods and compare CD8+ T cell estimates
methods_to_run <- c("quantiseq", "epic", "mcp_counter", "xcell")
all_res <- lapply(methods_to_run, function(m) {
  deconvolute(tpm_tumor, method = m)
})
names(all_res) <- methods_to_run

# Extract CD8 scores — each method uses slightly different cell type names
cd8_names <- c(
  quantiseq  = "T cell CD8+",
  epic       = "T cell CD8+",
  mcp_counter = "T cell CD8+",
  xcell      = "T cell CD8+"
)

cd8_mat <- sapply(methods_to_run, function(m) {
  res <- all_res[[m]]
  as.numeric(res[res$cell_type == cd8_names[m], -1])
})
rownames(cd8_mat) <- colnames(tpm_tumor)

# Cross-method Spearman correlation
cor(cd8_mat, method = "spearman", use = "pairwise.complete.obs")
# Expect rho > 0.5 between most pairs for CD8+ T cells
# rho < 0.3 for most pairs: check input format (log-transformed? wrong gene IDs?)

Visualization

library(ggplot2)

# --- Stacked bar plot (quanTIseq fractions) ---
qt_long <- tidyr::pivot_longer(qt, -cell_type,
  names_to = "sample", values_to = "fraction")

ggplot(qt_long, aes(x = sample, y = fraction, fill = cell_type)) +
  geom_bar(stat = "identity", position = "stack") +
  theme_classic(base_size = 12) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  labs(y = "Cell fraction", x = "Samples", fill = "Cell type")

# --- Boxplot by molecular subtype ---
score_df <- data.frame(
  sample = colnames(tpm_tumor),
  CD8 = as.numeric(mcp[mcp$cell_type == "T cell CD8+", -1]),
  subtype = sample_info$subtype
)

ggplot(score_df, aes(x = reorder(subtype, CD8, FUN = median), y = CD8, fill = subtype)) +
  geom_boxplot(outlier.size = 0.5) +
  geom_jitter(width = 0.15, alpha = 0.3, size = 0.5) +
  theme_classic(base_size = 12) +
  labs(y = "MCP-counter CD8+ T cell score", x = NULL) +
  theme(legend.position = "none")

# --- Cross-method correlation heatmap ---
library(ComplexHeatmap); library(circlize)
cor_mat <- cor(cd8_mat, method = "spearman", use = "pairwise.complete.obs")
col_fun <- colorRamp2(c(0, 0.5, 1), c("#2166AC", "white", "#B2182B"))

Heatmap(cor_mat, name = "Spearman\nrho", col = col_fun,
  cell_fun = function(j, i, x, y, w, h, fill) {
    grid.text(sprintf("%.2f", cor_mat[i, j]), x, y, gp = gpar(fontsize = 10))
  },
  column_title = "Cross-method concordance: CD8+ T cells")

# --- ESTIMATE purity vs immune score scatter ---
est_df <- data.frame(
  purity = as.numeric(est_scores["TumorPurity", ]),
  immune = as.numeric(est_scores["ImmuneScore", ]),
  subtype = sample_info$subtype
)
ggplot(est_df, aes(purity, immune, color = subtype)) +
  geom_point(alpha = 0.5, size = 1.5) +
  theme_classic(base_size = 12) +
  labs(x = "ESTIMATE tumor purity", y = "ESTIMATE immune score")

Output Specification

Validation Checks

After running deconvolution on TCGA-BRCA, verify:

Subtype ordering (consistent across methods):
  Basal-like should have highest CD8+ T cell scores
  Luminal A should have lowest overall immune infiltration
  HER2-enriched sits between Basal and Luminal B
  If Luminal A ranks highest for immune infiltration, something is wrong.

quanTIseq sanity:
  "uncharacterized cell" in solid tumors should be 0.5-0.9
  Sum of all fractions = 1.0 per sample (within floating point tolerance)
  Macrophage M2 typically > M1 in most solid tumors
  If all CD8 fractions are 0: gene symbols may not match quanTIseq signature

ESTIMATE sanity:
  TCGA-BRCA median purity: ~0.75
  Basal: lower purity (0.55-0.70), Luminal A: higher purity (0.80-0.90)
  If all purities cluster near 1.0 or near 0: wrong platform setting or wrong input scale

Cross-method concordance:
  CD8+ T cell Spearman rho > 0.5 between quanTIseq, EPIC, MCP-counter
  rho < 0.3 across most pairs: input is likely log-transformed or has wrong gene IDs

Red flags:
  CIBERSORT p > 0.05 for > 30% of samples: signature doesn't fit this cohort
  xCell returns all zeros for many cell types: gene overlap too low
  TIMER returns NA: cancer type not in the pre-built model set
  quanTIseq crashes: duplicate gene symbols in rownames (must deduplicate)
  Only 1 sample: TIMER and ConsensusTME need >= 2 samples per cancer type

Common Pitfalls

Input

1. Raw counts instead of TPM: quanTIseq, EPIC, xCell, MCP-counter, and CIBERSORT expect TPM. Feeding raw counts compresses low-expression signatures and inflates high-expression ones. Convert first: tpm <- counts / gene_lengths * 1e6 / colSums(counts / gene_lengths).
2. Log-transformed input: Most methods expect linear-scale TPM. Passing log2(TPM+1) flattens the dynamic range and produces wrong fractions. The immunedeconv docs say explicitly: "data must be on non-log scale."
3. Ensembl IDs as rownames: Every signature matrix uses HGNC gene symbols. Ensembl IDs will match nothing, and most methods fail silently with near-zero results rather than throwing an error.
4. Duplicate gene symbols: quanTIseq crashes on duplicate rownames. Deduplicate before calling any method — keep the row with highest mean expression.

Interpretation

5. Treating scores as fractions: xCell and MCP-counter output enrichment scores, not proportions. Saying "20% CD8 T cells" from an MCP-counter score of 3.7 is meaningless. Only quanTIseq, EPIC, and CIBERSORT return interpretable fractions.
6. Comparing values across methods: A CD8 estimate of 0.12 from quanTIseq and 4.8 from MCP-counter are on completely different scales. Compare ranks or standardized scores, never raw values.
7. Ignoring tumor purity: A low-purity sample will show high immune fractions simply because the denominator (total cells) contains fewer tumor cells. Control for purity when comparing groups with different purity distributions.
8. Treating deconvolution as ground truth: These are estimates. Cross-method concordance increases confidence; single-method results are hypotheses. For publication, validate key findings with IHC, flow cytometry, or spatial transcriptomics.

Method-specific

9. CIBERSORT without checking fit p-values: Samples with p > 0.05 have unreliable estimates. Excluding these is standard practice but often forgotten.
10. TIMER on unsupported cancer types: TIMER's regression models are pre-built for specific TCGA indications. Running it on a cancer type not in the training set gives uninterpretable results.
11. Single-sample runs with TIMER or ConsensusTME: These methods need at least 2 samples of the same cancer type. A single sample will either error or return NAs.
12. Mixing microarray and RNA-seq without batch correction: Deconvolution signatures are sensitive to platform effects. Don't pool microarray and RNA-seq samples without explicit normalization.

Advanced: Single-Cell Reference-Based Deconvolution

When matched scRNA-seq data is available for the same tissue type, second-generation methods outperform signature-based approaches.

# BayesPrism: top performer in 2024-2025 benchmarks
# install.packages("BayesPrism")  # or InstaPrism for faster runtime
library(BayesPrism)

# Requires: scRNA-seq reference (genes x cells) + cell type labels
# bulk: genes x samples matrix (counts, not TPM)
bp <- new.prism(
  reference = sc_counts,        # raw counts from scRNA-seq reference
  mixture = bulk_counts,         # raw counts from bulk samples
  input.type = "count.matrix",
  cell.type.labels = cell_labels,
  cell.state.labels = cell_states  # finer subtypes, optional
)
result <- run.prism(bp, n.cores = 4)
theta <- get.fraction(result, which = "final")  # samples x cell_types
# BayesPrism returns posterior mean fractions that sum to 1

When to use scRNA-seq-based methods vs signature-based:
  Have matched scRNA-seq reference for the same tissue type?
    -> BayesPrism or DWLS. Better accuracy than fixed signatures.
  No scRNA-seq reference, just bulk RNA-seq?
    -> quanTIseq + EPIC + MCP-counter (signature-based).
  Public scRNA-seq atlas exists for the tissue?
    -> Consider BayesPrism with the atlas as reference, but be aware
       of batch effects between your bulk data and the atlas.

Related Skills

• cancer-multiomics: TCGA data retrieval and expression matrix preparation (input for deconvolution)
• survival-analysis: Use deconvolution scores as covariates in Cox models or KM stratification

Public Datasets for Testing