单细胞测序分析软件包

20210829修改

本笔记主要参考seurat的网址https://satijalab.org/seurat/
主要是为了理清思路，记录一下，方便以后查找，希望也能够对各位同学有所帮助

Seurat软件介绍

seurat包是由Satija Lab开发，目前已经是v4版本

seurat安装

安装方便

install.packages('Seurat')

安装网址
https://satijalab.org/seurat/articles/install.html

导入分析需要的包

library(Seurat)
library(dplyr)
library(magrittr)

读入矩阵

读入的是10X网站中自带数据集，下载地址https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
该数据集包括以下：matrix.mtx, genes.tsv, and barcodes.tsv
#使用Read10X函数读取矩阵数据

pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_feature_bc_matrix/")

#使用CreateSeuratObject函数创建seurat对象，这里条件过滤条件是细胞中基因数目不能小于200，且基因至少在三个细胞中有表达

pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
#查看30个细胞中的基因表达情况
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]

其中，在矩阵“."代表没有检测到分子，因为在单细胞矩阵中大部分值是0，所以单细胞中用的是稀疏矩阵(sparse-matrix )。

一、数据前处理

数据前处理包括：基于QC度量的细胞选择与筛选，数据标化与缩放，及高度可变特征的检测。

1.质控与细胞选择

seurat列出了以下质控标准：

每个细胞中检测到的独特基因的数目
1.低质量的细胞或空的油滴通常会有很少的基因量
2.一个油滴中包含多个细胞基因表达将会特别高
细胞内检测到的分子总数（与独特的基因密切相关）
线粒体基因组的百分比
1.低质量或是死细胞会表现出线粒体污染；
2.使用PercentageFeatureSet()函数计算线粒体QC度量；
3.使用以“MT-”开头的作为线粒体基因

Seurat中引用了一篇质控文献
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4758103/

"[["这个操作符可以在列中添加到对象中，使用从MT-作为线粒体基因集

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

在创建对象 CreateSeuratObject() 的过程中，会自动计算细胞中独特基因与总基因数目，可以在目标对象中找到

head(pbmc@meta.data, 5)

QC指标可视化

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

从可视化结果进行细胞选择，过滤掉超过2500或少于200基因的细胞，同时过滤线粒体比例超过5%的细胞。

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

#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.

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

数据标准化

过滤细胞后，进行标准化数据，默认情况下，使用“lognormalize”全局缩放的归一化方法，通过总表达值对每个细胞的基因表达值归一化，并将其乘以缩放因子(默认为10,000)，最后对结果进行对数变换，标准化数据存储在pbmc@assays$RNA@data中。

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

不使用默认值，也可以使用以下代码

pbmc <- NormalizeData(pbmc)

特征性基因的选择

计算数据集中表现出高细胞间差异的基因子集(即，它们在一些细胞中高表达，而在另一些细胞中低表达)。在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。默认情况下，每个数据集选择2000个高变异基因用于下游分析。

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

#前10的高变基因

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

#将前10的高变基因在图中标注出来

plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

scale

接下来使用线性变换(缩放)处理数据，也是PCA的前处理阶段。
过程：

改变每个基因的表达，使细胞间的平均表达为0；
缩放每个基因的表达，使细胞间的差异为1（这一步在下游分析中给予相同的权重，使得高表达基因就不会影响后续分析）；
缩放的数据储存在pbmc[[“RNA”]]@scale.data / pbmc@assays$RNA@scale.data中。

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

线性降维

对缩放后的数据进行降维，默认降至50维。

PCA(Principal Component Analysis)，主成分分析，是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上，在尽可能保留原始数据信息的同时降低数据维度来加速数据分析。

过程就是从原始高维的空间中按顺序地找一组相互正交的坐标轴系统，新的坐标轴的选择与数据本身是密切相关的。

其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴。大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴，实现对数据降维。

默认情况下，只使用前面确定的变量特性作为输入，但如果使用一个不同的子集，可以使用feature参数定义变量特性

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

使用不同的方法检查和可视化PCA结果

#选择前5个维度进行查看
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
#使用DimPlot函数
DimPlot(pbmc, reduction = "pca")VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

