Zenn
😊

[R][Seurat][ggplot] scRNAseq plot tips

2023/02/16に公開・約13,800字
R
#
ggplot
#
seurat
#
bioinfomatics
#
scrnaseq
tech

230221 今後も定期的に追記していく予定

目次

シングルセルデータ解析ツールのSeuratには多彩なplotが用意されているが、各plotに用意されているオプションでは不十分に感じることがある。ここではSeurat plotをより自由に表現するtipsを紹介する。

必要なライブラリの読み込み
library(Seurat)
library(dplyr)
library(ggplot2)
library(patchwork)
library(pheatmap)
library(RColorBrewer)
library(viridis)

データの準備

10X GENOMICSが提供しているPBMC 1万細胞のデータを使用
https://support.10xgenomics.com/single-cell-gene-expression/datasets/3.0.0/pbmc_10k_v3?
上記ページのFeature/cell matrix HDF5(filtered)のリンクからh5ファイルをダウンロード。

mat <- Read10X_h5(filename = "pbmc_10k_v3_filtered_feature_bc_matrix.h5")
sc <- CreateSeuratObject(mat, min.cells = 3, min.features = 200)
sc <- PercentageFeatureSet(sc, pattern = "^MT-", col.name = "Percent_mito")

sc <- sc %>% 
  NormalizeData() %>% 
  FindVariableFeatures() %>% 
  ScaleData() %>% 
  RunPCA() %>% 
  RunUMAP(dims = 1:30) %>% 
  FindNeighbors(dims = 1:30) %>% 
  FindClusters(resolution = 0.5)

# resolution 0.5の時のcluster情報をmeta.dataスロットに保存
sc$cluster0.5 <- sc@active.ident

# cluster数が少ない版も保存しておく。
sc <- FindClusters(sc, resolution = 0.1)
sc$cluster0.1 <- sc@active.ident

軸やラベル

Serutaには凡例や軸を無くすような簡単な修飾コマンドが用意されている。

p1 <- DimPlot(sc)
p2 <- DimPlot(sc) + NoLegend() # 凡例を無くす
p3 <- DimPlot(sc) + NoLegend() + NoAxes() # さらに軸を無くす
p1|p2|p3

FontSize機能で、軸のラベルやタイトルのサイズを指定できる。plotのmain titleがある場合はmain=〇〇でサイズ指定可能。

DimPlot(sc) + 
	FontSize(x.text = 4, y.text = 8, x.title = 12, y.title = 16)

x軸ラベルを45度傾ける

DotPlot(sc, features = c("CD3E", "CD4", "CD8A", "MS4A1", "CD79B", "CD14", "CD68")) + 
	RotatedAxis()

ggplotの修飾コマンドを使用

Seuratでは上記のように簡単なplot修飾コマンドを用意してくれているが、それでは不十分であったり、Seurat版のコマンドを探すのが面倒だったりする。
基本はggplot2パッケージを使ったplotなので、ggplotの修飾コマンドが機能する。

以下は適当な例だが、ggplotの記法を使えばplotの各要素 (線やtext, legend,,,)などが自由度高く編集可能である。

DimPlot(sc) + 
  theme(axis.title = element_text(size = 10, face = "bold"),
        axis.text.x = element_text(size = 6, colour = "red"),
        axis.text.y = element_text(size = 6, colour = "blue"),
        axis.line = element_blank(),
        axis.ticks = element_line(colour = "gray"),
        legend.position = "bottom") +
  ggtitle("UMAP")

表示目盛りを細かくしてfiltering閾値を決めやすくする。

デフォルトの目盛り幅だと、フィルタリングの際の閾値を決めにくい。以下のFeatureScatterの例だと、nCount_RNAの0の次の目盛りが20000となってしまっている。

FeatureScatter(sc, feature1 = "nCount_RNA", feature2 = "Percent_mito", group.by = "orig.ident")

