Oprocha, Piotr

Schedae Informaticae

English

Jagiellonian Univeristy Repository

S C H E D A E I N F O R M A T I C A EVOLUME 15 2006Fast Solutions for DNA Code Words TestPiotr OprochaAGH University of Science and Technology, Faculty of Applied MathematicsAl. Mickiewicza 39, 30-059 Krako´w, Poland, andJagiellonian University, Institute of MathematicsReymonta 4, 30-059 Krako´w, Polande-mail: oprocha@agh.edu.plAbstract. An essential part of any DNA computation is to encode dataon DNA strands. Performing such computation we must be sure thatno undesirable hybridization will occur. The aim of this paper is topresent algorithms which can test if a given set of code words satisfiescertain coding and involutory properties, as it is necessary to preventundesirable DNA strands interactions. These algorithms are based onthe theory of codes and pattern matching methods.If we denote by m sum of the lengths of words in a given set, then allpresented algorithms can be realized with O(m2) complexity both fortime and space.Keywords: DNA, strings properties, code, complexity, algorithms.1. IntroductionEvery DNA molecule is an oriented sequence of nucleotides, so it can berepresented as a single word over alphabet ∆ = {A,C,G, T}. It is importantproperty of DNA that C is complementary to G, A is complementary to Tand two complementary sequences (single DNA strands) with opposite orien-tations joint together (hybridize) forming a double stranded DNA molecule(this property is called Watson-Crick complementarity).42It is interesting that DNA may be used in computational problems. AsDNA strands join together almost in a blink, it give the possibility to solvesome complicated problems (e.g. the existence the Hamilton Cycle problem)very fast. There are many technical problems with such experiments, howeverthere is hope for some improvements in future. Main ideas and results in DNAcomputing the reader may find in [6].In laboratory experiments it is very important to avoid some unnecessaryhybridizations (see Fig. 1). A situation, when a part of a given DNA sequence(longer than k) is complement to the other part of this sequence with the samelength, is undesired. For a given set of DNA sequences, a situation, when onesequence is a reverse complement of the subword of another sequence or it isa reverse subword of concatenation of other two sequences, is also unliked.Fig. 1. Unallowed hybridizationsPapers [9, 6] define properties of words and languages, which if fulfilled bya given set of DNA strands, prevent situations described above. Authors of[7, 10] present methods and try to implement theoretical results in softwareprogram. The most important part of this program is to perform tests ofgenerated sequences. The algorithm introduced there work in O(m4) time,where m denotes the sum of the lengths of DNA sequences in a given set.In this paper we describe how to construct algorithms with complexityO(m2), which speed up whole process of DNA code words generation. In thepresented algorithms the theory of codes is used as well as string matchingalgorithms like KMP or Aho-Corasick algorithm (see [3, 4]).2. Basic definitions and propertiesAn alphabet Λ is any nonempty finite set. We will denote by ∆ thealphabet {A,C,G, T} representing the DNA nucleotides.The word w ∈ Λ∗ is a prefix (resp. suffix) of the word x ∈ Λ∗, which wedenote w ⊏ x (resp. w ⊐ x), if there exists word u ∈ A∗ such that wu = x(resp. uw = x). The length of the word w = a1 . . . an we denote by |w| = n.43Definition 1. Let Λ be an alphabet. A subset X of the free monoid Λ∗is a code over Λ if, for all n,m ≥ 1 and for all x1, . . . , xn, x′1, . . . , x′m ∈ X,the conditionx1x2 . . . xn = x′1x′2 . . . x′mimpliesn = m and xi = x′i for all i = 1, . . . , n.The flower automaton of X is by the definition the automaton A(X) =(Q,E, I) containing a set of states Q, the labelling (transition) E ⊂ Q×Λ×Qand state I which is the initial and terminal state at once. If we denoteI = {∗} then we may define the set of states asQ = {∗} ∪ {qwi : w = a1 . . . an ∈ X, i = 2, . . . , n}and then the transitions in A(X) areE = {(∗, a1, qw1 ) | w = a1 . . . an ∈ X} ∪ {(qwn , an, ∗) | w = a1 . . . an ∈ X} ∪{(qwi , ai, qwi+1) | w = a1 . . . an ∈ X, i = 2, . . . , n − 1}.We also define labelling function λ : E → Λ such that λ((q1, a, q2)) = a. Thelabelling function can be extended to paths labelling in the usual way.The automaton A(X) is unambiguous if for all p, q ∈ Q and w ∈ Λ∗ thereis at most one path form p to q with label w.The square of A(X) is the automaton A′(X) = (Q′, E′, I ′) such thatQ′ = Q×Q, I ′ = {∗′} = {(∗, ∗)} andE′ ={((p, q), a, (p′, q′)) | (p, a, p′), (q, a, q′) ∈ E, a ∈ Λ}.Proposition 1. (see [2, chapt. IV thm. 2.1]) A finite set of words X isa code if and only if the flower automaton A(X) is unambiguous.The following proposition follows from the definitionsProposition 2. The flower automaton A(X) is unambiguous if and onlyif there is no path in A′(X) from a state (p, p) to a state (q, q) visiting a state(r, s) with r 6= s.Corollary 1. For any state (q, q) ∈ Q′ there is a path from ∗′ to (q, q)and a path from (q, q) to ∗′.44Corollary 2. A finite set of words X is a code if and only if there isno cycle in A′(X) from a state ∗′ visiting a state (r, s) with r 6= s.We recall that oriented graph G is a pair (V,E), where V is any finiteset, and E ⊂ V × V . For any u ∈ V we will denote by Adj(u) the setAdj(u) = {v ∈ V : (u, v) ∈ E}. If G = (V,E) then its transposition is thegraph GT = (V,ET ), where ET = {(v, u) : (u, v) ∈ E}.An involution θ : Λ → Λ of a set Λ is a mapping such that θ2 = idΛ.We will define two special involutions ν and ρ of the set ∆∗. The mappingν is defined by ν(A) = T , ν(T ) = A, ν(C) = G and ν(G) = C, and may beextended to an involution of ∆∗. The mapping ρ is defined inductively asρ(s) = s and ρ(us) = ρ(s)ρ(u) = sρ(u) for any s ∈ ∆ and u ∈ ∆∗.Definition 2. Let θ be an involution of a set Λ and let X ⊂ Λ∗ be afinite set.1. The set X is called θ(k,m1,m2)-subword compliant if for all u ∈ Λ∗such that |u| ≥ k we have Λ∗uΛmθ(u)Λ∗ ∩X = ∅ for m1 ≤ m ≤ m2.2. The set X is called θ-compliant if Λ∗θ(x)Λ∗ ∩X = ∅ for all x ∈ X.3. The set X is called θ-free if X2 ∩ Λ∗θ(x)Λ∗ = ∅ for all x ∈ X.4. The set X is called strictly θ if θ(x) /∈ X for all x ∈ X.Corollary 3. The set X is θ(k,m1,m2)-subword compliant if and onlyif for all u ∈ Ak we have A∗uAmθ(u)A∗ ∩X = ∅ for m1 ≤ m ≤ m2.Example 1. Let us consider following sets of wordsX1 = {AA,TT,CTT} , X2 = {AA,TAC,GT} .Obviously both defined sets are prefix codes. Let X˜i = θ(Xi), where θ denotesmorphism θ = νρ. In this caseX˜1 = {AA,TT,AAG} , X˜2 = {TT,GTA,AC} .Observe that X1 is not strictly θ as θ(AA) = TT ∈ X1. If we modifyX1 to X′1 = X1 \ {AA} then it becomes strictly θ but not θ-compliant asGθ(TT ) = GAA ∈ X ′1. Set X2 is strictly θ and θ-compliant but it is notθ-free. We have that (X2)2 = {AAAA,TACAA,GTAA, . . . , GTGT} andGTAA ∈ θ(TAC)A ⊂ A∗θ(TAC)A∗. Set X2 is θ-compliant so it is alsoθ(k,m1,m2)-subword compliant for any parameters k,m1,m2.45It is an important property of a single-stranded DNA molecule that itstwo ends are physically different. By convention, one end is called 5′ endand the other one is called the 3′ end. Let u ∈ ∆∗ denote DNA strand inits 5′ → 3′ orientation. If we denote by ←−u the Watson-Crick complement ofthe word u (i.e. oppositely-oriented complementary DNA strand which canbind to u forming double-stranded DNA), also in orientation 5′ → 3′, thenit is easily seen that ν(ρ(u)) = ρ(ν(u)) = ←−u (word ←−u consists of symbols(nucleotides) complement to symbols of u and its order is inverted).We would like to test if a given set of words over alphabet ∆ (DNAstrands) is νρ-free and νρ-subword compliant (fulfills conditions (1)–(4) forinvolution θ = νρ). These properties describe intermolecular hybridizationswhich are unlikely for us. We will try to construct fast algorithms checkingall these properties.3. AlgorithmsWe present algorithms checking if a given set is a code, θ-compliant andθ-free. We will also show that the upper bound for the complexity functionof presented algorithms is O(n2). All algorithms will be presented in thesimplified pseudocode.3.1. A test for a codeBy Corollary 2 it is enough to check if there exists path from ∗′ to ∗′ goingthrough state (r, s), where r 6= s. It is also important that paths throughstates (p, p) need not to be checked as the next state on the path is ∗′ or(q, q) for some q.The number of states n and the number of transitions m in the flowerautomaton A(X) for a code X ={x1, x2, . . . , x|X|}are equal ton = 1 +|X|∑i=1(|xi| − 1) and m =|X|∑i=1|xi|.We will assume that |x| > 1 for each x ∈ X. With this assumptionm < 2n and automaton is “thin” so we will represent the flower automatonby transition lists. In case of DNA we have only 4 symbols, so we will useone list for each symbol a ∈ Λ.46First, we present an algorithm which for a given set of words X builds aflower automaton A(X).Flower-Automaton(X)1: q ← 02: for each x ∈ X3: do Add-Trans(Q[0],x1,Q[q + 1])4: Add-Trans(Q[q + |x| − 1],x|x|,Q[0])5: for i← 2 to |x| − 16: do q ← q + 17: Add-Trans(Q[q],xi,Q[q + 1])8: return QStates of the A(X) are stored in the table Q and state ∗ is representedby Q[0]. Function Add-Trans(q, a, q′) add transition (q, a, q′) to the list oftransitions of the state q. This function may be realized in O(1), so the flowerautomaton is build in linear time O(n+m) = O(m).The next step is to build an A′(X) automaton. As we are only interestedin the existence of a special cycle visiting ∗′ it is enough for us to construct anoriented graph G with the vertex set equal to Q′. There is an edge betweentwo vertices q and p if there exists a ∈ Λ such that transition (q, a, p) ∈ E′. Infact there exists at most one such a. Vertices of graph G with its adjacencylists are stored in table V [i][j], where i, j = 1, . . . , n.Square-Graph(Q)1: n← |Q|2: for i← 0 to n− 13: do for j ← 0 to n− 14: do Add-Edges(i, j)5: return V47Function Add-Edges(i, j) works as follows. For each a ∈ Λ it takestransition list L(a, i) of the state Q[i] and transition list L(a, j) of the stateQ[j]. If transition (Q[i], a,Q[r]) ∈ L(a, i) and (Q[j], a,Q[s]) ∈ L(a, j) thenit adds an edge from V [i][j] to V [r][s] to the adjacency list of vertex V [i][j].Observe that every pair of transitions is considered only once. It implies thatcomplexity of the procedure Square-Graph(Q) is equal to O(n2 +m2) =O(m2).The last part of the algorithm is construction of the procedure whichdetermines if there exists a cycle starting in vertex V [0][0] and visiting ver-tex V [i][j] for some i 6= j. We will use some modification of the algorithmStrongly-Connected-Components presented in [3] to check this prop-erty. We will also use procedure Initialise-Colours which assigns colourwhite to vertices V [i][i] and gray to vertices V [i][j], where i 6= j. It obviouslymay be realized in O(m2) operations. Colours of vertices are stored in tablecolour. We will also use temporary table visited.Dfs-Visit(u, V )1: visited[u]← true2: for each v ∈ Adj(u)3: do if visited[v] = false4: then if colour[v] 6= black5: then if colour[v] = gray6: then colour[v]← black7: if Dfs-Visit(v, V ) = false8: then return false9: else return false10: return trueCode-Test(X)1: Q← Flower-Automaton(X)2: V ← Square-Graph(Q)3: Initialise-Colours(V)4: u← V [0][0]485: for each v ∈ V6: do visited[v]←false7: Dfs-Visit(u, V )8: change G to GT9: for each v ∈ V10: do visited[v]← false11: return Dfs-Visit(u, V )Function Dfs-Visit(u) is a simple modification of standard DFS. Whenit reaches gray vertex its colour is changed to black, which means that thereexists path from V [0][0] to this vertex. In the first use of DFS it is impossibleto reach a black unvisited vertex. White vertices are not considered as therecould happen anything “bad” on the path form such vertex to ∗′ or from∗′. If there exists a path in GT from V [0][0] to an unvisited vertex colouredblack, Dfs-Visit will return false and then given set X is not a code.For graph G defined by adjacency lists transposition graph GT may becalculated with complexity O(|V [G]| + |E[G]|) (see [3, chapt. 23]), whichin our case equals to O(n2 +m2) = O(m2). It means that function Code-Test(X) has complexity bounded by O(m2). Memory complexity is alsobounded by O(m2). Another approach to this problem (with the same com-plexity) is presented in [4]. Both approaches are faster than the algorithmintroduced in [7] which has time complexity O(n4) = O(m4).3.2. A test for θ-compliance and a strictly θ testAt the beginning we must recall KMP algorithm as it is necessary for usto match patterns in linear time.Suppose we are given text T = [1, . . . , t] and pattern P = [1, . . . , p]. Wewill use notation Pq = P [1..q]. The basic part of this algorithm is to computethe prefix function pi defined aspi[q] = max {k : k < q , Pk ⊐ Pq} .We will useCompute-Prefix-Function to compute pi[i] for i = 1, . . . , p.It may be done with complexity O(p) (see [3, chapt. 34, p. 975]). This func-tion is the main part of KMP-Matcher which decides if P is a substringof T . This function needs O(p+ t) time and O(p) space [3].49We will use the naive algorithm to check if a given set X is strictly θ, asit seems that there is no faster method.Strictly-θ-Test(X)1: for each x ∈ X2: do w ← θ(x)3: for each u ∈ X4: do if |u| = |w| and KMP-Matcher(w, u)5: then return false6: return trueIt is easy to see that complexity of this algorithm is bounded byO(|X|∑i=1(|xi|+m)) = O(|X| m) ≤ O(m2).To check θ-compliance we need only to remove condition |w| = |u| fromthe line 4 of the Strictly-θ-Test function.Next we will construct an algorithm testing if X is θ(k,m1,m2)-subwordcompliant. Main idea of this algorithm is given by Corollary 3. We will usemodified KMP-Matcher which sets M [i] with true if T [i . . . i+ |P |−1] =P and false otherwise, for all i = 1, . . . , |T |. We will call such table Mmatch table.KMP-Compute-Match-Table(T, P, pi,M)1: for i← 1 to |T |2: do M [i] = false3: q ← 04: for i← 1 to |T |5: do while q > 0 and P [q + 1] 6= T [i]6: do q ← pi[q]7: if P [k + 1] = T [i]8: then q ← q + 1509: if q = |P |10: then M [i] = true11: q ← pi[q]Subword-Compliance-Test(k,m1,m2,X)1: for each x ∈ X2: do for i← 1 to |x| − k3: do u← x[i . . . i+ k − 1]4: w ← θ(u)5: piu ← Compute-Prefix-Function(u)6: piw ← Compute-Prefix-Function(w)7: for each y ∈ X8: do KMP-Compute-Match-Table(y, u, piu,Mu)9: KMP-Compute-Match-Table(y,w, piw,Mw)10: i, j ← 111: while i ≤ |y| and j ≤ |y|12: do if m2 < j − i− k or Mu[i] = false13: then i← i+ 114: if m1 > j − i− k or Mw[j] = false15: then j ← j + 116: if m1 ≤ j − i− k ≤ m2 and17: Mu[i] = true and Mw[j] = true18: then return false19: return true51Function Subword-Compliance-Test returns true if a given set isθ(k,m1,m2)-subword compliant and false otherwise. Observe that lines 8to 18 of the procedure are realized in linear time O(|y|). We may assumethat time complexity of these lines is bounded by Cy and constant C isindependent of y. We may also assume that k < m otherwise procedurestops at line 2. Lines 3 to 6 need 4k operations. If we take C ′ = C + 4 thencomplexity of the whole function may be calculated as follows|X|∑i=1|xi|∑j=1(4k +∑y∈XCy) ≤|X|∑i=1|xi|∑j=1(4 + C)m =|X|∑i=1C ′|xi|m = C′m2 = O(m2).3.3. A test for θ freedomThe problem of θ freedom detection may be considered as the exactmatching problem with a set of patterns. Let us denote by P the set of allinvolutions of elements of the set X, i.e. P = {θ(x) : x ∈ X} = {p1, . . . , pk},and let T denote the set T = X2 = {uv : u, v ∈ X} = {t1, . . . , tl}, wherek = |X| and l = k2. Observe that X is θ-free iff pi is not a subword of tj forany i = 1, . . . , k and j = 1, . . . , l.To solve the problem, we will use Aho-Corasick algorithm. Main idea ofthis algorithm is to build a keyword tree of the set P (see Fig. 2).AAACCCC GGGTTTTTFig. 2. Keyword tree for P = {agtaa,ctgct,agct,agtac,ctggt}Let L(v) denote labelling of the path from root of K to node v. Definelp(v) as the length of the longest proper suffix of word L(v) that is prefixof some pattern in P. Then for node v let nv denote the unique node in Klabelled with the suffix of L(v) with length lp(v). When lp(v) = 0 then nv isthe root of K. We call ordered pair (v, nv) a failure link (see Fig. 3).In fact a keyword tree with failure links may be constructed in linear time(see [5, chapt. 3.4.]), i.e. the upper bound of time complexity is equal toO(∑pi) = O(m).52AAACCCC GGGTTTTTFig. 3. Keyword tree showing the failure linksGiven a keyword tree with failure links and word tj we need only O(|tj |)time to verify if there exists i such that pi is a subword of tj (see [5]). Weare searching for pattern matching in the set of texts T , so in our case timecomplexity of the searching process is bounded byO(l∑i=1|ti|) = O(|X|∑i,j=1|xi||xj |) = O(m2).Then θ freedom test complexity has O(m2) as the upper bound (O(m) fortree construction and O(m2) for words matching).4. References[1] Aho A.V., Hopcroft J.E., Ullman J.D.; The Design and Analysis of ComputerAlgorithms, Addison Wesley, 1974.[2] Bertsel J., Perrin D.; Theory of Codes, Academic Press, Inc. Orlando Florida,1985.[3] Cormen T.H., Leiserson C.E., Rivest R.L.; Introduction to Algorithms, Cam-bridge 1991.[4] Crochemore M., Rytter W.; Text Algorithms, Oxford University Press, 1994.[5] Gusfield D.; Algorithms on Strings, Trees and Sequences, Cambridge 1997.[6] Jonoska N., Mahalingham K.; Languages of DNA based code words, in: J. Chien,J. Reif, (eds.), Preliminary Proceedings of the 9th International Meeting on DNABased Computers, 2003, pp. 58–68.53[7] Jonoska N., Kephart D.E., Mahalingam K.; Generating DNA Code Words, Con-gressus Numerantium, 156, 2002, pp. 99–110.[8] Jonoska N.; Trends in Computing with DNA, J. Comput. Sci. Technol., 19, 2004.[9] Kari L., Konstantinidis S., Losseva E., Wozniak G.; Sticky free and overhangfree DNA languages, Acta Informatica, 40(2), 2003, pp. 119–157.[10] Kephart D.E., LeFevre J.; CODEGEN: The Generation and Testing of DNACode Words, to appear in IEEE.Received September 30, 2004

Fast solutions for DNA code words test

https://ruj.uj.edu.pl/xmlui/bitstream/handle/item/153725/oprocha_fast_solutions_for_dna_code_words_test_2006.pdf?sequence=1&isAllowed=y

Fast solutions for DNA code words test

Abstract

Similar works

Full text

Available Versions

Jagiellonian Univeristy Repository