本文參考Satija Lab的seurat 包，網址?https://satijalab.org/seurat

設置Seurat對象

在本教程中驶沼，我們將分析10X Genomics免費提供的外周血單核細胞（PBMC）數據集。在Illumina NextSeq 500上對2,700個單細胞進行了測序振亮。原始數據可以在這里找到挚歧。

我們首先閱讀數據摩梧。該Read10X函數從10X?讀取cellranger管道的輸出，返回唯一分子識別（UMI）計數矩陣卖漫。該矩陣中的值表示在每個細胞（列）中檢測到的每個特征（即基因;行）的分子數费尽。

接下來我們使用計數矩陣來創(chuàng)建一個Seurat對象。該對象充當容器羊始，其包含單細胞數據集的數據（如計數矩陣）和分析（如PCA或聚類結果）旱幼。有關Seurat對象結構的技術討論，請查看我們的GitHub Wiki突委。例如柏卤，計數矩陣存儲在pbmc[["RNA"]]@counts冬三。

library(dplyr)library(Seurat)

# Load the PBMC dataset

pbmc.data <- Read10X(data.dir ="../data/pbmc3k/filtered_gene_bc_matrices/hg19/")（路徑為下載的矩陣數據所在的位置，視各自下載路徑而定）

# Initialize the Seurat object with the raw (non-normalized data).

pbmc <- CreateSeuratObject(counts = pbmc.data, project ="pbmc3k", min.cells =3, min.features =200)

pbmc

## An object of class Seurat

## 13714 features across 2700 samples within 1 assay

## Active assay: RNA (13714 features)

計數矩陣中的數據是什么樣的缘缚？

# 讓我們檢查前30個細胞中的一些基因

pbmc.data[c("CD3D","TCL1A","MS4A1"),1:30]

## 3 x 30 sparse Matrix of class "dgCMatrix"

#### CD3D? 4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2? 3 . . . . . 3 4 1 5

## TCL1A . .? . . . . . . 1 . . . . . . . . . . .? . 1 . . . . . . . .

## MS4A1 . 6? . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .

矩陣中的值表示0（未檢測到分子）勾笆。由于scRNA-seq矩陣中的大多數值都是0，因此Seurat盡可能使用稀疏矩陣表示桥滨。這為Drop-seq / inDrop / 10x數據帶來了顯著的內存和速度節(jié)省窝爪。

dense.size <- object.size(as.matrix(pbmc.data))dense.size

## 709548272 bytes

sparse.size <- object.size(pbmc.data)sparse.size

## 29861992 bytes

dense.size/sparse.size

## 23.8 bytes

標準的預處理工作流程

以下步驟包括Seurat中scRNA-seq數據的標準預處理工作流程。這些代表了基于QC指標齐媒，數據標準化和縮放以及高度可變特征的檢測的細胞選擇和過濾蒲每。

QC并選擇細胞進行進一步分析

Seurat允許您根據任何用戶定義的標準輕松探索QC指標和過濾單元格。社區(qū)常用的一些QC指標包括

在每個細胞中檢測到的獨特基因的數量喻括。

低質量細胞或空液滴通常只有很少的基因

細胞雙峰或多重峰可能表現出異常高的基因計數

同樣邀杏，細胞內檢測到的分子總數（與獨特基因強烈相關）

映射到線粒體基因組的讀數百分比

低質量/垂死細胞通常表現出廣泛的線粒體污染

我們使用PercentageFeatureSet函數計算線粒體QC指標，該函數計算源自一組特征的計數百分比

我們使用所有基因MT-集作為一組線粒體基因

# The [[ operator can add columns to object metadata. This is a great place to stash QC stats

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern ="^MT-"

在下面的示例中唬血，我們可視化QC指標淮阐，并使用它們來過濾單元格。

我們過濾具有超過2,500或小于200的獨特特征計數的細胞

我們過濾線粒體計數> 5％的細胞

