RNA-Seq データ解析のための、正規化のプリセット・シナリオの解説

RNA-Seqのリードカウント(Count)のテーブルを見てみると、たくさんの 0 があることに気づくでしょう。0 は「マッピングされたリードが無い」ことを示していますが、即ちその遺伝子が「発現していない」とは言えないのが、RNA-Seqのデータ解析の難しいところです。

その理由は、繰返しのある実験データを見ると明らかです。例えば、繰返しのうちあるサンプルでは値があるのに、別のサンプルでは0ということがあります。これは、発現していてもシーケンサーにキャプチャーされない可能性があることを示しています。発現している遺伝子のCountが偶然 0 となるかどうかは確率的に決まり、発現が低いほど 0 となる確率が高くなることが予想されます。実際にデータを見ていると、繰り返しサンプルの一部のCountが0となる遺伝子は、ノイズ領域に集中していることが分かります。

ここで、グループAとグループBの平均値を比較して、ある遺伝子がグループBで発現量が下がって見えたとします。しかし、よく見るとグループBのサンプルにはCountが 0 のサンプルが含まれていて、測定されたCountだけ見るとグループAとあまり変わらない、ということがありえます。「Count 0」の扱い方がまずいと、偽の発現差を検出してしまうだろうということです。これを踏まえて、データ解析をどうするかを検討してみましょう。

Countデータを解析する場合

繰り返しサンプルがある場合、そのCountを散布図で見てみると、Countの高い領域では対角線上に点が収束していて、低い領域では収束していないように見えます。これがシグナル領域とノイズ領域です。シグナル領域とノイズ領域の境界は明確ではありません。また、データセットによってもまちまちですが、だいたいCountの値が100より高いところがシグナル領域、10より低いところがノイズ領域、そして10~100の間にシグナルとノイズの境界が見られることが多いです。上記のとおり、Countsが偶然0になっている遺伝子が多いのは、ノイズ領域です。(黒点)

the dynamic range of RNA-Seq - fig1

そこで、Subio Platformでは、RNA-Seq (Count) というプリセットのノーマライズシナリオを用意して、Countの値に対して次のような前処理・正規化を行います。

  1. Log Transformation: base 2
  2. Global Normalization: at 90th percentile
  3. Low Signal Cutoff: Replace <30 to 30
  4. Fill Missing Values: 4
  5. Centering

ステップ1で、Count 0は欠損値に置き換えられます。

ステップ2のデフォルトで90となっているpercentileの値は、ダイナミックレンジが広いほど低く、狭いほど高く調整してください。TPMはFASTQファイル処理の中で1kbあたりの正規化と、100万リードあたりの正規化が施された値です。そしてCountは、TPMから1kbあたりの正規化をキャンセルするよう逆算して得た値を代用しています。言い換えると、100万リードあたりの正規化は組み込まれており、本来Global Normalizationを重ねる必要はありません。しかし、私たちの経験からGlobal Normalizationを入れたほうが良い場合が多いので、シナリオに入れています。ヒストグラムを見て必要ないと判断される場合は除いて構いません。

ステップ3は、ノイズ領域の測定値が偽の発現差を検出するのを防ぐために、Countの下限値をシグナル領域の下限に合わせる措置です。

ステップ4は、Count 0だったところが欠損値になっているので、それを何かの値に置き換えます。ステップ3と同じ値か、それよりも少し低い値を設定するといいでしょう。ただし、指数で指定することに注意してください。4は、2の4乗=16の意味です。ステップ3で設定したデフォルトの下限値の約半分に相当します。

ステップ5は、全サンプルの平均値に対する比にするものです。このステップは、コントロールがある場合は、Ratio to Control Samplesブロックに置き換えることができます。

TPM・FPKM・RPKMデータを解析する場合

TPM・FPKM・RPKMデータを解析するのは、Countの場合よりも厄介です。なぜなら、ノイズ領域がシグナル領域混じりあって広がっている(黒点)ので、0を何らかの値に置換すると、冒頭の偽の発現差を検出する問題が顕在化してしまうからです。

the dynamic range of RNA-Seq - fig2

それでは、Subio Platformのプリセット・ノーマライズ・シナリオ RNA-Seq (TPM, FPKM, RPKM) を見てみましょう。

  1. Log Transformation: base 2
  2. Global Normalization: at 90th percentile
  3. Low Signal Cutoff: Replace <0.01 to 0.01
  4. Centering

ステップ2まではRNA-Seq (Count) のシナリオと同じですので、そちらの説明を参照してください。

ステップ3は、下限値を設定するものですが、Countの場合と異なり、ノイズ領域の影響をキャンセルすることはできません。下限値を設定する意味は最低限でしかありません。
Countの場合にはここでFill Missing Valuesブロックを適用していましたが、TPM・FPKM・RPKMの場合は使わないほうが良いでしょう。発現差解析をする際には、通常のFold値やP値を用いる方法だけでなく、グループAでは値が無いけどグループBでは値がある遺伝子群、あるいはその逆の遺伝子群も加える必要があります。この分手間がかかりますので、忘れないようにしてください。

ステップ4の説明も上記Centeringと同じですが、注意点は、Controlサンプルで測定値が欠損値となっている可能性があるので、Ratio to Control Samplesブロックへの置き換えはしないほうが良いということです。

ただ、この状態でクラスタリングなどの多変量解析を行うと、欠損値が邪魔になります。そこで、Centeringの後にFill Missing Valuesを追加して、-1など発現減少を表す何らかの値に置き換えるのはありでしょう。ということは、発現差解析のときと多変量解析のときで、Fill Missing Valuesの有無を切り替えることが必要となりますので、これも手間がかかります。

RNA-Seq データ解析 (TPM バージョン)

RNA-Seqのデータ解析は、Countsを使うよりTPMを使うほうが複雑になります。Countsを使ったデータ解析例もご覧ください。