1. サポート
2. Help - 解析理論とケーススタディ
3. RNA-SeqのGene Countsは比較できる絶対値ではない｜GEOデータ統合とバッチエフェクトの対処法

GEOには、多数のRNA-Seqやマイクロアレイデータが公開されています。 これらのデータを複数組み合わせれば、より多くのサンプルを使った解析ができるように思えるかもしれません。

しかし、実際には、複数のGEO Seriesをそのまま統合して解析することは、多くの場合うまくいきません。 疾患群と対照群の違いではなく、Seriesの違いがそのままクラスターの違いとして現れることが多いからです。

この問題は、GEOデータの再解析だけに限った話ではありません。 自分たちで取得するRNA-Seqデータでも、今年、来年、再来年と時間をあけて測定したデータが、 そのまま比較可能である保証はありません。

この記事では、複数のGEO Seriesを例に、 発現値そのものを統合するのではなく、 各データセット内で定義された比較を統合するという考え方を紹介します。 さらに、この考え方が長期にわたる前向き研究の実験計画にどのような示唆を与えるかを考えます。

RNA-Seqで得られるGene Countsは、比較可能な絶対値ではない

RNA-Seq解析では、各遺伝子にどれだけのリードが割り当てられたかを表す Gene Countsが得られます。 そのため、発現データ解析に慣れていない方にとっては、 「カウント値なのだから、異なるデータセット間でもそのまま比較できる絶対値ではないか」 と見えるかもしれません。

しかし、これはRNA-Seqデータ解析で最初につまずきやすい、大きな誤解です。 Gene Countsは、単純に「サンプル中にそのRNAが何個あったか」を表す値ではありません。 同じサンプルを測定していても、測定条件やデータ処理の違いによって、 同じ遺伝子のカウント値が大きく変わることがあります。
この違いは必ずしも「実験が下手だった」ために生じるわけではありません。 熟練した技術者が、同じプロトコルに従って丁寧に測定しても、 測定バッチ、試薬ロット、機器、ライブラリ調製、データ処理のわずかな違いによって、 カウント値のスケールや分布が変わることがあります。 これはRNA-Seqデータを扱ううえで、最初から想定しておくべき性質です。

したがって、異なるGEO Seriesに含まれるサンプル同士を、 Gene Countsや正規化後の発現値として直接比較できる可能性は極めて低いと考えるべきです。 これは、RNA-Seqデータ解析では基本的な前提です。
一方で、同一Series内でコントロールに対する比を取れば、 Series間でも比較できる場合があります。 発現値そのものは直接比較できないのに、比に変換すると比較できる場合がある。 ここがRNA-Seqデータのおもしろい性質です。 この記事では、その考え方を実際のGEOデータで確認します。

複数のGEO Seriesをそのまま統合すると何が起きるか

今回、Atopic Dermatitis（AD）の皮膚サンプルを対象とした複数のGEO Seriesを使いました。 RNA-Seqデータとして、GSE65832、GSE121212、GSE140227を使用しました。 これらはいずれも、AD患者から得られたnon-lesional skinとlesional skinのデータを含んでいます。

Gene Countsは、すべてGEO Generatedのデータを使用しました。 同じGEO Generated Gene Countsを使うことで、少なくともデータ処理パイプラインの違いはある程度そろえられます。 しかし、それでも各Seriesの発現値分布は大きく異なっていました。 シーケンサー、サンプル状態、read depthなどの違いにより、 発現値の分布は大きく変わっていたと考えられます。

特にGSE65832は、古いIllumina Genome Analyzer IIxで取得されたデータであり、 さらにFFPEサンプルに由来すると考えられます。 このSeriesでは、GSE121212やGSE140227に見られるような、 低カウント領域と高カウント領域に分かれる典型的な二山型の分布が明瞭ではありません。 この時点で、GSE65832を他のSeriesと発現値のまま統合することは 非常に難しいと予想されます。

実際に、これらのSeriesをそのまま統合してクラスタリングすると、 lesional / non-lesionalの違いよりも、Series間の違いが強く反映されました。 特にGSE65832は、予想どおり他の2つのSeriesから大きく離れており、 ほとんど逆転した発現プロファイルとして見えています。 同じRNA-Seqデータであっても、異なるSeriesを単純に結合すると、 生物学的な違いではなく、Series間の技術的な違いが強く反映されることを示しています。

「異なるSeriesのサンプル同士を直接比較できない」と言っても、 言葉だけでは分かりにくかったのではないかと思います。 しかし、この図を見ると、その意味が直感的に分かるのではないでしょうか。 この状態では、生物学的な違いよりも、Series間の技術的な違いの方が はるかに強くサンプル間の違いに反映されています。

各Series内のnon-lesionalに対する比にすると、Seriesを越えた比較が見えてくる

次に、各Series内のnon-lesionalサンプルの平均値を基準として、 各サンプルをlog2 ratioに変換しました。 つまり、GSE65832、GSE121212、GSE140227をまとめて一つの平均値で補正するのではなく、 それぞれのSeriesの中で測定されたnon-lesionalサンプルを基準にしました。

