The Wild West of spike-in normalization

Abstract

Spike-in normalization is a powerful approach to assess global changes in data obtained from genomic mapping of DNA-associated proteins by methods such as ChIP-sequencing (ChIP-seq)^, or CUT&RUN. While multiple spike-in methods provide detailed documentation, the implementation of these approaches often omit critical quality control steps and veer from the established procedures. Spike-in normalization typically makes use of a single scalar to normalize genome-wide data, making the approach particularly vulnerable to errors in implementation. Here, we show that proper application of spike-in normalization can increase quantification accuracy across a spectrum of conditions and outline how misuse of spike-in approaches can create erroneous biological interpretations. We conclude by providing guidelines to minimize pitfalls when applying this approach to normalize data from protein-DNA interaction results.

Fig. 1.. Demonstration of the ability of spike-in normalization to accurately capture signal variation over wide and narrow dynamic ranges.

(a) Reanalysis of data from Orlando et al., where untreated cells (high H3K79me2 signal) were mixed with DOT1 inhibited Jurkat cells (low H3K79me2 signal) in 5 ratios of treated/untreated (0/100, 25/75, 50/50, 75/25, 100/0). Dynamic range of H3K79me2 signal was measured by using Western blot results from Orlando et al. with ImageJ; the change in H3K79me2 between 0% treated and 100% treated samples was approximately 10-fold (Supplemental Fig. 1a). Maximum signal (10X) and minimum signal (1X) are labeled on x axes. The H3K79me2 ChIP-seq signal was quantified and plotted against the line of expected signal either by using standard read depth normalization (left) or spike-in normalization (right). Accuracy of the fit of the expected line was determined by R squared (the value is reported in top right of each plot). (b) The results of a similar titration experiment generated to focus on a narrower dynamic range. Mitotic-arrested cells were generated by treating with thymidine then S-trityl-cysteine (STC), hereafter termed “mitotic”. We estimate approximately 85-90% of cells arrested in prometaphase by this method. From previous mass-spectrometry data, the approximate fold change between mitotic H3K9ac (low H3K9ac signal) and interphase H3K9ac (high H3K9ac signal) was 3x, labeled on x axes (Supplemental Fig. 1b). We used interphase cells (high H3K9ac) mixed with mitotic-arrested cells (low H3K9ac) in six ratios (100/0, 95/5, 75/25, 50/50, 25/75, 0/100). To each of the samples, we spiked-in both D. melanogaster and S. cerevisiae chromatin. Quantification of the H3K9ac signal after read-depth normalization (left) or spike-in normalization using Drosophila (right) was plotted as in (a). The normalization using S. cerevisiae provided similar results (Supplemental Fig. 3). Within each plot, H3K79me3 or H3K9ac signal was min-max normalized according to the following equation: ${\mathsf{z}}_{\mathsf{i}}=({\mathsf{x}}_{\mathsf{i}}\text{-}\mathsf{min}(\mathsf{x}))∕(\max (\mathsf{x})\text{-}\mathsf{min}(\mathsf{x}))$. Here, $\mathsf{min}(\mathsf{x})$ is the average minimum signal (100% treated cells in (a) or 100% mitotic cells in (b)), $\max (\mathsf{x})$ is the average maximum signal (0% treated cells in (a) or 0% mitotic cells in (b)). ${\mathsf{X}}_{\mathsf{i}}$ is the signal for each sample and ${\mathsf{z}}_{\mathsf{i}}$ is the minmax normalized signal for each sample. QCs for the dataset are in Supplemental Fig. 2-6 and Supplemental Table 3.

Figure 2.. A schematic depicting the impact of misuse of spike-in normalization on downstream results.

Left, “True Enrichment” shows a hypothetical scenario of two conditions with variable histone acetylation profiles. Panels a-e show various scenarios observed in literature and their effect on both local signal (depicted as schematics of individual ChIP-seq peaks) and global signal (shown as a scatterplot of log-normalized signal for all peaks for each condition). The data represents 3 replicates, while for simplicity, the peaks and scatter plots show only the average of the replicates. Of note – these scenarios are not mutually exclusive, and combinations of these misuse of spike-in normalization can occur. (a) Traditional ChIP-seq with no spike-in added – the signal appears identical between the two conditions, not capturing the true patterns. (b) Properly performed normalization with spike-in ChIP (blue) yields results that accurately follow the original ground truth. (c) Variable proportion of spike-in added to the samples, this skews the analysis and condition #2 appears to wrongly have higher signal. (d) Low yield of spike-in data, either from input chromatin levels or technical issues. This precludes any QC of the spike-in. Also, the low number of reads that align to the spike-in genome could highly skew the results one way or the other (depicted as overlaying normalized peaks and two plots in both directions). (e) Exogenous chromatin from a species which is phylogenetically close to the sample species. This could skew the results in both directions, depending on the percent of misassigned reads in each sample (depicted as in (d)).

Figure 3.. Variations in the ratios of spike-in to sample chromatin in public datasets.

We examined 53 datasets. 51 of these were obtained from our survey. We additionally included the data we generated for Fig. 1b and Supplemental Fig. 2 (n=2). Of these only 27 had sufficient input samples for each condition (52.9%, pie chart; insert). One additional dataset plotted had > 1 input but was still missing inputs for some conditions. We aligned the inputs to a concatenated spike-in and target genomes (Methods). For each dataset, the input with lowest ratio of spike-in/sample reads was scaled down to 1, to capture the variation between input samples in the same dataset. Y-axis is shown on log2 scale to capture the diversity in variation within datasets. Alignment information for each individual sample is available in Supplemental Table 1, GEO accessions for each dataset are in Supplemental Table 5.

Figure 4.. An example of the ability of proper normalization strategy to correct for variations in the spike-in/sample ratio.

Mitotically arrested HeLa-S3 cells were treated with either the pan-HDAC inhibitor TSA (to increase global acetylation) or DMSO as a control. For each treatment, using the same chromatin preparation the amount of added spike-in was varied over five orders of magnitude (0.00025 to 2.5x spike-in/target ratio). The data were normalized by three approaches: (a) read-depth normalization – the change in the signal between the conditions cannot be observed; (b) normalization of relative IP ratios of spike-in and target – most of the change in the signal is not observable and the amount of spike-in added impacts the general signal; and (c) account for the ratio of spike-in/target by using the non-IP input controls and then normalize as in (b) using the relative spike-in/target IP – the known differences between the conditions can be detected in all cases but the one where the amount of spike-in is higher than the target. QCs for the dataset are in Supplemental Fig. 8-14 and Supplemental Table 4.