单细胞数据分析中的秩和检验与t检验

2023-02-24 18:14:02Python016

单细胞数据分析中的秩和检验与t检验,第1张

在单细胞数据分析的过程中，寻找差异基因的过程需要用到对基因统计的假设检验（例如函数FindAllMarkers中的test.use参数），我们这里来深入了解一下假设检验的方法和应用环境。

秩和检验适用于广泛的统计学环境，秩和检验是检验总体分布位置是否相同，因而称为非参数检验（Nonparametric test）。秩和检验（rank sum test）是一类常用的非参数检验。秩和检验首先将数据按从小到大或等级从弱到强转换成秩（也就是顺序），然后求秩和并计算秩和统计量，最后做出统计推断。本文简单介绍秩和检验的原理并基于R语言进行秩和检验的操作。

假设我们从总体A和总体B中分别采样n_a和n_b个样本构成样本集合a和b。通过将样本集a和b中的所有样本按从小到大顺序转化为秩之后我们可以通过绘图的方式对转换的结果进行展示，在图中我们使用“•”代表来自样本集a，使用“o”代表数据来自样本集b。

如果总体A和总体B总体分布位置分布相同（H_0:A=B），那么转换的结果如下图所示：

首先是python（范例）,借助于python模块scipy来实现。

其次是R的实现：(wilcox.test的函数)

这里可以发现，秩和检验仅仅和数据的总体分布有关，适用于一般的环境 ，在单细胞数据中寻找markergene 的过程中，大部分默认就是采用此方法，当然，这种检验只是一种很常规的检验，离我们真正的数据分析还很遥远。

t检验，亦称student t检验（Student's t test），主要用于样本含量较小（例如n <30），总体标准差σ未知的正态分布。 [1] t检验是用t分布理论来推论差异发生的概率，从而比较两个平均数的差异是否显著。它与 f检验、卡方检验并列。

这里我们需要注意一下：

（1）t检验对于大样本分布需要转换，而我们单细胞的数据分布属于大样本分布。

（2）数据分布为正态分布，单细胞数据是否为正态分布，在我的文章单细胞数据分析之PCA再认识与ScaleData函数做了详细的介绍，大家可以看一下。

t检验最常见的四个用途：

1、单样本均值检验（One-sample t-test）

用于检验总体方差未知、正态数据或近似正态的单样本的均值是否与已知的总体均值相等

2、两独立样本均值检验（Independent two-sample t-test）

用于检验两对独立的正态数据或近似正态的样本的均值是否相等，这里可根据总体方差是否相等分类讨论

3、配对样本均值检验（Dependent t-test for paired samples）

用于检验一对配对样本的均值的差是否等于某一个值

4、回归系数的显著性检验（t-test for regression coefficient significance）

用于检验回归模型的解释变量对被解释变量是否有显著影响。

单样本T检验用于比较一组数据与一个特定数值之间的差异情况。

应用场景：

某个医生检测40名从事铅作业工人的血红蛋白含量，其均数为130.83g/L，标准差为25.74g/L，试分析从事铅作业的工人血红蛋白含量是否不同于正常成年人平均值140g/L？

我们来看一下这个结果，以p=0.2696，以p=0.05为界，说明没有统计意义。

两独立样本t检验的目的是利用来自两个总体的独立样本，推断两个总体的均值是否存在显著差异。

2、使用的前提条件

（1）两个样本应该是相互独立的；

（2）样本来自的两个总体应该服从正态分布。

显然单细胞使用的就是两独立样本均值检验。

用于分析配对定量数据之间的差异对比关系。与独立样本t检验相比，配对样本T检验要求样本是配对的。两个样本的样本量要相同；样本先后的顺序是一一对应的。

配对样本t检验用于样品的两个相关组之间的比较手段。在这种情况下，同一样本有两个值（即一对值）。

举个例子，在1个月内有20只小鼠接受了治疗X。我们想知道处理X是否会对小鼠的体重产生影响。

