Seurat を駆使する会② - ばいばいバイオ

bioinfo.hatenablog.com

続き！

前回までのあらすじ
batch effectの除去
batch effect除去結果の確認
マーカー遺伝子の検出
マーカー遺伝子発現量の描画
新しいラベルの付与
マーカー遺伝子発現量の描画②
マーカー遺伝子発現量の描画③
オブジェクトの保存

前回までのあらすじ

インターフェロン刺激応答前後のPBMCデータセットを使用して、

データの入力
フィルタリング
正規化
Z-Scoreの算出
batch effectを除くためのCCA

を行ったのでした。

それでは今回は、実際にbatch effectを除くところから入っていきましょう。

CCAをなんのために実行したのかというと！
共通のcell typeをもつ細胞群を同定するためでした。

それではこの後何をしていくのか！
これら共通のcel typeをもつと推測される細胞同士の発現差異を小さくするような処理（batch effectの除去）を行います。

immune.combinedオブジェクトがある前提で話を進めます。

batch effectの除去

f:id:kimoppy126:20181024230059p:plain
Butler et al (2018) Integrated analysis of single cell transcriptomic data across conditions, technologies, and species, Nature Biotechnology volume 36, pages 411–420

AlignSubspace関数では、次元圧縮した結果を使って
DTW (Dynamic Time Warping)により、細胞発現パターンの類似度が最も近似する細胞の組み合わせを同定します。

具体的には、以下の手順を取ります。

各データセットについて各成分ベクトル（CC）ごとに遺伝子発現の相関を見て、どちらのデータセットにおいても発言変動に強く寄与している遺伝子を同定する（ Biweight midcorrelationを使用）。
各データセットについて、Biweight midcorrelationの高い順に、指定したnum.genes分の'metagene'を構成する。
95%基準範囲が一致するように、各'metagene'を線形変換し、2つの'metagene'の差異が10-90四分位範囲に収まるようにする（globalなアライメント）。
dynamic time warping (DTW) を使用して、大きいほうのデータセットに、小さいほうのデータセットの各細胞をマップする。
1 ~ 4を各CCについて行う。

うーーんちょっとよくわからない。
詳細は下記ドキュメントを参考にしてください。。 AlignSubspace function | R Documentation

とりあえずやっていることは、

CCごとに有意に発現している遺伝子の同定。
それら有意な細胞同士の発現パターンを揃える。

といったことかと思います。

immune.combined <- AlignSubspace(immune.combined, reduction.type = "cca", grouping.var = "stim", dims.align = 1:20)

options	description
reduction.type	次元削減の手法を選択
grouping.var	色分けに使用する変数
dims.align	使用するCCの次元数

batch effect除去結果の確認

それではCCAを用いてbatch effectを除去した結果をtSNE plotにより可視化してみましょう。
batch effect除去前同様RunTSNE関数を使用するのですが、reduction.use = "cca.aligned"を指定することにより、上記アライメント結果を使用していることに注意してください。

また、このtSNEによる次元圧縮結果においてクラスタリングを行います。
クラスタリングとは呼んで時の如く、似たような発現パターンを持った細胞を探してクラスター（集団）に分けることを指します。
Seuratではshared nearest neighbor (SNN) という手法をベースにしたクラスタリング手法を使用しています。
ちょっとここで詳細を描きたくはないので別記事にでもします。
今すぐ知りたい方は下記ドキュメントからどうぞ。
FindClusters function | R Documentation

immune.combined <- RunTSNE(immune.combined, reduction.use = "cca.aligned", dims.use = 1:20, 
    do.fast = T)

options	description
reduction.use	tSNEによる次元圧縮に使用する値
dims.use	次元圧縮に使用する次元数
do.fast	Barnes-hutのアルゴリズムを使用するか否か。これを使用すると実行速度が早くなる

immune.combined <- FindClusters(immune.combined, reduction.type = "cca.aligned", 
    resolution = 0.6, dims.use = 1:20)

options	description
reduction.type	クラスタリングに使用する値
resolution	クラスターをどの程度の解像度で作成するか
dims.use	使用する次元数

tSNE plotにより可視化します。

p1 <- TSNEPlot(immune.combined, do.return = T, pt.size = 0.5, group.by = "stim")
p2 <- TSNEPlot(immune.combined, do.return = T, pt.size = 0.5, do.label = T)
plot_grid(p1, p2)

