こんにちは、受託コンサルティングチームの hosor です。
もう12月ということで、時が過ぎるのが早くてなんだか悲しい気持ちです。
シングルセル解析において、様々なサンプルを統合した場合、各サンプルにどのクラスタの細胞がどれくらい含まれているのかを調べることがあると思います。
今回は、視覚的に分かりやすくするために、細胞の割合についての積み上げ棒グラフを描画していこうと思います。
作業環境
- CentOS Linux release 7.4.1708 (Core)
- R v4.3.1
- Seurat v4.3.0.1
- SeuratData v0.2.2
- ggplot2 v3.4.3
- reshape2 v1.4.4
- scales v1.2.1
使用データ
Human Cell Atlas によって作成された、8 人のドナーサンプルから得られたヒト骨髄単核球細胞 (BMNC) を使用しました。*1
サンプルごとに分割して、クラスタリング結果を UMAP 上で確認してみると 21 個のクラスタに分かれることが見て取れました。
MantonBM1, MantonBM2, ... が各ドナーサンプルを表しています。
例えば、MantonBM3 ではクラスタ 0 の細胞が少ない一方で、クラスタ 7, 8 の細胞が多いということが分かります。
こういった偏りを視覚的により分かりやすくするために、積み上げ棒グラフを描画していきます。
スクリプト
作成した Seurat オブジェクト hcabm40k を用いて、各サンプルにどのクラスタの細胞がどれくらい含まれているのか、といったことを割合で計算していきます。
※ 2023/12/15 追記
prop.table(cross_tabulation, margin = 1)の結果に誤りがありましたため、修正いたしました。
cross_tabulation <- table(hcabm40k$orig.ident, hcabm40k$seurat_clusters)
cross_tabulation
0 1 2 3 4 5 6 7 8 9 10 11 12
MantonBM1 1276 678 463 586 468 153 154 10 3 98 231 0 91
MantonBM2 609 1170 377 709 661 82 734 13 0 26 134 1 56
MantonBM3 9 378 210 120 414 125 145 1469 777 90 11 665 23
MantonBM4 362 453 1099 461 430 197 515 3 1 162 61 0 178
MantonBM5 126 645 771 267 482 576 93 380 16 179 25 18 57
MantonBM6 1045 460 601 689 191 825 253 0 0 82 45 1 80
MantonBM7 1195 548 386 635 347 483 401 5 1 45 193 1 133
MantonBM8 538 429 656 833 1089 415 126 42 30 105 42 9 64
13 14 15 16 17 18 19 20
MantonBM1 85 19 24 26 10 6 3 8
MantonBM2 72 108 8 9 7 20 12 0
MantonBM3 36 56 38 32 14 0 8 5
MantonBM4 73 152 78 72 19 9 8 4
MantonBM5 105 69 49 48 18 0 3 4
MantonBM6 49 20 36 31 15 4 7 5
MantonBM7 80 61 47 31 5 11 1 3
MantonBM8 51 10 47 40 5 0 5 1
# margin = 1: 各行の合計が 1.0 になるように計算
# margin = 2: 各列の合計が 1.0 になるように計算
cross_tabulation_ratio <- prop.table(cross_tabulation, margin = 1)
cross_tabulation_ratio
0 1 2 3 4 5
MantonBM1 0.290528233 0.15437158 0.10541894 0.13342441 0.10655738 0.03483607
MantonBM2 0.126663894 0.24334443 0.07841098 0.14746256 0.13747920 0.01705491
MantonBM3 0.001945946 0.08172973 0.04540541 0.02594595 0.08951351 0.02702703
MantonBM4 0.083467835 0.10445008 0.25340097 0.10629467 0.09914688 0.04542310
MantonBM5 0.032052913 0.16408039 0.19613330 0.06792165 0.12261511 0.14652760
MantonBM6 0.235413381 0.10362694 0.13539085 0.15521514 0.04302771 0.18585267
MantonBM7 0.259106678 0.11882047 0.08369471 0.13768430 0.07523851 0.10472680
MantonBM8 0.118580560 0.09455587 0.14458894 0.18360150 0.24002645 0.09147013
6 7 8 9 10
MantonBM1 0.03506375 0.0022768670 0.0006830601 0.022313297 0.052595628
MantonBM2 0.15266223 0.0027038270 0.0000000000 0.005407654 0.027870216
MantonBM3 0.03135135 0.3176216216 0.1680000000 0.019459459 0.002378378
MantonBM4 0.11874568 0.0006917224 0.0002305741 0.037353009 0.014065022
MantonBM5 0.02365810 0.0966675146 0.0040702111 0.045535487 0.006359705
MantonBM6 0.05699482 0.0000000000 0.0000000000 0.018472629 0.010137418
MantonBM7 0.08694709 0.0010841284 0.0002168257 0.009757155 0.041847355
MantonBM8 0.02777166 0.0092572184 0.0066122989 0.023143046 0.009257218
11 12 13 14 15
MantonBM1 0.0000000000 0.020719490 0.019353370 0.004326047 0.005464481
MantonBM2 0.0002079867 0.011647255 0.014975042 0.022462562 0.001663894
MantonBM3 0.1437837838 0.004972973 0.007783784 0.012108108 0.008216216
MantonBM4 0.0000000000 0.041042195 0.016831911 0.035047268 0.017984782
MantonBM5 0.0045789875 0.014500127 0.026710761 0.017552786 0.012465022
MantonBM6 0.0002252760 0.018022077 0.011038522 0.004505519 0.008109935
MantonBM7 0.0002168257 0.028837814 0.017346054 0.013226366 0.010190807
MantonBM8 0.0019836897 0.014106238 0.011240908 0.002204100 0.010359268
16 17 18 19 20
MantonBM1 0.005919854 0.002276867 0.0013661202 0.0006830601 0.0018214936
MantonBM2 0.001871880 0.001455907 0.0041597338 0.0024958403 0.0000000000
MantonBM3 0.006918919 0.003027027 0.0000000000 0.0017297297 0.0010810811
MantonBM4 0.016601337 0.004380908 0.0020751672 0.0018445930 0.0009222965
MantonBM5 0.012210633 0.004578988 0.0000000000 0.0007631646 0.0010175528
MantonBM6 0.006983555 0.003379139 0.0009011039 0.0015769317 0.0011263798
MantonBM7 0.006721596 0.001084128 0.0023850824 0.0002168257 0.0006504770
MantonBM8 0.008816399 0.001102050 0.0000000000 0.0011020498 0.0002204100
次に、データの整形をしていきます。
dataframe4barplot <- melt(cross_tabulation_ratio)
colnames(dataframe4barplot) <- c("sample", "cluster", "percentage")
dataframe4barplot$cluster <- factor(dataframe4barplot$cluster, levels = 0 : max(as.numeric(hcabm40k$seurat_clusters)))
# サンプル数 * クラスタ数のテーブルデータを、「あるサンプルにおける、あるクラスタの細胞割合」を行ごとに表すテーブルデータに変換した
head(dataframe4barplot)
sample cluster percentage
1 MantonBM1 0 0.290528233
2 MantonBM2 0 0.126663894
3 MantonBM3 0 0.001945946
4 MantonBM4 0 0.083467835
5 MantonBM5 0 0.032052913
6 MantonBM6 0 0.235413381
最後に、dataframe4barplot をもとに描画していきます。
# グラフの設定
theme <- theme(panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_text(colour = "black"),
text = element_text(size = 15),
axis.text.x = element_text(angle = 90, hjust = 1),
axis.title.x = element_text(size = 18),
axis.title.y = element_text(size = 18)
)
g <- ggplot(dataframe4barplot, aes(x = sample, y = percentage, fill = cluster)) # 描画する変数
g <- g + geom_bar(stat = "identity", width = 0.8) # 棒グラフ
g <- g + scale_y_continuous(expand = c(0, 0), limits = c(0, 1.02), labels = percent) # パーセント表示
g <- g + theme # 上記のグラフの設定
ggsave("BarPlot_cell_each_sample.pdf", width = length(unique(hcabm40k$orig.ident)) * 2)
各サンプルにおける、各クラスタの細胞割合を積み上げ棒グラフで描画することができました！
先述したように、MantonBM3 の各クラスタの細胞割合が、他のサンプルとは異なることが一目瞭然ですね。
少しスクリプトを変えるだけで、「各クラスタにおける、各サンプルの細胞割合」を表現することもできます。
cross_tabulation <- table(hcabm40k$orig.ident, hcabm40k$seurat_clusters)
cross_tabulation_ratio <- prop.table(cross_tabulation, margin = 2) # margin = 1 -> margin = 2
dataframe4barplot <- melt(cross_tabulation_ratio)
colnames(dataframe4barplot) <- c("sample", "cluster", "percentage")
dataframe4barplot$cluster <- factor(dataframe4barplot$cluster, levels = 0 : max(as.numeric(hcabm40k$seurat_clusters)))
# axis.text.x = element_text(angle = 90, hjust = 1) を消して、x 軸ラベルを回転させない
theme <- theme(panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_text(colour = "black"),
text = element_text(size = 15),
axis.title.x = element_text(size = 18),
axis.title.y = element_text(size = 18)
)
g <- ggplot(dataframe4barplot, aes(x = cluster, y = percentage, fill = sample)) # x と fill の変数を入れ替える
g <- g + geom_bar(stat = "identity", width = 0.8)
g <- g + scale_y_continuous(expand = c(0, 0), limits = c(0, 1.02), labels = percent)
g <- g + theme
ggsave("BarPlot_cell_each_cluster.pdf", width = max(as.numeric(hcabm40k$seurat_clusters))) # pdf の名前や幅を変更
R パッケージ Seurat が出力するグラフに似たデザインになるようにこだわったので、ご活用いただければ幸いです！
