The DNA Word Design Problem: A New Constraint Model and New Results

1 Introduction

This paper is about the design of error-correcting codes [MacWilliams and Sloane, 1977] for biological applications. The goal is to compute as large as possible a set S of length n code words over the alphabet {A, C, G, T} such that every pair of code words has Hamming distance (HD) at least d, and the set of additional constraints U on S is satisfied. The maximum cardinality of such a set S is denoted AU4 (n, d). The constraints in U typically include a subset of: the reverse complement constraint (RC) which requires that for each pair of words, the Hamming distance between the reverse of one and the Watson-Crick complement of the other is at least d; and a percentage constraint (P) which requires that 50% of letters in each word are from {C, G}.

The design of DNA code words has its motivation in tasks related to: information storage in DNA computing [Adleman, 1994; Lipton, 1995], storing information in DNA strands [Brenner and Lerner, 1992], DNA computing on surfaces [Frutos et al., 1997], probes in DNA micro-array technologies [Fodor et al., 1991; Schena et al., 1995], use as tracking tags [Qiu et al., 2003], and many others. DNA code

∗ Supported by the Israel Science Foundation, grant 182/13.

Vitaly Lagoon

Cadence Design Systems, USA lagoon@cadence.com

words are also used as molecular bar codes enabling to manipulate and identify individual molecules in complex chemical libraries [Shoemaker et al., 1996]

The impact of DNA codes is notable also today. Recent applications of DNA-encoded libraries to drug discovery are described, for example, in [Mullard, 2016; Goodnow Jr et al., 2017; Yuen and Franzini, 2017]. In [Buschmann and Bystrykh, 2013] the authors present recent results on errorcorrecting DNA codes. A recent application in the area of DNA data storage is described in [Blawat et al., 2016].

For the combinatorial constraints that arise in this type of problem, there are no known efficient algorithms. Techniques from coding theory have been applied [Brenner and Lerner, 1992; Frutos et al., 1997], stochastic local search is applied [Tulpan, 2000; Tulpan and Hoos, 2003; Chee and Ling, 2008], and genetic algorithms are also considered [Garzon et al., 1999; Ashlock et al., 2002; Ashlock and Houghten, 2005; 2009]. In [Montemanni et al., 2014], the authors present three different meta-heuristic approaches for these types of problems. They show improvements on some of the best known codes and present an extensive experimental comparison with previous works. Constraint and SAT programming [van Dongen, 1999; Metodi et al., 2011] have also been applied.

The instance A{P,RC} 4 (8, 4) occurs as problem CSPLib 033: Word Design for DNA Computing on Surfaces [van Dongen, 1999] and we focus special attention on this instance in the presentation of this paper. For this instance we mention the following results reported in the literature: In [Frutos et al., 1997] (1997) the authors report a solution with 108 DNA words. They search for solutions that have a specific form which they call template-map. In [Tulpan and Hoos, 2003] (2003) the authors present an algorithm based on stochastic local search. Initializing their algorithm with the 108 word set found by Frutos et al. [Frutos et al., 1997] they construct a set with 112 words in "less than one day of CPU time". Marc van Dongen also reports a solution with 112 words [van Dongen, 1999]. In [Mancini et al., 2008] (2008) the authors present a solution consisting of 87 DNA words, reported to be found in 554 seconds of CPU time using an OPL [Hentenryck et al., 1999] model. In [Metodi et al., 2011] (2011), the authors apply the BEE constraint solver and report finding a solution of 112 words "in a fraction of a second". They also show that this is the optimal code size for the DNA word design problem when using the template-map strategy of Frutos et al. [Frutos et al., 1997]. In his PhD thesis, Tulpan [Tulpan, 2000] presents a solution with 128 words. Prior to 2016, this solution was not known to the CSPLIB library of constraints. The technique presented in this paper finds a solution with 128 words in under 30 seconds of computation time. This is the first time that a SAT/CSP-based approach is able to match the best known solution for this instance.

In this paper we describe a new constraint model for the design of error-correcting codes for biological applications based on clustering of DNA words. We report solutions which improve the prior state-of-the-art.

2 The Problem Specification

For given parameters (n, d, U) (an instance of the DNA word problem) we seek a set S of cardinality |S| = m consisting of n-letter words over the alphabet Σ = {A, C, G, T } which satisfies the constraint (HD) as well as a subset U of the constraints {(RC), (P)} detailed below:

