2018年8月份的时候，我也使用过seurat来分析单细胞测序数据，然后最近也需要使用seurat包来分析实验室的单细胞测序数据，在R中安装完seurat的包后，我到网站上下载了pbmc3k_tutorial.Rmd的指导文件，下载下来按着这个markdown文件执行，发现好像跟我暑期的代码不大一样啊，然后我又返回到网站上重新看了下，发现原来是“On April 16, 2019 - we officially updated the Seurat CRAN repository to release 3.0!” 在2019年4月16的时候更新版本啦。好吧，R包更新的话，那么之前的代码就是以前的版本了，如果使用新的版本，那么代码就不能用，但是呢，原理是一样的，而且发现新版本更好用了，自己敲代码都少了哈哈哈。

一、总的分析流程

image.png

Note:自己的一个理解但是可能并不准确

二、具体的步骤

使用的数据是Seurat提供的pbmc_3k的10x单细胞数据，共有32738基因和2700x细胞

1.Setup the Seurat Object（创建Seurat对象）

library(dplyr)
library(Seurat)
rm(list = ls())
pbmc.data <- Read10X(data.dir = "~/seurat/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) # gene:13714,cells：2700
#比较稀疏矩阵和普通矩阵使用的内存大小
dense.size <- object.size(x = as.matrix(x = pbmc.data))
dense.size
sparse.size <- object.size(x = pbmc.data)
sparse.size
dense.size / sparse.size

1）先是用Read10x读入数据，10xGenomics得到的数据是用稀疏矩阵存放的，因为单细胞数据又很多0，也就是很多未检测到的数据，稀疏矩阵存放的话可以节省内存
2）接着用CreateSeuratObject创建Seurat对象，min.cells = 3, min.features = 200筛选的条件是min.cells = 3基因至少在3个细胞中表达， min.features=200，每个细胞中至少表达200个基因

2.Standard pre-processing workflow（数据预处理流程）

可以分为四个部分：

• Selection and filtration of the cells based on the QC metrics 基于QC的结果选择和过滤细胞

• Data normalization 数据标准化

• scaling 量化(中心化)

• Detection of the highly variable feature 高变异基因的检测

（1）Selection and filtration of the cells based on the QC metrics 基于QC的结果选择和过滤细胞

1)QC的标准

• The number of unique genes detected in each cell 每个细胞中的唯一基因的数目
• Similarly, the total number of molecules detected within a cell (correlates strongly with unique genes) 每个细胞中的UMI数目
• The percentage of reads that map to the mitochondrial genome 比对到线粒体基因中的read数的百分比（因为低质量/将死的细胞通常表现出广泛的线粒体污染）

2)代码实现

pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-") #线粒体基因以MT开头的选择出来
##nFeature_RNA就是每个细胞检测到的基因数目；nCount_RNA就是这些基因数目一共测到的count数目,也就是总的UMI数目；percent.mt（使用'PercentageFeatureSet`函数计算线粒体QC指标，该函数计算源自一组特征的计数百分比）
# Show QC metrics for the first 5 cells
head(x = pbmc@meta.data, 5)

#Visualize QC metrics as a violin plot 用小提琴图可视化QC特征
VlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3) 

# 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.
# FeatureScatter通常用于可视化特征 - 特征关系

plot1 <- FeatureScatter(object = pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt") 
plot2 <- FeatureScatter(object = pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") 
CombinePlots(plots = list(plot1,plot2))
#选择gene数目小于2500 以及大于200 以及 线粒体数目小于5%的细胞 #可以基于violin图来选择你要卡的阈值
#过滤完之后还剩余 genes:13714,cells：2638
pbmc <- subset(x = pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

3)结果

image.png
image.png

4)解释

对于小提琴图的几个参数的解释：nFeature_RNA就是每个细胞检测到的基因数目，就是以前版本的nGene；nCount_RNA就是这些基因数目一共测到的count数目,也就是以前版本的UMI数目；percent.mt（使用'PercentageFeatureSet`函数计算线粒体QC指标，该函数计算源自一组特征的计数百分比）。
对于QC标准的选择：根据小提琴图中的总的基因数目的分布以及线粒体所占的百分比，比如这边大多数基因落在200-2500之间，线粒体所占百分比小于5%这些是要保留下来的。

（2）数据标准化

1)对于数据标准化的理解

