5. シングルセルデータの正規化
前回はこちら。
解析の流れ
クオリティをチェックして、不要なデータをフィルタリングで除外したら、解析の準備は終了です。この後の解析の流れは、基本的に次のようなステップになります。
(1) データの正規化
(2) 複数サンプルであれば、データのマージ（インテグレーション）
(3) クラスタリング (tSNE or UMAP)
(4) 各クラスターのマーカー遺伝子を算出
(5) マーカーからクラスターの細胞種を推定
(6) 特定のクラスターと別のクラスターを比較して発現変動遺伝子を算出
(7) 発現変動遺伝子の機能解析（GO解析、エンリッチメント解析(GSEA)）
UMAP でクラスタリングしたとしても、クラスターの細胞の種類を同定することはできません。細胞がいくつぐらいのグループに分かれるかが、確認できるだけです。
そのため、各クラスターにおいて、どのような遺伝子が（ほかのクラスターより）多いのかという特徴（マーカー遺伝子）を見て、細胞の種類を推測する必要があります。
近年では、データベース化された過去の scRNA-seq のデータから、細胞の種類をアノテーションすることができるようになってきましたが、あくまで推測であることに注意する必要があります。詳細は後述するとして、まずは、「正規化」について解説します。
5.1 シングルセルデータの正規化
データ解析において、何かを比較して違いを見たいのであれば、（その違いが見えやすいように）データを補正する作業が必要になります。それが「正規化 (normalization)」です。
「正しい」という文字の響きから、正規化すれば、データが綺麗になって何でも正しい解析ができるような印象を与えるかもしれませんが、そうではありません。ばらつきを補正したりするのも正規化に含められることもありますが、それも「違い」を見やすくするための工夫の1つです。そういった補正によって、見えやすくなる違いもあれば、見えにくいままの違いもあります。
通常、正規化は、「何を見たいか」目的を持って行うものです。シングルセルデータの場合は、当然、細胞と細胞の違い（遺伝子ごとに）を見たいわけなので、そのための正規化を行います。
4. クオリティチェックとフィルタリング のバイオリンプロット(nCount_RNA)でも解説したように、シングルセルのデータでは、数百本のリードしか検出されていない細胞もあれば、数万本のリードが検出されている細胞もあります。
全体で1,000リード検出された細胞において、ある遺伝子が10本検出されていたとします。一方、全体で10,000リード検出された細胞において、その遺伝子が10本検出されていたらどうでしょう。同じ10本だからといって、この遺伝子については違いがないと考えて良いでしょうか？
このように、そのままのリードカウントの結果（raw カウント）で比較をすることは、遺伝子レベルで、細胞の違いを見る場合に好ましくありません。何らかの補正（正規化）が必要です。
5.2 Seurat で用いられる正規化
遺伝子発現データを正規化する計算方法（アルゴリズム）については、通常の RNA-seq の頃（その前のマイクロアレイの頃）から、さまざまな手法が提案されていますが、万能な方法はありません。RNA-seq については、 TMM という手法がよく使われますが、これは2グループ間の発現変動遺伝子を検出するための正規化方法で、シングルセルデータは想定されていません。
Seurat で使用される正規化方法（アルゴリズム）は、おもに2つあります。
- LogNormalize
- sctransform
LogNormalize
LogNormalize は、プログラム中で使用される正規化のための関数 NormalizeData() のデフォルトに設定されているアルゴリズムです。何も指定しない場合は、この LogNormalize が適用されます。
計算方法としては単純な手法です。細胞ごとに、個々の遺伝子のリードカウントをリード数の合計で割って、scale.factor (=10,000) をかけたのち、1を加えてから log 変換する（底は e ）というものです。
LogNormalize: Feature counts for each cell are divided by the total counts for that cell and multiplied by the ‘scale.factor’. This is then natural-log transformed using ‘log1p’
計算方法から分かるように、遺伝子ごとのカウントを1万リードあたりに換算したものです。log変換するのは、値の範囲を狭くするためです。また、その際に1を加えるのは、ゼロのカウントをlog変換後もゼロにするためで、統計的にはよく行われる処置です。
この正規化により、全体で1,000リード検出された細胞と、全体で10,000リード検出された細胞の間で、同じ遺伝子を比較できるようになります。
＊処理内容を確認すると、Seurat で正規化したシングルセルデータは、RNA-seq とは異なり、CPM や TPM 単位ではないことが分かります。また、トランスクリプトの長さを考慮した正規化は行われませんので、同じ細胞内の異なる遺伝子AとBを比較することは想定されていません。使用するプログラムによって、正規化方法は異なるので注意が必要です。
sctransform
LogNormalize の後に、Seurat の開発者らによって sctransform が提案されました。正規化のあとの解析において、変数選択や、次元削減のステップを考慮しているとされます。
しかし、最近 (Seurat v5 ) では、インテグレーションの際に、再び logNormalize を用いるようになっています。LogNormalize を完全に置き換えるものではないようです。インテグレーションは、刺激前と刺激後など、複数のサンプルがある場合の処理ですが、詳細は後述します。
＊選択肢として、sctransform を使ってインテグレーションする方法も選べます。
NormalizeData() の使用例
説明が長くなりましたが、Seurat 上で正規化を行う方法は簡単です。NormalizeData() 関数を使います。正規化に用いるアルゴリズムを引数に指定しない場合、上記の LogNormalize が適用されます。
pbmc <- NormalizeData(pbmc)
正規化した pbmc オブジェクトで、元の pbmc を上書きするような書き方です。しかしながら、実際には、元のリードカウント (raw カウント) が上書きされるわけではありません。
Seurat オブジェクトには、データを格納する場所（スロット）が複数あります。raw のカウントは、RNA@counts スロットに入っています。
> pbmc$RNA@counts[c("CD3D", "TCL1A", "MS4A1"), 1:3]
3 x 3 sparse Matrix of class "dgCMatrix"
AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1
CD3D 4 . 10
TCL1A . . .
MS4A1 . 6 .
一方、正規化したカウントデータは、 RNA@data スロットに入っています。
> pbmc$RNA@data[c("CD3D", "TCL1A", "MS4A1"), 1:3]
3 x 3 sparse Matrix of class "dgCMatrix"
AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1
CD3D 2.864242 . 3.489706
TCL1A . . .
MS4A1 . 2.583047 .
raw のカウントは、検出されたリードの本数そのままなので、整数値になっています。それに対して、正規化されたデータは、実数値であることが確認できます。
データの正規化ができたら、次のステップに進みます。
Discussion