# Visualize QC metrics as a violin plot

VlnPlot(pbmc, features =c("nFeature_RNA","nCount_RNA","percent.mt"), ncol =3)

# FeatureScatter is typically used to visualize feature-feature relationships, but can be used# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.

plot1 <- FeatureScatter(pbmc, feature1 ="nCount_RNA", feature2 ="percent.mt")

plot2 <- FeatureScatter(pbmc, feature1 ="nCount_RNA", feature2 ="nFeature_RNA")

CombinePlots(plots = list(plot1, plot2))

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

規(guī)范化數據

從數據集中刪除不需要的單元格后刁品，下一步是規(guī)范化數據泣特。默認情況下，我們采用全局縮放規(guī)范化方法“LogNormalize”挑随，通過總表達式對每個單元格的要素表達式度量進行標準化状您，將其乘以比例因子（默認為10,000），并對結果進行對數轉換兜挨。歸一化值存儲在pbmc[["RNA"]]@data膏孟。

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

為清楚起見，在前一行代碼（以及將來的命令）中拌汇，我們?yōu)楹瘮嫡{用中的某些參數提供默認值柒桑。但是，這不是必需的噪舀，可以通過以下方式實現相同的行為：

pbmc <- NormalizeData(pbmc)

識別高度可變的 features （ FindVariableFeatures ）

接下來魁淳，我們計算在數據集中表現出高細胞間差異的特征子集（即，它們在一些細胞中高度表達与倡，而在其他細胞中低表達）界逛。我們和其他人已經發(fā)現，在下游分析中關注這些基因有助于突出單細胞數據集中的生物信號纺座。

我們在這里詳細描述了Seurat3中的程序息拜，并通過直接建模單細胞數據中固有的均值 - 方差關系對先前版本進行了改進，并在FindVariableFeatures函數中實現。默認情況下少欺，我們?yōu)槊總€數據集返回2,000個要素喳瓣。這些將用于下游分析，如PCA赞别。

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

# Identify the 10 most highly variable genestop10 <- head(VariableFeatures(pbmc),10)

# plot variable features with and without labels

plot1 <- VariableFeaturePlot(pbmc)

plot2 <- LabelPoints(plot = plot1, points = top10, repel =TRUE)

CombinePlots(plots = list(plot1, plot2))

縮放數據

接下來畏陕，我們應用線性變換（'縮放'），這是在PCA等降維技術之前的標準預處理步驟氯庆。該ScaleData函數可以：

1.改變每個基因的表達蹭秋，使得跨細胞的平均表達為0

2.縮放每個基因的表達，以便跨細胞的方差為1

3.該步驟在下游分析中給予相同的權重堤撵，因此高表達的基因不占優(yōu)勢

4.結果存儲在?pbmc[["RNA"]]@scale.data

all.genes <- rownames(pbmc)

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

這個步驟需要太長時間仁讨！我可以加快速度嗎？

縮放是Seurat工作流程中必不可少的一步实昨，但僅限于將用作PCA輸入的基因洞豁。因此，默認值ScaleData僅為對先前標識的變量要素執(zhí)行縮放（默認為2,000）荒给。為此丈挟，省略features上一個函數調用中的參數，即

pbmc <- ScaleData(pbmc)

您的PCA和群集結果將不受影響志电。然而曙咽，Seurat熱圖（如下圖所示產生DoHeatmap）需要對熱圖中的基因進行縮放，以確保高表達的基因不會影響熱圖挑辆。為了確保我們以后不將任何基因遺留在熱圖之外例朱，我們正在擴展本教程中的所有基因。

?如何在Seurat v2中刪除不需要的變化來源鱼蝉？

在Seurat v2我們還利用ScaleData函數從單細胞集刪除變化的干擾源洒嗤。例如，我們可以“消退”與（例如）細胞周期階段或線粒體污染相關的異質性魁亦。這些功能仍然支持ScaleData中Seurat v3渔隶，即：

