研究報告

(1)

スーパーコンピュータ「京」における

融合遺伝子検出パイプラインの高速化

1 _{伊東聰，}

1 _{白石友一，}

2 _{島村徹平，}

1 _{千葉健一，}

1 _{宮野悟}

[email protected]

1

_{東京大学医科学研究所ヒトゲノム解析センター DNA情報解析分野}

2

_{名古屋大学大学院医学系研究科システム生物学分野}

平成26年度「京」における高速化ワークショップ

2014年12月19日

(2)

内容

• がん，DNAシークエンスとスーパーコンピュータ

• Genomon-fusion for 京の概要

– 京コンピュータへの移植

• CCLE RNA-seqの全計算

– 計算時間および解析結果などの概要

• blatのOpenMP化

• まとめと今後の課題

「京」における高速化ワークショップ

2 2014/12/19

(3)

がんと遺伝子異常

• がんは遺伝子異常の病気

– 加齢と共に発症率が増加

＝

後天的要素

による発症

– 日々の生活の中で徐々に蓄積

• 修復機構で治せなかったもの

– 遺伝子異常の種類

• 一塩基多型(SNP)

• 挿入/欠損

• コピー数変化

• 構造変異

• 融合遺伝子

• etc

Fig.1 BCR-ABL1 fusion gene

(http://watchingtheworldwakeup.blogspot.jp/200

9/07/strangeness-of-cancer.html)

DNA情報を読むことで，がんの原因

やメカニズムを解き明かせる！

(4)

ヒト17番染色体

約

8.1万文字のDNA情報

DNAの損傷認識と修復を行う遺伝

子。遺伝性乳がんとその変異との関

係がわかっている。

１ 81,189

DNA

Exon 1

Exon 2

Exon 24

Exon 3

Exon 23

約

8 万

文

字

の

領

域

に

エ

ク

ソ

ン

と

い

う

部

分

に

蛋

白

質

が

コ

ー

ド

RNA

へ

の

転

写

メ

カ

ニ

ズ

ム

RNA

mRNA

タ

ン

パ

ク

質

へ

の

翻

訳

エ

_ク

ソ

ン

部

分

が

切

り

出

さ

れ

る

ス

プ

ラ

イ

シ

ン

グ

と

い

う

プ

ロ

セ

ス

セントラルドグマ・アニメ

東京大学医科学研究所

宮野悟教授提供

Fig.2 Central dogma animation

(5)

DNAシークエンスのコスト推移

Fig.3 Trends of DNA sequencing cost