発現値のまま統合した場合と比べると、Seriesごとの分離は弱まりました。 さらに、non-lesionalに対してlesionalまたはchronic lesionでどのような遺伝子の発現が上昇し、 どのような遺伝子の発現が低下するのかについて、 Seriesを越えてある程度共通したパターンが見えるようになりました。

特に注目したいのはGSE65832です。 先ほどのクラスタリングでは、GSE65832は他の2つのSeriesから大きく離れ、 ほとんど逆転した発現プロファイルのように見えていました。 しかし、各Series内のnon-lesionalを基準にしたlog2 ratioでは、 GSE65832も他の2つのSeriesと似たような発現変化パターンとして見えるようになっています。 これは大きな進歩です。

もちろん、この方法でSeries間の違いが完全に消えるわけではありません。

ここからさらに、ComBatなどの補正アルゴリズムを使って、 Series間の違いをより小さく見せることもできます。 しかし、その結果が本来の生物学的状態に近づいているとは限りません。 見た目のクラスターがきれいになったとしても、 生物学的な差まで一緒に弱めてしまったり、 逆に人工的なパターンを作ってしまったりする可能性があります。

そのため、補正アルゴリズムを使う場合でも、 「きれいに混ざったから正しい」と判断するのではなく、 補正前後で何が変わったのかを確認しながら解釈する必要があります。

同じ患者のnon-lesionalを基準にすると、個体差を抑えられる

今回使用したSeriesでは、さらに同じ患者から得られたnon-lesionalサンプルと lesionalサンプル、またはchronic lesionサンプルのペアが含まれていました。 そこで次に、各患者ごとにnon-lesionalサンプルを基準として、 対応するlesionalまたはchronic lesionサンプルをlog2 ratioに変換しました。

この変換では、患者ごとの基礎的な発現プロファイルの違いが抑えられます。 そのため、患者間の個体差ではなく、 同じ患者内でnon-lesionalからlesional、またはchronic lesionへどのように変化したかに注目できます。

このクラスタリングでは、lesionalとchronic lesionのサンプルだけを表示しています。 各患者のnon-lesionalを基準にしたことで、 Series間の違いだけでなく、患者ごとの背景差もさらに抑えられ、 病変部に伴う発現変化をより直接的に比較できる状態になりました。

通常、複数のGEO Seriesを統合する場合には、 このように同じ患者から得られたnon-lesionalとlesionalのペアがそろっていることは多くありません。 しかし、もしペアがあれば、 患者ごとの基礎的な発現差を取り除き、 病変部に伴う変化だけに注目できます。 この結果を見ると、paired controlがいかに強力な補正手段になり得るかが分かります。

ペアサンプルがない場合には、各Series内に比較の基準として使えるコントロール群があるかが次の鍵になります。 ただし、複数Seriesで同じ意味を持つコントロール群がそろっているとは限りません。 そのようなコントロール群が存在すること自体が、統合解析では貴重な条件です。

マイクロアレイデータも、変化量としてなら同じ解析空間に置ける可能性がある

次に、Affymetrix HG-U133 Plus 2.0マイクロアレイデータであるGSE36842を追加しました。 GSE36842では、acute lesional skin（ALS）を「lesional」、 chronic lesional skin（CLS）を「chronic lesion」として扱いました。

マイクロアレイとRNA-Seqでは、発現値の単位、分布、ダイナミックレンジが大きく異なります。 そのため、マイクロアレイのシグナル値とRNA-SeqのGene CountsやCPMを、 発現値そのものとして一つの行列に統合することは適切ではありません。

しかし、各患者のnon-lesionalを基準にしてlesionalまたはchronic lesionへのlog2 ratioに変換すると、 比較対象は発現値そのものではなく、各患者内で病変部に向かってどのように変化したかになります。 この形であれば、RNA-Seqとマイクロアレイという測定技術の違いを越えて、 共通する発現変化を探索できる可能性があります。

実際に、GSE36842を追加してクラスタリングしても、 変動の方向はRNA-Seqデータとおおむね一致しているように見えます。 GSE36842はRMAにより数値化されたprocessed dataを用いているため、 RNA-Seq由来のlog2 ratioと比べると変動幅はかなり小さく見えます。 それでも、病変部で上昇する遺伝子群、低下する遺伝子群の大まかな方向は、 RNA-Seq側の結果と共通していることが分かります。

この性質は、RNA-Seqになってから新しく見つかったものではありません。 マイクロアレイ時代にも、異なるプラットフォームの発現値そのものを直接統合することは難しい一方で、 同じ比較から得られるfold changeやratioは、プラットフォームを越えて比較しやすい指標として扱われてきました。

実際、大規模な MicroArray Quality Control（MAQC）プロジェクト では、複数の参照RNAサンプルを、複数の施設・複数のマイクロアレイプラットフォームで測定し、 プラットフォーム内の再現性と、発現差遺伝子のプラットフォーム間一致性が検証されました。 この結果は、発現値そのものではなく、同じ比較から得られる変化量を見ることの有効性を示す代表的な例です。

