2014-01-21
次世代シークエンサーを使ったデータ解析〜オミックス統計学入門2014

ひとまず、次世代データ、Deep sequencing データの資料はこんなもの
「こうすれば、よい」という段階ではないので、コンセプト説明が重くなり、また、いざというときの汎用性重視になった…結局、「入門」としては役に立てにくいのだが…
まだ、1度も講義使用していないのでベータ版扱いです
90分で1-2回分相当
Rmdファイルです。html化、epub化できます(やり方はこちら)
html化、epub化が面倒くさければ、kindleで1米ドルでも(こちら)
次世代シークエンサーによるDeep Sequencing?オミックス統計学2014_3? 統計遺伝学
作者: ryamada
発売日: 2014/01/20
メディア: Kindle版
この商品を含むブログ (1件) を見る
# 次世代シークエンサーによるDeep sequencing

# はじめに
次世代シークエンサーによるDeep sequencingの概要に触れた後、
そのデータ産出特性に合わせて行われる確率的データ解析とはどんなことかを簡単にまとめ、
それに引き続いて、ゲノム・レファレンス配列を用いてバリアント検出をすることを対象として、データ処理・データ解析をするプロセスを確認する。
そのプロセスを確認しながら、そのデータ解釈やデータ処理のどういうところに確率や尤度やモデルが使われるのかに関して触れることで、確率的データ解析とその解釈の勘所を理解することを目標とする。
実務とは少しずれた話を中心にするが、それは、この分野のデータ解析はまだ実験系自体が発展途上であり、それに合わせたデータ解析手法も変化しており、その細部の知識が役に立つ期間は長くないと予想されるからである。逆に、どうして確率的にデータを扱うのか、その背景の考え方は何かを押さえておくことは、Deep sequencingにも有用でありそれ以外のデータ解析にも重要だからである。

# 次世代シークエンサーによるDeep Sequencingの基礎事項

## ゲノム・エピゲノム・トランスクリプトーム
次世代シークエンサーは核酸配列を読む。
従って、ゲノム(DNA)、トランスクリプトーム(RNA)の両方が対象となり、DNAの修飾状態によって対象を絞る技術と組み合わせることでエピゲノムの解析もできる(図)。

<img src="DeepSequencing.jpg" />

## 短いリードが大量に
次世代シークエンサーで読まれるリードは長くする方向での進展もあるが、現時点では短くてもよいから大量に高速に読むことでアドバンテージを稼ぐ技術となっている。
大量に読むと、ある特定の塩基箇所・範囲について複数回読むことになる。その複数の読み重ねを「深さ Depth」と言うのでDeep sequencingと言う。

そのため、図に示すように、データを活用するにあたって小規模実験データと異なる手順が必要となる
### 生データ処理

バルクで出てくるデータなので、それをひとまとめで扱うための処理がある。

### データのクオリティ・チェックとリードごとの処理

> データは全体としてどのようなクオリティなのかを、クオリティ指標の代表値と分布に照らして評価する。

> 同様に１本１本のリードごとに、そのリードのクオリティの評価をする。

> さらに個々のリードの個々の塩基のコールについても評価する。

### アセンブルかマッピングか

ヒトの全ゲノム配列のような、ある特定の配列(レファレンス配列)を地図として、その上のどこかしらにリード配列を割り当てることをマッピングと呼び、リード配列同士をつなぎ合わせて長い配列(contig)にすることをアセンブルと言う。

レファレンス配列がないときには図の左端のようにリード配列をアセンブルし、レファレンス配列があるときには、それ以外の場合のようにリード配列をマッピングする。

ジグソーパズルで言えば、完成図がない状態でピースをつなぐのがアセンブル、完成図に照らして絵を作るのがマッピングである。

## Deep sequensingの苦労
### 短いことの苦労
リードが短いと配列の特異性が低いので、区別のつかないリードがたくさんあり過ぎて、アセンブルもマッピングも難しい。
ジグソーパズルで言えば似過ぎたピースが多いのと、出来上がる絵のあちこちに似たような部分が多いのとで、完成が難しいのと同じようなものである。

### リードの精度による苦労
たくさんのリードが得られるがそれらの1つ1つが実験結果であり、きれいに読めているものとそうでないものとがある。高速に大量に読むので、個々のリードの精度には必然的にばらつきがある。不確か情報・偽情報もあるということである。
それらを用いることはジグソーパズルで言えばピースの絵柄に間違いがあったり汚れていて見えないところがあったりすること・ピースの輪郭の作り粗雑でピース同士(リード同市)をつなぐかつながないかの判断が難しかったりすること、に相当する。
これはアセンブル・マッピングの段階のことであるが、なんとか作り上げたアセンブリ・マッピング後情報を読み取る際にも苦労がある。同じくパズルに照らせば、２つのサンプルから作り上げた地図の違いを見つけるのは、間違いさがしであり、出来たパズルの絵に不確かさや汚れがあればその作業は難しい。
また、Deep sequencing では、リードというピースを同じ場所に積み上げ、その高さを場所ごとに比較したり(RNA-seqやChIP-seqの定量評価)、積み上がったところのアレルの比率を検討したり(バリアント検出、頻度推定)する作業も、同じく難しくなる。