使用DimHeatmap可视化函数查看样品中前500个细胞在前6个PCA中的热图：

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

确定数据集的维度

每个pc本质上代表一个“元特征”，它将相关特征集中的信息组合在一起。因此，越在顶部的主成分越可能代表数据集。然而，我们应该选择多少个主成分才认为我们选择的数据包含了绝大部分的原始数据信息呢？我们可以使用JackStraw函数来选择PC数目

在JackStraw()函数中, 使用基于零分布的置换检验方法。随机抽取一部分基因(默认1%)然后进行pca分析得到pca分数，将这部分基因的pca分数与先前计算的pca分数进行比较得到显著性p-Value。根据主成分（pc）所包含基因的p-value进行判断选择主成分。最终的结果是每个基因与每个主成分的关联的p-Value。保留下来的主成分是那些富集小的p-Value基因的主成分。JackStrawPlot()函数提供可视化方法，用于比较每一个主成分的p-value的分布，虚线是均匀分布；显著的主成分富集有小p-Value基因，实线位于虚线左上方。下图表明保留10个pca主成分用于后续分析是比较合理的。

pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)

另一种方法是“ElbowPlot”
基于每个主成分所解释的方差百分比的排序(ElbowPlot()函数)。在这个例子中，我们可以观察到PC9-10周围的一个“elbow"，这表明大部分真实信号是在前10台pc中捕获的。图中发现第9个主成分是一个拐点，后续的主成分(PC)变化都不大了。

ElbowPlot(pbmc)

注意：鉴别数据的真实维度不是件容易的事情；除了上面两种方法，Serat官方文档还建议将主成分(数据异质性的相关来源有关)与GSEA分析相结合。Dendritic cell 和 NK aficionados可能识别的基因与主成分 12 和 13相关，定义了罕见的免疫亚群 (i.e. MZB1 is a marker for plasmacytoid DCs)。如果不是事先知道的情况下，很难发现这些问题。

Serat文档因此鼓励用户使用不同数量的PC（10、15，甚至50）重复下游分析。然而也将观察到的，结果通常没有显著差异。因此，在选择此参数时，可以尽量选大一点的维度，维度太小的话对结果会产生不好的影响。

细胞聚类

Seurat v3应用基于图形的聚类方法，例如KNN方法。具有相似基因表达模式的细胞之间绘制边缘，然后将他们划分为一个内联群体。

在PhenoGraph中，首先基于pca维度中（先前计算的pca数据）计算欧式距离（the euclidean distance），然后根据两个细胞在局部的重合情况（Jaccard 相似系数）优化两个细胞之间的边缘权值。此步骤内置于FindNeighbors函数，输入时先前确定的pc数据。

为了聚类细胞，接下来应用模块化优化技术迭代将细胞聚集在一起。（the Louvain algorithm (default) or SLM [SLM, Blondel et al., Journal of Statistical Mechanics]），FindClusters函数实现这一功能，其中需要注意resolution参数，该参数设置下游聚类分析的“granularity”，更大的resolution会导致跟多的细胞类群。3000左右的细胞量，设置resolution为0.4-1.2是比较合适的。细胞数据集越大，需要更大的resolution参数, 会获得更多的细胞聚类。
查看细胞属于那个类群可以使用函数Idents。

#计算最邻近距离
pbmc <- FindNeighbors(pbmc, dims = 1:10)
#聚类，包含设置下游聚类的“间隔尺度”的分辨率参数resolution ，增加值会导致更多的聚类。
pbmc <- FindClusters(pbmc, resolution = 0.5)
#可以使用idents函数找到聚类情况：
#查看前5个细胞的聚类id
head(Idents(pbmc), 5)

非线性降维分析

Seurat提供了一些非线性降维度分析的方法，如tSNE和UMAP，以可视化和探索这些数据集；将上一步PC维度中的数据进一步降至二维，实现数据的可视化。之前pca是线性降维，能很好处理高维度数据，tsne和umap是非线性降维方式，对低维度数据处理速度快，但难以处理高维度数据。
相比tsne，umap能更好地反映原始数据结构，有着更好的连续性。

