Seurat-如何進行特定細胞的分群subcluster

前言

在許多單細胞數據研究的學術文章中不難發現,大多數研究人員將所有細胞群(cluster)做命名後,在後續的分析過程中,會再次將有興趣的細胞群獨立出來,再一次對該細胞群進行分群(subclustering),以便對該細胞群做一個更細緻的探討。不同的分析人員針對分群(subclustering)會有不同的分析方式,沒有所謂的對與錯,只要邏輯正確,想法正確,符合自身研究主題的需求,那我認為該分析方法與步驟就可以說是正確的。以下是我對於特定細胞分群(subclustering)的分析步驟與方法,如有不正確的地方,歡迎各位留言指教。

幫自己的文章做宣傳,如果大家有興趣的話,可以去看看我寫的Seurat教學:

範例數據

在這次的範例中,我們將分析來自10X Genomics的周邊血(PBMC)的兩組單細胞數據集,分別包含了2700和5025顆單細胞的數據集,原始數據的.rds儲存檔可以在我的雲端下載。如果對檔案有疑慮,可以到10X Genomics數據庫自行下載,然後按照單細胞數據分析:批次效應與Seurat Integration分析詳解這篇文章進行分析。

數據集下載位置:

3k PBMCs from a Healthy Donor

5k Peripheral Blood Mononuclear Cells (PBMCs) from a Healthy Donor (v3 chemistry)

1. 加載套件和載入.rds儲存檔

加載所需套件進入R環境,並將pbmc3k5K_final.rds數據載入R中,這裡使用的指令是readRDS()函數,代碼中提供的路徑只是一個示例路徑,需要根據電腦實際文件所在的路徑進行修改

# 加載套件
library(Seurat)
library(dplyr)
library(patchwork)

# 需用電腦中的正確路徑,須將"\"更改為"/"
pbmc.combined <- readRDS(file = "C:/Users/Administrator/Desktop/pbmc3k5K_final.rds")

# 檢查輸入結果
DimPlot(pbmc.combined, reduction = "umap", group.by = "orig.ident")
DimPlot(pbmc.combined, reduction = "umap")
DimPlot(pbmc.combined, reduction = "umap",label = TRUE, split.by = "orig.ident")

由於這套分析結果的解析度(Resolution)分得太細,不適合做後面的分群分析,因此我們先對結果做一些調整,將解析度降低,“假裝”我們得到的結果就是下面的結果:

# 降低解析度,以配合後續分析
pbmc.combined <- FindClusters(pbmc.combined, resolution = 0.1)
DimPlot(pbmc.combined, reduction = "umap",label = T)
recluster-1

再次提醒,這只是為了教學而故意降低解析度,在實際分析單細胞數據時,請根據自己的數據調整至最適合的解析度!!!

2. 提取想要分群的細胞群

接下來我們進一步提取有興趣想要往下分析的細胞,假設我們對第1群細胞群特別特別感興趣,可以用subset()函數將該細胞提取出來。

# 提取編號為0的細胞群
cluster1 <- subset(x = pbmc.combined,idents= c("1"))

# 檢查提取結果
DimPlot(cluster1, reduction = "umap", label = T)
subclustering

稍微檢查一下提取出來的細胞群是否正確,沒問題的話就能進行下一步了。

3. 進行數據整合(Seurat integration)以及PCA分析

接下來的分析與Seurat integration的基本流程一樣,先將assay轉到RNA assay,用原始RNA數據重跑第1群細胞群的數據整合(Seurat integration),進行數據縮放(scaling data)、PCA降維分析、確定數據維度以及進行細胞聚集分析(cell clustering)。

# 將assay轉為RNA,用原始RNA數據重跑一次數據整合 (Seurat integration)
DefaultAssay(cluster1) <- "RNA"
cluster1

