The SEQC2 epigenomics quality control (EpiQC) study

Abstract

Background: Cytosine modifications in DNA such as 5-methylcytosine (5mC) underlie a broad range of developmental processes, maintain cellular lineage specification, and can define or stratify types of cancer and other diseases. However, the wide variety of approaches available to interrogate these modifications has created a need for harmonized materials, methods, and rigorous benchmarking to improve genome-wide methylome sequencing applications in clinical and basic research. Here, we present a multi-platform assessment and cross-validated resource for epigenetics research from the FDA's Epigenomics Quality Control Group.

Results: Each sample is processed in multiple replicates by three whole-genome bisulfite sequencing (WGBS) protocols (TruSeq DNA methylation, Accel-NGS MethylSeq, and SPLAT), oxidative bisulfite sequencing (TrueMethyl), enzymatic deamination method (EMSeq), targeted methylation sequencing (Illumina Methyl Capture EPIC), single-molecule long-read nanopore sequencing from Oxford Nanopore Technologies, and 850k Illumina methylation arrays. After rigorous quality assessment and comparison to Illumina EPIC methylation microarrays and testing on a range of algorithms (Bismark, BitmapperBS, bwa-meth, and BitMapperBS), we find overall high concordance between assays, but also differences in efficiency of read mapping, CpG capture, coverage, and platform performance, and variable performance across 26 microarray normalization algorithms.

Conclusions: The data provided herein can guide the use of these DNA reference materials in epigenomics research, as well as provide best practices for experimental design in future studies. By leveraging seven human cell lines that are designated as publicly available reference materials, these data can be used as a baseline to advance epigenomics research.

Fig. 1

Sequencing and alignment metrics of whole methylome libraries, including all replicates across all cell lines. EM = EMSeq; MS = MethylSeq; SP=SPLAT; TS = TruSeq; TM = TrueMethyl. a Distribution of reference-based read alignment outcomes, including primary mapped reads (both mates mapped in correct orientation within a certain distance), multi-mapped reads (read pairs containing secondary or supplementary alignments), reads marked as PCR or optical duplicates, and unmapped reads. Ambiguous and duplicate reads can be a subset of properly aligned reads. b Median insert size distributions derived from distance between aligned paired end reads. c Percentage of bases trimmed per replicate, either due to low base quality, adapter content, or dovetailing reads. d Cumulative genomic coverage plot, averaged across cell line per assay. Coverage is cut off at 200× in this plot, but extends beyond for all assays. Dotted line indicates 20× mean coverage. e Nucleotide bias plot showing the log2 enrichment of covered versus expected mono- and di-nucleotides. f The relationship between the number of read pairs sequenced per assay and the mean depth of coverage per CpG dinucleotide, showing sequencing depth required to achieve a certain level of coverage. 20× CpG coverage is shown as the dotted line. g Same as f, but plotted using total bases sequenced, to include Oxford Nanopore sequencing, which produces variable read lengths

Fig. 2

Coverage of CpGs across the genome. All samples visualized here were downsampled to 20× mean coverage per CpG. a Empirical cumulative distribution functions for median coverage, averaged across samples for HG002-HG007. b Standard deviation between replicate beta values for HG002 as a function of average coverage. The expected curve (computed based on the assumption that replicate beta values are independent and identically distributed estimates of a common proportion p) is added as a solid black curve. c Intersection of CpG coverage (min 5×) across Chromosome 1. Exact values of CpGs covered per assay are shown on the right. d Count and genomic annotation for CpGs uniquely covered by an assay (left) and uniquely not covered by an assay (right). Up5kb = 5 kb upstream distance from promoter region; Promoter = within 1 kb upstream of transcript start site. e Distribution of coverage in CpG shelves, shores, and islands. EM = EMSeq; MS = MethylSeq; SP=SPLAT; TS = TruSeq; TM = TrueMethyl. f Mean coverage curves around transcript start sites (TSS)

Fig. 3

Estimates of methylation per CpG across the genome for HG002. All samples visualized here were downsampled to 20× mean coverage per CpG. a Methylation percentage distributions per assay. b Methylation bias (mbias) plots showing mean methylation per base for short-read assays (Nanopore excluded here). Dotted lines indicate recommended cutoffs for methylation calling for these data. Original top/bottom refer to mappings to bisulfite-converted strands in the reference genome. c Metagene plot showing mean methylation across genomic feature per assay. Promoter regions span 1 kb upstream of transcript start sites (TSS). d Mean methylation curves surrounding TSS across all genes. e Pearson correlation matrix of genome-wide methylation estimates. f Pearson correlation matrix of methylation estimates for sites where methylation was estimated to be between 20 and 80%. g Methylation percentage correlation between Oxford Nanopore and all other assays. Pearson correlation values shown on top. Marginal histograms show methylation curves per assay

Fig. 4

Mosaic plots illustrating agreement between assays for differentially methylated per assay (DMA) sites as coverage levels vary. Rows represent the number of the six assays for which each DMA site was also identified, with values ranging from 1 (indicating no other assays, shaded in red) to 6 (indicating all assays, shaded in purple). Columns indicate the median coverage across HG002-HG007, with values ranging between the 5th and 95th percentiles for each assay

Fig. 5

Microarray normalization and low-varying site definition. a Densities showing the percentage of DNA methylation variation explained by cell line across the epigenome (N = 677,520 overlapping CpG sites) for each normalization method. b Raw beta values at each of the 59 SNP probes on the Illumina EPIC arrays, with samples colored by lab. c Variance in methylation beta values (no normalization) within each genotype cluster at the 59 SNP probes, separated and colored by lab. The dotted vertical line represents the 95th percentile. d Variance in methylation beta values (normalized with funnorm + RCP) across the epigenome. Sites in the shaded area, which have less variation than 95% of SNP probe genotype clusters, are defined as low-varying sites. e Percentage of methylation (normalized with funnorm + RCP) variance explained by cell line across the epigenome, stratified by high-varying vs. low-varying sites

Fig. 6

a Density plots of sequencing/microarray concordance indicating the percent of variance explained (VE) by cell line, assay (sequencing or microarray), and residual variation for 841,833 CpG sites with complete information in all assays. b Distribution of percent variance explained by cell line in the sequencing/microarray variance partition analysis as a function of beta value variance (binwidth = 0.001) and median coverage (binwidth = 1) at each CpG site. 90% of the y-axis values fall between the outermost dotted lines for each bin along the x-axis