Seurat簡(jiǎn)介

??Seurat---幾乎是當(dāng)前單細(xì)胞RNA-seq分析領(lǐng)域的不可或缺的工具，特別是基于10X公司的cellrange流程得出的結(jié)果究孕，可以方便的對(duì)接到Seurat工具中進(jìn)行后續(xù)處理衩椒，簡(jiǎn)直是帶給迷茫在單細(xì)胞數(shù)據(jù)荒漠中小白的一眼清泉隆圆，相對(duì)全面的功能楣责，簡(jiǎn)潔的操作命令萌狂，如絲般順滑潭千。

seurat官網(wǎng)界面

更新記錄

??Seurat目前最新版為V3穿扳，第一版是為處理空間轉(zhuǎn)錄組而設(shè)計(jì)衩侥，第二版針對(duì)基因droplet的單細(xì)胞技術(shù)而開發(fā)的用于單細(xì)胞質(zhì)控，降維矛物，聚類茫死，鑒定細(xì)胞類型marker等多重功能，目前更新的第三版可以整合多組學(xué)履羞，多批次璧榄，多實(shí)驗(yàn)方法的數(shù)據(jù)，把握了后續(xù)數(shù)據(jù)分析需要多組學(xué)整合的熱點(diǎn)策略吧雹，符合下一步分析思路更全面骨杂，生物學(xué)解釋更透徹的科研需求。
以下為Seurat更新的記錄：
2019年4月16日發(fā)布3.0版
2018年11月2日版本3.0 alpha發(fā)布

• 發(fā)布的預(yù)印本描述了在單細(xì)胞數(shù)據(jù)集中識(shí)別“錨點(diǎn)”的新方法
  2018年3月23日 2.3版本發(fā)布
• 提高速度和內(nèi)存效率
• 用于分析來自Microwell-seq Mouse Cell Atlas數(shù)據(jù)集的~250,000個(gè)細(xì)胞的新視圖
  2018年1月10日： 2.2版本發(fā)布
• 支持多數(shù)據(jù)集對(duì)齊
• 評(píng)估對(duì)齊性能的新方法
  2017年10月16日： 2.1版發(fā)布
• 支持多模式單細(xì)胞數(shù)據(jù)
• 支持用于差異表達(dá)式測(cè)試的MAST和DESeq2包
  2017年7月26日： 2.0版本發(fā)布
• 發(fā)布預(yù)印本用于scRNA-seq數(shù)據(jù)集的綜合分析
• 數(shù)據(jù)集集成雄卷，可視化和探索的新方法
• 對(duì)代碼庫(kù)進(jìn)行重大重組搓蚪，以強(qiáng)調(diào)清晰度和清晰的文檔
  2016年10月4日： 1.4版本發(fā)布
• 為UMI計(jì)數(shù)數(shù)據(jù)添加了負(fù)二項(xiàng)回歸和差分表達(dá)式測(cè)試的方法
• 合并和下采樣Seurat對(duì)象的新方法
  2016年8月22日： 1.3版本發(fā)布
• 改進(jìn)的聚類方法 - 有關(guān)詳細(xì)信息，請(qǐng)參閱FAQ
• 所有函數(shù)都支持稀疏矩陣
• 消除不需要的變異來源的方法
• 一致的函數(shù)名稱
• 更新的可視化
  2015年5月21日： Drop-Seq手稿發(fā)布丁鹉。版本1.2發(fā)布
• 增加了對(duì)光譜t-SNE（非線性降維）和密度聚類的支持
• 新的可視化 - 包括pcHeatmap妒潭，dot.plot和feature.plot
• 擴(kuò)展包文檔，減少導(dǎo)入包負(fù)擔(dān)
• Seurat代碼現(xiàn)在托管在GitHub上揣钦，可以通過devtools包輕松安裝
• 小錯(cuò)誤修復(fù)
  2015年4月13日： 空間映射手稿發(fā)布雳灾。版本1.1發(fā)布

Seurat安裝

直接在CRAN加載安裝

# Enter commands in R (or R studio, if installed)
install.packages('Seurat')
library(Seurat)

同時(shí)SeuratV3支持使用UMAP算法進(jìn)行降維聚類可視化，需要配置對(duì)應(yīng)的pythonUMAP環(huán)境冯凹，方可執(zhí)行UMAP對(duì)應(yīng)功能谎亩。