其实原先我已经忘了数据标准化是干啥的了，我搞不明白为什么要用Log转化，与我们的FPKM有什么区别，之后到网上去看了下，发现数据标准化消除数据量纲之间的不同；消除数量级之间差别很大；避免数值问题：太大的数会引发数值问题；平衡各特征的贡献；一些模型求解的需要：加快了梯度下降求最优解的速度等需要
Seurat: 从数据集中删除不需要的cell后，下一步是标准化数据。 默认情况下，Seurat采用全局标准化方法“LogNormalize”，通过总表达式对每个cell的表达值进行标准化，将其乘以比例因子（默认为10,000），并对结果进行对数转换。 归一化值存储在pbmc [[“RNA”]] @ data中。

2)代码实现

#使用log转化，度量单位是10000
pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize", scale.factor = 1e4)

（3）Identification of highly variable features (feature selection) 识别高度变异基因(特征选择的过程)

1）解释

seurat: a subset of features that exhibit high cell-to-cell variation in the dataset 在数据集中表现出细胞间高变异的特征子集

2）代码实现

pbmc <- FindVariableFeatures(object = pbmc,selection.method = 'vst', nfeatures = 2000)
# Identify the 10 most highly variable genes 取前十个变化最大的基因
top10 <- head(x = VariableFeatures(object = pbmc), 10)
# plot variable features with and without labels 画出2000个高变异基因
plot1 <- VariableFeaturePlot(object = pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
CombinePlots(plots = list(plot1, plot2))

3) 结果展示

高变异基因

4）解释

q其实我不是很理解这边的高变异基因和差异基因的区别，变异基因是说在每个细胞之间变异比较大的基因吗？其实不是很理解，这边是什么意思，而且我觉得变异最大的基因会不会是系统误差导致的呢，并不是真正的生物学变异。我对这步选择高变异基因认为是在做特征选择的过程，选出变异比较大的基因用于后续分析。

（4）Scaling 中心化

1)理解

应用线性变换（'缩放'），这是在PCA等降维技术之前的标准预处理步骤(目的:所谓去中心化，就是将样本X中的每个观测值都减掉样本均值，这样做的好处是能够使得求解协方差矩阵变得更容易。)。 ScaleData函数
该步骤在下游分析中给予相同的权重，因此高表达的基因不占优势

2)代码

all.genes <- rownames(x = pbmc)
pbmc <- ScaleData(object = pbmc, features = all.genes)
#- `ScaleData` function to remove unwanted sources of variation from a single-cell dataset, 比如mitochondrial contamination线粒体污染
pbmc <- ScaleData(object = pbmc, vars.to.regress = 'percent.mt')

3)结果展示

image.png

3.Perform linear dimensional reduction线性降维 + Cluster the cells 对细胞进行聚类分析

（1） 线性降维：主成分分析

1）代码实现

#接下来，我们对中心化的数据执行PCA。 默认情况下，只有先前确定的变量要素用作输入，但如果要选择其他子集，则可以使用`features`参数定义。
pbmc <- RunPCA(object = pbmc, features = VariableFeatures(object = pbmc))

