2018年 解答

テスト対策と言うことで2018年の過去問（講義（Web: link）に載っている）を解いたので、その解答をまとめました。

問題1 Chaining

Chaining

• 文字列 \(\mathrm{Z}\) の開始位置 \(sx\) から終了位置 \(ex-1\) までの連続した部分文字列を、\(\mathrm{Z}[sx,ex)\) と表現する。
• ２つの文字列 \(\mathrm{P,Q}\) のなかで、類似した部分文字列\(\mathrm{P}[sx,ex),\mathrm{Q}[sy,ey)\)の組を、類似度を表すスコアとともに以下のように表現し、アラインメントと呼ぶ。
  $$\mathrm{A}=\left(\mathrm{P}[sx,ex),\mathrm{Q}[sy,ey),\text{score}\right)$$
  と表現し、アラインメントと呼ぶ。
• ２つのアラインメント \(\mathrm{A,B}\) 間で、以下が満たされている場合 \(\mathrm{A}<\mathrm{B}\) と定義する。（半開区間で定義しているので、\(\mathrm{A}\) と \(\mathrm{B}\) のアラインメントを重ねることができる。）
  • \(\mathrm{A}.ex\leq \mathrm{B}.sx\)
  • \(\mathrm{A}.ey\leq \mathrm{B}.sy\)
• アラインメント \(A_1<A_2<\ldots<A_k\) を鎖(chain)と呼ぶ。ここで、以下の二つを定義
  • 鎖の長さを \(k\)
  • 鎖のスコアを \(\sum_i\mathrm{A}_i.\text{score}\)

目的

動的計画法を用いて、最大スコアの鎖を求める。

アルゴリズム

#=== Initialization steps ===
X=[]
for A in Alignments:
    X.append((A.sx, A))
    X.append((A.ex, A))

# If x1 is equal to x2, they are sorted by whether it is "start point" or "end point".
X = sorted(X, key=lambda x:(x[0], -x[1].sx))

Y=[]

#=== Repeat steps ===
while(len(X)!=0):
    x,B = X.pop(0)

    #=== Case1 ===
    if x == B.ex:
        C_flags = [C.ey<=B.ey and C.score>=B.score for _,C in Y]
        # If alignment "C" which is better than alignment "B" is not in Y.
        if sum(C_flags)==0:
            Y.append((B.ey, B))
            Y = sorted(Y, key=lambda y:y[0])
            # If alignemd "D" is worse than alignment "B" remove "D" from Y.
            Y = [(x,D) for x,D in Y if not (B.ey<=D.ey and B.score>D.score)]

            #=== Print the Y's transition ===
            for y,A in Y:
                print("{}{}{}".format(A.score, A.name, y))
            print("-" * 20)

    #=== Case2 ===
    elif x == B.sx:
        # This flag indicates whether there is alignment "C" which alignment "B" can connect to.
        C_flags = [(C.ey, C) for _,C in Y if (C.ex<=B.sx and C.ey<=B.sy)]
        if C_flags:
            _,nearestC = sorted(C_flags, key=lambda c:c[0])[-1]
            newB = Alignment(B.name, nearestC.sx, nearestC.sy, B.ex, B.ey, nearestC.score+B.score)
            X = [(B.ex, newB) if A.name==B.name else (x,A) for x,A in X]

1

アラインメント \(\mathrm{B}1\) の開始位置 \(\mathrm{B}1.sx\) とアラインメント \(\mathrm{B}2\) の終了位置 \(\mathrm{B}2.ex\) が等しい時、②→①の順番だとスコアが高い終了位置のアラインメントが見落とされる可能性が出てしまう。

コメント

アルゴリズム（設問）の説明が悪いと感じた。 \(X\) の先頭を取り出す際（\(X\) をソートする際）に、同じ \(x\) であれば開始位置（ステップ②）のものより終了位置（ステップ①）のものを先に取り出すべきであり、ここではそのことに言及していると思われる。

2

ステップ②で、\(Y\) 内部から \(\mathrm{C}<\mathrm{B}\) を満たすもののうち \(\mathrm{B}.sy\) に \(ey\) が最も近いものを選択するため、スコアがより良く \(ey\) の小さい（＝\(\mathrm{B}.sy\) から遠い）\(\mathrm{C}\) が見過ごされてしまう可能性が出てしまう。

