Systematic bias in high-throughput sequencing data and its correction by BEADS

Abstract

Genomic sequences obtained through high-throughput sequencing are not uniformly distributed across the genome. For example, sequencing data of total genomic DNA show significant, yet unexpected enrichments on promoters and exons. This systematic bias is a particular problem for techniques such as chromatin immunoprecipitation, where the signal for a target factor is plotted across genomic features. We have focused on data obtained from Illumina's Genome Analyser platform, where at least three factors contribute to sequence bias: GC content, mappability of sequencing reads, and regional biases that might be generated by local structure. We show that relying on input control as a normalizer is not generally appropriate due to sample to sample variation in bias. To correct sequence bias, we present BEADS (bias elimination algorithm for deep sequencing), a simple three-step normalization scheme that successfully unmasks real binding patterns in ChIP-seq data. We suggest that this procedure be done routinely prior to data interpretation and downstream analyses.

Figure 1.

DNA fragments in high-throughput sequencing data are not uniformly distributed over the genome. (a) The patterns of raw sequencing signals of independent C. elegans input sequence extracts and genomic DNA samples (black) are similar to underlying GC content (red) and mappability (blue). Positions 11 075 000–11 098 000 of chromosome I of the C. elegans genome are shown. (b) GC frequency distributions of the C. elegans genome (solid line) and a set of input sequence reads (dashed line). (c) GC frequency ratio between input sequence data and the C. elegans genome.

Figure 2.

Bias in GC content (a, d, g, j), mappability (b, e, h, k) and raw input sequence signals (c, f, i, l) across internal exons and around transcript start sites in C. elegans and human. The error bars represent the 95% confidence intervals of the estimated mean GC, mappability or sequence signal values.

Figure 3.

BEADS normalization of high-throughput sequence reads of C. elegans input sequence (a, b), H3K4me3 ChIP (c, d), H3K36me3 ChIP (e, f) and human input sequence (g, h) libraries. Shown are the results following each step of correction: Raw (uncorrected), GC corrected, GC + mappability corrected, and GC + mappability + local corrected signals. Plots show signals across internal exons and around transcript start sites of all genes; except in d and f, only transcript start sites of highly expressed genes were used (see ‘Materials and Methods’ section). The normalized signals are plotted as relative normalized read counts (left-hand y-axis). The fully normalized signal is also plotted as fold-change relative to the genomic average (right-hand y-axis). Solid lines show average signal and shaded regions show 95% confidence intervals. Genomic read-count averages for the C. elegans input-control, H3K4me3, H3K36me3 and human input-control libraries are 33.0, 29.8, 14.6 and 0.4, respectively.