The Hamming distance constraint (HD): Each pair of distinct words in S differ in at least d positions.

The reverse complement constraint (RC): For every w1, w2 ∈ S (not necessarily distinct): wR 1 (the reverse of w1) and wc2 (the word obtained from w2 by replacing each A by T , each C by G, and vice versa) differ in at least d positions.

The percentage constraint (P): The number of letters from the set {C,G} in each word from S is between ⌊n/2⌋ and ⌈n/2⌉ (50%).

We denote by AU4 (n, d) the maximum cardinality of such a set S and we seek lower bounds of AU4 (n, d) together with corresponding sets S of words which satisfy the constraints.

For a letter x ∈ {A,C,G, T}, we denote by xc the complement letter of x: Ac = T , T c = A, C c = G, and Gc = C. We sometimes express the (RC) constraint in terms of two separate cases. (RC1): For every w ∈ S, wR and wc differ in at least d positions; and (RC2): Each pair of distinct words w1,w2 ∈ S, wR 1 and wc2 differ in at least d positions. For a word w we write P(w) and RC1(w) to specify that w satisfies the constraints (P) and (RC1). Similarly, for a pair u, v of words, we write HD(u,v) and RC2(u, v).

3 BaseWords and their Generated Images

Let Σ = {C, G, T, A} and denote Σn , the set of all n letter words over Σ. A word transformation is a mapping f : Σn → Σn . We are interested in (small) sets of word transformations which are self-inversive and which commute with respect to composition. A word transformation, f, is selfinversive if f ◦ f = id. A pair of word transformations, f, g is commutative with respect to composition if f ◦ g = g ◦ f. The composition of word transformations is always associative - a property inherited from the composition of relations. So for word transformations, f, g, h we always have (f ◦ g) ◦ h = f ◦ (g ◦ h). We say that a set of word transformations, T , is SCA if for all f, g, h ∈ T , f is self-inversive, f, g are commutative, and f, g, h are associative.

Let B ⊆ Σn and T be a set of word transformations. We denote by T (B) ⊆ Σn the closure of B under composition of transformations from T . These are all of the words that can be generated from B by repetitive application of transformations from T . Note that if T is SCA then T (B) consists of (at most) |B|×2|T | words. For simplicity, for w ∈ Σn , we write T (w) instead of T ({w}). We write T ′ (w) to denote T (w) \ {w}.

In this paper we focus on 4 simple word transformations introduced in the following example. In practice, any set of word transformations defined in terms of a self-inversive permutation of the alphabet Σ and a self-inversive permutation on the letters in a word is SCA.

Example 1 Consider the set T4 = {wc, fc, hs,ms} consisting of four word transformations defined by:

1. wc: (Watson-Crick complement) applies to a word swapping letters A ↔ T , C ↔ G.
2. fc: (full complement) applies to a word swapping letters A ↔ G, C ↔ T .
3. hs: (half swap) applies to a word by swapping the left ⌊n/2⌋ letters with the right ⌊n/2⌋ letters. If n is odd then the middle letter of the word remains in place.
4. ms: (middle swap) applies to a word by swapping letters in positions i ≤ ⌊n/2⌋ and n - i. If n is odd then the middle letter of the word remains in place.

One can check that (any subset of) T4 is SCA. We denote T3 = {wc, fc, hs}, T2 = {fc, wc}, T1 = {fc} and T0 = ∅.

Example 2 Consider the context where n = 8, the word w = C T A C G A A C and the set of word transformations T3 introduced in Example 1. The set of eight words generated from w under T3 is detailed below (to the right of each word the transformation under which it is obtained from w).

CTACGAAC w GAACCTAC hs(w) TCGTAGGT fc(w) GATGCTTG wc(w) AGGTTCGT hs(fc(w)) CTTGGATG wc(hs(fc(w))) AGCATCCA wc(fc(w)) TCCAAGCA hs(wc(fc(w)))

The following two examples illustrate that T3(w), for some words w might consist of less than eight words. In Example 3, each (individual) word in T3(w) is closed under the half swap transformation. In Example 4, the set of words T3(w) is closed under the half swap transformation.

