Top of Technical Information

マイクロアレイより少し複雑な、RNA-Seqのデータ解析手順

RNA-Seqでトランスクリプトームの測定を行うと、結果としてcountまたはFPKMのテーブルが得られます。これをエクセルで見てみると、0あるいは値が欠けているセルがたくさんあることがわかります。これらは「マッピングされたリードが無い」ことを表していますが、これが即ちその遺伝子が「発現していない」ことを表しているとは言えないのが、RNA-Seqのデータ解析の難しさの一つです。

その理由は、繰返し実験があると明らかです。例えば、2回の繰返しがあって、1回目のサンプルでは値があるのに、2回目では0ということがあります。一方のサンプルではcountの値がシグナル領域にあっても、もう片方では0となることもあります。これは、countが0でも「発現していない」とは言い切れず、「本当は発現しているにも関わらずシーケンサーにキャプチャーされなかった」可能性があることを示唆しています。データをよく見ると、発現量が少ないほど「キャプチャーされない」可能性が高く、発現量が高い遺伝子ほど「キャプチャーされる」可能性が高いというグラデーションがあることがわかります。

RNA-Seqでは、発現量の非常に少ない遺伝子は、数十回またはそれ以上の繰返しのうち数回で値があるように見えます。これはRNA-Seqがデジタルな手法であることに起因する現象で、マイクロアレイのようなアナログな手法では起こりません。非常に低い発現量のRNAを全て観測しようとしたら、アジレントのマイクロアレイであれば2,3回の測定で済むことが、RNA-Seqでは数十~数百回測定しなければならないのです。

いずれにせよ、RNA-Seqデータのこのような特徴を踏まえると、あるグループで平均値をとる場合、値があるものと0が混在していれば、0を含めて計算すると、そのグループにおいて発現が不当に下がって見えることになります。0の影響を除いてから平均値の計算や、検定を行う必要があります。

Subio Platformでは、このような問題に配慮して解析をするために、RNA-Seqのデータ解析では、下記のNormalization Scenarioを提案します。プリセットのNormalize Scenarioで、「Profiling (GDC RNA-Seq counts)」を選択すると、自動的にNormalize Blockは配置されます。ただし、各ブロックのオプション設定は、データに合わせて調整する必要があります。どのようなRNA-Seqデータにも適合するわけではありませんのでご注意ください。

  1. Log Transformation
  2. Global Normalization
  3. Low Signal Cutoff
  4. Centering あるいは Ratio to Control Samples

ノーマライズシナリオの解説

最初にLog Transformationブロックを持ってくることで、countあるいはFPKMが0だった場合、これらは欠損値になります。こうしておくことで、その後グループ間で比較するときに、欠損値は自動的に排除されて計算されることになります。

次のGlobal Normalizationブロックについては、CountもFPKMも、すでにノーマライズされている場合がほとんどですので、無くてもいい場合は外してください。しかしヒストグラムを見て、正規化が必要と判断される場合は適用してください。適用する場合は、ヒストグラムを見ながらシグナル領域の山のピークが揃うように、基準値を75 th ~ 90th percentileあたりで調整してください。シグナル領域にある遺伝子が少ないほど、揃える基準は高いところになります。インプットのRNAがごく微量の場合は、99th percentileという極端な値になることもあります。

Low Signal Cutoffブロックは、ノイズ領域で発現変動が不当に大きく見えるのを抑制する目的で使います。ノイズ領域とシグナル領域の境界あたりからその半分の値までの範囲で閾値を調整するとよいでしょう。

最後のブロックは、比にするものです。コントロールサンプルが無い場合はCenteringブロックを選び、コントロールサンプルがある場合はRatio to Control Samplesブロックを選んでください。

発現差のある遺伝子の集め方の注意点

ところで、このような前処理で生成されたProcessed Signalは、分母または分子の値が欠損値の場合には計算不能のため、欠損値になります。そうすると、fold値やT検定のP値も計算不可能となります。しかし、たとえばグループAのサンプルではすべてcountが0だったのに対し、グループBではある程度信頼できるcountの値があるという場合は、グループBで発現が上昇していると言えます。このような遺伝子群は、フィルターツールを使って、グループAでは100%のサンプルでCh1RawSignalが0の遺伝子群と、グループBでは50%以上のサンプルでCh1RawSignalがたとえば10より大きい遺伝子群を集め、この二つをVenn Diagramツールを使って掛け合わせて抽出してください。

RNA-Seq のデータ解析手順

データインポートから解析まで一通り学びたい方は、オンライントレーニングをお申し込みください。

クラスタリング解析

発現差解析は上記の前処理でいいのですが、クラスタリングには欠損値が邪魔になります。そこで、上記の前処理で発現差解析を行った後、欠損値の影響を避けたい場合は、Fill Missing Valuesブロックを追加してProcessed Signalを計算しなおしてから、クラスタリングを実行するほうがいいでしょう。おすすめの使い方は次の2通りです。

最後がCenteringになっている場合

Fill Missing Valuesブロックを 最後に 持ってきます。欠損値と入れ替える値は、-1(log2のスケールにおいて1/2に減少を表す)か-2にするといいでしょう。

最後がRatio to Control Samples になっている場合

F ill Missing ValuesブロックをRatio to Controls Samplesブロックの 前に 持ってきます。欠損値と入れ替える値は、Low Signal Cutoffで用いた閾値か、それよし少し低い値でいいでしょう。

ノーマライズ: Fill Missing Values ブロックの使い方」も合わせてご覧ください。