RNA-SeqデータをSubio Platformに取り込むときには、 入力データの種類に応じて、適切な正規化・前処理シナリオを選ぶ必要があります。 特に、Gene Counts、TPM、FPKM、RPKMでは、値の意味やノイズの見え方が異なるため、 同じ前処理をそのまま適用することはできません。
発現差解析ではGene Countsを起点にするのが基本です。 また、PCA、クラスタリング、ヒートマップなどでデータ構造を確認する場合には、 Gene Countsを適切に正規化・対数変換・フィルタリングしたデータを使うことで、 サンプル間の関係や発現パターンを確認しやすくなります。 GO・Pathway解析では、発現差解析で得られた遺伝子リストを用いて、 生物学的な意味づけを行います。
しかし、公開データや過去の解析結果では、 Gene Countsが手に入らない、あるいは非常に手間がかかることもあります。 そのため、TPM, FPKM, RPKMを使う場合の注意点についても解説します。
正規化と前処理の適切なやり方を見つける。
Subio Platformでは、下記のプリセットのシナリオを用意しています。データに合わせて選択し、それをベースにデータの特徴に合わせて調整していくと簡単にできます。
- Expression Microarray
- RNA-Seq (Count)
- RNA-Seq (FPKM, TPM, RPKM)
- Methylation Beta Values
- Pre-normalized Log2 Data - すでに正規化・対数比化までの処理が終わっている場合に選択します。
- Nothing - なにもしない状態に戻すときに使います。
ただ、解析の経験がないと適切な処理ができているかどうか判断するのが難しいかもしれません。 オンライントレーニング では、あなたが解析したいデータでインポート・正規化・解析までひととおり学んでいただけます。
前提となるRNA-Seq Gene Counts データの特徴
RNA-SeqのGene Counts のテーブルを見てみると、たくさんの 0 があることに気づくでしょう。 0 は「マッピングされたリードが無い」ことを示していますが、 それを単純に「発現していない」と読み替えることはできません。 これが、RNA-Seqのデータ解析の難しいところの一つです。
その理由は、繰返しのある実験データを見ると明らかです。 例えば、繰返しのうちあるサンプルでは値があるのに、別のサンプルでは0ということがあります。 これは、発現していてもシーケンサーにキャプチャーされない可能性があることを示しています。 発現している遺伝子のCountが偶然 0 となるかどうかは確率的に決まり、 発現が低いほど 0 となる確率が高くなることが予想されます。 実際にデータを見ていると、繰り返しサンプルの一部のCountが0となる遺伝子は、 低Counts領域に集中していることが分かります。
少し単純化して考えてみましょう。グループAとグループBの平均値を比較して、 ある遺伝子がグループBで発現量が下がって見えたとします。 しかし、よく見るとグループBにはいくつかCountが 0 のサンプルが含まれていて、 測定されたCountだけ見るとグループAとあまり変わらない、ということがありえます。 この場合、0が偶然測定されなかったのであれば、 実際の発現量は見かけほど下がっていないとみなすほうが正しいかもしれません。 RNA-Seqデータの0の扱いは、非常に難しいのです。 これを踏まえて、データ解析をどうするかを検討してみましょう。
Gene Counts データを解析する場合
繰り返しサンプルがある場合、そのCountを散布図で見てみると、 Countの高い領域では対角線上に点が収束していて、 低い領域では収束していないように見えます。 これがシグナル領域とノイズ領域です。 シグナル領域とノイズ領域の境界は明確ではありません。 また、データセットによってもまちまちですが、 十分に多いバルクRNA-Seqデータであれば、 だいたいGene Countsが100より高いところがシグナル領域、 10より低いところがノイズ領域、そして10~100の間に シグナルとノイズの境界が見られることが多いです。 そして、シグナルとノイズの境界は明確に切り替わるのではなく、 グラデーション状です。 上記のとおり、Countsが偶然0になっている遺伝子が多いのは、 ノイズ領域です。（黒点）
そこで、Subio Platformでは、 RNA-Seq (Count) というプリセットのノーマライズシナリオを用意して、 Countの値に対して次のような前処理・正規化を行います。
- Log Transformation: base 2
- Global Normalization: at 90th percentile
- Low Signal Cutoff: Replace <30 to 30
- Fill Missing Values: 4
- Centering
ステップ１で、Count 0は欠損値に置き換えられます。
ステップ２は、Global Normalizationを適用しています。 これは、サンプル間の全体的なスケールの違いを補正するための処理です。 総リード数で補正するCPMや、組成の偏りを考慮するTMM正規化とは計算方法が異なりますが、 サンプル間の全体的なスケール差を補正するという目的は共通していますので、 適用後の効果はだいたい似ています。
しかし実際には、100万で割っても正規化が中途半端なことがあるので、 たとえそのような正規化が適用されているデータに対しても、 Global Normalizationを重ねて適用するほうが安全です。 デフォルトで90となっていますが、これはデータを見て調整してください。 ヒストグラムで見ると、右側にある山の頂点付近の値を設定します。 75th percentileの位置を示す縦線を参考に値を決めてください。
ステップ３は、正規化後に値の下限値が、 ダイナミックレンジの違いを反映してズレるため、 これを揃えるための処理です。 実際のRNA-Seqのデータセットでは、 各サンプルの総リード数に数倍の差があることが普通で、 正規化後の値の分布を見ると、下限値側にズレが見られます。 これがマイクロアレイデータの正規化と大きく異なるポイントです。 ダイナミックレンジの狭いサンプルが、偶然ある群に多ければ、 その群では下限値がせり上がり、見せかけの発現上昇と見えてしまうので、 これを防ぐ必要があるのです。
そして、下限値を設定するときのポイントは、 最もダイナミックレンジの狭いサンプル群のシグナル領域の下限付近に設定することです。 ダイナミックレンジの広いサンプルがあるのにもったいないと思うかもしれませんが、 解析全体の信頼性を考えると、最も狭いダイナミックレンジに合わせて 信頼できる範囲をそろえる必要があります。 そうしなければ、狭いダイナミックレンジの信頼できない測定値を 解析に含めることになってしまうからです。
edgeRやDESeq2は、Gene Countsの分散構造を考慮して発現差を検定するための統計モデルです。 しかし、サンプルごとのダイナミックレンジの違いによって生じる 低Count領域の測定可能な下限値の違いをそろえるという考え方は、 これらの基本的なモデルには組み込まれていません。
そのため、低Count領域に多くの遺伝子が残ったまま解析すると、 本来は測定限界やサンプル間のダイナミックレンジの違いとして扱うべき変動が、 統計的な発現差として検出されることがあります。 実際のRNA-Seqデータでは、この影響によって、 edgeRやDESeq2で抽出された遺伝子リストに、 低Count由来の不安定な遺伝子が混入することがあります。
このため、edgeRやDESeq2を使う場合でも、 正規化後の下限値のズレに由来する遺伝子が 有意遺伝子リストに紛れていないかを確認することが重要です。
ステップ４は、Count 0だったところが欠損値になっているので、 それを何かの値に置き換えます。 ここにはステップ３と同じ値か、それよりも少し低い値を設定するといいでしょう。 ただし、指数で指定することに注意してください。 ４は、２の４乗＝16の意味です。ステップ３で設定した下限値と比べて、 あまり小さくしない方がいいです。その理由は、Gene Countsが0という意味が、 本当にステップ３で設定した下限値と比べて十分に低い保証がないからです。 であれば、よくわからない低い値を与えるよりも、下限値と比べて あまり変化しないように見えるようにするほうが安全だからです。
ステップ５は、全サンプルの平均値に対する比にするものです。 このステップは、コントロールがある場合は、Ratio to Control Samplesブロックに置き換えることができます。
以上の計算をエクセルで再現したファイルを用意しましたので、詳しくお知りになりたい方はダウンロードしてみてください。
TPM・FPKM・RPKMデータを解析する場合
TPM・FPKM・RPKMデータを解析するのは、Gene Counts の場合よりも厄介です。 なぜなら、ノイズ領域とシグナル領域が混じりあっているため、 0を何らかの値に置換すると、冒頭の偽の発現差を検出する問題が顕在化してしまうからです。 この低Counts由来の信頼できないばらつきは、 正規化や補正によって後から解消できるものではありません。 したがって、方針として次のどちらを取るかを選択する必要があります。
- 低Counts由来のばらつきを持つ遺伝子群をできるだけ排除する。ただし、解析可能な遺伝子が少なくなることを受け入れる。
- 解析可能な遺伝子をできるだけ残すことを優先する。その代わり、解析結果に信頼できない測定値に由来するノイズが含まれることを受け入れる。
それでは、Subio Platformのプリセット・ノーマライズ・シナリオ RNA-Seq (TPM, FPKM, RPKM) を見てみましょう。
- Log Transformation: base 2
- Global Normalization: at 90th percentile
- Low Signal Cutoff: Replace <0.01 to 0.01
- Centering
ステップ2まではRNA-Seq (Count) のシナリオと同じですので、そちらの説明を参照してください。
ステップ３は、下限値を設定するものですが、
Gene Counts の場合と異なり、ノイズ領域の影響をキャンセルすることはできません。
上記の二択の解析方針に従って、高い閾値を設定するか、低い閾値を設定するか、どちらかになります。
Gene Counts の場合には、信頼できない測定値が低Counts領域に固まっているために対処が可能ですが、
TPM・FPKM・RPKMの場合は、低Counts由来の信頼できない測定値が広く拡散しているため、
Fill Missing Valuesを使って適切に置換することができません。
そのため、TPM・FPKM・RPKMで log 変換後に欠損値となった値は、
無理に別の値へ置き換えず、欠損値として残すほうが安全です。
発現差解析をする際には、通常のFold値やP値を用いる方法だけでなく、
グループAでは値が無く、グループBでは値がある遺伝子群、
あるいはその逆の遺伝子群も加える必要があります。
この分手間がかかりますので、忘れないようにしてください。
ステップ４の説明も上記Centeringと同じですが、 注意点は、Controlサンプルで測定値が欠損値となっている可能性があるので、 Ratio to Control Samplesブロックへの置き換えはしないほうが良いということです。
ただ、この状態でクラスタリングなどの多変量解析を行うと、欠損値が邪魔になります。 そこで、Centeringの後にFill Missing Valuesを追加して、 -1など発現減少を表す何らかの値に置き換えるのはありでしょう。 ということは、発現差解析のときと多変量解析のときで、 Fill Missing Valuesの有無を切り替えることが必要となりますので、これも手間がかかります。