上一篇文章深入学习了Monocle2包，这次来继续学习Seurat包。
Seurat官方网站：https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html

Seurat包是个什么包？
Seurat是一个用于质量控制、分析和探索单细胞RNA-seq数据的R包。Seurat使用户能够鉴定和解释来自单细胞转录本测定的异质性来源，并可以整合成多种类型的单细胞数据。

安装Seurat

使用之前先安装。目前最新的是3.0版本，你可以直接在R里安装（要求R是3.4版本或更高），super easy：

> install.packages('Seurat')
> library(Seurat)

如果有些同学想安装低版本的怎么办呢？比如2.3.4这个版本：

> source("https://z.umn.edu/archived-seurat")

如果你还想安装更老的版本，就需要如下代码：

> install.packages('devtools')
# 替换成2.3.0版本
>devtools::install_version(package = 'Seurat', version = package_version('2.3.0'))
> library(Seurat)

建立Seurat对象

在这个教程里，我们将分析来自10X Genomics外周血单核细胞(PBMC)数据，这个dataset是用Illumina NextSeq 500平台进行测序，包含了2700个细胞。原始数据在这里下载：here.

这个原始数据来自CellRanger的处理，Cellranger返回一个UMI的count矩阵。矩阵里的列是细胞，行是基因。接下来我们用count矩阵创建Seurat对象。这个对象就好比一个容器，里面装着单细胞数据集，比如count矩阵，PCA，聚类结果等等。举个栗子：count matrix 储存在pbmc[["RNA"]]@counts里。

> library(dplyr)
> library(Seurat)

# 导入PBMC数据，用Read10×这个参数
> pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
# 起始的Seurat对象是原始数据，没有经过标准化的
> 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)

如果在R studio打开后，大概长这样：

那我想看看具体的某个基因在这个矩阵里每个细胞的表达情况怎么看呢?

> pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
3 x 30 sparse Matrix of class "dgCMatrix"
   [[ suppressing 30 column names 'AAACATACAACCAC', 'AAACATTGAGCTAC', 'AAACATTGATCAGC' ... ]]
                                                                   
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，意思是这些基因在这个细胞里没有被检测到。因为大部分基因在单细胞RNA-seq矩阵里都是0，所以Seurat尽可能使用稀疏矩阵（sparse-matrix）表示。这为Drop-seq/inDrop/10x数据节省了大量内存和速度。

标准的预处理工作流程

质量控制，选择需要分析的细胞

Seurat允许你简单的过滤一下你的细胞，质控的一般指标包括：

• 在每个细胞里检测到的基因
  • 低质量细胞或者空的droplets只有非常少量的基因
  • 细胞doublets 或者 multiplets 有非常高的基因count数
• 在一个细胞内检测到的分子数
• 读取到线粒体基因组的比例
  • 低质量/死细胞通常有很高的线粒体污染
  • 我们使用“PercentageFeatureSet”函数计算线粒体QC指标
  • 使用所有以“MT-”开头的基因作为线粒体基因

# 把表达矩阵里的线粒体QC单独存放在Seurat对象里，等于添加一列到metadata的对象里
> pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

添加后我们可以看一下，是否添加成功了：

#看metadata的前5行
> head(pbmc@meta.data, 5)
               orig.ident nCount_RNA nFeature_RNA percent.mt
AAACATACAACCAC     pbmc3k       2419          779  3.0177759
AAACATTGAGCTAC     pbmc3k       4903         1352  3.7935958
AAACATTGATCAGC     pbmc3k       3147         1129  0.8897363
AAACCGTGCTTCCG     pbmc3k       2639          960  1.7430845
AAACCGTGTATGCG     pbmc3k        980          521  1.2244898

我们还可以可视化QC指标，并用它们来过滤细胞：
（1）将unique基因count数超过2500，或者小于200的细胞过滤掉。
（2）把线粒体count数占5%以上的细胞过滤掉

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