### 場所ごとの読まれる回数のばらつきによる苦労
サンプルのどこがどれくらい読まれるかは、読まれる配列がどれくらいあるか、その配列がどれくらい読みやすいか、という要素に左右され、ばらつく。また、実験過程で配列を指数関数的に増幅する過程(PCR)を使う場合には、量のばらつきが増幅される。

元のサンプルがゲノム配列であるときは、2コピーずつあることが原則である。
この場合には、2コピーが均等に読まれていないことを考慮してバリアントの有無・ジェノタイプの推定をする必要が出る。また、部位ごとに得られる情報の量に違いが出るため、部位ごとにバリアント検出・アレル頻度推定の感度や特異度がばらつくという問題が生じる。
一方、読んだリードのデプスによって元のサンプルの量を定量的に評価したい場合(RNA-seq,ChIP-seq,プールしたゲノムを読んでアレル頻度定量をする)には、読まれる回数にばらつきがあることは、定量結果に直結する問題である。

# 確率的データ解析とは
## Deep sequencing における確率的データ解析の枠組み
以上のように次世代シークエンサーによるDeep sequencingでは複数の苦労がある。
これらはいずれも

> 不確かさのある情報が得られる

> 大量に得られるので分布が取れる

> 不確かさを分布によって補う

> 補うにあたって確率的判断をする

という枠組みで対処しようとしている。

## 確率・尤度、事前確率・事後確率
この確率的データ解析とは次のようにまとめられる。

> 知らないことがあって、それを知りたい(ある場所にバリアントがあるか、ある場所のアレル頻度はいくつか、ある座位のトランスクリプトはどれくらいあるか、など)

> 知らないことがもし、こうなっていたらと仮定する(バリアントがなかったら、バリアントがあったら、ある座位のトランスクリプトが別の座位に比べて●倍だったら、など)

> 仮定があったときに、ある結果がどれくらいの確率で生じるかを計算する

> 実際に得られた結果に照らして、どの仮定ではどれくらいの確率だったかを確認する。この確率を尤度と言う

> 仮定ごとに尤度の大小が出るので、それに応じて、どの仮定を信じるのがよいのかを判断する

これが、観察データに基づいて仮説を確率・尤度的に判断するプロセスである。

また、この確率・尤度の計算にあたって、色々と条件を想定する。
どのような配列は読まれやすいか(GC含量は読まれやすさに影響する要素)、読み間違いエラー率はどれくらいか・リードの個々の塩基はどれくらい信用できるか、などがそれにあたる。

このような条件を加味するときには、どのような値で取り込むかが問題になる。
ここで用いる条件の正確な値は知りえないからである。
条件に用いる値が変われば、それに基づいて計算する確率・尤度ももちろん変わり、最終結果ももちろん変わる。
従って、確率的データ解析の結果は、その真贋についても確率的に判断しなくてはならない。

また、大量のデータを取ってその分布を利用すると書いた。
この分布の利用にあたって、その分布がどのような形をしているのかにモデルとなる分布式を使うかもしれない。
モデルを入れると真実がそれとどれくらい合ってるのか、という点で、データ解析結果は影響を受ける。
モデルが複数考えられる場合には、どのモデルを使うかで結果が変わってくる。
従って、Deep sequencingデータの解析結果は、モデルによって結果が変わることも念頭に置いて解釈しなくてはならない。

# リシークエンシングによるバリアント検出

1人のDNAをレファレンス配列にマッピングし、バリアント箇所を検出する場合を考える。
検出結果の解釈に影響を与える点を中心にその過程をたどる。

## サンプル調整

読む個人を決めてDNAを採取・調整する。
全ゲノムにするか・エクソン部分だけにするか(エクソーム・シークエンシング)の選択肢があるが、

> エクソン部分だけにする場合には、その部分だけを取り出すというプロセスが入り、どのセグメントのどちらのアレル(母方・父方の２アレル)がどれだけsequencingに回るかに関してばらつきが入る。

> セグメントごとの取り出され方の差は、セグメントごとのDepthに影響が及ぶ。Depthが小さいセグメントは情報量が少なくなり、Depthが大きければ情報量が増えるので、セグメントごとにバリアントを見つけられるかどうかのパワーが異なり、バリアントであると判断した場合の信憑性が変わる

> バリアントは同一座位に１つのアレルがあるのか２つのアレルがあるのかを判断することで検出される。２アレルがサンプル調整の段階でどれだけ取り込まれるかにばらつきが出るので、そのばらつきがバリアントの検出感度に影響する。50:50で取り込まれた場合が最も感度が高くなる

> とりこみに差が出るのは、次のように分類できる

>> 片方がもう片方よりたまたま多く取り込まれることによって生じる

>> 片方が実験原理的に多く取り込まれる理由があることによって生じる

>>> バリアントのアレル自体が取り込まれに効率に影響を与える場合

>>> バリアントの周辺の配列差(近傍バリアント)が取り込みまれ効率に影響を与える場合


## マッピング前の評価

準備したサンプルをDeep sequencingする。

実験全体を評価する。

実験は生ものなので、その実験全体がうまく行っているのかどうかの評価が必要である。
その実験手法で「成功」とされる標準的な値、との範囲が量と質についてあるはずなので、それと比較する。

また、複数の実験をしているときには、実験単位で量的評価・質的評価をして、その異同を確認し、突出しているものについてはその原因を確認する。

