Sequence-dependent histone variant positioning signatures

Genomic sequence partially distinguishes target sites of double variant nucleosomes with those of nucleosomes bearing only H3.3 or H2A.Z

High-resolution studies on crystal structure of the nucleosome core particle have revealed that the DNA is wrapped around the histone octamer in a flat, left-handed superhelix [1]. Based on this finding, the notion of sequence-dependent affinity of the histone core is defined as the energetic cost to bend the DNA to accommodate the superhelical path. The observation that nucleosomes show higher affinity for particular DNA sequences [33] has fostered many efforts to find genomic signatures related to genome-wide nucleosome positioning [34–37]. Although the role of genomic sequence on nucleosome distribution is still a controversial topic, it has been confirmed, both in vivo and in vitro, that DNA sequence imparts nucleosome distribution and there actually exists genomic code for nucleosome positioning [38–40]. Regarding important functions of H2A.Z, many efforts have been spent on identifying genomic signatures that may affect its genome-wide distribution. Tolstorukov et al. [41] compared genomic sequences of H2A.Z-containing nucleosomes from yeast and human and concluded that human nucleosomal sequences do not show the pattern of 10-bp periodicity as observed in yeast. Despite trying a handful of approaches, Gervais et al. [42] could not find any specific DNA motif that can help to distinguish H2A.Z-containing nucleosomal sequences from those of nucleosomes containing canonical histones. Other computational approaches also failed to distinguish H2A.Z-containing nucleosomal sequences from linker sequences [37]. Recently, the work of Jin et al. [32] gives the evidence that most of H2A.Z-containing nucleosomes observed at important genomic regions, such as promoters, CTCF binding sites, etc. are actually double variant nucleosomes. Taken together, we speculate that at important genomic regions (e.g. promoters) although histone variant-containing nucleosomes do not occupy thermodynamically favorite locations, DNA sequence may help to distinguish target sites of double variant nucleosomes from those of nucleosomes bearing only H2A.Z or H3.3.

To verify this hypothesis, we applied the computational procedure proposed by Peckham et al. [36] to the problem of discriminating double variant nucleosomal sequences from H2A.Z-only or H3.3-only nucleosomal sequences. This method has shown competitive performance for the task of discriminating “nucleosome forming” sequences from “nucleosome inhibiting” sequences and was applied sucessfully on human data [43]. The SVM classifiers in our work were built using Radial Basis Function (RBF) kernel (see Methods) instead of linear (dot product) kernel as used in original work. Sequence datasets, called double-h2az and double-h33 for training SVM classifiers were created as described in Methods. Each sequence in the training sets was represented as a 2, 772-entry vector, in which each entry is a normalized count of the occurences of a particular k-mer or its reverse complement, with k=1 up to 6. These vectors were used to train SVM classifiers.

To evaluate the performances of the resultant classifiers we used 10-fold cross-validation procedure. According to this procedure, each dataset was divided randomly into 10 subsets. The classifiers were trained on 9 subsets and tested on the remaining one. This training-testing procedure was repeated 10 times using a different hold-out set at each time. To measure the performances of the classifiers, we utilized the receive-operator-characteristic (ROC) curve. The quality of the classifier can be evaluated by calculating the area-under-the-curve (AUC) (the “ROC score”), in which a random classifier achieves the ROC score of 0.5 and a perfect classifier achieves the ROC score of 1.0. The average ROC scores from 10-fold cross validation on two sequence datasets, double-h2az and double-h33, were 0.62 (SD ≈ 0.02) and 0.63 (SD ≈ 0.03), respectively. These are significantly higher than the performance of the random classifier (p = 0.3 × 10^–5 and p = 0.15 × 10^–4 (t-test), correspondingly). This result shows that, genomic sequence contributes in part to targeting double variant nucleosomes to sites different from those of H2A.Z-only and H3.3-only nucleosomes.

Characteristics of sequences wrapping diffenrent histone variants

The result above suggests us to search for sequence features characterizing target sites of different histone variants. Basically, there are two kinds of such features [44]: one is compositional discriminative motif, which may help to distinguish double variant nucleosomal sequences from H2A.Z-only and H3.3-only ones; and the other is periodic pattern, which may appear in the set of nucleosomal sequences.

The periodicity has been known as one of the fundamental features that may appear in a set of DNA sequences. For example, 3bp sequence period is known to characterize coding sequences [45], while ~10bp sequence period is known to affect the curvature, bendability [46] and establish rotational setting on the histone surface [31] of the DNA sequences. To identify periodic patterns that may appear in nuclesomal sequences wrapping different histone variants, we employed autocorrelation analysis, which was sucessfully applied to detect hidden sequence periodicities in sets of DNA sequences [47, 48](see Methods). Investigating the periodograms of DNA motifs which show strongest discriminative power, we found clear periods of 2, 4, 6 of dinucleotides CA/TG in the set of 2169 double variant nucleosomal sequences (Figure 1); and clear periods of 2, 4, 6, 8 in the set of 348 H3.3-only nucleosomal sequences (Figure 2) while the set of 894 H2A.Z-only nucleosomal sequences does not exhibit this feature (Figure 3).

