.
.

計測自動制御学会東北支部第 282回研究集会 (2013.
07.
17)

資料番号 282‑
10

ショートリード配列マッピンク.のためのFPGAアクセラレータの開発
DevelopmentofanFPGA‑BasedAcceleratorf
orDNA
Short‑ReadMappingUsingBurrows‑WheelerAlignment
0張山昌論ヘハシタムトウマラウイシディスーリヤヘ亀山充隆*

oNIasanoriHaxiyama*，HasithaMuthumalaWaidyasooriyaヘJ¥1ichitakaKameyama
本*東北大学大学院情報科学研究科
TohokuUniversity*GraduateSchoolo
fInformationSciences，
キーワード: ヘテロジニアスアーキテクチャ (Heterogeneousarchitecture)，高性能計算 (high‑performancesupercomputing)，FPGA(FPGA)
computing)，スーパーコンピューティング (連絡先: 干 980・8579 仙台市青葉区荒巻字青葉 ι6‑05
東北大学大学院情報科学研究科

張山昌論， Tel.
:(022)
7
95・7
153，Fax.
:(022)
263‑9167，E‑mail:hariyama0ecei.tohoku.ac・jp

1. はじめにトウェアツールでのショートリードのマッピン

次世代シーケンサーでは，数百万個以上のショートリード (DNAを分断した短い DNA鎖)を生成することができる.そのような，膨大な数のショートリードをマッピングすることが DNA解析の重要な基本処理となっている.ショートリー
ドのマッピング用ツールとしては，PCクラスタ上でのソフトウェアとして， Maq1)，Bowtie2)，BWA3)などのが存在する.ショートリー
ドのマッピングは，図 1に示すように，数 10から数 100個の塩基もつショートリードを， 1塩基だけが異なる場合 (SNP)，塩基が挿入される
場合 (Inse凶ion)，および，塩基が欠損する場合(Deletion)といった場合を考慮して，参照配列と呼ばれる既に解析済みの配列のどの箇所に対応するかを決定する.それにより，個人毎の遺伝子の違いを特定することが可能となる.ソフグの問題は，その膨大な処理時間である.大規
模な PCクラスタでの並列処理を用いたとしても数日 ‑‑1週間程度かかる.
この問題を解決するために， FPGA
(fieldprかgrammablegatearray)を用いたアクセラレー
タを提案する 4). FPGAでは回路構造を，コンフィグレーションと呼ばれるプログラムにより変更することにより，アプリケーション処理に適した専用プロセッサを構築できる.近年のFPGAは，半導体集積回路技術の進展により，数百万以上のプログラマプルロジック，DSPユニット，内蔵メモリ， DDR2および DD
R3などの外部メモリインタフェースを備えており，高性能な専用プロセッサ開発にも使えるようになってきた.実装するアルゴリズムとしては，ソフトウェアとして現在最も高速な手法としてしら

‑1‑

・・0

園園圏圏圏
囲園圏圏圏

Fig
. 1 ショートリードのマッピング
れる Burrows‑WheelerAlignment(BWA)3)をアルゴリズムを説明するために，図 3に示す用いる. BWAはショートリード毎の並列処理例を考える.図 3(a)に示すような参照配列 Xも可能でため， FPGA実装にも適する.提案のと，ショートリード W を考える.マッピング

FPGAペースアクセラレータでの処理性能は，処理の入力は， X の BWT，X の反転 (X')のソフトウェア実装と比較して 10倍以上高いといB¥;YT，生起配列 0(.，，)・0'(.，・)および C(.
)でう結果を得た.

ある.配列 C(.
)および X の Rotationを，そ(a)および 3(b)に示す.ソーとしれぞれ，図 3た結果を図 3(c)に示す.図 3(c)に示すように，2. BWAアルゴリズム
簡単に BWAアルゴリズム 3)について説明す

BWTBおよび生起配列 0(.，.)を得る.同様のやり方で， X'に対しでも B'および 0'(.，.)をる. BWAアルゴリズムは，文献 5)で提案されて得る.紙面の都合により，図 2のアルゴリズム