3

値の検索・挿入・削除を \(O(\log n)\) で計算する平衡木（2-3木、赤黒木）を使ってアラインメントのリスト \(Y\) を実装すれば、最悪計算量 \(O(n\log n)\) で最大スコアのchainが計算できる。

問題2 DNA修飾情報分析

1

真核生物の DNA は、核内において蛋白質複合体であるクロマチンの状態で存在する。クロマチンはヌクレオソームから構成され、4種類のヒストンタンパク質H2A,H2B,H3,H4がそれぞれ2つずつ合わさった8量体のコア・ヒストンに、DNAが巻きついた構造を取る。

ヒストンがメチル化やアセチル化などの化学的な修飾を受けると、クロマチンは

• DNAの凝集密度が低く転写因子が結合しやすく転写が活発に行われるユークロマチン
• DNAが凝集して転写因子が結合しにくいヘテロクロマチン

にその形状が変化する。

エンハンサーは遺伝子の転写量を増加させる作用をもつDNA領域のことで、遺伝子活性化因子と結合することで遺伝子の転写量を大幅に増大させる。

2

このアルゴリズムは、各記号の意味を正しく抑えることが重要である。

目的

実数値の列 \(\left\{L_i|i=1,2,\ldots,n\right\}\) 上の互いに交わらない区間の集合 \(C = \left\{[a_1,b_1],\ldots,[a_k,b_k]\right\}\) の中で、重み \(w(C) = \sum_{i\in I_j\in C}L_i\) を最大化するものを \(O(n)\) で計算するアルゴリズムの構築。

制約条件

この時、各区間の開始点 \(a_t\) と終点 \(b_t\) は、以下の不等式を満たす。

1. \(m_1\leq b_t-a_t+1\)：最終区間以外の各区間の長さは \(m_1\) 以上離れている。（\(t=1,\ldots,k-1\)）
2. \(m_0\leq a_t-b_{t-1}+1\)：隣接する区間は \(m_0\) 以上離れている。
3. \(a_1=1\) または \(a_1\geq m_0+1\)

※注意すべき点は、上記の赤字部分である。これらには生物学的な意味は全くなく、アルゴリズム設計を単純にするための条件である。

ここで、条件1を満たす区間集合を、最終区間 \([a_k,b_k]\) の状態で以下の2通りに場合分けする。

• \(C^0_{i,m}\)
  • \(b_k!=i\) つまり、位置 \(i\) で修了していない。
  • \(i\) が最終区間の終点 \(b_k\) から \(m\) 以上離れている。（\(m\leq i-b_k+1\)）
• \(C^1_{i,m}\)
  • \(b_k==i\) つまり、位置 \(i\) で修了している。
  • 最終区間の長さが \(m\) 以上ある。（\(m\leq b_k-a_k+1\)）

動的計画法によってこの問題を解くために、以下の4つの変数を導入する。

したがって、以下の漸化式を計算することで、重み \(w(C)\) を最大化することが可能である。

初期条件

漸化式

この時、\(w(\left\{[i-m_1+2,i]\right\})\) の計算が \(O(1)\) であるとすれば、各漸化式を計算するコストは \(O(1)\) これを \(i=[1,n]\) についてそれぞれ考えるから、計算量は \(O(n)\)

トレースバック \(\max\left\{W_{\text{long}}^1(n), W_{\text{short}}^0(n)\right\}\) から始めてトレースバック。

問題3 DNA3次元折りたたみ構造

参考論文：Comprehensive mapping of long range interactions reveals folding principles of the human genome

Hi-C解析とは

ゲノム配列そのものからは知ることができない染色体の立体構造を、ゲノム配列のシーケンシングによって明らかにする解析のこと。DNAの空間的な近接性は、制限酵素処理された断片間でのライゲーションの生じやすさで測れるという考え方が基になっている。

1. 実際に細胞の中で近接しているDNAとタンパク質の複合体をホルムアルデヒド処理することで固める。
2. 制限酵素を用いて断片化する。
3. 端にビオチン化ヌクレオチドを付加する。
4. ライゲーション処理をする。
5. 濃縮してビオチンプルダウンする。（ビオチンの付いたもの、つまりキメラ化したものだけを選択的に抽出したいから。）
6. ペアエンドでシーケンスする。（∵Forward, Reverseのリードがリファレンスゲノムへマッピングされた位置を調べ、それらのゲノム状の領域がもともと空間的に近接していた、と解釈する。）