为了回答这个问题，在治疗之前和之后测量了20只小鼠的体重。通过测量相同小鼠体重的两次，我们得到了治疗前的20组值和治疗后的20组值。

在这种情况下，可以使用配对t检验比较治疗前后的平均体重。

似然比(likelihood ratio, LR) 是反映真实性的一种指标，属于同时反映灵敏度和特异度的复合指标。即有病者中得出某一筛检试验结果的概率与无病者得出这一概率的比值。该指标全面反映筛检试验的诊断价值，且非常稳定。似然比的计算只涉及到灵敏度与特异度，不受患病率的影响。因检验结果有阳性与阴性之分，似然比可相应地区分为阳性似然比(positive likelihood ratio, ＋LR)和阴性似然比(negative likelihood ratio, －LR)。阳性似然比是筛检结果的真阳性率与假阳性率之比。说明筛检试验正确判断阳性的可能性是错误判断阳性可能性的倍数。比值越大，试验结果阳性时为真阳性的概率越大。阴性似然比是筛检结果的假阴性率与真阴性率之比。表示错误判断阴性的可能性是正确判断阴性可能性的倍数。其比值越小，试验结果阴性时为真阴阳性的可能性越大。

似然比检验（likelihood ratio test, LRT） 是一种检验 参数能否反映真实约束 的方法（分布或模型的某参数 θ 等于 θ 0 是否为真实约束）。似然比检验的思想是：“如果参数约束是有效的，那么加上这样的约束不应该引起似然函数最大值的大幅度降低。也就是说似然比检验的实质是在 比较有约束条件下的似然函数最大值与无约束条件下似然函数最大值 。” 可以看出，似然比检验是一种通用的检验方法（比 t 检验、 Κ 2 检验等具有更广的适用范围）。

这个有点难，我们不展开讨论了，主要就是检验分群结果结束以后，基因的表达分布是否是受到约束的

Identifies 'markers' of gene expression using ROC analysis. For each gene, evaluates (using AUC) a classifier built on that gene alone, to classify between two groups of cells. An AUC value of 1 means that expression values for this gene alone can perfectly classify the two groupings (i.e. Each of the cells in cells.1 exhibit a higher level than each of the cells in cells.2). An AUC value of 0 also means there is perfect classification, but in the other direction. A value of 0.5 implies that the gene has no predictive power to classify the two groups. Returns a 'predictive power' (abs(AUC-0.5) * 2) ranked matrix of putative differentially expressed genes.

关于roc的讲解在我的文章里深入理解R包AUcell对于分析单细胞的作用详细提到过，大家可以看一下。

1. Seurat对象查看当前的Assay

在进行了SCTransform操作后，矩阵默认会变成SCT矩阵，如果不加设置，后续的PCA等操作都是基于SCT矩阵。

修改DefaultAssay：

2. Seurat使用FindVariables找到高变基因后增删高变基因

3. 不同运行步骤中的文件来源和储存位置⚠️

⚠️：PCA的值是可以被覆盖的，使用三步法对矩阵进行标准化后进行PCA后再使用SCT矩阵进行标准化，PCA的矩阵变成了SCT的PCA矩阵，原有的PCA矩阵不会保留。后续的TSNE和UMAP降维图也和三步法不一样。

4. @data标准化矩阵和 @scale.data 归一化矩阵的区别

单细胞RNA 测序数据中，文库之间测序覆盖率的系统差异通常是由细胞间的cDNA 捕获或PCR 扩增效率方面的技术差异引起的，这归因于用最少的起始材料难以实现一致的文库制备。标准化旨在消除这些差异（例如长度，GC 含量），以使它们不干扰细胞之间表达谱的比较。这样可以确保在细胞群体中观察到的任何异质性或差异表达都是由生物学而不是技术偏倚引起的（批次矫正仅在批次之间发生，并且必须同时考虑技术偏差和生物学差异，标准化只需考虑技术差异）。

软件Seurat 提供了三种标准化的方法，分别为LogNormalize、CLR、RC，通常情况下我们采用LogNormalize 的方式进行标准化，计算公式为：log1p((Feature counts/total counts) ∗ scale. factor)

