F-Seq: a feature density estimator for high-throughput sequence tags

1 INTRODUCTION

High-throughput sequencing technologies that generate short sequence reads can be used to identify specific genomic features, such as transcription factor binding sites (Johnson et al., 2007; Robertson et al., 2007) and regions of open chromatin (Boyle et al., 2008) at a genome-wide level. In general, locations of biologically relevant features are defined by the presence of an enrichment of mapped sequence reads. To date, there is no standard means to summarize and visually display these data in an intuitive way. As the use of high-throughput sequencing becomes more prevalent, there is a growing need for a method to efficiently identify statistically significant genomic features based on sequence tags.

Published research using high-throughput sequencing data have employed histograms to calculate regions of dense sequence reads and make calls on sites of interest (Johnson et al., 2007; Robertson et al., 2007). Histograms are a non-parametric density estimator where the region covered is divided into equal-sized bins whose height is represented by the count of hits within that bin. These methods can be problematic as histograms are not smooth and can be strongly affected by the start/end points of the bins and the width of the bins (Fig. 1A and B).

To counteract bin boundary effects, one can instead calculate a kernel density estimate centered at each sequence allowing these estimates to overlap (Fig. 1D) (Parzen, 1962). Using a smooth kernel such as a Gaussian generates a smooth signal. This method does not alleviate the problem of bin width (or in the case of kernel density estimation, bandwidth) (Fig. 1C). Determination of an optimal bandwidth can present a problem, but this can be overcome by using the argument that minimizes the asymptotic mean integrated squared error (or other minimization techniques). However, the sparsity of data and size of whole-genomic sequences does not allow for estimating bandwidth with this method. Therefore, we suggest the use of a bandwidth based on the size of the feature being identified.

Although histogram methods have provided usable results, the dependency of resolution on bin size and the lack of statistical rigor in the treatment of the data begs for a new approach. We have developed an algorithm that uses kernel density estimation that can provide both a discrete and continuous probability landscape to better display genomic features of interest across the genome. These kernel density estimation-based probabilities, calculated at each base, are directly proportional to the probability of seeing a sequence read at that location.

2 F-SEQ DENSITY ESTIMATION

To generate the kernel density estimation, we consider the problem where we are given n sample points along a chromosome of length L. Our goal is to locate regions with high sample density. If we assume the points {x_i}_i=1ⁿ are sampled as , then an estimate of this probability density function (pdf) will provide a significance measure for high density regions. We use the univariate kernel density estimation (kde) to infer the pdf, written as

(1)

where b is a bandwidth parameter controlling the smoothness of the density estimates, and K() is a Gaussian kernel function with mean 0 and variance 1. Instead of explicitly setting b, a user provides a feature length parameter (default=600) which controls the sharpness of the pdf estimate. Larger features will naturally lead to smoother density estimates.

Computing the density at each point in the chromosome using all n points is computationally expensive and exceeds the precision available to common computing platforms. We therefore compute a default window size w as a function of the bandwidth parameter b and the Gaussian kernel such that

We expect that window sizes for typical bandwidth settings will be on the order of a few thousand, significantly less than the many millions of bases available.

We also compute a threshold level for evaluating the significance of density regions using the following background model:

1. Compute an average number of features for window w as n_w=nw/L.

2. Calculate the kernel density at a fixed point, x_c, within the window given a random and uniform distribution of the n_w features.

3. Repeat step 2 k times to obtain a distribution of the kernel density estimates for x_c. For large k the kdes become normally distributed.

4. The threshold is s SDs above the mean of this normal distribution.

Larger values of s reduce the false discovery rate and provide a natural statistical interpretation to the veracity of these density regions.

F-Seq takes an input a BED format file (http://genome.ucsc.edu/FAQ/FAQformat#format1) containing aligned sequence tags. Since calculation of kernel density estimation requires a point measure for each sequence, we use the estimated center of the DNA fragment being sequenced. In many cases, such as from ChIP-seq protocols, the aligned sequence represents only the 5′ end of a longer fragment and therefore should be extended to the average fragment size in the experiment. In the case of DNase-Seq protocols where the 5′ end of the sequence represents the point of enrichment, the alignment should be shortened to 1 bp in length. A perl script has been included to perform this task.

Output files can be created either as a continuous probability wiggle format (http://genome.ucsc.edu/goldenPath/help/wiggle.html) or as a discrete-scored regions BED format. The discrete regions are those where the continuous probability is above the threshold s SDs above the background mean. These output files are ready for immediate import into the UCSC Genome Browser (Kent et al., 2002) (http://genome.ucsc.edu).

3 EXAMPLE APPLICATIONS