Copy number profiles from DNA methylation arrays

If you run Illumina methylation arrays on tumors, you already have a copy number profile sitting in your IDAT files. The array was designed to measure methylation, but every probe also reports how much DNA hybridized to it, and that amount scales with copy number. This is why methylation-based brain tumor classification reports usually come with a genome-wide copy number plot next to the predicted class (Capper et al., 2018). One assay, two very different kinds of evidence, and the copy number plot is often what makes a pathologist comfortable with a call.

This post walks through how I generate those profiles in R with minfi and conumee, and how I read them. All the code runs on the example data in minfiData, so you can follow along before pointing it at your own IDATs.

Why total intensity tracks copy number

Each CpG on the array is measured by a methylated and an unmethylated signal. For Infinium I probes these come from two different bead types; for Infinium II probes they come from the two color channels of a single bead. The beta value, roughly M / (M + U), deliberately throws away the total amount of signal so that you get a proportion between 0 and 1.

For copy number you want exactly the part that beta discards. The sum M + U at a probe is (approximately) proportional to the number of DNA copies at that locus that made it onto the array. A region with an extra copy produces more total signal, and a deleted region produces less.

The catch is that raw totals are dominated by things that have nothing to do with copy number:

probe affinity, which varies a lot from probe to probe
GC content and probe design type (I vs II)
sample-level differences in DNA input, bisulfite conversion and hybridization

So you cannot read copy number off a single sample’s intensities. You compare the sample to a set of normal controls processed the same way, and most of the probe-specific effects cancel out in the ratio.

conumee does this in a few steps:

Fit. For each query sample, CNV.fit runs a multiple linear regression of the query’s combined intensities on the control samples’ intensities. The fitted values are the linear combination of controls that best resembles the query. It then computes a log2 ratio of query over fitted value for every probe.
Bin. CNV.bin groups neighboring probes into bins (by default at least 15 probes and at least 50 kb per bin) and takes the median log2 ratio per bin. It also computes a shift that moves the copy-neutral state close to zero, by minimizing the median absolute deviation of all bins from zero.
Detail. CNV.detail summarizes predefined regions of interest, typically oncogenes and tumor suppressors.
Segment. CNV.segment runs circular binary segmentation from the DNAcopy package on the binned values.

The one assumption that matters most is that the controls have a flat genome. Anything that is not flat in the controls will show up, inverted, in every query.

Setup

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install(c("minfi", "minfiData", "conumee", "CopyNeutralIMA"))

conumee pulls in the 450k and EPIC manifest and annotation packages it needs. Note that conumee (version 1.x) works in hg19 coordinates, so any downstream comparison with hg38 data needs a liftover.

Reading IDAT files with minfi

minfiData ships six 450k samples as raw IDATs together with a sample sheet, which makes it a convenient stand-in for a real experiment.

library(minfi)
library(conumee)

base_dir <- system.file("extdata", package = "minfiData")
targets <- read.metharray.sheet(base_dir)
rgset <- read.metharray.exp(targets = targets)
rgset

For your own data the pattern is the same. If you have an Illumina-style sample sheet, read.metharray.sheet finds it and builds the Basename column for you. If you just have a folder of IDATs, you can skip the sheet:

# Schematic: needs your own IDATs or a large GEO download, so it isn't run here
# Placeholder paths: replace with your own
my_dir <- "/path/to/your/idats"

targets <- read.metharray.sheet(my_dir)              # if you have a sample sheet
rgset_mine <- read.metharray.exp(targets = targets)

rgset_mine <- read.metharray.exp(base = my_dir, recursive = TRUE)  # without a sheet

Keep 450k and EPIC samples in separate objects. I process each array type on its own and deal with the overlap at the annotation step (more on that below).

Before anything else, check sample quality. Detection p-values are the quickest signal:

det_p <- detectionP(rgset)
colMeans(det_p > 0.01)   # fraction of probes failing detection, per sample

A sample with a large fraction of failed probes will usually give a noisy copy number profile too, and it is better to know that before you start interpreting bumps.

Preprocessing choices

conumee needs a MethylSet, because it uses the methylated and unmethylated intensities separately before summing them. That rules out preprocessing methods that only return ratios: preprocessQuantile returns a GenomicRatioSet, for example, and so does preprocessFunnorm with its default settings.

I use preprocessIllumina, which is also what the conumee vignette uses. It does background correction and normalizes to the array’s internal control probes, mimicking GenomeStudio. preprocessRaw and preprocessNoob also return a MethylSet and will work. The more important rule is that query and control samples go through exactly the same preprocessing. Mixing raw controls with normalized queries is a reliable way to make an interesting-looking but meaningless plot.

