生物情報科学科 情報基礎実験 ３年生冬学期 10/26 – 11/5 （森下研究室担当）
本演習では次世代シーケンサーデータを解析する基本的なプログラムを作成します。
(0) C/C++ および Java のプログラミング環境をインストール
l 環境のインストール（斉藤先生作成） http://www.xerial.org/wiki/lecture/2009/Setup
l C++プログラムのデバッグ法（斉藤先生作成） http://www.xerial.org/wiki/lecture/2010/CppDebug
l 使用言語は C/C++ を推奨。 ただしJava でも構いません。
(1) Suffix Array の実装
l
Induced Sorting を利用して線形時間アルゴリズムを実装
参考論文：Nong, Zhang, and Chan.
Two
Efficient Algorithms for Linear Suffix Array Construction.
Appendix に Cプログラムあります。自作が難しいようでしたら、このプログラムを解読して動かしてください。提出用のプログラムにはできるだけコメントを入れてください。なおこのプログラムの動きを例題を使って解説した資料は以下にあります。
supplementary_to_InducedSortingAlgorithm.pdf
l
チェック用のアルゴリズム（たとえば Larsson-Sadakane’s O(n log n)-time
algorithm）を実装して、同一の答えが出ることを確認してください。自分で書くことをお勧めしますが、大変な場合は、C++ で書いた Larsson-Sadakane’s O(n log n)-time algorithmの プログラムは以下にあります。
Larsson_Sadakane.cpp
Visual C++ 等のC++ プログラムを認識するツールで開いてください。無限ループする場合があるのでバグを修正しました（2010/11/10）。
l テスト用に、自分でランダムな配列を生成し、長さを 10^4, 10^5, 10^6, 10^7, 10^8 と伸ばしたとき、計算時間が長さに対してほぼ線形に増えてゆくことを確認してください。
l
実データとして、以下のページから最新のヒトゲノムデータ（2009年2月版 hg19はhuman
genomeの19版目を意味する）を入手します。
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/
たとえば chr1.fa.gz は、１番染色体(chr1) の塩基配列データを fasta 形式（最初の行に “>” からはじまるコメントあり）で格納したファイルで gzip で圧縮してあります。なお、 chr *** random.fa.gz および churn**** の名前をもつデータは染色体上での位置が確定していないデータなので、今回は利用しないでください。
データを読み込むには、有田先生の演習のときに教わった以下のページを参考にしてください。
http://www.metabolomics.jp/wiki/Aritalab:Lecture/Programming/Cpp/Genome
註：ゲノム中の N および大文字と小文字の扱い方
ゲノム配列中には大文字の A, C, G, T (コード領域)および小文字の a, c, g, t （コード領域以外）が現れていますが、今回の演習では区別しないので、すべて大文字に変換してください。その他、ゲノム配列には、A, C, G, T 以外に N という文字列が現れます。N は塩基を確定できなかったことを意味しています。たとえばテロメアやセントロメアのように繰返し配列が連続している領域では、解読することが諦められ N が続いています。
N の扱い方ですが、A, C, G, T に 1, 2, 3, 4 の番号を割り当てて、 N には 5 を割り当てる方法が自然です。また、多少変則的ですが、N に A, C, G, T を割り当てることで、4文字を 2 bit 表現して記憶領域を圧縮する方式もあります。主記憶や二次記憶装置が高価で潤沢に使えなかった時代に愛用されました。しかしN に A,C,G,Tを割り当てると、その位置に問合せ配列が誤って写像されることが起こります。しかしこのような悪影響はある程度避けられることもわかっています。ゲノム中には poly-A (A が連続して現れる配列)が多く存在します。この性質を利用して、 N に A を割り当てます。するともしpoly-Aを含む問合せ配列が問い合わされても、元からある多数の poly-A に写像されてしまい、位置を確定することはできません。したがって N に A を割り当てても、問合せ配列の位置を確定できない状況には変化が無いと言えるからです。もうひとつ愛用されているのは N に A,C,G,T をランダムに割り当てる考え方です。問合せ配列がある程度の長さ（長さ 20 以上）の時には、ランダムに生成された場所に誤って当たる確率は非常に低くなるため、間違いは起こりにくいからです。4文字に圧縮したい人は、自分の好きな方式で N を扱ってください。
l 染色体もしくはヒトゲノム全体で suffix array を構築します。
21もしくは22番染色体（ヒト常染色体では最も短いので）を対象に実験してください。suffix array が正しく構築できたことを、induced sorting と Larsson-Sadakane 双方で同じ結果が得られたということで確認できたということをレポートしてくれれば十分です。なお、ヒトゲノム全体を処理するのはオプション課題です。大きな主記憶が必要です。workstation (メモリ16GB) を利用してください。
(2) Burrows-Wheeler Transform (BWT) の実装 ミスマッチを許した問合せの実装
l suffix array から Burrows-Wheeler Transform を構築し、対象とするゲノム配列中に出現する問合せ配列の位置を枚挙するプログラムを作成します。長さn の配列（ヒトゲノム等）に対して、長さ k (= 10〜15) の問合せ配列をランダムに生成して、問合せ時間の平均時間を求め、k (および n) に関して線形比例するかどうか確認します。
l
Burrows-Wheeler
Transform のデータ構造である
Occ の主記憶容量を抑えるために、L 個毎に Occ の行を保存して、on-the-fly で他の行を補完する方法を実装します。L の値を 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024 ぐらいまで変化させたとき、記憶領域の変化と、長さ k (= 10〜15 等) のランダム配列に対する平均計算時間がどのように変化するか観察してください。記憶領域消費量の追跡には、Windows では Process Explorer が便利です。
http://technet.microsoft.com/ja-jp/sysinternals/bb896653.aspx
l 下記の論文を参考にして、ミスマッチおよびギャップを許して問合せ配列をゲノム上にアラインメントするプログラムを作成します。より多くのミスマッチおよびギャップを許容するに従って、速度が変化する様子を分析してください。
参考論文： Heng Li and Richard Durbin.
Fast
and accurate short read alignment with Burrows–Wheeler transform.
Bioinformatics 2009 25(14):1754-1760.
この論文にはいくつかの誤りがありますので、以下の資料で補います。
supplementary_to_Li_Durbin_Algorithm.pdf
(3) Chaining の実装と大規模ゲノムの比較
l ２つのゲノム中 G1, G2で対応する領域を探し出すために、G1 中の長さ k の塩基配列を、BWT
を使ってG2中にアラインメントし、その結果を Chaining
するプログラムを作成します。Chaining Algorithm は生物情報ソフトウエア論の講義ノートおよび配布資料を参照してください。
chaining.pdf (講義ノート)
k の値やミスマッチおよびギャップの許容範囲を工夫しないと、アラインメントを現実的な時間で探すのが難しくなります。まず、k の長さを 10〜20 程度の範囲で動かして、完全マッチするアラインメントを高速に見つけてみましょう。進化的に近い哺乳類内のヒト、マウス、ラットの配列比較だと、多くのアラインメントが見つかるはずですので、それらを chaining してみてください。また、k = 10〜20 のとき、chain の検出力がどのように変化し、計算時間がどのように増大するかを観察してみましょう。
G1 から長さ k の配列を取出すとき、1塩基毎にずらして kの長さの配列を抽出できますし、数を減らすなら重なり合わないようにして取り出す(non-overlapping)こともできます。後者のようにしても chaining には大きな影響がないかもしれません。
今回の演習では範囲外ですが、哺乳類と魚類を比較する場合は、k = 10〜20 でかなりミスマッチを許容しないとアラインメントは見つかってこないでしょう。たとえば
k = 20でミスマッチ数を 4 まで許容すると検出力は上がります。しかし計算時間も増大します。たとえば、いまわれわれはシーラカンスとチキンのゲノムを比較しているのですが（シーラカンスゲノムは解読されていないのですが、ゲノムが解読された一番近い種がチキンです）、このような工夫をしないと相同性の高い領域が見つかりません。時間がかかるので並列化して探索します。さらにヒトとショウジョウバエぐらい進化的に離れると、短いk で多くのミスマッチを許すなどして検出力を高める工夫が必要になります。
l 平衡木を実装する部分は、STL を使ってコードを書きましょう。
ü STLの使い方（斉藤先生作成）http://www.xerial.org/wiki/lecture/2010/STL
ü STL を使った Chaining のプログラム例（ただし方針だけで肝心なアルゴリズムの実装は伏せてあります）
chaining_sample_ver0.cpp
ü このプログラムで使っているテストデータ testdata_for_chaining.pdf
ü STL set, multimap, map のメソッドを解説したページ http://www.cppreference.com/wiki/start
l
例題１ 「ヒト、マウス、ラットのX染色体の比較」
ヒト、マウス、ラットゲノムのX染色体を比較し、下記の参考論文中の図６a「X染色体の比較図」に示す図を描画してみてください。
参考論文： Rat Genome Sequencing Project Consortium.
Genome
sequence of the Brown Norway rat yields insights into mammalian evolution.
Nature 428, 493-521 (1 April 2004)
図６ X染色体の比較
データは以下のサイトからダウンロードしてください。
http://hgdownload.cse.ucsc.edu/downloads.html
l
例題２ 「脊椎動物祖先での２回の全ゲノム重複の痕跡を探す」
例題1はかなり大きな染色体を比較するので大変かもしれません。代わりに、もう少し問題のサイズが小さいこちらの例を解いてもよいです。脊椎動物祖先での２回のゲノム重複は1970年に日本人遺伝学者 Susumu Ohno が予想し、多くの議論を巻き起こした説ですが、2005-7 年にかけてようやく実証されました。
http://genome.cshlp.org/content/17/9/1254/F4.large.jpg （２回全ゲノム重複後の脊椎動物の染色体進化）
ゲノム重複が２回起こったので、遺伝子が４倍にコピーされ、あらたな機能を獲得した遺伝子がうまれ、脊椎動物の繁栄をもたらしたと考えられています。以下の図は、HOX遺伝子群が４つの染色体に分かれて存在していることを示した図です。
two_round_whole_genome_duplications.pdf
HOX A グループはヒト７番染色体のおおよそ27,120,000-27,250,000塩基の位置（ヒトゲノムhg19版での位置）に広がっています。一方、HOX B, HOX C, HOX D のグループは各々 17, 12, 2 番染色体に存在しています。HOX A, B, C, D グループは脊椎動物の祖先でコピーされたので、その後の塩基変異によりだいぶ変化していますが、ある程度は類似しています。この類似度を頼りに HOX A グループの領域 chr7:27,120,000-27,250,000 と17, 12, 2 番染色体を比較し、chaining で HOX B, C, Dに対応する領域を見つけられるか？ というのが問題です。
ただし、長さ k の配列の完全マッチでは見つからないでしょう。ある程度ミスマッチを許容して探すことになります。またExon は類似していても、intron や遺伝子間領域は変異により類似性は失われています。ですからゲノム領域をそのまま chaining するのは、すこし荒っぽい探し方で、計算時間もかかります。そこでexon
領域だけを比較して複製された領域の対応関係を探してみます。たとえば HOXA1と HOXB1の遺伝子配列を chaining で比較してみると、 HOXA2 と HOXB1 の間よりは近いことがわかると思います。
HOX以外にも一般に４つの類似した領域がヒトゲノム内には見つかります。複製された領域であるという情報は遺伝子の機能を推定するのに役立っています。
l chaining の時間を軽減するには、ポリ A のようにゲノム上に多数出現する配列を問い合わせることを実行前にあらかじめ避けておく、もしくは実行時にさまざまな場所にマップされた問合せのアラインメントを使わない、等の工夫が必要になります。
レポートは課題毎に送ってください。〆切は次の演習がはじまる11月8日に一応しますが、11月末まで待ちます。11/10, 17, 24 の午後は理３号館３０９号室で質問に対応できます。
レポート提出先および質問等