pbmc <- ScaleData(pbmc, vars.to.regress ="percent.mt")

但是，特別是對于想要使用此功能的高級用戶洁奈，我們強烈建議使用新的規(guī)范化工作流程sctransform间唉。該方法在我們最近的預印本中進行了描述，這里使用Seurat v3單獨設置了一個小插圖睬魂。與此同時ScaleData终吼，該功能SCTransform還包括一個vars.to.regress參數。

執(zhí)行線性尺寸縮減

接下來氯哮，我們對縮放數據執(zhí)行PCA。默認情況下，只有先前確定的變量要素用作輸入喉钢，但features如果要選擇其他子集姆打，則可以使用參數定義。

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

提供定義可視化細胞和功能的幾種有用的方式PCA肠虽，包括VizDimReduction幔戏，DimPlot，和DimHeatmap

# Examine and visualize PCA results a few different ways

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")

DimPlot(pbmc, reduction ="pca")

特別是DimHeatmap允許容易地探索數據集中的異質性的主要來源税课，并且在試圖決定包括哪些PC以用于進一步的下游分析時可以是有用的闲延。細胞和特征都根據其PCA分數排序。設置cells為數字可以繪制光譜兩端的“極端”單元格韩玩，從而大大加快了大型數據集的繪圖速度垒玲。雖然顯然是一種監(jiān)督分析，但我們發(fā)現這是探索相關特征集的有用工具找颓。

DimHeatmap(pbmc, dims =1, cells =500, balanced =TRUE)

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

確定數據集的“維度”

為了克服scRNA-seq數據的任何單一特征中的廣泛技術噪聲合愈，Seurat基于其PCA分數對細胞進行聚類，每個PC基本上代表在相關特征集中組合信息的“元特征”击狮。因此佛析，頂部主要組件表示數據集的強大壓縮。但是彪蓬，我們應該選擇包含多少個組件寸莫？10？20档冬？100膘茎？

在Macosko?等人中，我們實施了一項受JackStraw程序啟發(fā)的重采樣測試捣郊。我們隨機置換數據的子集（默認為1％）并重新運行PCA辽狈，構建特征分數的“空分布”，并重復此過程呛牲。我們將“重要”PC識別為具有低p值特征的強大豐富的PC刮萌。

#NOTE:This process can take a long time for big datasets, comment out for expediency. More approximate techniques such as those implemented in ElbowPlot() can be used to reduce computation time

pbmc <- JackStraw(pbmc, num.replicate =100)

pbmc <- ScoreJackStraw(pbmc, dims =1:20)

該JackStrawPlot功能提供了一種可視化工具，用于比較每個PC的p值分布和均勻分布（虛線）娘扩∽湃祝“重要的”PC將顯示具有低p值的特征的強烈豐富（虛線上方的實線）。在這種情況下琐旁，看起來在前10-12個PC之后涮阔，重要性急劇下降。

JackStrawPlot(pbmc, dims =1:15)

另一種啟發(fā)式方法生成“肘圖”：基于每個（ElbowPlot函數）解釋的方差百分比對主成分進行排序灰殴。在這個例子中敬特，我們可以觀察PC9-10周圍的“肘部”，這表明大多數真實信號都是在前10PC中捕獲的。

ElbowPlot(pbmc)

識別數據集的真實維度 - 對用戶來說可能具有挑戰(zhàn)性/不確定性伟阔。因此辣之，我們建議考慮這三種方法。第一個是更多的監(jiān)督皱炉，探索PC以確定異質性的相關來源怀估，并且可以與GSEA一起使用。第二個實現基于隨機空模型的統(tǒng)計測試合搅，但對于大型數據集來說是耗時的多搀，并且可能不會返回明確的PC截止。第三種是常用的啟發(fā)式算法灾部，可以立即計算康铭。在這個例子中，所有三種方法產生了類似的結果梳猪，但我們可能有理由在PC 7-12之間選擇任何作為截止值的東西麻削。