From here on I use RGsetEx, the same six samples already loaded as an RGChannelSet, because its phenotype data has a status column marking normal and cancer samples.

data(RGsetEx, package = "minfiData")
mset <- preprocessIllumina(RGsetEx)

pData(mset)[, c("Sample_Name", "Sample_Group", "status")]

Building the annotation object

The annotation object defines which probes are used, how they are grouped into bins, and which regions get a detailed look. You only need to build it once per array type and parameter set, so I save it with saveRDS.

data(exclude_regions)
data(detail_regions)

anno <- CNV.create_anno(
  array_type      = "450k",
  exclude_regions = exclude_regions,
  detail_regions  = detail_regions
)
anno

A few arguments worth knowing:

array_type is one of "450k", "EPIC" or "overlap". Use "overlap" when the query and controls come from different array types; it restricts the analysis to probes present on both.
bin_minprobes (default 15), bin_minsize (default 50000) and bin_maxsize (default 5e6) control binning. The defaults were tuned for 450k data.
chrXY defaults to FALSE, so the sex chromosomes are excluded. If you turn it on, your controls need to be sex-matched to the query or you will “discover” X chromosome losses in every male sample.
exclude_regions removes known polymorphic regions that would otherwise produce recurring artifacts.
detail_regions lists regions to summarize individually. You can pass a GRanges object or a path to a BED file with a name column.

The bundled detail_regions covers a set of commonly altered cancer genes, including EGFR, CDK4, MDM2, MYC, MYCN, PTEN, RB1, TP53 and a combined CDKN2A/B region. Run detail_regions to see the full list. For brain tumor work it is worth adding a few regions of your own; PDGFRA, for example, is not in the default list.

Loading intensities and choosing controls

CNV.load sums the methylated and unmethylated signals per probe. For a MethylSet it takes sample names from the first phenotype column whose name contains “name” (here Sample_Name), so samples are called GroupA_1, GroupB_1 and so on.

d <- CNV.load(mset)
names(d)

It also checks the intensities and warns if the average is unusually low or high. Take those warnings seriously; they often mean a preprocessing mismatch or a failed sample.

Option 1: normals from the same experiment

minfiData includes three normal samples, which is enough to see how the pieces fit:

is_normal <- pData(mset)$status == "normal"
ctrl_inhouse <- d[is_normal]

Three controls is a demo, not a reference set. The regression has very little to work with, and any quirk in one control propagates into every profile.

Option 2: a public set of copy-neutral samples

The CopyNeutralIMA package provides samples from healthy individuals with no expected copy number changes (51 on the 450k array and 13 on EPIC, collected from public GEO series) through ExperimentHub, meant for exactly this use. It downloads the data the first time you run it.

library(CopyNeutralIMA)

ima <- annotation(mset)[["array"]]     # "IlluminaHumanMethylation450k"
rg_ctrl <- getCopyNeutralRGSet(ima)
mset_ctrl <- preprocessIllumina(rg_ctrl) # same preprocessing as the queries
ctrl <- CNV.load(mset_ctrl)

Older tutorials, including the conumee vignette, load controls from the CopyNumber450kData package. That package is not available in recent Bioconductor releases, which is why I point people to CopyNeutralIMA instead.

What makes a good reference set

In rough order of importance:

Same array type, or use array_type = "overlap".
Same preprocessing as the queries.
Copy-neutral. Normal tissue or blood from individuals without known constitutional copy number changes. Do not use “normal-looking” tumor samples.
Similar technical history. Controls run in the same lab, with similar DNA extraction and the same sample type (FFPE vs frozen), absorb more of the technical bias. The conumee authors recommend matched normal tissue profiled in the same experiment when you can get it.
Enough of them. More controls give the regression more flexibility to match each query. I would not go below a dozen or so if I had a choice.

Tissue matching helps but is less critical than you might expect, because total intensity is much less sensitive to methylation state than beta is. Normal brain is a natural reference for brain tumors. Flow-sorted blood datasets (for example FlowSorted.Blood.450k and FlowSorted.Blood.EPIC on Bioconductor) are also copy-neutral and are a common choice when no tissue normals are available.

Running conumee

With the annotation, query and controls in hand, the analysis is four calls. GroupB_1 is one of the cancer samples in minfiData.

x <- CNV.fit(query = d["GroupB_1"], ref = ctrl, anno = anno)
x <- CNV.bin(x)
x <- CNV.detail(x)
x <- CNV.segment(x)
x

