第七計(jì) 無中生有
本指本來沒有卻硬說有。現(xiàn)形容憑空捏造叁执。

設(shè)置Seurat對(duì)象

對(duì)于本教程茄厘，我們將分析可從10X Genomics免費(fèi)獲得的外周血單個(gè)核細(xì)胞（PBMC）數(shù)據(jù)集矮冬。在Illumina NextSeq 500上已對(duì)2700個(gè)單細(xì)胞進(jìn)行了測序〈喂可以在此處找到原始數(shù)據(jù)胎署。

我們從讀取數(shù)據(jù)開始。該[Read10X()](https://satijalab.org/seurat/reference/Read10X.html)函數(shù)從10X讀取cellranger管道的輸出窑滞，返回唯一的分子識(shí)別（UMI）計(jì)數(shù)矩陣琼牧。該矩陣中的值表示在每個(gè)單元格（列）中檢測到的每個(gè)特征（即基因；行）的分子數(shù)哀卫。

接下來巨坊，我們使用計(jì)數(shù)矩陣創(chuàng)建一個(gè)Seurat對(duì)象。該對(duì)象用作包含單個(gè)單元數(shù)據(jù)集的數(shù)據(jù)（如計(jì)數(shù)矩陣）和分析（如PCA或聚類結(jié)果）的容器此改。有關(guān)Seurat對(duì)象結(jié)構(gòu)的技術(shù)討論趾撵，請(qǐng)查看我們的GitHub Wiki。例如共啃，計(jì)數(shù)矩陣存儲(chǔ)在中pbmc[["RNA"]]@counts占调。

library(dplyr)
library(Seurat)
library(patchwork)

# 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, 0 variable features)

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

以下步驟涵蓋了Seurat中scRNA-seq數(shù)據(jù)的標(biāo)準(zhǔn)預(yù)處理工作流程。這些代表基于質(zhì)量控制指標(biāo)移剪，數(shù)據(jù)歸一化和縮放以及高度可變特征的檢測來選擇和過濾單元究珊。

質(zhì)控和選擇細(xì)胞進(jìn)行進(jìn)一步分析

Seurat允許您輕松瀏覽質(zhì)量控制指標(biāo)并根據(jù)任何用戶定義的標(biāo)準(zhǔn)過濾單元。社區(qū)常用的一些質(zhì)量控制指標(biāo)包括

• 每個(gè)細(xì)胞中檢測到的獨(dú)特基因的數(shù)量纵苛。
  • 低質(zhì)量的細(xì)胞或空液滴通常很少有基因
  • 細(xì)胞雙峰或多重峰可能顯示異常高的基因計(jì)數(shù)
• 同樣剿涮，細(xì)胞內(nèi)檢測到的分子總數(shù)（與獨(dú)特基因高度相關(guān)）
• 映射到線粒體基因組的讀段的百分比
  • 低質(zhì)量/瀕臨死亡的細(xì)胞通常表現(xiàn)出廣泛的線粒體污染
  • 我們使用[PercentageFeatureSet()](https://satijalab.org/seurat/reference/PercentageFeatureSet.html)函數(shù)計(jì)算線粒體QC指標(biāo)，該函數(shù)計(jì)算源自一組功能的計(jì)數(shù)的百分比
  • 我們使用所有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指標(biāo)幔虏，并使用它們來過濾單元格。

• 我們過濾具有超過2500或少于200的唯一特征計(jì)數(shù)的像元
• 我們過濾線粒體計(jì)數(shù)> 5％的細(xì)胞

# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

image

# 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")
plot1 + plot2

image

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

規(guī)范化數(shù)據(jù)

從數(shù)據(jù)集中刪除不需要的單元格后贝椿，下一步就是將數(shù)據(jù)標(biāo)準(zhǔn)化想括。默認(rèn)情況下，我們采用全局縮放歸一化方法“ LogNormalize”烙博，該方法將每個(gè)單元格的特征表達(dá)式測量結(jié)果與總表達(dá)式進(jìn)行歸一化瑟蜈，再乘以比例因子（默認(rèn)為10,000），然后對(duì)結(jié)果進(jìn)行對(duì)數(shù)轉(zhuǎn)換渣窜。規(guī)范化的值存儲(chǔ)在中pbmc[["RNA"]]@data铺根。

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

為了清楚起見，在前面的代碼行（以及以后的命令）中乔宿，我們在函數(shù)調(diào)用中提供了某些參數(shù)的默認(rèn)值位迂。但是，這不是必需的，并且可以通過以下方式實(shí)現(xiàn)相同的行為：

pbmc <- NormalizeData(pbmc)

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

接下來掂林，我們計(jì)算特征的子集臣缀，這些特征在數(shù)據(jù)集中顯示出較高的單元間差異（即，它們在某些單元格中高表達(dá)泻帮，而在其他單元格中低表達(dá)）精置。我們和其他人發(fā)現(xiàn)，在下游分析中關(guān)注這些基因有助于在單細(xì)胞數(shù)據(jù)集中突出顯示生物信號(hào)锣杂。

我們在Seurat中的過程在此進(jìn)行了詳細(xì)描述脂倦，并通過直接建模單細(xì)胞數(shù)據(jù)中固有的均值-方差關(guān)系對(duì)以前的版本進(jìn)行了改進(jìn)，并在[FindVariableFeatures()](https://satijalab.org/seurat/reference/FindVariableFeatures.html)函數(shù)中實(shí)現(xiàn)了該過程元莫。默認(rèn)情況下赖阻，我們?yōu)槊總€(gè)數(shù)據(jù)集返回2,000個(gè)要素。這些將用于下游分析柒竞，例如PCA政供。

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

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

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

image

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

接下來播聪，我們應(yīng)用線性變換（“縮放”）朽基，這是像PCA這樣的降維技術(shù)之前的標(biāo)準(zhǔn)預(yù)處理步驟。該[ScaleData()](https://satijalab.org/seurat/reference/ScaleData.html)函數(shù)：

• 移動(dòng)每個(gè)基因的表達(dá)离陶，以使整個(gè)細(xì)胞的平均表達(dá)為0
• 縮放每個(gè)基因的表達(dá)稼虎，從而使細(xì)胞之間的方差為1
  • 此步驟在下游分析中具有相等的權(quán)重，因此高表達(dá)的基因不會(huì)占主導(dǎo)地位
• 結(jié)果存儲(chǔ)在 pbmc[["RNA"]]@scale.data

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

<details style="box-sizing: border-box; display: block;"><summary style="box-sizing: border-box; display: list-item;">此步驟耗時(shí)太長招刨！我可以加快速度嗎霎俩？</summary></details> <details style="box-sizing: border-box; display: block;"><summary style="box-sizing: border-box; display: list-item;">與Seurat v2中一樣，如何刪除不需要的變化來源沉眶？</summary></details>

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

接下來打却，我們對(duì)縮放后的數(shù)據(jù)執(zhí)行PCA。默認(rèn)情況下谎倔，僅將先前確定的變量特征用作輸入柳击，但是features如果您希望選擇其他子集，則可以使用arguments進(jìn)行定義片习。

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

瑟拉提供可視化細(xì)胞和定義PCA捌肴，包括功能的幾種有用的方法VizDimReduction()，[DimPlot()](https://satijalab.org/seurat/reference/DimPlot.html)和[DimHeatmap()](https://satijalab.org/seurat/reference/DimHeatmap.html)

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

image

DimPlot(pbmc, reduction = "pca")

image

特別是[DimHeatmap()](https://satijalab.org/seurat/reference/DimHeatmap.html)可以輕松探索數(shù)據(jù)集中異質(zhì)性的主要來源藕咏，并且在嘗試確定要包括哪些PC以便進(jìn)行進(jìn)一步的下游分析時(shí)非常有用状知。單元和要素均根據(jù)其PCA分?jǐn)?shù)排序。設(shè)置cells為數(shù)字會(huì)在頻譜的兩端繪制“極端”單元孽查，從而極大地加快了大型數(shù)據(jù)集的繪制速度饥悴。盡管顯然是監(jiān)督分析，但我們發(fā)現(xiàn)這是探索相關(guān)特征集的有價(jià)值的工具。

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

image

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

image

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

為了克服scRNA-seq數(shù)據(jù)的任何單個(gè)特征中的廣泛技術(shù)噪聲西设，Seurat根據(jù)其PCA分?jǐn)?shù)對(duì)細(xì)胞進(jìn)行聚類起宽，每個(gè)PC本質(zhì)上代表一種“元特征”，該特征將跨相關(guān)特征集的信息進(jìn)行組合济榨。因此坯沪，最主要的主成分代表了數(shù)據(jù)集的穩(wěn)健壓縮。但是擒滑，我們應(yīng)該選擇包括多少個(gè)組件腐晾？10個(gè)？20嗎 100丐一？

在Macosko等人中藻糖，我們實(shí)施了受JackStraw程序啟發(fā)的重采樣測試。我們隨機(jī)置換數(shù)據(jù)的一部分（默認(rèn)為1％）库车，然后重新運(yùn)行PCA巨柒，以構(gòu)建特征分?jǐn)?shù)的“零分布”，然后重復(fù)此過程柠衍。我們將“重要”的PC識(shí)別為具有豐富的低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()](https://satijalab.org/seurat/reference/JackStrawPlot.html)功能提供了一個(gè)可視化工具，用于比較每個(gè)PC的p值分布與均勻分布（虛線）珍坊∥矗“重要”的PC將顯示出具有豐富的低p值功能（虛線上方的實(shí)線）。在這種情況下阵漏，似乎在前10到12臺(tái)PC之后驻民，重要性顯著下降。

JackStrawPlot(pbmc, dims = 1:15)

image

另一種啟發(fā)式方法會(huì)生成“肘形圖”：基于每個(gè)分量（[ElbowPlot()](https://satijalab.org/seurat/reference/ElbowPlot.html)函數(shù)）所解釋的方差百分比對(duì)主成分進(jìn)行排序履怯。在此示例中回还，我們可以在PC9-10周圍觀察到一個(gè)“彎頭”，這表明大多數(shù)真實(shí)信號(hào)是在前10臺(tái)PC中捕獲的叹洲。

ElbowPlot(pbmc)

image

識(shí)別數(shù)據(jù)集的真實(shí)維度–對(duì)用戶而言可能是挑戰(zhàn)/不確定柠硕。因此，我們建議考慮這三種方法疹味。第一個(gè)受到更多監(jiān)督仅叫，探索PC以確定異質(zhì)性的相關(guān)來源，例如可以與GSEA結(jié)合使用糙捺。第二種實(shí)現(xiàn)基于隨機(jī)空模型的統(tǒng)計(jì)測試诫咱，但是對(duì)于大型數(shù)據(jù)集而言非常耗時(shí)，并且可能無法返回清晰的PC截止值洪灯。第三種是常用的啟發(fā)式方法坎缭，可以立即進(jìn)行計(jì)算竟痰。在此示例中，所有三種方法都產(chǎn)生了相似的結(jié)果掏呼，但是我們可能有理由在PC 7-12之間選擇任何東西作為臨界值坏快。

我們在這里選擇10，但鼓勵(lì)用戶考慮以下內(nèi)容：

• 樹突狀細(xì)胞和NK愛好者可能認(rèn)識(shí)到與PC 12和13緊密相關(guān)的基因定義了罕見的免疫亞群（即MZB1是漿細(xì)胞樣DC的標(biāo)記）憎夷。但是莽鸿，這些組非常罕見，在沒有先驗(yàn)知識(shí)的情況下拾给，很難將它們與背景噪聲區(qū)分開來祥得。
• 我們鼓勵(lì)用戶使用不同數(shù)量的PC（10臺(tái)，15臺(tái)蒋得，甚至50臺(tái)＜都啊）重復(fù)進(jìn)行下游分析。正如您將看到的额衙，結(jié)果通常并沒有太大的不同饮焦。
• 我們建議用戶選擇此參數(shù)時(shí)偏高一點(diǎn)。例如窍侧，僅用5臺(tái)PC進(jìn)行下游分析確實(shí)會(huì)對(duì)結(jié)果產(chǎn)生不利影響县踢。

聚集細(xì)胞

Seurat v3在（Macosko等人）的初始策略的基礎(chǔ)上，應(yīng)用了基于圖的聚類方法疏之。重要的是殿雪，驅(qū)動(dòng)聚類分析的距離度量（基于先前確定的PC）保持不變。但是锋爪，我們將蜂窩距離矩陣劃分為群集的方法已得到了極大的改進(jìn)。我們的方法受到最近手稿的啟發(fā)爸业，這些手稿將基于圖的聚類方法應(yīng)用于scRNA-seq數(shù)據(jù)[SNN-Cliq其骄，Xu和Su，Bioinformatics扯旷，2015]和CyTOF數(shù)據(jù)[PhenoGraph拯爽，Levine等人，Cell钧忽，2015]毯炮。簡而言之，這些方法將單元格嵌入圖結(jié)構(gòu)（例如耸黑，K近鄰圖（KNN））桃煎，并在具有相似特征表達(dá)模式的單元格之間繪制邊，然后嘗試將該圖劃分為高度互連的“準(zhǔn)斜體”或“社區(qū)”大刊。

像在PhenoGraph中一樣为迈，我們首先基于PCA空間中的歐式距離構(gòu)造一個(gè)KNN圖，然后基于兩個(gè)相鄰像元在其局部鄰域中的共享重疊來細(xì)化任意兩個(gè)像元之間的邊緣權(quán)重（Jaccard相似性）。使用該[FindNeighbors()](https://satijalab.org/seurat/reference/FindNeighbors.html)功能執(zhí)行此步驟葫辐，并將先前定義的數(shù)據(jù)集維度（前10個(gè)PC）作為輸入搜锰。

為了對(duì)單元進(jìn)行聚類，我們接下來將應(yīng)用模塊化優(yōu)化技術(shù)（例如Louvain算法（默認(rèn)）或SLM [SLM耿战，Blondel等蛋叼，統(tǒng)計(jì)力學(xué)雜志]）將單元迭代地分組在一起，以優(yōu)化標(biāo)準(zhǔn)模塊化功能剂陡。該[FindClusters()](https://satijalab.org/seurat/reference/FindClusters.html)函數(shù)實(shí)現(xiàn)此過程鸦列，并包含一個(gè)分辨率參數(shù)，該參數(shù)設(shè)置下游集群的“粒度”鹏倘，其值越大薯嗤，導(dǎo)致集群的數(shù)量越多。我們發(fā)現(xiàn)纤泵，將此參數(shù)設(shè)置在0.4-1.2之間通陈娼悖可為大約3K個(gè)單元的單單元數(shù)據(jù)集返回良好的結(jié)果。對(duì)于較大的數(shù)據(jù)集捏题，最佳分辨率通常會(huì)增加玻褪。可以使用[Idents()](https://satijalab.org/seurat/reference/Idents.html)函數(shù)找到簇公荧。

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: 95965
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8723
## Number of communities: 9
## Elapsed time: 0 seconds

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

## AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 
##                2                3                2                1 
## AAACCGTGTATGCG-1 
##                6 
## Levels: 0 1 2 3 4 5 6 7 8

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

Seurat提供了幾種非線性降維技術(shù)带射，例如tSNE和UMAP，以可視化和探索這些數(shù)據(jù)集循狰。這些算法的目標(biāo)是學(xué)習(xí)數(shù)據(jù)的基礎(chǔ)流形窟社，以便將相似的單元格放置在低維空間中。上面確定的基于圖的聚類中的像元應(yīng)共位于這些降維圖上绪钥。作為UMAP和tSNE的輸入灿里，我們建議使用相同的PC作為聚類分析的輸入。

# If you haven't installed UMAP, you can do so via reticulate::py_install(packages =
# 'umap-learn')
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")

image

您可以在此時(shí)保存該對(duì)象程腹，以便可以輕松地將其重新加載匣吊，而不必重新運(yùn)行上面執(zhí)行的計(jì)算量大的步驟，也可以輕松地與協(xié)作者共享該對(duì)象寸潦。

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

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

Seurat可以幫助您找到通過差異表達(dá)定義聚類的標(biāo)記色鸳。默認(rèn)情況下，ident.1與所有其他單元格相比见转，它識(shí)別單個(gè)群集（在中指定）的正向和負(fù)向標(biāo)記命雀。[FindAllMarkers()](https://satijalab.org/seurat/reference/FindAllMarkers.html)自動(dòng)執(zhí)行所有集群的此過程，但是您也可以測試集群組之間的相互關(guān)系池户，或針對(duì)所有單元進(jìn)行測試咏雌。

該min.pct參數(shù)要求在兩組像元的任何一組中以最小百分比檢測特征凡怎，而thresh.test參數(shù)要求在兩組像元之間平均表達(dá)（平均）一定程度的特征。您可以將它們都設(shè)置為0赊抖，但是時(shí)間會(huì)大大增加-因?yàn)檫@將測試可能不太具有高度歧視性的大量功能统倒。作為加快這些計(jì)算速度的另一種選擇，max.cells.per.ident可以設(shè)置氛雪。這將降低每個(gè)身份類的采樣率房匆，使其沒有更多的單元格。盡管通常會(huì)斷電报亩，但速度的提高可能會(huì)非常明顯浴鸿，而差異化程度最高的功能可能仍會(huì)上升到最高位置。

# find all markers of cluster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster1.markers, n = 5)

##             p_val avg_log2FC pct.1 pct.2    p_val_adj
## IL32 2.593535e-91  1.2154360 0.949 0.466 3.556774e-87
## LTB  7.994465e-87  1.2828597 0.981 0.644 1.096361e-82
## CD3D 3.922451e-70  0.9359210 0.922 0.433 5.379250e-66
## IL7R 1.130870e-66  1.1776027 0.748 0.327 1.550876e-62
## LDHB 4.082189e-65  0.8837324 0.953 0.614 5.598314e-61

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

##                       p_val avg_log2FC pct.1 pct.2     p_val_adj
## FCGR3A        2.150929e-209   4.267579 0.975 0.039 2.949784e-205
## IFITM3        6.103366e-199   3.877105 0.975 0.048 8.370156e-195
## CFD           8.891428e-198   3.411039 0.938 0.037 1.219370e-193
## CD68          2.374425e-194   3.014535 0.926 0.035 3.256286e-190
## RP11-290F20.3 9.308287e-191   2.722684 0.840 0.016 1.276538e-186

# 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_log2FC)

## # A tibble: 18 x 7
## # Groups:   cluster [9]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene    
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>   
##  1 1.74e-109       1.07 0.897 0.593 2.39e-105 0       LDHB    
##  2 1.17e- 83       1.33 0.435 0.108 1.60e- 79 0       CCR7    
##  3 0\.              5.57 0.996 0.215 0\.        1       S100A9  
##  4 0\.              5.48 0.975 0.121 0\.        1       S100A8  
##  5 7.99e- 87       1.28 0.981 0.644 1.10e- 82 2       LTB     
##  6 2.61e- 59       1.24 0.424 0.111 3.58e- 55 2       AQP3    
##  7 0\.              4.31 0.936 0.041 0\.        3       CD79A   
##  8 9.48e-271       3.59 0.622 0.022 1.30e-266 3       TCL1A   
##  9 1.17e-178       2.97 0.957 0.241 1.60e-174 4       CCL5    
## 10 4.93e-169       3.01 0.595 0.056 6.76e-165 4       GZMK    
## 11 3.51e-184       3.31 0.975 0.134 4.82e-180 5       FCGR3A  
## 12 2.03e-125       3.09 1     0.315 2.78e-121 5       LST1    
## 13 1.05e-265       4.89 0.986 0.071 1.44e-261 6       GZMB    
## 14 6.82e-175       4.92 0.958 0.135 9.36e-171 6       GNLY    
## 15 1.48e-220       3.87 0.812 0.011 2.03e-216 7       FCER1A  
## 16 1.67e- 21       2.87 1     0.513 2.28e- 17 7       HLA-DPB1
## 17 7.73e-200       7.24 1     0.01  1.06e-195 8       PF4     
## 18 3.68e-110       8.58 1     0.024 5.05e-106 8       PPBP

Seurat可以使用test.use參數(shù)設(shè)置一些針對(duì)差異表達(dá)的測試（有關(guān)詳細(xì)信息弦追，請(qǐng)參見我們的DE小插圖）岳链。例如，ROC測試返回任何單個(gè)標(biāo)記的“分類能力”（范圍從0-隨機(jī)劲件，到1-完美）掸哑。

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

我們提供了幾種可視化標(biāo)記表達(dá)的工具。[VlnPlot()](https://satijalab.org/seurat/reference/VlnPlot.html)（顯示整個(gè)集群中的表達(dá)概率分布）零远，以及[FeaturePlot()](https://satijalab.org/seurat/reference/FeaturePlot.html)（在tSNE或PCA圖上可視化特征表達(dá)）是我們最常用的可視化苗分。我們還建議探索[RidgePlot()](https://satijalab.org/seurat/reference/RidgePlot.html)，[CellScatter()](https://satijalab.org/seurat/reference/CellScatter.html)和[DotPlot()](https://satijalab.org/seurat/reference/DotPlot.html)作為查看數(shù)據(jù)集的其他方法牵辣。

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

image

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

image

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

image

[DoHeatmap()](https://satijalab.org/seurat/reference/DoHeatmap.html)生成給定單元和特征的表達(dá)式熱圖摔癣。在這種情況下，我們將為每個(gè)聚類繪制前20個(gè)標(biāo)記（如果小于20纬向，則繪制所有標(biāo)記）择浊。

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

image

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

幸運(yùn)的是，對(duì)于該數(shù)據(jù)集罢猪，我們可以使用規(guī)范標(biāo)記輕松地將無偏聚類與已知單元格類型進(jìn)行匹配：

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

image

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