Aligning sequences with repetitive motifs?

Abstract

Pairwise sequence alignment is among the most intensively studied problems in computational biology. F We present a method for alignment of two sequences con- L taining repetitive motifs. This is motivated by biological studies of proteins with zinc nger domain, an important C H group of regulatory proteins. Due to their evolutionary his- Zn toofryd,i seeqreunetnczeisncof tnhgeeserspr(oshteoirnts scuobnsteaqiuneancveasriwaibtlhe nspuemcibecr H C H T symbols at each position). T G E K P YF G E K P Our algorithm uses two types of hidden Markov models (HMM): pair HMMs and pro le HMMs. Pro le HMMs Fig. 1. The structure of a zinc- nger. Highly variable sites describe the structure of sequence motifs. Pair HMMs as- are marked with black color. The most conserved amino sign a probability to alignment of two motifs. Combination acids are the four involved in binding the zinc ion [14]. of the these two types of models yields an algorithm that uses di erent score when aligning conserved vs. variable motif residues. The dynamic programming algorithm that computes the motif alignments is based on the well known positions are very conserved due to their importance Viterbi algorithm. We evaluated our model on sequences of in assuming desired function, while other positions are zinc nger proteins and compared it with existing alterna- highly variable, since they distinguish speci c DNA setives. quences where individual zinc ngers bind (Figure 1). We will focus our attention on the KRAB-ZNF 1 Introduction proteins that have a region encoding one or more Kruppel-associated box domains (KRAB, [2]) followed Pairwise sequence alignment is one of the most stud- by a zinc nger region (Fig. 2). The human genome enied problems in bioinformatics. We will concentrate on codes more that 600 of proteins from this family, and alignment of protein sequences, where a protein can a lot of e ort is dedicated to building and maintaining be represented as a string over the alphabet of 20 dif- their catalogues [9], [11], [3]. Complicated repetitive ferent amino acids. During the evolution, particular structure of these genes is a result of a dynamic evoamino acids in a protein can be substituted by an- lutionary history, full of sequence duplications [7, 12], other amino acid, or even get inserted or deleted. The and many mutations which help to gain new functions goal of sequence alignment is to compare two proteins, for duplicated copies. quantify their sequence similarity, and to identify pairs The repetitive nature of zinc nger protein of amino acids that have likely evolved from the same sequences complicates their sequence alignment. Traamino acid in the common ancestor. Over the years, ditional alignment methods based purely on sequence multitude of variations of this problem have been in- similarity frequently misalign individual zinc ngers, troduced and many practical software tools were de- or even align a single zinc nger in one sequence to veloped. parts of several di erent zinc ngers in the other seOur work is motivated by the study of zinc nger quence. Consequently, many studies of these proteins proteins. These proteins contain a variable number of limit their analyses and infer conclusions based only up to 40 zinc nger domains [18]. Zinc nger domain is on the the KRAB domains or sequences before the zinc a stretch of approximately 28 amino acids, the purpose nger region (e.g. [14], [7]), or dispute the relevance of of which is to bind DNA at speci c places. Comparison standard methods applied to genes with high variance of zinc ngers form di erent proteins reveals that some in the number of ngers [16].

-

? This work was supported by the European Community FP7 Marie Curie grants IRG-224885 to TV and IRG-231025 to BB, and by a grant from VEGA 1/1085/12.

In this work, we develop a new method for aligning sequences with repetitive motifs, such as zinc nger proteins. To overcome the problems outlined above, we combine the strength of pro le hidden Markov models

KRAB A

KRAB A KRAB B which are used to characterize the properties of these repetitive motifs, and pair hidden Markov models as a model of sequence alignments. We compare our work to MotifAligner that was previously used to align zinc

nger proteins [ 11 ], and we nd our new method to produce more accurate alignments on our testing set.

In the rest of this section, we introduce necessary background and notation and describe the MotifAligner approach to repetitive sequence alignment in more detail. In Section 2, we describe our new pro le-pro le-pair alignment method (PPP). We present the results of experimental comparison of PPP and MotifAligner in Section 3. 1.1

