Run HMMcopy on per-window reads for a single sample

Given a table of per-window counts with GC and mappability, run HMMcopy correction and segmentation, and return a GRanges of CNV calls for the selected sample/column.

runHMMCopy(CNV_RegionsWithReads, colname, plotDir = NULL)

Arguments

CNV_RegionsWithReads

data.frame / data.table / GRanges. Coercible table with columns: chr, start, end, width, strand, gc, map, and one sample-specific count column referenced by colname.

colname

character(1). Name of the sample count column to use (e.g., "sampleA").

plotDir

character(1) or NULL. Directory to write a per-sample PDF diagnostic plot named <colname>.pdf; if NULL, no plot is produced. Default: NULL.

Value

A GRanges with a metadata column named <colname> containing per-window CNV estimates (typically HMM state means on the log2 scale).

Details

Counts are GC/map corrected via HMMcopy::correctReadcount() and segmented with HMMcopy::HMMsegment(). The result is mapped back to the input windows and returned as a GRanges. Ensure that the binning (width) and covariate columns (gc, map) are consistent across the table.

See also

addHMMcopyCNV(), HMMcopy::correctReadcount(), HMMcopy::HMMsegment()

Other CNV: addHMMcopyCNV(), plotCNVheatmap(), removeCNV()

Examples

# Minimal synthetic example: build a table and pipe into runHMMCopy().
# Read depth doubles over 300 windows of chr2 (a copy-number gain).
set.seed(1)
n <- 2000L; w <- 5e4
cnv <- tibble::tibble(
    chr = factor(
        rep(c("chr1", "chr2"), each = n / 2),
        levels = c("chr1", "chr2")
    ),
    start = rep(seq(1, by = w, length.out = n / 2), 2),
    end = start + w - 1L,
    width = w, strand = "*",
    gc = runif(n, 0.05, 0.95), # wide spread => stable LOESS
    map = runif(n, 0.60, 0.98), # avoid extremes / exact 1.0
    sampleA = rpois(n, rep(c(500, 1000, 500), c(1500, 300, 200)))
) %>%
    runHMMCopy(colname = "sampleA", plotDir = NULL)
#> Warning: k-d tree limited by memory. ncmax= 438

# log2 copy-number estimate per window: about 1 over the gain, 0 elsewhere
table(round(cnv$sampleA))
#> 
#>    0    1 
#> 1700  300