Extraction of functional binding sites from unique regulatory regions:
The Drosophila early developmental enhancers

Dmitri A. Papatsenko, Vsevolod J. Makeev, Alex P. Lifanov, Mirelle Regnier, Anna G. Nazina and Claude Desplan.

• SITE MOTIFS
• CUTOFF VALUES
• REFINED AND CONSISTENT MAPS
• FORMULATION OF SCANSEQ ALGORITHM
• TRAINING PROCEDURE
• PREDICTED MAPS
• TEST ON IMPERFECT SEQUENCE DATA
• TEST ON RHODOPSIN PROMOTERS
• SOFTWARE TOOLS
• REFERENCES

1 Developmental enhancers: maps of binding site distribution.

1.1 Site motifs.

• Bicoid: Figure a1, alignment; motif
• Caudal: Figure a2, alignment; motif
• FTZ: Figure a3, alignment; motif
• Giant: Figure a4, alignment; motif
• Hunchback: Figure a5, alignment; motif
• Knirps: Figure a6, alignment; motif
• Krüppel: Figure a7, alignment; motif
• Paired HD: Figure a8, alignment; motif
• Paired PD: Figure a9, alignment; motif
• Tailless: Figure a10, alignment; motif
• Tramtrack: Figure a11, alignment; motif

1.2 Cutoff values.

• Bicoid maps for optimal cutoff=5.3, Table a2
• Caudal maps for optimal cutoff=4.7, Table a3
• FTZ maps for optimal cutoff=4.2, Table a4
• Hunchback maps for optimal cutoff=6.2, Table a5
• Knirps maps for optimal cutoff=3.9, Table a6
• Krüppel:
• maps for 20 best matches, Table a7
• maps for peak1, cutoff = 4.6, Table a8
• maps for peak2, cutoff = 5.0 (optimal), Table a9
• maps for peak3, cutoff = 5.4, Table a10

1.3 Refined and Consistent maps.

• abdominal-A enhancer: Table a11, maps; sequence [1]
• buttonhead cis element: Table a12, maps; sequence [2, 3]
• engrailed intron: Table a13, maps; sequence [4, 5]
• even-skipped stripe2 enhancer: Table a14, maps; sequence [6-8]
• even-skipped stripe3+7 enhancer: Table a15, maps; sequence [9]
• even-skipped stripe4+6 enhancer: Table a16, virtual map only; sequence [10]
• fushi-tarazu proximal enhancer: Table a17, maps; sequence [11-13]
• fushi-tarazu zebra region: Table a18, initial map only; sequence [14-16]
• gooseberry enhancer: Table a19, initial map only; sequence [17-19]
• hairy stripe5 enhancer: Table a20, maps; sequence [20]
• hairy stripe6 enhancer: Table a21, maps; sequence [21]
• hairy stripe7 enhancer: Table a22, maps; sequence [22, 23]
• kruppel region 730: Table a23, maps; sequence [24]
• orthodenticle enhancer: Table a24, initial map only; sequence [14, 25, 26]
• runt stripe 5 enhancer: Table a25, maps; sequence [27]
• spalt early enhancer: Table a26, maps; sequence [28-30]
• tailless enhancer: Table a27, maps; sequence [31, 32]
• ultrabithorax BRE enhancer: Table a28, initial map only; sequence [33, 34]
• ultrabithorax PBX enhancer: Table a29, initial map only; sequence [35, 36]

2 Formulation of Scanseq algorithm.

Consensus: The most general mathematical definition of a "motif" is a set of allowed words. In the simplest case, a motif can be represented by "the consensus", for example in the form of standard IUPAC abbreviations. For instance, the motif containing words AGCTAT and AGCTGT can be described by the consensus AGCTRT. All words that can be derived from the consensus by a set of allowed operation (substitutions, insertions, deletions etc) define a consensus neighborhood. If only substitutions are allowed (the so-called Hamming distance) the consensus neighborhood can be found using counting with a mismatch scheme.

Positional weighted matrix: A more advanced method of motif representation is a positional weighted matrix (PWM). In this case, the chance of a substitution in a given position is weighted. Therefore, in the case of PWM, a threshold is required each time to define a motif. The popular way of defining PWM is a probabilistic model, in which one uses probabilities of observing a specific letter at a specific position of the motif. The likelihood of observing the word in this probabilistic model is compared to the chosen threshold to conclude if the word does or does not belong to the motif.. The most practical way of evaluating potential matches is to calculate the log-odds ratio of the resulting probability to the corresponding probability yielded by a random model [37], in this case called the "standard".

Hidden Markov Model: In the HMM representation [38, 39], the positions are no longer independent. Position probabilities are linked with the Markov model and the transitional matrices are estimated from the word set. In the case where there are no insertions or deletions allowed [40, 41], searching for matches with HMM can be performed with the help of log-odds ratios, essentially in the same way as with PWM. Currently we do not use an HMM representation in our approach.

2.2 Count of motif occurrence.

Counting in double stranded DNA To count motif occurrences in a sequence one needs to count all words in the sequence that belong to the motif. This can be done with several efficient algorithms (see [42] for the review). Not always however it is relevant to count the complementary DNA strands separately. This motivates a counting scheme with the motifs counted in both complementary strands. For example consider the following double stranded DNA sequence:

3'-ACAACACGT
5'-TGTTGTGCA

The number of occurrences of AC is four, three in the forward strand, and one in the reverse (strands are read in opposite directions). This number is equal to the number of occurrences in any of the strands of {AC, GT} motif, where GT is the word complementary to AC. Observe now, that the self-complementary word ACGT is found twice in the double-stranded sequence. In this case, the problem of counting of such self-overlapping motifs in double-stranded DNA can be reformulated as a problem of weighted counting in a single strand. To achieve this, the motif set is modified by addition of all complementary words. Any word H that has its complement in the initial set (particularly any palindromic word) is assigned weight w(H) = 2. Other words are assigned weight w(H) = 1.