Seurat分析單細(xì)胞RNA-seq數(shù)據(jù)

官網(wǎng)提供了基于10X公司的PBMC樣本的數(shù)據(jù)作為示范；
首先讀入實(shí)例數(shù)據(jù)宇姚，利用Read10X函數(shù)讀取cellranger流程的輸出文件匈庭，返回唯一的分子識(shí)別（UMI）計(jì)數(shù)矩陣。該矩陣中的值表示在每個(gè)細(xì)胞（列）中檢測(cè)到的每個(gè)特征（即基因;行）的分子數(shù)浑劳。然后使用計(jì)數(shù)矩陣來創(chuàng)建一個(gè)Seurat對(duì)象阱持。該對(duì)象充當(dāng)容器，其包含單細(xì)胞數(shù)據(jù)集的數(shù)據(jù)（如計(jì)數(shù)矩陣）和分析（如PCA或聚類結(jié)果）

library(dplyr)
library(Seurat)
setwd("/Users/liuyang/Downloads/03project/02scRNA/10X_Seurat")
rm(list = ls())

# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "/Users/liuyang/Downloads/03project/02scRNA/10X_Seurat/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

標(biāo)準(zhǔn)的預(yù)處理工作流程基于QC指標(biāo)魔熏，數(shù)據(jù)標(biāo)準(zhǔn)化和縮放以及高度可變特征的檢測(cè)來選擇和過濾細(xì)胞衷咽。

QC并選擇細(xì)胞進(jìn)行進(jìn)一步分析

Seurat允許用戶根據(jù)自己定義的標(biāo)準(zhǔn)輕松探索QC指標(biāo)和過濾單元格鸽扁。常用的一些QC指標(biāo)包括

• 在每個(gè)細(xì)胞中檢測(cè)到的獨(dú)特基因的數(shù)量。
• 低質(zhì)量細(xì)胞或空液滴通常只有很少的基因
• 細(xì)胞雙峰或多重峰可能表現(xiàn)出異常高的基因計(jì)數(shù)
• 同樣镶骗，細(xì)胞內(nèi)檢測(cè)到的分子總數(shù)（與獨(dú)特基因強(qiáng)烈相關(guān)）
• 映射到線粒體基因組的讀數(shù)百分比（對(duì)于已經(jīng)組裝到染色體水平或有線粒體基因組的樣本）
• 低質(zhì)量/垂死細(xì)胞通常表現(xiàn)出廣泛的線粒體污染
• 我們使用PercentageFeatureSet函數(shù)計(jì)算線粒體QC指標(biāo)献烦，該函數(shù)計(jì)算源自一組特征的計(jì)數(shù)百分比
• 我們使用所有基因MT-集作為一組線粒體基因

pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
VlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

質(zhì)控

# 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(object = pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(object = pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

feature-feature relationships

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

從數(shù)據(jù)集中刪除不需要的單元格后，下一步是標(biāo)準(zhǔn)化數(shù)據(jù)卖词。默認(rèn)情況下，采用全局縮放規(guī)范化方法LogNormalize吏夯，通過總表達(dá)式對(duì)每個(gè)單元格的要素表達(dá)式度量進(jìn)行標(biāo)準(zhǔn)化此蜈，將其乘以比例因子（默認(rèn)為10,000），并對(duì)結(jié)果進(jìn)行對(duì)數(shù)轉(zhuǎn)換噪生。歸一化值存儲(chǔ)在pbmc[["RNA"]]@data裆赵。

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

識(shí)別高度可變的特征（特征選擇）

然后計(jì)算在數(shù)據(jù)集中表現(xiàn)出高細(xì)胞間差異的特征子集（即，它們?cè)谝恍┘?xì)胞中高度表達(dá)跺嗽，而在其他細(xì)胞中低表達(dá)）战授。作者和其他人研究者發(fā)現(xiàn)，在下游分析中關(guān)注這些基因有助于突出單細(xì)胞數(shù)據(jù)集中的生物信號(hào)桨嫁。