### 量的評価

どれくらい読めたかの総量に関する評価が必要である。

> 全部で何リードか $n.read$

> それはリード長を加味して何塩基対分か $\sum_{i=1}^{n.read} L_i$ (ただし$L_i$は第i番リードの長さ)

> 読む対象の全長(ゲノム全長はいくつか、エクソームなら対象の全長はいくつか) $L_{total}$に対して、平均して何回ずつ読まれているか $\frac{\sum_{i=1}^{n.read} L_i}{L_{total}}$

最後の指標はある一定以上のDepthがあって初めてバリアントの検出が可能であることから、大まかな指標としてわかりやすい。
Depthとバリアントの検出とに関する確率的実習は、後述の『確率的データハンドリング：Depthとバリアントの検出』を参照。




### 質的評価

リード属性によって実験全体の質的評価をする

リードの属性としては、量的属性として長さ、質的属性として塩基単位のクオリティ値がある。
多数のリードについて、その量的・質的属性の評価をする。

また、リードの属性には、読まれ方・データのクオリティに影響する因子として、
リード内の相対的位置(読み始め・読み終わり・その中間)と、塩基の構成比(GC割合)などがあるので、
その条件での評価も有用である。

> リード長の評価をする

> リードの各塩基のクオリティ値の評価をする

> リードの相対的位置(読み始め・読み終わり・その中間)とクオリティ値には関係があるので、リードの相対的位置ごとのクオリティ値の評価をする

> リードのGC割合を条件として、リードのクオリティ代表値を評価する。

なお、属性の評価をするとは、その分布をとることと、代表値を確認することである。

その分布・代表値は、実験手法の標準的なそれと比較する。
また、自前の実験群で同様に比較し、突出があれば原因を確認する。

塩基別のクオリティは、実験上の塩基が実際の塩基と一致しているかどうかを表している。
つまり、リードの塩基は誤っているかもしれないわけである。
この誤った塩基を持ったリードがマッピングされれば、そこに本来はないはずのバリアントが存在することになる。
従って、クオリティ値が低めなリードのデータを採用すると、バリアントの検出にあたってFalse positiveを生じることになる。

この問題の解決法としては、間違いの可能性が非常に低い(全くない、ということはない、という立場に立つのが確率的データ解析である)リードのみ、そのような塩基のみを解析に用いるという方法がある。
しかしながら、これは、活用する総リード数・総塩基数を減らすことであり、Depthを小さくする。
Depthが小さくなればFalse negativeが増えるわけで、むやみに不採用のリード・塩基を増やすことは適当ではない。

コールのエラー率とFalse positiveの関係については、後述の『確率的データハンドリング：コールエラーによるFalse positive』を参照


## マッピングの評価

ここまでは、実験全体、リードの集まり全体に関する評価について述べてきた。
リードをレファレンス配列にマッピングするステップを考える。

読まれるべきリードがあまり読まれなければ、完璧にマッピングが成功してもDepthは読まれた総本数以上にはならないし、読まれたリードのクオリティが完璧でなければ、マッピング自体がうまく行っても、リードの配列自体に偽情報があるわけであるから、やはりマッピングはうまくいかない。後者の場合は、

> マッピングできる位置が見つからない

> 間違ったところにマッピングされる

> マッピングが複数の場所に対応づく

という場合が多くなる、という形で現れる。

この節では、それらについては、いったん横に置き、リード自体は完璧であるが、マッピングに支障が生じる場合について検討する。

本数とクオリティはそれに重畳する形で影響を与えると考えればよい。

### レファレンス配列との不一致

レファレンスゲノム配列は、ある代表的な配列である。
従って、ある特定の個人のゲノム配列が、それと一致しないことは珍しいことではない。
逆に、バリアント検出のためにDeep sequencingをしている場合には、一致しない箇所を探しているのであるから、一致しない〜マッピングに困難がある〜場合にこそ興味があるわけである。

問題は、レファレンス配列と一致しない座位を有するリードと一致するリードに比べてマッピングされないかもしれないことと、もし、そうであるとすると、それはどういう場合であって、どの程度、どういう形で影響を及ぼすかである。

レファレンス配列とサンプル配列との違い

> リード上にバリアントが一つあれば、ない場合に比べて、リードがレファレンス配列とマッチするスコアは下がる

> 近傍にバリアントがあれば、さらに下がり、場合によっては、レファレンス配列にマッチしたとのスコアはさらに下がる

> ゲノム配列には民族特有なバリアントがあるので、サンプルを提供した個人とレファレンス配列の提供者とに民族的な違いが大きければ、リード配列とレファレンス配列との不一致確率はさらに上昇する

> ゲノム上には似た配列がある。いわゆるリピート配列もそれに相当し、コピー数多型もこれに該当する。また、pseudogeneのような相同性の高い配列もある。その場合には、より複雑な状況が出現する

### バリアントを含むリードのマッピング

#### バリアントの密度、ハプロタイプ

ある部分が長さ$N$のリードとして読まれるときに、レファレンス配列と同じ配列を読めば、完璧なリードの場合には$n_A=N$塩基が一致する。
$k$個の１塩基バリアントが含まれれば、$n_B=N-k$塩基が一致する。