The enrichment (and possibly periodic patterns) of CA/TG dinucleotides in the genomes of many eukaryotic organisms can be resulted from CpG-methylation-deamination process, which dominates point substitutions in vertebrates [49], in which highly-methylated CpG dinucleotide is deaminated and then mutated to CA or TG [50, 51]. This may lead to the possibility that the observed periodicities of CA/TG are caused by random deposition of H3.3-containing nucleosomes to CA/TG-enriched regions generated by this process. However, a recent work [52] has reported that H3.3 variant preferentially targets CpG-rich promoters and a large number of H3.3 promoters has low DNA methylation level. Taken together with this, our result suggests that the periodicities of CA/TG dinucleotides may play important roles in determining the target sites of nucleosomes containing H3.3 (both H3.3-only and double) variant. Also, from Table 2 we can see that the enrichment of CA/TG dinucleotides is one of the strongest discriminative signals distinguishing target sites of double variant nucleosomes from those of H3.3-only nucleosomes. No DNA motif, however, shows dominant effect on this target site selection. This suggests that the role of DNA sequence may be modest compared with that of other factors (such as chromatin remodeling complexes, PTMs) in identifying exact locations for double variant nucleosomes.

Previous work has shown that in both D. melanogaster and human cells, H3.3 is deposited into transcribed genes, promoters and gene regulatory elements [6]. It is incorporated into genes upon transcription induction and is associated with transcriptional elongation [24]. Recent in vitro study about the effects of histone variants on transcription reported, however, that the presence of H3.3 is transparent to the transcription and the effect of “hybrid” particles containing both H2A.Z and H3.3 (i.e. double variant nucleosomes) was entirely dominated by the presence of H2A.Z [53]. These results, together with what was reported in Jin et al. [32], have raised important questions that: why double variant nucleosomes are present at important regulatory regions and whether or not H3.3 simply marks sites of nucleosome destabilization or is itself important for mediating the process [54]. Our result here provides a complementary answer to these questions: the presence of H3.3 in double variant nucleosomes can help to position the particles at specific locations in promoters; based on that, H2A.Z, which is also included in the particles, can execute its functions on transcription process. This is consistent with the hypothesis that the process incorporating H3.3 into nucleosomes, which involves histone chaperones and nucleosome remodelers, also facilitates H2A.Z incorporation [54]; and with the in vitro result that H3.3 may be just a marker of chromatin regions in flux while the presence of H2A.Z directly affects the transcriptional properties of the particles [53].

Data preparation

Support vector machine (SVM) classifiers

Given a training set containing instance-class pairs (x _i , y _i ), i = 1, 2, ⋯, l where x _i ∈ R ^l and y _i ∈ {–1, 1} is a class label, an SVM classifier is a hyperplane w ^T ϕ(x _i ) + b, where ϕ(x _i ) is a function mapping x _i into a higher (maybe infinite) dimensional space, that best separates the two classes. The hyperplane is obtained by solving the following primal optimization problem:

(1)

Its dual is a quadratic optimization problem:

(2)

where e is an unit vector, C > 0 is an error penalty parameter, Q _ij = y _i y _j K (x _i ,x _j ), K (x _i ,x _j ) = ϕ(x _i ) ^T ϕ(x _j ) is a kernel function. In our work, we employed Radial Basis Function (RBF) kernel K (x _i ,x _j ) = exp(–(x _i – x _i )²) to build SVM classifiers. Discriminant value for a testing instance x given by a trained classifier is: f(x) = ∑ _i α _i y _i K(x,x _i ) + b. Gist software package [59] was utilized for the tasks of training and testing SVM classifiers.

Feature selection with Fisher criterion

Feature selection is a process of selecting a subset of relevant features available from the data that most contribute to distinguishing instances from different classes. In our method, sequence features related to two kinds of histone variants, double variant and H3.3 only or double variant and H2A.Z only, were identified and ranked by their Fisher scores (or F-score in short). This is one of statistical criteria that is simple, effective and independent of the choice of classification method. The discriminative strength of each feature is defined as following:

Given a dataset X with two classes, denote instances in class 1 as X¹, and those in class 2 as X². Assume is the average of the jth feature in X ^k , the F-score of the jth feature is:

(3)

where

(4)

The numerator indicates the discrimination between two classes, and the denominator indicates the scatter within each class. The larger the F-score is, the more likely this feature is more discriminative. Gist software package [59] was used to calculate F-scores for different DNA motifs.

Autocorrelation analysis

Correlation functions measure the enrichment of certain pairs of motifs at a distance of k bp. To calculate XX-autocorrelation fuction we followed the method described in [48]. Given a motif X and a DNA sequence S, we count in the entire S the number of pairs of two identical motifs X and X separated by k base pairs. There are L – k – l + 1 pairs in a sequence of length L, where l is the length of the motif. Consequently, the probability to find the pair X – X at the distance k can be estimated as:

(5)

The probability to find a single motif X, denoted by P _X (k), can be estimated as:

(6)

where is the number of motif X in the sequence S. If the pairs at a distance k are statistically independent we have: P _XX (k) = P _X (k) * P _X (k). Thus the difference, C _XX (k) = P _XX (k) – P _X (k) * P _X (k), measures the correlation at a distance of k base pairs. A positive peak of the covariance C _XX (k) implies that there are more X – X pairs at a distance of k than expected by chance. The mean covariance function for a set of sequences was calculated by averaging individual functions over all the sequences in the set.