我们还可以用FeatureScatter函数来可视化特征-特征之间的关系，可以使用Seurat对象里的任何东西，如对象中的列、PC分数等。

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

标准化数据

在去除掉不想要的细胞后，就可以标准化数据了。在默认情况下，我们使用global-scaling标准化方法，称为“LogNormalize”，这种方法是利用总的表达量对每个细胞里的基因表达值进行标准化，乘以一个scale factor（默认值是10000），再用log转换一下。标准化后的数据存放在pbmc[["RNA"]]@data里。

> pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

当然上面的代码中，我们用的参数都是默认值，在这种情况下，上面这行代码也可以写成更简单的形式：

> pbmc <- NormalizeData(pbmc)

鉴定高度变化的基因

接下来我们要计算在pbmc里细胞之间变化度很高的基因集（这些基因在一些细胞中高表达，在另一些细胞中低表达），这一步用FindVariableFeatures函数来执行。默认情况下，每个dataset返回2,000个基因。这些基因将用于下游分析，如PCA。

> pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

之后从这些基因里挑选10个变化程度最高的基因：

> top10 <- head(VariableFeatures(pbmc), 10)

下面是将这些基因可视化出来，并把top10的基因名称标记出来：

> plot1 <- VariableFeaturePlot(pbmc)
> plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
When using repel, set xnudge and ynudge to 0 for optimal results
> CombinePlots(plots = list(plot1, plot2))

scaling the data（数据比例缩放）

接下来，用一个线性变换(“缩放scaling”)。这是在降维前（例如PCA）一个标准的预处理步骤，用ScaleData函数:
（1）shift每个基因的表达，使细胞间的平均表达为0
（2）缩放每个基因的表达，使细胞间的差异为1
这一步给予下游分析同等的权重，这样那些非常高表达的基因就不会掩盖其他基因的变化，其结果存储在pbmc[["RNA"]]@scale.data中。

> all.genes <- rownames(pbmc)
> pbmc <- ScaleData(pbmc, features = all.genes)
Centering and scaling data matrix
  |====================================================================| 100%

scaling是Seurat流程中的一个重要步骤，但仅限于将作为PCA input的基因。因此，ScaleData中的默认值只是对前面确定的2000个基因执行缩放。要做到这一点，可以忽略前面函数调用中的features参数，用下面简单的代码就行：

pbmc < - ScaleData (pbmc)

这样你的PCA和聚类结果将不受影响。然而，Seurat heatmap需要对heatmap中的基因进行缩放，以确保高表达基因不会主导heatmap。为了确保我们在后面的热图中不会遗漏任何基因，我们将在本教程中缩放了所有基因。

线性降维

接下来我们对scaling后的数据进行PCA，在默认值下，只对前面选出的2000个基因基因降维，但是你也可以选择想要降维的基因集。

> pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
PC_ 1 
Positive:  CST3, TYROBP, LST1, AIF1, FTL, FTH1, LYZ, FCN1, S100A9, TYMP 
       FCER1G, CFD, LGALS1, S100A8, CTSS, LGALS2, SERPINA1, IFITM3, SPI1, CFP 
       PSAP, IFI30, SAT1, COTL1, S100A11, NPC2, GRN, LGALS3, GSTP1, PYCARD 
Negative:  MALAT1, LTB, IL32, IL7R, CD2, B2M, ACAP1, CD27, STK17A, CTSW 
       CD247, GIMAP5, AQP3, CCL5, SELL, TRAF3IP3, GZMA, MAL, CST7, ITM2A 
       MYC, GIMAP7, HOPX, BEX2, LDLRAP1, GZMK, ETS1, ZAP70, TNFAIP8, RIC3 
PC_ 2 
Positive:  CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1, HLA-DRA, LINC00926, CD79B, HLA-DRB1, CD74 
       HLA-DMA, HLA-DPB1, HLA-DQA2, CD37, HLA-DRB5, HLA-DMB, HLA-DPA1, FCRLA, HVCN1, LTB 
       BLNK, P2RX5, IGLL5, IRF8, SWAP70, ARHGAP24, FCGR2B, SMIM14, PPP1R14A, C16orf74 
