单细胞分析之Seurat分析教程（单样本）

由于Immujent最近比较忙，单细胞流程的教程就搁置了一段时间，今天小编继续拾起之前的剩下的内容来介绍一下Seurat包。很多单细胞的分析都是建立在确定了细胞分群（类型）以后，而常用来做细胞分群的软件就是Seurat，在2022年尹始，我们就一起来学习一下Seurat的分析流程吧~

1. 软件包的安装和加载

install.packages('Seurat') # 如果已安装则不需要运行
library(Seurat)
library(dplyr)

图片

2. 导入数据

这里用到的测试数据都是从Seurat官网下载的，大家有需要的也可以去官网下载（链接：https:///seurat/）或者通过本公众号获取，当然如果有自己的数据也可以用自己的数据来测试，具体得到矩阵的方法可参考单细胞分析流程之Cell Ranger和单细胞分析流程之Cell Ranger结果解读

setwd("J:/scRNA-seq/")  
#设置工作目录，这是事先创建好的目录，这样做的好处是方便进行数据整理
pbmc.data <- Read10X(data.dir = "hg19")    
#导入存放数据的目录，hg19是目录名
#如果是从GEO等网站中下载的单细胞表达矩阵，可以直接读入表达矩阵：pbmc.data <- read.table("FPKM.txt",header = TRUE)
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
# 创建seurat对象（10X数据和表达矩阵都可以通过这条命令创建对象），
# min.cells：一个基因至少在多少个细胞中表达（默认3）
# min.features：一个细胞中最少表达多少个基因（默认200）
# 一般在最开始创建对象的时候不会过滤掉太多细胞，建议使用默认参数即可
pbmc#查看pbmc中存的细胞数和基因数；另外可以用str(pbmc)查看pbmc中的包含的数据

图片

如上图可以看到，一共有13,714个基因，2,700个细胞

3. 质控并且挑选用于后续分析的细胞

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")     
# 计算线粒体基因的比例，人中线粒体的基因都是以MT开头，如果是小鼠的细胞则将MT替换成Mt即可
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3,pt.size = 0.1)  
# pt.size设置点的大小（展示的是基因数，UMI和线粒体比例的结果）` </pre>

图片

通过上面的命令可以先初步看一下数据质量，然后再挑选过滤的指标。

常用的质控信息：

a. 每个细胞检测到的基因数目的分布:基因数过低的细胞可能是捕获率低或细胞破碎或空beads,基因数过高可能是doublets或multiplets。

b. 每个细胞中总的UMI(reads)数的分布：与基因数一般是很强的正相关，拥有很少的reads或UMI分子数的样品可能是细胞破损或捕获失败，应该移除。

c. 线粒体基因表达所占的比例:一般低质量或死亡细胞线粒体污染比例高一些。

比较简单的质控方法如下图（去除异常极端值）：

图片

那么我们接下来就使用这种过滤异常极端值的方法来筛选细胞：

Q1 <- quantile(pbmc$nFeature_RNA)[2]
Q3 <- quantile(pbmc$nFeature_RNA)[4]
IQR <- Q3 - Q1
upper <- Q3+1.5*IQR
lower <- Q1-1.5*IQR
pbmc <- subset(pbmc, subset = nFeature_RNA > lower & nFeature_RNA < upper & percent.mt < 5)
pbmc  
# 查看过滤完以后查看还剩多少细胞和基因` </pre>

图片

可以看到，通过这种过滤方式一共还剩下13,714个基因，2,491个细胞。

过滤细胞这一步是需要反复调整的，主要看需要解决的科学问题是什么。如果一个课题如果想研究细胞是怎么凋亡的，那么在过滤的时候就不能把线粒体的比例作为筛选标准，因为一般凋亡的细胞都有较高的线粒体比例，一旦去掉以后则会把想重点研究的细胞都去掉，所以过滤这一步大家一定要反复测试。

4. 数据标准化

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

图片

可以把这一步理解为bulk RNA-seq中计算的TPM/100，因为每个细胞中标准化后的数据相加总和为10,000（scale.factor的值），标准化后的数据存放在pbmc[["RNA"]]@data，可以通过输入pbmc[["RNA"]]@data调用标准化后的值

如果有心的人会发现，如果直接计算pbmc[["RNA"]]@data（后面就简称为data）每一列的和，值并不是10,000，如下图：

图片

这里就是其中的一个重点啦，因为在做标准化的时候选择的标准化方法是normalization.method = "LogNormalize"，这样会在计算完'TPM'后对数据做log处理，所以直接去计算单个细胞中data的和是不等于10,000的。需要在计算之前对数据进行转换，如下命令：

TPM <- expm1(pbmc[["RNA"]]@data)
# 这样就能得到单细胞中的TPM矩阵（每一列的和为10,000）
head(colSums(TPM))

图片

5. 寻找高变基因

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# 默认情况会选择前2000个基因，nfeatures可以设置选择的基因数
top10 <- head(VariableFeatures(pbmc), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)  #标记出top10的var genes是哪些
plot1
plot2