※ビオチンの付いた部分であり、ホルムアルデヒドで固められた部分そのものをシークエンスしている訳ではない。制限酵素処理された時の配列の長さにも依存してしまう。
※精度の高いリファレンスゲノムがあることが前提。

1

定義

この時、\(C_{ij}\) は実対称行列なので、ある直交行列 \(P\) で対角化できる。

「固有値最大の固有ベクトル \(P_i\)」、主成分分析の言葉で言い換えれば「共分散行列 \(C_{ij}\) の第一主成分ベクトル」の正/負の値で分けられた2つの大域的な構造がA/B compartmentである。

生物学的な意義

同じコンパートメントに属している場合（AとBの場合）の方が、そうでない場合よりも、領域のスペースが近くなる傾向がある。また、上記論文（link）においては、コンパートメントAと遺伝子発現量や各ヒストンマーク（活性型のH3K36や抑制型のH3K27）のメチル化が正の相関があることを示している。

つまり、大域的に見れば「コンパートメント単位で遺伝子やエンハンサーの相互作用が制限されている」と考えることができる。

2

定義

※ 以下の変数 \(A,B\) は上記のABコンパートメントとは関係がない。

生物学的な意義

• TAD(Topologically Associated Domain)は、各コンパートメントの内のさらに局所的な構造である。
• TAD間の境界にはCTCFというタンパク質が結合していることが多い。コヒーシンがCTCF依存的に結合してインシュレーター（区切り壁）の役割を果たしていると考えられているため、TADごとにループ構造が（いくつか）あり、構造的に区切られていると考えることもできる。
• それぞれの細胞分化に必要な遺伝子群を、分化ステージや細胞系列に合わせて狭い領域にまとめることは難しい。このため全く異なる細胞で発現する遺伝子が隣接して存在することはゲノムでは普通に見られる。したがって、一つのTADの発現が、隣のTADにある遺伝子の発現に影響を及ぼさないよう構造化されているのであれば、それは特に初期胚での確立過程において極めて重要な役割を果たす。
• 各TADの間には、種間でよく保存された、CTCF分子の結合配列、house keeping遺伝子、tRNA遺伝子、そしてSINEトランスポゾンなどが集まった境界領域が存在している。

3

距離行列 \(D\) とグラム行列 \(G\) の関係は以下で表される。なお、距離行列は実験的にはHi-Cで観測する。

ここで、グラム行列の計算には、\(D\) の値に加えて \(\left|\vec{V}_{i}\right|^{2},\left|\vec{V}_{j}\right|^{2}\) の値が必要となる。そこで、点全体の重心が原点となるように各 \(\left|\vec{V}_{i}\right|\) を \(\left|\vec{V}_{0i}\right|\) へ平行移動する。（この時、\(D_{ij}\) は不変）

すると、上記の式よりグラム行列が \(D\) のみから計算することができる。

ここで、グラム行列 \(G\) は実対称行列なので、ある直交行列 \(P\)（各列が固有ベクトル）を用いて対角化できる。

定義より、\(G_{i j} =\vec{V}_{i} \cdot \vec{V}_{j}\) という関係があるので、

より、\(V_{i,j} = \lambda_j^{1/2}E_{j,i}\) と距離 \(D_{ij}\) からグラム行列 \(G\) を使って位置 \(\vec{V}_{i}\) を推定できることがわかる。

問題4 ヒトゲノムの配列構造

※この問題の解答は、Naoto-Yamaguchiに作成してもらいました。

1

• 物理的距離：遺伝子間の距離を塩基数で示したもの。ヒトゲノム解析で可能となった挿入・削除等により、個人差がある。
• 遺伝的距離：組み換え頻度から推定した遺伝子座の距離。１回の減数分裂で交叉が起きる回数の期待値が \(k\) となる２つの遺伝子座間の距離を \(k\ \mathrm{Morgan}\) と呼び、\(1\ \mathrm{cm} = 0.01\ \mathrm{Morgan}\) を単位として用いる。

両者の関係