If a query sample is also part of the reference set, CNV.fit notices (correlation above 0.99) and leaves it out of the fit with a message. That is a guard rail, not an invitation; keep your references separate.

To process every tumor in a cohort:

tumors <- names(d)[pData(mset)$status == "cancer"]

fits <- lapply(tumors, function(s) {
  CNV.segment(CNV.detail(CNV.bin(CNV.fit(d[s], ctrl, anno))))
})
names(fits) <- tumors

CNV.segment wraps DNAcopy, and its defaults (alpha = 0.001, nperm = 50000, min.width = 5, undo.splits = "sdundo", undo.SD = 2.2) were chosen for 450k data. Lowering undo.SD keeps more breakpoints; raising it merges more segments. I leave them alone unless I have a specific reason, and if I change them I change them for the whole cohort.

Plots

# Placeholder output path
pdf("GroupB_1_cnv.pdf", width = 12, height = 5)
CNV.genomeplot(x)                      # whole genome
CNV.genomeplot(x, chr = "chr7")        # single chromosome
CNV.detailplot(x, name = "EGFR")       # one detail region
CNV.detailplot(x, name = "CDKN2A/B")
CNV.detailplot_wrap(x)                 # all detail regions in one figure
dev.off()

In the genome plot, each dot is a bin, colored by its value (with the default colors, losses toward red and gains toward green), and segments are drawn as lines. Dashed vertical lines mark centromeres, and detail regions are labeled along the top. In the detail plot, individual probes are shown along with the bins that cover the region.

CNV.genomeplot accepts ylim (default c(-1.25, 1.25)), which is worth adjusting when a high-level amplification runs off the top of the plot.

Reading the plot

Most interpretation mistakes come from reading the y-axis as if it were an absolute copy number. It is not.

Zero is the most common state, not necessarily two copies. The shift computed in CNV.bin places the bulk of the genome at zero. In a near-diploid tumor that is two copies. In a near-triploid tumor with many changes, “zero” may be three copies, and a true diploid region can look like a loss.

Amplitudes are compressed. In theory a single-copy gain in a pure diploid tumor is log2(3/2), about 0.58, and a single-copy loss is log2(1/2) = -1. In practice you see much smaller shifts, because of contaminating normal cells, background signal and the array’s nonlinear response to DNA amount. Arm-level changes often appear as modest but clearly consistent offsets across hundreds of bins.

With that in mind, these are the patterns I look for:

Whole-chromosome and whole-arm changes. A whole arm sitting above or below the baseline, with breakpoints at or near the centromere.
Focal high-level amplifications. A small number of bins far above everything else. EGFR on 7p is the classic example in glioblastoma; MYCN, CDK4, MDM2 and PDGFRA are others you will meet in brain tumors.
Homozygous deletions. A narrow, deep drop, often at CDKN2A/B on 9p21, frequently sitting inside a broader single-copy loss of 9p. The signal never goes to zero, because normal cells in the sample still contribute two copies. Small homozygous deletions may be covered by only a handful of probes, so look at the detail plot and not just the genome plot.
Chromosome 7 gain with chromosome 10 loss. Combined whole-chromosome +7/-10 is a characteristic pattern of IDH-wildtype glioblastoma, often with EGFR amplification and CDKN2A/B deletion on top.
1p/19q codeletion. Loss of the entire 1p arm and the entire 19q arm is the defining alteration of IDH-mutant oligodendroglioma, classically caused by an unbalanced translocation. Whole-arm is the operative word. Partial 1p or 19q losses occur in other tumors, including IDH-wildtype glioblastoma, and are not the same thing.

A few cautions. Array-based copy number gives you relative gains and losses at bin resolution. It does not detect copy-neutral loss of heterozygosity, balanced rearrangements or small events between probes, and it struggles with subclonal changes present in a minority of cells. It is excellent supporting evidence and a good screening tool. For a result that drives a diagnosis, use whatever assay your lab has validated for that purpose (FISH, NGS-based copy number, or a validated array pipeline), and interpret copy number together with the methylation class and histology, not in isolation.

Purity and noise

Two things make profiles hard to read: low tumor content and noisy data. They look different.

Low purity compresses everything toward zero. A tumor with real arm-level changes can look nearly flat if it is mostly normal tissue. A flat profile therefore means “no detectable change,” not “normal genome.” Methylation-based purity estimates are a useful sanity check here.

