いろいろ倉庫

KNIME、EXCEL、R、Pythonなどの備忘録

【R】シングルセル解析を学んでみたい(Seurat)②

・お題:細胞一つ一つのオミクス情報を解析する手法をシングルセル解析と呼ぶらしい。Seuratというライブラリを使ってやってみたい。

・私は素人なので、以下のチュートリアルをなぞってみたい。

satijalab.org

・前回、チュートリアルのPCAのヒートマップを描くところまで行った。Seuratでは、PCAにより削減した次元を用いて続きの解析を進めていくらしい。このステップで大事になるのは、データセットの次元数。

・JackStraw法でやってみる。データのサブセットをランダムに並べ替えてPCAをやり直し、特徴スコアの帰無分布を構築し、この手順を繰り返すらしい。重要な主成分は、p値が低くなる。

>pbmc <- JackStraw(pbmc, num.replicate = 100)

・私の実行環境では5分程度かかった。次にスコアを算出し、図示してみる。

> pbmc <- ScoreJackStraw(pbmc, dims = 1:20)

>JackStrawPlot(pbmc, dims = 1:15)

・破線が一様分布なので、重要でないほど寝てくるらしい。上位6つの主成分がずぬけて大事っぽい。上位10個の主成分までが少し立っており、それ以降はずいぶん寝ているようにも見える。公式では、10-12程度が良さげっぽい。

・主成分分析の次元数の選び方といえば、エルボープロットも有名と思う。これもSeuratに実装されているので、やってみる。

> ElbowPlot(pbmc)

・見た感じPC8辺りで寝ているように見える。公式を見てみると、PC10ぐらいに肘があるとのこと。私の見方が厳しいのかも。。

・これらの結果として、7~12の主成分ならおおむね妥当なので、今回は10を設定する。ただし、差の小さい細胞集団を切り分けたい場合などは当然細かい違いを見ていくことになるので、目的によってどこまでとるか考えなければならないし、何なら主成分の数を振って解析するのが良いらしい。

・次に、クラスタリングしてみる。クラスタリングの考え方に関しては、公式サイトを参照。Jaccard 類似度に関しては、以下のサイトが分かりやすいと思う。

mieruca-ai.com

> pbmc <- FindNeighbors(pbmc, dims = 1:10)

Computing nearest neighbor graph
Computing SNN
> pbmc <- FindClusters(pbmc, resolution = 0.5)

・各細胞がどんな感じにクラスタ分けされたか見てみる。

> head(Idents(pbmc), 5)

・次に、図示する。よく論文で目にする七色のプロットがついに描ける…!次元削減手法はUMAPを選択する。次元数は先のクラスタ解析と同じ次元数にしておく。

> pbmc <- RunUMAP(pbmc, dims = 1:10)

> DimPlot(pbmc, reduction = "umap")

・次に、発現変動している特徴量を探してみる。クラスターを定義する遺伝子発現バイオマーカーを定義するともいえる。比較の仕方など、いろいろオプションがあるみたい。

・試しに、クラスター3のマーカーを探してみる。ident.1に差が見たいクラスタを指定すると、デフォルトで他のすべてのクラスタと比較してくれるっぽい。

>cluster3.markers <- FindMarkers(pbmc, ident.1 = 3, min.pct = 0.25)
>head(cluster3.markers, n = 5)

・結果を図示してみる。

> VlnPlot(pbmc, features = c("MS4A1", "CD79A", "CD79B"))

> FeaturePlot(pbmc, features = c("MS4A1", "CD79A", "CD79B"))

・確かに、クラスター3でのみ発現が変わっている遺伝子が釣れたっぽい。

・次に、クラスター5 vs クラスター0とクラスター3で差がある遺伝子を探してみる。この場合、ident.2に比較対象のクラスター番号を入れる。

> cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)

> head(cluster5.markers, n = 5)

> VlnPlot(pbmc, features = c("FCGR3A", "IFITM3", "CFD"))

・確かに、クラスター5とクラスター0&クラスター3で差があるっぽい遺伝子が出てきた。

・次に、各クラスターごとに特徴的な遺伝子を一気に調べてみた。ちなみに、only.posをTRUEにすると、positive markerだけ返してくれる(デフォルトはFALSE)。 

>pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

> pbmc.markers %>%
+   group_by(cluster) %>%
+   slice_max(n = 2, order_by = avg_log2FC)

・引数test.useで発現量の比較の仕方を変えられる(詳細はDE vignetteを参照)。例えば、ROCテストは、個々のマーカーのclassification powerを返すらしい。釣ってきた遺伝子をどのような用途に利用したいかによって弄ると良さそう。コードは割愛。
・各クラスタのマーカー遺伝子の発現を細胞ごとのヒートマップで表すこともできる。

> VlnPlot(pbmc, features = cluster0.markers %>% row.names %>% head(,2))
> pbmc.markers %>%
+   group_by(cluster) %>%
+   top_n(n = 10, wt = avg_log2FC) -> top10
> DoHeatmap(pbmc, features = top10$gene) + NoLegend()

・ここで「このクラスタ分類は、いったいどの細胞なのか?」という疑問がわいてきた。遺伝子発現パターン(マーカー遺伝子)が分かっている細胞なら、その発現に応じて細胞種のラベルを付与してやれば良い。今回はすでに細胞情報が与えられているとして、図示するときに細胞種を示してやる。

> new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")#クラスターごとのID(細胞種)を与える。クラスタの番号順に細胞種を与える。
> names(new.cluster.ids) <- levels(pbmc)#細胞種のベクトルの成分名にクラスタ番号を与える。

> new.cluster.ids

> pbmc <- RenameIdents(pbmc, new.cluster.ids)#identity classeを付けなおす。

> DimPlot(pbmc, reduction = "umap", pt.size = 0.5)

> DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()#凡例ではなく、プロット上に直接ラベルを記載。

・結果を保存するには、以下のコマンドを実行すれば良いらしい。file名は適宜変更。。
>saveRDS(pbmc, file = "../output/pbmc3k_final.rds")

 

・七色のお絵かきをできるのは楽しいけれど、ここから有用な情報を抽出できるように頑張りたい今日この頃。。

 

 

おわり。