Background and notation align them by inserting dashes to individual sequences so that they all have the same length and when we arrange them in a table, as in Figure 3, many columns contain the same or similar amino acids. Several consecutive dashes form a gap in the alignment, indicating that a part of the sequence was deleted or inserted during the evolution. The sequence alignment problem can be formulated as an optimization problem and solved by existing algorithms. For two sequences, the problem can be solved easily by Needleman-Wunsch dynamic programming algorithm [ 10 ], for multiple sequences it is NP-hard [ 6 ]. The scoring function for pairwise alignment is typically based on a substitution matrix scoring all pairs of aligned amino acids and on parameters for scoring gaps: gap opening penalty g for the rst dash in a gap and gap extension penalty e for each additional gap.

ZNF626_4799/12 YKC--EECGKAF-NQSSILTTHERIILERNZNF727_4861/2 YKC--EECGKDC--RLSDFTIQKRIHTADRS ZXDB_644/5 YQCAFSGCKKTF-ITVSALFSHNRAHFREQE LLNL1236_4814/2 SMC--PECSKTSATDSSCLLMHQRSHTGKRP ZNF23_141/15 FQC--KECGKAF-HVNAHLIRHQRSHTGEKP In this paper, we rely on several standard tools from Fig. 3. Alignment of ve sequences of zinc nger motifs computational biology, namely alignments, pair hid- from human proteins. den Markov models, and pro le hidden Markov models, which we brie y explain in this section.

We start by de ning hidden Markov models One way of systematically deriving a scoring func(HMMs). An HMM is a probabilistic nite state au- tion for pairwise alignments is to use pair HMMs [ 4 ]. tomaton. We can use it to generate a random sequence These models emit two sequences simultaneously. In over some alphabet as follows. We start in a desig- one step, the HMM can emit a single character in one nated start state B. In each step, we sample a charac- of the sequences or in both. The later case corresponds ter of the sequence from the emission probability dis- to two symbols aligned to each other, the former to tribution associated with the current state and then a symbol aligned to a dash. Figure 4 shows the pair randomly change the state according to the transition HMM used in our work. The match state M emits probability distribution. The process ends when we pairs of aligned characters, state X emits characters reach the designated nal state E. only in the rst sequence, and state Y emits charac

The sequence of states visited in the individual ters only in the second sequence. Given two sequences, steps is called a state path. We will denote the proba- we can nd the most probable state path that could bility of emitting x in state v as ev(x) and the proba- generate them and this will give us an alignment of bility of transition from state v to w as tv;w. The joint these two sequences. probability of emitting a sequence x = x1 : : : xn along To represent a typical sequence of a motif, we will the state path s = s1 : : : sn in a given HMM is use another kind of HMMs, called pro le HMMs [ 4 ].

A pro le HMM is typically constructed based on an n alignment of several motif instances, such as the one P (x; s) = es1 (x1) Y tsi 1si esi (xi): in Figure 3. Each position of the motif is represented i=2 by one state with emission probabilities set to the obA typical task solved with HMMs is to nd the most served frequencies of amino acids in the corresponding probable state path that could generate a given se- alignment column (possibly with some pseudocounts quence, i.e. to nd s = arg maxs P (x; s). This task added to avoid zero probabilities). These so called is solved by the Viterbi algorithm based on dynamic match states are arranged in a chain (see Figure 5). programming [ 19 ]. These states used alone would generate sequences of

The second important notion is sequence align- the same length. However, real sequences may have ment. Given a set of related protein sequences, we can various insertions and deletions compared to the conB

1-2δ-τ

HMMER in the original input sequences x and y, re

spectively. In the second step, MotifAligner computes scores of all gapless pairwise alignments of motifs tk; u`, for all 1 k a, 1 ` b: s[tk; u`] =

L X S[tki ; u`i ]; i=1 (1) where S[xi; yj ] is the score of aligning amino acids xi and yj (they use a standard BLOSUM85 substitution matrix [ 8 ]; motif occurrences are padded to have the same length).