両者の距離は、多型マーカーの代表例であるSNP間の距離を測るのに用い、ゲノム解読以前は、遺伝的距離、以降は物理的距離が主に用いられた。 しかし、交叉のHotspot（相同組換えの頻発部位）が存在することで、交叉の回数の期待値が場所に依存するため、両者は必ずしも比例しない

2

• アレル：相同の遺伝子座にあって、異なる遺伝情報を有する遺伝子のこと。対立遺伝子
• 遺伝子型：２つの対立遺伝子の組であり、同じアレルの組の時ホモ接合、異なるアレルの組の時ヘテロ接合と呼ぶ。

両者の関係

• アレルが優性遺伝形式に従う場合、そのアレルのホモ・ヘテロ接合で表現型として現れる。
• アレルが劣性遺伝形式に従う場合、そのアレルのホモ接合のみで表現型として現れる。
• 集団がハーディ・ワインバーグ平衡に従う場合、アレル頻度が \(A:a=p:q\) だとすると、\(AA:Aa:aa = p^2:2pq:q^2\) となる。

3

多型マーカーがHWEに従うことが重要なのは、以下のような平衡状態にない遺伝子型を排除することで、マーカーとして正しく機能することを保証できるからである。

• 近親婚などホモ接合が増える場合
• 遺伝子型の観測方法が不正確な場合

検定方法

ここで、遺伝子型ABを \(X\) 人のゲノムで観測した時、\(X\) 個以下となる確率の総和 \(\sum\left\{p(n_{AB})|n_{AB}\leq X\right\}\) が有意水準以下ならば、帰無仮説を棄却し、HWEにないとする。

4

ハプロタイプ

２倍体において片側の染色体上の遺伝子（or SNP）の列のこと。

連鎖不平衡

組換えがランダムな位置で多数起こった結果、ハプロタイプ上のアレルが、頻度率に従って並んだ状態のこと。

しかし、実際には近い遺伝子間など内部で組換えが起こりにくいハプロタイプである「連鎖不平衡ブロック」や交叉の起こりやすい「交叉のhot spot」が存在し、連鎖平衡でない連鎖不平衡の状態にあることが多い。

この不平衡状態を定量的に表す指標として、高頻度のアレルの連鎖不平衡係数 \(D_{AB}\) と、それを \(-1\sim+1\) に正規化した \(D_{AB}'\) を用いる。 この \(D_{AB}'=1\) とは、組換えを起こしておらず、完全に連鎖不平衡であることを示す。