# Examine and visualize PCA results a few different ways
print(x = pbmc[['pca']], dims = 1:5, nfeatures = 5)
VizDimLoadings(object = pbmc, dims = 1:2, reduction = 'pca')
DimPlot(object = pbmc, reduction = 'pca') #每个细胞在PC1 和 PC2中的位置图
# Doheatmap 画热图
DimHeatmap(object = pbmc, dims = 1, cells = 500, balanced = TRUE) #可以选择所要呈现的细胞数目
DimHeatmap(object = pbmc, dims = 1:3, cells = 500, balanced = TRUE) 
# Determine the 'dimensionality' of the dataset 决定选取多少主成分
# 使用jackStraw方法 来估计多少主成分比较合适
# We identify 'significant' PCs as those who have a strong enrichment of low p-value features. 显著的主成分是强烈富集 低P值的特征
pbmc <- JackStraw(object = pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(object = pbmc, dims = 1:20)
# JackStrawplot
JackStrawPlot(object = pbmc, dims = 1:15) #JackStrawPlot功能提供了一个可视化工具，用于比较每个PC的p值分布和均匀分布
# 碎石图
ElbowPlot(object = pbmc) #基于每一个解释的方差百分比（“ElbowPlot”函数）对主成分进行排序

2）结果展示

VizDimLoadings
DimPlot
heatmap
heatmap
JackStrawPlot
image.png

3）PCA中RunPCA时遇到的坑

pbmc <- RunPCA(object = pbmc, features = VariableFeatures(object = pbmc))
Error in irlba(A = t(x = object), nv = npcs, ...) : 
  max(nu, nv) must be strictly less than min(nrow(A), ncol(A))

image.png

心中纳闷啊，为啥报错啊,然后打开RunPCA的帮助文档

RunPCA

发现是npcs的选择是默认为50的。而我的矩阵是2000基因*24细胞的矩阵。我心中纳闷，我不是有2000个基因可以供选择吗？怎么我选取个主成分等于50就不行嘛？后来明白了PC的选择是小于你给的矩阵的行数或列数中的最小值的。因为PC计算过程是要计算奇异矩阵的，然后奇异矩阵是个对称矩阵哦。

image.png

(2)Cluster the cells 对细胞进行聚类分析

1)理解

• 在这里是对数据降维后在聚类，这里使用的是KNN聚类方法

• 与PhenoGraph一样，我们首先根据PCA空间中的欧氏距离构建KNN图，并根据其局部邻域中的共享重叠（Jaccard相似性）细化任意两个单元之间的边权重。 此步骤使用FindNeighbors函数执行，并将先前定义的数据集维度（前10个PC）作为输入。
• FindClusters函数实现了此过程，并包含一个分辨率参数，用于设置下游群集的“簇数目”，增加的值会导致更多的簇。
• 对于大约3000数目的单细胞数据集，我们发现将此参数设置在0.4-1.2之间结果较为良好。 较大的数据集通常会增加最佳分辨率。 可以使用Idents函数找到簇。

2)代码

pbmc <- FindNeighbors(object = pbmc, dims = 1:10)
pbmc <- FindClusters(object = pbmc, resolution = 0.5)
# Look at cluster IDs of the first 5 cells
head(x = Idents(object = pbmc), 5)

3)结果

The first 5 cells

4.Run non-linear dimensional reduction (UMAP/tSNE) 非线性降维： tSNE and UMAP

1）理解

• The goal of these algorithms is to learn the underlying manifold of the data in order to place similar cells together in low-dimensional space.
  -算法的目的：了解数据的多样性，以便在低维空间中将相似的细胞放在一起

• 作为UMAP和tSNE的输入，我们建议使用相同的PC作为聚类分析的输入。
• 降维有助于数据可视化
• UMAP（Uniform Manifold Approximation and Projection）是一种用于降维的新型流形学习技术。 UMAP由基于黎曼几何和代数拓扑的理论框架构成。 结果是适用于现实世界数据的实用可扩展算法。 UMAP算法与t-SNE竞争可视化质量，并且可以保留更多的全局结构，具有出色的运行时性能。 此外，UMAP对嵌入维度没有计算限制，使其可用作机器学习的通用降维技术

UMAP(https://github.com/lmcinnes/umap)

2)代码实现

reticulate::py_install(packages = "umap-learn") #安装UMAP
# UMAP
pbmc <- RunUMAP(object = pbmc, dims = 1:10)
DimPlot(object = pbmc, reduction = 'umap',label = TRUE )
# TSNE
pbmc <- RunTSNE(object = pbmc, dims = 1:10)
DimPlot(object = pbmc, reduction = 'tsne',label = TRUE )
# 您可以在此时保存对象，以便可以轻松地将其重新加载，而无需重新运行上面执行的计算密集型步骤，或者可以轻松地与协作者共享
saveRDS(pbmc, file = "G:/Ajitongjiuniversity/Bioinformatics/Single cell sequencing/scRNA_seq workflow/seurat/pbmc_tutorial.rds")

3)结果

UMAP图 TSNE图

4）收获

哈哈无意中指导了UMAP方法开心

5. Finding differentially expressed features 寻找差异表达的特征（聚类生物标志物）

1）理解