In this way we obtain a similarity score between each pair of motif occurrences. Next MotifAligner applies the Needleman-Wunsch algorithm [ 10 ] to T and U , treating motifs as sequence symbols and using matrix s as the substitution matrix. In this way we obtain pairs of aligned zinc ngers between the two proteins. 2

Pro le-pro le-pair alignment 1.2

MotifAligner approach

In this section, we present a new approach to align

ment of sequences with repetitive motifs. We adopt an approach similar to MotifAligner, however, we change Fig. 5. Example of a pro le HMM. States Mk are match the alignment algorithm and the scoring scheme to states, Ik are insert states and Dk are delete states. States take into account the structure of the repeated motif. B, E, and D1 : : : D3 are silent, which means that they do For example, the zinc- nger motif (Fig. 2) contains not generate any characters. several highly conserved positions, among them the four amino-acids binding the zinc ion (positions 3, 7, 20, 24). These four amino acids are crucial to the funcsensus motif; these are modeled by additional insert tion of the motif and as such should be used to anchor and delete states. Given a pro le HMM and a se- the whole alignment. However, the fact that these poquence, we can again nd the most probable state sitions match in the two aligned sequences should not path, which in this case gives us an alignment of the be very surprising and should not by itself contribute sequence to the motif represented by the pro le HMM. much to the resulting score. On the other hand, there Note, however, that the pro le HMM emits only a sin- are several variable positions, and the di erences at gle sequence; the motif itself is represented directly in these positions will be very informative of the evoluthe structure and parameters of the model. tionary distance.

To take these issues into account, we have developed a new pro le-pro le-pair method (PPP) for

To obtain high quality alignments even on sequences

with highly variable number of zinc nger motifs, Nowick et al. developed a pairwise alignment tool called MotifAligner [ 11 ]. To our knowledge, it is the only sequence alignment method designed speci cally to align sequences with variable number of repetitive motifs.

Part of our work was inspired by this algorithm.

MotifAligner rst uses a pro le HMM tool HMMER [ 5 ] and nds all canonical motif occurrences with statistically signi cant scores in both input sequences. Let T = (t1; : : : ; ta) and U = (u1; : : : ; ub) be the sequences of all motif occurrences found by Fig. 6. The pro le HMM of random 2000 human zinc ngers from the complete dataset, viewed as a HMM logo [ 15 ].

Thus our goal is to compute three paths sp; sx; sy through the pair HMM and the two pro le HMMs that would satisfy our constraints and the product of joint probabilities implied by all three models would be maximized: (sp; sx; sy) = arg max score(x; y; sp; sx; sy); spv;aslxid;sy (2) alignment of individual motifs. The method uses a combination of two pro le HMMs and a pair HMM for sequence alignment and aligns the two sequences by nding the best possible path through all three models simultaneously. To align the complete protein sequences containing these repeating motifs, we rst align each possible pair of motifs through PPP, compute their similarity score, and use a modi cation of a traditional global alignment algorithm, now operating on individual motif occurrences as a unit. We describe the details of the method in the remainder of this section. where score(x; y; sp; sx; sy) = Ppair(x; y; sp) Ppro le(x; sx) Ppro le(y; sy), where Ppair(x; y; s) is the joint probability of the state path s aligning sequences x and y in the pair HMM and Ppro le(x; s) is 2.1 Pairwise alignment of individual motifs the joint probability of the state path s and sequence x in the pro le HMM.

The input to PPP consists of two instances of the re- We obtain an optimal solution using dynamic propeating motif x = x1 : : : xLx and y = y1 : : : yLy , a pro- gramming similar to the Viterbi algorithm used to le HMM encoding the same motif, and a pair HMM compute the most probable state paths in individual characterizing the properties of a typical alignment. HMMs. Let S = (Sp; Sx; Sy) be the triplet of states Our goal is to align both x and y to a separate copy of of pair and x-pro le and y-pro le models satisfying pro le HMM and at the same time, use the pair HMM our conditions. We denote V [Sp; Sx; Sy; i; j] the score as a glue. of the highest scoring state path combination ending

