Seurat V4.9.9 – 一個強大的單細胞分析R套件

前言

Seurat是一個針對單細胞數據分析的R套件,最初由Satija Lab於2015年開發,是首批針對單細胞RNA-seq數據進行整合、分析和可視化的R套件之一,目前已經開發到V5 beta 版本。

如何下載Seurat套件?

首先,你必須在電腦上安裝R語言和RStudio,然後打開RStudio視窗後輸入安裝Seurat套件的指令:

# 下載 remotes 套件 
install.packages("remotes")
# 下載Rtools4.3 https://cran.r-project.org/bin/windows/Rtools/
# 下載 Seurat
remotes::install_github("satijalab/seurat", "seurat5", quiet = TRUE)

如果要下載早先版本的Seurat套件,可用以下方法進行下載:

# 下載 remotes package 
install.packages('remotes')
# 範例 Replace '2.3.0' with your desired version
remotes::install_version(package = 'Seurat', version = package_version('2.3.0'))
library(Seurat)

範例數據

在這次的範例中,我們將分析來自10X Genomics的周邊血(PBMC)單細胞數據集,這是一個包含2700個單細胞的數據集,採用Illumina NextSeq 500進行测序,原始數據可以在這裡下載

1. 加載套件和載入PBMC數據集

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

# 檢查版本
packageVersion("Seurat")

# 載入PBMC數據集,需用電腦中的正確路徑,須將"\"更改為"/"
pbmc.data <- Read10X(data.dir = "C:/Users/Administrator/Desktop/hg19")
# 創建Seurat對象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc

在這裡使用CreateSeuratObject()函數導入數據,其中counts參數是輸入的數據,project參數是用於標識Seurat對象的名稱,min.cells和min.features參數是為了過濾掉一些低質量細胞和基因。接著我們便可以繼續進一步的分析,這裡使用的分析方式是標準分析法。

2. 過濾掉低品質的細胞

低質量或者瀕死的細胞通常具有廣泛的粒線體污染,因此我們可以使用PercentageFeatureSet()函數來計算粒線體的指標,並進一步將異常的細胞移除,以免影響後續資料分析。至於移除參數的cut off值,則取決於研究人員

# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

[[ ]]符號可將列添加到對象數據里,在這個範例中,PBMC生成Seurat對象後,通過運行上述的指令,能將粒線體基因組成的QC指標添加到了PBMC的數據中。如果是人類數據集,pattern 為”^MT-“;如果是老鼠數據集,pattern 則設為”^mt-“。

接著執行QC指標視覺化指令:

# 使用小提琴圖將QC指標視覺化
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
  1. nFeature_RNA:每顆細胞的基因數目,大部分細胞基因數在1000 左右。
  2. nCount_RNA:每顆細胞的總 UMI count (Unique Molecular Identifier,是基因表現的量化標記),大部分細胞的count落在 2300 左右。
  3. percent.mt:粒線體基因所占的比例 (粒線體基因總 UMI count ÷ 細胞的總 UMI count),大部分細胞的粒線體比例落在 5%以下。

最後,Seurat範例篩選的條件是nFeature_RNA介於 200至2500 之間、粒線體基因占比小於 5% 的細胞做後續分析。

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc

執行完篩選條件後,使用在後續分析的細胞還剩下2638顆細胞。

3. 標準化 (Normalization)

在單細胞RNA數據分析中,由於不同細胞之間的RNA分子計數差異和批次效應等因素,需要對原始數據進行標準化以便於後續的分析。在Seurat中,常用的標準化方法是對每個細胞的基因表達數據(UMI count)進行總數標準化,即將每個細胞的所有基因的計數值除以該細胞的總計數,並將得到的商值乘以一個縮放因數(通常選擇10000)。這個過程可以通過Seurat中的NormalizeData()函數實現。進行標準化後,所有細胞的基因表達數據都處於相同的量級,方便後續分析。

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
# 若所有參數都是默認值,則可以簡寫成以下指令
# pbmc <- NormalizeData(pbmc)

4. 選擇高度變異的特徵 (Highly Variable Features, HVGs)

在單細胞RNA數據分析中,每個基因的表達模式都可能不同。有些基因的表達量可能在細胞之間變化較小,而有些基因的表達量則可能會在不同細胞之間有較大的差異。對於後者,這些基因的表達差異可能與不同細胞類型或狀態的差異相關聯。因此,在單細胞RNA數據分析中,需要識別具有高度差異性的基因,並將其用於下游的聚類和差異表達分析。在Seurat中,可以通過FindVariableFeatures()函數來選擇HVGs,這個函數通過計算每個基因的平均表達和離散程度(variance-to-mean ratio)來鑒定HVGs。通過選擇這些HVGs,可以減少對不差異的基因進行計算,從而提高下游分析的效率和準度。

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# 查看變異最高的10個基因
top10 <- head(VariableFeatures(pbmc), 10)

# 畫出所標示的10個高變基因
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