単純に、$t$箇所までの不一致は許容し、$t+1$箇所以上の不一致は許容しないとすれば、リードに含まれるほどの近傍にバリアントが$t-1$個までしかないバリアントは、100%マッピングされ、リードに含まれるほどの近傍に$t$個以上のバリアントがあって、自身を含めて$t+1$箇所以上が不一致となるバリアントは、まったくマッピングされない。

このような単純な割り切り方をする場合には、一定以上にバリアント密度が高い部位では、バリアントの検出感度が下がることがわかる。

また、この「近傍にバリアントがある」というときには、レファレンス塩基と異なる塩基が同一DNA分子上に並んで、複数のバリアントのハプロタイプを作っていることを指している。

従って、個人のゲノムにおいて、バリアントが同じように近傍に存在しているとしても、母方・父方で異なるDNA分子上にあるバリアントの場合より、バリアントが同一のDNA上にあってハプロタイプを形成している場合の方が、マッピングし損ねる率が上がる。

バリアントのゲノム上の密度と、ハプロタイプ構成とがマッピングに影響を与える。

バリアントの密度が高くなったり、民族性を反映して、そもそもレファレンス配列がマッピング先を含んでいない場合には、レファレンス配列を複数にしたり、レファレンス配列を三省せずにリード同士をつなぐ(アセンブル)の要素を取り込むなどを考慮する必要が生じる。

複数のレファレンス配列を用いる場合には、人工的に配列を重複させている側面もあるので、そのことに配慮したマッピング判定が必要となり、その判定は確率・尤度的な評価となる。
また、アセンブルを加味した場合には、それによって見出されるバリアントの感度と特異度は、アセンブルを加味しないそれとどのような関係にあるのかも、確率・尤度的に評価する必要がある。

#### リードのクオリティ

ここでは、簡単のために、塩基の一致数・不一致数による単純なマッピングを想定して述べたが、塩基コールのエラー率を加味してマッピングする場合にも、影響は及ぶ。

一致数・不一致数に関する閾値を設定する場合には、アレル別のリード数がバリアントで少なくなる傾向がある、という形で現れ、そのリード数は整数である。

他方、エラー率を加味すると、個々のリードがどこにどれくらいの確率でマッピングされるかの数字が細かい数字になる。
この確率は、エラー率を加味しても、バリアントを有する配列に由来するリードは、バリアントを有さない配列に由来するリードよりも低くなる。

マッピングにあたって、この確率を閾値に照らして、マッピングされるか否かを分けるのであれば、アレル別のリード数はやはり整数となり、バリアントを有する配列に由来するリード本数は少なくなる。
また、閾値を定めずに、マッピングされるかどうかの確率を用いて、重みづけ加算をするとしても、やはり、バリアントを有する配列に由来するリードの重みづけ加算値は小さくなる。

これは、バリアントの検出感度に影響する。
そして、その感度への影響は、個々の箇所についてはたいしたものではないが、多くの箇所をゲノムワイドに対象とする場合には、影響が無視しえない可能性が出てくる。

#### リードのクオリティとバリアントの存在

バリアント塩基を完璧なクオリティで読んだとき、そのリード配列とレファレンス配列との違いは、１塩基分になる。
同じ部分を読んだリードのバリアント塩基のクオリティが完璧でないときには、そのリードとレファレンス配列との違いは、1塩基分より小さいとみなせる。
単純に言うと、レファレンス配列との違いが小さい方が、マッピングの成功度は高いのであるから、クオリティが完璧でないリードの方が、マッピング成績が良くなる。

近傍にバリアントがあるときも同様である。
あるバリアントは完璧なクオリティで読まれているとする。
その近傍にあってハプロタイプを構成しているバリアントも完璧なクオリティで読まれれば、そのリードはかなりレファレンス配列から遠くなる。
一方、近傍バリアントの塩基クオリティが完璧ではないときには、かえってマッピング成績が上がるから、ここでも、クオリティが完璧でないリードの方が、マッピング成績が良くなる。

このように、「そもそも違っている」何かを検出するときに、「違いに曖昧さが入っている」方が、好成績でマッピングされたり、「曖昧さを持つ」リードを取り込んだ方が、検出されやすくなったりする可能性があることには注意が必要である。

クオリティ等の情報をすべて取り込んで、すべての可能性を確率的に扱えば、この問題のある部分は解消するかもしれないが、どこかしらのステップで閾値設定をすれば、その段差の影響が発生し、感度が下がったり上がったりすることになる。

また、すべての情報を取り込む、と簡単に書いたが、情報を取り込んだ上で、スコアリングするにあたっては、事前確率設定やモデル設定が影響を及ぼしてくるから、やはり、完璧ではない。

ここから言えることは、素朴な方法で臨んでも、工夫を凝らして情報を駆使しても、いずれにも何がしかの曖昧さや捉えきれない幅、false positive, false negativeが発生するということである。

その曖昧さの程度・大雑把さの程度についての感覚を得るには、実験系・処理フローとそのオプションが変わるごとに何がしかのシミュレーションをしてみることも役に立つ。


### リピートの存在、pseudogeneの存在

レファレンス配列には、多くの似通った配列が含まれる。
非コーディング遺伝子領域のリピート配列もそれに含まれるし、
遺伝子とそのpseudogeneもそうである。
また、複数のコーディング遺伝子が良く似た配列を共有している場合もある。
似通った配列が複数箇所にある場合の、マッピングの一般的な問題は、マッピング箇所が一意に決めにくいということである。