# 設定數據分割列表,進行標準化以及選擇高度變異的特徵的分析
cluster1.list<- SplitObject(cluster1, split.by = "orig.ident")
cluster1.list
cluster1.list <- lapply(X = cluster1.list, FUN = function(x) {
  x <- NormalizeData(x,normalization.method = "LogNormalize", scale.factor = 10000)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

# 選擇跨數據集重複變化的特徵進行整合
features <- SelectIntegrationFeatures(object.list = cluster1.list)
# 整合分析
cluster1.anchors <- FindIntegrationAnchors(object.list = cluster1.list, dims = 1:15, anchor.features = features)
# 創建一個integrated數據分析
cluster1.combined <- IntegrateData(anchorset = cluster1.anchors, dims = 1:15)

# 進行整合數據分析
# 設定使用先前整合的integrated assay進行分析。
# 原始未修改的數據仍然存在於RNA assay中,可以用DefaultAssay(cluster1.combined) <- "RNA"進行轉換
DefaultAssay(cluster1.combined) <- "integrated"
cluster1.combined <- ScaleData(cluster1.combined, verbose = TRUE)
cluster1.combined <- RunPCA(cluster1.combined, npcs = 50, verbose = TRUE)
#Dertermine the dimensionality of the dataset
ElbowPlot(cluster1.combined, ndims= 50)
# Cluster the cells (1:15)
cluster1.combined <- FindNeighbors(cluster1.combined, reduction = "pca", dims = 1:15)
cluster1.combined <- FindClusters(cluster1.combined, resolution = 0.2)
table(cluster1.combined@active.ident)
table(Idents(cluster1.combined), cluster1.combined$orig.ident)
subclustering

以解析度(Resolution)為0.2來分析的話,可以將第1群細胞群再次細分為3個子群。子群的數量可以經由resolution()函數的大小來調整,也可以用clustree來進行判斷,詳情請參考Clustree (0.5.0)-聚類(clustering)層次化的好物!

4. 非線性降維分析法(UMAP/tSNE)

接下來便是進行可視化非線性降維分析了,圖示以UMAP圖為例子。

# UMAP (1:15)
cluster1.combined <- RunUMAP(cluster1.combined, reduction = "pca", dims = 1:15)
DimPlot(cluster1.combined, reduction = "umap", group.by = "orig.ident")
DimPlot(cluster1.combined, reduction = "umap")
DimPlot(cluster1.combined, reduction = "umap", split.by = "orig.ident")
DimPlot(cluster1.combined, reduction = "umap",label = TRUE, split.by = "orig.ident")
# tSNE (1:15)
cluster1.combined <- RunTSNE(cluster1.combined, reduction = "pca", dims = 1:15)
DimPlot(cluster1.combined, reduction = "tsne", group.by = "orig.ident")
DimPlot(cluster1.combined, reduction = "tsne")
DimPlot(cluster1.combined, reduction = "tsne", split.by = "orig.ident")
DimPlot(cluster1.combined, reduction = "tsne",label = TRUE, split.by = "orig.ident")
subclustering-6

5. 差異表達基因分析

結果全部確認無誤後,就能進行下一步的差異表達基因分析(DEGs, Differentially expression genes analysis)了。因為integrated assay主要是針對整合所做的分析,因此我的習慣是在做DEGs分析時,會先將assay轉到RNA中,然後對RNA assay再次進行標準化、選擇高度變異的特徵、以及進行數據縮放,然後再以FindAllMarkers()函數找出各組細胞群的差異表達基因。

# 差異表達基因分析
DefaultAssay(cluster1.combined) <- "RNA"
cluster1.combined <- NormalizeData(cluster1.combined, normalization.method = "LogNormalize", scale.factor = 10000)
cluster1.combined <- FindVariableFeatures(cluster1.combined, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(cluster1.combined)
cluster1.combined <- ScaleData(cluster1.combined, features = all.genes)
cluster1.combined.markers <- FindAllMarkers(cluster1.combined, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
View(cluster1.combined.markers)
top10 <- cluster1.combined.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(cluster1.combined, features = top10$gene, label = F)
library(ggplot2)
ggsave(path = "C:/Users/Administrator/Desktop", filename = "Heatmaptop.jpeg", width=40, height=20, dpi = 300, units = 'cm',limitsize = FALSE)
subclustering-7

結語

在本篇文章中,我們重點討論了多組單細胞數據集的分群(subclustering)方法和步驟。分析時的各項函數數值的設定,都取決於分析人員對自身數據的考量以及判斷,分析人員應該根據數據的特點、實驗設計和研究目的來選擇合適的函數和參數設置,避免盲目跟從教學參數。

參考文獻

1. Malte D Luecken and Fabian J Theis, Current best practices in single‐cell RNA‐seq analysis: a tutorial, Mol Syst Biol. 2019 Jun; 15(6): e8746.

2. Tim Stuart, Andrew Butler, Paul Hoffman, Christoph Hafemeister et al., Comprehensive integration of single-cell data, Cell. 2019 Jun 13;177(7):1888-1902.e21.

3. Seurat 官網: https://satijalab.org/seurat/articles/integration_introduction.html

您的分享,我非常感激!!!
MillionQuesn
MillionQuesn

一個旅居台灣的異鄉人,分享突然靈光一現的亮點。

文章: 46

發佈留言

發佈留言必須填寫的電子郵件地址不會公開。 必填欄位標示為 *

error: Alert: Content selection is disabled!!