いるマッチング方法を用いている .αWは文字中の， Calculated処理の説明は省略する，ここ列α{，W}を表すとする.もし文字列 W が，文で， Calculated処理は W 中のミスマッチ数の字列 X の部分文字列で，かつ k(αW)壬l(αW)(.
)を与える処理である(より詳しくは文下限 D

で あ れ は 文 字 列 αWもX の部分文字列である献3)を参照).という性質がある.変数 kおよび lは，式 (1)お図 4 はショートリード W のマッピングを示よび (2)で与えられ，それぞれ X の suffixarrayしている.処理 InexR
βcur(l‑V，i，z，k，l)は， ミ(SA)インターパルの下限と上限を与える.

スマッチの下限 D(i)の計算の後に，呼び出される .Wの探索位置，許容されるミスマッチ数

k(αW)=C(α)+O(αル(α11')‑1)+1 (1)l(αW)=C(α)+O(α，l(αW))(2)(SNP，挿入，欠損の合計)，suffixarrayの下限，および suffixarrayの上限は，それぞれk，および tで与えられる.図 4に， InexRecur
処理の実行の様子を示す.ラベル

辞書式順序において αより小さい文字の数は，C(α)で与えられる.文字列 B[O，i]において αの生起する数は O(α，i)で与えられる.文字列 Xの BWT配列は B で与えられる.文字列 B の先頭から t個の文字を含む部分文字列は，B[O，i]で与えられる.
図 2 は BWAアルゴリズムを用いたショートリードのマッピングの方法 3)を示している.

i，z，(i，z，k，l)"は処理 InexRecur(W，i，ムk，l)"を表す.マッピン
グが終了した後に，参照配列中のショートリー
ドの位置， SNP，挿入および欠損(Indel)に関する情報を得る.この例では，図 5において SNPおよびIndelをミスマッチとして考慮した結果を得ている. 1個の挿入を持つ結果は SAイン
ターパル (5，5)を持っている(図 4).図3(c)に‑2‑。

2.3 4156
Posl
首cn 0:
11:

IDput: BWTstringBforreferences凶ngXArrayC(.
)叫 0(.，.
)fromB
BWTs
凶ngB'for曲erev，ぽ路 ofrefe.
renceX
Arrayぴ(.
，.
)fromB'ou
旬ut:Su節x‑arrayin
回v
als

CITGIAGI$IC

proωdures:

Calculated(W)b噸 n1ateD(山 内lowerbound山Dm叫諸rofmi
組 踊 凶I
民泊 W￨ωω

回叫(a)ReferencesequenceX，
short‑
readW andC(α)ofXInexRecur，w( i，z，k，1)b噸 nlt z<D(i)then￨ 関 知m ゆ2

TIGAIG$ICIC

3

GIAGIS CICIT

4AIGSICCITG

56

GI$ CICTIGA。

Occurrence(.，.
)
rrat0‑ c2

sccT G A G 010 011 。
4T G 010 012 。
。AcGcTsGcc
A G s 110 012 。

3C T G A G $IC 110 112

SAPosi
首on0:，1 :234 56d1=ゆ1=1UInex.
Recur(W，i‑1，z‑1，k，1)
foreachbE{A，C，G，T}dokb=C(的+O(b，k‑l)+1lb=C(b)+0(b，l)
ffkb~ lb血m1=IUIn飢Recur(W，i，z‑1，kbllb)I
fb=Wド]th阻
I1=IUIn
蹴Recur(W，i‑1，ztkb，lb}

SICICITIG生直(b)Rotat
ionofXBWT

聞 dlti<O也mIreturn [k，~阻

CICTIGAIGI$。

1GT

65

G $IC cT G A

53G AG62

TG A G

4$fA

111 112。
。

scC T 111 11211
scC 111 21211
固‑圃else(c)Sortingr
邸 ult

lI=ILJcur(w， i‑1，z‑1，kbtlb)

Ine.成e佃dv
lappingexample
Fig.3 Tendendm加 m 1

2.
1 FPGAペースアクセラレータのアー

end

main，w( z)
b噸 a

キテクチャ

￨ωc
叫a
凶 m

ショートリード配列マッピング問題において

Ine
点 目ur
(w， i，z，k，1)

end

は，参照ゲノム配列は通常同じものが用いられFig.2

る.そのため，本稿では， Burrows‑Wheeler変BWAアルゴリズム

換と生起配列をあらかじめオフラインで生成しSAインターパル (5，5)は参照文字列ておき，マッピング処理の際には FPGAに転送

X における位置 3に対応している.したがって，するだけにする.図 6に FPGAを用いたアクセ

W は位置 3にマッピングされることになる(図ラレータアーキテクチャを示す. FPGAアクセ5(a)). 1個の SNPを持つ結果は， SAインターラレータは，シンプルな 128個の PEと2個のパル (6，6)を持っている.したがって ，W は位D DR2メモりから構成されている.各 PEでは，

置 2にマッピングされることになる(図 5(b))同異なるショートリードに対するマッピング処理

様に， 1個の欠損はは S
Aインターパル (3，3)をが並列に行われる持ち，位置 1にマッピングされる(図 5(c))図 4は，それぞれ DDR2メモリと FPGA内のメモに示すように， Ine
xRe
cur処理は，深さ優先探リに転送される.

従うと，索，幅優先探索，そのハイプリッドにより行う

生起配列とショートリード

図 7に PEの構成を示す. PEは， 32ピット加算器，比較器，レジスタから構成され，式 (1)，ことができる.

および (2)を計算する. 1回の I
nxRecur処理︒
︒

i=2I i=1=‑1i=O I i

naMilloγ51(0，0， 1.1)、~(-1. ・1，5，5) in.
.
.
.
rtion
{
1.
1l，L 1 1 1 S N P目
唱 刷 出(‑
1，0，5，5)1insertion(1，0，5，5)

0，6，611SNP(‑
1，(1.
0，2，ヨ}6)1~I~~Ie~?n(‑1・，1，6，6)
nsertion(1.
0，4，5)
0，6，6) ((1，，‑1，3，3)
0Matd
司

園 圃 .

Fig.4

ーー争

SNP

山ー田争

Insertion1・1，0，
ヨ311deletion+
OeletionFig，6ーー

アクセラレータアーキテクチャ

ショートリード阿Fのマッピング(a)SAinterv
al:(5，5)，posi‑(b)SAinterval:(6，6)，pωi‑tion:3tion:2
SNP@ X(2)Insertion@ X(3)

Fig.7

PEの構成

データを消去するさらに，マッピング処理が行われている閑， CPUはショートリードデータをnterval:(3，3)，position:1(c)SAi
Deletion@ X(2}

少しずつ FPGAに転送している. FPGA内の一つのミスマッチを許した場合のマッピング結果

PEは，複雑な浮動小数点演算器や複雑な制御を終了した後に，レジスタファイルから新しいドのマッピングに特化した固定小数点演算のたデータか読み込まれる.図 2で示したように，各めの基づくシンプルな演算器である.そのため，

InxRecur処理では，新しい I
nxRecur処理が呼数多くの PEを 1個の FPGA上に搭載することばれる.この再帰的呼び出しのパラメータはレが可能となっている.

Fig. 5回路を有する CPUとは異なり，ショートリー

ジスタファイルに記憶され，再帰的呼び出しの図8，に，複数の PEによる DDR2メモリへ制御に使われる. PE内の ADD/SUBユニットのアクセスのタイムチャートを示す.本アーキ

は式 (1)，および (2)を計算するために使われテクチャでは，る.比較器と制御パスが条件分岐の制御を行う.

リへのアクセス性能で決まる.ショートリード

古いショートリードに対する処理が終了した後マッピングでのメモリアクセスはランダムアクに，新しいショートリードが PEに入力される.

セスであり，そのアクセスは入力されたショートFPGA内の出力メモリに記憶されたマッピング0GB/s
リードに依存する. CPUは理論的には 5

結果は， CPUから読み出される.出力メモリが以上の高い帯域を有しているが，ランダムアクオーバフローしないように， CPUは，出力バッ

セスに対しては実行メモリ帯域は大幅に低くなファから結果を読み出し，読み出しが終了した

る.一方， FPGAでは処理に特化したアドレス‑4‑トータルの性能は， DDR2メモ

~Tabl
e1 1個の PEのリソース使用量
リソース名

PEn (使用量(使用量)/(全利用可能リソース量)

x100%
.
20
LUTs
8
4
4 0Registers
1
01
5 0.
2
4.25kB 0
Onchipmemory 1.
05)笹野( ‑ ・ ・Fig.8 メモリアクセスのタイムチャート
Bandwidth(GB/s)
理論値￨ランダムアクセス時の値

CPU(i下3960x) I51.2 I
1.06FPGA (EP4SGXI
6.
4x2I
5.
15x2
530KH40C2)

域は大幅に低下することが分かる.ショートリー

Fig.9 構築したショートリード配列マッピン

ド配列マッピングにおけるメモリアクセスはラグシステム
生成田路により，ランダムアクセスであってもンダムであるため，実行メモリ帯域が大幅に低高いバンド幅を得ることができる.

下し，それにより CPUペースシステムの性能が大幅に低下する.一方， FPGAベースのシス
テムでは，対象とする処理に特化した無駄のな

3. 評価

いメモリアドレス生成ユニットを設計すること

評価には， DE4ボード 6)を用いた. DE4ボーができるので，ランダムアクセスであっても，ドには， Alter
aEP4SGX530KH40C2FPGAと5.15GB/s という理論的な上限値に近い高いバ4GBDDR2‑SDRAJ¥.
1が 2個搭載されている.図

ンド幅を得ることができる. FPGAには， 2個9に構築したシステムを示す. このシステムでの DDR2メモリを接続できるため， 10GB/s以orei
7‑396OxCPUを搭載したマザーボー
は， c上のバンド幅を得ることができる.その結果，

ドと， DE4ボードが PC1expressポートを経由

FPGAベースシステムでは， CPUベースシスで接続されている.アクセラレータの動作周波テムと比較して 10倍以上の性能を達成できる.

数は， 100tv
lHzである.表 1に， 1個の PEのさらに， 3個の FPGAボードを一つのホスト PC

リソース使用量を示す.この表から分かるように PC1expressポートを接続できるため，それに， PEサイズはレジスタ数により制限されてに比例して性能を向上できる.また，そのようおり， 1個の FPGAに 400個以上の PEを搭載なシステムを複数用いたり， FPGAボード同士できることが分かる.

を低レイテンシな専用 1/0により接続することにより，さらなる性能向上が期待できる.

表 2は， CPUペースシステムと， FPGAベースシステムのメモリ帯域の比較を示している.

我々は， 3000文字程度の参照配列に対して 100

CPUは理論的には高いバンド幅を持っている個のショートリードのマッピングが可能なシスが，ランダムアクセスに対しての実行メモリ帯テムの構築した.ソフトウェアでの処理と比較

υF円して， 10倍程度の高速化が達成された.現在，より大規模な参照配列に対応できるシステムを

6)htt
p:
jjw
ww
.
altera.
comjeducationjuniv
jmaterialsjboardsjde4
junv・de4‑board.
html

開発中である.

4. まとめ
本稿では， BWAアルゴリズムを高速化するための FPGAベースアクセラレータを提案した.
ソフトウェアでの処理と比較して 10倍程度の高速化が達成された.今回の開発で用いた FPGA
は，中規模の FPGAで高いグレードのものではない.現在利用可能な最新の FPGA(Altera社 StratixVなど)を用いることにより，さらなる高速化・高並列化が可能できることが期待される.

謝辞
本研究の一部は文部科学省科研費 12020735の助成を受けたものです.

参考文献1)H.L
i， Maq: Mappingandassemblywithqualiti回
ぺ http:jjmaq.sourceforge.netj，
2008.
. Tra
pnell，M. Pop，and2)B. Langmead，Cltr
afastandmemor
y‑efficient
S.Salzberg，Ualignmen
to
fs
hor
tdnas
equen
ce
stot
hehu
‑
mangenome
ヘGenomeBiology，Vol.10，No.
3，p.
R25，2009.
3) HengL
iandRichar
dDurbin，"Fastandac‑curateshor
treadalignmen
twithBurrows‑
ぺ Bioinfomatics，Vo.
125，Wheelertransfor
mNo.
14，pp.
1754‑1
760，2009.
4)HasithaMuthumalaWaidyasooriya，Masanoriy，出naandMichitakaKameyama，吋m‑
H釘 ie‑
plementation o
f a Custom Hardwar
AcceleratorforShort‑read Mapping Using"7 35もh An‑
Burrows‑
Wheeler Alignment，nualInternationalConferenceoft
he IEEE
Engineeringi
nMedicin
eandB
iologySociety(EJ¥'lBS)，
7pp.
651‑
654，?
2013.
5)P.Ferragin
aandG.:
Manzini，Opportunisticdatastructureswithapplications"，Proc.of
4
1stSymp.onFoun
dation
sofComputerSci‑pp.
390‑398，
2009.
en
ce，

‑6‑