研究報告

Size: px

Start display at page:

Download "研究報告"

たみえひのと
6 years ago
Views:

1 スーパーコンピュータ京における融合遺伝子検出パイプラインの高速化 1 伊東聰, 1 白石友一, 2 島村徹平, 1 千葉健一, 1 宮野悟 sito@hgc.jp 1 東京大学医科学研究所ヒトゲノム解析センター DNA 情報解析分野 2 名古屋大学大学院医学系研究科システム生物学分野平成 26 年度京における高速化ワークショップ 2014 年 12 月 19 日秋葉原 UDX4 階 UDX GALLERY NEXT-2

2 内容がん,DNAシークエンスとスーパーコンピュータ Genomon-fusion for 京の概要京コンピュータへの移植 CCLE RNA-seq の全計算計算時間および解析結果などの概要 blatのopenmp 化まとめと今後の課題京における高速化ワークショップ 2

3 がんと遺伝子異常がんは遺伝子異常の病気加齢と共に発症率が増加 = 後天的要素による発症日々の生活の中で徐々に蓄積修復機構で治せなかったもの遺伝子異常の種類一塩基多型 (SNP) 挿入 / 欠損コピー数変化構造変異融合遺伝子 etc Fig.1 BCR-ABL1 fusion gene ( 9/07/strangeness-of-cancer.html) DNA 情報を読むことで, がんの原因やメカニズムを解き明かせる! 京における高速化ワークショップ 3

4 セントラルドグマアニメヒト 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

