Inferring Copy Number Variation (inferCNV)

This tutorial is a prepared-input template: it validates the three inferCNV inputs, constructs an R analysis object, and guides inspection of expression patterns without treating them as DNA copy-number measurements.

Learning objectives and prerequisites

By the end of this lesson you will be able to:

• Validate prepared raw counts, cell annotations, and assembly-matched gene ordering, including the degree of gene overlap.
• Create an inferCNV object and run the core pipeline in R with reproducible commands.
• Locate and critically inspect the run heatmap and available outputs, while understanding caveats that prevent over-interpretation.
• Recognize common failure modes and apply practical troubleshooting steps.

Prerequisites:

• Completed basic single-cell RNA-seq (scRNA-seq) workflows: quality control, normalization, integration, clustering, and cell-type annotation.
• Familiarity with R (>= 4.0 recommended) and a working install of the inferCNV package (see References).
• Access to non-malignant reference cells within the same sample or a matched normal sample; using unrelated clinical samples as references is discouraged.

Context: why infer CNVs from scRNA-seq?

Copy number variation (CNV) refers to gains and losses of chromosomal segments and is a hallmark of many cancers. Direct measurement of CNV requires DNA-based assays (e.g., whole-genome sequencing, shallow WGS, or single-cell DNA-seq). However, broad CNV events often change the aggregate expression of genes across genomic regions. inferCNV is an R package that leverages single-cell transcriptomes to infer large-scale CNV patterns by comparing expression in "query" cells to a set of "reference" (presumed diploid) cells. This makes inferCNV useful for prioritising putative malignant or aneuploid populations in tumors, or for exploring chromosomal perturbations in heterogeneous tissues.

Important conceptual points:

• inferCNV infers broad-scale (megabase-level) copy number changes rather than precise small focal events.
• It is an expression-based, not DNA-based, method: differences in expression programs, cell cycle, or sample composition can create patterns that mimic CNV.
• Reference selection is the most critical factor: if the reference set contains cells with copy-number lesions, inferences will be biased.

Reproducible workflow

This is a run template, not an end-to-end reproducible dataset: it requires three prepared files and does not download or construct them. Before running, set the working directory to the folder containing those files, verify the identifier and genome-assembly contract below, and use the official file definitions to prepare missing inputs. The commands read inputs and write files in that directory.

Check the installed R and package versions before analysis with R --version and, in R, packageVersion("infercnv"). The examples assume raw, unnormalized gene-by-cell counts; do not substitute normalized or scaled expression. Record these versions with the assembly and identifier mapping.

Install and verify the JAGS prerequisite

The R package needs native JAGS 4.x as well as the R binding rjags. Installing R packages alone does not install JAGS. With R and its build tools already configured, Ubuntu 24.04 users can install JAGS in a terminal:

sudo apt-get update
sudo apt-get install jags

On other operating systems, use the official installation instructions to install JAGS first. Then run this separate setup step in the R environment you will use for analysis. Bioconductor selects the release matching your R version; update=FALSE avoids a blanket update of existing packages.

# Setup only: do this before opening the analysis script below.
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager", repos = "https://cloud.r-project.org")
}
BiocManager::install("infercnv", ask = FALSE, update = FALSE)
needed <- c("rjags", "infercnv")
missing <- needed[!vapply(needed, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing)) stop("Setup incomplete: ", paste(missing, collapse = ", "))
message("JAGS: ", as.character(rjags::jags.version()))
message("infercnv: ", as.character(packageVersion("infercnv")))

Proceed only when both packages load and the native JAGS version prints. Package installation can warn about failures without stopping R, so do not treat the return from the installer as proof of readiness. Keep setup and analysis in the same R library/environment.

Verified execution and limits

The original R route was executed on Ubuntu 24.04.4 with R 4.3.3, infercnv 1.18.1 and JAGS 4.3.2 using 170 genuine public cells from a CELLxGENE-curated GSE84465 subset and 20,374 exact GRCh37 gene matches. The run retained 14,319 genes and all 170 cells, with HMM and denoising enabled. The returned object was reloaded in a fresh R process; gene coordinates, cell identities, HMM report links and opaque white-background heatmaps were checked. Reference labels, input paths and a two-thread limit were adapted for this bounded test.