options	description
do.return	ggplot2のオブジェクトを返すか否か
pt.size	プロットに使用するドットのサイズを調整する値
grouping.var	グループ分けに使用する変数
do.label	クラスター番号のラベルを表示するか否か

f:id:kimoppy126:20181102234911p:plain

マーカー遺伝子の検出

セルタイプ特異的に発現変動を示すような遺伝子のことをマーカー遺伝子と呼んだりします。
インターフェロン刺激前後でクラスター特異的に発現しているような遺伝子をマーカー遺伝子として検出してみましょう。
この処理を行うには、FindConservedMarkers関数を使用します。

因みに、ここでは2状態での有意性（p-value）の統合解析を行っているのですが、
Seuratでは内部的にMetaDEというパッケージで実装されているmeta-analysis methodを使用しています。
例えばクラスター7におけるマーカー遺伝子を検出するには、以下のように行います。

> cluster1.markers <- FindConservedMarkers(immune.combined, ident.1 = 1, grouping.var = "stim", 
    print.bar = FALSE)
Testing 1_CTRL vs 0_CTRL, 2_CTRL, 3_CTRL, 4_CTRL, 5_CTRL, 6_CTRL, 7_CTRL, 8_CTRL, 9_CTRL, 10_CTRL, 11_CTRL, 12_CTRL
Testing 1_STIM vs 0_STIM, 2_STIM, 3_STIM, 4_STIM, 5_STIM, 6_STIM, 7_STIM, 8_STIM, 9_STIM, 10_STIM, 11_STIM, 12_STIM

options	description
ident.1	マーカー遺伝子を検出する対象のクラスター番号
grouping.var	グループ分けに使用する変数

結果の上位6遺伝子を確認してみましょう。

> head(cluster1.markers)
       CTRL_p_val CTRL_avg_logFC CTRL_pct.1 CTRL_pct.2 CTRL_p_val_adj
SELL            0      1.4666075      0.649      0.127              0
GIMAP7          0      1.3515324      0.853      0.274              0
LTB             0      1.2678666      0.806      0.256              0
RPS6            0      0.8749207      1.000      0.992              0
RPS18           0      0.8567551      1.000      0.987              0
RPS14           0      0.8561853      1.000      0.985              0
          STIM_p_val STIM_avg_logFC STIM_pct.1 STIM_pct.2 STIM_p_val_adj
SELL    0.000000e+00      1.3684940      0.861      0.344   0.000000e+00
GIMAP7  0.000000e+00      1.2279994      0.961      0.443   0.000000e+00
LTB    1.881262e-307      1.1502014      0.732      0.253  2.643737e-303
RPS6    0.000000e+00      0.9697095      0.999      0.938   0.000000e+00
RPS18   0.000000e+00      0.8860020      1.000      0.969   0.000000e+00
RPS14   0.000000e+00      0.9231849      1.000      0.963   0.000000e+00
            max_pval minimump_p_val
SELL    0.000000e+00              0
GIMAP7  0.000000e+00              0
LTB    1.881262e-307              0
RPS6    0.000000e+00              0
RPS18   0.000000e+00              0
RPS14   0.000000e+00              0

クラスター1では、SELLというCD4ナイーブT細胞のマーカー遺伝子が有意に発現していることがわかります。

マーカー遺伝子発現量の描画

同様に他のクラスターにおいても有意に発現しているマーカー遺伝子を検出した結果を参考に、FeaturePlot関数により、
各遺伝子の発現量を載せたヒートマップを描画してみましょう。

FeaturePlot(object = immune.combined, features.plot = c("CD3D", "SELL", "CREM", 
    "CD8A", "GNLY", "CD79A", "FCGR3A", "CCL2", "PPBP"), min.cutoff = "q9", cols.use = c("lightgrey", 
    "blue"), pt.size = 0.5)

f:id:kimoppy126:20181103001103p:plain

options	description
features.plot	プロットする遺伝子
min.cutoff	プロットする遺伝子発現量のカットオフ値を分位数で指定
cols.use	使用する色の組み合わせ
pt.size	プロットに使用するドットのサイズを調整する値

新しいラベルの付与

先ほどのヒートマップ結果から、各マーカー遺伝子に対応するセルタイプ名でラベルを上書きします。

