单细胞分析入门——一套经典的Seurat分析流程
前几天看到周运来老师的一句话，放在这儿，与诸君共勉：
看到才能想到，想到才能做到，做到才能得到，得到才能失去，只有失去了才知道适不适合自己
好了，切入正题，今天想给大家分享的是入门级的经典的Seurat分析流程，采用的数据为10X Genomics官方提供的外周血单个核细胞(peripheral blood mononuclear cell, PBMC)的测序数据。放一个下载链接：点击即可下载，建议电脑操作!。这里我们下载的实际上就是可以直接利用Seurat进行分析的数据。经过层层解压，你会得到这三个文件：
文件内容.jpg
当然，这三个文件是什么，可以看看我前面的文章我眼中的barcodes.tsv.gz/features.tsv.gz/matrix.mtx.gz。
下面就正式开始利用Seurat包进行分析。我个人将Seurat的上游分析分为几个部分，分别是数据质控、特征选择、降维聚类、差异基因筛选与可视化。
Seurat包目前已经更新到4.0的版本了，还好它已经托管到了CRAN下面，可以直接用install.packages()进行下载。
install.packages('Seurat')
library(Seurat)
1. 首先使用Read10X()读入数据：
#保证那三个文件位于这个工作目录下
pbmc <- Read10X(data.dir = 'C:/Users/DELL/Downloads/filtered_feature_bc_matrix')
当然对于这个函数，你?Read10X一下，你可能会发现一些你以后可能会用到的参数，例如其中的gene.column和cell.column参数，对着你的features.tsv和barcodes.tsv文件，看看它是啥意思。
2. 然后就是创建SeuratObject对象：
pbmc <- CreateSeuratObject(counts = pbmc, project = "pbmc3k", min.cells = 3, min.features = 200)
实际上在这里是有一小步的质控的，就在于min.cells和min.features两个参数上，这两个参数的意思是：过滤掉表达基因种类小于200个的细胞和被表达在不足3个细胞中的基因。
3. 对线粒体基因进行统计
看你自己项目的实际需求，比如你做单细胞核测序，一般来说就应该要做这一步，当然你也可以用同样的正则匹配去做核糖体基因等。
mt.genes <- grep(pattern = "^MT-", x = rownames(pbmc), value = T)
注意，我们做的人的数据，人的线粒体基因一般是以MT开头的，我们在pbmc这个文件的行名中匹配这个模式，找打线粒体基因，value这个参数，设置为T就是返回线粒体基因的名称，否则返回所在的行号。
对每个细胞计算其线粒体基因的含量，在这里我列出3中方法，当然也只是在看教程而已：
pbmc[['percent.mt']] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
pbmc <- PercentageFeatureSet(object = pbmc, pattern = "^MT-", col.name = "percent.mt")
pbmc <- PercentageFeatureSet(object = pbmc, features = mt.genes, col.name = "percent.mt")
你会发现，其实我们也可以自己给它一个基因的集合，让它去进行计算，特别是对于那些不太方便用正则表达式进行展示的基因集合。
4. 数据质控
这个部分需要根据你自己的数据选择合适自己的参数，这里我和官方文档设置的一样。
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nCount_RNA < 2500 & percent.mt < 5)
核心意思在于过滤后只剩下：表达基因超过200种，少于2500种和线粒体基因含量小于5%的细胞。
5. 降维聚类
pbmc <- NormalizeData(object = pbmc)
#=============findvariablefeatures==============#
pbmc <- FindVariableFeatures(object = pbmc, nfeatures = 2000)
#=================scale data=================#
pbmc <- ScaleData(object = pbmc, features = rownames(pbmc))
#=================PCA analysis===================#
pbmc <- RunPCA(object = pbmc, features = VariableFeatures(pbmc), verbose = F)
至于什么是NormalizeData和ScaleData，我在单细胞分析中的NormalizeData()与ScaleData()区别在哪儿？写过一点自己的感悟。这里为什么ScaleData()是针对于全部的基因，官方也给出了解释，确实这样会在这一步消耗更多的时间，特别是在大数据集的时候，但是如果我们在这里不对高变基因以外的基因进行Scale，就会在绘制heatmap时出现不必要的麻烦。
下一个重点在于我们选择多少个主成分进行下游的分析，Seurat给了两种图形显示让我们对我们的数据有直观的感受。
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)
ElbowPlot(pbmc)
进行UMAP和TSNE聚类分析，两者的差异在我的Literature Review(1): A short but comprehensive comparison about tSNE, UMAP and PCA这篇帖子里有过分享，当然在这里我都做了，结果会储存在SeuratObject中的不同位置，彼此并不会相互影响，这也是SeuratObject的一个突出优势。
#=============UMAP & TSNE analysis================#
pbmc <- FindNeighbors(object = pbmc, reduction = "pca", dims = 1:10)
pbmc <- FindClusters(object = pbmc, resolution = 0.5)
pbmc <- RunUMAP(object = pbmc, dims = 1:10)
pbmc <- RunTSNE(object = pbmc, dims = 1:10)
注意，这里FindClusters()函数中的resolution参数非常重要，这个值越大，聚出来的cluster也就越多。
6. 类型注释
一般来说，单细胞的上游分析，最费劲的地方就在于调整你的cluster数，聚出合适的类，进行比较准确的注释。在这里我就直接套用官方的聚类注释结果，直接进行注释了。
cluster_cell_type <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono",
"NK", "DC", "Platelet")
names(cluster_cell_type) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, cluster_cell_type)
pbmc[['cell_type']] <- pbmc@active.ident
注意，pbmc[[]]这种形式是直接在meta.data中添加元素（列），是在细胞的水平上，这在我们前面的统计线粒体基因处也用到了。当然有很多种不同的将注释结果映射到每个细胞上的方法，后面有时间再一起总结、学习！
#TSNE
DimPlot(object = pbmc,
reduction = 'tsne',
cols = c('#FF0033','#CC3399','#993366','#009933', '#0066CC','#003399','#9966CC','#99CC33','#003399'),
label = T,
label.box = T,
pt.size = 1)
#UMAP
DimPlot(object = pbmc,
reduction = 'umap',
cols = c('#FF0033','#CC3399','#993366','#009933', '#0066CC','#003399','#9966CC','#99CC33','#003399'),
label = T,
label.box = T,
pt.size = 1)
tSNE.png
UMAP.png
看来，两种不同的聚类方式绘制出的图还是很不一样的哈！最后，我想说，Seurat最好的教程是它的官网Seurat，没有谁比创造它的人更清楚它是怎么工作的了，我们不过在重复它自己给出的教程罢了。
今天也是摸鱼的一天！