Noise shows up as scatter between neighboring bins and as wavy patterns that segmentation may turn into many short, low-amplitude segments. Common causes are degraded DNA from FFPE, low input, poor bisulfite conversion and a reference set that does not match the query well. conumee stores a simple noise estimate in the fit, the root mean squared difference between log2 ratios of neighboring probes:

x@fit$noise

sapply(fits, function(f) f@fit$noise)   # compare across the cohort

I use it to rank samples within a cohort rather than as an absolute cutoff, because the scale depends on the array type, preprocessing and reference set. Samples at the noisy end get looked at carefully before anyone interprets a small segment.

Two more checks I find useful:

Recurring artifacts. If the same small segment appears in most samples at the same position, suspect the reference set or a polymorphic region, not biology.
Known truth. If you have a few samples with copy number from another assay (or cell lines with well-known karyotypes), run them through the same pipeline. It is the fastest way to calibrate what a real single-copy change looks like in your hands.

450k vs EPIC

The EPIC array has roughly 850,000 probes. It keeps most of the 450k content, drops some 450k probes and adds many probes in enhancer regions. For copy number, the practical consequences are:

Use array_type = "EPIC" for EPIC queries with EPIC controls.
Use array_type = "overlap" when mixing array types, for example EPIC queries against 450k controls. You lose the EPIC-only probes but keep comparability.
Bins depend on probe density, so the same bin_minprobes and bin_minsize settings give a different bin layout on EPIC than on 450k. Do not compare bin-level values across arrays directly; compare segments or arm-level calls.
CopyNeutralIMA has EPIC controls as well: annotation(mset)[["array"]] returns "IlluminaHumanMethylationEPIC" for an EPIC MethylSet, and getCopyNeutralRGSet accepts that value. The EPIC set is smaller, so in-house EPIC normals are worth collecting if you run that array routinely.

Exporting results

CNV.write returns a data frame, or writes a file if you give it a path.

seg <- CNV.write(x, what = "segments")
head(seg)

# Placeholder output paths
CNV.write(x, what = "segments", file = "GroupB_1.CNVsegments.seg")
CNV.write(x, what = "bins",     file = "GroupB_1.CNVbins.igv")
CNV.write(x, what = "detail",   file = "GroupB_1.CNVdetail.txt")
CNV.write(x, what = "probes",   file = "GroupB_1.CNVprobes.igv")

The segment table has one row per segment with columns ID, chrom, loc.start, loc.end, num.mark, bstat, pval, seg.mean and seg.median, with the shift from CNV.bin already applied. It is in SEG format, so you can drop it straight into IGV, and it is a convenient starting point for cohort-level summaries such as arm-level gain and loss frequencies. The detail table gives one value per detail region, which is the easiest way to tabulate EGFR or CDKN2A/B status across a cohort. Any threshold you put on these values should be calibrated against samples with known status, because the right cutoff depends on purity, noise and your reference set.

A checklist for your own data

Read IDATs with read.metharray.exp, one array type per object.
Check detection p-values and drop or flag failed samples.
Preprocess queries and controls identically with a method that returns a MethylSet (preprocessIllumina is a good default).
Build the annotation once with the right array_type, exclude_regions and your own detail_regions, and save it.
Use a copy-neutral reference set of the same array type, as technically similar to your samples as you can manage.
Run CNV.fit, CNV.bin, CNV.detail and CNV.segment with the same settings for the whole cohort.
Look at every genome plot yourself, with the noise estimate and an idea of purity next to it.
Export segments and detail values, and calibrate any calls against samples with known status.

Note added later: the Hovestadt lab has since released conumee 2.0 (conumee2, on GitHub as hovestadtlab/conumee2), which adds support for the EPIC v2 array and mouse arrays and a separate step for detecting focal amplifications and homozygous deletions. The overall logic described here (normalize against copy-neutral controls, bin, segment) carries over, but check its documentation for the current function arguments.

References

Aryee MJ, et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics 30, 1363-1369 (2014).
Capper D, et al. DNA methylation-based classification of central nervous system tumours. Nature 555, 469-474 (2018).
Hovestadt V, Zapatka M. conumee: Enhanced copy-number variation analysis using Illumina DNA methylation arrays. Bioconductor package.
Olshen AB, et al. Circular binary segmentation for the analysis of array-based DNA copy number data. Biostatistics 5, 557-572 (2004).
Pastor Hostench X, Przybilla M. CopyNeutralIMA. Bioconductor package.
Sturm D, et al. Hotspot mutations in H3F3A and IDH1 define distinct epigenetic and biological subgroups of glioblastoma. Cancer Cell 22, 425-437 (2012).

← All writing