# Efficient Serial and Parallel Algorithms for Selection of Unique Oligos in EST Databases

Authors: Manrique Mata-Montero; Nabil Shalaby; Bradley Sheppard
Journal: Advances in Bioinformatics
 (2013)

---

## Abstract

Obtaining unique oligos from an EST database is a problem of
great importance in bioinformatics, particularly in the discovery of
new genes and the mapping of the human genome. Many algorithms
have been developed to find unique oligos, many of which are much
less time consuming than the traditional brute force approach. An
algorithm was presented by Zheng et al. (2004)
which finds the solution of the unique oligos search problem efficiently.
We implement this algorithm as well as several new algorithms based
on some theorems included in this paper. We demonstrate how, with
these new algorithms, we can obtain unique oligos much faster than
with previous ones. We parallelize these new algorithms to further
improve the time of finding unique oligos. All algorithms are run on
ESTs obtained from a Barley EST database.


---

## Body

## 1. Introduction

Expressed Sequence Tags (or ESTs) are fragments of DNA that are about 200–800 bases long generated from the sequencing of complementary DNA. ESTs have many applications. They were used in the Human Genome Project in the discovery of new genes and are often used in the mapping of genomic libraries. They can be used to infer functions of newly discovered genes based on comparison to known genes [1].An oligonucleotide (or oligo) is a subsequence of an EST. Oligos are short, since they are typically no longer than 50 nucleotide bases. Oligos are often referred to in the context of their length by adding the suffix “mer”. For example, an oligo of length 9 would be referred to as a 9-mer. The importance of oligos in relation to EST databases is quite significant. An oligo that is unique in an EST database serves as a representative of its EST sequence. The oligonucleotides (or simply oligos) contained in these EST databases have applications in many areas such as PCR primer design, microarrays, and probing genomic libraries [2–4].In this paper we will improve on the algorithms presented in [2] to solve the unique oligos search problem. This problem requires us to determine all oligos that appear in one EST sequence but not in any of the others. In addition, we will consider two oligos to be virtually identical if they fall within a certain number of mismatches from each other. In the  appendix we include all the algorithms used and developed in this paper.


## 2. The Unique Oligos Search Problem