This is software-execution evidence, not DNA or clinical validation. The input is a Smart-seq2 GBM subset: cutoff=0.1 was retained to test the displayed template, not established as the appropriate threshold for these full-length data. Choose the assay-appropriate cutoff before your own analysis. Candidate reference labels were not independently verified as diploid. Setup reused existing R libraries, so this was not a fresh whole-stack installation or a test of inferCNV on other operating systems or hardware.

Required files (inputs)

• Raw counts matrix in R format or as an on-disk matrix: counts_matrix.rds (genes x cells). Name rows with consistent gene identifiers, such as approved HGNC symbols or stable Ensembl IDs; use exactly matching identifiers in the gene-ordering file.
• Cell annotations: a table with cell_id and cell_type columns, one row per cell. IDs must match the count-matrix column names; the package input written below is headerless.
• Gene ordering: gene_ordering.tsv with gene, chromosome, start, and end coordinates from a declared assembly (the example metadata table has headers). Use the same gene identifiers as the count-matrix row names and document the source, assembly, and any identifier conversion; only overlapping genes are analyzed.

Example R session

# Start R in the directory containing the three prepared input files.
# Do not install inside this analysis script; follow official Bioconductor
# installation instructions, then record R.version.string and packageVersion().
stopifnot(requireNamespace("infercnv", quietly = TRUE))
message("R: ", R.version.string, "; infercnv: ", as.character(packageVersion("infercnv")))
library(infercnv)

# 1) Load inputs. Replace filenames as needed; record assembly and mapping.
# gene_ordering.tsv coordinates must match the identifiers/assembly of counts.
counts <- readRDS("counts_matrix.rds")             # gene x cell matrix (raw counts)
annotations <- read.delim("cell_annotations.tsv", header=TRUE, stringsAsFactors=FALSE)
gene_order <- read.delim("gene_ordering.tsv", header=TRUE, stringsAsFactors=FALSE)

# Check cell IDs and coordinate metadata before writing headerless package inputs.
stopifnot(!anyDuplicated(colnames(counts)), !anyDuplicated(annotations$cell_id))
stopifnot(setequal(colnames(counts), annotations$cell_id))
stopifnot(all(c("cell_id", "cell_type") %in% names(annotations)))
stopifnot(!anyNA(annotations$cell_type), all(nzchar(annotations$cell_type)))
stopifnot(is.matrix(counts) || inherits(counts, "Matrix"))
stopifnot(all(is.finite(counts)), all(counts >= 0), all(counts == floor(counts)))
stopifnot(all(c("gene_symbol", "chromosome", "start", "end") %in% names(gene_order)))
stopifnot(!anyDuplicated(gene_order$gene_symbol))
stopifnot(!anyNA(gene_order$gene_symbol), all(nzchar(gene_order$gene_symbol)))
stopifnot(!anyNA(gene_order$chromosome), all(nzchar(gene_order$chromosome)))
stopifnot(all(is.finite(gene_order$start)), all(is.finite(gene_order$end)))
stopifnot(all(gene_order$start >= 1), all(gene_order$start <= gene_order$end))
matched_genes <- intersect(rownames(counts), gene_order$gene_symbol)
coverage <- length(matched_genes) / nrow(counts)
message(sprintf("Gene-order overlap: %d / %d count genes (%.1f%%)",
                length(matched_genes), nrow(counts), 100 * coverage))
if (length(matched_genes) == 0L) stop("No shared genes: inspect identifiers, mapping, and assembly before proceeding")
unmatched_genes <- setdiff(rownames(counts), matched_genes)
if (length(unmatched_genes) > 0L) {
    print(utils::head(unmatched_genes, 10))
    warning("Unmatched count genes remain: review the reported IDs and document exclusions before analysis")
}
# Review chromosome distribution and why genes are absent before continuing.
# Neither the overlap percentage nor a fixed gene count establishes adequate
# genomic coverage or biological validity; no universal cutoff is imposed here.