Accounting for background distribution. It is insufficient to count simply the number of occurrences of some motif to infer whether it is "too abundant" in the sequence. If the input sequence exhibits strong composition biases the most frequent words found would be the words that contain the most frequent nucleotides. Therefore, the nucleotide frequencies (in the simplest case) must also be taken into account. A popular way to account for sequence composition [43] is to assess how often a motif found would appear in the random sequence with the same nucleotide content. In this case, Z-score (difference between the observed O and the expected E number of the motifs, normalized by the standard deviation V) is used to count the number of motif occurrences in the random sequence:
 

(1)

To calculate the Z-score, one needs to calculate the expected number of motif occurrences in the random sequence with a given nucleotide distribution, and estimate the variance for this quantity.

Counting with overlaps Calculation of the expectation and the variance of the number of motif occurrences in a random sequence is not a trivial problem, because one should take into account possible overlaps between words belonging to the same motif. Consider for instance the simplest case of a sequence with n=3 over a two-letter alphabet {AB} with equal probabilities of the letters. Inspection of the table below demonstrate that words AA and AB have the same expectation (0.5), but their variances are different:

Table a30. Overlap effect on expectation and variance.
 

	AA	AB
AAA	2	0
AAB	1	1
ABA	0	1
ABB	0	1
BAA	1	0
BAB	0	1
BBA	0	0
BBB	0	0
Expectation	4/8=0.5	4/8=0.5
Variance	6/8-0.25=0.5	4/8-0.25=0.25

Obviously this happens because the motif AA can overlap itself. This is also true in a general case [44]: whereas the expectation of the motif that contains many overlapping words does not change, the variance is changed, sometimes substantially.

Poisson approximation. If the probability of observation for two overlapping motifs in the sequence is much smaller than that for observation of one motif, the Poisson approximation  can be implemented. The applicability of the Poisson approximation depends on the nature of a motif.

Calculation of variance and expectation. It was found that both the variance and the expectation for the number of occurrences of a motif (a set of words) can be exactly calculated analytically in a random Bernoulli type sequence using the generation function technique [45, 46]. Bernoulli type sequence must contain no correlation between probabilities of letter occurrences at neighboring positions. Let H be the motif of m-words , . Let S be an n-letter sequence over the four-letter DNA alphabet with the vector of letter probabilities p. Then for each given word H, the probability will be, where i is the ith character of the word. For the entire motif one can calculate: .

The overlapping and correlation sets. It was observed that the variance depends on the structure of the words. Given two strings H and F, the overlapping is the H-suffixes that are F prefixes. Associated F suffixes form the correlation set A_H,F. Each word in A_H,F, when concatenated to the right of H creates an occurrence of F. For example, let H = AACAA, and F= AAAT, then A_H,F = {AT, AAT} and A_F,H is the empty set. The autocorrelation set A_H,H={,CAA, ACAA} also contains the empty string  and denotes the associated probability.

Variance and expectation for the single strand counting. It was recently demonstrated [45] that the expectation of the number of occurrences of a set of patterns H in a random sequence of length n generated by a Bernoulli model, N_H, is:
 

(2)

And the variance is:
 

(3)

Where

and  , but also can be analytically calculated.

Variance and expectation for the double strand counting In the case of counting at double strands, the word weights appear in formulae (2) and (3):

Let is the weighted number of occurrences of the weighted pattern H. Then its expectation writes:
 

(4)

and the linearity constant c₁ for the expectation writes:
 

(5)

An example Let H contain all words H that differ from some fixed word h for no more than k mismatches. If h is palindrome, then any word in its k-neighborhood has its complementary word in this k-neighborhood; hence all the words in the k-neighborhood are assigned weight 2, and  while . In the case of a true random sequence, Z-score does not change as compared to the case of single letter counting. However, the Poisson approximation V~E obviously does not satisfy.

3 Parameters and predictions.

3.1 Training procedure.

3.2 Predicted maps for 12 enhancers.

• abdominal-A enhancer, Figure a13 A, B
• engrailed intron, Figure a14 A, B
• even-skipped stripe2 enhancer, Figure a15 A, B
• even-skipped stripe3+7 enhancer, Figure a16 A, B
• even-skipped stripe4+6 enhancer, Figure a17 A, B
• fushi-tarazu proximal enhancer, Figure a18 A, B
• hairy stripe6 enhancer, Figure a19 A, B
• hairy stripe7 enhancer, Figure a20 A, B
• kruppel region 730, Figure a21 A, B
• runt stripe 5 enhancer, Figure a22 A, B
• spalt early enhancer, Figure a23 A, B
• tailless enhancer, Figure a24 A, B

4 Software tools.

• Programs to refine motif and to search for binding sites with PWM:
• MOTIF_CUT: finds the informational content for the user-defined alignment and generates PWM.
• MOTIF_MATCH: searches with the user defined PWM for matches in a given sequence.
• Programs to search for redundant motifs in a single sequence:
• SCANSEQ: finds related word families (motifs) for each word in a given sequence and calculates statistical significance (Z-score) of the found motifs.
• MAP_GENERATOR: builds a map of redundant motifs on the base of Scanseq data.
• Program to train parameters of the Scanseq:
• TRAIN_SCAN: finds best parameter combination (m_max, m_min, k_max, and c) for defined distribution of binding sites on a given sequence.
• Required accessory tools:
• GNUPLOT
• math library

5 References.