這里詳細(xì)描述了SeuratV3中的程序植兰，并通過直接建模單細(xì)胞數(shù)據(jù)中固有的均值 - 方差關(guān)系對(duì)先前版本進(jìn)行了改進(jìn)，并在FindVariableFeatures函數(shù)中實(shí)現(xiàn)璃吧。默認(rèn)情況下楣导，為每個(gè)數(shù)據(jù)集返回2,000個(gè)要素。這些將用于下游分析畜挨，如PCA筒繁。

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

# Identify the 10 most highly variable genes
top10 <- head(x = VariableFeatures(object = pbmc), 10)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(object = pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
CombinePlots(plots = list(plot1, plot2))

特征選擇

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

接下來，應(yīng)用線性變換（'縮放'）巴元，這是在PCA等降維技術(shù)之前的標(biāo)準(zhǔn)預(yù)處理步驟毡咏。該ScaleData函數(shù)：
改變每個(gè)基因的表達(dá)，使得跨細(xì)胞的平均表達(dá)為0
縮放每個(gè)基因的表達(dá)逮刨，以便跨細(xì)胞的方差為1
該步驟在下游分析中給予相同的權(quán)重呕缭，因此高表達(dá)的基因不占優(yōu)勢(shì)結(jié)果存儲(chǔ)在 pbmc[["RNA"]]@scale.data

執(zhí)行線性降維聚類

對(duì)縮放數(shù)據(jù)執(zhí)行PCA。默認(rèn)情況下修己，只有先前確定的變量要素用作輸入臊旭，但features如果要選擇其他子集，則可以使用參數(shù)定義箩退。
瑟拉提供定義可視化細(xì)胞和功能的幾種有用的方式PCA离熏，包括VizDimReduction，DimPlot戴涝，和DimHeatmap

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

DimPlot(object = pbmc, reduction = "pca")      

PC

特別是DimHeatmap允許容易地探索數(shù)據(jù)集中的異質(zhì)性的主要來源滋戳，并且在試圖決定包括哪些PC以用于進(jìn)一步的下游分析時(shí)可以是有用的钻蔑。細(xì)胞和特征都根據(jù)其PCA分?jǐn)?shù)排序。設(shè)置cells為數(shù)字可以繪制光譜兩端的“極端”單元格奸鸯，從而大大加快了大型數(shù)據(jù)集的繪圖速度咪笑。雖然顯然是一種監(jiān)督分析，但我們發(fā)現(xiàn)這是探索相關(guān)特征集的有用工具娄涩。

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

PC熱圖

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