Negative:  NKG7, PRF1, CST7, GZMB, GZMA, FGFBP2, CTSW, GNLY, B2M, SPON2 
       CCL4, GZMH, FCGR3A, CCL5, CD247, XCL2, CLIC3, AKR1C3, SRGN, HOPX 
       TTC38, APMAP, CTSC, S100A4, IGFBP7, ANXA1, ID2, IL32, XCL1, RHOC 
PC_ 3 
Positive:  HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1, HLA-DPA1, CD74, MS4A1, HLA-DRB1, HLA-DRA 
       HLA-DRB5, HLA-DQA2, TCL1A, LINC00926, HLA-DMB, HLA-DMA, CD37, HVCN1, FCRLA, IRF8 
       PLAC8, BLNK, MALAT1, SMIM14, PLD4, LAT2, IGLL5, P2RX5, SWAP70, FCGR2B 
Negative:  PPBP, PF4, SDPR, SPARC, GNG11, NRGN, GP9, RGS18, TUBB1, CLU 
       HIST1H2AC, AP001189.4, ITGA2B, CD9, TMEM40, PTCRA, CA2, ACRBP, MMD, TREML1 
       NGFRAP1, F13A1, SEPT5, RUFY1, TSC22D1, MPP1, CMTM5, RP11-367G6.3, MYL9, GP1BA 
PC_ 4 
Positive:  HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1, CD74, HLA-DPB1, HIST1H2AC, PF4, TCL1A 
       SDPR, HLA-DPA1, HLA-DRB1, HLA-DQA2, HLA-DRA, PPBP, LINC00926, GNG11, HLA-DRB5, SPARC 
       GP9, AP001189.4, CA2, PTCRA, CD9, NRGN, RGS18, GZMB, CLU, TUBB1 
Negative:  VIM, IL7R, S100A6, IL32, S100A8, S100A4, GIMAP7, S100A10, S100A9, MAL 
       AQP3, CD2, CD14, FYB, LGALS2, GIMAP4, ANXA1, CD27, FCN1, RBP7 
       LYZ, S100A11, GIMAP5, MS4A6A, S100A12, FOLR3, TRABD2A, AIF1, IL8, IFI6 
PC_ 5 
Positive:  GZMB, NKG7, S100A8, FGFBP2, GNLY, CCL4, CST7, PRF1, GZMA, SPON2 
       GZMH, S100A9, LGALS2, CCL3, CTSW, XCL2, CD14, CLIC3, S100A12, CCL5 
       RBP7, MS4A6A, GSTP1, FOLR3, IGFBP7, TYROBP, TTC38, AKR1C3, XCL1, HOPX 
Negative:  LTB, IL7R, CKB, VIM, MS4A7, AQP3, CYTIP, RP11-290F20.3, SIGLEC10, HMOX1 
       PTGES3, LILRB2, MAL, CD27, HN1, CD2, GDI2, ANXA5, CORO1B, TUBA1B 
       FAM110A, ATP1A1, TRADD, PPA1, CCDC109B, ABRACL, CTD-2006K23.1, WARS, VMO1, FYB 

Seurat 提供几种方法可以对PCA分析后的细胞和基因进行可视化：VizDimReduction, DimPlot, 和DimHeatmap。

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

> DimPlot(pbmc, reduction = "pca")

下面的热图里，可以帮你简单的探索一下dataset里异质性的来源，在你决定下游用哪些PC进行分析时会非常有用。细胞和基因根据它们的PCA值进行排序。用cells这个参数可以定义在排序里最高的和最低的一些细胞：

#只画第一个主成分
> DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)

也可以画很多个PC:

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

只截了一部分图

决定“维数”

为了克服scRNA-seq数据中大量的技术噪声，Seurat基于它们的PCA值对细胞进行聚类。那么怎么才能知道我们需要几个主成分进行分析呢？
有两种方法：