バリアントが存在する場合には、事情が込み入ってくる。

二箇所に似通った配列があるとする。

今、あるリードがこの二箇所のどちらか片方を読んだものだとすると、マッピングではどちらにもマッピングされうる。
二箇所の配列が全く同一であれば、優劣は付かない。
二箇所の配列が異なっていれば、優劣が付く。

今、二箇所のレファレンス配列は同一であるとする。
片方にだけバリアントがあるとすると、
バリアントのない配列を読んだリードは、二箇所に同確率でマッピングされる。
バリアントを含む配列を読んだリードの、二箇所へのマッピングスコアは、バリアントを含まない場合よりも低くなるが、二箇所で同一の値となる。

二箇所のレファレンス配列が１塩基だけ違っている場合には、ぞれぞれの座位から読まれたリードは、それぞれのレファレンス配列に最も高スコアとなるが、もう片方の座位にも、一塩基バリアントがある場合と同じ程度にはマッピングスコアが出る。

バリアント探索をしているときには、このようなマッピングの良さは捨てたくないレベルであるから、気がもめる。

わずかな差であっても一致の良い方に迷わずマッピングをするのであれば、紛れることはないが、リードのクオリティが下がると、二箇所へのマッピングの良さの差が縮まるから、問題はそれほど簡単ではない。

さらには、二箇所のレファレンス配列に一塩基の違いがあり、その塩基箇所について片方の座位にバリアントがあり、そのバリアントのアレルが、もう片方の座位のレファレンス配列と一致している場合には、バリアント配列からのリードは、本来の座位とは別の座位にマッピングされることになる。

また、コーディング遺伝子とそれに似たpsuedogeneとがあった場合で、両方が増幅されリードを生んでいるにも関わらず、マッピング先をコーディング遺伝子に限定したとすると、pseudogeneに由来するリードはマッピング先がないので、次善のマッピング先として、コーディング遺伝子にマッピングされることになるが、これによって見出されるリード配列とレファレンス配列の不一致はバリアントではない。

したがって、バリアント探索においては

> リード配列とレファレンス配列の異同だけではなく

> レファレンス配列上にありえる、複数の似た配列のセットと、リード配列との異同とを、併せて比較し、

> リードがバリアント配列を反映していたら、どのようなマッピングになるかを考慮する

必要がある。

## マッピング後の評価

### バリアント検出

バリアントの検出とは、サンプルのアレル別配列数が(レファレンスアレル):(バリアントアレル)の比として、$2:0$か$1:1$か$0:2$かのいずれであるかを識別することである。

以下の『確率的データハンドリング』に、諸条件がその検出にどのような影響を与えるかを述べた。

プールサンプルでのバリアント検出では、$M$人の場合にアレル別配列数が$2M:0,(2M-1):1,(2M-2):2,...,M:M,...,0:2M$を識別することになる。

### 定量比較

#### 同一サンプル内での、座位上の比率の定量

同一座位にマッピングされたリードの情報から、何かしらの定量をすることは、バリアントのジェノタイプ推定・プールサンプルのアレル頻度推定と基本的には同じである。
ジェノタイプのバリアントの方が整数比であることを前提に推定しているのに対して、こちらの定量は連続値であると考える点が異なる。

#### 同一サンプル内での座位間の定量比較

座位別の定量比較をする場合には、素朴にはマッピングされたリード数を比較する。

座位ごとに読まれやすさ$\beta$に違いあることを想定した上で、座位間の定量比較をするのであれば、$\beta$を用いたモデルに照らして、サンプル中のテンプレート量の推定をし、その量を比較することもできる。

また、２座位のテンプレート量の比を問題にするのであれば、同様に座位ごとの読まれやすさを加味した上で、比を推定することができる。

#### 異なるサンプル間での定量比較

絶対量の比較か、(標準配列に対する)相対量の比較かで推定するべき値は変わるが、基本的には、実験が異なるので、両者が同じ土俵に乗るように、標準化するか、総量確認をした上での比較となる。

# 確率的データハンドリング
## Depthとバリアントの検出~理想的な条件~

ここで、Depthがどのようにバリアント検出に影響を与えるかに関する確率的評価をしてみる。

今、サンプルに２つのアレルA,Bが$n_A,n_B$本ずつあるとする。
ゲノムの場合には理想的には$n_A=n_B$である。

この$n_A,n_B$本がそれぞれ、平均Depth $\beta_A,\beta_B$になるような読まれ方をしたとする。
そのときに、実験データとしてA,Bのアレルが$N_A,N_B$本のリードとして得られる確率について考える。

ある確率でぽつりぽつりと起きるような現象があって、ある時間内に平均$\beta$回、起きるとき、実際に$0,1,2,...$回起きる確率はポアソン分布に従うということが知られている。

実際Deep sequencingをして、座位ごとのDepthの分布をとると、ポアソン分布とみなしても悪くなさそうな分布が得られる。

