2018年第3湖_____________________________________________________
文章编号=1009 -2552 (2018)03 -0067 -06 DOI：10. 13274/jki.hdzj.2018. 03. 014
"f s息技术基于G PU的BW A序列比对算法分析与加速
海玲\刘俊霞\海志民2，刘岩\杨嘉鹏\刘智勇1
(1.新疆工程学院电气与信息工程系，乌鲁木齐830052; 2.大唐环境产业集团股份有限公司，新疆呼图壁831200 )
摘要：文中实现了 G PU平台加速的BW A-M EM算法，将BW A-M EM算法中的两个热点模块：SM EM查和ch a in生成模块利用G PU平台进行加速，通过重构算法流程、精简需要向CUDA设 备传输的数据结构，采用合理的任务划分方式来提升BW A-M EM在G P U平台的性能。论文对BW A-M EM算法的特点进行了深入分析，总结了 BW A-M EM算法在G PU平台加速效果受到限制的原因。
关键词：全球气候模式；谱方法；图形处理器；并行计算
中图分类号：TP391 文献标识码：A
Analysis and accelerating of BWA sequeuce comparison algorithm
based on GPU
HAI Ling1，LIU Jun-xia1，HAI Zhi-min2，LIU Yan1，YANG Jia-peng1，LIU Zhi-yong1
(1. Department of Electrical and Information Engineering，Xinjiang Institute of Engineering，Urumqi 830052，China；
2. Datang Environmental Industry Group Co.，Ltd.，Hutubi 831200，Xinjiang，China) Abstract:This paper implements the BWA-MEM algorithm on CUDA platform.SMEM finding and
chain generating,two of the hotspots，are accelerated with GPU devices.It improves the performance of GPU-based BWA-MEM by reconstructing the algorithm，reducing the data that has to be transferred between host and device and developing a reasonable task partition method.It analyzed the characteristics of the BWA-MEM algorithm，and summed up the reason why BWA-MEM algorithm is lim ited in GPU platform.
Key words:community atmosphere model;spectral element;GPU;parallel computing
0引百
生物序列比对是计算机领域经典的文本比对问 题在生物领域的应用[|]。随着不断涌现的新型分 子生物学技术的发展，随之而来的比如基因变异、R N A表达、蛋白质和基因相互作用等分子生物学研 究需要研究人员采用高通量的方法去解释。这对高 性能计算提出了新的挑战。下一代测序又称高通量 测序，它的一大特点是高通量，如最新的Illu m in a测 序仪HiSeq2500在“高输出”（High Output)模式下在 11天能一次性产生600G B数据，而在‘快速’（Rapid Run)模式下能在约27小时内产生120G B数据[2]。这种规模的数据超出了中小型实验室的数据处理能 力。高通量测序时代的挑战不再是数据的产生，而是数据存储，处理和分析。第三代测序技术的发展将进 一步加快测序速度，同时产生较长的测序片段，对序 列比对技术的发展提出了更高的要求。
与此同时，计算机硬件平台也在近几年飞速发展, 不断推陈出新，新的多核平台和众核平台不断出现，众 核协处理器加速设备如GPU使高性能计算机的性能 快速提高,GPU异构平台正成为构建高性能计算机的 重要方式。这也使得计算机体系结构变得越来越复
收稿日期：2017-07 -07
基金项目：新疆维吾尔国际科技合作计划项目（20156007)；新疆工程学院博士基金项目（2015B Q J021812)；新疆维吾尔自
治区高校科研计划自然科学基金项目（XJE D L2014S074,
XJE D L2016S106)
作者简介：海玲（1982-)，女，硕士研究生，研究方向为信号处理与检测，并行计算。
—67 —
^ ^
s l s s l i H ^S I O B W A
S 胀驾商油_^勹戍'辟黯
数据结构与算法论文_靖3许
s $^l , , s m
S B W A -M E M «»^
胀驾s i ?
嗤簦项>$1许3^ )N #带#^s f f 黯，裥许#簦项翁筠s  $
虽«辂^丨。羅_*&_簦项漭讳S -M S 姊I  S S S 0^$,)N #
»m s 芻
將：t ^〇M I t r v ,*)4$-H 泠爾浮瞄骓 卡卟，
X 4B W A
雜湼裥BWA -M EM «
辂掛泠爾卡卟§y s 0030^0^$s ^ I «n «辂諷涪m ^^s f f 黯。1
BW A
犍硌知冷導L l »鉼B W A
裥丨->茂窟姊齊项当民笃里丨->汁嗤# 姗$屈降(民宫洚A
滌$M l a )h r 姙S 粦羊电r 5.M l M ^^a H _
^^^r a ^r f  (r a 日 3W
S -W h e e le r  T r a n s -f o m o s ^H ^E  鄯娜c li l  雒t I 〕。^&11|->«辂 laa:B W A -backtrack, B W A -S W  fP  B W A -M E M 。
BWA-backtrack
M m ^
-M e s f t loobpsillumina
簦项 r e a d s ^i -f s 〇^^a ->M m ^m ^s ^s ,
^ 70bp
M  I M b p m i ^i -f s 。B W A -M E M  fp  B W A -S W
#[s l s i ^5s ^^ 3—{M
裥 B W A -M E M ,m M m _i ?s «»,M m i i #^s #a  来鸯鲥^_s i u s ,M ^[v }m ^,m f »S A a w 70b p 〜
loobpsillumina  read s h r , B W A -M E M i&tb  B W A -
b a
c k t r a c k #m ^s f f t o ^BWA -backtrack
)t 油驾鑛葛燁攝雒菡菊$M I B ^s ,)t a b a c k w a r d  search®^
碎来|^課 r e a d s
掛 #<i *^s h r s f <l ;w 。^^讁漭H I e , B W A -b a c k t r a c k l ^l >«»^a l s >>^,^^s ^->l s ^ /J N I ^S r e a d s ^g ?w ^<i *^s F M -In d e x
h r s  v a l ,N t M i n t e r v a l ^^,^s l ^i t j s ^s w ^<i *^s i n
t e r v a l ^M l ^i t I ^s w ^<i *^s ^s ^s 0f #f <l ;s 〇M i  ,a ^s i t j ^s 胀t t r f <l ;s , c(
a )^#<i *^ s 丑•f t ^/J N -H a s y l w +l p o ,, 04a
掛 t a o h 丑圧芻s 涔燁。l >3(s 3=c (a )+o (s ,l >3(3-l )+l  >3l (a s =c s +o (a ,>3l (s )函一 B W A M
臶回連 i i -M >>^驾^带雜#H m ,^a b o u n d e d  search^
碎u p «
辂1洚邻。p r a 5
c p lc u lp t io n
«
«C a lc u la t e  B W T  s trin g t a f o r  r e f e r e n c e  s t r in g  X  C a lc u la t e  a r r a y  c( •) a n d o (
•，•) fro m t a  C a lc u la t e  B W T  s trin g t a 、f o r  t h e  r e v e r s e  r e f e r e n c e
C a lc u la t e  arrayo、( •，•) fro m t a 、P r o c 05
d llr a 5s:
l —^H X A C T c /)H A I l c f f i  (5
Z )
n A L C U L A T u (
$)r e t u r n  I  麵REi
(勻_与| -
1 ,z,l
,
1
X 1 丨
1)
n A L C U L A T u ( w)
k <—1
T m  —1
for 、= 0 t o  |与| — 1 d o
A :T C ?〔5) +0、(^〔5 V I I ) +1
z
丨 c (^7〕) +0、(与〔5b  i H -h A : > Z  t h e n  ktl
T m  —1
D i i Y Z
^(5r -N v b if  z  <b s  t h e n  r e t u r n  0 if  ~ < 0 t h e n  r e t u r n  { I  A :, Z
二T 07个7 u  INEXREi(^k-rz-lv,z) f o r  e a c h  6"A,c,G,r} d o
T s ) +0(M -1) +1
z 丨c(6) +o (M )if  A :/A z  ths
7个7 U r E X  R
-e c u r  (^v l l r ^o  if  6 H ^s t h e n r ~N  c r x ^d E C U R  ( $? 一—一 , z, k, 1)e ls e
T 7 U  I _R E l w ,i —r z —l r l )
r e t u r n  7
»辂 1 BW
A W 辂穿ft®
B W A -M E M »{#^
賺笤函 2驾邻，»知啤滞
E  屮；1->^4。(1
) i M A »c l i l  著 r e a d s :
掛r -r s U F  练驾31卿涂#<i *^s ®^»c li l , FM D-index^E
掛丨->»$
画 B 4®^f r s f p ^s ^s s »c l i l  ,^E )t s  洚 _讁漭H I e ( M E M )^4m t #o o 0%s &铒。b w a ^
图2 B W A -MEM 算法流程
reads 不是一次性全部读人的，而是分块读取。分块
大小是一个可由用户指定的参数，默认为
100000000个碱基。这样可以避免reads 过大导致 无法放人内存。B W A 完成对当前分块的比对并写
回比对信息之后，再开始读人下一个分块。
(2) 查SMEM :BW A-M EM 算法采用经典的
seed -and-extend  范式，所谓 seed -and-extend  范式即
先到一段比较精确的匹配子段作为seed ，然后采 用如Smith-Waterman 等算法向两侧进行扩展，直到 达到比对条件的边界，得到的即为一个候选比对结
果。BW A-M EM 采用的seed 是“超级最长精确匹配 (super-maximal  exact  matches  ) ”，简称 “ SMEM  ”，其 定义为，读段序列中与参考基因组能够精确匹配的 最长的子序列，并且这个序列不能向两边进一步地 拓展，同时，这个子序列不能被其他SM EM 包含。
BW A-M EM 所采用的寻seed 的方法是在文献[4]
中提出的基于FM D -IN D E X 的方法，将在后续章节 中进一步介绍该算法。在某些情况下，最优的比对
结果可能并不包含任何一个“ SMEM ”，针对这种情 况，BW A-M EM 算法引人了 re-seeding 环节。如果得 到的SM EM 在in d e x 中出现了&次，并且长度过长 (默认参数下为28个bp )，则进人re-seeding 模块， 在re-seed 的过程中，条件加强为覆盖了上述SMEM 的中间碱基，并且在in d e x 中出现了& + 1次的最
SMEM 。通过这种方法能够避免由于SMEM 过长而
导致不到最优匹配的情况。
(3) 生成Chain ：在BW A-M EM 算法中定义了一 个C hain 的概念，由部分重叠在一起的seed (不存在 覆盖的seed )或距离足够近的seed 组成一条Chain 。
C hain 的意义是一方面将部分重叠在一起的seed 合
并后可以减少搜索空间，提高效率，另一方面可以允 许在最终比对中存在少量的错配和gap ，提高比对 质量。并不是每一个chain 在扩展后都能够得到一 个比对结果。为了提高算法效率，ch a in 生成后，还 会对这些chain 进行一次过滤操作。过滤的目的是
去除那些在更长的c h a in 中被大量覆盖住的、并且 质量得分及长度都远不如更长的c h a in 的那些
chain 。这里的chain 并不要求一定是精确的。
(4) 基于Smith-Waterman 算法拓展Chain ：在
一步中生成的C hain 的基础上，为了让整个读段能 够在参考序列上进行比配，需要对C hain 在两个方 向上进行拓展。在BW A-M EM 算法中，seed 拓展使 用的是基于Sm ith-W aterm an 的动态规划方法。
BW A-M EM 算法的seed 拓展方法与标准的Smith - Waterman 拓展算法主要有两点区别。第一，BWA -
M EM 加人了启发式剪枝，这个启发式剪枝算法能够
提高比对速度。第二，在seed 拓展过程中，BWA -
M EM 会记录每一个reads 末端的最优得分，如果这
个得分与局部最优得分的差值小于一个阈值，那么 即使局部最优得分更高，也不会被选取。BWA -
M EM 通过这种策略实现在局部最优与全局最优比 对的自动选择。(5) 生成并输出比对结果：经过上述步骤后以得到比对结果的元数据。接下来BW A 会根据这些 元数据，并结合用户指定的参数，来生成比对结果信 息，包括匹配结果、位置、注释等信息，输出比对结果。
1.3算法热点分析
使用GPU 加速BW A-M EM 算法，首先需要对算 法进行热点分析。BW A-M EM 算法适用于70b p 以 上的reads 比对，但是对于不同长度的reads ，算法的 热点有一定区别。将BW A-M EM 算法按模块分解， 并分别用100b p 和250b p 的reads 数据进行热点测 试。图3-4中可以看出，BW A-M EM 算法中有两个 最主要的热点模块，m em _chain 和chain 2aln 。其中
m em _chain 负责从reads 中出超级最长精确匹配
一
69
—
的作用是对生成的C hain进行双端拓展，利用改进 的Smith-Waterman算法，寻整个reads在参考序 列上的最佳匹配。
图4对长度250b p re a d s的热点对比
随着reads长度的增加，chain2a ln模块的时间 占比逐渐超过了mem_ c h a in模块。这是由于chain2a ln中使用的Smith-Waterman算法的时间复 杂度是0 (n2 )，而mem_chain模块的时间复杂度是 〇(n)，显然随着n的增长，chain2a ln模块运行时间 的增长速度要高于mem_chain模块。因此，对于较 长的reads，chain2a ln模块时间占比更大。
根据测试结果，BW A-M EM算法加速最重要的 就是对mem_chain模块和chain2a ln模块进行加速。2B W A-M E M的G P U加速设计与实现
2. 1整体设计
mem_chain和chain2a ln这两个模块是整个程序 的两个最重要的热点。如前文所述，mem_chain算法 主要完成SMEM的查和C hain的生成，而chain2aln 模块则负责将C hain的两端利用Smith-Waterman算 法进行拓展，到最佳比对位置。Smith-Waterman算 法的GPU加速已经有科研人员进行了深人的研究，
取得了很好的加速效果。而SMEM查和C hain的生成目前未见有GPU加速相关研究出现。因此重点 关注mem_chain模块的GPU加速。mem_chain模块 的整体流程如图5所示。
m em_chain模块中包含两个重要模块：mem_ collect_intv和 generate_chain。其中 mem_collect_in-tv模块从rea d中查最大精确匹配（MEM)，而gen-erate_chain模块将距离符合条件的mem构成chain。
由于整个算法流程要经过多个模块，为了能 够使程序在GPU上有更好的性能，设计了图6所 示的多kernel并行框架。G PU加速程序主要由三 一 70 —
图5m e m_chain模块算法流程
个kernel组成，分别对应前文中提到的三个热点模块。其中，基于Smith-W aterman的G PU加速算 法已经有很多相关研究，Ernst Joachim Houtgast 在中完成了对ch a in拓展环节的G PU加速，取 得了 3倍的加速效果。在本文中，主要测试GPU 对生成M E M和生成C hain模块的加速效果。
图6B W A-MEM多kernel并行框架
为了完成多kernel模式的移植，
需要对原程序
中的多线程模式进行重构。在原程序中，多线程分 支从程序完成数据读人后就开始执行，每个线程执 行整个比对流程。这种多线程模式不适用于GPU 设备上进行并行。因此将其优化为图7的多线程模 式。用更小的，多个kernel完成比对任务。
图7B W A-MEM重构后的多线程模式
2.2任务划分
任务划分涉及到两方面的工作，分别是G rid、B lock维度划分和序列与线程的划分。
对于G rid和B lo ck维度划分，根据C U D A编程 模型，G rid中的各个b lo ck会被分配到G PU的各个 流多处理器（stream multiprocessor，SM)中执行。在 维度选择时优先考虑B lo c k维度设计，G rid—般越 大越好。B lo ck维度设计需要根据实际计算中每个 S M中使用的寄存器、共享存储器数量与硬件每个 SM上可用的资源量来确定。通过编译工具可以观 察到mem_collect_intv这个ke rn e l占用了 66个寄存 器，496字节的常量存储器。对于generate_chain这 个k e rn e l占用了 46个寄存器，416字节常量存储器。
接着由硬件型号确定每个S M的资源，如每SM 中w arp总数上限、每S M上b lo c k总数上限、每SM 寄
存器数量、每S M上共享存储器数量。最后根据 二者计算b lo c k的维度。一般维度选择的原则是每 个b lo c k中的线程数量是32的整数倍，线程数量保 持在64 ~ 256之间。由于BW A-M EM算法对全局存 储器合并访问较困难，因此g rid和b lo c k维度对性 能的影响也更加复杂。再结合需要在全局存储器中 为每个线程预留足够的运算结果存储空间，选定b lo ck大小为64，g rid大小为64，均为一维。
接下来是reads与线程的划分。在比对过程中 reads具有天然的无关性，即不同reads之间完全独 立，只需要参考序列信息，不需要其他reads的信 息。因此采用最直观的处理方法，每个reads分配 一个线程进行运算。
2.3存储设计与优化
CUDA设备与主机端进行数据传输时，只能够对 支持位拷贝的数据进行传输。也就是说，对于包含指 针变量的结构体，指针的指针等数据结构，无法直接 在设备和主机端进行传输。因此需要对数据结构进 行调整，将需要在设备与主机端传输的数据选取出 来，将这些数据独立分配为一块连续的存储空间，加 快数据传输速度。以保存生成ch a in的数据结构为 例，优化后的chain存储模式分别如图8所示。
图8优化后的chain存储格式
G PU内部的6种存储器存储空间有限，因此需 要根据数据的数据量和访问特点来确定其存储位 置。BW A-M EM的ke rn e l中使用的数据表1所示。
表1B W A-MEM传输到GPU数据的存储方式变量含义长度/属性存储位置
s e q s查询序列查询序列长度，只读共享存储器b w t- >b w t参考序列B W T很大，只读纹理存储器
b w t bwt数据结构信息很小，只读常量存储器cu d a_o p t比对参数很小，只读常量存储器b w t- >s a参考序列后缀数组很大，只读纹理存储器global_chain生成的
c h a in较大，读写全局存储器g lo b a l_s e e
d s ch a in中的s
e e d s较大，读写全局存储器
一71
—

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。