表达量矩阵差异分析决定上下调基因的阈值

生信技能树

发布于 2022-06-27 20:41:14

1.2K0

接近十年前，我作为药企小部门的小喽啰，那个时候辅助化学和生物研发团队做非酒精性脂肪肝的药物研发，但是关于这个非酒精性脂肪肝疾病的基本上就十几个相关表达量芯片研究公布在GEO数据库里面。

我们需要批量下载它们并且进行最简单的表达量矩阵差异分析，并且给出来统计学显著的上下调基因。现在看来当然是非常简单了，公众号推文在：

但是那个时候的我第一次接触表达量矩阵芯片，不同芯片产商的不同探针，不同数量值范围，还有差异分析后决定上下调基因的阈值都让我头疼无比。

尤其是logFC的阈值问题，不同文章完全不一样，有使用1，1.2，1.5，甚至2的。最开始跟部门小伙伴讨论希望定下来阈值，但是发现同一个阈值没办法适用于不同数据集，导致汇报给领导的时候有一些数据集就十几个上下调基因，有一些确实好几千个，非常的尴尬。

现在想起来，当然知道为什么了，因为表达量矩阵形式不一样，而且不同数据集里面的两个分组的组间差异和组内差异很不一样。我在生信技能树的教程：《你确定你的差异基因找对了吗？》提到过，必须要对你的转录水平的全局表达矩阵做好质量控制，最好是看到标准3张图：

左边的热图，说明我们实验的两个分组，normal和npc的很多基因表达量是有明显差异的
中间的PCA图，说明我们的normal和npc两个分组非常明显的差异
右边的层次聚类也是如此，说明我们的normal和npc两个分组非常明显的差异

如果分组在3张图里面体现不出来，实际上后续差异分析是有风险的。比如超级多的统计学显著的上下调基因或者说超级少的基因。

神奇的 |logFC| > [mean(|logFC|) + 2sd(|logFC|)]

它依托于一个很常见的统计学理念，就是正态分布，也叫做钟形曲线，有一个概念是置信区间。置信区间是在预先确定好的显著性水平下计算出来的，显著性水平通常称为α(希腊字母alpha)，如前所述，绝大多数情况会将α设为0.05。置信度为(1-α)，或者100×(1-α)%。于是，如果α=0.05，那么置信度则是0.95或95%，如下所示：

钟形曲线

实际的求解步骤也非常简单：

第一步：求一个样本的均值。
第二步：计算出抽样误差。
第三步：用第一步求出的“样本均值”加、减第二步计算的“抽样误差”，得出置信区间的两个端点。

也就是说，0.05这样的阈值其实在钟形曲线上面就换算成了 1.96个标准差（sd值），但是我当时为了方便就直接写成了2个sd，所以就沿用至今。

所以大家如果最近看到我们的火山图，会发现里面有各种稀奇古怪的阈值：

稀奇古怪的阈值

学徒作业

熟读文献：《Identification of transcriptomic signatures and crucial pathways involved in non-alcoholic steatohepatitis》，定位到里面的gse数据集，并且每个数据集独立的差异分析，看看各自的变化倍数的范围，是不是动态变化很大，如果是统一上下调阈值比如1，1.2，1.5，甚至2，多个数据集看看是否有交集。如果不统一阈值，而是各自数据集内部的 |logFC| > [mean(|logFC|) + 2sd(|logFC|)] 作为阈值，也就是说正态分布的95%之外的那些基因，再次看看是否交集会更多。

尤其是跟文献里面的Robust rank aggregation (RRA)算法得到的多个数据集的上下调基因，对比一下。