归一化的目的则是使特征具有相同的度量尺度

参考： Seurat的normalization和scaling

5. 关于有些细胞属于同一个cluster但是在umap或者tsne图上相聚较远的问题：UMAP和TSNE是各自的算法在PCA降维的基础上再进行非线形降维，在二维图上把其各自算法认为相近的细胞聚在一起。但是FindClusters输入的不是UMAP或TSNE降维的数据，而是FindNeighbors的数据，而FindNeighbors输入的数据是PCA降维数据，是用另外一种算法计算的细胞之间的距离。因此会出现有些细胞被认为是同一个cluster，但是在umap或者tsne图上相聚较远（尤其是一些散在的，脱离主群的细胞）

6. marker基因鉴定，查看marker基因的表达是使用RNA矩阵还是sct矩阵？

这是一个争议性问题，两个都可以，目前建议最好使用RNA矩阵。

sct的到的count并不是真实的基因表达值，而是通过scaledata倒推出来的，它是一个回归，运算之后的残差。

7. 关于FindAllMarks找到的基因

如下图，先看cluster0的Marker基因：cluster0的差异基因是cluster0的细胞和剩下的所有的cluster合在一起的细胞做对比得到的。pct.1是这个基因在cluster0中的表达比例（S100A8在cluster0的细胞中的表达比例是100%），pct.2是这个基因在除了cluster0以外的所有细胞中的表达比例（S100A8在除cluster0以外的细胞中的表达比例是51.2%）。avg_log2FC是表达差异倍数，p_val_adj是校正后的p值。

8. 在提取Marker基因时比较好的办法：因为单细胞矩阵算出来的结果，p_val_adj==0的有很多，所以可以先把p_val_adj==0先提出来。再把p_val_adj<0.01的按差异倍数靠前的20/30/500...（按需要）个基因提出来，然后把两个矩阵合在一起（取交集）用来做细胞鉴定。（结合p值和fc来做筛选效果更好）

⚠️：提取没有核糖体和线粒体的marker基因更好。（这些基因对鉴定没有帮助）

有些基因比如Foxp3，对细胞鉴定很重要，但常常在筛选Marker基因的时候筛选不出来。非负矩阵分解可能更好。

参考：过滤线粒体核糖体基因

9. 提取亚群

⚠️ 新提取的亚群需要重新进行降维聚类 （和大群相比，标记基因发生了变化），并重新寻找marker基因，重新分群，注释。❗️subset提取子集后，不同样本间批次校正的信息也被去除，需要重新进行批次矫正

参考： Seurat取子集时会用到的函数和方法 ⚠️⚠️⚠️

10. 取子集后如何去除空子集（还存在这个level，但里面包含的细胞为0，如何去除） as.factor(as.character())

11. 双细胞的预测和去除如DoubletFinder建议单样本进行，不建议双样本一起预测。除此之外，其他步骤都可以多样本一起做，质控的时候也可以多样本一起做，但是建议每个样本都单独看一看。

12. 单细胞多样本整合：merge()；多样本拆分：SplitObject()

13. 在做多组数据整合，每个组又有多个样本的时候，最好把单独的每个样本当成批次，而不是把不同的组当成批次。

14. 多核运算

参考：单细胞数据分析中future包的使用

15. pbmc3k.final@commands$FindClusters 可以查看FindClusters运行时间和记录。Seurat是记录其分析过程的，也可查看command下其他操作

16. 关于质控标准：同一组织的最好用同样的标准，不同组织的可以不一样。不同组织线粒体含量等可能存在差异。

17. 可视化的方法总结

参考： https://www.jianshu.com/p/0d1e2c7d21a4

18. circos图绘制

19. 单细胞数据思维导图，有利于查看单细胞数据格式。

https://www.jianshu.com/p/7560f4fd0d77

20. 对于旧版本Seurat对象的更新

scRNA <- UpdateSeuratObject(scRNA)