我們在這里選擇了10，但鼓勵用戶考慮以下內容：

1.樹突細胞和NK研究者可以認識到與PC12和13強烈相關的基因定義了罕見的免疫亞群（即MZB1是漿細胞樣DC的標記）春弥。然而呛哟，這些組是如此罕見，如果沒有先驗知識匿沛，很難將這種大小的數據集與背景噪聲區(qū)分開來扫责。

2.我們鼓勵用戶使用不同數量的PC（10,15甚至50！）重復下游分析逃呼。正如您將觀察到的鳖孤，結果通常沒有顯著差異。

3.在選擇此參數時抡笼，我們建議用戶在較高的一側犯錯苏揣。例如，僅使用5臺PC進行下游分析確實會對結果造成嚴重影響推姻。

聚類細胞

Seurat v3采用基于圖形的聚類方法平匈，建立在（Macosko?等人）的初始策略之上。重要的是藏古，驅動聚類分析的距離度量（基于先前識別的PC）保持不變增炭。然而，我們將細胞距離矩陣分成簇的方法已經大大改進拧晕。我們的方法受到最近手稿的啟發(fā)隙姿，這些手稿將基于圖形的聚類方法應用于scRNA-seq數據[SNN-Cliq，Xu和Su厂捞，Bioinformatics输玷，2015]和CyTOF數據[PhenoGraph队丝，Levine?等，Cell饲嗽，2015]炭玫。簡而言之奈嘿，這些方法將單元格嵌入圖形結構中 - 例如K-最近鄰（KNN）圖形貌虾，在具有相似特征表達模式的單元格之間繪制邊緣，然后嘗試將此圖形劃分為高度互連的“準集團”或“社區(qū)”裙犹。

與PhenoGraph一樣尽狠，我們首先根據PCA空間中的歐氏距離構建KNN圖，并根據其局部鄰域中的共享重疊（Jaccard相似性）細化任意兩個單元之間的邊權重叶圃。使用該FindNeighbors函數執(zhí)行該步驟袄膏，并將先前定義的數據集維度（前10個PC）作為輸入。

為了聚類細胞掺冠，我們接下來應用模塊化優(yōu)化技術沉馆，例如Louvain算法（默認）或SLM?[SLM，Blondel?等德崭，Journal of Statistical Mechanics]斥黑，將細胞迭代分組在一起，目的是優(yōu)化標準模塊化功能眉厨。該FindClusters函數實現此過程锌奴，并包含一個分辨率參數，用于設置下游群集的“粒度”憾股，增加的值會導致更多的群集鹿蜀。我們發(fā)現將此參數設置在0.4-1.2之間通常會返回大約3K單元格的單細胞數據集的良好結果。較大的數據集通常會增加最佳分辨率服球≤钋。可以使用該Idents功能找到群集。

pbmc <- FindNeighbors(pbmc, dims =1:10)

