Database-Indexed Protein Sequence Search on Multicore CPUs

14 Pages·2016·1.57 MB·English

by  Jing Zhang

Zhangetal.BMCBioinformatics (2016) 17:443
DOI10.1186/s12859-016-1302-4
SOFTWARE OpenAccess
muBLASTP: database-indexed protein
sequence search on multicore CPUs
JingZhang1*,SanchitMisra2,HaoWang1andWu-chunFeng1
Abstract
Background: The BasicLocalAlignmentSearchTool(BLAST)isafundamentalprograminthelifesciencesthat
searchesdatabasesforsequencesthataremostsimilartoaquerysequence.Currently,theBLASTalgorithmutilizesa
query-indexedapproach.Althoughmanyapproachessuggestthatsequencesearchwithadatabaseindexcan
achievemuchhigherthroughput(e.g.,BLAT,SSAHA,andCAFE),theycannotdeliverthesamelevelofsensitivityasthe
query-indexedBLAST,i.e.,NCBIBLAST,ortheycanonlysupportnucleotidesequencesearch,e.g.,MegaBLAST.Dueto
differentchallengesandcharacteristicsbetweenqueryindexinganddatabaseindexing,theexistingtechniquesfor
query-indexedsearchcannotbeusedintodatabaseindexedsearch.
Results: muBLASTP,anoveldatabase-indexedBLASTforproteinsequencesearch,deliversidenticalhitsreturnedto
NCBIBLAST.OnIntelHaswellmulticoreCPUs,forasinglequery,thesingle-threadedmuBLASTPachievesuptoa
4.41-foldspeedupforalignmentstages,anduptoa1.75-foldend-to-endspeedupoversingle-threadedNCBIBLAST.
Forabatchofqueries,themultithreadedmuBLASTPachievesuptoa5.7-foldspeedupsforalignmentstages,andup
toa4.56-foldend-to-endspeedupovermultithreadedNCBIBLAST.
Conclusions: WithanewlydesignedindexstructureforproteindatabaseandassociatedoptimizationsinBLASTP
algorithm,were-factoredBLASTPalgorithmformodernmulticoreprocessorsthatachievesmuchhigherthroughput
withacceptablememoryfootprintforthedatabaseindex.
Keywords: BLAST,Databaseindex,Localalignment,Multicore
Background information requires substantial computing power. Con-
The Basic Local Alignment Search Tool (BLAST) [1] is sequently, significant research effort has been invested
afundamentalalgorithminlifesciencesthatcomparesa intoacceleratingtheBLASTsearchalgorithm.
query sequence to a database of sequences, i.e., subject Much of this research effort has focused on the paral-
sequences, to identify sequences that are the most sim- lelizationofBLASTondifferentparallelarchitecturesdue
ilar to the query sequence. The similarities identified by toitscompute-anddata-intensivenature.NCBIBLAST+
BLASTcanbeusedtoinferfunctionalandstructuralrela- [2] uses pthreads to speed up BLAST on a multicore
tionships between the corresponding biological entities, CPU. On CPU clusters, TurboBLAST [3], ScalaBLAST
forexample. [4], and mpiBLAST [5] have been proposed. To achieve
Withtheadventofnext-generationsequencing(NGS), higher throughput on a per-node basis, BLAST has also
whether at the outset or downstream from NGS, the been mapped and optimized onto various accelerators,
exponentialgrowthofsequencedatabasesisarguablyout- includingFPGAs[6,7]andGPUs[8–13].However,there
stripping our ability to analyze the data. Specifically, the are few recent studies that focus on improving the per-
increasingdemandstominesequencedatabasesforuseful formance of CPU implementations of the widely-used
BLASTalgorithm.
Mostpreviousstudies[1,4,14,15]adoptqueryindexing
*Correspondence:[email protected] for sequence search. Query indexing uses a lookup table
1DepartmentofComputerScience,VirginiaTech,225StangerStreet,24060 torecordpositionsofeachwordintheinputquery.These
Blacksburg,VA,USA
BLAST algorithms then scan each database sequence to
Fulllistofauthorinformationisavailableattheendofthearticle
©TheAuthor(s).2016OpenAccessThisarticleisdistributedunderthetermsoftheCreativeCommonsAttribution4.0
InternationalLicense(http://creativecommons.org/licenses/by/4.0/),whichpermitsunrestricteduse,distribution,and
reproductioninanymedium,providedyougiveappropriatecredittotheoriginalauthor(s)andthesource,providealinktothe
CreativeCommonslicense,andindicateifchangesweremade.TheCreativeCommonsPublicDomainDedicationwaiver
(http://creativecommons.org/publicdomain/zero/1.0/)appliestothedatamadeavailableinthisarticle,unlessotherwisestated.
Zhangetal.BMCBioinformatics (2016) 17:443 Page2of14
findshortmatches,extendthesematchestooptimalalign- two-hit ungapped extension in protein sequence search.
ments, and then calculate the final similarity scores. In Thus,databaseindexingthatalignsaquerytomillionsof
contrast,otherapproachessuggestthatdatabaseindexing databasesequencesinsteadofasingledatabasesequence
canyieldmuchfasterspeedthanqueryindexing[16,17]. iterativelywillsuffermorefromsuchirregularities,lead-
Examples of such tools include BLAT [1], SSAHA [18], ingtoseriousperformancedegradation.
MegaBLAST [19], and CAFE [20]. However, these tools To overcome these challenges of database indexing
cannotprovidethesamelevelofsensitivityastheBLAST for protein sequence search, we propose muBLASTP
algorithm [17, 21, 22], or support nucleotide sequence (i.e., microprocessor-based BLASTP), a novel BLASTP
search. algorithm that includes an advanced index data struc-
SSAHAandBLAT,forexample,aresignificantlyfastfor ture for sequences of the database and a set of opti-
findingnear-identicalmatches.However,toreducemem- mizations for the BLASTP algorithm. The experimental
ory footprint and search space, both tools build indexes results show that on a modern multicore architecture,
ofnon-overlappingwordsfromthedatabase,whichleads namely Intel Haswell, for a single query, the single-
to extremely fast search but compromised sensitivity. threaded muBLASTP can deliver up to a 4.41-fold
More specifically, BLAT, for example, builds database speedupforalignmentstages,anduptoa1.75-foldend-
index with non-overlapping words of length W. With to-end speedup over the single-threaded NCBI BLAST.
this approach, the size of database index is significantly For a batch of queries, the multithreaded muBLASTP
reduced, roughly 1 the size of an index with overlap- can achieve up to a 5.7-fold speedup for alignment
W
ping words. However, it requires a matching region of stages, and 4.56-fold end-to-end speedup over the mul-
2W − 1 letters between two sequences for guarantee- tithreaded NCBI BLAST using 24 threads. The exper-
ing to detect it. CAFE is another search tool supporting imental results also shows that on a older generation
protein sequence with database index, but the search multicore architecture, namely Intel Nehalem, for a sin-
method and scoring phase are substantially changed. gle query, muBLASTP still can deliver up to a 3.8-fold
MegaBLASTistheonlyBLASTvariantbasedondatabase speedupforalignmentstages,anduptoa1.94-foldend-
index.MegaBLASTacceleratesthesearchforhighlysim- to-end speedup over the single-threaded NCBI BLAST.
ilar sequences by using a large word size (W = 28) to For a batch of queries, the multithreaded muBLASTP
reducethesearchworkloadandthememoryusage.How- can achieve up to a 8.59-fold speedup for alignment
ever, according to the previous studies [23–25], increas- stages, and 3.85-fold end-to-end speedups over the mul-
ing word size could sacrifice the sensitivity and accu- tithreaded NCBI BLAST using 12 threads. In addition
racy.Furthermore,MegaBLASTonlysupportsnucleotide toimprovingperformancesignificantly,muBLASTPpro-
sequences,astheauthorsclaimedthatitisverychalleng- duces identical hit returned to NCBI BLAST, which is
ingtosupportproteinsequencebasedontheirdesign. importanttothebioinformaticscommunity.
Becausequeryindexingusuallycontainsahighpercent-
age of empty slots due to few letters in a query, most of Implementation
the optimizations of query indexing seek to reduce the
Databaseindex
sparsity of the index, e.g., the thick backbone and the
The most challenging component of muBLASTP is the
positionarrayinNCBIBLAST[26]andthedeterministic
designofthedatabaseindex.Theindexshouldincludethe
finiteautomaton(DFA)inFSA-BLAST[14].Fordatabase
positionsofoverlappingwordsfromallsubjectsequences
indexing,whichisfullofpositionsfrommillionsofsubject
of the database. Thus, each position contains the infor-
sequencesfromadatabase(e.g.,about6millionsequences
mation for the sequence id and the offset in the subject
in env_nr database, and over 85 million sequences in nr
sequence, i.e., subject offset. For the protein sequence
database), the major challenges differ substantially from
search, the BLASTP algorithm uses the small word size
query indexing. First, the size of the database index can (W = 3),largealphabetsize(22letters),andneighboring
be prohibitive, especially for the protein database, which
word comparisons. Because these factors may make the
hastheincreasedalphabetandtheshortwordlength.Sec-
database index very large, we design our database index
ond,unlikenucleotidesequencesearch,proteinsequence
with the following techniques: index blocking, sorting,
search needs to search the hits of similar words, i.e.
andcompression.
the neighboring words rather than merely and exactly
matched words. Including neighboring words increases Indexblocking
the size of the index by one or two orders of magni- Figure 1a illustrates the design of index blocking. We
tude.Third,BLASTemploysinput-sensitiveheuristicsto first sort the database by the sequence length; partition
quickly eliminate unnecessary search spaces. However, thedatabaseintosmallblocks,whereeachblockhasthe
thisheuristicintroducessignificantirregularitiesinmem- samenumberofletters;andthenbuildtheindexforeach
oryaccesspatternsandincontrolflowpaths,e.g.during block separately. In this way, the search algorithm can
Zhangetal.BMCBioinformatics (2016) 17:443 Page3of14
Fig.1Anexampleofbuildingacompresseddatabaseindex.Thefigureshowstheflowfromtheoriginaldatabasetothecompressedindex.
aIndexblockingphasepartitionsthesorteddatabaseintoblocks.bBasicindexingphasegeneratesbasicindex,whichcontainspositionsofall
wordsinthedatabase.cIndexsortingsortspositionsofeachwordbysubjectoffsets.dIndexcompression-mergemergespositionswiththesame
subjectoffset.eIndexcompression-incrementdoneonthemergedpositionsgeneratesincrementsofsubjectoffsetsandsequenceids
go through the index blocks one by one and merge the fewerbitsforsubjectoffsetsintheblockshavingshortand
high-scoringresultsofeachblockinthefinalstage.Index mediumsequences,andmorebitsonlyfortheblockshav-
blocking can enable the database index to fit into main ingextremelylongsequences.(Thisisthereasonwhywe
memory,especiallyforlargedatabaseswhosetotalindex sortthedatabasebythesequencelength).
size exceeds the size of main memory. By configuring Furthermore,indexblockingallowsustoparallelizethe
thesizeoftheindexblock,wecanachievebetterperfor- BLASTP algorithm via the mapping of one block to a
mance.Forexample,iftheindexblockissmallenoughto thread on a modern multicore processor. For this block-
fitintotheCPUcache,thehitdetectionandgappedand wiseparallelmethodtoachievetheidealloadbalance,we
ungappedextensionmayachievebetterdatalocality. partitionindexblocksequallytomakeeachblockhavea
Anotherbenefitofusingindexblockingistoreducethe similarnumberofletters,insteadofanidenticalnumber
indexsize.Withoutindexblockingandassumingatotalof ofsequences.Toavoidcuttingasequenceinthemiddle,if
Msequencesinthedatabase,weneedlog Mbitstostore thissequencereachesthecapoftheblocksize,weputit
2
sequence ids. After dividing the database into N blocks, intothenextblock.
each block contains M sequences on average. Thus, we Afterthedatabaseispartitionedintoblocks,eachblock
N
only need log (cid:2)M(cid:3) bits to store sequence ids. For exam- isindexedindividually.AsshowninFig.1b,theindexcon-
2 N
ple, if there are 220 sequences in a database, we need 20 sistsoftwoparts:thelookuptableandthepositionarray.
bitstorepresentthesequenceids.With28blocks,ifeach Thelookuptablecontainsawentries,whereaisthealpha-
block contains 212 sequences, then we only need a max- betsizeofaminoacidsandwisthelengthofthewords.
imum of 12 bits to store the sequence ids. In addition, Each entry contains an offset to the starting position of
because the number of bits for storing subject offsets is thecorrespondingword.Inthepositionarray,aposition
determinedbythelongestsequencesineachblock,after ofthewordconsistsofthesequenceidandthesubjectoff-
sorting the database by the sequence length, we can use set. For protein sequence search, the BLASTP algorithm
Zhangetal.BMCBioinformatics (2016) 17:443 Page4of14
not only searches the hits of exactly matched words, but
it also searches the neighboring words, which are simi-
larwords.ThequeryindexusedinexistingBLASTtools,
e.g., NCBI BLAST, includes the positions of neighbor-
ingwordsinthelookuptable.However,forthedatabase
index in muBLASTP, if we store the positions for the
neighboring words, the total size of the index becomes
extraordinarilylarge. To addressthis problem, insteadof
storing positions of the neighboring words in the index,
weputtheoffsets,whichpointtotheneighboringwords
of every word, into the lookup table. The hit detection
stage then goes through the positions of neighbors via
theoffsetsaftervisitingthecurrentword.Inthisway,we
useadditionalstridememoryaccessestoreducethetotal
memoryfootprintfortheindex.
Indexcompression
As shown in Fig. 1b, a specific subject offset for a word
mayberepeatedinmultiplesequences.Forexample,the
word“ABC”appearsinposition0ofsequence1and3.In
lightofthisrepetition,itispossibletocompresstheindex
byoptimizingthestorageofsubjectoffsets.Next,wesort
thepositionarraybythesubjectoffsettogroupthesame Fig.2Anexampleofresolvingoverflowsinthecompressedindex.
subject offsets together, as shown in Fig. 1c. After that, aResolvingtheoverflowinthenumberofpositions.bResolvingin
theincrementalsubjectoffsets
we reduce the index size via merging the repeated sub-
jectoffsets:foreachword,westorethesubjectoffsetand
thenumberofpositionsonceandstorethecorresponding
the second has the incremental subject offset 0 for an
sequence ids sequentially, as shown in Fig. 1d. After the
additional45positions.
indexmerging,weonlyneedasmallarrayforthesorted
subject offsets. Furthermore, because the index is sorted
OptimizedBLASTPalgorithmwithdatabaseindex
bysubjectoffsets,insteadofstoringtheabsolutevalueof
BecausetheBLASTPsearchalgorithmintroducesamore
subjectoffsets,westoretheincrementalsubjectoffsets,as
irregular memory access pattern when using a database
notedinFig.1e,andonlyuseeight(8)bitsfortheincre-
index(ratherthanaqueryindex),weproposeandrealize
mental subject offsets. Because the number of positions
hitreorderingwithtwo-levelbinning inordertomitigate
for a specific subject offset in one block is generally less
theirregularmemoryaccesspatternandirregularcontrol
than256,wecanalsouseeight(8)bitsforthenumberof
flow,especiallyforthetwo-hitungappedextension.
positions.Thus,intotal,weonlyneeda16-bitintegerto
storeasubjectoffsetanditsnumberofpositions. Hitreorderingwithtwo-levelbinning
However,thiscompressedmethodpresentsachallenge. The two-hit ungapped extension in protein sequence
Whenweuseeight(8)bitseachfortheincrementalsub- searchrequiressearchingfortwo-hitpairs,wheretwohits
ject offset and the number of repeated positions, there are on the same diagonal and close together, to trigger
still exist a few cases that the increment subject offsets ungappedextensions.Thetraditionalmethod,namelythe
or the number of repeated positions is larger than 255. last-hitarray-basedmethod,iscommonlyusedinquery-
Whensuchsituationsareencountered,wesplitoneposi- indexed BLAST. The last-hit array method uses an array
tionentryintomultipleentriestomakethevaluelessthan torecordthelasthitofeachdiagonal.Whenanewhitis
255. For example, as shown in Fig. 2a, if the increment detected, the algorithm checks the distance between the
subject offset is 300 with 25 positions, then we split the newly found hit and the last hit in the same diagonal of
subjectoffsetintotwoentries,wherethefirstentryhasthe thelast-hitarrayandupdatesthelasthitwiththenewhit.
incrementalsubjectoffset255andthenumberofrepeated Althoughthealgorithmscansthesubjectsequencefrom
position0,andthesecondentryhastheincrementalsub- thebeginningtotheend,thediagonalaccessforanewhit
ject offset 45 for the 25 positions. Similarly, as shown in canberandom.Therandommemoryaccessesonlast-hit
Fig. 2b, for 300 repeated number of positions, the sub- arrays is a critical problem for database-indexedBLAST,
ject offset is split into two entries, where the first entry whichalignsaquerytothousandsofsubjectsequencesat
hastheincrementalsubjectoffset2for255positions,but once (rather than aligning a subject sequence to a single
Zhangetal.BMCBioinformatics (2016) 17:443 Page5of14
query, as is done in query-indexed BLAST). Therefore, from different sequences are put into different bins and
to improve the performance of finding two-hit pairs, we sorted by diagonal ids. After that, we can quickly detect
propose a new method that reorders hits with two-level two-hitpairsbyscanningeverysecond-levelbin.
binning. To improve the performance of the two-hit ungapped
As shown in Fig. 3, each bin is mapped to a diagonal extension further, we filter out the hits that cannot
inthefirstlevelofbinning,andthehitsaregroupedinto be used to trigger the ungapped extension (instead of
binsbydiagonalids,whicharecalculatedbysubjectoff- directly putting all the hits into the second-level bins).
sets minus query offsets. Because query offsets can be Thisoptimization,ascapturedinFig.4,candramatically
calculated by subject offsets minus diagonal ids, we only reduce processing overhead by reducing memory usage,
storethesequenceidsandsubjectoffsetsdirectlyfromthe andinturn,improvingperformance.
indexinordertotominimizememoryusage. Specifically,beforewritingahitintoasecond-levelbin,
Afterthefirst-levelbinning,hitshavingthesamediago- wecheckitsdistancetothelasthitinlast-hitarray.Only
nalidsareplacedintothesamebins.However,ineachbin, if the distance of the current hit to the last hit satisfies
the hits from different sequences are interleaved. Thus, the distance thresholds, i.e., less than threshold_A and
wedesignasecondlevelofbinningtoreorderthehitsby greaterthanorequaltooverlap,thehitcanbeputintothe
sequenceids.Incontrasttofirst-levelbinning,wherethe second-levelbins.Asthenumberofsequencesinaindex
binidisequaltothediagonalid,second-levelbinningsets blockcanbeadjustedbyconfiguringthesizeoftheindex
thebinidtothesequenceid.Becausewescanthebinsof block, the size of the last-hit array may be small enough
thefirst-levelbinningonebyone,thehitsinasecond-level to fit in the cache: not only in the last-level cache (LLC)
bin are sorted naturally by the diagonal id. As shown in ontheHaswellCPUinourevaluationbutalsointheL2
Fig.4,ahitinthesecond-levelbincontainsthesubjectoff- cache. As a result, this optimization to ungapped exten-
setandthediagonalid.Withthesecond-levelbinning,hits sion exhibits excellent data locality when accessing the
Fig.3Anexampleoftwo-levelbinningwithoutfiltering.First-levelbinninggroupshitsintobinsaccordingtotheirdiagonalids.Second-level
binningscanshitsinthefirst-levelbinsbinbybin,andregroupshitsintosecond-levelbinsbytheirsequenceids
Zhangetal.BMCBioinformatics (2016) 17:443 Page6of14
Fig.4Anexampleoftwo-levelbinningwithfiltering.Whilescanninghitsinfirst-levelbins,wecheckthedistanceofeachhittothelasthitinthe
last-hitarray.Onlyifthedistancefitsintothethreshold,thehitcanbeputintothesecond-levelbins
reordered hits, thus improving performance. Moreover, alignment of each query to a thread. However, this
because our optimization filters out the majority of hits, approach results in different threads potentially access-
we also significantly reduce the time spent on memory- ingdifferentindexblocksatthesametime.Inlightofthe
write operations, and in turn, improve performance limited cache size, this approach results in severe cache
further. contention between threads. To mitigate this cache con-
Ifthesubjectoffsetsareunsortedinthedatabaseindex, tention and maximize cache-sharing across threads, we
as shown in Fig. 5a, the binning method can intro- exchangeexecutionorder,asshowninAlgorithm1.That
duce random memory accesses, which would adversely is, the first two stages, i.e., hit detection and ungapped
impact performance. However, sorting the subject off- extension, which share the same database index, access
sets in the database index, as shown in in Fig. 1c, can the same database block for all batch query sequences
resolvethisproblem.Oncetheindexsortingiscomplete, (fromLine6to10).So,weapplytheOpenMPpragmaon
asshowninFig.5b,boththereadsonthedatabaseindex theinnerlooptomakedifferentthreadsprocessdifferent
and the writes on the first-level binning are contiguous, inputquerysequencesbutonthesameindexblock.Then,
thus improving the binning performance via better data threads on different cores may share the database index
locality. that is loaded into memory and even cache. The aligned
results for each index block are then merged together
Optimizationsviamultithreading for the final alignment with traceback, as shown on
In BLAST algorithm, the query sequence is aligned to Line9.
eachsubjectsequenceinthedatabaseindependentlyand For better load balancing, and in turn, better per-
iteratively.Thus,wecanparallelizetheBLASTalgorithm formance, we leverage the fact that we already have a
withOpenMPmultithreadingonthemulticoreprocessors sorted database with respect to sequence lengths. We
inacomputenode,e.g.,ourpairof12-coreIntelHaswell thenpartitionthisdatabaseintoblocksofequalsizeand
CPUs or 24 cores in total. However, achieving robust leverageOpenMPdynamicscheduling.
scalability on such multicore processors is non-trivial,
particularly for a data-/memory-intensive program like Discussion
BLAST, which also introduces irregular memory access In muBLASTP, we use the composition-based statistics
patternsaswellasirregularcontrolflow.Atahighlevel, presented in [27], which is also the default method used
twomajorchallengesexistforparallelizingBLASTwithin in NCBI BLAST. For other composition-based statistics
acomputenode:(1)cacheandmemorycontentionamong methodsinNCBIBLAST,suchas[28],ourcurrentcode
threads on different cores and (2) load balancing among basedoesnotsupportit.Weleavethisworkforthefuture
thesethreads. versions.
Because the alignment on each query is independent, Moreover, the current version of muBLASTP can only
a straightforward approach to parallelization maps the producetheidenticalresultstoNCBIBLASTwhenboth
Zhangetal.BMCBioinformatics (2016) 17:443 Page7of14
Algorithm1OptimizedmultithreadedmuBLASTP
1: functionMUBLASTP(Q,S,I)
2: Input: Q: queries, S: subject sequences, I:
databaseindex
3: Output:Print_Result(R)
4: forIndexblockIbinI do
5: #pragmaompparallelforschedule(dynamic)
6: forQueryQiinQdo
7: H =Hit_Detection(Ib,Qi)
8: U =Ungapped_Extension(H,Sb,Qi)
9: Gi =Gi+Gapped_Extension(U,Sb,Qi)
10: endfor
11: endfor
12: #pragmaompparallelforschedule(dynamic)
13: forQueryQiinQdo
14: Ri =Traceback(Gi,Qi,S)
15: endfor
16: endfunction
use the default output format (i.e., “pairwise” format)
andthedefaultcomposition-basedstatisticsmethod.Asa
result, our software can only generate the similar results
toNCBIBLASTifanyotherparameterisset.Inthefuture
updatesofthissoftware,wewilladdthesupportsfordif-
ferentformats,makingmuBLASTPtobeacomprehensive
toolasNCBIBLAST.
Results
We conducted our experimental evaluations on two dif-
ferentmulticoreCPUplatforms—Haswellplatformand
Nehalemplatform.HaswellplatformconsistsoftwoIntel
Haswell Xeon CPUs (E5-2680v3), each of which has 12
Fig.5Anexampleoffirst-levelbinninghitswithunsortedindexand
cores, 30MB shared L3 cache, and 32KB L1 cache and
sortedindex.Intheexample,thehitsaregeneratedforthewordin
256KB L2 private cache on each core. Haswell plat- queryoffset1.aFirst-levelbinningwithunsortedindex.bFirst-level
form also has 128GB of 2133-MHz DDR main memory. binningwithsortedindex
Nehalem platform consists of two Intel Nehalem Xeon
CPUs(E5645),eachofwhichhas6cores,12MBsharedL3
197aminoacids(orletters),respectively.Thenrdatabase
cache,and32KBL1cacheand256KBL2privatecacheon
consistsofabout85,000,000sequenceswiththetotalsize
eachcore.Nehalemplatformalsohas24GBof1600-MHz
at53GBandwhosemedianlengthandaveragelengthare
DDR main memory. In the experiments, all programs
292and366aminoacids(orletters),respectively.
are compiled by an Intel C/C++ compiler 15.3 with the
Figure6showsthedistributionofsequencelengthsfor
compileroption-O3 -fopenmp.Intheexperiments,all
the uniprot_sprot, env_nr and nr databases. The sizes of
performancenumbersareaveragevaluesofmultipleruns.
most sequences from the two databases lie in the range
from 60 amino acids to 1000 amino acids and with only
a handful of sequences longer than 1000 amino acids.
Databases WeusedthreetypicalproteinNCBIdatabases
Similar observations are also reported in other studies
from GenBank [29]: uniprot_sprot, env_nr and nr. The
[17,30,31].
uniprot_sprot database includes approximately 300,000
sequenceswithatotalsizeof250MBandwhosemedian
length and average length are 292 and 355 amino acids Queries The performance of BLAST depends in part
(or letters), respectively. The env_nr database consists of on the query length. Based on the length distribution
about 6,000,000 sequences with the total size at 1.7 GB shown in Fig. 6, we evaluated the performance of our
andwhosemedianlengthandaveragelengthare177and single-threadmuBLASTPusingthreesetsofquerieswith
Zhangetal.BMCBioinformatics (2016) 17:443 Page8of14
./configureCC=iccCXX=icpc\
—without–gui—without–debug
make
makeinstall
We formatted database, and ran NCBI BLAST with
defaultparameters,asnotedbelow.
makeblastdb<–indatabase><–dbtypeprot>
blastp<–queryquery><–dbdatabase>\
[–num_threadsnumber_of_threads]
Fig.6Sequencelengthdistributionofuniprot_sprot,env_nrandnr As the usage of indexdb program shown above, the
database
index block size is a configurable variable. By default, its
value is set to 128K amino acids (or letters), making the
index block size around 256 KB and fitting into the L2
different lengths — around 100, 500 and 1000 — where cache (256 KB on both Haswell and Nehalem). The rea-
each query set contains 50 queries. For the evaluation son to set the index block size based on the L2 cache is
of our multithreaded muBLASTP, we built three query thatsincetheL2cacheisprivateforeachcore,wecould
batches,eachcontaining100querieswithlengthsaround avoidheavycachecontentionsacrossdifferentthreadsin
100, 500, and 1000, respectively. In addition, we con- themultithreadingmodeiftheindexdatacanbelocated
structedamixed-lengthbatchofsequencesbyrandomly fromtheL2cache.IfmuBLASTPisrunningwithasingle
selecting100queriesofarbitrarysizeinordertoevaluate thread,wecouldincreasetheindexblocksizeandtryto
therealworldperformanceofmultithreadedmuBLASTP, fullyutilizetheL3cacheaswellastheL2cache.Because
especially with respect to scalability and load balancing. increasingindexblocksizemaygeneratemuchmorehits
Table 1 captures the statistical profile of query lengths ineachblock,thepracticalvaluesare2048Kaminoacids
fromourmixed-lengthquerybatchesoftheuniprot_sprot, (letters)onHaswelland1024KaminoacidsonNehalem
env_nr and nr databases, respectively. The details for inourexperimentsforthesinglethreadmode.
queriesaregiveninAdditionalfile1.
To align queries with muBLASTP, as the following
Indexsize
commands, we first formatted and sorted the database
Table2showstherawfile(FASTAformat)sizeforthecor-
using the formatdb and sortdb program, respec-
respondingdatabase(“Database”row),thecorresponding
tively. And then, we indexed the database with a config-
index file size with neighboring words (“Index w/ neigh-
urableblocksizeusingtheindexdbprogram,andfinally
bors”row),theindexfilesizewithoutneighboringwords
alignedqueriesagainstthedatabaseusingthemublastp
(“Index w/o neighbors” row), and the compressed index
program.
file size (“Compressed index” row). Except “Database”
row,thelatterthreerefertothedifferentindexingmech-
formatdb<–idatabase>
anisms presented in this paper. According to Table 2,
sortdb<–idatabase><–osorted_database>
indexdb<–isorted_database>\ the database index with neighboring words, when com-
pared to the database index without neighboring words,
[–sblock_size(Kletters),default128]
mublastp<–iquery><–dsorted_database>\ can be on the order of 20 times larger. Index com-
pression achieves 1.47-fold compression rate for the
[–tnumber_of_threads]
uniprot_sprotdatabase,1.46-foldcompressionrateforthe
In experiments, we compared muBLASTP with NCBI env_nr database, and 1.47-fold compression rate for the
BLAST (version 2.3.0), which was configured and built nrdatabase.Asaresult,thecompressed(database)index
withthefollowingcommands. for the uniprot_sprot database is 2 times the size of the
Table1Statisticsofquerylengths(aminoacids)inmixed-length Table2Sizeofdatabaseandindexfilesingigabytes(GB)
querybatches uniprot_sprot env_nr nr
Targetdatabase Averagelength Medianlength Maximumlength Database 0.25 1.89 52.4
uniprot_sprot 333 289 1187 Indexw/neighbors 18.1 116.6 N/A
env_nr 191 175 504 Indexw/oneighbors 0.76 5.82 122.3
nr 312 263 1127 Compressedindex 0.51 3.97 83.1
Zhangetal.BMCBioinformatics (2016) 17:443 Page9of14
originaldatabasewhileitis1.8timesthesizeoftheorig- Single-threadedmuBLASTPvs.single-threadedNCBIBLAST
inal env_nr database, and it is 1.6 times the size of the Figure 7 shows the speedups of singled-threaded
original nr database. Because we embedded the offsets muBLASTP over single-threaded NCBI BLAST on
to neighboring words into the database index, our index Haswell platform, using different query lengths.
without neighboring words can achieve identical results muBLASTP achieves 2.22∼3.35-fold, 1.17∼1.7-fold,
astheindexwithneighboringwordsbutwithsignificantly and 1.06∼1.3-fold speedups over NCBI BLAST on
lessmemoryusage. the uniprot_sprot database with queries of length 100,
500, and 1000, respectively. For the env_nr database,
Performancecomparisonforalignmentstages muBLASTP achieves 2.24∼3.51-fold, 1.3∼1.77-fold, and
To evaluate the performance improvement with index 1.26∼1.39-fold speedups with queries of length 100, 500
structure and re-factored BLAST algorithm, we used and 1000, respectively. For the nr database, muBLASTP
gettimeofday() functions to measure the execution achieves 2.3∼4.41-fold, 1.34∼1.5-fold, and 1.21∼1.26-
timeofallfouralignmentstagesforbothmuBLASTPand fold speedups with queries of length 100, 500 and 1000,
NCBIBLASTwithoutI/O. respectively. muBLASTP achieves higher speedup on
Fig.7Speedupforalignmentstagesofsingle-threadedmuBLASTP Fig.8Speedupforalignmentstagesofsingle-threadedmuBLASTP
oversingle-threadedNCBIBLASTonHaswellplatformwithdifferent oversingle-threadedNCBIBLASTonNehalemplatformwithdifferent
querylengthsonuniprot_sprotdatabase(a),env_nrdatabase(b)and querylengthsonuniprot_sprotdatabase(a),env_nrdatabase(b)and
nrdatabase(c) nrdatabase(c)
Zhangetal.BMCBioinformatics (2016) 17:443 Page10of14
Table3SpeedupforalignmentstagesofmultithreadedmuBLASTPovermultithreadedNCBIBLASTonHaswellplatform(with24
threads)andNehalemplatform(with12threads)withquerybatchesofdifferentquerylength
Haswellplatform Nehalemplatform
Querylength uniprot_sprot env_nr nr uniprot_sprot env_nr nr
Aver SD Aver SD Aver SD Aver SD Aver SD Aver SD
100 5.7 0.013 2.67 0.013 1.94 0.021 8.2 0.013 4.52 0.017 8.59 0.015
500 3.22 0.012 1.44 0.012 1.39 0.014 2.25 0.018 1.69 0.018 2.69 0.016
1000 2.85 0.016 1.35 0.013 1.20 0.015 1.76 0.018 1.41 0.016 1.72 0.015
Mixed 4.16 0.013 1.41 0.012 1.49 0.012 2.82 0.018 2.34 0.013 4.1 0.015
Averistheaveragevalueofthreeruns,andSDisthestandarddeviationofthreeruns
the larger database because the BLAST algorithm on a uniprot_sprot database, up to a 4.52-fold speedup when
large database needs to process significantly more hits, using the env_nr database, and up to a 8.59-fold
i.e., spending more time on hit detection and two-pair speedup when using the nr database. In this case,
hit ungapped extension, which are the stages that our muBLASTP achieves much higher speedups on the
optimizationsfocuson. smaller uniprot_sprot database, which indicates that
Figure 8 shows the speedups of singled-threaded muBLASTPdeliversbetterscalabilitythanNCBIBLAST
muBLASTP over single-threaded NCBI BLAST on onsmallerdatabases.
Nehalem platform, using different query lengths. We also tested muBLASTP performance with query
muBLASTP achieves 1.6∼3.17-fold, 1.33∼1.47-fold, batchesofmixedlengths.Table3showsthatonHaswell
and 1.09∼1.3-fold speedups over NCBI BLAST on platform muBLASTP achieves a 4.16-fold speedup over
the uniprot_sprot database with queries of length 100, NCBI BLAST on uniprot_sprot database, a 1.41-fold
500, and 1000, respectively. For the env_nr database, speedup over NCBI BLAST on env_nr database, and
muBLASTP achieves 2.38∼3.8-fold, 1.3∼1.53-fold, and a 1.49-fold speedup on nr database. Table 3 also
1.14∼1.25-fold speedups with queries of length 100, 500 shows that on Nehalem platform muBLASTP achieves a
and 1000, respectively. For the nr database, muBLASTP 2.82-fold speedup over NCBI BLAST on uniprot_sprot
achieves 2.0∼3.21-fold, 1.1∼1.49-fold, and 1.00∼1.25- database, a 2.34-fold speedup over NCBI BLAST
fold speedups with queries of length 100, 500 and 1000, on env_nr database, and a 4.1-fold speedup on nr
respectively. database.
MultithreadedmuBLASTPvs.multithreadedNCBIBLAST MultithreadedmuBLASTPvs.single-threadedmuBLASTP
When using query batches of different lengths, Table 3 We also evaluated parallel efficiency of multithreaded
shows that our multithreaded muBLASTP on Haswell muBLASTP. Table 4 shows that multithreaded
platform achieves up to a 5.7-fold speedup over mul- muBLASTP using 24 threads on Haswell platform can
tithreaded NCBI BLAST when using the uniprot_sprot achieve 19.8∼21.6-fold speedups over single-thread
database,uptoa2.67-foldspeedupwhenusingtheenv_nr muBLASTP with query batches of different lengths on
database,anduptoa1.94-foldspeedupwhenusingthenr differentdatabases.Table4alsoshowsthatmultithreaded
database. muBLASTP using 12 threads on Nehalem platform can
Table 3 shows that our multithreaded muBLASTP on achieve 10.7∼11.6-fold speedups over single-thread
Nehalem platform achieves up to a 8.2-fold speedup muBLASTP with query batches of different lengths on
over multithreaded NCBI BLAST when using the differentdatabases.
Table4SpeedupforalignmentstagesofmultithreadedmuBLASTPoversingle-threadedmuBLASTPonHaswellplatform(with24
threads)andNehalemplatform(with12threads)withquerybatchesofdifferentquerylength
Haswellplatform Nehalemplatform
Querylength uniprot_sprot env_nr nr uniprot_sprot env_nr nr
Aver SD Aver SD Aver SD Aver SD Aver SD Aver SD
100 19.8 0.013 20.2 0.016 20.0 0.013 10.7 0.014 11.2 0.012 11.6 0.012
500 20.9 0.022 20.6 0.011 21.4 0.017 10.9 0.012 11.5 0.013 11.5 0.012
1000 21.4 0.013 21.2 0.013 21.5 0.012 10.8 0.012 11.3 0.013 11.6 0.018
Mixed 21.4 0.012 21.5 0.017 21.6 0.013 10.9 0.018 11.3 0.018 11.5 0.015
Averistheaveragevalueofthreeruns,andSDisthestandarddeviationofthreeruns