ggplotのscale_x_continuousのbreaks引数で表示目盛りを細かくすると良い。これなら2000-3000ぐらいを下限の閾値に設定すればいいことがわかる。

FeatureScatter(sc, feature1 = "nCount_RNA", feature2 = "Percent_mito", group.by = "orig.ident") + 
  scale_x_continuous(breaks = seq(0,max(sc$nCount_RNA), by=1000)) + 
  scale_y_continuous(breaks = seq(0,max(sc$Percent_mito), length.out = 20)) + 
  theme(axis.text.x = element_text(size = 6, angle = 90, hjust = 1),
        axis.text.y = element_text(size = 6))

seqのlength.out引数は指定範囲を指定要素数で分割してくれる。

離散値の色 (VlnPlotやDimPlot)

デフォルトの色

DimPlot(sc, group.by = "cluster0.1", label = T)

各クラスターの色を変更するにはDimPlotのcols引数で色を指定する。group.by引数で指定した情報で群分けした場合の群の数だけ色の指定が必要。色数が足りないとエラーになる。

色ベクトルの用意に関しては、私のスクラップにまとめているのでそちらを参照して欲しい。
https://zenn.dev/rchiji/scraps/0895ffa42d68ba

連続値の色 (FeaturePlotやDotPlot, DoHeatmap)

FeaturePlotの色

遺伝子発現量を可視化するFeaturePlotにもcols引数があるが、困ったことにcolsの色数が増えると凡例に表示される色の範囲がおかしくなる。正規化して対数変換をした遺伝子発現量をplotしているのでこの値が変わるのはよろしくない。

default
FeaturePlot(sc, features = c("CD3D", "MS4A1", "CD14"), ncol = 3)

2 colors
FeaturePlot(sc, features = c("CD3D", "MS4A1", "CD14"), ncol = 3, 
	cols = c("lightgray", "darkred"))

over 3 colors
FeaturePlot(sc, features = c("CD3D", "MS4A1", "CD14"), ncol = 3, 
	cols = c("lightgray", "gray", "orange", "red", "darkred", "black"))

これを解決するには、ggplotによる色修飾機能を使用する。

scale_color_gradientnで複数色を連続値に割り当てる。
FeaturePlot(sc, features = c("CD3D", "MS4A1", "CD14"), ncol = 3) & 
	scale_color_gradientn(colours = c("lightgray", "gray", "orange", "red", "darkred", "black"))

scale_color_gradientnで複数色を連続値に割り当てる。
FeaturePlot(sc, features = c("CD3D", "MS4A1", "CD14"), ncol = 3) & 
	scale_color_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))

viridisのパレットもscale_color_viridis機能から使用可能。(※ scale_color_brewerでRColorBrewerのパレットが使えるが、これは離散値用)

scale_color_viridis
FeaturePlot(sc, features = c("CD3D", "MS4A1", "CD14"), ncol = 3) &
	scale_color_viridis(option = "D")

DotPlotの色

DotPlotではAverage Expressionの値は変わらないが、複数色指定しても最初の2色しか使われてない。

default
DotPlot(sc, features = c("CD3D", "MS4A1", "CD14"))

over 3 colors
DotPlot(sc, features = c("CD3D", "MS4A1", "CD14"), 
	cols = c("lightgray", "orange", "red", "darkred", "black"))

DotPlotもggplotの色修飾機能で変色可能。

scale_color_gradientnで複数色を連続値に割り当てる。
DotPlot(sc, features = c("CD3D", "MS4A1", "CD14")) & 
	scale_color_gradientn(colours = c("lightgray", "orange", "red", "darkred", "black"))

scale_color_viridis
DotPlot(sc, features = c("CD3D", "MS4A1", "CD14")) & 
	scale_color_viridis()

DoHeatmapの色

default
DoHeatmap(sc, features = c("CD3D", "MS4A1", "CD14"))

デフォルトでscale.dataスロットが使用するため0を平均に正負の値が表示される。