今回、AffymetrixマイクロアレイデータをRNA-Seqデータに追加しても、 各患者のnon-lesionalに対するlesionalまたはchronic lesionのlog2 ratioとして表現すれば、 変動幅には違いがあるものの、発現変化の方向はおおむね共通して見えました。

測定技術がマイクロアレイからRNA-Seqへ変わっても、 「発現値そのものよりも、同じ比較から得られる変化量の方が再現しやすい」 というトランスクリプトミクスデータの性質は受け継がれていると考えられます。

発現値を統合するのではなく、比較を統合する

ここまでの結果から、複数データセットを統合するときに重要なのは、 発現値そのものを無理に同じ尺度にそろえることではないと分かります。

各データセット内で信頼できる比較を作り、 その比較どうしを統合する必要があります。 言い換えれば、統合すべきなのは発現値そのものではなく、 実験内で定義された変化です。

たとえば、次のような比較です。

• disease / control
• treatment / untreated
• after / before
• time X / time 0
• differentiated / undifferentiated
• stimulated / unstimulated

これらの比較が各データセット内で適切に定義できる場合、 そのlog2 ratioや差分を使うことで、 Seriesやプラットフォームを越えた比較が可能になる場合があります。

ただし、この方法がどんな場合でもうまくいくわけではありません。 比に変換しても、Seriesごとに分かれてしまうことは現実にあります。 それでも、各データセット内で共通のコントロールとして使えるサンプルがあるなら、 試してみる価値はあります。

バッチエフェクトへの対処法として、比の基準を用意する

ここまで、複数のGEO Seriesを統合する例を見てきました。 しかし、ここで得られた知見は、複数のデータセットの統合にとどまりません。 より一般的な、バッチエフェクトへの対処法として考えることができます。

バッチエフェクトとは、疾患状態や処理条件といった本来見たい生物学的な違いではなく、 測定時期、シーケンサー、ライブラリ調製キット、試薬ロット、サンプル調製条件、 解析パイプラインなどの違いによって生じる系統的な差です。 この差が大きいと、サンプルは生物学的な状態ではなく、 測定されたSeriesやバッチごとに分かれて見えてしまいます。

このようなバッチエフェクトに対して、補正アルゴリズムを使う方法もあります。 しかし、バッチエフェクトを完全に補正することは、現実には簡単ではありません。 そのため、補正アルゴリズムに頼るより、 データを取得する前、つまり実験計画の段階から対策を組み込む方が安全です。

具体的には、 各バッチに「比の基準」となるサンプルを含めておくことです。 ここでいうコントロールは、必ずしも疾患群に対する健常者、 処理群に対する未処理群という意味ではありません。 異なるバッチや時期のデータをあとから比較するための、 統合解析用の基準サンプルという意味です。

この基準サンプルには、解析で注目したい遺伝子がある程度の強さで発現している必要があります。 基準サンプルでほとんど発現していない遺伝子では、 比を計算したときに値が不安定になり、 わずかな測定誤差でも大きなlog2 ratioとして見えてしまうためです。

同じバッチ内に基準となるサンプルがあれば、 各サンプルをその基準に対する比や差分として表現できます。 絶対的な発現値をバッチ間で直接比較するのではなく、 「同じバッチで測定した基準サンプルに対して、どのように変化しているか」 として比較できるようになります。

前向き研究では、将来の比較可能性を実験計画に組み込む

この考え方は、長期にわたってサンプルを収集する前向き研究で特に役立ちます。 患者登録やサンプル収集が数年にわたる場合、バッチエフェクトは必ず顕在化します。

今年測定したサンプルと、来年、再来年に測定するサンプルは、 そのまま比較できないという前提に立って計画を立てる必要があります。 シーケンサー、ライブラリ調製キット、試薬ロット、read length、read depth、 解析パイプライン、測定施設、担当者など、時間とともに変わる要因が必ず出てくるからです。

まず、各測定バッチに同じ基準サンプルを含めることが基本になります。 そうしておけば、各バッチ内で新しいサンプルを基準サンプルに対する比や差分として表現できます。

長期研究では、基準サンプルそのものの一貫性も問題になるので、さらに対策が必要です。 1種類の基準サンプルだけに依存すると、そのサンプルの保存状態、調製ロット、 細胞状態などが変わった場合に、比の基準自体が不安定になる可能性があります。 そのため、可能であれば、性質の異なる複数の基準サンプルを用意しておく方が安全です。

また、基準サンプルでは、解析で注目したい遺伝子がある程度の強さで発現している必要があります。 研究対象がまだ絞り込まれていない場合は、 できるだけ幅広い遺伝子がある程度の強さで発現している状態を作った人工的なコントロールサンプルが有効です。

長期研究で避けたいのは、 数年かけて集めたサンプルが、測定時期の違いによって比較できなくなることです。 発現値を後から無理にそろえるのではなく、 将来比較できる形でデータを取得しておく。 これが、前向き研究におけるバッチエフェクトへの現実的なリスク管理になります。