16S rRNA解析のクオリティフィルタリングについて

最近は仕事で微生物の 16S rRNA解析用のツールを作っています。
そのツールには入力データとして、NGS で読んだリードデータを入れるのですが、その前にリードのフィルタリングを行わなければなりません。
NGS のデータにはシーケンスエラーが含まれていて、それをなるべく除去しておかないと、解析結果に影響が出てしまいます。
というわけで、シーケンスのクオリティ情報をもとにフィルタリングを行うんですが、これをクオリティフィルタリングと言います。*1

ところで、アセンブリやらマッピングやらの場合と、16S rRNA解析の場合とでは、クオリティフィルタリングの意味合いが異なってきます。
アセンブリなどの場合、シーケンスエラーが多少あっても、リードがたくさんあった方がコンセンサス配列に厚みができ、結果としてコンティグの品質は上がります。
一方、16S rRNA解析の場合、一本のリードがひとつの微生物に対応し、数塩基の違いが生物種の違いとなるため、リードひとつひとつのクオリティは高くなければなりません。
というわけで、16S rRNA解析では、通常の NGS 解析より厳しいクオリティでのフィルタリングが必要とされます。

Reducing the Effects of PCR Amplification and Sequencing Artifacts on 16S rRNA-Based Studies では、Human Microbiome Project のデータを使って、様々なフィルタリング手法を比較しています。
今回はこれらのフィルタリング手法とその評価を紹介してみたいと思います。

フィルタリング手法の評価指標

フィルタリング手法の評価指標としては、エラーレートとフィルタリング後に残ったリード数の二つを用います。
フィルタリングとして望ましいのは、結果として残ったリードの「エラー率はなるべく低く」かつ「リード数はなるべく多く残る」ような手法です。
ちなみにフィルタを全く適用していない状態では、エラー率 0.0061、リード数は 100% 残っている状態です。

	エラー率	リード数[%]
Raw	0.0061	100

これを基準として順番に見ていきましょう。

Basicフィルタ

まずは、基本的なフィルタリングを適用してみます。下記の表をご覧ください。

バーコード配列のミスマッチ	0〜1	2
	0.0056	0.0101
プライマー配列のミスマッチ	0〜2	3
	0.0056	0.0073
Ambiguous base call(N)	なし	あり
	0.0058	0.0111
リード長	200bp以上	200bp未満
	0.0059	0.0163
Homopolymer*2 の長さ	8未満	8以上
	0.0061	0.0663

これは、リードを二つのグループに分け、それぞれのグループでのエラー率を算出しています。
左側のグループと右側のグループでは、右側のグループの方が有意に高いエラー率となっています。
なので、基本的なフィルタリングとして、右側のグループを削除してみます。
これにより、全体のエラー率は 0.0056、リード数は 84.1% となります。

	エラー率	リード数[%]
Raw	0.0061	100
Basic	0.0056	84.1

First 250bp フィルタ

次に、単純でばかげたフィルタを考えてみます。
シーケンスのクオリティは、リードの読み始めの部分では高いのですが、読み進んでいくとだんだん低くなります。
なので、最初の高いほうだけ抜き出せば、エラー率を低く抑えることができます。
First 250bp フィルタは、リードの最初の 250bp のみを抜き出すフィルタです。
このフィルタの結果は次のようになります。

	エラー率	リード数[%]
Raw	0.0061	100
Basic	0.0056	84.1
First 250bp	0.0017	≒90

単純なわりにはエラー率が大幅に改善されています。
ただ、シーケンスのクオリティスコアを考慮していないので、クオリティの高い部分を捨ててしまうことや、逆に低いのに通してしまうという欠点があり、実用的ではありません。*3

Hard Cutoff フィルタ

というわけで、ここからはシーケンスのクオリティスコアを考慮した手法を考えていきたいと思います。
Hard Cutoff フィルタは、リードの先頭からクオリティスコアを見ていって、閾値よりも低くならない範囲を採用するフィルタです。
クオリティスコアの閾値を 25 とした場合の結果を次に示します。

	エラー率	リード数[%]
Raw	0.0061	100
Basic	0.0056	84.1
First 250bp	0.0017	≒90
Hard Cutoff(Q25)	0.0007	36.5

大幅にエラー率が下がりましたが、リード数も大幅に減ってしまいました。*4
この手法の良い点は、実装が簡単であり、なおかつエラー率を大幅に改善できるところです。
NGS はリードを大量に出力するので、リードの減少が気にならないというのであれば、この方法を採用することも考えられます。

Rolling Average フィルタ

Hard Cutoff フィルタは、採用する範囲すべてにおいて高いクオリティスコアを要求する厳しいフィルタです。
しかし、実際のクオリティスコアは単調減少するわけではなく、多少のゆれがあります。
クオリティスコアが一時的に閾値より低くなったとしても、その先に高い状態が続くということもまれではありません。
この部分を捨てるのはもったいないということで、Rolling Average フィルタは、リードの先頭からのクオリティの平均値を計算し、それが閾値より低くならない範囲を採用します。

	エラー率	リード数[%]
Raw	0.0061	100
Basic	0.0056	84.1
First 250bp	0.0017	≒90
Hard Cutoff(Q25)	0.0007	36.5
Rolling Average(Q35)	≒0.0030	≒50

この手法は、それほどエラー率が改善されないわりにはリード数が減少してしまいます。
この手法はぜんぜん実用的ではありませんが、アイデアは次に説明する Sliding Window に引き継がれます。

Sliding Window フィルタ

Rolling Average では、リードの先頭からのクオリティスコアの平均値を算出しました。
しかし、リードの先頭はクオリティが高く、その影響がのちのちまで残ってしまうと、クオリティの低い部分を採用してしまうことにつながります。
そこで、50bp〜100bp ぐらいの着目領域（Window）を考え、リードの先頭からスライドさせていきます。
スライドするごとに Window 内だけでのクオリティスコアの平均値を計算します。
その平均値が閾値を下回ったとすると、そこはクオリティスコアがいっせいに下がったポイントであると分かるわけです。
これが Sliding Window フィルタのアイデアです。
結果は次のようになります。

	エラー率	リード数[%]
Raw	0.0061	100
Basic	0.0056	84.1
First 250bp	0.0017	≒90
Hard Cutoff(Q25)	0.0007	36.5
Rolling Average(Q35)	≒0.0030	≒50
Sliding Window(Q35,50bp)	0.0010	60〜65

エラー率の改善もいいし、リードも結構残っています。
実装も簡単なので、Hard Cutoff ではリードが減りすぎるという場合に使えそうです。

shhh.flows

最後にこの論文の著者らが開発した手法(shhh.flows)の結果だけ。

	エラー率	リード数[%]
Raw	0.0061	100
Basic	0.0056	84.1
First 250bp	0.0017	≒90
Hard Cutoff(Q25)	0.0007	36.5
Rolling Average(Q35)	≒0.0030	≒50
Sliding Window(Q35,50bp)	0.0010	60〜65
shhh.flows(350flows)	0.0006	≒80

詳しく見ていないので、詳細は分からないのですが、非常に良いです。
もちろん著者らの主張はこれ(shhh.flows)がいいと言っています。
著者らの開発した mothur というツールに入っているそうです。
時間があったらもう少し調べてみようと思っています。

まとめ

16S rRNA解析用のクオリティフィルタリング手法とその評価を見てきました。
16S ではフィルタリング手法は解析結果に大きな影響を与えるので重要です。
現状、Hard Cutoff してリードが足りない場合は OTU のクラスタ数を減らしてごまかすみたいなことをやってるので、ここらへんはもうちょっと調べたいところです。