```{r}
beta_A <- 3
N_A <- 0:30
x_A <- dpois(N_A,beta_A)
plot(N_A,x_A,type="h")
```
平均して$beta_A=3$回だとすると、1度も起きない確率もそれなりにあることも含めて、図のように示される。
今、ヘテロ接合体であるとして、$beta_A=beta_B$であると仮定する。
この仮定はどういうことかと言うと、個人のゲノムDNAの段階でA,Bの２つのアレルが等量ずつあり、そのDNA抽出・サンプル調整の過程で、２つのアレルがサンプルに等量ずつ入ることに成功し、さらに、そこから等しい効率でリードが読まれるという状況を設定している。

このとき、2アレル併せて$N=N_A+N_B$本のリードが得られたとする。
そのときの$(N_A,N_B)$の内訳は$(N_A,N_B)={(0,N),(1,N-1),...,(N,0)}$の$N+1$通りある。
このそれぞれがどのくらいの確率で起きるかを計算してみる。
そうすると、本当はヘテロ接合体であるのに、観察リードからはホモ接合体のように見える確率が計算できる。

```{r}
N <- 4
# アレルBはアレルAと同条件にする
beta_B <- beta_A
N_B <- N_A
x_B <- dpois(N_B,beta_B)
# A,Bが0:N,N:0で読まれる確率を積で計算
prob.N <- x_A[1:(N+1)]*x_B[(N+1):1]
# これが全体なので、和が1となるように補正
prob.N <- prob.N/sum(prob.N)
plot(0:N,prob.N,type="h",ylim=c(0,max(prob.N)))
# A,Bどちらかのリードしかない確率
prob.N[1]+prob.N[N+1]
```

では、この理想的な条件で総本数(Depth)が変わると誤ってホモ接合体とみなしてしまう確率が計算できる。

```{r}
Ns <- 1:50
prob.Ns <- rep(0,length(Ns))
N_A <- N_B <- 0:max(Ns)
x_A <- dpois(N_A,beta_A)
x_B <- dpois(N_B,beta_B)
for(i in 1:length(Ns)){
  N <- Ns[i]
  prob.N <- x_A[1:(N+1)]*x_B[(N+1):1]
  prob.N <- prob.N/sum(prob.N)
  prob.Ns[i] <- prob.N[1]+prob.N[N+1]
}
plot(Ns,prob.Ns,type="l")
# 確率を提示
round(prob.Ns,5)
```
$N=5$程度で、誤ってホモ接合体とみなす確率は0.05を切っていますから、Depth=5程度で信頼ができるとも考えられます。
しかしながら、ゲノム全体であれば30億塩基対($3\times 10^9$)あり、その100〜数100塩基に1箇所のバリアントがあるとすれば、$10^7$オーダーのバリアント箇所があり、それをエクソン部分(ゲノム全体の1%程度)に限定しても、まだ$10^5$箇所のバリアントが検出されるべきです。5%を見つけそこなうとすると、見つけそこなうバリアント数は$10^3$オーダーになります。
この見つけそこね数を3桁減らすには、「誤ってホモ接合体とみなす確率」を3桁小さくする必要が出ます。そうすると
```{r}
# 10^t で表示
round(log10(prob.Ns),2)
```
これで$-5$程度の値がほしいことになり、それは$N=19$程度であることがわかります。

実際、「このくらいのDepthがほしい」と思って実験計画をしても、実際にはDepthの小さい座位、大きい座位が生じる。
したがって、得られたDepthに応じて、
バリアントを見つけそこねている確率(False Negative)の率の概数を織り込みながら、データを解釈することが必要となる。

## Depthとバリアントの検出~理想的でない条件~

上のシミュレーションでは、ヘテロ接合体ならば、サンプル中の2アレルの割合が等しく、そのリード産生効率も等しいと仮定した。
また、実験エラーはないものとした。

サンプル中のアレル比、そのリード産生効率の不均等は、ヘテロ接合体であるのに、ホモ接合体であるように見えるデータを生む、と言う意味で、感度を下げる働きをする。

逆に、実験エラーがゼロではないことを想定すると(実際、ゼロではない)、バリアントはないのに、バリアントであるかのような実験結果となる(False positive)。

この２つのタイプの「非理想的条件」について、実習する。

なお、ここでは、細かいことよりも、条件によって結果の解釈がどのように変わるかを実感することを主目的とする。

### 不均等

サンプル調整段階でアレル別のDNA分子数の不均等とDeep sequencingにおけるリード産生の不均等に分けられる。

2種類のものから構成されている集まりを考える。
それがたくさんちょうど同数あるところとする。
そこから、適当に抜き出すとき、抜き出したものはおよそ同数であり、抜き出す数が多ければ多いほどその比率は半々に近くなる。
逆に言えば、少ないサンプルでは比率がばらつく。

大容量の実験系・十分な濃度の標本を用いた実験系・よく混ぜ合わせてある実験系では、偏りが生じにくく、比率が半々になりやすいが、その反対の条件では半々からずれる。

Deep sequeincingはさまざまなところでミクロ化しているので、偏りを生じる要素は少なくない。

さらに、それを増幅する過程があることで偏りは大きくなる。

ここでは、どのような過程がどの程度の偏りを生じるかについては触れず、偏りが生じた条件でDepthとFalse negative(ヘテロ接合体なのにホモ接合体のようなデータが得られること)とがどのような関係にあるかを確認する。

### 不均等状態でのDepthとFalse negativeの関係