slot引数で"data"を指定すれば、FeaturePlotやDotPlotと同様の尺度で表示されるが、色合いの悪さが目立つようになる。

data slot
DoHeatmap(sc, features = c("CD3D", "MS4A1", "CD14"), slot = "data")

DoHeatmap機能にはそもそも色を指定するオプションが無いが、ggplotの色修飾で編集可能である。なお、FeaturePlot/DotPlotではscale_color_シリーズであったが、DoHeatmapではscale_fill_シリーズを使用する。

2 colors
DoHeatmap(sc, features = c("CD3D", "MS4A1", "CD14")) + 
	scale_fill_gradient(low = "grey", high = "red")

3 colors
DoHeatmap(sc, features = c("CD3D", "MS4A1", "CD14")) + 
	scale_fill_gradient2(low = "blue",mid = "white",high = "red", midpoint=0)

over 3 colors
DoHeatmap(sc, features = c("CD3D", "MS4A1", "CD14")) + 
	scale_fill_gradientn(colors = RColorBrewer::brewer.pal(7, "RdYlBu"))




ちなみにheatmapの上部に表示されるcluster名がデフォルトだとfontsizeが大きくて、細胞数が少ないクラスターの所では文字が重なってしまう。DoHeatmapのsize引数でここのfontsizeだけを指定することができる。

font size
DoHeatmap(sc, features = c("CD3D", "MS4A1", "CD14"), size = 4)

値の範囲

外れ値のような発現を持つ細胞がいるとそっちに色が取られて、どこまでが発現細胞なのか見づらい時がある。

FeaturePlot(sc, features = c("CD3D","CD4","CD8A"), ncol = 3) & 
	scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))

色を割り当てる範囲をlimitsで狭める。

上限より大きい値を持つ細胞には色が割り当てられないので灰色に表示されてしまうが。

 FeaturePlot(sc, features = c("CD3D","CD4","CD8A"), ncol = 3) & 
	 scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")), 
		 limits = c(0,1))

与える色ベクトルの最後の色を繰り返して、ある値以上は同じ色を持つようにする。

col <- rev(brewer.pal(n = 11, name = "RdBu"))
col <- c(col, rep(col[length(col)], times = 11)) # timesの回数は好きに決定
 FeaturePlot(sc, features = c("CD3D","CD4","CD8A"), ncol = 3) & 
	 scale_colour_gradientn(colours = col)

FeaturePlotのorder引数をTRUEにすると陽性細胞を前面に持って来れる。

FeaturePlot(sc, features = c("CD3D","CD4","CD8A"), ncol = 3, order =T) & 
	scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))

VlnPlot

細胞数が多い時のVlnPlotはpt.size = 0を指定

細胞数が多いデータの場合、点が重なりすぎてviolin plotが見づらくなる。

VlnPlot(sc, features = c("nFeature_RNA", "nCount_RNA", "Percent_mito"), group.by = "orig.ident")


VlnPlot(sc, features = c("nFeature_RNA", "nCount_RNA", "Percent_mito"), pt.size = 0, group.by = "orig.ident")

本当は透明度などもいじれれば細胞数が多くても点を残しても良いのだが、VlnPlot機能には無いので、ggplotのgeom_jitterを使用するのもあり。

VlnPlot(sc, features = c("nFeature_RNA", "nCount_RNA", "Percent_mito"), pt.size = 0, group.by = "orig.ident") &
  geom_jitter(size = 0.1, height = 0, alpha = 0.2)

見たい遺伝子数が多い時はstack引数が有用

VlnPlot(sc, features = c("CD3D","CD4","CD8A","PTPRC","MS4A1","CD19","CD79A","CD14","CD68","ITGAX"), stack = T)

さらにflip引数もTRUEにした方が見栄えが良い。

VlnPlot(sc, features = c("CD3D","CD4","CD8A","PTPRC","MS4A1","CD19","CD79A","CD14","CD68","ITGAX"),
        stack = T, flip = T)