Example 3 Consider the word w = C T T G C T T G and the set of word transformations T3 introduced in Example 1. The set T3(w) consists of the following four words (to the right of each word we detail the transformation under which it is obtained from w). Note that each word v ∈ T3(w) is closed under the half swap transformation. Namely, hs(v) = v.

CTTGCTTG w AGGTAGGT wc(fc(w)) GAACGAAC wc(w) TCCATCCA fc(w)

Example 4 Consider the word w = C A C T G T G A and the set of word transformations T3 introduced in Example 1. The set T3(w) consists in the following four words (to the right of each word we detail the transformation under which it is obtained from w). Note that T3(w) is closed under the half swap transformation. Namely, hs(T3(w)) = T3(w).

CACTGTGA w TGTCACAG wc(fc(w)) GTGACACT wc(w) ACAGTGTC fc(w)

In this paper we represent elements of Σn as length 2n bit vectors. We adopt the convention that each letter of Σ = {A, C, G, T} is represented as follows by two bits:

C = 00 G = 01 (both have msb=0) A = 10 T = 11 (both have msb=1)

Example 5 The eight words detailed in Example 2 are represented respectively by the following eight bit vectors:

00 11 10 00 01 10 10 00 01 10 10 00 00 11 10 00 11 00 01 11 10 01 01 11 01 10 11 01 00 11 11 01 10 01 01 11 11 00 01 11 00 11 11 01 01 10 11 01 10 01 00 10 11 00 00 10 11 00 00 10 10 01 00 10

For a given word w ∈ Σn in bit vector representation, and a set of word transformations T ⊆ T4 of Example 1 the set T (w) can be specified by a corresponding set of literals. Thus, when encoding to CNF a constraint model to represent the set T (w) no CNF clauses are required - all of the information is in the Boolean variables.

Example 6 Consider w ∈ Σ8 represented by the bit vector 〈x1x2 x3x4 x5x6 x7x8 x9x10 x11x12 x13x14 x15x16〉. Then, the following are the bit-vector representations of w, hs(w), fc(w), wc(w), hs(fc(w)), wc(hs(fc(w))), wc(fc(w)), and hs(wc(fc(w))), respectively:

1. 〈x1x2 x3x4 x5x6 x7x8 x9x10 x11x12 x13x14 x15x16〉 2. 〈x9x10 x11x12 x13x14 x15x16 x1x2 x3x4 x5x6 x7x8〉 3. 〈¯x1¯x2 ¯x3¯x4 ¯x5¯x6 ¯x7¯x8 ¯x9¯x10 ¯x11¯x12 ¯x13¯x14 ¯x15¯x16〉 4. 〈¯x9¯x10 ¯x11¯x12 ¯x13¯x14 ¯x15¯x16 ¯x1¯x2 ¯x3¯x4 ¯x5¯x6 ¯x7¯x8〉 5. 〈x1¯x2 x3¯x4 x5¯x6 x7¯x8 x9¯x10 x11¯x12 x13¯x14 x15¯x16〉 6. 〈x9¯x10 x11¯x12 x13¯x14 x15¯x16 x1¯x2 x3¯x4 x5¯x6 x7¯x8〉 7. 〈¯x1x2 ¯x3x4 ¯x5x6 ¯x7x8 ¯x9x10 ¯x11x12 ¯x13x14 ¯x15x16〉 8. 〈¯x9x10 ¯x11x12 ¯x13x14 ¯x15x16 ¯x1x2 ¯x3x4 ¯x5x6 ¯x7x8〉

The interested reader can double check that for w = C T A C G A A C, represented by 00 11 10 00 01 10 10 00 these eight bit vectors correspond precisely to those given as Examples 2 and 5.

Let T ⊆ T4 and U ⊆ {(P), (HD), (RC)}. We say that w ∈ Σn is a base word with respect to T and U if the words in T (w) satisfy the constraints U, and w is the minimal element of T (w) (under the lexicographic order on the bit vector representation of the words).

Example 7 The word w1 = C T A C G A A C is a base word with respect to T3 and U = {(P), (HD), (RC)}. The eight words of T3(w1) detailed in Example 2 satisfy the constraints (P), (HD), and (RC), and the bit vector representation of w1 is minimal in the lexicographic order of the bit vector representations detailed in Example 5. The word w2 = C T T G C T T G is also a base word with respect to T3 and U. The four words of T3(w2) detailed in Example 3 satisfy the constraints, and w2 is minimal (in the lexicographic order of the respective bit vector representations).