pbmc <- FindClusters(pbmc, resolution =0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

#### Number of nodes: 2638

## Number of edges: 96033

#### Running Louvain algorithm...

## Maximum modularity in 10 random starts: 0.8720

## Number of communities: 9

## Elapsed time: 0 seconds

# Look at cluster IDs of the first 5 cellshead(Idents(pbmc),5)

## AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG AAACCGTGTATGCG

##? ? ? ? ? ? ? 1? ? ? ? ? ? ? 3? ? ? ? ? ? ? 1? ? ? ? ? ? ? 2? ? ? ? ? ? ? 6

## Levels: 0 1 2 3 4 5 6 7 8

運行非線性降維（UMAP / tSNE）

Seurat提供了幾種非線性降維技術斩熊，如tSNE和UMAP往枣，可視化和探索這些數據集。這些算法的目標是學習數據的基礎流形座享，以便在低維空間中將相似的單元放在一起蛋济。上面確定的基于圖的聚類內的單元應該共同定位這些降維圖。作為UMAP和tSNE的輸入华望，我們建議使用相同的PC作為聚類分析的輸入希停。

# If you haven't installed UMAP, you can do so via reticulate::py_install(packages =# 'umap-learn')

（V3 版本這里有一個bug，裝了umap-learn以后還是會報錯 cant find umap淳衙，解決方法為：

1.restart R

2.重新run一次或兩次下面的代碼）

pbmc <- RunUMAP(pbmc, dims =1:10)

# note that you can set `label = TRUE` or use the LabelClusters function to help label individual clusters

DimPlot(pbmc, reduction ="umap")

您可以在此時保存對象蘑秽，以便可以輕松地將其加載回來饺著，而無需重新運行上面執(zhí)行的計算密集型步驟，或者可以輕松地與協(xié)作者共享肠牲。

saveRDS(pbmc, file ="../output/pbmc_tutorial.rds")

尋找差異表達的特征（聚類生物標志物）

Seurat可以幫助您找到通過差異表達定義聚類的標記幼衰。默認情況下，它標識單個群集（指定ident.1）的正負標記缀雳，與所有其他單元格進行比較渡嚣。FindAllMarkers為所有群集自動執(zhí)行此過程，但您也可以測試群集彼此之間或所有單元格肥印。

該min.pct論證要求在兩組單元中的任一組中以最小百分比檢測特征识椰，并且thresh.test參數要求在兩組之間以一定量差異地表示（平均）特征。您可以將這兩個設置為0深碱，但時間會急劇增加 - 因為這將測試大量不太可能具有高度不同的功能腹鹉。作為加速這些計算的另一種選擇，max.cells.per.ident可以設置敷硅。這將對每個標識類進行下采樣功咒，使其不再具有比設置的更多的單元格。雖然通常會出現功率損失绞蹦，但速度增加可能是顯著的力奋，并且最高度差異表達的功能可能仍然會升至頂部。

# find all markers of cluster 1

cluster1.markers <- FindMarkers(pbmc, ident.1 =1, min.pct =0.25)

head(cluster1.markers, n =5)

# find all markers distinguishing cluster 5 from clusters 0 and 3

cluster5.markers <- FindMarkers(pbmc, ident.1 =5, ident.2 = c(0,3), min.pct =0.25)

head(cluster5.markers, n =5)

# find markers for every cluster compared to all remaining cells, report only the positive ones

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_logFC)

Seurat有幾個差分表達式測試坦辟，可以使用test.use參數設置（有關詳細信息刊侯，請參閱我們的DE插圖）。例如锉走，ROC測試返回任何單個標記的“分類能力”（范圍從0 - 隨機滨彻，到1 - 完美）。

cluster1.markers <- FindMarkers(pbmc, ident.1 =0, logfc.threshold =0.25, test.use ="roc", only.pos =TRUE)

我們提供了幾種用于可視化標記表達的工具挪蹭。VlnPlot（顯示跨群集的表達概率分布）亭饵，以及FeaturePlot（在tSNE或PCA圖上可視化特征表達）是我們最常用的可視化。我們還建議您探索RidgePlot梁厉，CellScatter以及DotPlot查看數據集的其他方法辜羊。

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

# you can plot raw counts as well

VlnPlot(pbmc, features = c("NKG7","PF4"), slot ="counts", log =TRUE)

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

DoHeatmap為給定的單元格和要素生成表達式熱圖。在這種情況下词顾，我們正在為每個群集繪制前20個標記（或所有標記八秃，如果小于20）。

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n =10, wt = avg_logFC)

DoHeatmap(pbmc, features = top10$gene) + NoLegend()

將單元類型標識分配給集群

幸運的是肉盹，在這個數據集的情況下昔驱，我們可以使用規(guī)范標記輕松地將無偏的聚類與已知的細胞類型相匹配：

new.cluster.ids <- c("Naive CD4 T","Memory CD4 T","CD14+ Mono","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()

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