In this paper we use the notationHD(x,y) to denote the Hamming Distance between the strings x and y. Given an EST database D={x1,x2,…,xk}, where xi is a string over the alphabet {A,C,G,G}, integers d and l, and l-mer y, we say that y occurs approximately in D if there exists a substring z of some EST xi such that HD(y,z)≤d. We also say that an m-mutant list of a string s is a list of all possible strings, s*, of length |s| over the alphabet {A,C,G,T} such that HD(s,s*)≤m. Such a string s* is referred to as an m-mutant of s. A unique oligo of D is defined as an l-mer u such that u occurs exactly in one EST and does not occur approximately in any other EST. The unique oligos search problem is the problem of finding all unique oligos in an EST database.Many algorithms have been presented to solve this problem [5, 6]. The algorithm presented in [2] relies on an observation that if two l-mers agree within a specific Hamming Distance, then they must share a certain substring. These observations are presented in this paper as theorems.Theorem 1.
Suppose one has twol-mers l1 and l2 such that HD(l1,l2)≤d. If one divides them both into ⌊d/2⌋+1 substrings, l11l12⋯l1⌊d/2⌋+1 and l21l22⋯l2⌊d/2⌋+1, and each lji, except possibly lj⌊d/2⌋+1, has length ⌈l/(⌊d/2⌋+1)⌉, then there exists at least one i0∈{1,2,…,⌊d/2⌋+1}, such that HD(l1i0,l2i0)≤1.Proof.
Suppose by contradiction that for anyi∈{1,2,…,⌊d/2⌋+1}, l1i and l2i have at least 2 mismatches. Then HD(l1,l2)≥d+2 which is a contradiction to the fact that HD(l1,l2)≤d.Using this observation, an algorithm was presented in [2] which solves the unique oligos search problem in time O((l-q)qr24q). The algorithm can be thought of as a two-phase method. In the first phase we record the position of each q-mer in the database into a hash table of size 4q. We do so in such a way that for each q-mer x over the alphabet {A,C,G,T} we have that hashtable[hashfunction[x]]={{s1,p1},{s2,p2},…,{sn,pn}} whereby si is an EST sequence, pi is the position of x within that sequence, and n is the number of occurrences of x in the database. In the second phase, we extend every pair of identical q-mers into l-mers and compare these l-mers for nonuniqueness. We also do the same for pairs that have a Hamming Distance of 1. If they are nonunique, we mark them accordingly. Theorem 1 guarantees that if an l-mer is nonunique, then it must share a q-mer substring that differs by at most one character with another q-mer substring from another l-mer. Hence, if an l-mer is nonunique, it will be marked during phase two.Assuming there aren symbols in our EST database, the filing of the q-mers into the hash table takes time Θ(qn). In phase two, we assume that the distribution of q-mers in the database is uniform; in other words, that each table contains r≈n/4q entries. Thus we have O(r2) comparisons within each table entry. Each q-mer also has a 1-mutant list of size 3q, so, we have O(qr2) comparisons for each entry in the table. Also, the time required to extend each pair of q-mers to l-mers is 2(l-q+1). Given that we have 4q entries in the hash table, we have a total time complexity of
(1)O((l-q)qr24q)=O((l-q)q(n4q)24q)=O((l-q)qn24q),
where
(2)q=l⌊d/2⌋+1.In [7], several variations of Theorem 1 are presented. We can use these theorems to generate similar algorithms with slightly different time complexities.Theorem 2.
Suppose one has twol-mers l1 and l2 such that HD(l1,l2)≤d. If one divides them both into d+1 substrings, l11l12⋯l1d+1 and l21l22⋯l2d+1, and each lji, except possibly ljd+1, has length ⌈l/(d+1)⌉, then there exists at least one i0∈{1,2,…,d+1}, such that l1i0=l2i0.Proof.
Suppose by contradiction that we cannot find anyi0∈{1,2,…,d+1} such that l1i0=l2i0. Then there exists at least one mismatch between l1i and l2i for each i∈{1,2,…,d+1}, and thus we have at least d+1 mismatches which contradicts the fact that HD(l1,l2)≤d.Based on Theorem2 we can design a second algorithm that works in a similar way to Algorithm 1. The major difference between these algorithms is that in Algorithm 2 we are not required to do comparisons with each hash table entries mutant list. This means we have O(r2) comparisons within each table entry which yields a total time complexity of
(3)O((l-q)r24q)=O((l-q)(n4q)24q)=O((l-q)n24q),
where
(4)q=ld+1.Algorithm 1:Algorithm for the unique oligos problem.
Require: EST database D={x1,x2,…,xk}, integer l (length of unique oligos) and integer d
(maximum number of mismatches between non-unique oligos)
Ensure: All unique l-mers in D
(1)q←l/(⌊d/2⌋+1)
(2)posi← findqmers(q) (hashtable of positions of all qmers in D)
(3)for  i←1 to 4q{split loop iterations among processors}  do
(4)x←i as a base 4 integer of length q
(5)mismatchlist← list of base 4 integers of length q mismatching x by 1 digit
(6)modifiedmismatchlist← the numbers in mismatchlist in base 10
(7)mut← list of each hashtable[i] for all i∈modifiedmismatchlist
(8) goo2(q,l,d,posi[i],mut)
(9)end forAlgorithm 2:Algorithm for the unique oligos problem.
Require: EST database D={x1,x2,…,xk}, integer l (length of unique oligos) and integer d
(maximum number of mismatches between non-unique oligos)
Ensure: All unique l-mers in D
(1)q←l/(d+1)
(2)posi← findqmers(q) (hashtable of positions of all qmers in D)
(3)for  i←1 to 4q{split loop iterations among processors}  do
(4) goo(q,l,d,posi[i])
(5)end forA third theorem was also briefly mentioned [7]; however, it was not implemented in an algorithm. We use this theorem to create a third algorithm to solve the unique oligos search problem.Theorem 3.
Suppose one has twol-mers l1 and l2 such that HD(l1,l2)≤d. If one divides them both into ⌊d/3⌋+1 substrings, l11l12⋯l1⌊d/3⌋+1 and l21l22⋯l2⌊d/3⌋+1, and each lji, except possibly lj⌊d/3⌋+1, has length ⌈l/(⌊d/3⌋+1)⌉, then there exists at least one i0∈{1,2,…,⌊d/3⌋+1}, such that HD(l1i0,l2i0)≤2.Proof.
Suppose by contradiction that for anyi∈{1,2,…,⌊d/3⌋+1}, l1i and l2i have at least 3 mismatches. Then HD(l1,l2)≥d+3 which is a contradiction to the fact that HD(l1,l2)≤2.The algorithm is somewhat similar to Algorithm1. The main difference is that we compare every q-mer to q-mers in its corresponding 2-mutant list, rather than its 1-mutant list. Each q-mer has 9(q2)+3q=9q(q-1)/2+3q 2-mutants, so we have O(q2r2) comparisons for each entry in the hash table yielding a total time complexity of
(5)O((l-q)q2r24q)=O((l-q)q2(n4q)24q)=O((l-q)q2n24q),
where
(6)q=l⌊d/3⌋+1.It is important to note the4q term in the denominator of our time complexity expressions. Since this term is exponential, it will have the largest impact on the time taken to run our algorithms. Based on this observation, we expect Algorithm 3 to run the fastest, followed by Algorithm 1 and then Algorithm 2.Algorithm 3:Algorithm for the unique oligos problem.
Require: EST database D={x1,x2,…,xk}, integer l (length of unique oligos) and integer d
(maximum number of mismatches between non-unique oligos)
Ensure: All unique l-mers in D
(1)q←l/(⌊d/3⌋+1)
(2)posi← findqmers(q) (hashtable of positions of all qmers in D)
(3)for  i←1 to 4q{split loop iterations among processors}  do
(4)x←i as a base 4 integer of length q
(5)mismatchlist← list of base 4 integers of length q mismatching x by at most 2 digits
(6)modifiedmismatchlist← the numbers in mismatchlist in base 10
(7)mut← list of each hashtable[i] for all i∈modifiedmismatchlist
(8) goo2(q,l,d,posi[i],mut)
(9)end for


