Seurat单细胞分析常见代码-02

Posted 2023-04-03

参考技术A 默认是500 * 1024 ^ 2 = 500 Mb

future在seurat的具体应用详见： https://satijalab.org/seurat/archive/v3.0/future_vignette.html

CaseMatch()为seurat包内部函数， https://github.com/satijalab/seurat/blob/master/R/utilities.R

如果细胞太密，有重叠，可以设置透明度，我一般alpha.use设置0.8或者0.9

一般我常用的配色是：c("lightgrey", "red")或者c("lightgrey", "blue")

为了对比细胞间的基因表达差异，常使用不同的配色；
使用viridis调色板

可参考的配色有

我们常常用一些小的特征基因（3-5 个基因）使用 FeaturePlot() 来辅助celltype的鉴定，文献中也有此操作( https://www.nature.com/articles/s42003-020-0922-4/figures/1 )

AddModuleScore()的算法跟上面的计算基因集均值的方法不同，后面细究下源码。
AddModuleScore()为Seurat包中自带函数，需要先计算基因集中所有基因的平均值，再根据平均值把表达矩阵切割成若干份，然后从切割后的每一份中随机抽取对照基因（基因集外的基因）作为背景值。因此，在整合不同样本的情况下，即使使用相同基因集为相同细胞打分，也会产生不同的富集评分；

DotPlot的features为list时，list组间会有间隔，便于比较不同celltype的marker基因差异，可视化效果更好；

计算表达基因的细胞比例的函数，在 https://github.com/satijalab/seurat/issues/371 中有网友给出了一个自定义函数，如下：

同样，也可以计算在cluster的占比

Seurat软件使用操作教程

参考技术A Peripheral Blood Mononuclear Cells (PBMC) 是 10X Genomics dataset page 提供的一个数据，包含2700个单细胞，出自Illumina NextSeq 500平台。

PBMCs是来自健康供体具有相对少量RNA（around 1pg RNA/cell）的原代细胞。在Illumina NextSeq 500平台，检测到2700个单细胞，每个细胞获得69000 reads。

matrix.mtx：matrix.mtx 是 MatrixMarket格式文件；更多内容见： http://math.nist.gov/MatrixMarket/formats.html

总共得到2700个细胞和13714个基因。
count matrix长什么样呢？我们可以看

其中的.表示0，即no molecules detected。当然，这个地方还有另外一种含义就是这个基因是真的没有表达。

由于单细胞测序数据中大多数的值都为0，因此，seurat使用一个稀疏矩阵来保存测序得到的count matrix，这样有利于数据存储空间的节省。

我们来看看使用稀疏矩阵和使用0来存储两种方式的大小对比。

dense为转换为0后的matrix存储大小，709591472 bytes，
sparse为.即稀疏矩阵的大小，29905192 bytes。
两者比值为23.7 倍，即使用系数矩阵来存储单细胞水平的基因表达值非常节省空间。
使用pbmc数据初始化Seurat对象

预处理主要包括基于QC指标的细胞和基因过滤，数据标准化和归一化，高变基因选择。

一般筛选条件有三个：

我们将使用以下标准进行基因过滤：
We filter cells that have unique feature counts over 2,500 or less than 200
We filter cells that have >5% mitochondrial counts
过滤前三个指标可视化

我们就出了小提请图。
查看基因数目, 线粒体基因占比与UMI数目的关系

筛选检测到基因数目超过2500或低于200的细胞
单个细胞中线粒体基因数目占比超过>5%

默认使用数据标准化方法是LogNormalize, 每个细胞总的表达量都标准化到10000，然后log取对数；结果存放于pbmc[["RNA"]]@data。
标准化前，每个细胞总的表达量

鉴定在细胞间表达高度变化的基因，后续研究需要集中于这部分基因。Seurat内置的FindVariableFeatures()函数，首先计算每一个基因的均值和方差，并且直接模拟其关系。默认返回2000个基因。

线性转换缩放数据，ScaleData()函数可以实现此功能。

最终每个基因均值为0，方差为1。

结果存放于pbmc[["RNA"]]@scale.data。

设置参数features是因为ScaleData默认处理前面鉴定的差异基因。这一步怎么做都不会影响到后续pca和聚类，但是会影响做热图。

移除影响方差的因素。

对缩放后的数据进行PCA分析，默认使用前面鉴定表达变化大的基因。使用features参数可以重新定义数据集。

VizDimReduction, DimPlot, 和 DimHeatmap可以从基因或细胞角度可视化pca结果
查看对每个主成分影响比较大的基因集

查看前五的可变基因

可视化对每个主成分影响比较大的基因集

两个主成分的展示

DimHeatmap绘制基于单个主成分的热图，细胞和基因的排序都是基于他们的主成分分数。对于数据异质性的探索是很有帮助的，可以帮助用户选择用于下游分析的主成分维度。

为了避免单个基因影响，Seurat聚类细胞时使用pca结果。首先需要确定的是使用多少个主成分用于后续分析。常用有两种方法，一种是基于零分布的统计检验方法，这种方法耗时且可能不会返回明确结果。另一种是主成分分析常用的启发式评估。
JackStraw()
在JackStraw()函数中, 使用基于零分布的置换检验方法。随机抽取一部分基因(默认1%)然后进行pca分析得到pca分数，将这部分基因的pca分数与先前计算的pca分数进行比较得到显著性p-Value，。根据主成分（pc）所包含基因的p-value进行判断选择主成分。最终的结果是每个基因与每个主成分的关联的p-Value。保留下来的主成分是那些富集小的p-Value基因的主成分。
处理大数据时会花费大量时间；ElbowPlot()内置了一些其它的方法可以减少运行时间。

JackStrawPlot()函数提供可视化方法，用于比较每一个主成分的p-value的分布，虚线是均匀分布；显著的主成分富集有小p-Value基因，实线位于虚线左上方。下图表明保留10个pca主成分用于后续分析是比较合理的。

ElbowPlot

启发式评估方法生成一个Elbow plot图。在图中展示了每个主成分对数据方差的解释情况（百分比表示），并进行排序。根据自己需要选择主成分，图中发现第9个主成分是一个拐点，后续的主成分(PC)变化都不大了。

注意：鉴别数据的真实维度不是件容易的事情；除了上面两种方法，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。

添加细胞标签

此时可以保存数据，方便下次直接导入数据修改图形。

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

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

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

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

Seurat有多种方法可视化标记基因的方法
VlnPlot: 基于细胞类群的基因表达概率分布
FeaturePlot：在tSNE 或 PCA图中画出基因表达情况
RidgePlot，CellScatter，DotPlot

Seurat有多种方法可视化标记基因的方法
VlnPlot: 基于细胞类群的基因表达概率分布
FeaturePlot：在tSNE 或 PCA图中画出基因表达情况
RidgePlot，CellScatter，DotPlot

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

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

12 最后保存这个运行过得文件