Accele a ing Exhaus i e Pai wise Me agenomic
Compa isons
Es eban P´e ez-Wohl eil1, Osca To eno1, and Oswaldo T elles1

?
1Depa men o Compu e A chi ec u e, Uni e si y o Malaga,
Boule a d Louis Pas eu 35, Malaga, Spain
{es ebanpw,osca ,o elles}@uma.es
Abs ac . In his manusc ip , we p esen an op imized and pa allel
e sion o ou p e ious wo k IMSAME, an exhaus i e gapped aligne o
he pai wise and accu a e compa ison o me agenomes. Pa alleliza ion
s a egies a e applied o ake ad an age o mode n mul ip ocesso a chi-
ec u es. In addi ion, sequen ial op imiza ions in CPU ime and mem-
o y consump ion a e p o ided. These algo i hmic and compu a ional en-
hancemen s enable IMSAME o calcula e nea op imal alignmen s which
a e used o di ec ly assess simila i y be ween me agenomes wi hou e-
qui ing e e ence da abases. We show ha he o e all e iciency o he
pa allel implemen a ion is supe io o 80% while e aining scalabili y as
he numbe o pa allel co es used inc eases. Mo eo e , we also show ha
sequen ial op imiza ions yield up o 8x speedup o scena ios wi h la ge
da a.
Keywo ds: High Pe o mance Compu ing ·Pai wise Compa ison ·Pa -
allel Compu ing ·Nex Gene a ion Sequencing ·Me agenome Compa i-
son
1 Backg ound
A me agenome is de ined as a collec ion o gene ic ma e ial di ec ly eco e ed
om he en i onmen . In pa icula , a me agenome is composed o a la ge num-
be o eads (DNA s ings) d awn om he species p esen in he o iginal popu-
la ion. To his day, he ield o compa a i e me agenomics has become big-da a
d i en [1] due o new echnological imp o emen s in high- h oughpu sequenc-
ing. Howe e , he analysis o la ge me agenomic da ase s ep esen s a compu a-
ional challenge and poses se e al p ocessing bo lenecks, specially o sequence
compa ison algo i hms.
T adi ional me agenomics compa ison in ol e in e media e pai wise (and in-
di idual) compa isons agains a e e ence da abase. This p ocedu e allows o ex-
ac a mapping dis ibu ion be ween eads and species, and hus enables o la e
on compa e hese dis ibu ions. A simila i y measu e can hen be compu ed om
he wo dis ibu ions. Howe e , due o he unknown and complex composi ion
?Co esponding au ho .
o me agenomes, adi ional compa isons based on a e e ence equi e da abases
o be la ge, which o en in oduce bias and d as ically inc ease execu ion imes.
In his line, di ec compa isons be ween me agenomes o assess o e all simila i y
gain in e es as unning imes can be sho ened and bias a oided. Fu he mo e,
he e is ye no accep ed consensus on how simila i y should be assessed, and in
se e al scena ios ce ain y comes a he expense o exhaus i e and op imal align-
men s. S ill, op imal alignmen o la ge da ase s is no easible wi hou making
use o pa allel in as uc u e and op imiza ion echniques.
Nex Gene a ion Sequencing pla o ms a e gene a ing la ge amoun s o da a
pe un, o highe quali y and a a lowe p ice. Howe e , he pe o mance o
compu a ional app oaches used o p ocess me agenomic da a su e in e sely
p opo ional o ha o he size o gene a ed samples. High Pe o mance Com-
pu ing echniques can be applied in o de o o e come he p ocessing bo lenecks
and accele a e unning imes.
Se e al pa allelism s a egies ha e been al eady applied o so wa e o bo h
compa a i e genomics and me agenomics, such as he mul ipu pose BLAST [2]
amily, whe e di e en ypes o a chi ec u es ha e been exploi ed (e.g. mpi-
BLAST [3] o dis ibu ed memo y o TERABLAST [4] o i s use on FPGAs).
Pa alleliza ion o sequence alignmen algo i hms ha e also been applied o GPUs
such as GSWABE [5] o CUSHAW2-GPU [6]. FPGAs ha e also been employed
o accele a e sequence compa isons (e.g. SWAPHI [7]). Howe e , GPUs and FP-
GAs a e expensi e and hei speci ici y o en o ce he use o a educed sub-
se o p og ams due o pla o m dependence es ic ions. O he gene al pa allel
app oaches which make use o CPU mul i h eading such as BOWTIE [8] o
PARALLEL-META3 [10] use POSIX h eads [11] o UNIX-based en i onmen s
in a sha ed memo y a chi ec u e.
Howe e , he abo e men ioned sequence aligne s a e no speci ically designed
o compa e ead o eads o con igs and pa icula ly no o assess simila i y be-
ween me agenomes. Fo ins ance, BOWTIE wo ks bes when aligning sho
eads o la ge e e ence genomes. Mo eo e , in [12] he e ec o in oducing
in e media y agen s o ul ima ely assess simila i y be ween me agenomes was
a gued. In his line, al e na i e app oaches ha we e capable o di ec compa -
isons we e discussed (e.g. MASH [13], SIMKA [14]) and when possible (BLAST,
COMMET [15]), compa ed o IMSAME. Addi ionally, i was shown ha coa se-
g ained app oaches o me agenomics compa ison could lead o esul s ha we e
highly dependen on hype pa ame e s (such as he ini ial seed size).
To ul ill he gap, IMSAME was p esen ed as a pa allel, ine-g ained, and
exhaus i e gapped ead- o- ead (including con igs and sca olds) aligne . In his
manusc ip , we p esen an op imized e sion o IMSAME which is able o com-
pu e as e while using a linea and con olled amoun o memo y. Mo eo e ,
High Pe o mance Compu ing echniques ha e been applied o balance he wo k-
load among h eads o educe h ead synch oniza ion.
2
2 Me hods
IMSAME (”Inc emen al Mul i-S age Alignmen o ME agenomes”) is in ended
o compa e eads o eads di ec ly, i.e. wi hou using a e e ence da abase, and
o assess simila i y be ween hem while p o iding a con iden le el o ce ain y. I
p oceeds by combining a di e en se o alignmen - ee, gapped- ee and gapped
alignmen s. Each o hese p ocedu es is in ended o yield a di e en le el o
speed and sensi i i y, depending on he alignmen s age. Fo ins ance, he ini ial
de ec ion o seeds be ween eads ( o be e e ed as hi s) is pe o med using k-
me s (wo ds o leng h k), whe eas p obabilis ic il e ing is applied when hi s a e
ex ended in o High-sco ing Segmen Pai s (HSPs). HSPs wi h su icien ly small
p obabili y o belonging by chance o he unde lying dis ibu ion a e kep and
used as ancho s o a bounded Needleman-Wunsch (NW in ad ance) global align-
men [16]. Figu e 1 shows he o e all a chi ec u e o IMSAME. The ollowing
sec ions illus a e each o he me hods employed in IMSAME.
The compu a ion o hi s, ex ended agmen s (HSPs) and gapped alignmen s
will o en ep esen mo e han 85% o he compu a ion ime. The e o e he
pa allel s a egy in IMSAME is ocused in hese s ages, whe eas loading he
da abase and wo kload gene a ion and dis ibu ion is pe o med sequen ially.
2.1 Compu a ion o alignmen s
This sec ion depic s he in e nal p ocedu e ollowed by IMSAME o compu e
pai wise alignmen s be ween he sequences con ained wi hin he inpu s.
Hash able gene a ion and diagonal il e ing A hash able is buil o
wo ds o size 12 (i.e. 12-me s) o he e e ence me agenome. This p ocedu e
s a s by linea ly scanning he e e ence me agenome and adding an en y in
he hash able o each 12-me . The posi ion in he ile and he ead numbe o
which i belongs is s o ed. Each en y in he hash able will hold a linked lis
in o de o handle collisions. Since he e e ence me agenome is conside ed as a
la ge sequence (i.e. he coo dina es o eads is global in espec o he ile), hen
he inse ion o 12-me s in he hash- able is so ed in e ms o he diagonal (as
in [17]) be ween que y and e e ence me agenome.
Hi s de ec ion and ex ension o agmen s Once he hash able is buil
o he e e ence me agenome, he algo i hm p oceeds by loading he que y
me agenome and ma ching 12-me wo ds o hose s o ed in he hash able. These
hi s se e as seeds o ex end he alignmen s. Fo e e y hi , a linea , ungapped
ex ension which allows mu a ions -bu no indels- is pe o med in bo h di ec-
ions, o wa d and backwa d espec i e o he sequence. This ex ension wo ks
by op imizing a sco ing unc ion ha akes in o accoun he leng h and numbe
o sha ed iden i ies.
3
Fig. 1. T adi ional me agenomics compa ison (le ) and o e all diag am o he wo king
p ocedu e o IMSAME ( igh ). Le : (a) Bo h me agenomes a e compa ed indi idually
agains a chosen e e ence (b) da abase. (c) Indi idual esul s a e me ged in o a uni-
ied (d) esul which p opaga es indi idual biases. Righ : (a) and (b) ep esen he
inpu me agenomes. (c) Compu a ion o he hash able o ini ial seeds. (d) Wo kload
dis ibu ion o h eads. (e) Each h ead asks o a job (i.e. block o eads). ( ) Each
h ead de ec s hi s be ween eads. (g) A e de ec ion o hi s, an HSP is compu ed by
ex ending he hi linea ly. (h) I he compu ed HSP has no been gene a ed by chance,
hen a gapped alignmen be ween he wo eads is pe o med.
Ancho ed gapped alignmen In o de o apply a bounded compu a ion o
bo h CPU ime and memo y equi emen s, i is necessa y o use heu is ic me h-
ods o explo e a educed subspace o he whole sea ch space. Hence, he quali y
o he esul s will be di ec ly a ec ed by he used heu is ic me hod. IMSAME
uses a simple ye powe ul ancho ing p ocedu e, which is illus a ed as ollows:
Once a hi has been de ec ed and ex ended, he expec ed alue o he esul ing
ex ended agmen is compu ed. I he expec ed alue is su icien ly small, hen
i is used as ancho o he alignmen . A s aigh line is compu ed be ween he
global s a (0,0) o he wo sequences and he ancho ed agmen . Ano he line
is compu ed om he ending o he ancho ed agmen o ha o he wo se-
quences. These wo lines a e hen disc e ized using B esenhams [18] algo i hm.
The B esenhams algo i hm will de ine a disc e e succession o numbe s ha
ep esen he guides o he alignmen p ocedu e. A window o a iable size is
used o explo e a subse o cells in he NW ma ix o he le and o he igh
4
espec i e o he cen e o he guides (and hus con o ming a window). This p o-
cedu e enables a much as e compu a ion based on he educed sea ch space. A
he same ime, high quali y alignmen s a e s ill p oduced due o he ancho ed
compu a ion o e egions ha a e known o be simila .
Bounded compu a ion in CPU ime and memo y Since only a subse o
he sea ch space is explo ed, less memo y is equi ed o s o e he able compu ed
in he dynamic p og amming algo i hm. The e o e, he size o he bounded win-
dow will de e mine he educ ion in memo y and CPU ime o e he algo i hm.
Gene ally, a NW algo i hm will equi e O(n2) ime and space in he size o he
inpu . In his case, o sequences o leng h nand m,O(nm) will be equi ed,
which g ows quad a ically. Using a bounded window in he compu a ion educes
one o he a iables o a cons an . The e o e i we ake nas he size o one o he
inpu s and kas he size o he window, we will ha e ha k << n. Hence, he
space and ime complexi y d ops om quad a ic o linea in he size o he inpu .
Howe e , in ce ain cases he di e ence in leng h be ween he sequences o be
aligned can be conside ably la ge. Fo example, conside wo sequences whose
leng h di e ence is o one o de o magni ude. I his di e ence is no aken in o
accoun , he algo i hm will p obably explo e ma ix cells ha a e ou side he
bounda ies o he ancho ing, and hus esul ing in mo e compu a ion ime. In
his sense, we p opose using he geome ic mean (√nm) o p oduce a mapping
be ween sequences leng h a io and window size. Mo eo e , kis u he adjus ed
applying a use -de ined pa ame e . No ice ha he geome ic mean will no as-
sign he same window size o wo di e en se s o sequences s1and s2o leng h
10l+1,10l+1 and s3,s4o leng h 10l+2,10l, al hough he size o he sea ch space
would be equal.
2.2 Dynamic wo kload pa i ioning and dis ibu ion o h eads
IMSAME uses a modi ied Guided Sel Scheduling (GSS) [19] o handle wo kload
assignmen . In his line, he que y me agenome ( he e e ence me agenome is
only p ocessed o gene a e he hash able) is sepa a ed in o Mpa i ions each o
which will con ain misubpa i ions, wi h i anging om 1 o M. Each o hese
subpa i ions will hold a numbe o eads ha will be dynamically assigned o
h eads using a h ead-sa e queue when hese un ou o wo k. Mo eo e , he
numbe o pa i ions is use -de ined, and can be adjus ed depending on he size
o he inpu s. Each subpa i ion is likewise di ided by he numbe o h eads
in o blocks o eads. Thus, he numbe o eads con ained in each block is
de e mined by he cu en le el o pa i ioning i, he o al numbe o eads R
and he numbe o h eads . The exp ession o calcula e he numbe o eads
assigned o a block a pa i ion iis as ollows:
Bi=
R
M
∗i=R
M∗ ∗i(1)
Whe e Ris he o al numbe o eads in he me agenome, Mis he numbe
o pa i ions, is he numbe o a ailable h eads and iis he cu en pa i ion.
5