DimPlot

cells.highlight

特定の細胞のみを色付きでplotすると、どこに分布しているのかを確認しやすい。
DimPlotのcells.highlight引数に細胞バーコードのベクトルを指定する。

DimPlot(sc, cells.highlight = c("AAATGGATCGTAGGGA-1"))

WhichCells

特定の条件を満たす細胞バーコードベクトルを取得するにはWhichCells機能が有用である。

cluster 1と3の細胞バーコードをハイライト表示
DimPlot(sc, 
	cells.highlight = WhichCells(sc, idents = c("1","3")))

WhichCellsのidents引数は、active.identスロットに保存された名前を対象にする。
active.identの内容はlevels関数などで確認可能である。

active.identの確認
levels(sc)

[1] "0" "1" "2" "3" "4" "5" "6" "7" "8" "9"

active.identの変更

meta.dataスロットに保存している他のメタ情報をactive.identに設定するには以下のようにIdents(sc)にmeta.dataスロットの列名を指定するか、SetIdent機能を使用する。

Idents機能。resolution=0.5でFindClusterした結果を上書き
Idents(sc) <- "cluster0.5"
levels(sc)

[1] "0" "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15" "16" "17" "18" "19" "20"

SetIdent機能。
sc <- SetIdent(sc, value = "cluster0.1")
levels(sc)

[1] "0" "1" "2" "3" "4" "5" "6" "7" "8" "9"



WhichCellsのexpression引数を使えば、active.identを変えずにmeta.dataスロットを参照できる。
expression = meta.data列名 %in% 欲しいlevel と書く。meta.data列名は""が要らないのと、levelが複数ある場合は==ではなく%in%とする必要がある。

active.ident以外のmeta.dataスロットの情報を使用する。
Idents(sc) <- "cluster0.5"
DimPlot(sc, 
	cells.highlight = WhichCells(sc, expression = cluster0.1 %in% c("1","3")))

またexpression引数を使えば、遺伝子発現量の閾値で細胞バーコードを取り出すことができる。

CD8A発現を有する細胞をハイライト表示
DimPlot(sc, cells.highlight = WhichCells(sc, expression = CD8A>0))

CD8A発現を有する細胞かつGZMA発現が1より大きい細胞をハイライト表示
DimPlot(sc, cells.highlight = WhichCells(sc, expression = CD8A>0 & GZMA>1))

複数ハイライトを色分けして表示

WhichCellsで複数条件を満たす細胞バーコードを取り出しても、1つのハイライトになってしまう。
各条件で複数ハイライトを作りたければcells.highlight引数に細胞バーコードのリストを指定する。
listの要素に名前をつけていなければ、legendがGroup_1, Group_2,,,になる。
cols.highlight引数でlist数に対応した色ベクトルを用意しなければ全て同じハイライト色になるので注意。

DimPlot(sc, 
        cells.highlight = list(cluster1 = WhichCells(sc, expression = cluster0.1 %in% "1"),
                               cluster3 = WhichCells(sc, expression = cluster0.1 %in% "3")),
        cols.highlight = c("darkred", "darkblue")
        )

DoHeatmap

plot marginを調整

ヒートマップ上部にラベルされるクラスター名が長いと、途中で途切れてしまう。

Seurat PBMCチュートリアルで作成したSeuratオブジェクトとDEGを使用
DEGs <- FindAllMarkers(pbmc, only.pos = T, logfc.threshold = 1)
DEGs %>%
     group_by(cluster) %>%
     top_n(n = 5, wt = avg_log2FC) -> top5
DoHeatmap(pbmc, features = top5$gene)

ggplot2のplot margin調節を使用すると改善可能。

DoHeatmap(pbmc, features = top5$gene) + 
	theme(plot.margin = unit(c(1,0,0,0),"cm"))

Discussion

ログインするとコメントできます