in the bioinformatics and computational biology ( bcb ) domain , biological sequence alignment is a quite common task .
it includes comparison of dna ( nucleotide ) , rna ( nucleotide ) , and protein ( amino acid ) sequences . in such comparison
, a query sequence is aligned to subject sequences from a large database to find similarities .
scientists are used to adopting sequence alignment to infer biological information of a newly discovered sequence from a set of previously known sequences . for instance , if a recently discovered sequence is similar to a known disease gene , then the biological information about the functionality of the new sequence can be inferred .
in addition , biological sequence alignment plays an important role in the study of evolutionary development and the history of species .  
being one of the most important tools applied in finding biological sequence alignment , blast ( basic local alignment search tool )   has been widely implemented on commodity pc clusters .
but with the exponential growth of the biosequence databases ( see figure 1 ) , meeting the computational requirements using current platforms is becoming a difficult task .  
compared to general - purpose pc microprocessors , field programmable gate arrays ( fpgas ) are able to provide a much higher degree of bit - level data parallelism which leads to higher computational speed .
scientists have been investigating parallel architecture design for the blast algorithm on fpgas to achieve higher computational efficiency .  
this paper proposes a design and implementation of a parallel architecture on fpga for blast , namely , systolic array - based parallel architecture for the blast algorithm .
the design is implemented with very high speed integrated circuit hardware description language ( vhdl ) , making it compatible across a number of fpga architectures ( e.g. , xilinx , altera ) .
we first present the overview of the blast algorithm as well as the related work .
after that , architecture performance including storage resource requirements , fpga resource consumption , and estimated speedups against other architectures will be analyzed based on synthesizing , functional , and timing simulation .
although heuristic algorithms produce local alignments that are not always optimal , they are normally much quicker than exhaustive dynamic programming algorithms .
the execution time of biological sequence alignment is crucial to bioinformatics scientists , especially due to the fact that the size of the biosequence databases has grown exponentially over years .
the three main steps of the blast algorithm are finding hits , ungapped extension , and gapped extension .
( 1 ) finding hitsin this step , we use a protein sequence which consists of eight residues ( amino acids ) as an example shown below :  pqgefgvy . 
  in this step ,
we use a protein sequence which consists of eight residues ( amino acids ) as an example shown below :  pqgefgvy .
the query sequence is divided into overlapping k - words as illustrated below with k = 3 :  word 0 : pqg  word 1 : qge  word 2 : gef  word 3 : efg  word 4 : fgv  word 5 : gvy.six words are extracted from this sequence . generally speaking , ( m  k ) +
1 words would be extracted when the size of the sequence is m. after dividing the query sequence , these k - words plus a score matrix are used to find the possible matching words which have a score with at least the threshold value of t. finally , subject sequences in the database are scanned one by one to find exact matches to those possible matching words .
previous work shows that this step could be finished by high - level application software running on a pc . by designing an architecture which can implement every step of blast on a fpga , the interface issue for software and hardware connection
can be avoided , saving execution time . instead of finding exact matches to those possible matching words , in this paper
, the implementation searches for the exact matches of those k - words directly on the fpga .
( 2 ) ungapped extensionin this step , each hit recorded in the last step is extended by aligning with a subject sequence in both directions without any gap .
a score matrix , for example , blosum50 , is used to compute residue pair alignment scores .
the extension process is terminated when the total score of the alignment is below a threshold value .
these results from ungapped extension are called high - scoring segment pairs ( hsps ) . 
  in this step , each hit recorded in the last step is extended by aligning with a subject sequence in both directions without any gap .
a score matrix , for example , blosum50 , is used to compute residue pair alignment scores .
the extension process is terminated when the total score of the alignment is below a threshold value .
these results from ungapped extension are called high - scoring segment pairs ( hsps ) . 