（1）用JackStrawPlot函数

它可以比较每个PC的p-值的分布，显著的PC会有很低的p value，：

#这一步需要很长的时间
> pbmc <- JackStraw(pbmc, num.replicate = 100)
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04m 34s
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s  
> pbmc <- ScoreJackStraw(pbmc, dims = 1:20)

可视化：

> JackStrawPlot(pbmc, dims = 1:15)

在这图里可以看出，显著性从第10个PC到12个PC就开始急剧下降了

（2）Elbow plot法

弯道图主要是看画出来的点在哪一个PC处“转弯”：

> ElbowPlot(pbmc)

这里我们选择了10，但是我们建议在用户分析自己的数据时，选择不同的PC数来重复下游分析（10,15，甚至50），这里也不要选择过少的PC数，会显著影响下游的分析结果。

细胞聚类

首先我们在PCA空间里根据欧氏距离构建一个KNN图，并在任意两个细胞间要确定它们的边缘权重，这个过程用FindNeighbors函数，这里的input就是前面我们定义的dataset的维数。
接下来再用FindClusters函数，该函数有一个“分辨率”的参数，该参数设置下游聚类的“粒度”，值越高，得到的聚类数越多。这个参数设置在0.4-1.2之间，对于3千个左右的单细胞数据通常会得到比较好的结果。对于较大的数据集，最佳分辨率通常会增加。

> pbmc <- FindNeighbors(pbmc, dims = 1:10)
Computing nearest neighbor graph
Computing SNN
> 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...
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8720
Number of communities: 9
Elapsed time: 0 seconds

#用Idents函数看一下前5个细胞的聚类情况
> head(Idents(pbmc), 5)
AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG AAACCGTGTATGCG 
             1              3              1              2              6 
Levels: 0 1 2 3 4 5 6 7 8

非线性降维

Seurat提供了几种非线性的降维技术，如tSNE和UMAP。这些算法的目标是在低维空间中将相似的细胞放在一起。上面所确定的基于图的集群中的单元应该在这些降维图上共同定位。作为UMAP和tSNE的input，我们建议使用与聚类分析的输入相同的PCs作为输入。

> pbmc <- RunUMAP(pbmc, dims = 1:10)
15:28:31 UMAP embedding parameters a = 0.9922 b = 1.112
15:28:31 Read 2638 rows and found 10 numeric columns
15:28:31 Using Annoy for neighbor search, n_neighbors = 30
15:28:31 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:28:32 Writing NN index file to temp file /tmp/RtmpoZ1H8s/file30314c4df890
15:28:32 Searching Annoy index using 1 thread, search_k = 3000
15:28:33 Annoy recall = 100%
15:28:33 Commencing smooth kNN distance calibration using 1 thread
15:28:34 Initializing from normalized Laplacian + noise
15:28:34 Commencing optimization for 500 epochs, with 105386 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:28:43 Optimization finished
> DimPlot(pbmc, reduction = "umap")

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

寻找差异表达基因

Seurat可以帮助你定义cluster的差异表达，默认情况下，它定义单个cluster的阳性和阴性的marker（与其他细胞群比较）。FindAllMarkers可以执行这个过程。

#在cluster1里寻找marker
#min.pct参数定义了一个基因在两组细胞任意一组里的最低可检测到的百分率
> cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04s  
> head(cluster1.markers, n = 5)
            p_val avg_logFC pct.1 pct.2    p_val_adj
IL32 1.894810e-92 0.8373872 0.948 0.464 2.598542e-88
LTB  7.953303e-89 0.8921170 0.981 0.642 1.090716e-84
CD3D 1.655937e-70 0.6436286 0.919 0.431 2.270951e-66
IL7R 3.688893e-68 0.8147082 0.747 0.325 5.058947e-64
LDHB 2.292819e-67 0.6253110 0.950 0.613 3.144372e-63

#寻找cluster5与cluster0和3之间的差异marker
> cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=10s  
> head(cluster5.markers, n = 5)
                      p_val avg_logFC pct.1 pct.2     p_val_adj
