Seurat.jpg

后臺有讀者翻到了一年前發(fā)的文獻(xiàn)解讀，請教了一下文章的圖的做法。正好前段時間剛做過單細(xì)胞轉(zhuǎn)錄組分析，今天就給大家介紹一下常用工具Seurat的用法神汹。

Seurat 4.0 使用指南

設(shè)置Seurat對象

示例數(shù)據(jù)

10X Genomics免費提供的外周血單核細(xì)胞(PBMC)數(shù)據(jù)集。通過Illumina NextSeq 500測序的2700個單細(xì)胞古今。示例數(shù)據(jù)下載:https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz屁魏。(或者后臺回復(fù)seurat,領(lǐng)取完整代碼及示例數(shù)據(jù))

setwd(".../filtered_gene_bc_matrices/hg19") #設(shè)置工作環(huán)境到數(shù)據(jù)所在文件夾
#安裝和加載所需包
BiocManager::install("Seurat") 
BiocManager::install("dplyr")
BiocManager::install("patchwork")
library(dplyr)
library(Seurat)
library(patchwork) 

#導(dǎo)入示例數(shù)據(jù)
pbmc.data <- Read10X(data.dir = ".../filtered_gene_bc_matrices/hg19/")#自行填寫數(shù)據(jù)所在文件夾
#創(chuàng)建Seurat對象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
#過濾檢測少于200個基因的細(xì)胞（min.features = 200）和少于3個細(xì)胞檢測出的基因（min.cells = 3）
pbmc

#參數(shù)解釋
CreateSeuratObject(
  counts, #未標(biāo)準(zhǔn)化的數(shù)據(jù)，如原始計數(shù)或TPMs
  project = "CreateSeuratObject",#設(shè)置Seurat對象的項目名稱
  assay = "RNA", #與初始輸入數(shù)據(jù)對應(yīng)的分析名稱
  names.field = 1,#對于每個cell的初始標(biāo)識類捉腥，從cell的名稱中選擇此字段氓拼。例如，如果cell在輸入矩陣                       #中被命名為BARCODE_CLUSTER_CELLTYPE抵碟，則設(shè)置名稱披诗。字段設(shè)置為3以將初始標(biāo)識設(shè)置為                    #CELLTYPE。
  names.delim = "_", #對于每個cell的初始標(biāo)識類立磁，從cell的列名中選擇此分隔符。例如剥槐，如果cell命名                         #為bar - cluster - celltype唱歧，則將此設(shè)置為“-”，以便將cell名稱分離到其組成部分                     #中粒竖，以選擇相關(guān)字段颅崩。
  meta.data = NULL, #要添加到Seurat對象的其他單元級元數(shù)據(jù)。應(yīng)該是data.frame蕊苗，其中行是單元格名稱沿后，列                    #是附加的元數(shù)據(jù)字段。
  ...
  min.cells #包含至少在這些細(xì)胞檢測到的features朽砰。
  min.features #包含至少檢測到這些features的細(xì)胞
)

> pbmc
An object of class Seurat 
13714 features across 2700 samples within 1 assay 
Active assay: RNA (13714 features, 0 variable features)
#1個數(shù)據(jù)集尖滚，包含2700個細(xì)胞喉刘，13714個基因。

數(shù)據(jù)矩陣

# 查看這三個基因的前三十行矩陣
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
> pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
3 x 30 sparse Matrix of class "dgCMatrix"
   [[ suppressing 30 column names ‘AAACATACAACCAC-1’, ‘AAACATTGAGCTAC-1’, ‘AAACATTGATCAGC-1’ ... ]]
                                              
CD3D  4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2
TCL1A . .  . . . . . . 1 . . . . . . . . . . .
MS4A1 . 6  . . . . . . 1 1 1 . . . . . . . . .
                          
CD3D   3 . . . . . 3 4 1 5
TCL1A  . 1 . . . . . . . .
MS4A1 36 1 2 . . 2 . . . .

.在矩陣中的值表示0(未檢測到分子)漆弄。由于scRNA-seq矩陣中的大多數(shù)值為0,因此Seurat在任何可能的情況下都使用稀疏矩陣表示睦裳。這為Drop-seq/inDrop/10x數(shù)據(jù)節(jié)省了大量內(nèi)存和速度。

標(biāo)準(zhǔn)的預(yù)處理流程

下面的步驟包含了Seurat中的scRNA-seq數(shù)據(jù)的標(biāo)準(zhǔn)預(yù)處理流程撼唾。包括QC（質(zhì)控）廉邑、數(shù)據(jù)歸一化以及細(xì)胞的選擇和過濾。

QC和細(xì)胞篩選

常用的質(zhì)控指標(biāo)：

