こんにちはtsuyuhです。

最近ますます寒くなってきて、1月・2月は一体どうなってしまうんだろうと戦々恐々としています。

今回は、シングルセルRNA-seq解析によく用いられるRパッケージ "Seurat" のFindMarkers() 関数についての話題です。
FindMarkers() は二群間での発現量の有意差を調べることができる関数です。
その解析結果はVolcano plotやMA plotにしたいところですね。
しかしながら、この関数はデフォルトのオプションではプロットを描くのに適していません。
それについて少しだけ解説します。

やりたいこと

Seuratオブジェクト "seurat_obj" に含まれる "cell" という種類の細胞について、"sampleA" と "sampleB" 間での二群間比較をおこない、その結果をVolcano plotで可視化したいと思います。

なお、Volcano plotを綺麗に描画する方法は色々ありますが、今回は単にplot() で描きます。

作業環境

R v4.0.3
Seurat v4.0.5

結論

> library(Seurat)

# 細胞種を推定するまでは省略

# 全細胞から細胞種 "cell" を抽出
> cell <- subset(seurat_obj, celltype == "cell")

# FindMarkers関数でsampleAとsampleBの二群間比較を実行
> result <- FindMarkers(cell, ident.1 = "sampleA", ident.2 = "sampleB", group.by = "orig.ident", min.pct = 0.25, logfc.threshold = 0)

# 二群間比較の結果を使ってVolcano plotを描く
> volcano <- cbind(result$avg_logFC, -log10(result$p_val_adj))
> colnames(volcano) <- c("logFC", "-logPval")
> pdf("volcano_plot.pdf")
> plot(volcano)
> dev.off()

解説

FindMarkers関数は、デフォルトでlogFC >= 0.25の遺伝子だけを抽出してからDEG計算をおこなっているので、普通にやるとresult変数に格納される遺伝子はlogFC >= 0.25のものだけになります。
なので、この結果をそのまま使っても、全遺伝子でvolcano plotは描けません。

そこで、logfc.threshold = 0 というオプションを追記し、logFCではフィルタリングしないようにします。

※ logFC = log2{ident.1の各細胞の(発現を正規化した値)の平均} - log2{ident.2の各細胞の(発現を正規化した値)の平均}」

オプション詳細

min.pct = xx
少なくともどちらか一方のサンプル群において、xxの割合以上の細胞で検出された遺伝子をDEG解析対象にする。デフォルトは0.1。これを設定しないと、両方のサンプル群で発現が稀な遺伝子をDEGとして検出してしまい、生物学的に謎な結果になる。
logfc.threshold = yy
logFCがyy以上の遺伝子をDEG解析対象にする。デフォルトは0.25。 logfc.threshold = 0 にしないと、全遺伝子がVolcano plotに反映されない。ただし、解析にかかる時間は爆増する。

他の色々なオプションは公式のページをご参照ください↓
https://satijalab.org/seurat/reference/findmarkers