• Seurat可以帮助您找到通过差异表达定义聚类的标记。
• 默认情况下，它标识单个cluster的正marker和负marker（在ident.1中指定），与所有其他cluster进行比较。
• FindAllMarkers为所有cluster自动执行此过程，但您也可以测试集群彼此之间或所有cluster。
• FindAllMarkers参数：min.pct参数在两组cluster中的任意一组中以最小百分比检测到特征；ident.1：用于定义标记的标识类；ident.2：用于和ident.1比较的标识类；logfc.threshold: Log FC的倍数值
• FindAllMarkers（参数是test.use）中用于找两个cluster之间的差异基因的方法有：wilcox（默认）；bimod；roc；t；negbinom；poisson；LR；MAST；DESeq2
• %>% R中的管道传输函数

2)代码实现

#- (1)找出所有cluster中的marker
# find all markers of cluster 1 找出cluster1 的所有marker
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
head(x = cluster1.markers, n = 5)
# find all markers distinguishing cluster 5 from clusters 0 and 3 找到cluster 5 和cluster(0,3)比较的所有marker
cluster5.markers <- FindMarkers(object = pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(x = cluster5.markers, n = 5)
# find markers for every cluster compared to all remaining cells, report only the positive ones 找所有cluster的marker
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC) #将pbmc.markers传给group_by，group_by按cluster 排序，再将结果传给top_n，top_n按avg_logFC排序，显示每个类中的前两个

#- (2)可视化这些marker在这些集群中的表达情况

#--1）用VlnPlot(显示跨集群的表达概率分布)
#--2）FeaturePlot（在tSNE或PCA图上可视化特征表达）这两种是比较常用的
#--3）还可以用`RidgePlot`，`CellScatter`和`DotPlot`这几种方法查看数据集情况

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

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

#RidgePlot
RidgePlot(object = pbmc, features = c("MS4A1", "CD79A"))

#CellScatter
CellScatter(object = pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"),cell1 = "CATTACACGGAGTG", cell2 = "CATTACACCAACTG")

# DotPlot
DotPlot(object = pbmc, features = c("MS4A1", "CD79A"))

# DoHeatmap 热图

# 观察每个cluster中的TOP10的基因在每个cluster中的表达情况
pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC) -> top10
DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()

3)结果

VlnPlot 图1
VlnPlot 图2
FeaturePlot
RidgePlot
CellScatter
DotPlot
DoHeatmap

6.Assigning cell type identity to clusters 鉴定每个cluster的身份

(1)使用经典的marker

1）代码实现

#我们可以使用经典的标记轻松地将无偏聚类与已知细胞类型相匹配
new.cluster.ids <- c("Memory CD4 T", "CD14+ Mono", "Naive CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Mk")
names(x = new.cluster.ids) <- levels(x = pbmc)
pbmc <- RenameIdents(object = pbmc, new.cluster.ids)
DimPlot(object = pbmc, reduction = 'umap', label = TRUE, pt.size = 0.5) + NoLegend()

library(ggplot2)
plot <- DimPlot(object = pbmc, reduction = "umap", label = TRUE, label.size = 4.5) + xlab("UMAP 1") + ylab("UMAP 2") + 
  theme(axis.title = element_text(size = 18), legend.text = element_text(size = 18)) + 
  guides(colour = guide_legend(override.aes = list(size = 10)))
ggsave(filename = "G:/Ajitongjiuniversity/Bioinformatics/Single cell sequencing/scRNA_seq workflow/seurat/pbmc3k_umap.png", height = 7, width = 12, plot = plot)
saveRDS(pbmc, file = "G:/Ajitongjiuniversity/Bioinformatics/Single cell sequencing/scRNA_seq workflow/seurat/pbmc_tutorial.rds")

2）结果展示

Dimplot 好漂亮

(2)全自动细胞注释-SingleR(暂未跑通)

7.总结

• 今天下午13.30到晚上20.26的收获吧，算是把seurat搞清楚了，比之前刚用得时候清楚。
• 收获的点在于:再写一遍笔记的时候会比较着急，我自己是用notepaid写了一遍了，然后想着用简书也整理一遍，才发现原来自己很着急哦，很多笔记都是直接copy自己写在notepaid中的，希望自己下次能够把简书和Notepaid同步起来，开心
• 不理解的地方：高度变异基因和差异基因的区别？还有SingleR是如何注释的还是需要搞清楚，这个还是很有用的