# 2) Write annotation file required by inferCNV (if not already in file form)
# inferCNV expects cell_id and group, without a header row.
write.table(annotations[,c("cell_id","cell_type")],
            file="infercnv_annotations.tsv", sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)

# 3) Write gene ordering file in inferCNV expected format:
# headerless columns: gene, chromosome, start, end
write.table(gene_order[,c("gene_symbol","chromosome","start","end")],
            file="infercnv_gene_ordering.tsv", sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)

# 4) Create the inferCNV object.
# Choose reference groups: names in 'cell_type' that represent normal/diploid cells
reference_groups <- c("T_cell","B_cell","Endothelial")  # example; adjust to your annotations

infercnv_obj <- CreateInfercnvObject(
    raw_counts_matrix = counts,
    annotations_file = "infercnv_annotations.tsv",
    delim = "\t",
    gene_order_file = "infercnv_gene_ordering.tsv",
    ref_group_names = reference_groups
)

# 5) Run inferCNV.
# - cutoff: 0.1 for droplet-based 10x (sparse); 1 for Smart-seq2-like full-length counts.
# - out_dir: path to write outputs (creates directory if needed).
out_dir <- "infercnv_output"
if (dir.exists(out_dir) && length(list.files(out_dir, all.files = TRUE, no.. = TRUE))) {
    stop("Choose an empty output directory for a new run; do not mix old and new results")
}
dir.create(out_dir, showWarnings = FALSE)

infercnv_obj <- infercnv::run(
    infercnv_obj,
    cutoff = 0.1,
    out_dir = out_dir,
    cluster_by_groups = TRUE,
    denoise = TRUE,
    HMM = TRUE
)

# 6) Save the returned object and verify it can be read back.
final_path <- file.path(out_dir, "returned_final_infercnv.rds")
saveRDS(infercnv_obj, final_path)
final_check <- readRDS(final_path)
stopifnot(methods::is(final_check, "infercnv"))
methods::validObject(final_check)
stopifnot(nrow(final_check@expr.data) > 0L, ncol(final_check@expr.data) > 0L)
stopifnot(all(is.finite(final_check@expr.data)))
stopifnot(all(colnames(final_check@expr.data) %in% colnames(counts)))
stopifnot(identical(rownames(final_check@expr.data), rownames(final_check@gene_order)))
message("Retained ", nrow(final_check@expr.data), " genes and ",
        ncol(final_check@expr.data), " cells; saved object: ", final_path)

# 7) HMM=TRUE must produce readable HMM region reports, not just a heatmap.
# Prefixes vary by release; the documented region-report suffix is used here.
hmm_reports <- list.files(out_dir, pattern = "\\.pred_cnv_regions\\.dat$",
                          full.names = TRUE, recursive = TRUE)
if (!length(hmm_reports)) stop("No HMM region report: inspect the run log and installed-version output documentation")
hmm_type <- final_check@options$HMM_type
if (length(hmm_type) != 1L || !hmm_type %in% c("i6", "i3")) {
    stop("Unrecognized HMM type: review installed-version documentation before checking state codes")
}
legal_states <- if (hmm_type == "i6") 1:6 else 1:3
for (path in hmm_reports) {
    report <- read.delim(path, header = TRUE, check.names = FALSE)
    stopifnot(all(c("cell_group_name", "cnv_name", "state", "chr", "start", "end") %in% names(report)))
    if (nrow(report)) {
        stopifnot(all(report$state %in% legal_states))
        stopifnot(all(is.finite(report$start)), all(is.finite(report$end)))
        stopifnot(all(report$start >= 1), all(report$start <= report$end))
    }
    message(basename(path), ": ", nrow(report), " reported regions")
}
writeLines(capture.output(sessionInfo()), file.path(out_dir, "session-info.txt"))
print(list.files(out_dir, recursive = TRUE))
# Inventory is not proof of completion. Inspect the heatmap and reports below.