# install UMAP： reticulate::py_install(packages ='umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)
pbmc <- RunTSNE(pbmc, dims = 1:10)
#画图
DimPlot(pbmc, reduction = "umap")
DimPlot(pbmc, reduction = "tsne")#添加细胞类群xiba的标签
DimPlot(pbmc, reduction = "umap",label = TRUE)
LabelClusters(DimPlot(pbmc, reduction = "umap"),id = 'ident')#可以在这里保存为以下形式，就不必再执行以上步骤
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")

差异分析

seurat可以通过差异表达找到每个聚类的marker gene，差异分析可以有多种形式，如找到所有聚类的marker gene（如cluster1中所有的markgene是指cluster1相对于其余所有cluster是差异的）、两个cluster之间的差异分析、某个cluster中两个样品之间差异分析等。

seurat可以通过差异表达分析寻找不同细胞类群的标记基因。FindMarkers函数可以进行此操作，但是默认寻找单个类群(参数ident.1)与其他所有类群阳性和阴性标记基因。FindAllMarkers函数会自动寻找每个类群和其他每个类群之间的标记基因。

min.pct参数:设定在两个细胞群中任何一个被检测到的百分比，通过此设定不检测很少表达基因来缩短程序运行时间。默认0.1

thresh.test参数：设定在两个细胞群中基因差异表达量。可以设置为0 ，程序运行时间会更长

max.cells.per.ident参数：每个类群细胞抽样设置；也可以缩短程序运行时间。

# find all markers of cluster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)# find all markers distinguishing cluster 5 from clusters 0 and 3(找到cluster5和cluster0、3之间的markgene)
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)# find markers for every cluster compared to all remaining cells, report only the positive ones(找到每个cluster相比于其余cluster的markgene，只报道阳性的markgene)
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)

Seurat可以通过参数test.use设定检验差异表达的方法(详情见 DE vignett/ https://satijalab.org/seurat/archive/v3.0/de_vignette.html)。

cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
head(cluster1.markers, n = 5)

添加细胞注释信息

如果你对每个cluster已经鉴定出细胞类型，就可以在可视化中将细胞类型标注在图中。

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

可视化Seurat有多种方法可视化标记基因的方法

VlnPlot: 基于细胞类群的基因表达概率分布
FeaturePlot：在tSNE 或 PCA图中画出基因表达情况
RidgePlot，CellScatter，DotPlot

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

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

#查看特定基因在聚类图中的表达量分布情况

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

#查看特定基因在cluster中的表达山脊图

features <- c("LYZ", "CCL5")
RidgePlot(pbmc, features = features, ncol = 2)

#查看特定基因在cluster中的分布点图，点的大小代表表达该基因的细胞比例，颜色代表平均表达水平

DotPlot(pbmc, features = features) + RotatedAxis()

DoHeatmap为指定的细胞和基因画表达热图。每个类群默认展示top 20标记基因

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

Assigning cell type identity to clusters

根据已知细胞类型与基因标记的对应关系，可以为细胞类群找到对应的细胞类型。

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

参考文章：
https://www.jianshu.com/p/75c6d55a63bb?from=timeline
https://www.jianshu.com/p/44c1f526e450
https://www.bilibili.com/read/cv6411199/
https://www.jianshu.com/p/fef17a1babc2

默认情况下，我们的seurat对象中是一个叫RNA的Assay。

在我们处理数据的过程中，做整合（integration），或者做变换（SCTransform），或者做去除污染（SoupX），或者是融合velocity的数据等，我们可能会生成新的相关的Assay，用于存放这些处理之后的矩阵。

在之后的处理中，我们可以根据情况使用指定Assay下的数据。不指定Assay使用数据的时候， Seurat给我们调用的是Default Assay下的内容。

可以通过对象名@active.assay查看当前Default Assay，通过DefaultAssay函数更改当前Default Assay。

Assay数据中，counts为raw，data为normalized，scale为scaled

调用Assay中的数据的方式为，以调取一个名为PBMC的Seurat对象中Assay integrate中的nomalized数据为例：

PBMC@assays$RNA@data