Our strategy to solve the DNA word design problem is to seek a solution S constructed in terms of a much smaller set, B, of base words such that S = T (B) for a set of word transformations T . In this sense, one can view T as inducing a clustering of the words in Σn . However, it is not just about the size of the set B. Expressing the constraints of the DNA word design problem in terms of the set B turns out to be much more compact as well. We demonstrate this in the next three lemmas where we focus on the case where T ⊆ T4.

Let x1 . . . x2n be a bit vector representing a word w ∈ Σn and B = T (w) represented as a set of bit vectors. A model that constrains w to be a base word must constrain w to be the lexicographic minimum in T (w). That means imposing at most 16 constraints (assuming that T ⊆ T4). A model that constrains w to be a base word must impose constraints (P) and (RC1) on each word of T (w). The following lemma implies that it suffices to impose them only on w.

Lemma 1 Let w ∈ Σn . If w satisfies Constraint (P), then so does T (w). If w satisfies Constraint (RC1), then so does T (w).

Proof: Let w ∈ Σn and f ∈ {hs,wc, fc,ms}. We show that if w satisfies the constraint, then so does f(w). Denote w = 〈x1, . . . , xn〉 and f(w) = 〈y1, . . . , yn〉 and assume the constraints hold for w. The proof for constraint (P): If f ∈ {hs,wc,ms} then w and f(w) have the same number of letters from {C,G}. If f = fc and n is even, then the number of {C,G} in f(w) equals the number of {A, T} in w, which by constraint (P), equals the number of {C,G} in w. For odd n the proof follows with a similar argument. The proof for constraint (RC1): consider the Boolean values bi ↔ xi ̸= xcn-i+1 and di ↔ yi ̸= ycn-i+1 for 1 ≤ i ≤ n. For each choice of f , direct inspection reveals that ∑ni=1 bi = ∑ni=1 di. So (RC1) holds for f(w). □

A model that constrains w to be a base word must impose constraints (HD) and (RC2) on each pair of distinct words from T (w) i.e., two constraints for each pair of distinct words from T (w). In the case that T (w) contains k words that means at most (k2) constraints in total. The following lemma implies that it suffices to consider only the pairs including the base word w i.e., at most k - 1 constraints only.

Lemma 2 Let w ∈ Σn . If all pairs (u, w), where u ∈ T (w) and u ̸= w satisfy Constraints (HD) and (RC2), then T (w) satisfies Constraints (HD) and (RC2).

Proof: We detail the proof for constraint (HD). The proof for constraint (RC2) is similar. Assume that w is a base word and that w differs from every other word w′ ∈ T4(w) in at least n/2 positions. We argue (∗) that for any pair of distinct words u, v ∈ T4(w), also u and v differ in at least 4 positions. By definition of T4(w), u = f(w) where f is some composition of functions from {fc,wc, hs,ms}. The pair u, v differ in the same number of positions as the pair f(u), f(v). The argument (∗) follows because the functions composing to f are all inversive and so f(u) = w. But we already know that w and f(v) differ in at least four positions. □

Now consider a set B of base words. A model that constrains S = T (B) to be a solution of the DNA word design problem must impose constraints (HD) and (RC2) on pairs of words: for each w1,w2 ∈ B (viewing B as a sequence we can assume that w1 occurs before w2) all pairs u, v where u ∈ T (w1) and v ∈ T (w2). The following lemma implies that it suffices, for the pair (w1,w2) to impose the constraints for pairs of the form (w1, v) with v ∈ T (w2) and ignoring the symmetric ones of the form (u,w2). That means checking a linear number of constraints, two for each pair (u,w2), instead of a quadratic number.

Lemma 3 Let w1,w2 ∈ Σn . If all pairs (w1, v) with v ∈ T (w2) satisfy Constraints (HD) and (RC2), then all pairs u, v in T (w1) × T (w2) satisfy Constraints (HD) and (RC2).

Proof: The proof is similar to that of Lemma 2. □

4 The Constraint Model

Let (n, d, U) specify an instance of the DNA word problem. To ease presentation, we introduce the model in three steps.