UpdateSeuratObject {SeuratObject} ：Update old Seurat object to accommodate new features

21. 对有些操作需要用到python设置的情况

22. 单细胞数据做pooling的好处：可以尽量的降低dropout的问题。（dropout就是矩阵中的zero，这些zero实际上并不是0，而是每个液滴里面起始反应量太低了。而一般的反转录效率只能到30%左右，70%的转录本实际上在反转录那一步是被丢掉的，这是单细胞测序一个比较大的问题）。

但是一旦做了pooling，你必须要证明pooling对结果是没有影响的（下图的右面三个图）。

23. Seurat的VlnPlot中的combine参数，在如下画三个基因的情况下，设置成T就画一张图，设置成False，会将三个基因各画一张图。

24. rev()这一步是将横坐标的基因反过来排序

这两个画出来的图横坐标基因的顺序是相反的（见NicheNet）

25. 堆叠小提琴图的绘制

完成这个需求有以下几种实现方法：

1. Seurat包直接就可以实现（stack = T）

2. 通过Seuart->scanpy来实现，第一张是Seurat包VlnPlot函数画的图，第二张是scanpy中stacked_violin函数画的图，那么现在问题就变成为Seurat对象到scanpy对象的转换

3. 用R原生函数实现StackedVlnPlot的方法

4. 使用基于scanpy包衍生的scanyuan包来说实现

5. 使用R包MySeuratWrappers来实现

最简单的方法1如下：

如果不设置level，会按字母顺序排列，case会自动排在con前面。

使用Seurat的 RenameIdents 函数也可以

x: table

margin: a vector giving the margins to split by. E.g., for a matrix 1 indicates rows, 2 indicates columns, c(1, 2) indicates rows and columns. When x has named dimnames, it can be a character vector selecting dimension names.

得到的HC_1样本的orig.ident默认是样本名中第一个_号的前一部分。所以要保证矩阵的列名是样本名_细胞barcode 这样的格式。

如果有多个分组，例如两个样本矩阵中细胞分别命名为HC_1_barcode，HC_2_barcode，在直接通过如下方法得到两个Seurat对象，再对其进行merge之后，两个样本会被合并成一个。也就是样本信息只保留了第一个_号之前的HC，没有保留_号之后的1和2。

为了避免这种情况，可以在构建Seurat对象时通过参数进行设置

⚠️PC数的选择：Seurat官网提供的三种方法只能给出PC数的粗略范围，选择不同PC数目，细胞聚类效果差别较大，因此，需要一个更具体的PC数目。作者提出一个确定PC阈值的三个标准：

一般先选默认分辨率(0.8)，大概可能会分出十几个群。因为最终都是要注释到每一个barcode，所以首先可以看大类marker的分布(不受分辨率影响)，可以根据marker基因的分布来调整分辨率。是否需要精细的分群得看精细的分群对研究有没有决定作用，还有很重要的一点是 看分出的各个cluster在Findallmarkers给出的结果中marker的热图是不是能明显分开 。精细划分的细胞本来就很类似，如果有部分小群的热图明显分不开或者非常类似，就可以考虑把分辨率调小。

这实际上是没有必要的必须保持一致的。下游的都是用pca之后的，pca是为了压缩数据。

umap和tsne是为了可视化(仅仅是可视化)，但是FindNeighbor是计算细胞间距离矩阵。找类群数目和可视化可以说没有关系。

map函数：

R语言循环第三境界：purrr包map函数！

浅析R语言中map(映射)与reduce(规约)

参考： monocle2

查看不同细胞群的中位基因也是一样

查看不同样品的中位基因也是一样

或者也可以

时代的洪流奔涌而至，单细胞技术也从旧时王谢堂前燕，飞入寻常百姓家。雪崩的时候，没有一片雪花是无辜的，你我也从素不相识，到被一起卷入单细胞天地。R语言和Seurat已以势如破竹之势进入4.0时代，天问一号探测器已进入火星乌托邦平原了，你还不会单细胞吗？那么为了不被时代抛弃的太远，跟着我们一起通过学习seurat系列教程入门单细胞吧。