不均等状態でA,Bのアレルから平均$\beta_A \ne \beta_B$本が観察されるような状況になったとする。
このとき、ヘテロ接合体であるのに観察リードのすべてが同一アレルになる確率(False Negative率)を計算してみる。
上の例で$\beta_A=\beta_B$の場合をやってあるから、$\beta_B = \beta_A \times (1+\alpha); \alpha =0,0.5,1,1.5,2,...,10$として計算する。

```{r}
Ns <- 1:50
beta_A <- 3
N_A <- N_B <- 0:max(Ns)
x_A <- dpois(N_A,beta_A)
alpha <- seq(from=0,to=10,by=0.5)
prob.Ns <- matrix(0,length(alpha),length(Ns))
for(j in 1:length(alpha)){
  beta_B <- beta_A*(1+alpha[j])
  x_B <- dpois(N_B,beta_B)
  for(i in 1:length(Ns)){
   N <- Ns[i]
   prob.N <- x_A[1:(N+1)]*x_B[(N+1):1]
   prob.N <- prob.N/sum(prob.N)
   prob.Ns[j,i] <- prob.N[1]+prob.N[N+1]
  }
}

matplot(Ns,t(prob.Ns),type="l",xlab="Depth")
```

両者の比が大きくなるほど、Depthが大きくないとFalse Positiveになることがわかる。
ここで、両者の比を最大で10としたが、これは、指数関数的増加プロセスが何かの拍子で作用すれば、生じるのにそれほど難しい数字ではない。

また、プールサンプルにおけるアレル頻度推定をする場合には、A,Bアレルが10% vs. 90 % などは、見出したいと考えているだろう。このときのA,Bアレルの比は10である。
これは5人をプールして、10本の染色体のうち1本がバリアントである場合である。
このことに鑑み、より大きなA,Bアレル比、$\beta_B = \beta_A*2^{\alpha};\alpha = 0,1,2,...,6$も実施してみる。これは
```{r}
print(paste("アレル頻度比",2^(0:6),"倍"))
      
print(paste("アレル頻度",1/(1+2^(0:6))))
```
のことである。
```{r}
Ns <- 2^(0:10)
beta_A <- 3
N_A <- N_B <- 0:max(Ns)
x_A <- dpois(N_A,beta_A)
alpha <- 0:6
prob.Ns <- matrix(0,length(alpha),length(Ns))
for(j in 1:length(alpha)){
  beta_B <- beta_A*2^alpha[j]
  x_B <- dpois(N_B,beta_B)
  for(i in 1:length(Ns)){
   N <- Ns[i]
   prob.N <- x_A[1:(N+1)]*x_B[(N+1):1]
   prob.N <- prob.N/sum(prob.N)
   prob.Ns[j,i] <- prob.N[1]+prob.N[N+1]
  }
}

matplot(Ns,t(prob.Ns),type="l",xlab="Depth")
```

### ポアソン分布仮定よりも負の二項分布仮定
以下は、番外編。
部位によらず一定のパラメタで表されるポアソン分布に従うと考えてきたが、部位によって$\beta$の値がばらつくのが実情である。

場所によって$\beta$がばらつくと、どこでも$\beta$が一定であると仮定したときと比べて、Depthの分布のばらつきが大きくなる。

実際、部位ごとに$\beta$の値をいくつにするのかについて、$\beta$値に影響を与える因子で補正してやると、
部位ごとに$\beta$の値のより真実に近いと思われる値が推定できることが確かめられている。

その$\beta$を基にすると、実際のDeep sequencingデータでの、部位別のDepthをポアソン分布とみなすモデルのあてはまりは、$\beta$値をゲノム全体で一様にした場合よりも、実データのあてはまりが良くなる。
しかしながら、まだあてはまりは十分ではない。それをさらに改善するために、負の二項分布を用いることもできる。

これはどういうことかというと、平均してDepthが$\beta$になるような座位で、実際に何本が読まれるかが、ポアソン分布よりもばらつくことを意味している。読まれる本数がばらつくということは、たまたまあるアレルだけ読まれない、という確率が高くなる、ということである。

従って、現実には、ポアソン分布を仮定するよりも、False negativeの起きる確率は高く、その高いFalse negative率を推定するには、負の二項分布モデルを使った計算がより正確だ、ということである。

以下は、負の二項分布をモデルとした際の計算という細部である。

その場合には、さらに、False neagtive率が上がることをシミュレーションしたのが下記である。

ここでのポイントは、False negativeをもたらしている要因を取り込み切らなければ、汲み取っていないFalose positive率があるから、宝探しをしているときには、その見落としの可能性に注意を払い、どうしても確実に見つけたい部分があるならば、そこについては、単純なモデルのFalse negative率で安心せずにDepthを上げる必要がある、ということである。
```{r}
my.dnegbinom <- function(x,k,r){
  choose(k+r-1,k)*(1-x)^r*x^k
}
r <- 20
p <- r*beta_A/(1+r*beta_A)
x_A.NB <- my.dnegbinom(p,N_A,r)
x_B.NB <- my.dnegbinom(p,N_B,r)
prob.N.NB <- x_A.NB[1:(N+1)]*x_B.NB[(N+1):1]
# これが全体なので、和が1となるように補正
prob.N.NB <- prob.N.NB/sum(prob.N.NB)
plot(0:N,prob.N.NB,type="l",ylim=c(0,max(prob.N.NB)))
# A,Bどちらかのリードしかない確率
prob.N.NB[1]+prob.N.NB[N+1] 
```