Step 1: Consider the search for a solution S given in terms of m base words with respect to a single set of word transformations F1 ⊆ T4. Let M be a m × 2n Boolean matrix. Each row in M represents an n-letter base word (in 2n bits) and represents 2|F1| different words. We view a matrix as the set of its rows. We write w ∈ M to indicate that w is a row of M. If w,w′ ∈ M we write w <M w′ to indicate that row w occurs before row w′ in M. We denote by ≺ the lexicographic order on Boolean vectors. The solution S represented by matrixM is S = ∪{ F1(w) | w ∈M }.

For each base word w ∈ M, each single word u in F1(w) is required to satisfy constraints P(u) and RC1(u). Moreover, each pair of distinct words u and v in F1(w) is required to satisfy constraints HD(u, v) and RC2(u, v). Here, according to Lemma 1 and 2, it is sufficient to require P and RC1 only for the base word w, and it is sufficient to require HD and RC2 only for pairs that include w.

Observe that Constraint (1) implies that for each base word w, a row in M, the words in F1(w) are pairwise distinct.

The next constraint ensures that words u, v generated from different base words w, w′ ∈ M satisfy constraints HD(u, v) and RC2(u, v). Here, according to Lemma 3, it is sufficient to consider pairs which include w.

The next constraint in our model is a symmetry breaking constraint. It ensures that words inM are minimal in the sets they generate under F1, and that the matrix of base words is sorted.

Example 8 Consider the instance (12, 8,∅), in which we seek a solution S that contains 12-letter words over the alphabet {A,C,G, T} such that every pair w1,w2 ∈ S has Hamming distance at least 8. Consider a clustering based on T4 as in Example 1 where each cluster consists of 16 DNA words. The following 7 base words (translated from the corresponding 5 × 24 Boolean matrix) provide a solution with 7 × 16 = 112 DNA code words, improving the previous best known result of 64 due to [Bogdanova et al., 2001].

The constraint model includes Constraints (1-3) on each of the two matrices, and additional constraints on pairs of words u ∈ F1(w1) and v ∈ F2(w2) where w1 ∈ M1 and w2 ∈ M2 satisfy constraints HD(u, v) and RC2(u, v).

In some cases we can do better as indicated by the following generalization of Lemma 3.

Lemma 4 Let T ′ ⊆ T be sets of word transformations and w1,w2 ∈ Σn . If all pairs (w1, v) with v ∈ T (w2) satisfy Constraints (HD) and (RC2), then all pairs u, v in T ′ (w1) × T (w2) satisfy Constraints (HD) and (RC2).

So, if F1 ⊆ F2, we can introduce a refined constraint to the model:

Example 9 Consider the setting of Example 8 and a clustering based on T4 and T3 as in Example 1 where each cluster consists of 16 or 8 DNA words respectively. The following base words in M1 and M2 (translated from the corresponding 5 × 24 and 2 × 24 Boolean matrices) provide a solution with 7×16+2×8 = 128 DNA code words, doubling the previous best known result of 64 due to [Bogdanova et al., 2001].

5 Implementing the Constraint Model

The computations reported in this paper are performed using the finite-domain constraint compiler BEE [Metodi et al., 2013] which compiles constraints to a CNF, and solves it applying an underlying SAT solver. In the configuration for this paper we use CryptoMiniSAT v2.5.1 as the underlying SAT solver. All computations were performed on a cluster of Intel E8400 cores, each clocked at 2 GHz. Each of the cores in the cluster has computational power comparable to a core on a standard desktop computer. Each SAT instance is run on a single thread.

A solution is a tuple of matrices M1, . . . ,Mk as described in step 3 of Section 4, representing m1, . . . ,mk base words respectively. These are then applied to generate a set S of corresponding DNA words which are a solution to the DNA word design problem.

Several optimizations are applied in the implementation. Constraint (3), introduced above, is a first symmetry break. Base words are selected to be the minimal words in the corresponding clusters of words generated by F1, . . . ,Fk, and the rows of the matrices M1, . . . ,Mk are constrained to be sorted. Columns in the matrices cannot be arbitrarily reordered. However there is a symmetry that allows to swap the letters at positions i and n + 1 - i in the n-letter base words and we have introduced this symmetry break.