FCGR3A        7.583625e-209  2.963144 0.975 0.037 1.040018e-204
IFITM3        2.500844e-199  2.698187 0.975 0.046 3.429657e-195
CFD           1.763722e-195  2.362381 0.938 0.037 2.418768e-191
CD68          4.612171e-192  2.087366 0.926 0.036 6.325132e-188
RP11-290F20.3 1.846215e-188  1.886288 0.840 0.016 2.531900e-184

你还可以寻找每一个cluster里与其他所有细胞相比之后的差异marker:

> pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
Calculating cluster 0
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s  
Calculating cluster 1
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  
Calculating cluster 2
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=08s  
Calculating cluster 3
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s  
Calculating cluster 4
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s  
Calculating cluster 5
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=11s  
Calculating cluster 6
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=09s  
Calculating cluster 7
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=09s  
Calculating cluster 8
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=08s  
> pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
# A tibble: 18 x 7
# Groups:   cluster [9]
       p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene    
       <dbl>     <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>   
 1 1.96e-107     0.730 0.901 0.594 2.69e-103 0       LDHB    
 2 1.61e- 82     0.922 0.436 0.11  2.20e- 78 0       CCR7    
 3 7.95e- 89     0.892 0.981 0.642 1.09e- 84 1       LTB     
 4 1.85e- 60     0.859 0.422 0.11  2.54e- 56 1       AQP3    
 5 0.            3.86  0.996 0.215 0.        2       S100A9  
 6 0.            3.80  0.975 0.121 0.        2       S100A8  
 7 0.            2.99  0.936 0.041 0.        3       CD79A   
 8 9.48e-271     2.49  0.622 0.022 1.30e-266 3       TCL1A   
 9 2.96e-189     2.12  0.985 0.24  4.06e-185 4       CCL5    
10 2.57e-158     2.05  0.587 0.059 3.52e-154 4       GZMK    
11 3.51e-184     2.30  0.975 0.134 4.82e-180 5       FCGR3A  
12 2.03e-125     2.14  1     0.315 2.78e-121 5       LST1    
13 7.95e-269     3.35  0.961 0.068 1.09e-264 6       GZMB    
14 3.13e-191     3.69  0.961 0.131 4.30e-187 6       GNLY    
15 1.48e-220     2.68  0.812 0.011 2.03e-216 7       FCER1A  
16 1.67e- 21     1.99  1     0.513 2.28e- 17 7       HLA-DPB1
17 7.73e-200     5.02  1     0.01  1.06e-195 8       PF4     
18 3.68e-110     5.94  1     0.024 5.05e-106 8       PPBP   

Seurat有几种方法可以检测差异表达，用test.use参数来定义你使用的是哪一种，比如说ROC检验，ROC检验返回的是每一个marker的 ‘classification power（分类能力）’（值从0-1，1是最好的）：

> cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s  
> head(cluster1.markers, n = 5)
      myAUC  avg_diff power pct.1 pct.2
RPS12 0.822 0.5029988 0.644 1.000 0.991
RPS27 0.822 0.5020359 0.644 0.999 0.992
RPS6  0.820 0.4673635 0.640 1.000 0.995
RPL32 0.815 0.4242773 0.630 0.999 0.995
RPS14 0.807 0.4283480 0.614 1.000 0.994

我们还可以可视化这些marker。VlnPlot （展示在不同cluster里的表达可能分布），FeaturePlot（根据tSNE和PCA结果可视化基因表达）是最常用的可视化方法。你也可以用其他一些方法，例如：RidgePlot, CellScatter, DotPlot。

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

你也可以展示原始的count值：

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

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

我们还可以用DoHeatmap来画细胞和基因的热图，在这里我们画每一个cluster里的top20的marker。如果不够20个，则画所有的marker。

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

给每一个cluster标记名称

很幸运的是，我们可以很清楚的知道每一个cluster的marker是什么，并可以为它们标记：

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