5. 特徵縮放 (Scaling data)

特徵縮放是一個重要步驟,它用於調整細胞之間基因表達的差異大小,以便在後續分析中更好地捕捉細胞間的差異。特徵縮放是以z-score進行轉換,這可以使得每個基因的平均值為0且標準差為1。在Seurat中,用scale.data() 函數可以實現這個步驟。

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

6. 線性降維分析 (linear dimensional reduction)

線性降維(linear dimensional reduction)是指將高維度數據轉換為低維度數據的一種數據降維方法,通過將原始高維度特徵空間映射到低維度特徵空間,可以去除冗餘特徵,同時保留有用的信息。在單細胞RNA數據分析中,線性降維可以用來探索數據的主要變異方向,發現細胞之間的相似性和差異性。在Seurat中,我們可以使用 RunPCA() 函數(PCA, Principle Component Analysis)來進行PCA降維,可以將 verbose 參數設置為 TRUE 來查看更多有關PCA計算的信息,並且使用 VizDimReduction()、DimPlot()以及DimHeatmap()函數來可視化PCA的結果。

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

# 查看PCA的結果,這裡以DimHeatmap為例子
# VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
# DimPlot(pbmc, reduction = "pca")
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
# 叫出PC1至PC15的結果
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)

7. 確定數據的維度

Seurat基於其PCA分數對細胞進行聚類,其中每個PC基本上表示一個主成分,它結合了相關特徵集中的信息。接下來,我們應該選擇多少個PCA主成分來進行分析呢?Seurat套件中有兩種方法來確定PCA主成分的選擇,一種是JackStraw()函數,而另一種為Elbowplot()函數,在這裡我將會以Elbowplot()函數當作例子:

ElbowPlot(pbmc)

從圖中結果中可以看出,PC9-10附近有一個拐點(elbow),這表明大部分真實信號是在前10個PC中被捕獲,因此可以選擇10個PCA主成分作為參數用於後續分析。

8. 細胞聚集 (Cell Clustering)

細胞聚集是指將單細胞RNA-seq數據中相似的細胞群聚在一起的分析過程。通過細胞聚集,可以將大量單細胞數據分為不同的亞群,從而更好地理解細胞的異質性。

pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
# 查看前五個細胞的ID
head(Idents(pbmc), 5)

dims = 1:10 意思是選取前十個主成分來進行細胞分群,千萬千萬記得加 “:” 號,這才表示選擇了1-10個PCA!!!

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

非線性降維分析法是一種將高維度資料轉換為低維度資料的方法,且盡可能保留資料的局部結構和相對距離。在單細胞RNA數據中,UMAP(Uniform Manifold Approximation and Projection)和tSNE(t-Distributed Stochastic Neighbor Embedding)是常用的非線性降維分析法,用於將數據映射到二維或三維空間,以便於視覺化和理解數據間的關係。

# 使用UMAP聚類
pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")
# 顯示cluster標記
DimPlot(pbmc, reduction = "umap", label = TRUE)
# 使用TSNE聚類
pbmc <- RunTSNE(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "tsne")
# 顯示cluster標記
DimPlot(pbmc, reduction = "tsne", label = TRUE)
圖例為UMAP結果圖

10. 差異表達基因分析

差異表達基因分析(DEGs, Differentially expression genes analysis)是指將細胞分成不同群集(cluster)後,比較這些群集中基因的表現量,以此來確定哪些基因在不同的群集中有顯著的差異表達。

# 尋找cluster 2中所有的特徵基因
cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n = 5)
# 找出區分cluster 5與cluster 0和cluster 3的所有特徵基因
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
# 找出所有cluster的特徵基因,並且只報告陽性特徵基因而已(only positive markers)
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)

11. 可視化結果

VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
# 用counts畫圖並且取log值
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
# 熱圖(heatmap)
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()

12. 鑑定每個群集的細胞類型

如果有判別各類細胞的特定基因,那麼接下來就可以直接執行鑑定細胞類型的步驟。不同細胞型態的特定基因,可以使用資料庫(database)、文獻查詢或者R套件來協助細胞類型的註解。整個單細胞數據分析過程中,就屬這一步最消耗時間,因為一旦判別錯誤,基本上就等於要重來一遍。

當判別完細胞類型後,就可以標註細胞的名稱了。

new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono",
 "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

13. 保存

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

結語

Seurat是一款非常強大的單細胞分析R套件,可以說是單細胞分析中的大哥也不為過。但是在分析的過程中,了解各項參數的背景、原理和設定考量是研究人員必須要做的功課,這樣才不會分析到最後才驚覺自己做了個寂寞

參考文獻

1. Butler, A., Hoffman, P., Smibert, P., Papalexi, E., & Satija, R. (2018). Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nature Biotechnology, 36(5), 411–420.

2. Seurat 官網 https://satijalab.org/seurat/index.html

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

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

文章: 46

發佈留言

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

error: Alert: Content selection is disabled!!