new.ident <- c("CD14 Mono", "CD4 Naive T", "CD4 Memory T", "B", "CD16 Mono", 
    "T activated", "CD8 T", "NK", "DC", "B activated", "Mk", "pDC", "Eryth")
for (i in 0:12) {
    immune.combined <- RenameIdent(object = immune.combined, old.ident.name = i, 
        new.ident.name = new.ident[i + 1])
}

TSNEPlot(immune.combined, do.label = T, pt.size = 0.5)

options	description
old.ident.name	rename前のラベル
new.ident.name	rename後のラベル

f:id:kimoppy126:20181103001502p:plain

マーカー遺伝子発現量の描画②

インターフェロン刺激前後において保存されたセルタイプのマーカー遺伝子を検出する際にはSplitDotPlotGG関数を使用します。
発現細胞の割合と発現量の2つの側面から遺伝子発現を見ることができます。
以下ではセルタイプごとに2-3のマーカー遺伝子をプロットしています。

immune.combined@ident <- factor(immune.combined@ident, levels = (c("pDC", "Eryth", 
    "Mk", "DC", "CD14 Mono", "CD16 Mono", "B activated", "B", "CD8 T", "NK", 
    "T activated", "CD4 Naive T", "CD4 Memory T")))
markers.to.plot <- c("CD3D", "CREM", "HSPH1", "SELL", "GIMAP5", "CACYBP", "GNLY", 
    "NKG7", "CCL5", "CD8A", "MS4A1", "CD79A", "MIR155HG", "NME1", "FCGR3A", 
    "VMO1", "CCL2", "S100A9", "HLA-DQA1", "GPR183", "PPBP", "GNG11", "HBA2", 
    "HBB", "TSPAN13", "IL3RA", "IGJ")
sdp <- SplitDotPlotGG(immune.combined, genes.plot = rev(markers.to.plot), cols.use = c("blue", 
    "red"), x.lab.rot = T, plot.legend = T, dot.scale = 8, do.return = T, grouping.var = "stim")

options	description
genes.plot	プロットする遺伝子
cols.use	使用する色の組み合わせ
x.lab.rot	x軸のラベルを縦にするか否か
plot.legend	レジェンドを表示するか否か
dot.scale	プロットに使用するドットのサイズを調整する値
do.return	ggplot2のオブジェクトを返すか否か
grouping.var	グループ分けに使用する変数

f:id:kimoppy126:20181103001610p:plain

マーカー遺伝子発現量の描画③

マーカー遺伝子発現量を見る別の方法として、FeatureHeatmap関数を使用するとクラスターごと、状態ごとの遺伝子発現のヒートマップを描くことができます。

FeatureHeatmap(immune.combined, features.plot = c("CD3D", "GNLY", "IFI6", "ISG15", 
    "CD14", "CXCL10"), group.by = "stim", pt.size = 0.25, key.position = "top", 
    max.exp = 3)

options	description
features.plot	プロットする遺伝子
group.by	グループ分けに使用する変数
pt.size	プロットに使用するドットのサイズを調整する値
key.position	レジェンドを表示する場所
max.exp	スケーリングした発現量の最大値

f:id:kimoppy126:20181103002534p:plain CD3DやGNLY 等といったT 細胞やNK細胞/CD8 T細胞の標準的な細胞マーカーは、インターフェロンの刺激応答によってほとんど発現分布が変化していないことがわかります。
一方、 IFI6やISG15といったインターフェロン応答遺伝子においては、すべてのセルタイプにおいて発現量が増加していることが示唆されます。

また、CD14、CXCL10は、細胞特異的な発現変動をする遺伝子となっています。
CD14はCD14単球のクラスターにおいて発現変動が低下していることがわかり、
CXCL10はインターフェロン刺激後の単球とB細胞のクラスターにおいて、顕著に遺伝子発現量が増幅していることがわかります。

オブジェクトの保存

saveRDS関数を使えば、Rで使用できるオブジェクトを、バイナリ形式で保存しておくことができます。 immune.combinedという変数に代入したオブジェクトをimmune_combined.rdsという名前で保存する場合、以下のようにします。

saveRDS(immune.combined, file = "~/Projects/datasets/immune_combined.rds")

因みに読み込むときはloadRDS関数を使用します。

immune.combined <- loadRDS(file = "~/Projects/datasets/immune_combined.rds")