For experimentation, an instance takes the form (n, d,U) together with a tuple of pairs 〈(m1,F1), . . . , (mk,Fk)〉 where F1, . . . , Fk ⊆ T4 are sets of word transformations and we seek a solution with m1 base words with respect to F1 and m2 base words with respect to F2 etc.

Table 1 summarizes our experimentation and details the best lower bounds for AU4 (n, d) that we have found using our approach. The first three columns detail the instance (n, d,U). The forth column ("lit") details the best previous result that we found in the literature together with a reference indicating where this result can be found in literature as detailed in the legend of Table 1. The fifth column ("new") details our result. The sixth column ("clustering") details the word transformations that were used and the number of base words per set of transformations (e.g, 〈(T4, 4), (T0, 3)〉 means that 4 base words were found for the word transformations in T4 and 3 base words were found for the word transformations in T0, summing up to a total of 4×16+3×1 = 67 DNA code words). The last three columns describe the CNF encoding (number of clauses and variables) and the SAT solving time required to find the corresponding DNA code.

In Table 1 results are for a given type of clustering. We have not addressed in this paper the issue of how to find a suitable type of clustering. For each instance of the DNA word problem we have experimented with a wide range of possible choices applying a mixture of intuition and brute force search. We only present here the instances for which we have improved on the best previous result, with the exception of instance (8, 4, {hd, rc, p}) for which we obtained a solution with 128 words, equal in size to the result by Tulpan [Tulpan, 2000]. In this case, we report the result as it is obtained in less than 30 seconds of computation, as opposed to the result by Tulpan which takes orders of magnitude longer.

Tables 2, and 3 illustrate some of the results in the search process for a solution for the instances (7, 3, {hd, rc, p}) and (9, 5, {hd, rc, p}). For these instances, there is no solution that contains a cluster of type T4. The first four columns detail the number of clusters of type T3, T2, T1 and T0. The fifth column details the total number of words represented by the corresponding clustering. The final column indicates the time (in seconds) to find a solution with the corresponding clustering. Here we assumed a timeout of 24 hours. Table 4 illustrates the search for the instance (12, 8, {hd}). Here solutions include also clusters of type T4.

6 BottomLine and Discussion

We have improved on several lower bounds AU4 (n, d) in the context of error-correcting codes for DNA word problems. Our approach is constraint based and expressed in terms of a partition of the class of all n-letter words into clusters of words as determined by subsets of a set T of word transformations. Each cluster is required to satisfy the constraints of the word problem. A pair of clusters that contain conflicting words cannot both be in a solution.

To evaluate the quality of our proposed clustering of DNA words let us consider the instance (8, 4, {(P), (HD), (RC)}) with a clustering based on T3. There are 48 = 65,536 DNA words (of length eight). From these, there are only 17,440 words which satisfy constraints (P) and (RC1) and hence can appear in a solution to the DNA word design problem. These 17,440 words partition into a total of 2262 clusters when applying T3. From these 2262 clusters only 1544 satisfy constraints (HD) and (RC2). From these 1544, 1424 are clusters of size eight, and 120 are of size four. The 1544 clusters that could appear in a solution comprise 11,872 DNA words. This is 68% of the 17,440 that could appear in a solution. We also note that there are 2262 - 1544 = 718 clusters that satisfy constraints (P) and (RC1) but does not satisfy constraints (HD) or (RC2). We observe that these "bad" clusters, are not so due to a single contradicting pair of words: by applying Lemma 2 it happens that every word in a "bad" cluster is in conflict with at least one other word in that cluster. In clusters of size four there is an average of 2 contradicting pairs. In clusters of size eight there is an average of 6.06 contradicting pairs. Furthermore, when a pair of "good" clusters (from the 1544) does not satisfy either (HD) or (RC), then this is on average because of 8.68 contradicting pairs. Indeed, by applying Lemma 3 it is obtained that in this case every word from one cluster is in conflict with at least one word from the other cluster.

As can be seen in Table 1, most of the new codes obtained applying our technique are found in under 10 minutes of computation time. Several of the codes require considerable more time. All of the codes are available on request from the authors and will be made publicly available,

Finally, we note that our approach is incomplete. There may be (larger) solutions which cannot be found using our approach as they cannot be expressed in terms of the specific clusters we have selected.