(http://www.genome.gov/sequencingcosts/)

(6)

DNAシークエンス解析の概要

• 次世代シークエンサーが出してくれるデータ

– RNA (20~30GB) < Exome < Whole genome (200~300GB)

– フォーマットは大体Fastq

– リード長：50bp～100bp

• DNAは細胞膜を壊して取り出す際に断片化されている

• bp: base pair＝塩基対，ATGCで表される文字一つ分のこと

「京」における高速化ワークショップ

Fig.4 Sample of Fastq data

@CCLE_253J_BV_RNA_08.CELL_1/1

CGGGGCTAAGAGACAGCAAATACACATGAACAGAAAGAAGAGGTCAAAGAAAAGGCTGACGGCAAGTTAACAAAAAGAAAAATGGAAAATAATACCCGTCC

+

@??DDDDDB?3AAEAGGF;AA9<2:::FDGB@??F??DGFFI=0/?DFFICB;8((..67;E89,5;>ADCA:?B5'5<5<<9@#################

@CCLE_253J_BV_RNA_08.CELL_614/1

GTGGGGGTGGTGTTAGTACCCCATCTTGTAGGTCTTGAGAGGCTCGGCTACCTCAGTGTGGAAGGTGGGCAGTTCTGGAATGGTGCCCGGGGCCGAGGGGG

+

CCCFFFFAFHFHHIJJIJJJJJJJJJJJIJJJEHIJIIJJJJIJJ6A?=?AE@@;A;@;@?5@>@,:?B################################

@CCLE_253J_BV_RNA_08.CELL_1227/1

GTGCGGGGTGGGCCCAGTGATATCAGCTGCCTGCTGTTCCCCAGATGTGCCAAGTGCATTCTTGTGTGCTTGCATCTCATGGAACGCCATTTCCCCAGACA

+

??=?DDDDHHA?BGADDHICGBGIGIHIIBFFEH97..7?ACC?B;C;ACAC>,>C;(;;55>:>C9A5;>@A4::+>:@C@:@9@99&8>>C:A?#####

@CCLE_253J_BV_RNA_08.CELL_1840/1

GGCAGAAGAGGGGCGGGGAGCTGTGTGCCCTAAGATCTCATTGCCTTTTTATGCCGATTAACATGCTTTTAGCCCCTACTGAGCTTATAGTTAACAGAAGT

+

C@@FFFFFGGGHHFIIIF>?8=A@ACBCCDCCCCCDDDDDADCCDDDDDD@DCCCBDBBBCCDDC:>CCDCDCDD58?BD<:@CCA@CD:>BA((:>ACDC

6 2014/12/19

(7)

DNAシークエンスとスーパーコンピューティング

• DNAデータの解析

– アラインメント（マッピング）

• 各リードが，DNA鎖のどの部分であるかを探すこと

– アノテーション/フィルタリング/etc

• DNAのどこに変異（異常）があるかを探す

• 我々のミッション：がんに関連するDNA変異を見つける！

– 遺伝子の不均一性

• 健常人と患者，臓器による変化，人種依存性，時間変化，…

数千万～数億リードに対し，そのポジションを同定す

る作業（アラインメント）が必要であり，その計算コス

トが非常に大きい．

多数検体，検体内多部位など，多くのシークエンス

データの統計的検定により得られる知見が重要．

(8)

Genomon-fusion

• “

Whole Transcriptome

シークエンスの結果である

FASTQ

ファイルのマッピング，

fusion gene

を検出し，

fusion gene

の候補の一覧を出力するソフトウェアです．”

http://genomon.hgc.jp/rna/

より抜粋

• オープン/フリーソフトウェアを利用したパイプライン

– Bowtie

– Blat

– CAP3

– fasta36

– SAMtools

– Picard

– bedtools

「京」における高速化ワークショップ

8 2014/12/19

(9)

概要と特徴

• シークエンスのアラインメントを丁寧に行う．

– 丁寧なアラインメントによりスプライシングを高感度に検出することで，スプラ

イシングの検出ミスによるアーティファクトを除く．

• ２段階のアラインメント

1. トランスクリプトームの配列に

bowtieでアラインメント

．

（その後ゲノム配列の座標に変換）

• 既知のスプライシング情報を用いる（既知のトランスクリプト配列に張り付ける）ことで，スプ

ライシングの検出を感度良く行うことができる！！

• アラインメントされた配列はゲノム配列の座標に変換．

2. １でアラインメントされなかったリードを，

blatを用いて

，ゲノム配列にアライ

ンメントを行う．

• blatにより，データベースに登録のないスプライシングも検出しつつアラインメントを行うこと

ができる．

• 計算時間はちょっとかかる（

→スパコンがないと厳しいか）

東京大学医科学研究所

白石友一氏提供

(10)

Genomon-fusionのフロー（アラインメント部分）

Fastq

split

fastq.aaa

bowtie

aligned.sam

unaligned.

fastq

blat

aligned2.sam

merged.sam

fastq.???

bowtie

aligned.sam

unaligned.

fastq

blat

aligned2.sam

merged.sam

Result.sam

「京」における高速化ワークショップ

10 グリッドエンジン

逐次プロセスをグリッドエン

ジンがまとめて一つのジョブ

としてシステムに投入する

仕組み．

2014/12/19

(11)

京コンピュータへの移植

• 移植の目的：多数検体を同時に処理する

• Genomon-fusionのHGCスーパーコンピュータでの運用

– 1検体のFastqファイルをsplitして並列処理

– グリッドエンジン（SGE/UGE）を利用したマルチジョブ

投入で実現

• 京コンピュータでは…

– MPI/OpenMPを用いた

並列プログラムのみ

受け付け

• シリアルプログラムは動くが，グリッドエンジンはない

– 独自のファイルシステム（ステージング）への対応

• ステージング：ログインノードと計算ノードのディスク間でのファイル授受

(12)

Genomon-Fusion for 京（アラインメント部分）

• GFKの大雑把な流れ

split

INPUT.0

aligned.sam

unaligned.

fastq

aligned2.sam

merged.sam

INPUT.?

aligned.sam

unaligned.

fastq

aligned2.sam

merged.sam

Result.sam

Fastq

Result.sam

_sam.0

bowtie

blat

bowtie

blat

前処理

(GFKpre.sh)

本処理

(GFKmain.sh)

全体を一つの

MPIプロセス化

後処理

(GFKpost.sh)

「京」における高速化ワークショップ

12 2014/12/19

(13)

CCLE RNA-seq 全計算

• CCLE

– Cancer Cell Line Encyclopedia

http://www.broadinstitute.org/ccle/home

• RNA-seq: 780 検体, 20TByte

• Calc. Summary:

– Total core time: 7,027,933,648 sec (976,102ノード時間)

• 分野1課題4の

2013年度割り当て分の12.7%

– # of core usage: 499,401

• 2 core/Node, メモリ使用量（8GB/proc）を確保するため

– 2 jobs=16 hours

(14)

blatのプロファイル

• テストケース：3ケース

– Case1: 平均値 (round-robin split)

– Case2: 最遅

– Case3: 最速

• Input data:

– Fastaフォーマット

• bowtieでマッピングされなかったリード

– ファイルサイズ（行数）

• Case1: 52,222

• Case2: 68,922

• Case3:

6,776

ー＞ケース3では，ほとんどのリードが

bowtieによりマッピングされている！

2014/12/19

「京」における高速化ワークショップ

14 Fig.7 Sample: CCLE_253J_BV_RNA_08

1File=5,000,000 lines (125k pair-reads)

(15)

blatのプロファイル

Timer start = main(0) Timer start = blat(0)

Timer start = blat(part1)(0)

Loaded 3095693983 letters in 25 sequences Timer stop = blat(part1)(0) 76.049922 Timer start = blat(part2)(0)

Timer start = blat(part2 region1)(0)

Timer stop = blat(part2 region1)(0) 0.000524 Timer start = blat(part2 region2)(0)

Timer stop = blat(part2 region23)(0) 0.000766 Timer stop = blat(part2 region2)(0) 272.296611 Timer start = blat(part2 region3)(0)

Timer start = searchOneIndex(1)

Timer start = searchOneIndex region 4(1)

Timer stop = searchOneIndex region 4(26111) 6213.238720

Searched 2637211 bases in 26111 sequences Timer stop = searchOneIndex(1) 6213.241863 Timer stop = blat(part2 region3)(0) 6213.243220 Timer start = blat(part2 region4)(0)

Timer stop = blat(part2 region4)(0) 0.000504 Timer stop = blat(part2)(0) 6485.543722 Timer stop = blat(0) 6561.628977

Timer stop = main(0) 6561.679652 Finished! nest = 0

Timer start = main(0) Timer start = blat(0) type = static, mod=0 num_threads = 1

Timer start = blat(part1)(0)

Loaded 3095693983 letters in 25 sequences Timer stop = blat(part1)(0) 11.941816 Timer start = blat(part2)(0)

Timer stop = blat(part2 region23)(0) 0.000004 Timer stop = blat(part2 region2)(0) 256.131548 Timer start = blat(part2 region3)(0)

Timer start = searchOneIndex(1)

Timer start = searchOneIndex region 4(1) Timer start = searchOneIndex region 41(1)

Timer stop = searchOneIndex region 41(26111) 0.036332 Timer start = searchOneIndex region 42(1)

Timer stop = searchOneIndex region 42(0) 1651.952223

Timer stop = searchOneIndex region 4(0) 1651.988864 Searched 2637211 bases in 26111 sequences

Timer stop = searchOneIndex(1) 1651.992685 Timer stop = blat(part2 region3)(0) 1651.992701 Timer start = blat(part2 region4)(0)

Timer stop = blat(part2 region4)(0) 0.000002 Timer stop = blat(part2)(0) 1908.124283 Timer stop = blat(0) 1920.315536

Timer stop = main(0) 1920.316067 Finished! nest = 0

Fig.8 Calc. time of blat for case1(averaged)

Left: Fujitsu FX10, Right Shirokane1

(16)

OpenMPによるスレッド並列化

2014/12/19

「京」における高速化ワークショップ

16

void searchOneIndex(int fileCount, char *files[], struct genoFind *gf, char *outName,

boolean isProt, struct hash *maskHash, FILE *outFile, boolean showStatus) /* Search all sequences in all files against single genoFind index. */

{ int i;

char *fileName; int count = 0;

long long totalSize = 0; gfOutputHead(gvo, outFile); for (i=0; i<fileCount; ++i)

{

fileName = files[i]; if (nibIsFile(fileName))

{

struct dnaSeq *seq; if (isProt)

errAbort("%s: Can't use .nib files with -prot or d=prot option¥n", fileName);

seq = nibLoadAllMasked(NIB_MASK_MIXED, fileName); freez(&seq->name);

seq->name = cloneString(fileName);

searchOneMaskTrim(seq, isProt, gf, outFile, maskHash, &totalSize, &count); freeDnaSeq(&seq);

}

else if (twoBitIsSpec(fileName)) {

struct twoBitSpec *tbs = twoBitSpecNew(fileName); struct twoBitFile *tbf = twoBitOpen(tbs->fileName); if (isProt)

errAbort("%s is a two bit file, which doesn't work for proteins.", fileName);

if (tbs->seqs != NULL) {

struct twoBitSeqSpec *ss = NULL;

for (ss = tbs->seqs; ss != NULL; ss = ss->next) {

struct dnaSeq *seq = twoBitReadSeqFrag(tbf, ss->name, ss->start, ss->end); searchOneMaskTrim(seq, isProt, gf, outFile,

maskHash, &totalSize, &count); dnaSeqFree(&seq);

} } else

{

struct twoBitIndex *index = NULL;

for (index = tbf->indexList; index != NULL; index = index->next) {

struct dnaSeq *seq = twoBitReadSeqFrag(tbf, index->name, 0, 0);

searchOneMaskTrim(seq, isProt, gf, outFile, maskHash, &totalSize, &count); dnaSeqFree(&seq); } } twoBitClose(&tbf); } else {

static struct dnaSeq seq;

struct lineFile *lf = lineFileOpen(fileName, TRUE);

while (faMixedSpeedReadNext(lf, &seq.dna, &seq.size, &seq.name)) {

searchOneMaskTrim(&seq, isProt, gf, outFile, maskHash, &totalSize, &count); } lineFileClose(&lf); } } carefulClose(&outFile); if (showStatus)

printf("Searched %lld bases in %d sequences¥n", totalSize, count); }

Fig.9 Source code of searchOneIndex

(blat/blat.c)

Region 1

Region 2

Region 3

Region 4

(17)

• faMixedSpeedReadNext

– タイプ：読み込み

– 内容

• インプットファイルからの

リード(k-mer)読み込み

• searchOneMaskTrim

– タイプ：メイン計算

– 内容

• リード(k-mer)のアラインメント，

最も重い部分

• lineFileClose

– タイプ：ファイル書き込み

– 内容

• アラインメント結果をファイルに書き出す

OpenMPによるスレッド並列化

searchOneMaskTrim(&seq, isProt, gf, outFile, maskHash, &totalSize, &count); }

lineFileClose(&lf);

(18)

• searchOneMaskTrimのOpenMP化

OpenMPによるスレッド並列化

2014/12/19

「京」における高速化ワークショップ

18

searchOneMaskTrim(&seq, isProt, gf, outFile, maskHash, &totalSize, &count); }

lineFileClose(&lf);

#pragma omp parallel private(i, ii, thread_num)

{

thread_num = omp_get_thread_num();

#pragma omp for // Mail loop

for( ii = 0; ii < lcount; ii++ ) {

searchOneMaskTrim(&seq[ii], isProt, gf, outFile[thread_num], maskHash, &totalSize, &count, thread_num); } // End of for

} // End of parallel region

Fig.11 Source code of region 4 with OpenMP

スレッド数取得

(19)

• OpenMP化のポイント

– インプットファイルの一括読み込み

• whileループのforループ化のために必要

– アウトプットファイルの分割

• スレッド処理の独立性を担保するため．全処理終了後にファイルをマージする．

– cat blat.psl.* > blat.psl

– 行ごとにデータが独立しているので，順番を問わない点を活用．

– gvo

• ファイル出力に関するグローバル変数

– 関数(searchOneIndex)内での変更がないので，全スレッドでコピーを持たせる

– メモリ確保

• 関連する配列のメモリはスレッドごとに確保・解放する．

OpenMPによるスレッド並列化

(20)

OpenMPによるスレッド並列化

2014/12/19

「京」における高速化ワークショップ

20

Case1 Case2 Case3

# of threads 1 2 4 1 2 4 1 2 4 Elapse time 6811.76 4719.19 2923.65 17137.47 9459.88 6788.9 495.67 386.66 387.09 Time (main) 6380.9 4375.55 2491.25 16791.53 9116.41 6441.93 65.27 43.08 40.13 Speed-up (main) 1 1.46 2.56 1 1.84 2.6 1 1.51 1.62 0 2000 4000 6000 8000 10000 12000 14000 16000 18000

T

ime

(s

ec

)

Calculation time of blat on FX10

(21)

まとめ

• がん，DNAシークエンスとスーパーコンピューティング

• Genomon-fusion for 京の概要

• CCLE RNA-seq 全計算

– 約800検体をおよそ97万ノード時間で計算終了．

– アラインメント部分（blat）が非常に計算コストがかかっていた．

• blatのOpenMP化

– メモリ使用量の関係から遊んでいたコアの有効活用．

– 4コアで2.5倍程度の高速化を実現した．

• 今後の課題

– 実用上の問題として，データマネジメントが見えてきた．

研究報告

スーパーコンピュータ「京」における

融合遺伝子検出パイプラインの高速化

1

伊東 聰，

1

白石友一，

2

島村徹平，

1

千葉健一，

1

宮野 悟

[email protected]

東京大学 医科学研究所 ヒトゲノム解析センター DNA情報解析分野

名古屋大学大学院 医学系研究科 システム生物学分野

平成26年度「京」における高速化ワークショップ

2014年12月19日

内容

• がん，DNAシークエンスとスーパーコンピュータ

• Genomon-fusion for 京の概要

– 京コンピュータへの移植

• CCLE RNA-seqの全計算

– 計算時間および解析結果などの概要

• blatのOpenMP化

• まとめと今後の課題

「京」における高速化ワークショップ

2

2014/12/19

がんと遺伝子異常

• がんは遺伝子異常の病気

– 加齢と共に発症率が増加

＝

後天的要素

による発症

– 日々の生活の中で徐々に蓄積

• 修復機構で治せなかったもの

– 遺伝子異常の種類

• 一塩基多型(SNP)

• 挿入/欠損

• コピー数変化

• 構造変異

• 融合遺伝子

• etc

Fig.1 BCR-ABL1 fusion gene

(http://watchingtheworldwakeup.blogspot.jp/200

9/07/strangeness-of-cancer.html)

DNA情報を読むことで，がんの原因

やメカニズムを解き明かせる！

ヒト17番染色体

約

8.1万文字のDNA情報

DNAの損傷認識と修復を行う遺伝

子。遺伝性乳がんとその変異との関

係がわかっている。

１

81,189

DNA

Exon 1

Exon 2

Exon 24

Exon 3

Exon 23

約

8

万

文

字

の

領

域

に

エ

ク

ソ

ン

と

い

う

部

_{伊東聰，}

_{白石友一，}

_{島村徹平，}

_{千葉健一，}

_{宮野悟}

_{東京大学医科学研究所ヒトゲノム解析センター DNA情報解析分野}

_{名古屋大学大学院医学系研究科システム生物学分野}

_ク