Tha is, he me agenome is i s ly di ided by he numbe o pa i ions, and each
o hese is hen di ided by he numbe o h eads mul iplied by he subpa i ion
dep h. Figu e 2 shows how he que y me agenome is decomposed pe pa i ion.
Fig. 2. Que y decomposi ion in o wo kload blocks. Ini ially, he que y me agenome is
di ided in o pa i ions ( ed solid lines). Each pa i ion is again di ided by he cu -
en pa i ioning le el ( e ical dashed lines). Finally, each subpa i ion is di ided in o
blocks o equal numbe o eads depending o he numbe o h eads a ailable.
The wo kload dis ibu ion unc ion belongs o he amily o 1/x unc ions
and shows a decay in he size o blocks in he ea ly pa i ions in o de o assign
smalle jobs o he h eads as hese a e consumed. Addi ionally, he unc ion
shows a ho izon al asymp o e ha gua an ees blocks o eads o a minimum size
despi e he numbe o pa i ions chosen.
3 Resul s and discussion
Two sepa a e compa isons we e ca ied ou in o de o es he wo di e en
aspec s ha ha e been imp o ed o e he o iginal IMSAME e sion. All se-
quences used in he compa isons belong he o he Human Mic obiome P ojec 1
(HMP), in pa icula o he Illumina WGS Assemblies. The un iden i ie s a e
p o ided a each compa ison pe o med in o de o allow ep oducibili y. The
es ing scena io is se up as ollows:
1. A single compa ison in ol ing sca olds o es sequen ial wi h op imiza ions.
2. A compa ison in ol ing h ee di e en da ase s a ying in sizes, om small
o la ge. These compa isons will se e o accoun o he scalabili y o he
pa allel imp o emen s in espec o he sequen ial e sion.
A las , he speedup om bo h op imiza ion pe spec i es a e discussed and
add essed.
1h p://hmpdacc.o g/HMASM/
6
3.1 In as uc u e
The Picasso supe compu e loca ed a he Uni e si y o Malaga (Malaga, Spain)
[20] was used o es he pa alleliza ion s a egies. The compu a ion was pe -
o med using only he a nodes which con ain 8 In el E7-4870 p ocesso s and 2
TB o RAM each. The s o age is managed by a Lus e ile sys em suppo ed by
a DDN s o age ack wi h i e h ee-dimensional disk enclosu es and wo edun-
dan SFA10000 con olle s. The execu ions desc ibed in his manusc ip ange
om 1 o 32 co es inc easing by s eps o powe s o wo. Run ime execu ions we e
measu ed using he ime command om UNIX-based en i onmen s.
3.2 Compa ison be ween he o iginal and bounded IMSAME
Besides he imp o ed pa alleliza ion s a egies, IMSAME has addi ionally been
imp o ed wi h bounded CPU ime and memo y usage. In o de o pe o m com-
pa isons o measu e he speedup p oduced by he pa alleliza ion echniques, he
imp o emen s a e es ed be ween he o iginal IMSAME and he bounded in a
sequen ial ashion. Fo his pu pose, wo uns composed o sca olds, namely
SRS016105 and SRS017451 we e aken om he HMP da abase and compa ed
using bo h e sions. Since sca olds a e longe han eads, highe penal ies we e
used o a ine gap model (−8 o inse ion o gap and −4 o ex ension o
gap). Resul s emained equal o bo h execu ions, al hough a ia ions can be ob-
se ed i la ge-scale ea angemen s ake place (e.g. long ange ansposi ions).
The o iginal e sion o IMSAME ook 7 minu es and 29 seconds, whe eas he
bounded e sion ook 55 seconds, ep esen ing a speedup o app oxima ely 8x.
This speedup is mos ly p oduced by wo ac s (1) he diagonal il e ing p o-
cedu e educing he numbe o linea ex ensions pe o med p io o a gapped
alignmen and (2) he bounded window applied o he NW algo i hm, which
subs an ially educes he space sea ch. Howe e , i is impo an o no e ha he
la e speedup is di ec ly p opo ional o he size o he eads (i.e. smalle eads,
smalle speedup).
3.3 Speedup e alua ion o he pa alleliza ion s a egy
The speedup in oduced by he pa alleliza ion s a egy is measu ed by using
h ee da ase s: (1) a small-sized one, (2) a medium-sized one and (3) a la ge-
sized one. Table 1 summa izes he h ee da ase s. This p ocedu e enables us o
e alua e he s abili y o he speedup as a unc ion o he inpu size.
Table 2 shows he execu ion imes o each o he da ase s along wi h he
speedup and e iciency o each execu ion. In he same line, Figu e 4 shows he
speedup e alua ion plo . The speedup is calcula ed as he ime needed by he
algo i hm un using only one co e di ided by he ime needed using mo e co es.
The e iciency is calcula ed as he a io be ween he achie ed speedup and he
op imal speedup (equal o he numbe o co es). As can be seen in Table 2 and
Figu e 4, he speedup is nea ly op imal in he scena io o enough da a, i.e. he
la ge da ase . Howe e , in he case o he small and mid-sized da ase , a decay
7
Table 1. Summa y o he da ase used o he speedup e alua ion. F om le o igh :
(1) Me agenome pai s compa ed, (2) Sum o eads om bo h me agenomes, (3) Sum
o he size o bo h me agenomes in megaby es and (4) A e age size o eads in base
pai s.
Da ase (Run ID) Numbe o eads Size (MB) A e age ead leng h (bp)
SRS017697
SRS019119 613,983 77 91
SRS064376
SRS065347 3,401,514 463 100
SRS018359
SRS057022 14,295,910 1809 99
in e iciency can be obse ed due o he size o da a no being la ge enough o
he numbe o co es. Mo eo e :
1. Small da ase (in pu ple in Figu e 4): The peak o e iciency is achie ed using
2 and 4 co es ( eaching 84% e iciency). When using mo e han 4 co es, he
da ase size becomes oo small and some h eads become inac i e while o h-
e s a e s ill p ocessing. To imp o e e iciency on he small da ase , a highe
numbe o pa i ions should be used o emo e h ead balance synch oniza-
ion a he end o he compu a ion.
2. Medium da ase (in ed in Figu e 4): The medium da ase ep esen s a 143%
inc ease in size o e he small da ase , and shows a much highe e iciency,
wi h a peak a 86% using 8 co es. Howe e , simila ly o he smalle case, he
e iciency (and hus he speedup) decays when using 32 co es.
3. La ge da ase (in g een in Figu e 4): The la ge da ase shows he o e all
bes e iciency and speedup, wi h a peak a 8 co es and slow decay up o 32
co es, whe e 81% e iciency is achie ed. Howe e , he ac ha he speedup
is op imal a 2 and 4 co es indica es ha p obably s ill mo e pa i ioning
le els a e equi ed in o de o a oid h ead synch oniza ion.
Table 2. Execu ion imes, speedup and e iciency o he execu ions o IMSAME using
om 1 o 32 co es. The ows indica e he numbe o co es whe eas he columns e e
o ime consump ion (in seconds), speedup and e iciency pe each o he da ase s.
Small Medium La ge
Co es Time (s) Speedup E iciency Time (s) Speedup E iciency Time (s) Speedup E iciency
1 425 1.00 1.00 8,634 1.00 1.00 76,253 1.00 1.00
2 252 1.69 0.84 4,770 1.81 0.91 38,009 2.01 1.00
4 126 3.37 0.84 2,571 3.36 0.84 18,867 4.04 1.00
8 76 5.59 0.70 1,262 6.84 0.86 9,838 7.75 0.97
16 42 10.12 0.63 683 12.64 0.79 5,282 14.44 0.90
32 30 14.17 0.44 388 22.25 0.70 2,937 25.96 0.81
8
Fig. 3. The speedup is shown o he di e en da ase s (in pu ple, g een and ed) along
wi h he op imal speedup (in blue). The x-axis shows he numbe o co es used pe
compa ison, whe eas he y-axis shows he calcula ed speedup in espec o he numbe
o co es used.
4 Conclusions
In his manusc ip , we ha e shown an op imized e sion o IMSAME in which
we applied wo pa alleliza ion s a egies, namely (1) a dynamic schedule o
he dis ibu ion o wo k and (2) an n-le el pa alleliza ion in he compu a ion o
alignmen s using POSIX h eads in a sha ed-memo y en i onmen . We ha e also
applied se e al sequen ial imp o emen s o e he o iginal e sion, which ha e
imp o ed he o e all algo i hm complexi y and e iciency. Addi ionally, we ha e
ca ied ou wo sepa a e compa isons o p o e he pe o mance o IMSAME,
ha is, i s ly, one o alida e he sequen ial imp o emen s o e he o iginal
e sion and secondly, ano he one using h ee di e en da ase s anging in sizes
o e alua e he achie ed speedup and he pa allel e iciency.
In o de o keep de eloping IMSAME, we a e cu en ly wo king on:
1. Pa alleliza ion o he loading s age.
2. Imp o e wo kload dis ibu ion by building a eg ession model o au oma i-
cally se he numbe o pa i ioning le els.
3. Use ROC cu es [21] o se he op imal pe cen age h esholds.
4. Use gene ic algo i hms o de e mine op imal scheduling.
Acknowledgmen s
This wo k has been pa ially suppo ed by he Eu opean p ojec ELIXIR- EX-
CELERATE (g an no. 676559), he Spanish na ional p ojec s Pla a o ma de Re-
cu sos Biomolecula es y Bioin o m icos (ISCIII-PT13.0001.0012) and RIRAAF
(ISCIII-RD12/0013/0006) and he Uni e si y o Malaga.
9