该图的X轴是基因在所有细胞中的均值，Y轴是变异程度，Y轴越大变成程度越大，默认会根据Y轴选择前2,000个基因作为高变基因（红色的点）。

值得注意的是，这里的高变基因也可以根据课题的需要自己搜集一些关键基因用来替换这里的VariableFeatures，在后面涉及到VariableFeatures的地方替换成搜集的基因集即可（这也是根据课题的需要）

6. 数据归一化

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
# 对所有基因进行归一化，在做下一步PCA之前需要对数据进行归一化处理
# 数据缩放后的结果保存在pbmc[["RNA"]]@scale.data` </pre>

7. 线性降维（PCA）

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc),verbose = FALSE) 
# 通过设置 features来选择用于降维的基因，也可以输入自己搜集的基因集（见 5. 寻找高变基因）
DimPlot(object = pbmc, reduction = "pca")
# 绘制细胞PCA的结果` </pre>

ElbowPlot(pbmc)
# 选择用于做后续分析中的数据（PC），一般选取趋紧平缓的点，即可选择10-15 PC.
# 注意：如果选择的PC数太少，会影响数据。如果在不知道怎么选择时，可以多测试一下

8. 细胞聚类

pbmc <- FindNeighbors(pbmc, dims = 1:10)
# dims：通过上一步确定的PC，这里选择的是1~10 PC，大家在做的过程中可以多调整一下，可以尝试1:12,1:15等多个组合，看看PC的增加对结果的影响
pbmc <- FindClusters(pbmc, resolution = 0.7)  
# resolution设置分辨率，即分群的数目，需不断调试；同样的数据resolution越大，分的群越多

Number of communities：分群的个数，从上图可以看到一共分了10群，

9. 可视化（TSNE/UMAP）

pbmc <- RunTSNE(pbmc, dims = 1:10)
# dims中的参数需要和FindNeighbors中选择的PC保持一致
DimPlot(pbmc, reduction = "tsne")
# TSNE的结果

# If you haven't installed UMAP, you can do so via reticulate::py_install(packages = 'umap-learn')
pbmc <- RunUMAP(object = pbmc, dims = 1:10)
DimPlot(object = pbmc, reduction = "umap")
# UMAP的结果

从上面两个图大家可以看出，UMAP的结果会相对来说更集中一些，会把同一个群中的细胞集中在一起。当然这个也要看具体的数据，现在主要是推荐用UMAP进行可视化，不过小编觉得TSNE和UMAP都是可视化的方法，是为数据服务的，所以当然是哪个图好看就用哪个图啦 ~不用一味的追求一定要用某种方法。小编在做课题的时候就会两种都跑一下，然后看哪个图好看就会在文章中放哪个图（希望这点能帮到你

10. 寻找每个cluster的Marker gene并鉴定细胞类型

pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
# min.pct：基因在最少多大比例的细胞中检测到才用于后续的检测分析；
# logfc.threshold：高于给定的最小的变化倍数(log2)的基因用于后续统计分析。数值越大，计算会越快，差异基因可能会少；注意的是seurat v4中计算的差异倍数是log2，而seurat v3中计算的差异倍数是loge,也就是ln
# 计算每个cluster相对于其他cluster 的marker gene。
# 使用FindAllMarkers就能一步把所有cluster的marker都计算完，比较简单的方式

计算完每个群的差异基因就可以用这些差异基因去鉴定细胞类型啦，鉴定细胞类型的方法这次就不介绍了。

11. 基因表达可视化

下面就介绍一些展示基因表达的方法，可以放到文章中哦~

1. 小提琴图

library(ggplot2)``VlnPlot(object = pbmc, features = c("MS4A1")) + theme(plot.title = element_text(face = "bold.italic"))
# theme(plot.title = element_text(face = "bold.italic")) 使得基因名斜体，不过这个命令只能对单个基因绘图时有用，文章中的基因名要斜体，作图时候斜体了就不怕后面忘记修改啦~

2. 染色图

FeaturePlot(object = pbmc,features = c("MS4A1"),cols = c("#C0C0C0", "#E41C12"),reduction = "umap") + theme(plot.title = element_text(face = "bold.italic"))
# reduction：选择在pca,tsne,umap中对基因进行染色
# cols：颜色的设置
# pt.size：设置点的大小，当细胞比较稀疏时，可以将点设置较大一些

3. 气泡图

markers <- c("CD3D", "CD8A", "CD8B", "SELL", "CD4", "MS4A1", "CD79A", "MIR155HG", "NME1", "FCGR3A",  "VMO1", "CCL2", "S100A9", "HLA-DQA1", "GPR183", "PPBP", "GNG11", "TSPAN13", "IL3RA", "IGJ")
DotPlot(object = pbmc, features = markers,cols = c("#C0C0C0", "#E41C12")) + RotatedAxis
# 可用于展示每个cluster中的marker表达情况，每一列的点是可以直接比较的（对列进行scale）

图片

4. 表达分布图

RidgePlot(object = pbmc, features = c("LDHB", "FCGR3A"), ncol = 2)` </pre>

图片