また、各遺伝子間の \(D'\) を解析することで、連鎖不平衡ブロック、組換えのhot spotがわかった他、「組換え位置の特徴」として、組換えのhot spotのモチーフ配列や「結合タンパク質のPRDH9」が発見された。

連鎖不平衡ブロックの代表点を安価に調査し、多様な疾患の罹患リスク予測やハプロタイプ別祖先推定を行うサービスが23andMeから発表され話題を呼んだ。

問題5 Genome Wide Association Study(GWAS)

多重性は、「数撃てば当た」ってしまうことが問題である。例えば、サイコロを一回投げた場合、\(1\) の目が出る確率は \(1/6\) である。しかし、サイコロを \(n\) 回数投げた場合、少なくとも \(1\) 回 \(1\) の目が出る確率は \(1-\left(1-\left(\frac{1}{6}\right)\right)^n\) である。これは、次第に大きくなってしまう。

1

例えば \(n=5\times10^5\) 個の各マーカーを有意水準 \(\alpha=0.05\) で検定すると、有意になるマーカー数の期待値は \(2.5\times10^4\) 個となり、非常に多くなってしまう。

2

Bonferroniの補正は、「\(n\) 個のすべての位置の遺伝子型が対象疾患に相関しない」を有意水準 \(\alpha\) で検定する作業を「各位置において遺伝子型が対象疾患に関連しない」を有意水準 \(\alpha/n\) で検定する操作を繰り返すことで近似する

この時、

• 「全マーカーで帰無仮説『各位置において遺伝子型が対象疾患に関連しない』が棄却されない確率 \((1-\alpha/n)^n\) 」
• \(1-\)「帰無仮説『\(n\) 個のすべての位置の遺伝子型が対象疾患に相関しない』が棄却される確率」\(=\alpha\)

の２つが等しくなれば上手く近似されたことになるが、\(\alpha/n\approx0\) を仮定してテイラー展開を用いることで以下が成り立つので、ある程度の正当性は見込まれる。

ただし、(1)の例であれば有意水準が \(\alpha/n = 1\times10^{-7}\) となり、少し厳しすぎてしまう。

3

FDR(False discovery rate)は、真の帰無仮説を誤って棄却した割合の期待値を表している。

ここで、有意水準 \(\alpha\) で棄却されるマーカーの中で、疾患に関連しないが誤って棄却されるマーカーの数の期待値は \(n\alpha\) で表される。

よって、以下の手順で設定したFDRを満たす有意水準を見つけることができる。

1. \(\mathrm{FDR}\) を決める。
2. \(n\) 個のマーカーそれぞれに対して検定を行い、p値を計算する。
3. \(n\) 個の \(p\) 値を昇順に並べる。（\(p_1\leq p_2\leq\ldots\leq p_n\)）
4. \(i=m\) とおく。
5. \(np_i\leq i\mathrm{FDR}\) を満たす場合、\(k=i\) として6に進む。そうでなければ \(i\leftarrow i-1\) として再びこの操作を繰り返す。なお、\(i=1\) に達した場合、どの仮説も棄却することなく終了する。
6. 仮説 \(i=1,2,\ldots,k\) に対応するマーカーを全て棄却する。

※ ただし、一般に \(i \rightarrow i-1\) とした際に \(np_i/i\) が小さくなるとは限らない（例：\(p_i=p_{i-1}\)）ため、この方法にも限界がある。

補足

以下の式の意味を考える。

• \(np_i\) は、有意水準 \(p_i\) で棄却したマーカーの中で、疾患に関連しないが誤って棄却されるマーカーの数の期待値
• \(\mathrm{FDR}\) は、棄却されたマーカーの中で、疾患に関連しないが誤って棄却されるマーカーの数の期待値

したがって、上の不等式を満たす場合、関連しないマーカーを棄却してしまう確率を予め設定したFDRによって上から抑えることが出来る。

問題6 連鎖解析

※この問題の解答は、Naoto-Yamaguchiに作成してもらいました。

1

IBD vectorとは、ある家系の「３世代」の継承関係を元に孫の順序付遺伝子型を数値に変換してベクトル表現したもの。

子の各「アレル」が祖父由来なら \(0\)、祖母由来なら \(1\) と定義することで、4次元ベクトルで \(2^4=16\) 通りの継承関係を表現できる。

ここで、「罹患同胞対(りかんどうほうつい, affected sib-pair, 同じ疾患に罹患した兄弟で観察された共有する同祖遺伝子,IBD)」を分析することで、「劣性遺伝形式」ならIBD=2,「優性遺伝形式」ならIBD=1,2を説明できる。

2

\(X_i = \mathrm{AABC},\ \text{IBD vector} = 1100\)

• 子の順序付き遺伝子型は \(\mathrm{AABC},\mathrm{AACB}\)
• 親の順序付き遺伝子型は \(\mathrm{BACA},\mathrm{CABA}\)
• 求める確率 \(P(X_i,I_i) = 1/16 p_{\mathrm{A}}^2p_Bp_C\times2\)

\(X_i = \mathrm{ABAC},\ \text{IBD vector} = 0010\)

• 子の順序付き遺伝子型は \(\mathrm{ABAC},\mathrm{ABCA},\mathrm{BAAC},\mathrm{BACA}\)
• 親の順序付き遺伝子型は、子の順序付き遺伝子が \(\mathrm{BACA}\) の時のみ矛盾なく決めることができ、\(\mathrm{BCA?}\)
• 求める確率 \(P(X_i,I_i) = 1/16 p_{\mathrm{A}}p_Bp_C\)

（以下省略）

3

Hidden Markov Model を用いて、最もらしい（\(P(X_1\cdots X_M|I_i)\) を最大化する）IBD vector \(I_i\) の列を推定する。 →推定したIBD vectorから組み換え位置を予測できる!!

ここで、\(P(X_1|I_1) = 1,P(X_M|I_M) = 1\) と初期化して、以下の漸化式を回すことで、以下を最大化するIBD vector \(I_i\) の列を推定できる。

漸化式

※ 2つのIBD vector \(I_j\) と \(I_{j-1}\) の間で一致しない要素の数を \(k\) と定義している。

4

常染色体劣性形式の場合、IBD=2の場合でしか劣性遺伝（疾患の原因変異が劣性）を説明できないが、これは全16パターンあるIBD vectorのうち4パターンのみである。したがって、\(4/16=1/4\) の領域まで責任遺伝子領域を絞ることが出来るから。

問題7 AdaBoost

\(\sum_{i=1}^{N} w_{i}^{t+1} \leq\left(\sum_{i=1}^{N} w_{i}^{t}\right) \times 2 \varepsilon_{t}\) の証明

• \(\ast1\) \(\alpha>0\) で \(\gamma=0,1\) の時、以下が成り立つことは自明
  $$\alpha^{\gamma} = 1-(1-\alpha)\gamma$$
• \(\ast2\) 定義より、
  $$ \begin{aligned} 1-\varepsilon_{t} &=\sum_{i=1}^{N} p_{i}^{t}-\sum_{i=1}^{N} p_{i}^{t}\left|h_{t}\left(\vec{x}_{i}\right)-y_{i}\right| \\ &=\sum_{i=1}^{N} p_{i}^{t}\left(1-\left|h_{t}\left(\vec{x}_{i}\right)-y_{i}\right|\right) \\ &=\sum_{i=1}^{N} w_{i}^{t}\left(1-\left|h_{t}\left(\vec{x}_{i}\right)-y_{i}\right|\right) /\left(\sum_{i=1}^{N} w_{i}^{t}\right) \end{aligned} $$
• \(\ast3\) 定義より、\(\beta_t = \varepsilon_t/(1-\varepsilon_t)\) なので、
  $$1-\beta_t = \frac{1-2\varepsilon_t}{1-\varepsilon_t}$$

\(h_{f}\left(\vec{x}_{i}\right) \neq y_{i}\) ならば \(\prod_{t=1}^{T} \beta_{t}^{1-\left|h_{t}\left(\vec{x}_{i}\right)-y_{i}\right|} \geq\left(\prod_{t=1}^{T} \beta_{t}\right)^{1 / 2}\) の証明

※ 最終クラス分類器 \(h_f\) の出力は以下である。

• \(h_{f}\left(\vec{x}_{i}\right)=1, y_{i}=0\) の時 \(\sum_{t=1}^{T}\left(\log \beta_{t}\right)\) を \(\left(\ast\right)\) の両辺に加算すると、
  $$ \sum_{t=1}^{T}\left(\log \beta_{t}\right)\left(1-h_{t}\left(\vec{x}_{i}\right)\right) \geq \sum_{t=1}^{T}\left(\log \beta_{t}\right) \frac{1}{2} $$
  ここで、\(1-h_{t}\left(\vec{x}_{i}\right)=1-\left|h_{t}\left(\vec{x}_{i}\right)-y_{i}\right|\) が成り立つから、題意を満たす。
• \(h_{f}\left(\vec{x}_{i}\right)=0, y_{i}=1\) の時 \(\left(\ast\right)\) の \(\text { otherwise }\) であることから、
  $$\sum_{t=1}^{T}\left(-\log \beta_{t}\right) h_{t}\left(\vec{x}_{i}\right)<\sum_{t=1}^{T}\left(-\log \beta_{t}\right) \frac{1}{2}$$
  ここで、\(h_{t}\left(\vec{x}_{i}\right)=1-\left|h_{t}\left(\vec{x}_{i}\right)-y_{i}\right|\) が成り立つから、題意を満たす。

\(\sum_{i=1}^{N} w_{i}^{T+1} \geq \varepsilon\left(\prod_{t=1}^{T} \beta_{t}\right)^{1 / 2}\) の証明

※ここで、\(\varepsilon=\sum_{\left\{i | h_{f}\left(\vec{x}_{i}\right) \neq y_{i}\right\}} \frac{1}{N}\)

ここまでに証明した各種不等式を用いてBoosting Propertyを示す!!

\(\varepsilon \leq \prod_{t=1}^{T} 2 \sqrt{\varepsilon_{t}\left(1-\varepsilon_{t}\right)}=\prod_{t=1}^{T} \sqrt{1-\left(1-2 \varepsilon_{t}\right)^{2}}\) の照明

定義より \(\beta_{t}=\varepsilon_{t} /\left(1-\varepsilon_{t}\right)\) だから、