Notes on parameters:

• cutoff: filters genes by mean counts across cells. The Broad workflow documentation suggests 0.1 for sparse droplet counts and 1 for Smart-seq-like data; treat these as starting points, not universal thresholds.
• denoise: enables inferCNV denoising procedures. Inspect the preliminary and processed views; whether denoising helps depends on the data and chosen filters.
• HMM: runs an optional Hidden Markov Model to assign discrete copy-number states. Treat HMM outputs as smoothed calls, not definitive DNA copy numbers.

Alternative: infercnvpy in Python

# infercnvpy is a distinct Python implementation. Unlike the R workflow above,
# its documented input is normalized, log-transformed expression.
# This is an API sketch, not a runnable preparation recipe: adata must already
# contain filtered cells and normalized, log-transformed X, coordinate columns
# in adata.var, and cell labels in adata.obs. See the official input guide:
# https://infercnvpy.readthedocs.io/en/latest/infercnv.html

import infercnvpy as cnv

# Inputs:
# adata: filtered AnnData with normalized, log-transformed expression in X
# reference_cat: list of cell-type names to use as diploid reference
cnv.tl.infercnv(
    adata,
    reference_key="cell_type",
    reference_cat=["T_cell", "B_cell", "Endothelial"],
    window_size=100
)

# With inplace=True (the default), results are stored at adata.obsm["X_cnv"].
# Confirm the installed version and coordinate metadata before downstream use.
cnv_matrix = adata.obsm["X_cnv"]

Output inspection and interpretation

Primary outputs to inspect in the inferCNV output directory:

• The documented final heatmap is commonly named infercnv.png; the preliminary view is infercnv.preliminary.png. Read the plot labels and package version when interpreting axes and colors.
• If HMM prediction is enabled, inspect its generated files using the installed version documentation. These are model-derived states, not direct measurements of DNA dosage.
• InferCNV also saves stage-specific R objects and text matrices; consult the output directory and package documentation rather than relying on an assumed filename.

The checkpoint saves the object returned by this run rather than guessing which stage file is final. It also checks the HMM region-table schema and the state codes for the recorded HMM type. A report with zero retained regions can be a valid outcome after neutral-state or posterior-probability filtering: it is not proof that every cell is diploid. Inspect the unfiltered and filtered reports and the run log. If another release changes the output schema, consult that release's documentation instead of renaming files or inventing calls to satisfy the check. This template uses a new, empty output folder so stale artifacts cannot masquerade as new results; deliberate checkpoint resumption needs separate review of the package's resume options.

What to look for:

• Look for broad, contiguous expression shifts across genes in genomic order. A regional pattern is a candidate signal to investigate, not a CNV call by itself.
• Ask whether a broad pattern recurs across cells and is compatible with the reference and sample design. Shared expression patterns can still arise from cell state or other confounding, so they do not establish clonality.
• Compare HMM states with continuous heatmap patterns where both are present. Agreement is a consistency check within related model outputs, not independent validation; disagreement warrants review of assumptions and parameters.

What a plausible positive result means:

• Coordinated amplification or deletion across many genes in a contiguous genomic region is consistent with a copy-number lesion and can be used to prioritise cells/clusters for DNA validation.
• Shared CNV patterns across cells that also express malignant markers strengthens the hypothesis of malignancy.

What the result does NOT prove:

• An expression-based CNV call is not a diagnostic DNA measurement. Confounding factors (cell cycle, strong transcriptional programs that co-regulate nearby genes, or an inappropriate reference) can produce false positives.
• Small focal events and balanced rearrangements are outside the method’s reliable scope; use an appropriate DNA assay to assess them.
• Quantitative copy number (exact integer ploidy) should not be assumed from inferred signal alone without orthogonal DNA data.

Troubleshooting

Below are common symptoms, likely causes, and practical fixes.

Symptom: No clear signal (noisy or flat heatmap)