In particular, we are simultaneously seeking the with the triplet S and covering the pre xes x1 : : : xi, three paths through the three HMMs that satisfy the y1 : : : yj of the two sequences. following constraints: The computation of V [Sp; Sx; Sy; i; j] depends on the types of states Sp; Sx; Sy. For example, if Sp is the Constraint 1 (Pro le match states constraint) match state M of the pair HMM and Sx is the match If xi and yj are emitted by the same match state Mk state Mk of the pro le HMM, then according to our in their pro le models then the pair model has emit xi constraints Sy must be the same match state Mk and and yj together in the match state M . we have the following recurrence: V [M; Mk; Mk; i; j] = eM (xi; yj )eMk (xi)eMk (yj )

Constraint 2 (Pair match state constraint) If

the pair model emits xi and yj together in the match state M then both pro le models emit xi and yj in the same match state Mk or in the same insert state Ik.

In other words, if the pair model is in the state X or Y (which is interpreted as a gap in one of the sequences), the two pro le models should not be in the same match state: symbols belonging to the same consensus column should be aligned. However, if both pro le models are in the same insert state they can max either be evolutionarily related, in which case they issopfthreaonttthuhdeelseedynoptfbcrltyooeh,ueawlldelphigahminciarheovdmdeweobuoldseuseielnadln.rgeCcinoMonsrnererseistttsterhpaadeotirnneintdino2ttfhoatetlhuhsseoseeiqisnmpaugamepiXnrleicemesam,nionatddhtdecaYehlt-, >>>>>>>>>>>>>>>>>>>>tttYYXMMMtttffMMooIkrr``MM001 MkkttkIMtk``In;k<nM1M1k<kMkkkVV;V[nY[Y[;X6=;MM;`I`;k`;MI1kn;I;k1i; i1;1i;1j;1j; j1]11]] satraetien Mcokmnpolerteinly tdhie esraemntecionlusemrtnsstoarteinItkhe(is.ea.meeitchoelr- >>>>>>>>tY M tfIokr 01Mk t`M<`Mkk V [Y; Ik 1; M`; i 1; j 1] umn k, but di erent states Mk and Ik), which means :tY M tIk 1Mk tIk 1Mk V [Y; Ik 1;Ik 1; i 1; j 1] that the symbols being emitted are unrelated, then The value at V [M; Mk; Mk; i; j] has to include the the pair model should not be in the match state. These emission probabilities of xi and yj in all three models. constraints thus ensure that the sequence and the pro- Then we take a maximum over all choices of previous le alignment can be interpreted in a consistent man- cells from which the current cell value could be comner. puted. Every value considered in the maximum is the 8>tMM tM`Mk tM`Mk V [M; M`; M`; i 1; j 1] >>>> for 0 ` < k > >>>tMM tIk 1Mk tIk 1Mk V [M; Ik 1; Ik 1;i 1; j 1] >>>>tXM tM`Mk tMnMk V [X; M`; Mn; i 1; j 1] >>>> for 0 `; n < k; n 6= ` > >>>tXM tM`Mk tIk 1Mk V [X; M`; Ik 1; i 1; j 1] >>>> for 0 ` < k > >><tXM tIk 1Mk tM`Mk V [X; Ik 1; M`; i 1; j 1]

for 0 ` < k product of the value of the predecessor cell and tran- Number of Finger Motifs sition probabilities in all three models. All the other Genome Genes Variants Total Average Median cases can be derived analogously; we omit the deriva- hg19 612 1071 13363 12.48 12 tions due to the space constraints. mm9 302 513 5226 10.19 10

Every time we compute a value for any cell, we keep rchaneMFaamc22 457778 1802180 129124539 1121..0128 1121 a pointer to the cell from which the value was derived.

We use those pointers later to trace back the resulting Table 1. The complete dataset, based on genes from the state paths. Of particular importance is the path in the whole human genome. One gene can have multiple variants pair model, since it de nes the alignment of x and y. that di er in organization of zinc ngers. For each of the resulting state path, we also compute its joint probability its respective model, obtaining values Ppair(x; y; sp), Ppro le(x; sx); and Ppro le(y; sy). 2.2

Alignment of complete motif arrays

and the length of motifs is almost the same, so we can say that Lx; Ly; L = O(m) and hence the time required to compute the alignment of one motif pair is O(m6). From the same observation, one can easily see that the space complexity is O(n2 + m4). The running time and memory is practical, since values of n and m tend to be small in real proteins (for zinc- nger arrays, both n and m are less than 30).

We use the same procedure as MotifAligner for alignment of complete motif arrays. We compute all pairwise alignments of individual motifs, where the score of a pairwise motif alignment is based on joint probabilities of motif sequences and state paths in all three models as described below. Since we perform the motif alignment for all pairs of motifs and assign a score to each such alignment, we get a scoring system simi- 3 Experiments and evaluation lar to a scoring matrix. Treating motifs as symbols and using this scoring matrix, we obtain the full alignment Gold standard data set. We evaluated our approach of input motif arrays using Needleman-Wunsch algo- on human zinc- nger genes and their counterparts in rithm. related species macaque, mouse, and dog. We down

More formally, for motif arrays Ax = (x1; : : : ; xn) loaded the set of annotations of KRAB zinc nger and Ay = (y1; : : : ym), we calculate n m matrix S, genes from the Human KZNF Catalog [ 9 ] and where remapped the annotation to the current human genome assembly hg19 using liftOver tool. To S(xi; yj ) = ln Ppair(xi; yj ; sp;i;j ) ; (3) obtain the sequences of these zinc nger genes in other

Ppro le(xi; sx;i;j )Ppro le(yj ; sy;i;j ) species, we used the whole genome alignments from where sp;i;j , sx;i;j and sy;i;j are the three state paths the UCSC genome browser [ 17 ] as a mapping between computed when aligning motifs xi, yj by PPP. This the human (hg19) and the macaque (rheMac2), score compares the hypothesis that the two motif se- mouse (mm9), and dog (canFam2) genomes. quences are related (given by probability from the pair The resulting genomic sequences were translated HMM) to the hypothesis that these are simply two into amino acid sequences and cleaned for apparent independent sequences following the same pro le (as artifacts. In particular, we removed genes that condetermined by scores from the two pro le HMMs). tained ngers shorter than 10 amino acids. Summary statistics of the resulting dataset is shown in Table 1.

Because of relatively high time complexity of the PPP 2.3 Algorithm complexity algorithm, alignment of genes with high number of ngers takes a lot of time. For that reason, we prepared a subset of the complete dataset, omitting genes from human chromosome 19 and their putative orthologs in other genomes. These genes contain the highest numbers of repeating motifs (30 or more). Summary statistics for this restricted dataset are shown in the Table 2.

The time complexity of the PPP algorithm on two motif arrays with O(n) motifs, each of length O(m) is O(n2m6). There are O(n2) individual motif alignments. The time needed to compute one such alignment is O(LxLyL4), where Lx; Ly are the lengths of motifs and L is the number of columns in the pro

le HMM. This follows from the observation that in the recurrent step of individual motif alignment we

ll 3 L L Lx Ly matrix, and time required to compute each cell is at most O(L2), the upper bound on the number of values considered in the recurrence.

Typically, the number of columns in the pro le HMM

Model parameters. The emission probabilities of the

pair HMM used in our experiments were based on the BLOSUM85 substitution matrix. This particular matrix was chosen in order to compare our results to MotifAligner [ 11 ]. In particular, we used the probability

Genome hg19 mm9 canFam2 rheMac2 distributions p and q from which the matrix was derived, as supplied in EMBOSS software package [ 13 ].

The transition probability parameters (see Figure 4) were set as follows: = 0:0345, so that the expected length of an alignment is 28, which is the length of a typical human C2H2 zinc nger motif; = 0:05185 so that the expected length of a match region is 13:45, because the most variable region of a zinc nger motif spans positions 12-15; " = 0:4769 so that the expected length of a gap is 1:1.

The complete parameter set of the pro le model was acquired from the Pfam database entry for the ZNF C2H2 family [ 1 ]. The length of the pro le is 23, which is shorter than a typical human zinc nger motif. The reason is that the model is based on a more diverse set of sequences from various species.

Histogram of ppp.related[ , 2 ] − (ppp.related[ , 3 ] + ppp.related[ , 4 ]) Histogram of ppp.unrelated[ , 2 ] − (ppp.unrelated[ , 3 ] + ppp.unrelated[ , 4 ]) where the related set was treated as positives and the random set as negative examples. The classi cation 3.1 The PPP score distribution performance of the PPP is clearly better, demonstrating that our scheme is more suitable as a score for To compare the scoring function of the PPP model classi cation of paired motifs from random pairs. with the scores used by MotifAligner, we created two sets of zinc- nger motif pairs. The related set contained 1000 ngers from the human genome, each 3.2 Alignment accuracy paired with the corresponding nger from macaque, Next we use the PPP and the simpler MotifAligner mouse or dog. The random set contained 1000 ran- method for scoring pairs of zinc ngers as building dom pairs of ngers; we assume that these ngers are blocks in the whole motif array alignment. The on average more distantly related to each other than Needleman-Wunsch algorithm for the whole motif arpaired ngers in the rst set. ray alignment has three parameters: the gap opening

We have computed a PPP alignment of each se- penalty g, the gap extension penalty e, and the subquence pair in both samples. The score distributions stitution matrix s that scores individual motif alignare shown in Figure 7. Both distributions resemble the ments. For MotifAligner, we have used the original normal distribution, with mean of the related set close parameters [ 11 ], in particular the BLOSUM85 substito 20 and mean of the random set at around 5. tution matrix and the gap penalties set to g = 84 and

For comparison purposes, we have reimplemented e = 75:6. In the PPP model, the matrix s is deterMotifAligner algorithm as described in [ 11 ]. Figure 7 mined by the equation 3, and we have tested several shows the score distributions of the MotifAligner ap- di erent settings of the parameters g and e. proach to alignment of individual motifs, based on the We carried out three tests. In the rst one, we BLOSUM85 substitution matrix. These score distribu- aligned all zinc nger arrays of orthologous proteins in tions do not resemble the normal distribution. In par- the complete dataset. The second and the third experticular, the distribution for the related set has a heavy iments simulated a loss of ngers during the evolution tail, which is clearly not desirable. { we created two arti cial datasets with 1/5 and 1/3

The important property of the scoring scheme is of the total number of ngers removed in each zinc how well it is able to distinguish positive examples nger array in all four genomes, and we aligned the from negative ones. Figure 9 shows a ROC curve, original human dataset with the four reduced sets. the only existing program speci cally designed to align sequences with repetitive motifs.

There is still a room for improvement of our work.

Apart from obvious upgrades, like a more e cient implementation, the underlying model can be enhanced in several ways. For example, an interesting question is whether some other scoring function of individual motif alignments would perform better. Such a function might be based on di erent properties of the underlying models, e.g. the full probability of a sequence, instead of the probability of the Viterbi path.

To alleviate problems caused by the computational complexity of the algorithm, various heuristics could be applied, especially methods avoiding exhausting computations of the whole dynamic programming matrix. In order to apply our model to other protein families with repeating motifs, a more robust procedure for parameter estimation should be established. In addition, a method for assessment of statistical signi cance of alignments may be helpful when computing alignments of large datasets where random similarities are more likely to occur.

The model we have implemented is not the only way of doing sequence alignment with repetitive motifs. It is very appealing to use a monolithic probabilistic model instead of multiplying probabilities of three separate models. We have tried to develop such a model, but we were not able to overcome some of e t a R e v iits o P e u r T

PPP

MotifAligner

False Positive Rate We have designed and implemented an algorithm for alignment of sequences with repetitive motifs. The algorithm is built on top of two types of hidden Markov models. It utilizes positional information from two copies of a pro le HMM and uses a pair HMM to align the motif sequences. We were able to apply our model on real world data, and obtained better results than

References

1. A. Bateman , S. Boehm , E. L. L. Sonnhammer , F. Gago: Mulitple sequence alignment of zinc nger C2H2 type family . Pfam Family: zf-C2H2 (PF00096) , 2011 . Online. http://pfam.sanger.ac.uk/family/PF00096.

2. E. J. Bellefroid , D. A. Poncelet , P. J. Lecocq , O. Revelant , J. A. Martial : The evolutionarily conserved Krppel-associated box domain de nes a subfamily of eukaryotic multi ngered proteins . Proc. Natl. Acad. Sci . U.S.A. 88 ( 9 ), 1991 , 3608 { 3612 .

3. G. Ding , P. Lorenz , M. Kreutzer , Y. Li , H. J. Thiesen: SysZNF: the C2H2 zinc nger gene database . Nucleic Acids Res . 37 ( Database issue ), 2009 , D267{ 273 .

4. R. Durbin , S. Eddy , A. Krogh , G. Mitchison: Biological sequence analysis . 1st edition . Cambridge University Press. 1998 , 356 p., ISBN : 978 - 0521629713 .

5. S. R. Eddy : Accelerated pro le HMM searches . PLoS Comput. Biol . 7 ( 10 ), 2011 , e1002195 .

6. I. Elias: Settling the intractability of multiple alignment . Journal of Computational Biology 13 ( 7 ), 2006 , 1323 { 1339 .

7. A. T. Hamilton , S. Huntley , M. Tran-Gyam , D. M. Baggott , L. Gordon , L. Stubbs: Evolutionary expansion and divergence in the ZNF91 subfamily of primate-speci c zinc nger genes . Genome Res . 16 ( 5 ), 2006 , 584 { 594 .

8. S. Heniko , J. G. Heniko : Amino acid substitution matrices from protein blocks . Proc. Natl. Acad. Sci . U.S.A., 89 ( 22 ), 1992 , 10915 { 10919 .

9. S. Huntley , D. M. Baggott , A. T. Hamilton , M. TranGyam , S. Yang , J. Kim , L. Gordon , E. Branscomb , L. Stubbs: A comprehensive catalog of human KRABassociated zinc nger genes: insights into the evolutionary history of a large family of transcriptional repressors . Genome Res. , 16 ( 5 ), 2006 , 669 { 677 .

10. S. B. Needleman , C. D. Wunsch : A general method applicable to the search for similarities in the amino acid sequence of two proteins . J. Mol. Biol ., 48 ( 3 ), 1970 , 443 { 453 .

11. K. Nowick , C. Fields , T. Gernat , D. Caetano-Anolles , N. Kholina , L. Stubbs: Gain, loss and divergence in primate zinc- nger genes: a rich resource for evolution of gene regulatory di erences between species . PLoS ONE 6 ( 6 ), 2011 , e21553 .

12. K. Nowick , A. T. Hamilton , H. Zhang , L. Stubbs: Rapid sequence and expression divergence suggest selection for novel function in primate-speci c KRABZNF genes . Molecular Biology and Evolution 27 ( 11 ), 2010 , 2606 { 2617 .

13. P. Rice , I. Longden , A. Bleasby: EMBOSS: the European molecular biology open software suite . Trends Genet ., 16 ( 6 ), 2000 , 276 { 277 .

14. D. Schmidt , R. Durrett: Adaptive evolution drives the diversi cation of zinc- nger binding domains . Mol. Biol. Evol . 21 ( 12 ), 2004 , 2326 { 2339 .

15. B. Schuster-Bockler , J. Schultz , S. Rahmann: HMM Logos for visualization of protein families . BMC Bioinformatics 5 , 2004 , 7 .

16. J. H. Thomas , R. O. Emerson: Evolution of C2H2-zinc nger genes revisited . BMC Evol. Biol. 9 , 2009 , 51 .

17. UCSC Human Genome Feb . 2009 ( hg19, GRCh37) Pairwise Alignments . Online. http://hgdownload.cse.ucsc.edu/downloads.html.

18. R. Urrutia: KRAB-containing zinc- nger repressor proteins . Genome Biology 4 ( 10 ), 2003 , 231 .

19. A. J. Viterbi : Error bounds for convolutional codes and an asymtotically optimum decoding algorithm . IEEE Transactions on Information Theory IT-13 , 1967 , 260 { 267 .