5 DNA シークエンスのコスト推移 Fig.3 Trends of DNA sequencing cost ( 京における高速化ワークショップ 5

6 DNA シークエンス解析の概要次世代シークエンサーが出してくれるデータ RNA (20~30GB) < Exome < Whole genome (200~300GB) フォーマットは大体 Fastq リード長 :50bp~100bp DNAは細胞膜を壊して取り出す際に断片化されている bp: base pair= GTGGGGGTGGTGTTAGTACCCCATCTTGTAGGTCTTGAGAGGCTCGGCTACCTCAGTGTGGAAGGTGGGCAGTTCTGGAATGGTGCCCGGGGCCGAGGGGG + GTGCGGGGTGGGCCCAGTGATATCAGCTGCCTGCTGTTCCCCAGATGTGCCAAGTGCATTCTTGTGTGCTTGCATCTCATGGAACGCCATTTCCCCAGACA GGCAGAAGAGGGGCGGGGAGCTGTGTGCCCTAAGATCTCATTGCCTTTTTATGCCGATTAACATGCTTTTAGCCCCTACTGAGCTTATAGTTAACAGAAGT + C@@FFFFFGGGHHFIIIF>?8=A@ACBCCDCCCCCDDDDDADCCDDDDDD@DCCCBDBBBCCDDC:>CCDCDCDD58?BD<:@CCA@CD:>BA((:>ACDC 京における高速化ワークショップ Fig.4 Sample of Fastq data 6

7 DNA シークエンスとスーパーコンピューティング DNA データの解析アラインメント ( マッピング ) 各リードが,DNA 鎖のどの部分であるかを探すことアノテーション / フィルタリング /etc DNA のどこに変異 ( 異常 ) があるかを探す数千万 ~ 数億リードに対し, そのポジションを同定する作業 ( アラインメント ) が必要であり, その計算コストが非常に大きい. 我々のミッション : がんに関連する DNA 変異を見つける! 遺伝子の不均一性健常人と患者, 臓器による変化, 人種依存性, 時間変化, 多数検体, 検体内多部位など, 多くのシークエンスデータの統計的検定により得られる知見が重要. 京における高速化ワークショップ 7

8 Genomon-fusion Whole Transcriptome シークエンスの結果である FASTQ ファイルのマッピング,fusion gene を検出し,fusion gene の候補の一覧を出力するソフトウェアです. より抜粋オープン / フリーソフトウェアを利用したパイプライン Bowtie Blat CAP3 fasta36 SAMtools Picard bedtools 京における高速化ワークショップ 8

9 概要と特徴東京大学医科学研究所白石友一氏提供シークエンスのアラインメントを丁寧に行う. 丁寧なアラインメントによりスプライシングを高感度に検出することで, スプライシングの検出ミスによるアーティファクトを除く. 2 段階のアラインメント 1. トランスクリプトームの配列に bowtie でアラインメント. ( その後ゲノム配列の座標に変換 ) 既知のスプライシング情報を用いる ( 既知のトランスクリプト配列に張り付ける ) ことで, スプライシングの検出を感度良く行うことができる!! アラインメントされた配列はゲノム配列の座標に変換でアラインメントされなかったリードを,blat を用いて, ゲノム配列にアラインメントを行う. blat により, データベースに登録のないスプライシングも検出しつつアラインメントを行うことができる. 計算時間はちょっとかかる ( スパコンがないと厳しいか ) 京における高速化ワークショップ 9

10 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

11 京コンピュータへの移植移植の目的 : 多数検体を同時に処理する Genomon-fusionのHGCスーパーコンピュータでの運用 1 検体のFastqファイルをsplitして並列処理グリッドエンジン (SGE/UGE) を利用したマルチジョブ投入で実現京コンピュータでは MPI/OpenMP を用いた並列プログラムのみ受け付けシリアルプログラムは動くが, グリッドエンジンはない独自のファイルシステム ( ステージング ) への対応ステージング : ログインノードと計算ノードのディスク間でのファイル授受京における高速化ワークショップ 11

12 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

13 CCLE RNA-seq 全計算 CCLE Cancer Cell Line Encyclopedia 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

14 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

15 blat のプロファイル Timer start = main(0) Timer start = blat(0) Timer start = blat(part1)(0) Loaded letters in 25 sequences Timer stop = blat(part1)(0) Timer start = blat(part2)(0) Timer start = blat(part2 region1)(0) Timer stop = blat(part2 region1)(0) Timer start = blat(part2 region2)(0) Timer start = blat(part2 region21)(0) Timer stop = blat(part2 region21)(0) Timer start = blat(part2 region22)(0) Timer stop = blat(part2 region22)(0) Timer start = blat(part2 region23)(0) Timer stop = blat(part2 region23)(0) Timer stop = blat(part2 region2)(0) Timer start = blat(part2 region3)(0) Timer start = searchoneindex(1) Timer start = searchoneindex region 4(1) Timer stop = searchoneindex region 4(26111) Searched bases in sequences Timer stop = searchoneindex(1) Timer stop = blat(part2 region3)(0) Timer start = blat(part2 region4)(0) Timer stop = blat(part2 region4)(0) Timer stop = blat(part2)(0) Timer stop = blat(0) Timer stop = main(0) 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 letters in 25 sequences Timer stop = blat(part1)(0) Timer start = blat(part2)(0) Timer start = blat(part2 region1)(0) Timer stop = blat(part2 region1)(0) Timer start = blat(part2 region2)(0) Timer start = blat(part2 region21)(0) Timer stop = blat(part2 region21)(0) Timer start = blat(part2 region22)(0) Timer stop = blat(part2 region22)(0) Timer start = blat(part2 region23)(0) Timer stop = blat(part2 region23)(0) Timer stop = blat(part2 region2)(0) 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) Timer start = searchoneindex region 42(1) Timer stop = searchoneindex region 42(0) Timer stop = searchoneindex region 4(0) Searched bases in sequences Timer stop = searchoneindex(1) Timer stop = blat(part2 region3)(0) Timer start = blat(part2 region4)(0) Timer stop = blat(part2 region4)(0) Timer stop = blat(part2)(0) Timer stop = blat(0) Timer stop = main(0) Finished! nest = 0 京における高速化ワークショップ 15

16 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

17 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

18 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

19 OpenMP によるスレッド並列化 OpenMP 化のポイントインプットファイルの一括読み込み while ループの for ループ化のために必要アウトプットファイルの分割 gvo スレッド処理の独立性を担保するため. 全処理終了後にファイルをマージする. cat blat.psl.* > blat.psl 行ごとにデータが独立しているので, 順番を問わない点を活用. ファイル出力に関するグローバル変数関数 (searchoneindex) 内での変更がないので, 全スレッドでコピーを持たせるメモリ確保関連する配列のメモリはスレッドごとに確保解放する. 京における高速化ワークショップ 19

20 Time (sec) OpenMP によるスレッド並列化 Calculation time of blat on FX Case1 Case2 Case3 # of threads Elapse time Time (main) Speed-up (main) Fig.12 Result of OpenMP blat 京における高速化ワークショップ 20

21 まとめがん,DNAシークエンスとスーパーコンピューティング Genomon-fusion for 京の概要 CCLE RNA-seq 全計算約 800 検体をおよそ97 万ノード時間で計算終了. アラインメント部分 (blat) が非常に計算コストがかかっていた. blat の OpenMP 化メモリ使用量の関係から遊んでいたコアの有効活用. 4コアで2.5 倍程度の高速化を実現した. 今後の課題実用上の問題として, データマネジメントが見えてきた. Sparc64VIIIfxの整数演算性能についてポスト京への期待京における高速化ワークショップ 21

22 謝辞 HPCI 戦略プログラム分野 1 予測する生命科学医療および創薬基盤課題 4 大規模生命データ解析 (hp140230) 東京大学医科学研究所ヒトゲノム解析センタースーパーコンピュータ ( 財 ) 高度情報科学技術研究機構山本秀喜様, 野口孝明様京における高速化ワークショップ 22

nagasaki_GMT2015_key09

nagasaki_GMT2015_key09 Workflow Variant Calling 03 長崎は遺伝研大量遺伝情報研究室の所属です国立遺伝学研究所生命情報研究センター 3F 2F 欧州EBIと米国NCBIと密接に協力しながら DDBJ/EMBL/GenBank国際塩基配列データベースを構築しています私たちは塩基配列登録を支援するシステムづくり登録データを活用するシステムづくり高速シーケンス配列の情報解析を行なっています