• Possible causes:
  • Reference groups contain contaminated or aneuploid cells (non-diploid).
  • Too aggressive filtering of genes removed informative genes.
  • Incorrect gene ordering or mismatched genome build (e.g., gene symbols matched to different coordinates).
  • cutoff parameter too high for sparse data.
• Fixes:
  • Confirm reference groups are high-confidence normal cells (e.g., stromal or immune cells not expected to be malignant). If unsure, restrict to well-characterised cell subsets.
  • Check the gene_ordering file for completeness and that gene symbols match rownames(counts). Use Ensembl-to-HGNC mapping if necessary.
  • Check whether the chosen cutoff excludes too many genes, using the installed version documentation as a guide. A run with HMM=FALSE can help separate continuous expression patterns from discrete HMM states; it does not guarantee improved signal.

Symptom: Many small fragmented gains and losses scattered across genome

• Possible causes:
  • Batch effects or strong cell-type-specific transcriptional programs dominate the signal.
  • Over-smoothing window size or inappropriate denoising parameters producing spurious local excursions.
  • Counts were normalized prior to inferCNV (inferCNV expects raw counts).
• Fixes:
  • For the R inferCNV workflow shown here, use raw counts. Do not pass log-normalized or scaled values as raw counts.
  • cluster_by_groups=TRUE organizes cells by annotation for visualization; it does not correct batch effects or validate a regional signal.
  • Experiment with disabling HMM and examining continuous signal; adjust denoise and smoothing parameters.

Symptom: Reference appears shifted relative to expectation

• Possible causes: reference set contains a subpopulation with true CNV; reference composition differs across samples.
• Fixes: build a stricter reference set (e.g., only well-characterised immune subsets), or use matched normal tissue when available. Run inferCNV with reference groups from the same sample where possible.

Limitations, version and reference-data caveats, and responsible use

Limitations and caveats:

• Genome build consistency: gene coordinates in the gene_order file must match the gene symbols and the genome assembly you intend to represent (e.g., hg19/GRCh37 vs hg38/GRCh38). Mismatches will produce incorrect genomic ordering and misleading heatmaps.
• Resolution depends on gene coverage, expression, smoothing, and data quality. This expression-based approach is intended for broad patterns, not single-exon or small focal CNV detection.
• Noise model: sequencing depth, protocol differences (10x vs Smart-seq2), and sample composition influence sensitivity. Parameter choices (cutoff, denoise) must reflect data type.
• Reference fidelity: a diploid reference is assumed. If references are not diploid, inferred relative shifts may be inverted or dampened.

Version and reproducibility:

• Different inferCNV releases may change default parameters and output formats. Record the package version (e.g., sessionInfo() in R) and document genome build and gene ID mapping used.
• Keep the input raw counts, annotation file, and gene_ordering file in an archival folder alongside the output directory for reproducibility.

Responsible use and interpretation:

• Treat inferCNV results as hypothesis-generating. Any claim about malignancy or diagnostics should rely on orthogonal DNA-based validation and expert pathology review.
• Do not report expression-inferred CNVs as clinical diagnoses. Label results explicitly as inferred from transcriptome data in any downstream communication.

Next-step checklist

1. Confirm raw counts and gene-ordering are consistent (matching gene symbols and genome build).
2. Define and document reference groups; re-check that these cells are likely diploid.
3. Run inferCNV with denoise=TRUE and appropriate cutoff for your platform; inspect both continuous heatmap and HMM outputs.
4. Compare CNV patterns to cluster annotations and known marker genes to prioritise candidate malignant clusters.
5. Validate prioritized candidates with orthogonal DNA assays (shallow WGS, targeted DNA panels, FISH, or single-cell DNA-seq) before making biological or clinical claims.

References

• inferCNV GitHub repository (Broad Institute)
• Broad Institute inferCNV: running the analysis and documented outputs
• Bioconductor: infercnv package and current release documentation
• infercnvpy API reference and input preparation and method notes

Record the R and inferCNV versions, genome assembly, gene identifier mapping, reference groups, and parameter values with each run.