Seuratを駆使する会 ①

きっかけ

いつも以上に丁寧目を心がけて。Let's Seurat!

Seuratってなあに

single cell RNA seqデータの解析に用いられるRのパッケージです。

QC(Quality Control)、統計解析、可視化全部入り。Seurat独自のオブジェクト(SeuratObject)を作って解析を進めていきます。

以下のような特徴があります。

  • クラスタリング、細胞のタイプ、状態の決定に教師データを必要としない。
  • 細胞の状態、プロトコル、種が異なるサンプルの統合的な(batch effectを除く)解析が可能。
  • 明瞭で、魅力的かつ視覚的に解釈しやすく、dryのラボにとっても、wetのラボにとっても容易に使用可能。

Satija氏の研究室が管理しています。
scRNA解析に使用できるパッケージの中では比較的扱いやすく、擬似系譜解析などの高度な解析を必要としないのであればこのパッケージ1つで事足りてしまいます。

用意するもの

f:id:kimoppy126:20180814060307p:plain

  • 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

使用するデータについて

このチュートリアルでは、以下のデータを使用しています。

出典
Multiplexed droplet single-cell RNA-sequencing using natural genetic variation | Nature Biotechnology

サンプル
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ファイルと呼ぶことにします。

f:id:kimoppy126:20180609102950p:plain:w300

ダウンロードしたデータを解凍しましょう。解凍後の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.datastim.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 総細胞数の下限を指定する。5を指定することで、5細胞以上の細胞のみにフィルタリングする。
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の先頭6要素を抽出して確認しています。

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%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

出力 f:id:kimoppy126:20180813090748p:plain

> 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%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

出力 f:id:kimoppy126:20180813090944p:plain

g.1 <- head(rownames(ctrl@hvg.info), 1000)
g.2 <- head(rownames(stim@hvg.info), 1000)

g1g2で共通して発現変動の大きい遺伝子を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を除去する際には、以下の手順が取られます。

  1. CCAを行い、共通のcell typeに由来する細胞を同定する。
  2. それらの発現差異を小さくするような処理を行う。

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オブジェクト(ctrlstim)は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

出力 f:id:kimoppy126:20180813212224p:plain

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にする。

出力 f:id:kimoppy126:20180813212323p:plain

上記二つを左右に並べてプロット

plot_grid関数を用いてプロットを左右に並べます。

plot_grid(p1, p2)

出力 f:id:kimoppy126:20180813212409p:plain

各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

出力 f:id:kimoppy126:20180814070022p:plain

図を見てみると、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 正に有意に発現している遺伝子と、負に有意に発現している遺伝子を同じ数プロットする

出力 f:id:kimoppy126:20180813212526p:plain

続きは

いつか書くよ。