### コールエラーによるFalse positive

実験データであるから、不確かなデータ、誤った塩基コールは必ずある。
クオリティ値は、塩基ごとに、どのくらいのエラーの可能性があるかについて、数値として提供している。

今、ある座位についてホモ接合体のサンプルを読み$N$リードが得られたとすれば、コールエラー率が$e=0,10^{\alpha};\alpha=-5,-4,-3,-2,-1$のとき、$N$本のすべてにエラーが入らない確率は$1-(1-e)^N$である。


```{r}
# Depth
Ns <- 1:100
# エラー率
es <- c(0,10^((-5):(-1)))
X <- matrix(0,length(Ns),length(es))
for(i in 1:length(Ns)){
  X[i, ] <- 1-(1-es)^Ns[i]
}
matplot(X,type="l")
```

では、エラーコール率が$e=0.01$で一定のとき、読まれる総本数が$N=10$のとき、エラーリードの数が$Ne=0,1,2,...$の確率はいくつかを確認する。

```{r}
e <- 0.01
N　<- 10
X <- dbinom(0:N,N,e)
plot(0:N,X,type="h")
```

今、Ne=0の確率は
```{r}
X[1]
```
と、確実と言うには頼りないが、
今、Ne=0,1の場合の確率を合算すると
```{r}
sum(X[1:2])
```
とかなり高いので、$Ne=0,1$の場合ならば、ホモ接合とみなしてもよいと考えるかもしれない。

では$e$が増えるとどうなるだろうか。

```{r}
e <- 0.05
N　<- 10
X <- dbinom(0:N,N,e)
plot(0:N,X,type="h")
cumsum(X)
```

示した数字の4番目か5番目当たり($Ne=3,4$)までは、ホモ接合体でも起きうると考えないといけないかもしれない。

特に、ゲノム上の大多数の塩基はホモ接合なので、Deep sequencingで大量の塩基対($10^6$〜$10^9$)を読むきには、もっと厳しい値$0.99999999$くらいまではホモ接合とみなしておかないとFalse positiveの確率は低くても、False positiveの件数は馬鹿にならなくなる。

$e=0.01$に戻して、$N$を増やしてみよう

```{r}
e <- 0.01
N <- 50
X <- dbinom(0:N,N,e)
plot(0:N,X,type="h")
cumsum(X)
```

$N=50$リードのうち、5,6リードはエラーとして許容しないといけないことがわかる。

### False negative とFalse positiveとのバランスを取る

バリアント検出にあたっては、若干のエラーを許容し、それに伴って、ある割合のリード数がレファレンスと異なっていても、ホモ接合とみなさなくていけない、という要請と、ヘテロ接合であっても、リードの2アレル比は半々からずれるが、どれくらいまでの偏りはヘテロとみなすか、という要請とのせめぎ合いとなる。

その様子を確認する。

```{r}
e <- 0.01
N <- 10
X.homo <- dbinom(0:N,N,e)
beta_A <- beta_B <- N/2
N_A <- N_B <- 0:N
# ヘテロについては理想的なポアソン分布モデルを使用
x_A <- dpois(N_A,beta_A)
x_B <- dpois(N_B,beta_B)
X.hetero <- x_A*x_B[(N+1):1]
X.hetero <- X.hetero/sum(X.hetero)
matplot(0:N,cbind(X.homo,-X.hetero),type="h",main="e=0.01,N=10")
```

上向きの黒いバーはホモ接合のときにバリアントに見えるリード本数が0,1,2,...,N本、観察される確率であり、
下向きの赤い破線バーはヘテロ接合のときにあるアレルが観察される本数が0,1,2,...,N本である確率である。

$N=10$本のリードが観察されるときに、少ないアレルのリード本数が1本のとき、2本のときには、ホモ接合・ヘテロ接合のどちらにしても、目に見える程度のFalse positive, False negativeがあることがわかる。

Nの値が大きくなれば、2つの山の分離は良くなる。
その様子は以下の通り。

```{r}
e <- 0.01
N <- 30
X.homo <- dbinom(0:N,N,e)
beta_A <- beta_B <- N/2
N_A <- N_B <- 0:N
# ヘテロについては理想的なポアソン分布モデルを使用
x_A <- dpois(N_A,beta_A)
x_B <- dpois(N_B,beta_B)
X.hetero <- x_A*x_B[(N+1):1]
X.hetero <- X.hetero/sum(X.hetero)
matplot(0:N,cbind(X.homo,-X.hetero),type="h",main="e=0.01,N=30")
```

また、eの値を小さくすれば、黒の山はごく左端に片寄る。

しかしながら、Nの値が大きい場合に限定すれば、判断できない座位が増えるし、eの値を小さくすれば、活用できるリード数・塩基数が減るので、そのバランスを取ることとなる。
もしくは、バランスを取らずに、各座位のDepth・塩基のクオリティごとに、
False positive,False negativeの確率を算出しながら、情報をすべて活用する、ということになる。

なお、False positive, False negativeの値は、エラー率、ポアソン分布モデルに依存していることも、再度、指摘する。
その想定を変更すれば、おのずと、判断に用いるべき、黒・赤ピークの分離閾値や、False positive, False netative率も変化する。