RNA-Seqによる発現データの解析には、リードカウントではなくTPM、FPKM、またはRPKMを用いるべきだと信じている人が多いです。その理由の一つに、サンプル間の系統誤差を回避することがあげられます。しかし、100万リードあたりで均すやり方は、一種のグローバルノーマライズに過ぎず、これで除去できるのは線形の系統誤差だけです。これは、現実のオミクスデータ解析に対して甘すぎる前提と言わざるを得ません。特に大規模な患者集団のデータなどでは必ずデータ品質の問題が紛れ込むので注意が必要です。
それでは、GSE159751を例としてみてみましょう。このデータセットは、2群に分かれる26サンプルで構成されています。Supplementary fileよりGSE159751_RAW.tar をダウンロードできますが、このファイルはGSMのレコードによるとFPKM値のテーブルのようです。Fig1をの右下のヒストグラムを見てください。各サンプルのFPKM値の分布が視覚化されていますが、大まかに言うと、一山型と二山形の二つのタイプがあるように見えます。これは、系統誤差の存在を強く示唆しています。
ということは、このデータセットでクラスタリングを実行すると、生物学的文脈ではなく、非線形バイアスに従ったクラスターが形成されることが予想されます。Fig2 のPCAの結果をご覧ください。予想どおり、青と黄色のDisease Stateは第一主成分（横軸）とあまり関係ないように見えます。（Fig2 上図）そこで、PC1スコアの高い群と低い群にピンクと緑の印をつけてみました。（Fig2 下図）
そして、階層型クラスタリングを実行してみました。もちろんサンプルは、Disease Stateにより分かれるのではなく、予想通りPC1スコアによってクラスターができています。さらに興味深いことに、PCAスコアによる分類は、ヒストグラムの形と関係がありそうです。（Fig3）
正規化の手法を学んだ方であれば、もっと洗練された正規化法を適用すればいいじゃないかと思うかもしれません。そこで、このデータをquantile normalizationによる正規化をして、再度クラスタリングを実行してみました。しかし、結果はそう変わりません。PC1スコアの高いサンプルと低いサンプルが再びクラスターを形成しています。（Fig4）
実は、「たとえ分布の形を似せるような強力な正規化手法を用いたとしても、非線形の系統誤差がある場合、それを補正することはできない」ということは、大量のマイクロアレイデータの解析から既に分かっていることですが、これはRNA-Seqになっても変わらないと言えます。オミクスデータの非線形系統誤差を補正するバイオインフォマティクス的手段は存在しないのです。悪いことに、そのような「洗練された」正規化手法は、実験者の目をデータの質の問題から逸らす役割を果たしてきました。（詳細は、RMAの問題や、Z-score normalizationに関するトピックをご覧ください。）
もう一つの実験をお見せしましょう。FASTQファイルをダウンロードして、まったく異なる数値化パイプライン（fastp, HISAT2, StringTie）を通して算出された発現量のデータを解析してみました。異なるアルゴリズムで数値化されていますから、当然ながらヒストグラムの形はまったく異なって見えます。しかしそれでも、PC1スコアの高いサンプルと低いサンプルは、またしてもクラスターを形成しています。（Fig5）つまり、この差はFASTQファイルに本質的な違いとして存在しているのです。そしてその差は、その後のデータ解析における手法の選択や工夫を凌駕するほど、圧倒的に強いのです。
データ解析は魔法ではありません。比較不能なデータセットの前では無力なのです。従って、第一のルールは、均質なサンプルを得ることです。しかし、それがほとんど不可能なケースもあります。特に、ある程度以上のサンプル数を集める実験や、ことなる時点でとられるサンプルを集める実験の時です。ですから、私たちは第二のルールをアドバイスています。それは、非線形系統誤差が含まれることは不可避という前提で、それに対して頑健性を持たせるよう実験デザインを工夫することです。
そして、上記と同じ理由で、様々な研究者が出したRNA-Seqのデータを統合して解析することは、特殊なケースを除いてほとんど不可能です。もしそのようなデータ解析をしようと思っているのでしたら、どうぞSubioにご相談ください。