• 每個細(xì)胞在檢測到的特異基因數(shù)
  • 低質(zhì)量細(xì)胞或空液滴通常只能檢測到非常少的基因
  • 兩個或多個細(xì)胞被同時捕獲通常會有很高的基因數(shù)
• 每個細(xì)胞檢測到的分子總數(shù)（與基因密切相關(guān)）
• 每個細(xì)胞的線粒體基因比例
  • 低質(zhì)量/瀕死細(xì)胞常表現(xiàn)出廣泛的線粒體污染
  • 使用PercentageFeatureSet()函數(shù)計算線粒體QC指標(biāo)
  • 使用所有以MT-開頭的基因作為一組線粒體基因

#向pbmc新增一列percent.mt數(shù)據(jù)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

QC指標(biāo)儲存在哪倒谷？

每個細(xì)胞基因數(shù)和總分子數(shù)在建立seurat對象時就已經(jīng)自動計算好了蛛蒙。

#展示前5個細(xì)胞的QC指標(biāo)
head(pbmc@meta.data, 5)
> head(pbmc@meta.data, 5)
                 orig.ident nCount_RNA
AAACATACAACCAC-1     pbmc3k       2419
AAACATTGAGCTAC-1     pbmc3k       4903
AAACATTGATCAGC-1     pbmc3k       3147
AAACCGTGCTTCCG-1     pbmc3k       2639
AAACCGTGTATGCG-1     pbmc3k        980
                 nFeature_RNA percent.mt
AAACATACAACCAC-1          779  3.0177759
AAACATTGAGCTAC-1         1352  3.7935958
AAACATTGATCAGC-1         1129  0.8897363
AAACCGTGCTTCCG-1          960  1.7430845
AAACCGTGTATGCG-1          521  1.2244898

#使用小提琴圖可視化QC指標(biāo)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

• nFeature_RNA代表每個細(xì)胞測到的基因數(shù)目。
• nCount_RNA代表每個細(xì)胞測到所有基因的表達(dá)量之和渤愁。
• percent.mt代表測到的線粒體基因的比例牵祟。

QC

#FeatureScatter通常用于可視化 feature-feature 相關(guān)性，
#nCount_RNA 與percent.mt的相關(guān)性
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
#nCount_RNA與nFeature_RNA的相關(guān)性
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2 #合并兩圖

SCATTER

過濾線粒體基因表達(dá)比例過高的細(xì)胞猴伶，和一些極值細(xì)胞（可以根據(jù)小提琴圖判斷课舍，查看兩端離群值）。

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
#濾掉 2500 > nFeature_RNA >200 和percent.mt < 5的數(shù)據(jù)

數(shù)據(jù)標(biāo)準(zhǔn)化

默認(rèn)情況下他挎，使用全局縮放歸一化方法“LogNormalize”筝尾，用總表達(dá)量對每個細(xì)胞的基因表達(dá)式進(jìn)行歸一化，再乘以一個縮放因子(默認(rèn)為10,000)办桨，然后對結(jié)果進(jìn)行log轉(zhuǎn)換筹淫。標(biāo)準(zhǔn)化的數(shù)值存儲在pbmc[["RNA"]]@data中。

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

若所有調(diào)用的參數(shù)都是默認(rèn)值呢撞，則可省去损姜。

pbmc <- NormalizeData(pbmc)

鑒定高變基因

接下來，計算數(shù)據(jù)集中表現(xiàn)出高細(xì)胞間變異的特征基因(即殊霞，它們在某些細(xì)胞中高表達(dá)摧阅，而在其他細(xì)胞中低表達(dá))。這些基因有助于突出單細(xì)胞數(shù)據(jù)集中的生物信號绷蹲。

用FindVariableFeatures()函數(shù)實現(xiàn)棒卷。默認(rèn)情況下，每個數(shù)據(jù)集返回2000個features 祝钢。這些將用于下游分析比规，如PCA。

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

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

# 畫出不帶標(biāo)簽或帶標(biāo)簽基因點圖
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

variabel

數(shù)據(jù)縮放

線性變換(“縮放”)拦英，是在PCA降維之前的一個標(biāo)準(zhǔn)預(yù)處理步驟蜒什。ScaleData()函數(shù)功能:

• 轉(zhuǎn)換每個基因的表達(dá)值，使每個細(xì)胞的平均表達(dá)為0

• 轉(zhuǎn)換每個基因的表達(dá)值疤估，使細(xì)胞間的方差為1

  • 此步驟在下游分析中具有相同的權(quán)重灾常，因此高表達(dá)的基因不會占主導(dǎo)地位

• 結(jié)果存儲在pbmc[["RNA"]]@scale.data中

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

線性降維

接下來霎冯，對縮放的數(shù)據(jù)執(zhí)行PCA。默認(rèn)情況下岗憋，只使用前面確定的變量特性作為輸入肃晚，但是如果想選擇不同的子集，可以使用features參數(shù)來定義仔戈。

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