為了克服scRNA-seq數(shù)據(jù)的任何單一特征中的廣泛技術(shù)噪聲窗怒，Seurat基于其PCA分?jǐn)?shù)對(duì)細(xì)胞進(jìn)行聚類，每個(gè)PC基本上代表在相關(guān)特征集中組合信息的“元特征”蓄拣。因此扬虚，頂部主要組件表示數(shù)據(jù)集的強(qiáng)大壓縮。但是球恤，我們應(yīng)該選擇包含多少個(gè)組件辜昵？10？20咽斧？100堪置？在Macosko 等人中，實(shí)施了一項(xiàng)受JackStraw程序啟發(fā)的重采樣測(cè)試张惹。我們隨機(jī)置換數(shù)據(jù)的子集（默認(rèn)為1％）并重新運(yùn)行PCA舀锨，構(gòu)建特征分?jǐn)?shù)的“空分布”，并重復(fù)此過程宛逗。我們將“重要”PC識(shí)別為具有低p值特征的強(qiáng)大豐富的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(object = pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(object = pbmc, dims = 1:20)

該JackStrawPlot功能提供了一種可視化工具，用于比較每臺(tái)PC的p值分布和均勻分布（虛線）拧额”撸“重要的”PC將顯示具有低p值的特征的強(qiáng)烈豐富（虛線上方的實(shí)線）。在這種情況下侥锦，看起來在前10-12PC之后进栽，重要性急劇下降。

JackStrawPlot(object = pbmc, dims = 1:15)

JackStrawPlot

另一種啟發(fā)式方法生成“肘圖”：基于每個(gè)（ElbowPlot函數(shù)）解釋的方差百分比對(duì)主成分進(jìn)行排序恭垦。在這個(gè)例子中快毛，我們可以觀察PC9-10周圍的“肘部”，這表明大多數(shù)真實(shí)信號(hào)都是在前10臺(tái)PC中捕獲的番挺。

ElbowPlot(object = pbmc)

識(shí)別數(shù)據(jù)集的真實(shí)維度 - 對(duì)用戶來說可能具有挑戰(zhàn)性/不確定性唠帝。因此，我們建議考慮這三種方法玄柏。第一個(gè)是更多的監(jiān)督襟衰，探索PC以確定異質(zhì)性的相關(guān)來源，并且可以與GSEA一起使用粪摘。第二個(gè)實(shí)現(xiàn)基于隨機(jī)空模型的統(tǒng)計(jì)測(cè)試瀑晒，但對(duì)于大型數(shù)據(jù)集來說是耗時(shí)的绍坝，并且可能不會(huì)返回明確的PC截止。第三種是常用的啟發(fā)式算法苔悦，可以立即計(jì)算轩褐。在這個(gè)例子中，所有三種方法產(chǎn)生了類似的結(jié)果玖详，但我們可能有理由在PC 7-12之間選擇任何作為截止值的東西把介。
我們?cè)谶@里選擇了10，但鼓勵(lì)用戶考慮以下內(nèi)容：
樹突細(xì)胞和NK愛好者可以認(rèn)識(shí)到與PC12和13強(qiáng)烈相關(guān)的基因定義了罕見的免疫亞群（即MZB1是漿細(xì)胞樣DC的標(biāo)記）蟋座。然而拗踢，這些組是如此罕見，如果沒有先驗(yàn)知識(shí)蜈七，很難將這種大小的數(shù)據(jù)集與背景噪聲區(qū)分開來。
我們鼓勵(lì)用戶使用不同數(shù)量的PC（10,15甚至50Ｄ！）重復(fù)下游分析飒硅。正如您將觀察到的，結(jié)果通常沒有顯著差異作谚。
在選擇此參數(shù)時(shí)三娩，我們建議用戶在較高的一側(cè)犯錯(cuò)。例如妹懒，僅使用5臺(tái)PC進(jìn)行下游分析確實(shí)會(huì)對(duì)結(jié)果造成嚴(yán)重影響雀监。

聚類細(xì)胞

Seurat v3采用基于圖形的聚類方法，建立在（Macosko 等人）的初始策略之上眨唬。重要的是会前，驅(qū)動(dòng)聚類分析的距離度量（基于先前識(shí)別的PC）保持不變。然而匾竿，我們將細(xì)胞距離矩陣分成簇的方法已經(jīng)大大改進(jìn)瓦宜。我們的方法受到最近手稿的啟發(fā)，這些手稿將基于圖形的聚類方法應(yīng)用于scRNA-seq數(shù)據(jù)[SNN-Cliq岭妖，Xu和Su临庇，Bioinformatics，2015]和CyTOF數(shù)據(jù)[PhenoGraph昵慌，Levine 等假夺，Cell，2015]斋攀。簡(jiǎn)而言之已卷，這些方法將單元格嵌入圖形結(jié)構(gòu)中 - 例如K-最近鄰（KNN）圖形，在具有相似特征表達(dá)模式的單元格之間繪制邊緣淳蔼，然后嘗試將此圖形劃分為高度互連的“準(zhǔn)集團(tuán)”或“社區(qū)”悼尾。
與PhenoGraph一樣柿扣，我們首先根據(jù)PCA空間中的歐氏距離構(gòu)建KNN圖，并根據(jù)其局部鄰域中的共享重疊（Jaccard相似性）細(xì)化任意兩個(gè)單元之間的邊權(quán)重闺魏。使用該FindNeighbors函數(shù)執(zhí)行該步驟未状，并將先前定義的數(shù)據(jù)集維度（前10個(gè)PC）作為輸入。
為了聚類細(xì)胞析桥，我們接下來應(yīng)用模塊化優(yōu)化技術(shù)司草，例如Louvain算法（默認(rèn)）或SLM [SLM，Blondel 等泡仗，Journal of Statistical Mechanics]戈咳，將細(xì)胞迭代分組在一起，目的是優(yōu)化標(biāo)準(zhǔn)模塊化功能贴彼。該FindClusters函數(shù)實(shí)現(xiàn)此過程泛粹，并包含一個(gè)分辨率參數(shù)，用于設(shè)置下游群集的“粒度”截亦，增加的值會(huì)導(dǎo)致更多的群集爬泥。我們發(fā)現(xiàn)將此參數(shù)設(shè)置在0.4-1.2之間通常會(huì)返回大約3K單元格的單細(xì)胞數(shù)據(jù)集的良好結(jié)果。較大的數(shù)據(jù)集通常會(huì)增加最佳分辨率崩瓤∨鄯龋可以使用該Idents功能找到群集。

