damidseq_pipeline: an automated pipeline for processing DamID sequencing datasets

1DamID is a well-established technique for discovering regions of DNA bound by or associated with proteins. It has been used to map the genome-wide binding of transcription factors, chromatin proteins, nuclear complexes associated with DNA and RNA pol II (for e.g.;;; ). The technique can be performed in cell culture, whole organisms or with cell-type specificity, and requires no fixation or antibody purification.

2DamID involves the fusion of a bacterial DNA adenine methylase (Dam) to any DNA-associated protein of interest. The bacterial Dam protein methylates adenine in the sequence GATC and, given that higher eukaryotes lack native adenine methylation, the DNA-binding footprint of the protein of interest is uniquely detectable through isolating sequences flanked by methylated GATC sites. However, a major consideration with DamID is that any Dam protein within the nucleus will non-specifically methylate adenines in GATC sequences at accessible regions of the genome. For this reason, DamID is always performed concurrently with a Dam-only control, and the final DNA-binding profile is typically presented as a log2(Dam-fusion/Dam-only) ratio.

3Although the majority of published DamID experiments have used tiling microarrays for data analysis, next-generation sequencing (NGS) allows greater sensitivity and higher accuracy. Although several recent studies have used NGS with DamID, these have relied upon a comparison of peak binding intensities between read-count-normalized Dam-fusion and Dam samples. Depending on the characteristics of the Dam-fusion protein (see later) this approach may lead to real signal being lost, and correct normalization of the datasets is required to detect all binding by many Dam-fusion proteins. Here, we describe a software pipeline for the automated processing of DamID-sequencing (DamID-seq) data, including normalization and background reduction algorithms.

4Although DamID-seq data can be aligned and binned as per all NGS data, two issues arise that are specific to DamID. The first major consideration is the correct normalization of the Dam-fusion and Dam-control samples. The greatest contribution to many Dam-fusion protein datasets is the non-specific methylation of accessible genomic regions (e.g. Fig. 1B), with a mean correlation between Dam alone and Dam-fusion datasets of 0.70 (n = 4, Spearman’s correlation). Representing the data as a (Dam-fusion/Dam) ratio in theory negates such non-specific methylation. However, strong methylation signals at highly bound regions in the Dam-fusion dataset will reduce the relative numbers of reads present at accessible genomic regions in this dataset (see, for example, the occupancy of Dam-RNA Pol II over the eyeless gene in Fig. 1), and normalizing the data based on read counts alone can therefore produce a strong negative bias to the ratio file [Fig. 1B (iii), Supplementary Fig. S5A]. Depending on the characteristics of the fusion protein, this negative bias can lead to real signal being lost (Fig. 1). Although microarray data inadvertently overcame this bias through the manual adjustment of laser intensities during microarray scanning, until now no method has existed for correctly normalizing DamID-seq datasets.

5In order to correct for this negative bias we use the read counts from accessible genomic regions—as determined from the Dam-only dataset—as the basis for normalization, while avoiding regions likely to contain real signal in the Dam-fusion sample. We use the following algorithm to adjust the Dam-fusion dataset. Given the GATC-site resolution of DamID, we divide the read counts into GATC fragments.All GATC fragments lacking read counts are excluded. The remaining GATC fragments are divided into deciles.Given the high probability that the highest 10% of Dam-fusion read counts represent bound signal rather than background signal, we exclude fragments that have scores in this decile.The first three deciles of the Dam sample can generate inconsistent normalization values if included (Supplementary Table S2), so we exclude fragments that lie within this range.The distribution of the log2(Dam-fusion/Dam) ratio for all remaining fragments is determined via the Gaussian kernel density estimate, where h is the bandwidth, estimated via the method of: (where σ is the standard deviation of the sample and IQR the interquartile range). For speed considerations, we estimate kernel density over 300 equally spaced points within the interval.The point of maximum kernel density represents the point of maximum correspondence between Dam-fusion and Dam values; if both samples are correctly normalized this value should equal 0. We therefore normalize all Dam-fusion values by.In addition to ensuring correct normalization, a second important consideration is the reduction of background noise. Regions without specific methylation will have randomly distributed background counts that, when a ratio file is generated, will generate a large degree of noise. Such noise can potentially obscure peak detection. In order to mitigate this effect we add pseudocounts to both datasets. In order to maintain equivalence between replicates with differing numbers of reads (assuming that genomebound ≪ genomeunbound) the number of pseudocounts added is proportional to the sequencing coverage, thus, where c is a constant. (Supplementary Table S1 for a comparison of gene calls with different read-depths). Adding pseudocounts increases the number and the total genomic coverage of detected peaks and increases the signal:noise ratio (Supplementary Figs S1–S4).

6Given the GATC-site resolution of DamID, we divide the read counts into GATC fragments.

7All GATC fragments lacking read counts are excluded. The remaining GATC fragments are divided into deciles.

8Given the high probability that the highest 10% of Dam-fusion read counts represent bound signal rather than background signal, we exclude fragments that have scores in this decile.

9The first three deciles of the Dam sample can generate inconsistent normalization values if included (Supplementary Table S2), so we exclude fragments that lie within this range.

10The distribution of the log2(Dam-fusion/Dam) ratio for all remaining fragments is determined via the Gaussian kernel density estimate, where h is the bandwidth, estimated via the method of: (where σ is the standard deviation of the sample and IQR the interquartile range). For speed considerations, we estimate kernel density over 300 equally spaced points within the interval.

11The point of maximum kernel density represents the point of maximum correspondence between Dam-fusion and Dam values; if both samples are correctly normalized this value should equal 0. We therefore normalize all Dam-fusion values by.

12The combination of these two methods compares favorably with previously published microarray data [Fig. 1B (iv)] or DamID-seq data (Supplementary Figs S1–S4; Supplementary Fig. S5).

13The damidseq_pipeline software is implemented in Perl, and will process multiple single-end read sequencing files in FASTQ or BAM format. The pipeline can match sequencing adaptors to sample names, automatically identifies the Dam-only control, and performs alignment, read-length extension, normalization, background reduction and ratio file generation. (Supplementary Methods for details).

14A large number of user-configurable options are provided, including the ability to adjust the normalization algorithm parameters, generate read-count normalized files and add a user-specified number of pseudocounts. Parameters specified on the command-line can be saved as defaults if the user desires.

15The damidseq_pipeline software is open-source and freely available at http://owenjm.github.io/damidseq_pipeline. A detailed set of installation and usage instructions are provided at the above website, along with a small example dataset.