标题：Seurat - Guided Clustering Tutorial
网址：https://satijalab.org/seurat/v3.2/pbmc3k_tutorial.html
部分内容转载来自：caokai001 链接：https://www.jianshu.com/p/7e4c205a576e
部分内容转载来自： 刘小泽 链接：https://www.jianshu.com/p/b46b6b6d344f
一、设置Seurat对象
对于本教程，我们将分析可从10X Genomics免费获得的外周血单个核细胞（PBMC）数据集。Illumina NextSeq 500上已对2700个单细胞进行了测序。原始数据可在此处找到。
我们从读取数据开始。该Read10X函数从10X读取cellranger管道的输出，返回唯一的分子识别（UMI）计数矩阵。该矩阵中的值表示在每个单元格（列）中检测到的每个特征（即基因；行）的分子数。
接下来，我们使用计数矩阵创建一个Seurat对象。该对象用作包含单个单元格数据集的数据（如计数矩阵）和分析（如PCA或聚类结果）的容器。有关Seurat对象结构的技术讨论，请查看我们的GitHub Wiki。例如，计数矩阵存储在中pbmc[["RNA"]]@counts。
library(dplyr)
library(Seurat)
library(patchwork)
# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "./10x_Rawdata/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
#每个基因至少在3个细胞中表达，每个细胞中至少有200个基因表达。
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)
# Lets examine a few genes in the first thirty cells
>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数据节省大量内存并节省速度。
>dense.size <- object.size(as.matrix(pbmc.data))
>dense.size
## 709591472 bytes
>sparse.size <- object.size(pbmc.data)
>sparse.size
## 29905192 bytes
>dense.size/sparse.size
## 23.7 bytes
二、标准的预处理流程
以下步骤包含Seurat中scRNA-seq数据的标准预处理工作流程。这些代表基于质量控制指标，数据归一化和缩放以及高度可变特征的检测来选择和过滤细胞。
2.1 质控和选择细胞进行进一步分析
通过Seurat，您可以轻松地探索质量控制指标并根据任何用户定义的条件过滤单元。社区常用的一些质量控制指标包括
- 每个细胞中检测到的独特基因的数量。
- 低质量的细胞或空液滴通常很少有基因
- 细胞双峰或多重峰可能显示异常高的基因计数
- 同样，细胞内检测到的分子总数（与独特基因高度相关）
- 映射到线粒体基因组的读段的百分比
- 低质量/濒死细胞通常表现出广泛的线粒体污染
- 我们使用
PercentageFeatureSet函数计算线粒体QC指标，该函数计算源自一组功能的计数的百分比 - 我们使用所有
MT-以线粒体基因开头的基因
> pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
> head(pbmc@meta.data,5)
orig.ident nCount_RNA nFeature_RNA percent.mt
AAACATACAACCAC-1 pbmc3k 2419 779 3.0177759
AAACATTGAGCTAC-1 pbmc3k 4903 1352 3.7935958
AAACATTGATCAGC-1 pbmc3k 3147 1129 0.8897363
AAACCGTGCTTCCG-1 pbmc3k 2639 960 1.7430845
AAACCGTGTATGCG-1 pbmc3k 980 521 1.2244898
在下面的示例中，我们可视化QC指标，并使用它们来过滤单元格。
- 过滤 unique feature counts 超过2500或少于200的细胞
- 过滤线粒体 > 5％的细胞
> VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
#CombinePlots(plots = list(plot1,plot2))
过滤
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
>pbmc
An object of class Seurat
13714 features across 2638 samples within 1 assay
Active assay: RNA (13714 features, 0 variable features)
之前是13714 features across 2700 samples ，过滤掉了62个细胞
三、 预处理之归一化
关于Seurat归一化原理，可以看这一篇：https://www.biorxiv.org/content/biorxiv/early/2019/03/18/576827.full.pdf
标准化是一个比较简单的过程，使用的是"logNormalize", 就是将每个基因的表达量对该细胞总表达量进行平衡，然后乘以一个因子(scale factor, 默认值为10,000), 然后中进行对数转换。
默认是进行一个全局的LogNormalize操作：log1p(value/colSums[cell-idx] *scale_factor) ，其中log1p指的是log(x + 1) ，得到的结果存储在：pbmc[["RNA"]]@data
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
四、鉴定差异基因HVGs（highly variable features）
我们在Seurat3的流程详细描述在这里，并通过直接建模在单细胞数据固有的均方差关系改进了以前的版本，并在实现FindVariableFeatures功能。V3默认选择2000个差异基因。这些将用于下游分析，例如PCA。
如果手动计算，相当于计算变异系数，而不是标准差。因为标准差对于数据呈现对称分布，及单峰情况体现差异比较好，但是对于双峰情况，变异系数cv 可能更好。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
SLM <<- FindVariableFeatures(SLM, selection.method = "mvp", mean.cutoff = c(0.1,100),dispersion.cutoff=c(-1,Inf))
##"mvp"这种写法，最后鉴定出来了4000多个差异基因，不受默认值影响。
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)
- V3横坐标范围设定参数改成：mean.cutoff；纵坐标改成：dispersion.cutoff
- 鉴定差异基因的算法包含三种：vst(默认)、mean.var.plot、dispersion
- vst：首先利用loess对 log(variance) 和log(mean) 拟合一条直线，然后利用观测均值和期望方差对基因表达量进行标准化。最后根据保留最大的标准化的表达量计算方差
- mean.var.plot: 首先利用mean.function和 dispersion.function分别计算每个基因的平均表达量和离散情况；然后根据平均表达量将基因们分散到一定数量（默认是20个）的小区间（bin）中，并且计算每个bin中z-score
- dispersion：挑选最高离散值的基因
- V3默认选择2000个差异基因，检查方法也不同（V3用VariableFeatures(sce)检查，V2用sce@var.genes检查）
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
> plot1 + plot2
Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto), :
Viewport has zero dimension(s)
此外: Warning messages:
1: Transformation introduced infinite values in continuous x-axis
2: Transformation introduced infinite values in continuous x-axis
有问题，未实现图片：
红色的HVGs 有2000个，从上面 FindVariableFeatures ，命令参数就可以知道.
五、标准化
从上图中，我们可以看到很多基因都表现出很大的差异性，但是这些差异性并不一定都是事实，有些可能是因为实验，或者数据采集的方式导入的。因此人们想尽量的消除这些因素带来的差异。这些差异有可能是因为batch-effect， 或者cell alignment rate，或者其它。这些因素都可以写入metadata中，然后通过ScaleData函数来消除。
ScaleData 可以通过metadata信息，消除很多其他因素影响。
5.1 缩放的标准
应用线性变换缩放数据，这是一个标准的预处理步骤，比PCA等降维技术更重要。使用ScaleData函数:
（1）使每个基因在所有细胞间的表达量均值为0
（2）使每个基因在所有细胞间的表达量方差为1
缩放操作给予下游分析同等的权重，这样高表达基因就不会占主导地位，其结果存储在pbmc[["RNA"]]@scale.data中，这里默认ScaleData 进行了中心化，标准化。
## 全局缩放
all.genes <- rownames(pbmc)
pbmc <- ScaleData(object = pbmc, features = all.genes)
> head(pbmc[["RNA"]]@scale.data[, 1:3])
## AAACATACAACCAC AAACATTGAGCTAC AAACCGTGCTTCCG
## AL627309.1 -0.06452497 -0.06452497 -0.06452497
## AP006222.2 -0.03726409 -0.03726409 -0.03726409
## RP11-206L10.2 -0.04153370 -0.04153370 -0.04153370
## RP11-206L10.9 -0.03734170 -0.03734170 -0.03734170
## LINC00115 -0.08343360 -0.08343360 -0.08343360
## NOC2L -0.32126699 -0.32126699 -0.32126699
缩放是Seurat工作流程中必不可少的步骤，但仅限于将用作PCA输入的基因。因此，默认设置ScaleData是仅对先前确定的变量特征执行缩放（默认为2,000）。为此，请忽略features上一个函数调用中的参数，即
>pbmc <- ScaleData(pbmc) #默认为2,000个高边基因
这样操作的结果中降维和聚类不会被影响，但是只对HVGs基因进行标准化，下面画的热图依然会受到全部基因中某些极值的影响。所以为了热图的结果，还是对所有基因进行标准化比较好。
ScaleData的操作过程是怎样的?
看一下帮助文档就能了解：
ScaleData这个函数有两个默认参数：do.scale = TRUE和do.center = TRUE，然后需要输入进行scale/center的基因名（默认是取前面FindVariableFeatures的结果）。
scale和center这两个默认参数应该不陌生了，center就是对每个基因在不同细胞的表达量都减去均值；
scale就是对每个进行center后的值再除以标准差（就是进行了一个z-score的操作）
再看一个来自BioStar的讨论：https://www.biostars.org/p/349824/
问题1：为什么在NormalizeData()之后，还需要进行ScaleData操作？
前面NormalizeData进行的归一化主要考虑了测序深度/文库大小的影响（reads *10000 divide by total reads, and then log），但实际上归一化的结果还是存在基因表达量分布区间不同的问题；而ScaleData做的就是在归一化的基础上再添zero-centres （mean/sd =》 z-score），将各个表达量放在了同一个范围中，真正实现了”表达量的同一起跑线“。
另外，新版的ScaleData函数还支持设定回归参数，如（nUMI、percent.mito），来校正一些技术噪音
问题2： FindVariableGenes() or RunPCA() or FindCluster()这些参数是基于归一化Normalized_Data还是标准化Scaled_Data ？
所有操作都是基于标准化的数据Scaled_Data ，因为这个数据是针对基因间比较的。例如有两组基因表达量如下:
g1 10 20 30 40 50
g2 20 40 60 80 100
虽然它们看起来是g2的表达量是g1的两倍，但真要降维聚类时，就要看scale的结果：
> scale(c(10,20,30,40,50))
[,1]
[1,] -1.2649111
[2,] -0.6324555
[3,] 0.0000000
[4,] 0.6324555
[5,] 1.2649111
attr(,"scaled:center")
[1] 30
attr(,"scaled:scale")
[1] 15.81139
> scale(c(20,40,60,80,100))
[,1]
[1,] -1.2649111
[2,] -0.6324555
[3,] 0.0000000
[4,] 0.6324555
[5,] 1.2649111
attr(,"scaled:center")
[1] 60
attr(,"scaled:scale")
[1] 31.62278
二者标准化后的分布形态是一样的（只是中位数不同），而我们后来聚类也是看数据的分布，分布相似聚在一起。所以归一化差两倍的两个基因，根据标准化的结果依然聚在了一起
问题3： ScaleData()在scRNA分析中一定要用吗？
实际上标准化这个过程不是单细胞分析固有的，在很多机器学习、降维算法中比较不同feature的距离时经常使用。当不进行标准化时，有时feature之间巨大的差异（就像👆问题2中差两倍的基因表达量，实际上可能差的倍数会更多）会影响分析结果，让我们误以为它们的数据分布相差很远。
很多时候，我们会混淆normalization和scale：那么看一句英文解释：
Normalization "normalizes" within the cell for the difference in sequenicng depth / mRNA throughput. 主要着眼于样本的文库大小差异
Scaling "normalizes" across the sample for differences in range of variation of expression of genes . 主要着眼于基因的表达分布差异
另外，看scale()函数的帮助文档就会发现，它实际上是对列进行操作：
那么具体操作中是要先对表达矩阵转置，然后scale，最后转回来
t(scale(t(dat[genes,])))
问题4：RunCCA()函数需要归一化数据还是标准化数据？
通过看这个函数的帮助文档：
RunCCA(object, object2, group1, group2, group.by, num.cc = 20, genes.use,
scale.data = TRUE, rescale.groups = FALSE, ...)
显示scale.data = TRUE，那么就需要标准化后的数据
删除不需要的数据
使用ScaleData函数从单细胞数据集中删除不需要的变量。例如，我们可以“regress out”与细胞周期阶段或线粒体污染相关的异质性（去除不需要的变量，即校正协变量，校正线粒体基因比例的影响）。如下：
# 删除不需要的数据，校正线粒体基因引起的影响
# 这个步骤非常耗时
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
我们对矩阵进行了细胞水平质控，但是如何对基因进行质控，当然我们创建seurat 对象时候，已经设定一个cutoff. 基因必须早多少个细胞中表达. 但是如果直接用这些基因进行聚类和降维（T-sne/UMAP 可视化 可能引起偏差，所以需要先进行一下PCA 初步选取主要的主成分，对应着权重不同的基因，再用这些基因进行聚类，降维可视化效果更好）
六、PCA
接下来，我们对缩放后的数据执行PCA。默认情况下，仅将先前确定的变量特征用作输入，但是features如果您希望选择其他子集，则可以使用arguments进行定义。默认使用标准化的数据
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
# 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
Seurat提供了多种手段来查看PCA分析的结果，其中包括PrintPCA, VizPCA, PCAPlot, 以及PCHeatmap
【一维可视化】 对每一个维度，权重较大的基因进行可视化
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
【二维可视化】 将细胞降维到二维空间
DimPlot(object = pbmc, reduction = "pca")
这个图中每个点就是一个样本，它根据PCA的结果坐标进行画图，这个坐标保存在了cell.embeddings中：
> head(pbmc[['pca']]@cell.embeddings)[1:2,1:2]
PC_1 PC_2
AAACATACAACCAC -4.7296855 -0.5184265
AAACATTGAGCTAC -0.5174029 4.5918957
其实还支持定义分组：根据pbmc@meta.data中的ident分组：
> table(pbmc@meta.data$orig.ident)
pbmc3k
2638
# 当然这里只有一个分组
Tips: 下图以PC_1 来看，是解释细胞分布最大的基因组合方向，所有的基因是维度，有些基因对这个方向共享很大，也就是说更容易分开细胞，这些基因可能更加重要。下面是热图1展示这些基因对于区分细胞的效果，可以看到将很多细胞区分成紫色和黄色.效果不错.
可以看到5,6维度貌似区分不开了。
探索异质性来源—DimHeatmap
每个细胞和基因都根据PCA结果得分进行了排序，默认画前30个基因（nfeatures设置），1个主成分(dims设置)，细胞数没有默认值（cells设置）
> DimHeatmap(pbmc, dims = 1:12, cells = 500, balanced = TRUE)
# 其中balanced表示正负得分的基因各占一半
降维的真实目的是尽可能减少具有相关性的变量数目，例如原来有700个样本，也就是700个维度/变量，但实际上根据样本中的基因表达量来归类，它们或多或少都会有一些关联，这样有关联的一些样本就会聚成一小撮，表示它们的”特征距离“相近。最后直接分析这些”小撮“，就用最少的变量代表了实际最多的特征
既然每个主成分都表示具有相关性的一小撮变量（称之为”metafeature“，元特征集），那么问题来了：
降维后怎么选择合适数量的主成分来代表整个数据集？
七、确定PCA的维数
Suerat使用了jackStraw方法[4]来估计应该使用多少个principal components的基因来进行下游分析。
# 注意: 这可能是一个非常费时的过程。
pbmc <- JackStraw(pbmc, num.replicate = 100) # 重复一百次
pbmc <- ScoreJackStraw(pbmc, dims = 1:20) # 选择20个PC
JackStrawPlot以及PCElbowPlot函数提供给我们两种可视化手段来查看jackStraw的结果。使用JackStrawPlot可以查看每个PC的p value，这样方便我们通过p值进行筛选。PCElbowPlot可以方便的查看从哪个PC开始，差分度就进入的平台期。
- 1.
JackStrawPlot
JackStrawPlot(pbmc, dims = 1:15)
-
2.
PCElbowPlot: 拐点图 # 最大可设置到50ElbowPlot(pbmc,ndims = 50, reduction = "pca") # 最大可设置到50
它是根据每个主成分对总体变异水平的贡献百分比排序得到的图，我们主要关注”肘部“的PC，它是一个转折点（也即是这里的PC9-10），说明取前10个主成分可以比较好地代表总体变化
但官方依然建议我们，下游分析时多用几个成分试试（10个、15个甚至50个），如果起初选择的合适，结果不会有太大变化；另外推荐开始不确定时要多选一些主成分，也不要直接就定下5个成分进行后续分析，那样很有可能会出错
先择PC，对于得到正确的结果具体关键性的作用。除了上面的方法以外，还可以使用其它的一些办法，比如说使用GSEA等富集分析的办法来选择PC，或者说使用一个随机模型来查看具体哪些PC会比较有效果。但是后者这一方法在实现起来需要消耗过多的资源。
在本例中，因为是Seurat挑选的例子，所以通过上面的JackStraw方法，只要把cut.off值设置在7-10之间，就可以得到差不多的结果。
Tips: 确定进行聚类和降维的基因和细胞后，下一步就是聚类了，在这里需要搞情况聚类和降维作用是不一样的。特别是T-sne是降维可视化作用，而不是用它来聚类.为什么看到T-sne 图是不同颜色的色块呢？常见的聚类包括kmeans/hclust .
八、细胞聚类
当决定了使用哪些PC中的基因对细胞进行分类之后，就可以使用FindClusters来对细胞聚类了。分群后的结果会保存在object@ident slot中。
Seurat 3版本提供了graph-based 聚类算法，参考两篇文章：[SNN-Cliq, Xu and Su, Bioinformatics, 2015]和 [PhenoGraph, Levine et al., Cell, 2015]. 简言之，这些graph的方法都是将细胞嵌入一个算法图层中，例如：K-nearest neighbor (KNN) graph 中，每个细胞之间的连线就表示相似的基因表达模式，然后按照这个相似性将图分隔成一个个的高度相关的小群体（名叫‘quasi-cliques’ or ‘communities’）；
（Step-1）在PhenoGraph中，首先根据PCA的欧氏距离结果，构建了KNN图，找到了各个近邻组成的小社区；然后根据近邻中两两住户（细胞）之间的交往次数（shared overlap）计算每条线的权重（术语叫Jaccard similarity） 。这些计算都是用包装好的函数FindNeighbors() 得到的，它的输入就是前面降维最终确定的主成分。
（Step-2）得到权重后，为了对细胞进行聚类，使用了计算模块的算法（默认使用Louvain，另外还有包括SLM的其他三种），使用FindClusters() 进行聚类。其中包含了一个resolution的选项，它会设置一个”间隔“值，该值越大，间隔越大，得到的cluster数量越多。一般来说，这个值在细胞数量为3000左右时设为0.4-1.2 会有比较好的结果
pbmc <- FindNeighbors(pbmc, dims = 1:10) ## dims 作为输入的维度
pbmc <- FindClusters(pbmc, resolution = 0.6) ## 分辨率参数，和cluster 个数有关系，如何你想获得粗糙的分类，需要调成大一些.
# 查看一下前几个的cluster,结果聚成3类.可能发现新的细胞系.
head(Idents(pbmc), 5)
## AAACATACAACCAC AAACATTGAGCTAC AAACCGTGCTTCCG AAACCGTGTATGCG AAACGCTGGTTCTT
## 4 3 0 6 4
## Levels: 0 1 2 3 4 5 6 7 8
九、降维后，用聚类的类别可视化【UMAP/tSNE】
使用UMAP或者tSNE的好处在于将多维的PC图降维至2维图像中来，这样方便分集细胞的可视化。 UMAP皱缩成团分不开，而tSNE图则分散得多。
pbmc <- RunUMAP(pbmc, dims = 1:10)
P1 <- DimPlot(pbmc, reduction = "umap", label = TRUE)
pbmc <- RunTSNE(pbmc, dims = 1:10)
P2 <- DimPlot(pbmc, reduction = "tsne",label = TRUE)
P1 + P2
您可以在此时保存该对象，以便可以轻松地将其装入，而不必重新运行上面执行的计算密集型步骤，也可以轻松地与协作者共享该对象。
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
得到聚类分组结果，通常需要继续分析，在每组里面那些标记基因。由于10x genomic 对每一个细胞都进行标签化，不知道每一个标签对应那种细胞，可以通过每组的标记基因对每组进行细胞类型判别（根据文献说明的已知marker）,进而发现新的标记基因.
十、寻找分群marker基因
生成了tSNE图像后，我们还需要确认这一堆堆的细胞都分别是些什么细胞，这时，可以使用特定细胞的marker来对不同的细胞群进行标记。但是想要对它们进行标记，首先需要得到每个集群差异表达的基因。
在比较时，Seurat使用一比多的办法。如果需要对cluster 1的细胞(ident.1)进行差异分析时，它比较的对象将是除了cluster 1以外的其它所有的细胞。也可以使用ident.2参数来传递需要比较的对象。 在FindMarkers函数中，min.pct是指的marker基因所占细胞数的最低比例，这里使用1/4。 使用max.cells.per.ident可以让函数对细胞进行抽样，以加速计算。 FindAllMarkers可以一次性找到所有的高表达的marker。
FindAllMarkers()函数，它默认会对单独的一个细胞群与其他几群细胞进行比较找到正、负表达biomarker（这里的正负有点上调、下调基因的意思；正marker表示在我这个cluster中表达量高，而在其他的cluster中低）。
# find all markers of cluster 1
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
## p_val avg_logFC pct.1 pct.2 p_val_adj
## RPS12 4.371171e-118 0.5264969 1 0.990 5.994624e-114
## RPS6 1.681008e-114 0.5234111 1 0.994 2.305334e-110
## RPS27 1.719993e-107 0.5242477 1 0.991 2.358799e-103
## RPS14 1.546910e-105 0.4500125 1 0.993 2.121432e-101
## RPS25 5.848778e-98 0.5427189 1 0.973 8.021015e-94
上面这个代码中min.pct的意思是：一个基因在任何两群细胞中的占比最低不能低于多少，当然这个可以设为0，但会大大加大计算量.
还可以计算cluster 5与（cluster0+cluster3）的差异基因
# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(object = pbmc, ident.1 = 5,
ident.2 = c(0,3), min.pct = 0.25)
head(cluster5.markers, n = 5)
## p_val avg_logFC pct.1 pct.2 p_val_adj
## FCGR3A 3.849393e-139 2.723391 0.981 0.098 5.279058e-135
## RHOC 1.149445e-91 1.665090 0.865 0.130 1.576350e-87
## CDKN1C 6.799344e-75 1.061123 0.519 0.022 9.324621e-71
## IFITM2 1.347125e-73 1.560680 1.000 0.644 1.847448e-69
## HES4 3.616785e-69 1.107357 0.596 0.052 4.960059e-65
计算所有类和其他类的标记基因,（对所有cluster都比较一下，only.pos = TRUE为只挑出来正表达marker）：
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE,
min.pct = 0.25, logfc.threshold = 0.25)
#挑每个cluster的前2个
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 0\. 3.72 0.969 0.136 0\. 0 S100A8
## 2 6.30e-298 3.77 0.996 0.236 8.63e-294 0 S100A9
## 3 7.44e- 94 0.807 0.938 0.587 1.02e- 89 1 LDHB
## 4 3.14e- 61 0.848 0.426 0.105 4.31e- 57 1 CCR7
## 5 1.22e- 80 0.975 0.981 0.635 1.67e- 76 2 LTB
## 6 5.02e- 69 0.956 0.783 0.315 6.88e- 65 2 IL7R
## 7 0\. 2.95 0.931 0.045 0\. 3 CD79A
## 8 5.62e-214 2.46 0.617 0.024 7.71e-210 3 TCL1A
## 9 3.86e-164 2.39 0.973 0.22 5.30e-160 4 CCL5
## 10 2.97e-142 2.12 0.595 0.052 4.07e-138 4 GZMK
## 11 1.66e-171 2.34 0.981 0.135 2.27e-167 5 FCGR3A
## 12 1.68e-106 1.90 1 0.368 2.30e-102 5 LST1
## 13 2.24e-211 3.48 0.972 0.063 3.07e-207 6 GZMB
## 14 1.84e-127 3.52 0.943 0.138 2.52e-123 6 GNLY
## 15 1.11e-180 2.68 0.806 0.013 1.52e-176 7 FCER1A
## 16 5.01e- 20 1.93 1 0.548 6.87e- 16 7 HLA-DPB1
## 17 4.08e-171 5.02 1 0.012 5.60e-167 8 PF4
## 18 6.09e- 99 5.96 1 0.026 8.35e- 95 8 PPBP
>top2 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
> top2$gene
[1] "LDHB" "CCR7" "LTB" "AQP3" "S100A9" "S100A8" "CD79A" "TCL1A"
[9] "CCL5" "GZMK" "FCGR3A" "LST1" "GZMB" "GNLY" "TBCB" "NDUFA2"
[17] "FCER1A" "HLA-DPB1" "PF4" "PPBP"
还可以使用不同的算法来寻找差异表达的基因。在test.use选项中设置（详见：DE vignette ）例如：
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 0,
thresh.use = 0.25, test.use = "roc",
only.pos = TRUE)
想了解更多关于查找差异表达基因的信息，可以查看Seurat的vignette。
【一维可视化】标记基因在不同cluster 表达量
拿到差异表达基因后，可以使用VlnPlot来查看结果。
VlnPlot(object = pbmc, features = c("MS4A1", "CD79A"))
# you can plot raw UMI counts as well
VlnPlot(pbmc, features = c("MS4A1", "CD79A"), slot = "counts", log = TRUE)
## RidgePlot
## RidgePlot
RidgePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
## DotPlot
DotPlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
【二维可视化】标记基因在不同细胞表达量
FeaturePlot(object = pbmc,
features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A",
"FCGR3A", "LYZ", "PPBP", "CD8A"),
cols = c("blue", "red"))
提取所有cluster 的top10 标记基因可视化
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
# 设置slim.col.label=TRUE将避免打印每个细胞的名称，而只打印cluster的ID.
DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()
目前我们通过聚类知道了，将细胞分成8个类别，但是并不知道具体是那种细胞。所有需要进一步通过已知的特定细胞标记基因对cluster 打上细胞类型标签.
十一、对分群进行细胞类型标记
对于示例数据，我们已经有一些比较常见的细胞marker了，很方便的就可以对分集的细胞进行标记了。
| Cluster ID | Markers | Cell Type |
|---|---|---|
| 0 | IL7R, CCR7 | Naive CD4+ T |
| 1 | IL7R, S100A4 | Memory CD4+ |
| 2 | CD14, LYZ | CD14+ Mono |
| 3 | MS4A1 | B |
| 4 | CD8A | CD8+ T |
| 5 | FCGR3A, MS4A7 | FCGR3A+ Mono |
| 6 | GNLY, NKG7 | NK |
| 7 | FCER1A, CST3 | DC |
| 8 | PPBP | Platelet |
将原来的cluster 0-9 ，修改成对应的细胞名称
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T",
"CD14+ Mono", "B", "CD8 T",
"FCGR3A+ Mono",
"NK", "DC", "Platelet","NULL")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
P1 <- DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
P2 <- DimPlot(pbmc, reduction = "tsne", label = TRUE, pt.size = 0.5) + NoLegend()
P1 + P2
saveRDS(pbmc, file = "./output/pbmc3k_final.rds")
当我们修改了聚类的参数，结果会不同。比如分辨率参数.
# 生成一个快照, 其实就是把pbmc@ident中的数据拷贝到pbmc@meta.data中。
# 当我们想把pbmc@meta.data中的值赋值给pbmc@ident的话，可以使用SetAllIdent函数。
# 比如pbmc <- SetAllIdent(object=pbmc, id = 'orig.ident')
pbmc <- StashIdent(object = pbmc, save.name = "ClusterNames_0.6")
# 现在我们使用不同的resolution来重新分群。当我们提高resolution的时候，
# 我们可以得到更多的群。
# 之前我们在FindClusters中设置了save.snn=TRUE，所以可以不用再计算SNN了。
# 下面就直接提高一下区分度，设置resolution = 0.8
pbmc <- FindClusters(pbmc, resolution = 0.8)
番外
分辨率左图是0.8，右图0.6
# 并排画两个tSNE,右边的是用ClusterNames_0.6 (resolution=0.6) 来分组颜色的。
# 我们可以看到当resolution提高之后，CD4 T细胞被分成了两个亚群。
plot1 <- DimPlot(object = pbmc, reduction = "umap", label = TRUE) + NoLegend()
plot2 <- DimPlot(object = pbmc, reduction = "umap", label = TRUE,
group.by = "ClusterNames_0.6") +NoLegend()
plot1 + plot2
类别不同，我们又可以按照这一次聚类结果，寻找marker基因
# 我们再一次寻找不同分集中的marders。
tcell.markers <- FindMarkers(object = pbmc, ident.1 = 0, ident.2 = 1)
# Most of the markers tend to be expressed in C1 (i.e. S100A4).
# However, we can see that CCR7 is upregulated in
# C0, strongly indicating that we can differentiate memory from naive CD4 cells.
# cols.use demarcates the color palette from low to high expression
FeaturePlot(object = pbmc, features = c("S100A4", "CCR7"),
cols = c("blue", "red"))