( 3 ) gapped extension in this stage , the smith - waterman or needleman - wunsh algorithm is applied to perform the gapped alignment in both directions for those collected hsps . after this step ,
, the smith - waterman or needleman - wunsh algorithm is applied to perform the gapped alignment in both directions for those collected hsps . after this step ,
we choose the needleman - wunsch algorithm for gapped extension implementation because of its ability to find optimal global gapped alignment between the query sequence and the subject sequence .
suppose that there are two sequences , q and s , whose lengths are m and n , respectively .
q[1 : i ] denotes the first i characters in sequence q and s[1 : j ] denotes the first j characters in sequence s. the similarity score between q[1 : i ] and s[1 : j ] can be obtained from the similarity scores computed for q[1 : i  1 ] and s[1 : j  1 ] .
a dynamic programming score matrix d is built , where each element of the matrix denotes the best alignment score between the q[1 : i ] segment of q and the s[1 : j ] segment of s.   the first step of needleman - wunsch algorithm is to create the score matrix d with m + 1 columns and n + 1 rows .
m and n are lengths of sequences q and s. the first row and the first column are assigned the value zero .
the computation complexity of the initialization step is o(m + n ) .   the second step is comprised of filling the matrix .
di , j is used to express the cell ( i , j ) of the scoring matrix d. di , j is defined as the maximum of the similarity scores computed in positions ( i  1 , j  1 ) , ( j , j  1 ) , and ( i  1 , j ) associated with gap penalty or si , j .
the following equation is used to compute the value of the cell ( i , j ) : 
 ( 1)d(i , j)=max{d(i1,j1)+si , j , d(i1,i)+d , d(i , j1)+d . in this equation , si
si , j = 1 if the ith character in q is equal to the jth character in s ( match score ) ; otherwise si , j = 0 ( mismatch penalty ) .
d(i  1 , j ) + 1 is the score of an alignment between xi and a gap in y. d(i , j  1 ) + d is the score of an alignment between a gap in x and yj .
trace back can form the best global alignment between q and s. it starts from the bottom - right value in the matrix and ends at the up - left value . for each step ,
the next alignment pair which is added to the front of the alignment segment accords to 3 rules : ( 1 ) if the cell ( i , j ) was derived from cell ( i  1 , j  1 ) , the pair of alignment xi and yj is added to the front of the alignment segment .
( 2 ) if cell ( i , j ) was derived from cell ( i  1 , j ) , xi and a gap in y are added to the front of the alignment segment .
( 3 ) if cell ( i , j ) was derived from cell ( i , j  1 ) , a gap in x and yj are added to the front of the alignment segment . in some cases
it 's hard to decide which cell cell ( i , j ) comes from .
the diagonal cell is chosen under this situation to shorten the length from the bottom - right cell to the upper - left cell .
it corresponds to the following alignment : t_g g c c a g t c c g and t t g_c_a_t_g .
fpgas not only provide high speedups and high bit - level data parallelism compared to regular processors , but also include the convenience of reprogrammability . more and more bioinformatics scientists have designed parallel fpga architectures for the blast algorithm in order to solve the biological sequence alignment issue .
many architectures have been proposed already such as mercury blastn [ 1214 ] , tree - blast , mercury blastp , rc - blast , fpga / flash accelerator , multiengine blastn accelerator [ 18 , 19 ] , and many industry applied systems like bee2   and mitrion .
most of them adopt the word position record - based search ( wprbs ) method , which constructs all storage tables to record the position of each word in query sequence , then sends those words to the accelerator one by one where the storage table is searched to find the hit .
although this approach is widely used , it still suffers from time and storage limitations . using wprbs , at most one hit can be found in one - clock cycle because only one word can be searched per clock cycle .
fpga memory port number also limits wprbs searching capacity no matter where the storage tables are stored : internal chip memory or external chip memory .
take storing tables in an external memory for instance , mercury blastn   and mitrion   built a hash table from the query . during the searching process , an accelerator checks words in the subject sequence database against the hash table one by one .
an external sram attached to the fpga is used to store the hash table since the internal block rams do not have enough storage space to hold tables for long query sequences .
this situation obviously causes a timing issue : accessing delay to external sram leads to long pipeline cycle time .
the database is formatted as an index structure in which every word is associated with its position in the sequence and its adjacent environment so that short , ungapped alignments could be computed simultaneously , avoiding many random accesses to the database . but the size of the database index is very large , a 150  gb index for the human genome is needed when storing a 40 amino acid substring environment .
this requires 50 times more storage than the raw data . in the future , with the exponential growth of the database , the storage resource requirement will become a serious issue . to get better searching efficiency , multiengines blastn   adopted 64 identical computing units in a single chip to compare the query with 64 subject sequence in database at the same time .
mercury blastp   proposes a two - hit generator for accelerating the first stage of blastp .
one is used as a subject sequence database and the other is used as a results container which receives hsps lists from the basic layer .
this layer scans all subject sequences from database for the query sequence in order to get a hsps list . after that
, the list is sent to the host computer software that assigns statistical significance specified by bioinformatics scientists . in the top layer of this system , a dell precision t1500 desktop computer with an intel core i7 microprocessor and 4  gb memory runs the host software . in the middle layer , two 1  gb
the basic layer is the systolic array - based parallel architecture that implements the blast algorithm .
this layer is mapped to one fpga chip ( xilinx virtex-5 xc5vlx110 t ) which is embedded on the xilinx ml509 evaluation platform .
the initialization set and data control block is used to initialize the status of registers for the hit finder array and to forward subject sequence as well as query sequence to the hit finder array .
for example , it can find adjoined overlapping hits  tkel  and  kelp  then combine them into  tkelp  .
ungapped extension blocks extends ungapped hits in both directions until the score meets a thread value t. gapped extension block extension gapped extension on hsps through exhaustive algorithm like needleman - wunsch .
multiple hits detection module is used to detect 3-word hits and record the hits ' address in the query and the subject sequence . compared to wprbs method which could detect at most one hit in one clock cycle ,
as the figure illustrated , there is a systolic array with 32 processing units , every 3 units are connected to one 3-input and gates .
the value of each register is sent to the corresponding hits information extraction units for recording the hits address in the query and the subject sequence . the systolic array and the hits information extraction unit
clock_en ( array_clk in figure 4 ) drives hits information extraction units ; clock_out is the clock source of the systolic array .
these two clocks work together driving a pipelined architecture without any pause . at each clock_out rising edge , query or subject sequence moves forward for one processing unit ;
at each clock_en rising edge , a bit value at a certain position in each register is used to detect the location of hits .
the whole architecture works as follows : first , a query sequence with 32 characters is forwarded into the systolic array so that each processing unit holds a character from the query sequence .
then the subject sequence is driven in to the systolic array by each internal clock rising edge .
meanwhile , the incoming subject character and the query character which are held by the unit are compared , if they are identity , the logic 
1  would be generated ; otherwise , the logic  0  would be generated .
so , the systolic array with 3-input and gates can detect multiple hits at one internal clock rising edge . the architecture of the processing unit in the systolic array is illustrated in figure 6 . in the implementation , 20 3-bit long binary numbers (  00001  to  10100  ) are defined to present 20 amino acids in a protein .
shift registers are able to shift characters to the adjacent processing unit so that the query and the subject sequence can go through the systolic array .
when a new character from the query or the subject sequence is shifted into the unit , the corresponding i d register would increase its value by 1 .
each processing unit of the systolic array will compare two characters in it at each internal clock rising edge .  
as shown in figure 4 , outputs of 32 3-input and gates goes into 2 16-bit registers .
the systolic array with 32 processing units cooperates with 32 3-inputs and gates to detect hits in both sequences .
if there is high similarity between query and subject sequence , the multiple hits detection module may output a large amount of hits per clock cycle .
ungapped extension block is relatively slower than multiple hits detection module and the systolic array in this situation . as a result
previous research [ 23 , 24 ] suggested that in many cases there is a high similarity between the query sequence and the subject sequence .
processing those hits not only wastes the execution cycle for multiple hits detection module , but also consumes additional extension time for ungapped extension block because one hit is extended twice .
for instance , two adjacent hits  atk  and  tkp  are found but they are actually one hit 
if they are not combined into one hit , they would have been recorded twice in multiple hits detection module and extended twice in ungapped extension block . in this section , we implement hits combination block to address this issue
. it can detect the overlapping hits and merge them to reduce verbose hits and maintain the sensitivity of blast .  
this block contains a hits ( first in first out ) fifo buffer which is used to store hits location addresses from both query and subject sequence .
the algorithm implementation engine executes hits combination algorithm then forward the hit addresses and hit length to the ungapped extension block .
the order of the hit addresses in the hits fifo buffer follows the direction of the stream flow ( figure 5 ) . for example , the address of  akl  should be stored in the fifo buffers first , then the  klp  .
the combination algorithm is triggered when there is more than one record in the hits fifo buffer . here
, the hit addresses in subject sequence are recorded and forwarded to the following operations .
first , the algorithm implementation engine calculates the ending address of the ith and ( i + 1)th hit in subject sequence . in this paper ,
the starting address is the recorded address minus 1 , because the recorded address is the address of the middle word in a 3-word hit .
the ending address is the address of the last character in the sequence , so ( ending  address = starting  address + ( seqence  length  1 ) ) . here , the sequence length is equal to 3
. then , a comparison between ith hit 's ending address and ( i + 1)th hit 's starting address is performed . if ith hit 's ending address is greater than or equal to the ( i + 1)th hit 's starting address , there would be an overlap between two hits .
the algorithm can determine new addresses in both subject sequence and query sequence of the merged hit and calculate its length .
in the next step , the algorithm records the new address and hit length at the ith place in hit 's fifo buffer .
if ith hit 's ending address is less than the ( i + 1)th hit 's starting address , the ( i + 1)th hit address will be chosen to be compared with ( i + 1)th hit address .
this operation continues until the hits fifo buffer is empty . because multiple hits detection module can detect multiple hits in one internal clock cycle , multiple modules working together in parallel speeds - up hit detection . in order to construct a longer array ,
as figure 8 shows , the large systolic array contains many multiple hits detection modules .
every module matches one hits combination block , each of which records and combines hits detected by its mapped module .
then fifo combination modules combine hits fifos and deliver hits information to larger hits fifos . in this architecture ,
multiple hits detection modules are connected in a serial way to extend the array comparison size .
each hits combination block maps to a multiple hits detection module so that each block can detect overlapping hits and merge them simultaneously .
the score matrix blosum50 [ 10 , 25 ] is applied to find the alignment score for a pair of words .
if this block detects a hit in a fifo , the address of the hit word in query sequence , and subject sequence and the hit length would be extracted from the fifo to compute both extension points .
after this step , the ungapped extension algorithm implementation engine can finish the ungapped extension step .
the data flow in the ungapped extension block is shown in figure 9 . the extension procedure in algorithm implementation engine
after computing the addresses of characters  k  and  p ,  the starting addresses of the extension , the algorithm drives the extension in both directions .
each pair of the residues own one alignment score which can be found in the blosum50 score matrix . in figure 13 , the alignment score 
l  and  d  is 3 , and the score of the pair 
then endpoints of this high - scoring segment in both directions on both query and subject sequences are recorded and forwarded to the gapped extension block for gapped alignment . in this step ,
hits which are sent from hits fifo of the hits combination block are extended to either side to identify a hsp .
mercury blastn prefilters the hit 's characters because only one hit can be searched at one clock cycle .
fpga / flash could get the subsequence directly because it constructs the index for each word .
compared to methods applied in mercury blastn and fpga / flash , multiple hits detection module is more efficient since it could detect multiple hits in one clock cycle .
however , the serial ungapped extension process is slower than multiple hits detection module . to speedup this step , we proposed two methods in the above sections : adding hits combination block between multiple hits detection module and ungapped extension block ; dividing the long components array into several parallel components array groups .
as shown in figure 11 , we propose one more parallel architecture for ungapped extension block . in this architecture , each hits fifo of hits combination block connects to an ungapped extension block .
each block connects to the query sequence memory as well as initialization set & data control block because the context of the hits is necessary for extension .
system may slow down if hits from all hits fifos are queued in only one ungapped extension block . in this architecture , hits information could be distributed into different paths and save processing time . in this block , hsps from the fifo buffer are used to perform gapped extension .
the gapped algorithm implementation engine will be triggered once the block detects that there is a pair of segments in fifo .
this paper focuses on the first two steps of blast algorithm including finding hits and ungapped extension because they consume most execution time in blast algorithm . in order to implement the needleman - wunsch algorithm on the fpga
we have mapped multiple hits detection module , hits combination block , and ungapped extension block architecture on the xilinx xc5vlx110 t fpga which is embedded on ml509 education board .
the architecture was designed under xilinx ise design suite 12.1 in very high - speed integrated circuit hardware description language ( vhdl ) .  
in our system , the xilinx ml509 education board has a pci express(x1 ) interface which could provide 500  mb / s throughput .
the board is then connected to the host computer : dell precision t1500 pc with a 2.80  ghz intel core i7 microprocessor and 4  gb memory .
ncbi blast software runs on the host . in our experiment , a systolic array - based parallel architecture with 2048 processing units was implemented on the xilinx ml509 . except for this implementation
, we also mapped systolic array based architectures with 1204 , 2048 , and 3072 processing units to different fpga devices under xilinx ise design suite .
the clock values came from the timing analyze tool called postplace and route which integrated in the ise .
we will analyze the architecture performance from two aspects : one is comparing our architecture with previous architectures on hardware , the other is comparing our architecture against blast software versions . for hardware comparison
, we focus on 3 parts including synthesis performance , storage requirement , and word scanning speed . for software versions comparison
, a hardware performance comparison is made between ours and other fpga architectures from 3 aspects including synthesis performance , storage requirement , and word scanning speed .
tree - blast and mercury blastn are wprbs architectures , families of fpga - based accelerators is a systolic array - based architecture .
all data which reflects the capability of architectures in the comparison from previous publication .  
because every four components need a bram to index the scoring matrix , the input array size the system could process is restricted by bram . as a result , the query size is limited up to 600 on xc2vp70 and 1024 on xc4vlx160 .
xia et al .   implemented a fpga - based accelerator which is also a systolic array - based one for blast algorithm called families of fpga - based accelerators .
in addition , the data stream in the array needs to be shut down to record hits addresses . as a result , it needs more registers to store the systolic array status in order to resume it . in our architecture ,
wprbs is replaced by multiple hits detection module to find hits so that the bram is not the bottleneck of array size anymore but luts are .
we can build 2048 components on xc5vlx110 t which is two times the array size of tree - blast . based on synthesis report
our design has longer array size and consumes less resource , comparing with tree - blast built on xc4vlx160 with 78% of the slices . for synthesis performance comparison , we also implement our architecture with 2048 components on fpga xc2vp70 and fpga xc4vlx160 .
we constructed two systolic arrays with 1024 and 3072 processing units to make a comparison in synthesis performance with the families of fpga - based accelerators .
the memory usage is 3408 kbits in total which is 20.3% of the memory capacity in xc5vlx110 t .
since many wprbs architectures need to record the position of each word in query sequence , they generally consume much more storage .
for instance , rc - blast spends 64  k x 64bits memory to store the query index in xilinx 4085xla .
mercury and mitrion have external sram to store the hash table because internal memory is not enough for the record .
compared to these methods , our architecture saves storage , avoids the complexity of memory access , and simplifies the fpga layout and routing .
the speed of multiengine blastn accelerator which applies 64 identical parallel engines reaches 6400  m matches / s . the word scanning speed of our multiple hits detection module approach achieves 14450  m matches / s which is over 2 times than related work .  
our design is mapped to the xilinx xc5vlx110 t fpga to accelerate the first two steps of ncbi blast family programs .
as table 3 illustrated , swiss - prot from ebi   and drosoph.net from ncbi blast are two databases in experiments .
swiss - prot is a protein sequence database while drosoph.net is a database for dna sequences .
the swiss - prot in our experiments is uniprotkb / swiss - prot protein knowledgebase release 2012_11 .
queries which are picked from drosoph.nt translated into six - frame protein sequence before searching against swiss - prot .
tblastx is sued for searching six - frame translation of dna sequence against six - frame translation of dna database . for each algorithm
, we built 6 arrays with different sizes ( 128 - 3k ) to carry different queries with different lengths .
for each length , 15 queries with the same length are picked up randomly from the relevant database .
the final execution time for a certain array is the average of running time of 15 queries .  
the execution time of the software versions in this table comes from the previous work .
the hardware execution time ( ht ) is less than software execution time ( st ) for all cases .
blastp software version only spends 1901  ms in searching database with a 128 long query while 25132  ms is taken when it is searching database with a 3-k long query .
index construction not only needs more time but also increases the size of searching object . in our architecture , however , the execution time does not increase so rapidly .
this is because the searching time is equal to the sum of the time of database stream flowing through the array and the extension time of valid hits .
the time of database stream flowing is determined by the size of database ; the extension time of valid hits is influenced by the architecture design . in this experiment ,
the hits combination block and the parallel design reduce the number of valid hits and their extension time .
for the same reason , the execution time of tblastn , blastx , and tblastx which are implemented on the fpga does not increase sharply . from table 3
, it can be seen that tblastn software version is slower than blastp software version for the same query length .
the reason is that all dna sequences are translated into the 6 possible potential proteins in tblastn software version before the searching .
the translation time is not counted for the hardware version in our experiments because the translation is a one - time job , and results can be used for many other applications . for tblastx software version ,
both queries and the database are translated into the six possible potential proteins before searching . the translation time is not counted for the same reason as in tblastn .
our architecture could implement 3 steps of blast algorithm simultaneously because the fpga is a good platform to parallel the blast algorithm . paralleling blast algorithm on the fpga
patternhunter has 48  m/4s word scan speed at most while ours could get 14450  m matches / s .
in this paper , we presented a systolic array - based fpga parallel architecture for blast algorithm .
the architecture contains four different kinds of blocks which are multiple hits detection module , hits combination block , ungapped extension block , and gapped extension block .
we designed the multiple hits detection module to improve the efficiency of the bottleneck step .
the multiple hits detection module based on the systolic array could find more than one hit in only one clock cycle .
the performance is faster than wprbs searching approach used in other proposed architectures cited in this paper . in order to speedup hits searching , hits
combination block between multiple hits detection module and ungapped extension block merges adjacent overlapping hits into one hit in order to speedup the extension step .
lastly , we proposed a parallel architecture for ungapped extension block as figure 12 illustrated .
we have implemented multiple hits detection module , hits combination block as well as ungapped extension block .
based on the ise synthesis report , we compared our architecture performance to the performance of other related work .
first , compared to tree - blast , a twice longer array can be built in our architecture with less fpga resource .
thirdly , the word scanning speed of our architecture faster than mercury blastn 's word scanning speed .
last , compared to blastp software , our architecture is more suitable for dealing with longer sequences . in the future ,
we hope to design an efficient architecture for gapped extension block . to satisfy more users ' requirement
the multiple hits detection module needs to be changed , for instance , 3-input and gate should be changed to 11-input and gate .
the clock_division block needs to be redesigned to guarantee the pipeline works correctly , more details will be addressed in the further design .
our architecture is only available for users who use xilinx xc5vlx110 t fpga because we only test it on that fpga to guarantee the correctness of its behavior and timing . after enough tests on different devices ,
 a design of systolic array - based field programmable gate array ( fpga ) parallel architecture for basic local alignment search tool ( blast ) algorithm is proposed . 
 blast is a heuristic biological sequence alignment algorithm which has been used by bioinformatics experts . 
 in contrast to other designs that detect at most one hit in one - clock - cycle , our design applies a multiple hits detection module which is a pipelining systolic array to search multiple hits in a single - clock - cycle . 
 further , we designed a hits combination block which combines overlapping hits from systolic array into one hit . 
 these implementations completed the first and second step of blast architecture and achieved significant speedup comparing with previously published architectures . 