解析の詳細
宮本 真理 Ph.D.
シニアフィールドバイオインフォマティクスサイエンティスト
㈱CLCバイオジャパン
はじめに
• 今日のセミナーでお話しすること
– データ解析の流れ
–
内部でのデータ処理の流れとその原理
• 今日のセミナーでお話ししないこと
– 詳細な使い方はお話ししませんが、デモにてどのように実行可能か
お話しします。パラメータについては、必要な個所はスライドに含めて
います。
• データフォーマットとクオリティ
• クオリティトリミング
• Duplicate Removal
• マッピング
• ローカルリアライメント
• 変異検出
• 変異検出解析例
• De Novoアセンブリ
解析の流れ
PCR Duplicate除去※ 変異検出 ローカルリアライメント※ マッピング アノテーション付け クオリティチェック インポートデータ準備
マッピング関連
変異検出と
アノテーション
データインポート
• 次世代シーケンサー解析には主に つの情報を使うためイン
ポートが必要です。
– リード配列
– 参照配列（ゲノムやコンティグ）
– アノテーションファイル
•
– 行が リード分
データフォーマット：
リード配列
@S1_89:1:1:672:654/1
AATTACATGACTCCAGTACCCAACCTCCACTTACCATTTATTGTCTCTGTCTACCCTGCCCAATGTCTCCCAGAAA
+
BBBBBCCBBCBBBBBB@CBBBCBBB<BBBBBBBABBABBBAAB@BABA>>AB@??=<>?<AA@@A>@@@@?@@6==
@S1_89:1:1:672:654/2
TAGACAGAGACAATAAATGGTAAGTGGAGGTTGGGTACTGGAGTCATGTAATTGAATCTGAAGGGAGAAGCTGCTA
+
BBBACBABBBBBBB=AABBB9BBB9>A@BA9@@A7+::7;:.71<5;<19?=913>==::6:7==6503>895263
@S1_89:1:1:657:649/1
ACAGAGACAATAAATGGTAAGGGGAGGTTGGGCACTGGCGTCATGTAATTGAATCTGAAGGGAGAAGCTGCTACAG
+
;;)95;5;2;;60;;;8+;02-;;3;;68;76&)614,)4'4)2411,0/4112$40/&/&4/0+/,&/1))++'-
•
– 行目：リードの
– 行目：配列
– 行目：＋
データフォーマット：
リード配列
@S1_89:1:1:672:654/1
AATTACATGACTCCAGTACCCAACCTCCACTTACCATTTATTGTCTCTGTCTACCCTGCCCAATGTCTCCCAGAAA
+
BBBBBCCBBCBBBBBB@CBBBCBBB<BBBBBBBABBABBBAAB@BABA>>AB@??=<>?<AA@@A>@@@@?@@6==
@S1_89:1:1:672:654/2
TAGACAGAGACAATAAATGGTAAGTGGAGGTTGGGTACTGGAGTCATGTAATTGAATCTGAAGGGAGAAGCTGCTA
+
BBBACBABBBBBBB=AABBB9BBB9>A@BA9@@A7+::7;:.71<5;<19?=913>==::6:7==6503>895263
@S1_89:1:1:657:649/1
ACAGAGACAATAAATGGTAAGGGGAGGTTGGGCACTGGCGTCATGTAATTGAATCTGAAGGGAGAAGCTGCTACAG
+
;;)95;5;2;;60;;;8+;02-;;3;;68;76&)614,)4'4)2411,0/4112$40/&/&4/0+/,&/1))++'-
データフォーマット：
リード配列
クオリティスコア
•
シーケンサーにて読まれたリードの塩基ごとのクオリティを示す値。
ファイル
では、記号にて表記されている。
•
へインポートされたデータはすべて
フォーマットのクオリ
ティスコアへ変換される。
•
リードの各塩基がエラーである確率を とすると、
のクオリティスコアは以下
𝑄
𝑆𝑎𝑛𝑔𝑒𝑟
= −10 log
10
𝑝
Phred Quality Score
エラー確率
Base Callの精度
10
1/10
90%
20
1/100
99%
30
1/1000
99.9%
40
1/10000
99.99%
データフォーマット：
リード配列
• その他、各社シーケンサーメーカーにより、提供するフォー
マットが異なる場合があります。
データフォーマット：
参照配列
•
•
データフォーマット：
アノテーション
•
遺伝子名や、その場所などを示している 列からなるファイル。
–
列目：染色体名や、コンティグ名
–
列目：ファイルを作成したプログラム名
–
列目：
や
などアノテーションのタイプ
–
列目：アノテーションの開始ポジション
–
列目：アノテーションの終了ポジション
–
列目：０から
までのスコア。ただし最後の列で
を指定している場合のみ。そ
れ以外は
ではブラウザでのグレーの濃さを表す。
–
列目：センス鎖かアンチセンス鎖か
–
列目：コーディングエクソン上にあるかどうか
–
列目：グループ情報
データフォーマット：
アノテーション
•
ターゲットリシーケンスなどを行った際のターゲット情報を保存したりするた
めに使われる
•
– 変異情報を保存するために使われる。
データフォーマット：その他
•
– マッピング後のデータについて、マッピングの状態を含めたファイル
–
は
のバイナリ形式のファイル
データフォーマット
• まとめ
– バイオインフォマティクス、
データ解析では様々なフォーマットが利
用されている。
– 拡張子で見たことのないものがあれば、とりあえず
で検索「
フォーマット」などでおおよそ検索可能
– フォーマットの変換が必要になることは比較的多く発生する。その時
も
で検索することで、様々な方法を調べることが出来るので、自
分が出来そうなものからトライしていく。
各社メーカーのクオリティとデータに見られる特性について
•
– 末端へ向かってみられるクオリティの低下
•
– ホモポリマー領域でのミスコール
–
エラーの可能性
–
リッチな領域でのバイアス
•
– ホモポリマー領域でのミスコール
目的別クオリティトリミングの考え方
•
アッセンブリ
– クオリティの低いリードは間違った
作成の原因に加え、グラフが
多く作成されるため、メモリの消費、処理時間増加につながる。
• 変異検出
– 明らかに低い
を取り除く。
– 変異検出では、変異の検出の際に、再度フィルターを個別にかけら
れるため、一部悪いデータが残っていてもあとでフィルタリングできる
ことを覚えておく。
–
も間違った変異検出の原因となるため、
の可能性があ
る場合は、取り除いておく。
クオリティチェック流れ
作成
– インポートしたリードのクオリティがどのぐらいか、その後のトリミング
や、
の状況などを確認するためにレポートを作成。
の除去 （オプション）
– フラグメント作成の過程で
が異常にかかってしまったものを補正。
トリミング （オプション）
– アダプターの除去、クオリティスコアによる除去、長さを指定した除去
などを選択・組み合わせてトリミング。
上記処理の後に再度
を作成すると処理前と処理後でのリー
ドのクオリティを比較でき、便利です。
レポート
•
• GC-content
長さの分布
リードごとに含まれるGC含有量を全
レポート
•
• Quality distribution
あいまいな塩基としてコールされたものがど
レポート
•
• Nucleotide contributions
塩基ごとのカバレッジ分布
レポート
•
• Ambiguous base-content
レポート
•
レポート
•
QCレポートは数字で見ることも可能。これに
より、具体的なリードの数や塩基数が把握
できる。これらのレポートはエクスポートする
対策
プラグイン
により２つのタイプの
除去が可能。
•
– リードファイルから
を除去。
– 参照配列がない場合にも利用可能。
– おもに
の前に利用される。
•
– マッピング後のファイルから
を
除去。
– 参照配列があることで、リードの向きを
考慮して
を取り除くことができる。
• PCRが異常にかかてってしまった影響を除去するためのプラグイン。
• リードの先頭20塩基、または20塩基が均一に4箇所一致する箇所があり、残りの
塩基が完全一致のアライメントスコアの80％以上あれば、PCR Duplicate と判断し、
そのリードを除去。
•
プラグインのインストールが必要です。
※プラグインのインストールは、管理者権限で行う必要があります。
•
リードの先頭 塩基以上、または 塩基以上が均等に 箇所において一致する箇所
がある。
•
かつ残りの塩基がコンセンサス配列のアライメントスコアの ％以上あれば、
と判断し、そのリードを除去。
注意点！
• リードのスタート地点や、アライメントスコアの高い一致をもっ
て、
の可能性としているため、トリミング前に実施して
トリミング
• あらかじめ登録されているアダプターの除去
• 新規で独自の配列を登録することも可能
アダプター除去
• Quality Score を使い、Quality の低い配列が連続
するようになる箇所からカット
• 正確に読めていない塩基をいくつ許容するか
クオリティトリミング
• 塩基数を指定して、5末端、3末端をカット
• Quality Scoreでカット後、短くなりすぎた配列を
カット
長さによる除去
クオリティトリミング原理
•
では
を使い、累積の
がある一定の値より大
きいものが続いた場合に、その箇所を取り除く、という処理を行います。
• 具体的には以下：
を 値へ変換
中に設定するパラメータ（
）と 値の差を計算
差の累積和を計算。このとき、 以下の値は とする
後のリード開始点は累積和がはじめて 以上になった点。
後のリード終
了点は累積和が最大の点
クオリティトリミング原理
0 20 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 リード配列 G C C C A T G T T C G A T G C Phred score 4 8 15 30 32 23 10 31 31 20 15 11 10 10 9 p値 0.40 0.16 0.03 0.00 0.00 0.01 0.10 0.00 0.00 0.01 0.03 0.08 0.10 0.10 0.13 Limit - p値 (D) -0.35 -0.11 0.02 0.05 0.05 0.04 -0.05 0.05 0.05 0.04 0.02 -0.03 -0.05 -0.05 -0.08 (D)の累積和 0.00 0.00 0.02 0.07 0.12 0.16 0.11 0.16 0.21 0.25 0.27 0.24 0.19 0.14 0.06 スタート点： 累積和が0より大 きくなった塩基 終了点： 累積和が最大を 示す塩基 Phred score の棒グラフ グラフより、ある程度クオリティが高くなった場所からリードを使い、クオリティが 連続して悪くなっている箇所からリードをトリムしていることがわかる。マッピング
• マッピングとは
– マッピングとはシーケンサーから算出されたリード配列を参照配列
（ゲノム）へマッピングする手法を指します。
Mapping
特徴
•
Suffix Array を使い参照配列をインデックス化し、高速なマッピングを可能
にしています。
•
ローカルアライメント・グローバルアライメントによるスコア計算が可能。
•
異なるシーケンステクノロジー、ペアエンド、シングルエンド、をあわせて
マッピング可能。
•
カラースペースによる配列のエラー補正も可能。
37マッピングの詳細
インデックスファイル作成
スコア計算
マッピングの詳細
• インデックスファイルの作成
– ？インデックス
Genome
Read
どこが似てるかな・・・？
Genomeの端から端まで順番に調べていては、膨大な時間がかかる
マッピングの詳細
•
にインデックスと言う辞書の索引のようなものを作成し、
検索効率を上げる。
Genome
Read
どこが似てるかな・・・？
Index
•
i
1
2
3
4
5
6
7
S[i] A
G
T
T
C
G
$
Suffix
i
AGTTCG$
1
GTTCG$
2
TTCG$
3
TCG$
4
CG$
5
G$
6
$
7
Suffix
i
$
7
AGTTCG$
1
CG$
5
G$
6
GTTCG$
2
TCG$
4
TTCG$
3
i
1
2
3
4
5
6
7
A[i] 7
1
5
6
2
4
3
Suffix Array, A
Sort
i
1
2
3
4
5
6
7
A[i]
7
1
5
6
2
4
3
1
$
A
C
G
G
T
T
2
G
G
$
T
C
T
3
T
$
T
G
C
4
T
C
$
G
5
C
G
$
6
G
$
7
$
i
1
2
3
4
5
6
7
S[i] A
G
T
T
C
G
$
マッピングの詳細
スコアリング
最適なマップ場所を
で探索
リード配列（
）が全て一致した場合
CGTATCAATCGATTACGCTATGAATG
||||||||||||||||||||
ATCAATCGATTACGCTATGA
20
マッピングの詳細
CGTATCAATCGATTACGCTATGAATG
|||||||||||||||||||
TTCAATCGATTACGCTATGA
19
CGTATCAATCGATTACGCTATGAATG
|||||| ||||||||||||
TTCAATCAATTACGCTATGA
16
CGTATCAATCGATTACGCTATGAATG
|||||| ||| |||||||
TTCAATCAATTGCGCTATGC
12
マッピングの詳細
•
マッピングのプロセスでは、各リードがもっとも高いアライメントスコア（参照配列との一致度
を示すスコア）を示す場所にマッピングをしています。しかしながら、時には近傍のリードの
マッピングの状況から、最も高いアライメントスコアではなくとも、もっともらしいマッピング結
果が考えられる場合があります。
•
たとえば上記例では、
は左横にずれることで、他のリードのマッピングとも一致し
もっともらしいマッピングになると考えられます。マッピングの段階では、各々のリードのアラ
イメントスコアのみを考えているため、このような状況が発生します。
•
では、このような状況を修正するため、マッピングを部分的にやり直します。
この際、通常のマッピングの段階とは異なり、他のリードのマッピング状況を考慮するため、
先ほどのマッピングは以下のように変化します。
•
先ほどのマッピングよりも、こちらの方がもっともらしい結果であることが直感的に分かりま
原理
原理
•
グラフにして書き直し、それぞれのパスを通るリードのカバレッジを記入すると以下のように
• Toolbox > NGS Core Tools >
Local Realignment
•
2種類のLocal Realignmentsがありま
す。さらにGuided にはNo forceと
Forceの2種類があります。
– Non guided
– Guided
• No force
• Force
• Guided Local Realignment
– ガイドとなるような変異（InsertionやDeletion）の情報をあらかじめ与えておく
ことで、その領域のInsertion、Deletionを考慮してリアライメントを行う。
– ガイドとなる変異情報がない場合、Local Realignment では、少なくとも1本の
リードがInsertionやDeletionを支持している必要がある。このような場合、ガイ
ドとなる変異情報を与えることで、InsertionやDeletion を効率的に検出できる
ようになる。
• マッピングされた後のリードから
を取り除きます。
注意点！
• リードのスタート地点を起点として考えているため、
を行い、リード
の 末端がトリムされると、正しく
が行われない可能性があります。
• クオリティの高いデータの場合、 末端がカットされることは非常に稀なた
め、トリミングによる影響が出ることが考えにくく、トリミング後のご利用に
大きな問題は考えられません（トリミングの設定のもよりますが）。
• どうしてもクオリティの低いデータで実施したい場合、マッピング後、マッ
プしたリードを抜出し、トリミングを行う方が安全です。
種類の検出方法
•
：
クオリティと、変異の見られる頻
度から変異のサイトを検出
以前の
。
•
：
確率モデルを使い、変異のサイト
を検出。
使い分け：
変異の見られる頻度が、その領域において ％以下のような場合は、
それよりも多い場合は、
をご利用ください。
Mapping後のデータに対し、を設定し、許容するミスマッ
チや、gap、またQuality ScoreによりSNP detectionに含
めるデータのフィルタリングを行う。
：結果
Count: クオリティのフィルターをパスしたリードの数
Coverage: クオリティのフィルターをパスしたリードの数
Frequency: 変異が見られた頻度
Forward reads: その領域に見られたForwardリードの数
Reverse reads:その領域に見られたReverseリードの数
Forward/reverse: Forward/Total reads または Reverse/Total reads のうち小さい方の値。Forwardと Reverseが同じなら、0.5となる。
Average quality: 該当する領域の平均リードクオリティ。
詳細
• 確率モデル（
）を使った変異検出
与えられるリードから、そのポジションの
を推定
と推定した
が異なる場合、変異として結
果
A
A
A
T
T
C
?
?
: Site type (ex) A/A, A/T, A/C ... ?
Reference
詳細
A
B
A∩B
P(A)
P(B)
P(A∩B)
)
(
)
|
(
)
(
)
|
(
)
(
)
|
(
)
(
)
(
)
|
(
)
(
B
P
B
A
P
A
P
A
B
P
B
P
B
A
P
B
A
P
A
P
A
B
P
B
A
P
)
(
)
|
(
)
|
(
P
B
A
P
A
B
P
B
ベイズの定理
事前確率
Prior
)
(
)
(
)
|
(
)
|
(
R
P
S
P
S
R
P
R
S
P
Reads
:
type
Site
:
R
S
A
A
A
T
T
C
?
?
: Site type (ex) A/A, A/T, A/C ... ?
Reference
)
|
(
R
S
P
)
(
P
S
: Error Model を使って推定
: Genome Model を使って推定
詳細
•
–
が のとき、
の大部分は になると仮定し、初期の確率を以下のように設定
し、 アルゴリズムを使ってそれぞれの確率を推定する。
• アルゴリズム（ ）は、得られたデータから推定したい現象が観察できない場合に、 その確率を推定する、一般的な統計の手法。Site Type Initial Probability
A/A 0.2475 A/C 0.001 A/G 0.001 A/T 0.001 T/C 0.001 T/G 0.001 T/T 0.2475 G/C 0.001 C/C 0.2475 G/G 0.2475 G/- 0.001 A/- 0.001 C/- 0.001
詳細
•
– リードに含まれるエラーを考慮するため、尤度のところにエラーを考慮した確
率を推定する。初期値を以下のように設定し、 アルゴリズムにて確率を推
定する。
Reference
A
C
G
T
-
Reads
A
0.90
0.025
0.025
0.025
0.025
C
0.025
0.90
0.025
0.025
0.025
G
0.025
0.025
0.90
0.025
0.025
T
0.025
0.025
0.025
0.90
0.025
-
0.025
0.025
0.025
0.025
0.90
詳細
変異コール
•
モデルと
モデルにより事後確率が計算できました。この時、リ
ファレンスと同じアレルである場合も計算されます。
•
： －＞
と考えます。
の事後確率が
と計算できたとし
ます。
• ウィザード中のパラメータで、 参照配列と異なる確率 を指定しています。
これを
とすると、
の確率は ％以下であるということになります。
•
の確率が ％という事は、指定した閾値を満たさないため、このポジショ
ンは変異としてコールされません。
A
Reference
それぞれの事後確率
A/A = 0.15
A/T = 0.8
?
詳細
変異コール
•
参照配列と異なる確率を ％とすると、
が ％の場合、そのポジションは変異
があるとされ、リファレンスと異なるアレル（
）のうち、最も事後確率が高いも
のを変異のアレルとして返します
。
A
Reference
それぞれの事後確率
A/A = 0.15
A/T = 0.8
A/C = 0.6
A/G = 0.01 .. etc.
?
活用例
• アプリケーションノート
健常デンマーク人２００名とデンマーク
人結腸癌患者の比較による癌体細胞
変異の検出
健常者解析フロー
マッピング
変異検出
アミノ酸置換
クオリティの悪いものを除去
Local Realignment
癌患者解析フロー
マッピング
変異検出
アミノ酸置換
各種フィルタリングLocal Realignment
体細胞変異を除く
原理
•
では
グラフというネットワーク理論に基づいた方
法で
アセンブリを実行します。
• 各リードからさらに短い長さの配列のセットを作成し、グラフを作成。
•
を利用しているオープンソースの方法では
が有名です。
ライブラリ配列
リード
Word セット
原理
•
グラフではリードを短い配列に分断し（
）、グラフを作成します。
（例） リード長
の場合は 個の
ができる。
リード AGTTGATCTTACTAGAGGAA 1 AGTTGATCTT 2 GTTGATCTTA 3 TTGATCTTAC 4 TGATCTTACT 5 GATCTTACTA 6 ATCTTACTAG 7 TCTTACTAGA 8 CTTACTAGAG 9 TTACTAGAGG 10 TACTAGAGGA 11 ACTAGAGGAA原理
• グラフ作成 簡単な例として
で考える
AACGT ACGTC CGTCA GTCAA TCAAGAACGT – ACGTC – CGTCA – GTCAA - TCAAG
AACGT ACGTC CGTCA CGTCG
GTCAA
CGTCA - GTCAA – TCAAG
AACGT – ACGTC
CGTCG
AACGTCAAG
AACGTCAAG
原理
CGTCA - GTCAA – TCAAG - CAAGT - AAGTC
AACGT – ACGTC AGTCC - GTCCA CGTCG - GTCGA - TCGAG - CGAGT - GAGTC
CGTCA - GTCAA – TCAAG AACGT – ACGTC
CGTCG