Seurat提供了幾種有用的方法來可視化細(xì)胞和定義PCA的特性关串，包括VizDimReduction、DimPlot和DimHeatmap

#查看PCA結(jié)果
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
> print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
PC_ 1 
Positive:  CST3, TYROBP, LST1, AIF1, FTL 
Negative:  MALAT1, LTB, IL32, IL7R, CD2 
PC_ 2 
Positive:  CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1 
Negative:  NKG7, PRF1, CST7, GZMB, GZMA 
PC_ 3 
Positive:  HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1 
Negative:  PPBP, PF4, SDPR, SPARC, GNG11 
PC_ 4 
Positive:  HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1 
Negative:  VIM, IL7R, S100A6, IL32, S100A8 
PC_ 5 
Positive:  GZMB, NKG7, S100A8, FGFBP2, GNLY 
Negative:  LTB, IL7R, CKB, VIM, MS4A7 

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

VizDimLoadings

DimPlot(pbmc, reduction = "pca")

DimPlot

DimHeatmap()可以方便地探索數(shù)據(jù)集中異質(zhì)性的主要來源监徘，并且可以確定哪些PC維度可以用于下一步的下游分析晋修。細(xì)胞和基因根據(jù)PCA分?jǐn)?shù)來排序。

DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE) #1個PC 500個細(xì)胞

pc1

DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE) #15個PC

pc1-15

確定數(shù)據(jù)的維度

主成分分析的原理非常簡單凰盔，概括來說就是選擇包含信息量大的維度（features）墓卦，去除信息量少的“干擾”維度。所以這里會有個問題——如何知道保留幾個維度是最佳的呢户敬？我們希望通過保留盡可能少的維度來留存盡可能多的信息落剪。Seurat有兩種方法來確定維度。

• JackStraw

pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)

JackStraw

可以看出尿庐，在10-12個PC之后忠怖，顯著性大幅下降，也就是前10-12個維度包含了大部分的樣本信息抄瑟。

• Elbow plot

ElbowPlot(pbmc)

Elbow

可以看出凡泣，PC9-10附近有一個拐點（“elbow”），這表明大部分真實信號是在前10個pc中捕獲的皮假。

綜合以上方法鞋拟，選擇10個主成成分作為參數(shù)用于后續(xù)分析。

細(xì)胞聚類

Seurat使用KNN算法進(jìn)行聚類惹资。

pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
#dims = 1:10 即選取前10個主成分來分類細(xì)胞贺纲。
#查看前5個細(xì)胞的分類ID
head(Idents(pbmc), 5)

非線性降維（UMAP/tSNE）

pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")
# 顯示在聚類標(biāo)簽
DimPlot(pbmc, reduction = "umap", label = TRUE)
# 使用TSNE聚類
pbmc <- RunTSNE(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "tsne")
# 顯示在聚類標(biāo)簽
DimPlot(pbmc, reduction = "tsne", label = TRUE)

#保存rds，用于后續(xù)分析
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")

找差異表達(dá)基因(聚類標(biāo)志cluster biomarkers)

利用 FindMarkers 命令褪测，可以找到找到各個細(xì)胞類型中與其他類別的差異表達(dá)基因猴誊，作為該細(xì)胞類型的生物學(xué)標(biāo)記基因。

• dent.1參數(shù)設(shè)置待分析的細(xì)胞類別

• min.pct參數(shù)汰扭，在兩組細(xì)胞中的任何一組中檢測到的最小百分

• thresh.test參數(shù)，在兩組細(xì)胞間以一定數(shù)量的差異表達(dá)（平均）

• max.cells.per.ident參數(shù)福铅，通過降低每個類的采樣值萝毛，提高計算速度

# cluster 1的標(biāo)記基因
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
#找出區(qū)分cluster 5與cluster 0和cluster 3的所有標(biāo)記
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
# 找出每個cluster的標(biāo)記與所有剩余的細(xì)胞相比較，只報告陽性細(xì)胞
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)

可視化

VlnPlot(pbmc, features = c("MS4A1", "CD79A"))

violin

# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)

violin_2

FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", 
    "CD8A"))

feauture

#每個聚類前10個差異基因表達(dá)熱圖(如果小于10滑黔，則繪制所有標(biāo)記)
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()

heat

鑒定細(xì)胞類型

這個數(shù)據(jù)集的markers與已知細(xì)胞的marker可以輕松配對笆包。也可以通過查閱相關(guān)文獻(xiàn)人工注釋环揽，或者利用singleR（挖個坑，有空再來填）自動注釋庵佣。

#加上注釋
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()

umap

#保存
saveRDS(pbmc, file = "../output/pbmc3k_final.rds")

往期內(nèi)容：

R繪圖實戰(zhàn)|GSEA富集分析圖