## 3. Implementation

We implement these algorithms using C on a machine with 12 Intel Core i7 CPU 80 @ 3.33 GHz processors and 12 GB of memory. The datasets we use in this implementation are Barley ESTs taken from the genetic software HarvEST by Steve Wanamaker and Timothy Close of the University of California, Riverside (http://harvest.ucr.edu/). We use two different EST databases, one with 78 ESTs and another with 2838. In our experiments we search for oligos of lengths 27 and 28 since they are common lengths for oligonucleotides. As we increase the size of the database, we see that Algorithm 3 is the most efficient as anticipated (data shown in Tables 1 and 2).Table 1
Results of serial algorithms.
Algorithm
l
d
q
Dataset
Time taken (secs)
Non-unique oligos
Algorithm2
28
6
4
1 (78 ESTs)
163
46,469
Algorithm1
28
6
7
1 (78 ESTs)
131
46,469
Algorithm3
27
6
9
1 (78 ESTs)
231
46,564
Algorithm2
28
6
4
2 (2838 ESTs)
197, 500
1,611,241
Algorithm1
28
6
7
2 (2838 ESTs)
117, 714
1,611,241
Algorithm3
27
6
9
2 (2838 ESTs)
94, 317
1,614,235Table 2
Results of parallel algorithms on 12 processors.
Algorithm
l
d
q
Dataset
Time taken (secs)
Non-unique oligos
Algorithm2
28
6
4
1 (78 ESTs)
33
46,469
Algorithm1
28
6
7
1 (78 ESTs)
29
46,469
Algorithm3
27
6
9
1 (78 ESTs)
66
46,564
Algorithm2
28
6
4
2 (2838 ESTs)
40, 420
1,611,241
Algorithm1
28
6
7
2 (2838 ESTs)
22, 848
1,611,241
Algorithm1
27
6
9
2 (2838 ESTs)
18, 375
1,614,235One important thing to note about all of these algorithms is the fact that the main portion of them is a for loop which iterates through each index of the hash table. It is also obvious that loop iterations are independent of each other. These two factors make the algorithms perfect candidates for parallelism. Rather than process the hash table one index at a time, our parallel algorithms process groups of indices simultaneously. Ignoring the communication between processors, our algorithms optimally parallelize our three serial algorithms.There are many APIs in different programming languages that aid in the task of parallel programming. Some examples of this in the C programming language are OpenMP and POSIX Pthreads. OpenMP allows one to easily parallelize a C program amongst multiple cores of a multicore machine [8]. OpenMP also has an extension called Cluster OpenMP which allows one to parallelize across multiple machines in a computing cluster.A new trend in parallel programming is in the use of GPUs. GPUs are the processing units inside computers graphics card. C has several APIs which allow one to carry out GPU programming. The two such APIs are OpenCL and CUDA [9, 10].In the second implementation of our algorithms we use OpenMP to parallelize our algorithms throughout the 12 cores of our machine. We can easily see that we achieve near optimal parallelization with our parallel algorithms; that is, the time taken by the parallel algorithms is approximately that of the serial algorithms divided by the number of processors.


## 4. Conclusion

In this paper we used three algorithms to solve the unique oligos search problem which are extensions of the algorithm presented in [2]. We observed that we can achieve a significant performance improvement by parallelizing our algorithms. We can also see that Algorithm 3 yields the best results for larger databases. For smaller databases, however, the time difference between each pair of algorithms is negligible, but results in Algorithm 3 being the slowest, and this is due to the time required to compute the mismatches of each q-mer. Other algorithms can be obtained by setting q to different values. See Algorithms 1, 2, 3, 4, 5, 6, 7, and 8.Algorithm 4:Findqmers (q).
Require: EST database D={x1,x2,…,xk}, integer q
Ensure: A hashtable of all qmer positions.
(1)hashtable← a hashtable of all qmer positions in D
(2)for  i←1 to k  do
(3)for  j←1 to length(D[i])-q+1  do
(4)hashedqmer← map(D[i],j,j+q-1)
(5)hashtable[hashedqmer] ← Append(hashtable[hashedqmer], {i,j})
(6)end for
(7)end forAlgorithm 5:Map (string s,i,j).
(1)r← substring(s, i, j)
(2)t←r under the transformation {A,C,G,T}→{0,1,2,3}
(3)return  tAlgorithm 6:Substring (string s,i,j).
(1)r← substring of s from character i to character j
(2)return  rAlgorithm 7:goo2(q,l,d,posi,mut).
(1)posi← a list of positions of a specified qmer in D
(posi={{x1,y1},{x2,y2},…} where {x,y} corresponds to position y of sequence x)
(2)mut← a list of positions of qmers in D that mismatch this qmer by either 1 or 2 characters
(depending on the filtration algorithm using this function)
(3)for  i←1 to length(posi)  do
(4)for  j←i+1 to length(posi)  do
(5)if  posi[i][1]≠posi[j][1]  then
(6)lq1← list of l-mers generated from the extension of the qmer in position posi[i]
(7)lq2← list of l-mers generated from the extension of the qmer in position posi[j]
(8)for  x←1 to length(lq1)  do
(9)for  y←1 to length(lq2)  do
(10)if  HD(lq1[x],lq2[y])≤d  then
(11)       mark thelmers as non-unique
(12)end if
(13)end for
(14)end for
(15)end if
(16)end for
(17)for  k←1 to length(mut)  do
(18)if  posi[i][1]≠mut[k][1]  then
(19)lq1← list of l-mers generated from the extension of the qmer in position posi[i]
(20)lq2← list of l-mers generated from the extension of the qmer in position mut[k]
(21)for  x←1 to length(lq1)  do
(22)for  y←1 to length(lq2)  do
(23)if  HD(lq1[x],lq2[y])≤d  then
(24)       mark thelmers as non-unique
(25)end if
(26)end for
(27)end for
(28)end if
(29)end for
(30)end forAlgorithm 8:goo(q,l,d,posi).
(1)posi← a list of positions of qmer in D
(posi={{x1,y1},{x2,y2},…} where {x,y} corresponds to position y of sequence x)
(2)for  i←1 to length(posi)  do
(3)for  j←i+1 to length(posi)  do
(4)if  posi[i][1]≠posi[j][1]  then
(5)lq1← list of l-mers generated from the extension of the qmer in position posi[i]
(6)lq2← list of l-mers generated from the extension of the qmer in position posi[j]
(7)for  x←1 to length(lq1)  do
(8)for  y←1 to length(lq2)  do
(9)if  HD(lq1[x],lq2[y])≤d  then
(10)      mark thelmers as non-unique
(11)end if
(12)end for
(13)end for
(14)end if
(15)end for
(16)end for



---

*Source: 793130-2013-04-08.xml*