定義：ゲノムアセンブリ
※ インデックスは0-originである。
| 記号・用語 | 定義 |
|---|---|
| $$\Sigma=\{A,C,G,T\}$$ | 塩基の集合 |
| $$\Sigma^{\ast}$$ | 全ての配列 |
| $$^-:\Sigma^\ast\rightarrow\Sigma^\ast$$ | 配列の相補塩基対（通常のワトソンクリックの相補対）を取るという演算 |
| $$\mathbf{G}$$ | ゲノム。塩基上の文字列 |
| $$G=|\mathbf{G}|$$ | ゲノムの文字列長 |
| $$r$$ | リード。ゲノム $\mathbf{G}$ からサンプリングされた、長さ $l \in \mathbb{N},l < G$ の配列。より正確には、$[0\cdots G-1]$ 上の一様分布に従う確率変数 $X$ の関数 $\mathbf{G}[X\cdots X+l)$ なお、スライスがインデックスの範囲を超えた場合は、環状ゲノムのことを示していると考える。$$\mathbf{G}[X\cdots X+l] = \mathbf{G}[X\cdots G] + \mathbf{G}[0\cdots X-G+l] $$ |
| $$R$$ | 独立にサンプリングした $N\in\mathbb{N}$ 本のリードの集合 $\{r_1,r_2,\cdots,r_N\}$ 各リードの長さは固定長 $l$ とする。 |
| $$c=lN/G$$ | 平均深度（データサイズ/ゲノムサイズ）。$\times10.4$ などと記されている。ここで、$\mathbf{G}$ がわかっていないにも関わらず $G$ がわかるのは、染色体の大きさから推定することができるからである。 |
| 深度 | 何回読まれたか |
| 被覆度 | どのくらいの割合が、最低1回読まれるか。 |
| 強連結成分 $V^{\prime}$ | 頂点の部分集合で（$V^{\prime}\subset V(G)$）、「$V^{\prime}$ の頂点の中では互いに行き来できる」が、「$V^{\prime}$ から外には行き来できない（行くだけ、来るだけは可能でも良い。）」もの。 |
Lemma 1 (Coverage)
Proof
Lemma 2
Proof
以上より、「十分多くのリードを読めば、全ての箇所が読まれる確率が1に近づく」ことがわかる。
またここで、ゲノム $\mathbf{G}$ がわかっていて、2つのリード $r_i,r_j$ が $l\ast\theta\ (0<\theta\leq1)$ だけ重なっていれば、その2つは「重なった部分から取ってこられた」と検知できるとする。
Lemma 3
Proof
なお、これはアセンブリが完全に実現できた場合の解析であり、ある種のアセンブリの上界を与えていると捉えることができる。
ところで、ここにおける確証のためのパラメータ $\theta$ はどのくらいだろうか？これはゲノムの構成に依存する。例えば、ゲノムの中に、長い $A$ の連続が2箇所以上あった場合、ポリAのリードに対しては、$\theta=1$ でも、私達が思っている意味での確証はできないように思われる。敷衍して言えば、このパラメータ $\theta$ は、ゲノム配列とゲノム上の位置に依存するパラメータであって、どちらもアセンブリの前には手に入らないものだ。
さて、では「これらのリード $R$ から、元のゲノム $\mathbf{G}$ を復元する方法」について考える。ここでは、以下の2つを挙げる。
OLC¶
最も愚直な手法は、以下のOLCである。
- (Overlap) 各リード $r_i$ と $r_j$ が、ゲノム上で連続した領域にあるかを判断する。
- (Layout) これらの情報を元に、リード $r_i$ を並べる。
- (Consensus) 並んだリードにたいして、各位置で、何らかの手法（例えば多数決）で塩基を確定させる。
OLCは、ゲノムアセンブリ初期に行われ、現在でも、ロングリードシークエンサーの発達とともに再び日の目を見ている。しかし、
- Overlapは計算が重い： 「全リードvs.全リード」の計算を行う必要があり、Suffix Array等で前計算を行い（$O(n\ast l)$）、SIMD(Single Instruction Multiple Data)等で並列にDPを行うことで $O(n\ast l)$ で"近似的に"実行することができるが、それでも重い。
- Overlapが真であるかを確かめることが難しい： 「$r_i,r_j$ がアラインメントをすると良いスコアでアラインできる」からといって「$r_i,r_j$ がゲノム上で重複する部分から読まれている」とは言えない。
- Layoutは、厳密解を求めるのが難しい： 「各リードをノードとして、overlapしている組をエッジで結んだグラフ」において、「全てのノードを一回だけ使って、全てのノードを回るパスを発見すること」であるが、この問題はNP完全に属するハミルトン路問題と同一である。
といった問題から、OLCのパラダイムに基づいて、ゲノムアセンブリを実行する際には、発見的方法、制約条件や作業仮説をつけることが必須になる。
de Bruijn Graph¶
定義：de Bruijn Graph
正の整数 $k$ に対して、$k$-de Bruijn Graphとは、次のような有向グラフ $G=(V,E)$ のことである。（文字に対しては半開区間でスライスを作ることに注意）
- $V\subset\Sigma^k$
- $(u,v)\in E\Longleftrightarrow u[1,\ldots,k) = v[0,\ldots,k-1)$
簡単に言えば、$k$-mer の集合を、『後ろ $k-1$ 文字と先頭 $k-1$ 文字が一致したら辺を引く』というルールに基づいて辺を引けるだけ引いたグラフである。
ここで、リードの集合 $R$ から作られるde Bruijn Graphを、$R$ 中に出現する $k$-merをノードとするde Bruijn Graphだとする。つまり、$k+1$ より長いリードに関しては、短い断片にバラバラにしてしまう。（この時 $R$ 中に何回その $k$-merが現れたかを数えてノードに計数としてのせておく。）
それでは、このde Bruijn Graphを用いてアセンブリを行う手法について考えるが、これは「ある $k$ を固定して、与えられたリードから、全ての $k$-merを抽出してde Bruijn Graphを作り、オイラー閉路（環状ゲノム）/オイラーパス（線状ゲノム）を構築する」ことに他ならない。また、このアルゴリズムは明らかに $O(|E|)$ で計算できる。
一般に、グラフ $G=(V,E)$ のオイラー閉路/パスを探すアルゴリズムを実装する際は、次のようにパス $P$、開始頂点 $v$ を持っておいて、深さ優先探索を一度だけ、グラフ $G$ を破壊しながら、再帰呼び出して行う擬似コードが用いられることが多い。
1. For each nodes in {u | (v,u) ∈ E}
2. (a) Remove (v,u) from E
(b) Call DFS(G,u)
3. Push v to Pこのプログラムは、スタック数の上限にさえ気をつければ、辺の数に関して線形時間で解ける。本来は難しいはずの問題が、辺の数に関して線形時間で解けるのは、どこかでおかしいことが起きているということだと思われる。
- 「ここまでの定式化のどこがおかしかったのだろうか？」
- 「現実世界で目にするインスタンスの中で、先程のオイラー閉路/オイラーパスとして解けないものには、どのようなものがあるだろうか？」
詳しく見ていく。
相補鎖¶
シークエンシングを行う際に、「相補鎖・鋳型鎖」のどちらかを選択することは困難であり、またあるリードの対に対して、それらが同じ鎖から出てきたのか、それとも相補鎖から出てきたのかを知ることもできない。このとき、
- 2つの鎖が共通の $k$-merを持っていない：2つの弱連結成分からなり、オイラー閉路/パスが2つ獲得できる。
- 2つの鎖が共通の $k$-merを持っている：グラフを巡回する間に別の鎖に乗り移ってしまうことがあり、これに対処するために以下の面倒な実装が必要となる。
- 各ノードに $k$-merとその相補鎖を対応させる。
- ノード間のエッジが必ず両方向に伸びる。
- 自分自身への多重ループができうる。
- エッジに「鋳型鎖-相補鎖」というようなラベル付けをする。
リピート（反復配列）¶
- 2つの異なったゲノム上の位置を占める配列が、全く同じ $k$-merを持ってしまう場合（$\ldots ARB\ldots CRD\ldots$）：
- ここから構築されるde Bruijn Graohから、グラフが「$\ldots ARB\ldots CRD\ldots$」から読んだリードから構成されたものなのか、「$\ldots ARD\ldots CRB\ldots$」から構成されたものなのかを判別することができない。
- このようなグラフから、オイラー閉路/オイラーパスを構築できるかも定かではない。（パスが合流するところ（$R[0]$）では「出次数<入次数」となるが、分流する場所ではその逆となる。）
- paired endというシークエンシング法を使えば、元のリードの情報から、オイラー閉路を選んだり・ノードを複製したり・エッジを複製したりしてアセンブリ結果を良くすることができる。
- タンデムリピートの場合（$\ldots RRRR\ldots RRRRR\ldots$）
- ループ構造になり、$R$ の繰り返し回数を推定することが極めて難しい。
※ これらの難しい問題に対処するために、現代的なアセンブラでは、リピートに関しては、解決できない部分は諦めて、確信が持てる– 必ず正しいような– 配列のみを出力するようになっている。
シークエンシングエラー¶
- de Bruijn Graphは、「全てのリードが完全に正しいこと」を要求しており、一塩基の誤りが、$k$ 個の誤った頂点と、$k − 1$ 個の誤った辺を誘導する。
- （辺の重みを持たない）de Bruijn Graph そのものには、このようなエラーを訂正する機能がないため、オイラー閉路/パスを辿る過程でも、エラーかどうか判定できない。
- 被覆度を上げるためには十分な量のリードを読む（スループットを上げる）必要があるが、これに伴って、エラーの総数が、リードの個数に線形に増えていってしまう。
- 対処法としては、以下が挙げられる。
- Fastqファイルの信頼度を利用し、
- 何らかの閾値を用いて塩基の精度が低いものを取り除く。
- ほとんど現れない $k$-merに関しては、単に無視する。
- de Bruijn Graphの実装において、
- ある頂点の近くだけで探索することによって、エラーを検知して、他の辺にマージする。
- 深さ優先探索をしていく中で、極端にすぐに行き止まる（本当は無い道を進んでしまった）場合は、直前の分岐までをなかったことにする。
- Fastqファイルの信頼度を利用し、
実装における工夫：BloomFilter¶
$k = 50$ 程度でde novo アセンブリを行うと、あり得る $k$-mer の数も多くなり、シークエンス深度によっては、数百ギガバイトものメモリを要求することもある。現在では、多数の実装が開発されており、その最適化の方法、実際の問題に対する性能もまちまちだが、今回はBloom Filter を用いた確率的 de Bruijn Graphについて記述する。
BloomFilterとは、1970 年にBloomらによって開発された、確率的データ構造である。このデータ構造は、ハッシュが計算できるような要素の集合に対する表現を与えるもので、空間効率が非常に高くなるほか、集合に対する所属、追加、そして集合同士の和を高速に計算できることが特徴である。
演算¶
BloomFilter $B,B^{\prime}$、要素 $e\in U$ に対して、次のように演算を定義する。
| # | 演算 | Return | 説明 |
|---|---|---|---|
| 1 | add:(BloomFilter:B, U:e) → BloomFilter |
$$B^{\prime} \text { such that } B^{\prime}[i]=B[i] \wedge h_{1}(e) \wedge h_{2}(e) \cdots \wedge h_{p}(e) \text { for all } i \in\{1, \cdots, M\}$$ | 各ハッシュ値を計算して、その番地にtrueを立てる。 |
| 2 | has:(BloomFilter:B, U:e) → {t,f} |
$$B\left[h_{1}(e)\right] \wedge B\left[h_{2}(e)\right] \cdots \wedge B\left[h_{p}(e)\right]$$ | 各ハッシュ値を計算して、それらの番地の場所に全てtrueが立っているか返す。 |
| 3 | union:(BloomFilter:B, BloomFilter:B′) → BloomFilter |
$$B\mid B^{\prime}$$ | 2つのBloomFilterのBitORを取って返す。 |
- 1つの要素 $e\in U$ に対して、$p$ 個のハッシュ値を計算し、それらが「全て
trueであれば含まれている」/「どれか一つでもfalseであれば含まれていない」とする。 - 一度BloomFilter $B$ に要素 $e$ を挿入(
B.add(e)すれば、以降、この要素が存在するかどうかのクエリ(B.has(e)は、必ず正しい答え(true) を返す。unionがBitORを取っていることから、この条件はunion後も満たされる。 - 挿入していないある要素 $e^{\prime}$ については、
B.has(e′)はほとんどの場合falseになるが、偶然、ハッシュ値 $h_1(e^{\prime})\cdots h_p(e^{\prime})$ がすでに全て埋められているときには、間違えてtrueを返してしまう。
※ つまり、BloomFilterには偽陰性はないが、偽陽性はある。直感的には、
- ハッシュ函数が少なすぎると、同じハッシュ函数値の集合を持った要素がたくさん出てきてしまう
- ハッシュ函数が多すぎると、すぐにBitVectorがすべて
trueで埋まってしまう。 - $N$ 個の要素を挿入した、$M$ 長のBitVectorと $p$ 個のハッシュ函数からなるBloomFilterの偽陽性は、簡単な計算から $\left(1-e^{-pN/M}\right)^p$ と近似できる。
- この近似から、最適なハッシュ函数の個数が $p=\frac{N}{M}\ln2$ の近くの整数になることがわかる。
このBloomFilterを用いた「ゆるい」de Bruijn Graohがprobabilistic de Bruijn Graphである。
- de Bruijn Graphにおいて、「2つのノードに対してエッジが引かれること」と「それらの $k-1$ 長の接頭辞と接尾辞が同じこと」は同値であるため、ノードを確定すれば、エッジは自動的に決まる。
- ゆえに、de Bruijn Graph を表現するためには、ある $k$-mer に対して、「それが含まれるか否かを答えるメソッド」が定義できればよいことになる。
- このとき問題となるのは、偽陽性であり、そのうちの「偽陽性であるとしてしまった $k$-merが、本当の $k$-mer間にパスを構成する時」である。（小さな引き込み路しか作らないなら、リードに混入するエラーを取り除くアルゴリズムを流用すれば対処可能。）
- うその $k$-merの集合を予め全て列挙（Black List）しておけば間違いは無くなるが、これだと情報論理の観点から、一切の圧縮ができなくなる。
- そこで、偽陽性のうちの重大なもの（`B.has(e)`とクエリされる $e$ のみ）を列挙しようと考える。
定義：Critical False Positive
| 記号・用語 | 定義 | definition |
|---|---|---|
| $$\mathcal{S}$$ | リード中の $k$-merの集合 | true positive nodes |
| $$\mathcal{E}$$ | $\mathcal{S}$ の要素と直接辺で結ばれている $k$-mer。言い換えれば、グラフを巡回する際に、クエリになる $k$-merの集合 | extensions of nodes from $\mathcal{S}$ |
| $$\mathcal{P}$$ | $e\in\mathcal{E}$ で BloomFilter $B$ が B.has(e) = true を返す $k$-merの集合 |
all elements of $\mathcal{E}$ for which the Bloom filter answers yes. |
ここで、我々が気をつけるべき偽陽性は、『$\mathcal{P}$ のうち、$\mathcal{S}$ には無いようなもの』である。このことをCritical False Positiveと呼ぶ。$\left(cFP = \mathcal{P}\setminus\mathcal{S}\right)$
課題1: Bloom Filter(One of two)
次の三題から二題を選んで解け。
- 集合 $U$ から $N$ 個の要素をランダムにサンプリングして、$M$ 長の BitVector を使った BloomFilter で表現する時、偽陽性が最も小さくなるようなハッシュ関数の個数を求めよ。ハッシュ関数は $U$ の要素を、$\{1,\cdots,M\}$ までの値に偏りなくマップすると考えてよい。
- $M$ を $3$ 以上の素数とする。$A,B$ を、$\{1,\cdots,M\}$ 上の一様分布に独立に従うものとする。確率変数 $Y_i=A\ast i+B, i=1,\cdots,t$ について、『$i\neq j\Rightarrow\text{$Y_i$ and $Y_j$ are independent.}$』であることを示せ。 ただし、素数 $M$ の剰余類で、連立方程式 $y_i=a\ast i+b, y_j=a\ast j+b$ は、$a,b$ について必ず一意に解けることを用いてもよい。ここから、BloomFilter において、$p$ 個のハッシュ関数をどのように構築すればよいか示せ。ハッシュ函数についてより知りたい場合は "Adam Kirsch and Michael Mitzenmacher. Less hashing, same performance: Building abetter Bloom filter.Random Structures and Algorithms, 33(2):187–218, sep 2008." などを参考にせよ。
- k-merの集合 $S$ に対して、CriticalFalseNegative を $O(\|S\|k)$ で構成するアルゴリズムを与えよ。ランダムアクセスメモリを仮定してよく、メモリ空間は必要な分だけ好きに確保できるとしてよい。
解答
1¶
Bit Vectorの $i$ 番目の位置に true を立てるハッシュ値の数を $B[i]$ で表す。ここで、$Np$ 個のハッシュ値が一様分布に従っていると仮定すると、
なる確率分布に従う。従って、偽陽性となる確率(False Positive Rate; FPR)は、以下で近似できる。
$$ \begin{aligned} \operatorname{FPR}(p) &= \left(1-\operatorname{Pr}\left\{B[i]=0\right\}\right)^p\\ &= \left(1-e^{-\lambda}\right)^p\\ &= \left(1 - e^{-pN/M}\right)^p \end{aligned} $$ここで、
$$ \begin{aligned} \frac{\partial\ln\left(\operatorname{FPR}(p)\right)}{\partial p} &= \frac{\partial}{\partial p}\left(p\ln\left(1 - e^{-pN/M}\right)\right)\\ &=\ln\left(1 - e^{-pN/M}\right) + \frac{pN}{M}\frac{e^{-pN/M}}{1 - e^{-pN/M}}\\ \therefore\left.\frac{\partial\ln\left(\operatorname{FPR}(p)\right)}{\partial p}\right|_{p = \frac{M}{N}\ln2} &=\ln\left(1 - \frac{1}{2}\right) + \ln2\frac{1/2}{1 - 1/2}\\ &= -\ln2 + \ln2 \\ &= 0 \end{aligned} $$となるため、最適なハッシュ函数の個数は整数値であるから、$p = \frac{M}{N}\ln2$ 付近の整数
2¶
素数 $M$ の剰余類で、連立方程式
$$\begin{cases} y_i&=a\ast i+b\\ y_j&=a\ast j+b \end{cases}$$は、$a,b$ について必ず一意に解けることを用いれば、$y_i$ が与えられた時、
$$y_j = a\ast \left(j-i\right) + y_i$$となり、これは $i\neq j$ の時、$y_i$ と独立になる。
従って、
$$h_i = A\ast i + B\quad\operatorname{mod} M$$というハッシュ関数を $p$ 個用意すれば良いことになる。
3¶
Critical False Positive は、「『リード中のk-merの集合 $\mathcal{S}$ の要素と直接辺で結ばれている $k$-merで、BloomFilter $B$ が B.has(e) = true を返すもの』 のうち、$\mathcal{S}$ に無いようなもの」を指す。
したがって、以下のアルゴリズムで、$O(\|S\|k)$ によってCritical False Positive の構築が可能である。
cFP = []
P = []
for kmer in S:
for e in kmer.extensions:
if B.has(e): # O(1)
# P.append(e)
if e not in S: # O(k): Using Trie Tree
cFP.append(e)
課題2: Information theory
- $n$ 通りのものから一つを選び、それを表現するときは、$\log n$ bit が必要になる（$n$ 個を『辞書式順序』で並べて、何番目かを伝えれば良い）。
- $n$ 通り値を取り得るものから、重複無しに $m$ 個選んだ集合を表す（i.e, $V\in U, |U|=n,|V|=m$ なる $V$ を表す）場合には、$\Omega\left(\log \left(\begin{array}{c}{n} \\ {m}\end{array}\right)\right)$ bit 必要になる。
ここで、$k$-de Bruijn Graph を表現するときは、頂点、すなわち $\Sigma^k$ の任意の部分集合を表現すれば良い。したがって、$k$-de Bruijn Graph を表現するには、$k$、ノードの個数 $n$ に続けて、どの $k$-mer を選んだかを伝えなければならず、最低でも以下のbitが必要になる（桁数を表すのに追加でbitが使われる）。
$$ \Omega\left(\log k+\log n+\log \left(\begin{array}{l}{\Sigma^{k}} \\ {n}\end{array}\right)\right) $$次の問いに答えよ。
- 一般に、$n$ 通りの値を取り得るものから、重複無しに $m$ 個選んだ集合を Bloom Filter で表現すると、上記で示した下界より小さくできる。これは下界を破っているように見えるが、そうではない理由を説明せよ（つまり、この BloomFilter は対象の集合を表せていないことを示せ）
- $k$-de Bruijn graph を、二つの BloomFilter を使って表現することを考える。このデータ構造は複数の BloomFilter からなり、さらに、de Bruijn graph をたどるようなクエリに対しては偽陽性はない。しかも、必要なbit数は上記の下界よりも小さい。一見、この状態は下界を破っているように見えるが、実はそうではない理由を説明せよ。
解答
raise KeyboardInterrupt(
"I didn't understand what the problem statement means."
)
課題3: de Bruijn Graph(Optional)
以下で用いるデータセットは、人為的に作成されており、次のような性質を満たす。
- カバレッジは十分量ある（被覆度は1になるように設計してある）
- 各リードは必ずゲノム中の位置からサンプリングされているし、エラーはない。
- どちらか一方の鎖からしかリードはサンプリングされていない。つまり、相補鎖について考える必要はない。
- あるkがあり、このkを用いて作った k-mer de Bruijn Graph は "ほとんど"リニアになっている。つまり、ほとんどのノードの入次数と出次数は1になる。
（確率的であれ、別のものであれ）de Bruijn Graph を実装し、.fasta fileと整数 $k$ を入力として、何らかの .fasta fileを出力するようなプログラムを書け（注：頑張りすぎないこと。ノードを選んで探索すればよい）。
余裕がある人は、2つのデータセット tiny と large に対して、何か $k$ を決めて、リードから $k$-deBruijn Graph を構築して、その上で何かしらのパスを求めてみよ。
参考のために、リードをサンプリングした配列もデータセットに含めておいた。
解答
from kerasy.utils import read_fastseq
true_seq = read_fastseq("dataset/tiny/assumed_genome_size1K_readlength_100.fasta")[0]
reads = read_fastseq("dataset/tiny/tiny.fa")
from kerasy.search.debruijn import kmer_deBruijnGraph
k = 89
model = kmer_deBruijnGraph(k=k)
model.build(reads)
print(f"num sequences : {model.num_reads}")
print(f"average length: {model.ave_read_length}")
model.export_graphviz(f"deBruijnGraph_{k}mer.png")
from kerasy.Bio.alignment import SmithWaterman
score_model = SmithWaterman()
score_model.load_params()
print(f"length of the true sequence: {len(true_seq)}")
print(true_seq)
for init in model.inits:
result_trail = model.assemble(init)
for result in result_trail:
score_model.align(true_seq, result, verbose=-1)
print()
print("*"*60)