スーパーコンピュータ 京 における 融合遺伝子検出パイプラインの高速化 1 伊東聰, 1 白石友一, 2 島村徹平, 1 千葉健一, 1 宮野悟 sito@hgc.jp 1 東京大学医科学研究所ヒトゲノム解析センター DNA 情報解析分野 2 名古屋大学大学院医学系研究科システム生物学分野 平成 26 年度 京 における高速化ワークショップ 2014 年 12 月 19 日秋葉原 UDX4 階 UDX GALLERY NEXT-2
内容 がん,DNAシークエンスとスーパーコンピュータ Genomon-fusion for 京の概要 京コンピュータへの移植 CCLE RNA-seq の全計算 計算時間および解析結果などの概要 blatのopenmp 化 まとめと今後の課題 京 における高速化ワークショップ 2
がんと遺伝子異常 がんは遺伝子異常の病気 加齢と共に発症率が増加 = 後天的要素による発症 日々の生活の中で徐々に蓄積 修復機構で治せなかったもの 遺伝子異常の種類 一塩基多型 (SNP) 挿入 / 欠損 コピー数変化 構造変異 融合遺伝子 etc Fig.1 BCR-ABL1 fusion gene (http://watchingtheworldwakeup.blogspot.jp/200 9/07/strangeness-of-cancer.html) DNA 情報を読むことで, がんの原因やメカニズムを解き明かせる! 京 における高速化ワークショップ 3
セントラルドグマ アニメ ヒト 17 番染色体 タンパク質への翻訳 mrna スエプクソラインシ部ン分グがと切いりう出プさロれセるス RNA DNA Exon 1 Exon 2 Exon 3 RNA への転写メカニズム Exon 23 Exon 24 約 8.1 万文字の DNA 情報 1 DNA の損傷認識と修復を行う遺伝子 遺伝性乳がんとその変異との関係がわかっている 81,189 京 における高速化ワークショップ 4 約 8 万文字の領域にエクソンという部分に蛋白質がコード 東京大学医科学研究所宮野悟教授提供 Fig.2 Central dogma animation
DNA シークエンスのコスト推移 Fig.3 Trends of DNA sequencing cost (http://www.genome.gov/sequencingcosts/) 京 における高速化ワークショップ 5
DNA シークエンス解析の概要 次世代シークエンサーが出してくれるデータ RNA (20~30GB) < Exome < Whole genome (200~300GB) フォーマットは大体 Fastq リード長 :50bp~100bp DNAは細胞膜を壊して取り出す際に断片化されている bp: base pair= 塩基対,ATGCで表される文字一つ分のこと @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 京 における高速化ワークショップ Fig.4 Sample of Fastq data 6
DNA シークエンスとスーパーコンピューティング DNA データの解析 アラインメント ( マッピング ) 各リードが,DNA 鎖のどの部分であるかを探すこと アノテーション / フィルタリング /etc DNA のどこに変異 ( 異常 ) があるかを探す 数千万 ~ 数億リードに対し, そのポジションを同定する作業 ( アラインメント ) が必要であり, その計算コストが非常に大きい. 我々のミッション : がんに関連する DNA 変異を見つける! 遺伝子の不均一性 健常人と患者, 臓器による変化, 人種依存性, 時間変化, 多数検体, 検体内多部位など, 多くのシークエンスデータの統計的検定により得られる知見が重要. 京 における高速化ワークショップ 7
Genomon-fusion Whole Transcriptome シークエンスの結果である FASTQ ファイルのマッピング,fusion gene を検出し,fusion gene の候補の一覧を出力するソフトウェアです. http://genomon.hgc.jp/rna/ より抜粋 オープン / フリーソフトウェアを利用したパイプライン Bowtie Blat CAP3 fasta36 SAMtools Picard bedtools 京 における高速化ワークショップ 8
概要と特徴 東京大学医科学研究所白石友一氏提供 シークエンスのアラインメントを丁寧に行う. 丁寧なアラインメントによりスプライシングを高感度に検出することで, スプライシングの検出ミスによるアーティファクトを除く. 2 段階のアラインメント 1. トランスクリプトームの配列に bowtie でアラインメント. ( その後ゲノム配列の座標に変換 ) 既知のスプライシング情報を用いる ( 既知のトランスクリプト配列に張り付ける ) ことで, スプライシングの検出を感度良く行うことができる!! アラインメントされた配列はゲノム配列の座標に変換. 2. 1 でアラインメントされなかったリードを,blat を用いて, ゲノム配列にアラインメントを行う. blat により, データベースに登録のないスプライシングも検出しつつアラインメントを行うことができる. 計算時間はちょっとかかる ( スパコンがないと厳しいか ) 京 における高速化ワークショップ 9
Genomon-fusion のフロー ( アラインメント部分 ) 逐次プロセスをグリッドエンジンがまとめて一つのジョブとしてシステムに投入する仕組み. Fastq split fastq.aaa bowtie グリッドエンジン fastq.??? bowtie aligned.sam unaligned. fastq aligned.sam unaligned. fastq blat blat aligned2.sam aligned2.sam merged.sam merged.sam Result.sam 京 における高速化ワークショップ 10
京コンピュータへの移植 移植の目的 : 多数検体を同時に処理する Genomon-fusionのHGCスーパーコンピュータでの運用 1 検体のFastqファイルをsplitして並列処理 グリッドエンジン (SGE/UGE) を利用したマルチジョブ投入で実現 京コンピュータでは MPI/OpenMP を用いた並列プログラムのみ受け付け シリアルプログラムは動くが, グリッドエンジンはない 独自のファイルシステム ( ステージング ) への対応 ステージング : ログインノードと計算ノードのディスク間でのファイル授受 京 における高速化ワークショップ 11
Genomon-Fusion for 京 ( アラインメント部分 ) GFK の大雑把な流れ Fastq Fastq Fastq Fastq split 前処理 (GFKpre.sh) INPUT.0 bowtie 本処理 (GFKmain.sh) 全体を一つの MPI プロセス化 INPUT.? bowtie aligned.sam unaligned. fastq unaligned. fastq aligned.sam blat blat aligned2.sam aligned2.sam merged.sam merged.sam Result.sam Result.sam Result.sam sam.0 後処理 (GFKpost.sh) 京 における高速化ワークショップ 12
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 実際には, プリプロセスの都合上,100 検体ずつ 8job で処理 京 における高速化ワークショップ 13
blat のプロファイル テストケース :3 ケース Case1: 平均値 (round-robin split) Case2: 最遅 Case3: 最速 Input data: Fasta フォーマット bowtie でマッピングされなかったリード ファイルサイズ ( 行数 ) Case1: 52,222 Case2: 68,922 Case3: 6,776 ー > ケース 3 では, ほとんどのリードが bowtie によりマッピングされている! Fig.7 Sample: CCLE_253J_BV_RNA_08 1File=5,000,000 lines (125k pair-reads) 京 における高速化ワークショップ 14
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 start = blat(part2 region21)(0) Timer stop = blat(part2 region21)(0) 9.632589 Timer start = blat(part2 region22)(0) Timer stop = blat(part2 region22)(0) 262.455965 Timer start = blat(part2 region23)(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 Fig.8 Calc. time of blat for case1(averaged) Left: Fujitsu FX10, Right Shirokane1 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 start = blat(part2 region1)(0) Timer stop = blat(part2 region1)(0) 0.000003 Timer start = blat(part2 region2)(0) Timer start = blat(part2 region21)(0) Timer stop = blat(part2 region21)(0) 3.424628 Timer start = blat(part2 region22)(0) Timer stop = blat(part2 region22)(0) 252.706811 Timer start = blat(part2 region23)(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 京 における高速化ワークショップ 15
OpenMP によるスレッド並列化 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); Region 1 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, Region maskhash, &totalsize, 2 &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, Region isprot, gf, 3outFile, 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); } Region 4 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) Target block is region4! 京 における高速化ワークショップ 16
OpenMP によるスレッド並列化 famixedspeedreadnext タイプ : 読み込み 内容 インプットファイルからのリード (k-mer) 読み込み searchonemasktrim タイプ : メイン計算 内容 リード (k-mer) のアラインメント, 最も重い部分 linefileclose タイプ : ファイル書き込み 内容 アラインメント結果をファイルに書き出す 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); Fig.10 Source code of region 4 京 における高速化ワークショップ 17
OpenMP によるスレッド並列化 searchonemasktrim の OpenMP 化 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); while -> for ループ #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 スレッド数取得 スレッドごとにファイル出力 京 における高速化ワークショップ 18
OpenMP によるスレッド並列化 OpenMP 化のポイント インプットファイルの一括読み込み while ループの for ループ化のために必要 アウトプットファイルの分割 gvo スレッド処理の独立性を担保するため. 全処理終了後にファイルをマージする. cat blat.psl.* > blat.psl 行ごとにデータが独立しているので, 順番を問わない点を活用. ファイル出力に関するグローバル変数 関数 (searchoneindex) 内での変更がないので, 全スレッドでコピーを持たせる メモリ確保 関連する配列のメモリはスレッドごとに確保 解放する. 京 における高速化ワークショップ 19
Time (sec) OpenMP によるスレッド並列化 18000 Calculation time of blat on FX10 16000 14000 12000 10000 8000 6000 4000 2000 0 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 Fig.12 Result of OpenMP blat 京 における高速化ワークショップ 20
まとめ がん,DNAシークエンスとスーパーコンピューティング Genomon-fusion for 京の概要 CCLE RNA-seq 全計算 約 800 検体をおよそ97 万ノード時間で計算終了. アラインメント部分 (blat) が非常に計算コストがかかっていた. blat の OpenMP 化 メモリ使用量の関係から遊んでいたコアの有効活用. 4コアで2.5 倍程度の高速化を実現した. 今後の課題 実用上の問題として, データマネジメントが見えてきた. Sparc64VIIIfxの整数演算性能について ポスト京への期待 京 における高速化ワークショップ 21
謝辞 HPCI 戦略プログラム分野 1 予測する生命科学 医療および創薬基盤 課題 4 大規模生命データ解析 (hp140230) 東京大学医科学研究所ヒトゲノム解析センタースーパーコンピュータ ( 財 ) 高度情報科学技術研究機構山本秀喜様, 野口孝明様 京 における高速化ワークショップ 22