首先我们学习经典的标准流程，这个教程我当初入门时候曾经花1000购买过，也曾花6000购买过，不同的单细胞培训班，买的是一样的教程。现在免费送给你，别说话，开始学吧！

Seurat4.0系列教程1：标准流程 - (jianshu.com)

Seurat4.0系列教程2：多模式数据联合分析 - (jianshu.com)

Seurat4.0系列教程3：合并两个样品的10x数据集 - (jianshu.com)

Seurat4.0系列教程4：整合分析 - (jianshu.com)

Seurat4.0系列教程5：交互技巧 - (jianshu.com)

Seurat4.0系列教程6：常用命令 - (jianshu.com)

Seurat4.0系列教程7：数据可视化方法 - (jianshu.com)

Seurat4.0系列教程8：细胞周期评分和回归分析 - (jianshu.com)

Seurat4.0系列教程9：差异表达检测 - (jianshu.com)

Seurat4.0系列教程10：降维 - (jianshu.com)

Seurat4.0系列教程11：使用sctransform - (jianshu.com)

Seurat4.0系列教程12：大数据集整合的方法 - (jianshu.com)

Seurat4.0系列教程13：使用互惠 PCA （RPCA）快速整合数据 - (jianshu.com)

Seurat4.0系列教程14：整合scRNA-seq and scATAC-seq数据 - (jianshu.com)

Seurat4.0系列教程15：映射和注释查询数据集 - (jianshu.com)

Seurat4.0系列教程16：多模式参考映射注释细胞 - (jianshu.com)

Seurat4.0系列教程17：Mixscape - (jianshu.com)

Seurat4.0系列教程18：Weighted Nearest Neighbor Analysis（WNN)） - (jianshu.com)

Seurat4.0系列教程19：多线程并行策略 - (jianshu.com)

Seurat4.0系列教程20：单细胞对象的格式转换 - (jianshu.com)

Seurat4.0系列教程21：结合Cell Hashing分析双细胞 - (jianshu.com)

Seurat4.0系列教程22：空间转录组的分析、可视化与整合 - (jianshu.com)

Seurat4.0系列教程告一段落，但这决不是终点。这个系列教程只是给大家打开一扇窗，知道Seurat4.0有这些功能可用，少走弯路，起到一个抛转引玉的作用，后续还要自己深入研究。不要像我当初入门单细胞之时，为了找到整合方法耗费大量时间及不必要的金钱。君子生非异也，善假于物也。学，然后知不足。加油吧！少年！用好单细胞测序及分析这个技术，为人类健康研究做出自己的贡献！不足之处，恳请批评指正！------以此纪念粉丝突破1000。

CellChat三部曲1：使用CellChat对单个数据集进行细胞间通讯分析 - (jianshu.com)

CellChat 三部曲2：使用CellChat 对多个数据集细胞通讯进行比较分析 - (jianshu.com)

CellChat 三部曲3：具有不同细胞类型成分的多个数据集的细胞通讯比较分析 - (jianshu.com)

RNAvelocity系列教程1：RNA速率简介及scVelo安装 - (jianshu.com)

RNAvelocity系列教程2：使用 Seurat 和 scVelo 估计 RNA 速率 - (jianshu.com)

RNAvelocity系列教程3：使用Seurat和velocyto估算RNA速率 - (jianshu.com)

RNAvelocity系列教程4：velocyto.R的使用教程 - (jianshu.com)

RNAvelocity系列教程5：scVelo 的基本知识 - (jianshu.com)

RNAvelocity系列教程6：# scVelo用于RNA 速率基本流程 - (jianshu.com)

RNAvelocity系列教程7：# scVelo实战-胰腺数据集 - (jianshu.com)

RNAvelocity系列教程8：# scVelo实战-齿状回数据集 - (jianshu.com)

RNAvelocity系列教程9：# scVelo应用-动力学模型 - (jianshu.com)