きっかけ
Single-Cell2018 大阪セミナー シングルセル2018 株式会社エー・イー企画2018 に参加した際、dry人口の少なさに驚愕した。
てかそもそもscRNA-seqの日本語の記事ほとんどなくない？？
いつも以上に丁寧目を心がけて。Let's Seurat!
Seuratってなあに
single cell RNA seqデータの解析に用いられるRのパッケージです。
QC(Quality Control)、統計解析、可視化全部入り。Seurat独自のオブジェクト（SeuratObject）を作って解析を進めていきます。
以下のような特徴があります。
- クラスタリング、細胞のタイプ、状態の決定に教師データを必要としない。
- 細胞の状態、プロトコル、種が異なるサンプルの統合的な（batch effectを除く）解析が可能。
- 明瞭で、魅力的かつ視覚的に解釈しやすく、dryのラボにとっても、wetのラボにとっても容易に使用可能。
Satija氏の研究室が管理しています。
scRNA解析に使用できるパッケージの中では比較的扱いやすく、擬似系譜解析などの高度な解析を必要としないのであればこのパッケージ１つで事足りてしまいます。
用意するもの
- R (v >= 3.3)
- (Rstudio)
Satija labでは v 3.4.1が推奨されていますが、v 3.4.3, v 3.4.4, v 3.4.5, v 3.5.0でも動作確認済みです。
Rstudioを使うと可視化・フィルタリングを繰り返しながら可視化することができるので便利です。
インストール
CRANで配布されています。
CRANはRのパッケージを配布しているwebサイトです。ミラーサイトが全世界各地に存在します。
RまたはRstudioを開いて以下のコマンドを打ちましょう。
install.packages('Seurat') library(Seurat)
するとバーーっとインストールされます。途中gccのバージョンによっては警告メッセージが出ますが恐れることはありません。
Windows環境では正常にインストールされましたが、CentOS、Ubuntuではライブラリ依存のエラーが出ました。
そういった場合、以下のように不足しているライブラリをインストールする必要があります。
CentOS
yum install [足りないライブラリ名]
Ubuntu
apt-get install [足りないライブラリ名]
実践！！
2018年3月22日にSatija labのホームページで公開された以下のチュートリアルを実践してみましょう。
https://satijalab.org/seurat/immune_alignment.html
使用するデータについて
このチュートリアルでは、以下のデータを使用しています。
サンプル
PBMC（末梢血単核細胞）
詳細
control と stimulatedに分け、stimulatedグループに対しては、インターフェロンβを添加。
インターフェロンβとは
インターフェロンαとβはリンパ球（T細胞、B細胞）、マクロファージ、線維芽細胞、血管内皮細胞、骨芽細胞など多くのタイプの細胞で産生され特に抗ウイルス応答の重要な要素である（詳しくはI型インターフェロンの項を参照）。インターフェロンαとβはマクロファージとNK細胞をともに刺激し、腫瘍細胞に対しても直接的に増殖抑制作用を示す。 インターフェロン - Wikipedia
補足
インターフェロンに対する応答により、cell type特異的な遺伝子発現変化が起きると、異なるstimulation condition同士を統合解析することは困難になります（異なるstimulation condition間で共通のcell typeを同定することが困難になるため）。
そこで、アライメントにButler et al, 2017で用いられているCCA（Canonical Correlation Analysis）と呼ばれる次元圧縮手法を利用し、共通のcell typeの同定（stimulation condition間のbatch effectの除去）を促進します。 異なる手法を使用して取得された4つの膵島のデータセット(ここにある)を統合し、
- データセット中に含まれる細胞のcell typeの同定
- control、stimulatedで保持されるcell type マーカーの取得
- データセットを比較することでcell type特異的な刺激応答を比較
を行います。
解析を行うディレクトリの作成
まずは適当な名前を付けたディレクトリを作成しましょう。
今回はimmune_alignment_matrixというディレクトリを作成し、その中に入ります。
$ mkdir immune_alignment_matrix $ cd immune_alignment_matrix
発現matrixデータのダウンロード
今回のデータセットはここからダウンロードできます。wgetコマンドでダウンロードしましょう。
$ wget https://www.dropbox.com/s/79q6dttg8yl20zg/immune_alignment_expression_matrices.zip
このデータは発現データのmatrix (行が遺伝子、列が細胞) となっており、scRNA-seqでは、基本的にこのフォーマットに変換してから解析を行います。
DropSeqの論文中では、このフォーマットをdge (digital expression matrix)ファイルとか呼んでたりします。 この記事でもdgeファイルと呼ぶことにします。
ダウンロードしたデータを解凍しましょう。解凍後のzipファイルは削除しましょう。
$ unzip immune_alignment_expression_matrices.zip
$ rm immune_alignment_expression_matrices.zip
解凍後、出力されたデータを見てみましょう。
$ ls
immune_alignment_expression_matrices.zip* immune_control_expression_matrix.txt.gz immune_stimulated_expression_matrix.txt.gz
controlのdgeファイルがimmune_control_expression_matrix.txt.gz
stimulatedのdgeファイルがimmune_stimulated_expression_matrix.txt.gzです。
さらに解凍します。
$ gunzip immune_stimulated_expression_matrix.txt.gz
$ gunzip immune_control_expression_matrix.txt.gz
$ ls
immune_alignment_expression_matrices.zip* immune_control_expression_matrix.txt immune_stimulated_expression_matrix.txt
Seuratオブジェクトの作成
データを読み込み、ctrl.data、stim.dataというオブジェクトにそれぞれ格納します。
オブジェクトを作成する際に、発現している総細胞数、発現している総遺伝子数等でデータをフィルタリングすることができます。
> ctrl.data <- read.table("immune_control_expression_matrix.txt", sep = '\t', row.names = 1) > stim.data <- read.table("immune_stimulated_expression_matrix.txt", sep = "\t", row.names = 1)
seuratオブジェクトの作成は、CreateSeuratObject関数で行います。
> ctrl <- CreateSeuratObject(raw.data = ctrl.data, project = "IMMUNE_CTRL", min.cells = 5)
| options | description |
|---|---|
| project | プロジェクト名を指定する。ここでえ与えた情報は、メタデータとしてオブジェクト内に保持される。 |
| min.cells | 総細胞数の下限を指定する。５を指定することで、５細胞以上の細胞のみにフィルタリングする。 |
| min.genes | 総遺伝子数の下限を指定する。使用していない。 |
metadataの付与
メタデータの情報は、作成したseuratオブジェクト中metadataスロットのctrl列にも保持させておきます。
> ctrl@meta.data$stim <- "CTRL" > head(ctrl@meta.data$stim) [1] "CTRL" "CTRL" "CTRL" "CTRL" "CTRL" "CTRL"
上では、head()関数でmetadataスロットのctrlの先頭６要素を抽出して確認しています。
class()関数を使えば、オブジェクトの型を調べることができます。
> class(ctrl) [1] "seurat" attr(,"package") [1] "Seurat"
確かにseuratオブジェクトが作成されていることがわかります。
フィルタリング
FilterCells()関数を使えば、フィルタリング条件、上限、下限を指定して細胞をフィルタリングすることができます。
今回は、発現している総遺伝子数で細胞をフィルタリングします（500 =< nGene =< ∞）。
> ctrl <- FilterCells(ctrl, subset.names = "nGene", low.thresholds = 500, high.thresholds = Inf)
| options | description |
|---|---|
| subset.names | フィルタリング条件を指定する。"nGene"を指定することで、総遺伝子数でのフィルタリングを行う。 |
| high.thresholds | フィルタリングの上限を指定する。 |
| high.thresholds | フィルタリングの下限を指定する。 |
正規化
次に発現量の正規化を行います。
デフォルトでは、各細胞について総発現量で割った値にSize Factorを掛け合わせ、対数値変換を行っています。
Size Factor（scale.factor）はデフォルトで10,000 となっています。
> ctrl <- NormalizeData(ctrl, scale.factor = 10000) Performing log-normalization 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| **************************************************|
正規化後のデータはhvg.info スロットに格納されます。
head(ctrl@hvg.info) gene.mean gene.dispersion HBB 2.7100010 8.221465 HBA2 1.6422454 7.115376 HBG1 0.2289593 6.932304 HBA1 1.1739589 6.571082 FTH1 6.7435404 6.478593 FTL 5.8811076 6.320713 gene.dispersion.scaled HBB 4.343644 HBA2 6.552497 HBG1 17.496792 HBA1 10.889386 FTH1 0.000000 FTL 0.000000
Z-Scoreの算出
Seuratでは、発現量の一般的な統計指標としてZ-Scoreを使用します。
Z-Scoreとは
分布の平均値からのずれを示す値．注目している標本値と分布の平均値の差を分布の標準偏差で割った値で定義される．Z-scoreの絶対値が大きければ大きい程，分布の平均値からのずれが大きいことを示している． Z-score：バイオキーワード集｜実験医学online：羊土社
Z-Scoreの算出には、ScaleData関数を使用します。
batch effectやconfounding factorと呼ばれる、technicalな影響に由来するデータ変動を取り除く必要がある場合には、
vars.to.regress引数に除きたい変数を与えてあげることでtechnicalな影響を取り除くこともできます（percent.mitoなど）。
今回は特に引数は指定せず、単にZ-Scoreを求める処理を行いましょう。
スケーリング後のデータは、scale.dataスロットに格納されます。
ScaleDataを行う前はこのスロットは空になっています。
> ctrl@scale.data NULL
ScaleData()関数でデータのスケーリングを行うと、、
> ctrl <- ScaleData(ctrl) Scaling data matrix |==============================================================================| 100%
> ctrl@scale.data[0:6,0:6] AAACATACATTTCC.1 AAACATACCAGAAA.1 AAACATACCTCGCT.1 AAACATACCTGGTA.1 AL627309.1 -0.02458278 -0.02458278 -0.02458278 -0.02458278 RP11-206L10.2 -0.02445482 -0.02445482 -0.02445482 -0.02445482 LINC00115 -0.07280831 -0.07280831 -0.07280831 -0.07280831 NOC2L -0.28284480 -0.28284480 -0.28284480 -0.28284480 KLHL17 -0.04240574 -0.04240574 -0.04240574 -0.04240574 PLEKHN1 -0.05282801 -0.05282801 -0.05282801 -0.05282801 AAACATACGATGAA.1 AAACATACGGCATT.1 AL627309.1 -0.02458278 -0.02458278 RP11-206L10.2 -0.02445482 -0.02445482 LINC00115 -0.07280831 -0.07280831 NOC2L -0.28284480 -0.28284480 KLHL17 -0.04240574 -0.04240574 PLEKHN1 -0.05282801 -0.05282801
このように全発現量のZ-Scoreがctrl@scale.data に格納されます。
stimulatedに関しても
ここまでの処理を、stimulatedのデータに関しても行います。
# Setup stimulated object stim <- CreateSeuratObject(raw.data = stim.data, project = "IMMUNE_STIM", min.cells = 5) stim@meta.data$stim <- "STIM" # Filtering stim <- FilterCells(stim, subset.names = "nGene", low.thresholds = 500, high.thresholds = Inf) # Normalization stim <- NormalizeData(stim) # Convert expression value to Z-Score stim <- ScaleData(stim)
使用する遺伝子の選定
FindVariableGenes関数で、分散の大きい（発現変動の大きい）遺伝子を抽出することができます。
以下では、g1にcontrol内で発現変動の大きい遺伝子1000遺伝子を、
g2にstimulated内で発現変動の大きい遺伝子1000遺伝子を、それぞれ抽出しています。
> ctrl <- FindVariableGenes(ctrl) Calculating gene means 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| **************************************************| Calculating gene variance to mean ratios 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| **************************************************|
出力
> stim <- FindVariableGenes(stim) Calculating gene means 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| **************************************************| Calculating gene variance to mean ratios 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| **************************************************|
出力
g.1 <- head(rownames(ctrl@hvg.info), 1000) g.2 <- head(rownames(stim@hvg.info), 1000)
g1、g2で共通して発現変動の大きい遺伝子をunique関数で抽出します。
その後、intersecit関数で、各scale.dataスロット中に含まれるもののみにフィルタリングします。
genes.use <- unique(c(g.1, g.2)) genes.use <- intersect(genes.use, rownames(ctrl@scale.data)) genes.use <- intersect(genes.use, rownames(stim@scale.data))
これで前処理は終了です。次にこれらの遺伝子の発現データを可視化します。
batch effect除去の実行
始めのほうに述べた通り、異なるstimulation condition同士を統合解析するために、CCA（Canonical Correlation Analysis）をベースにしたbatch effect除去を行います。
batch effectを除去する際には、以下の手順が取られます。
- CCAを行い、共通のcell typeに由来する細胞を同定する。
- それらの発現差異を小さくするような処理を行う。
CCAによるオブジェクトのマージ
まず、CCAを行います。
CCAとは
正準相関分析（Canonical Correlation Analysis：CCorA，ときどきCCAと略されることがあるが，我々は CCA を 正準コレスポンデンス分析（Canonical Correspondence Analysis）の略として使用する）は，2つの変数の集合の間の関係性を調査できるたくさんの統計手法の中の1つである．これは，変数の2つの集合間の相関を調査し，両方の表にできるだけ相関していて，お互いが直交している正準変数の集合をこれらの表から抽出する． (https://www.xlstat.com/ja/solutions/features/canonical-correlation-analysis-ccora)
要するに相関要素を組み込んだPCAといった具合です。
SeuratでCCAを行うには、RunCCA関数を用います。
> immune.combined <- RunCCA(ctrl, stim, genes.use = genes.use, num.cc = 30) Running CCA Merging objects Performing log-normalization 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| **************************************************| Scaling data matrix |==============================================================================| 100%
| options | description |
|---|---|
| genes.use | 使用する遺伝子の選択 |
| num.cc | canonical vectorの数 |
CCAの結果は、drスロット中の$cca に格納されます。
この関数を使うと、2つあったseuratオブジェクト（ctrl、stim）は1つのseuratオブジェクト（immune.co
mbined）へとマージされます。
どちらのサンプルに由来するか（stimulated か、controlか）という情報はmeta.dataスロットに保存されます。
データの可視化
CC1, CC2を軸にしたDimention plot
次元削減後の変数をオブジェクト内に格納した後は、それらの値を用いてbiplotを行うことができます。
> p1 <- DimPlot(object = immune.combined, reduction.use = "cca", group.by = "stim", pt.size = 0.5, do.return = TRUE) > p1
| options | description |
|---|---|
| features.plot | 次元削減に用いる手法を指定 |
| group.by | 色分けするgroup名を指定 |
| do.return | 変数に格納する場合はTrue |
出力
CC1をy軸に取ったViolin plot
1つのCCについて、Violin plotを行います。
p2 <- VlnPlot(object = immune.combined, features.plot = "CC1", group.by = "stim", do.return = TRUE)
| options | description |
|---|---|
| features.plot | y軸に取るCCを指定する。"CC1"を指定している。 |
| group.by | 色分けするgroup名を指定する。セルタイプの情報は"stim"という変数に格納しているため、これを指定する。 |
| do.return | 変数に格納する場合はTrueにする。 |
出力
上記二つを左右に並べてプロット
plot_grid関数を用いてプロットを左右に並べます。
plot_grid(p1, p2)
出力
各CCで有意に発現している遺伝子を検出
PrintDim 関数を使用します。 発現変動の有意性については、分散/平均の大きさから算出しています。
> PrintDim(object = immune.combined, reduction.type = "cca", dims.print = 1:2, + genes.print = 10) [1] "CC1" [1] "RPS6" "RPS18" "RPL21" "RPL13" [5] "RPL13A" "RPS14" "RPL3" "RPS3" [9] "RPS2" "RPL32" [1] "" [1] "TYROBP" "FCER1G" "C15orf48" [4] "SOD2" "FTL" "ANXA5" [7] "CST3" "CD63" "TIMP1" [10] "TYMP" [1] "" [1] "" [1] "CC2" [1] "NKG7" "GNLY" "GZMB" "CCL5" [5] "PRF1" "CST7" "KLRD1" "GZMH" [9] "CLIC3" "GZMA" [1] "" [1] "RPS9" "RPL8" "HLA-DRA" [4] "RPS23" "RPL11" "RPS8" [7] "RPL18A" "RPL32" "RPL12" [10] "RPL26" [1] "" [1] ""
| options | description |
|---|---|
| reduction.type | 次元削減の手法を選択する。 |
| dim.print | プロットする次元数を指定する。 |
| genes.print | 各CCについて、出力する遺伝子数を選択す。 |
有意に発現している遺伝子を見てみると、 CC1では骨髄球特異的な遺伝子が有意に発現しており、 CC2ではリンパ球特異的な遺伝子が有意に発現していることがわかります。
使用するCCの選定
下流解析で使用するCCの数を決定するため、
CC内での相関の強さを縦軸、各CCを横軸取ったプロットを行います。
> p3 <- MetageneBicorPlot(immune.combined, grouping.var = "stim", dims.eval = 1:30, + display.progress = FALSE) Rescaling group 1 Rescaling group 2
出力
図を見てみると、CC20より先においてCC内の相関の強さが下落していることがわかるため、
CC1-20を下流の解析に使用することにします。
各CCにおけるヒートマップの描画
下流解析で使用するCC数の決定には、各CCにおけるヒートマップも有用です。
DimHeatmap関数で描画することができます。
DimHeatmap(object = immune.combined, reduction.type = "cca", cells.use = 500, dim.use = 1:9, do.balanced = TRUE)
| options | description |
|---|---|
| reduction.type | 次元削減の手法を選択 |
| cells.use | 使用する細胞数を選択 |
| dim.use | プロットする次元数 |
| do.balanced | 正に有意に発現している遺伝子と、負に有意に発現している遺伝子を同じ数プロットする |
出力