pbmc <- FindNeighbors(object = pbmc, dims = 1:10)
pbmc <- FindClusters(object = pbmc, resolution = 0.5)

運(yùn)行非線性降維（UMAP / tSNE）

Seurat提供了幾種非線性降維技術(shù)却桶，如tSNE和UMAP境输，可視化和探索這些數(shù)據(jù)集。這些算法的目標(biāo)是學(xué)習(xí)數(shù)據(jù)的基礎(chǔ)流形颖系，以便在低維空間中將相似的單元放在一起嗅剖。上面確定的基于圖的聚類內(nèi)的單元應(yīng)該共同定位這些降維圖。作為UMAP和tSNE的輸入嘁扼，我們建議使用相同的PC作為聚類分析的輸入窗悯。

# If you haven't installed UMAP, you can do so via reticulate::py_install(packages = 'umap-learn')
pbmc <- RunUMAP(object = pbmc, dims = 1:10)
# note that you can set `label = TRUE` or use the LabelClusters function to help label individual clusters
DimPlot(object = pbmc, reduction = "umap")

聚類圖

尋找差異表達(dá)的特征（聚類生物標(biāo)志物）

Seurat可以幫助您找到通過差異表達(dá)定義聚類的標(biāo)記。默認(rèn)情況下偷拔，它標(biāo)識(shí)單個(gè)群集（指定ident.1）的正負(fù)標(biāo)記蒋院，與所有其他單元格進(jìn)行比較。FindAllMarkers為所有群集自動(dòng)執(zhí)行此過程莲绰，但您也可以測(cè)試群集彼此之間或所有單元格欺旧。
該min.pct論證要求在兩組單元中的任一組中以最小百分比檢測(cè)特征，并且thresh.test參數(shù)要求在兩組之間以一定量差異地表示（平均）特征蛤签。您可以將這兩個(gè)設(shè)置為0辞友，但時(shí)間會(huì)急劇增加 - 因?yàn)檫@將測(cè)試大量不太可能具有高度歧視性的功能。作為加速這些計(jì)算的另一種選擇，max.cells.per.ident可以設(shè)置称龙。這將對(duì)每個(gè)標(biāo)識(shí)類進(jìn)行下采樣留拾，使其不再具有比設(shè)置的更多的單元格。雖然通常會(huì)出現(xiàn)功率損失鲫尊，但速度增加可能是顯著的痴柔，并且最高度差異表達(dá)的功能可能仍然會(huì)升至頂部。

cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
head(x = cluster1.markers, n = 5)

Seurat有幾個(gè)差分表達(dá)式測(cè)試疫向，可以使用test.use參數(shù)設(shè)置

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

我們提供了幾種用于可視化標(biāo)記表達(dá)的工具咳蔚。VlnPlot（顯示跨群集的表達(dá)概率分布），以及FeaturePlot（在tSNE或PCA圖上可視化特征表達(dá)）是我們最常用的可視化搔驼。我們還建議您探索RidgePlot谈火，CellScatter以及DotPlot查看數(shù)據(jù)集的其他方法。

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

小提琴圖

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

特征圖

DoHeatmap為給定的單元格和要素生成表達(dá)式熱圖舌涨。在這種情況下糯耍，我們正在為每個(gè)群集繪制前20個(gè)標(biāo)記（或所有標(biāo)記，如果小于20）囊嘉。

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()

熱圖

將單元類型標(biāo)識(shí)分配給集群

根據(jù)每個(gè)類群中高表達(dá)的特異基因温技，我們可以利用已知細(xì)胞類型特定的標(biāo)記基因，對(duì)每個(gè)類群確定細(xì)胞類型哗伯。

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

定細(xì)胞類型

參考材料：

https://satijalab.org/seurat/