哈佛课程笔记-理解DEseq2差异分析原理

前言：网上差异分析的教程一搜一大把，大家都知道了如何走流程，但是对原理模模糊糊。刚学完哈佛大学chan bioinformatics core的差异分析课程，跟大家分享一下学到的DEseq2差异分析原理：https://hbctraining.github.io/DGE_workshop_salmon_online/schedule/links-to-lessons.html

一． RNA-seq数据分析的基础知识

1. RNA-seq测序数据count数的分布模型

我们平常拿到的差异分析数据都长这样：

横轴是基因，纵轴是基因在每个样本中测到的read数，也就是count数。

想知道RNA-seq的每个基因count值是如何分布的？我们来画个柱状图：

纵轴是基因，横轴是count数，这个图揭示了某个样本中所有基因的count数分布。

上图体现了RNA-seq 基因count数分布的3个基本特点

1. count数量少，但分布范围广 –-因为没有最大值限制
2. 非正态分布 （第三个埋藏在后文中，看看谁能找到）

对于这种类型的数据，很明显正态分布模型就不适用了（注意，这个与芯片数据不同哦），那我们应该用什么模型去拟合RNA-seq gene的count数的分布呢？

泊松分布，二项分布都在一定程度上能拟合这种偏态分布，但是目前的观点认为，负二项分布在拟合RNA-seq的count数分布上表现最佳，因为RNA-seq counts数还具有这样的分布规律：

随着基因表达量的增加，基因count数分布的方差快速增加。低表达基因中count数分布的方差大小不一（学名：具有方差异质性）。泊松分布和二项分布表示，这么奇怪的分布，我干不来，找别人吧！哼！但是负二项分布对于这种分布拟合度较高，它的公式长这样：（给有需要的人吧）

1. 生物学重复: 我们做RNA-seq的时候抓耳挠腮：多少个样本才够哇？技术重复能不能代替生物学重复？生物学重复对结果的影响有多大？

别急，哈佛大佬速速来解答疑问：

首先，为啥要那么多重复呢?因为重复会让差异分析更加精确，可信！基因表达差异也有能因为一些无关的原因，实验污染，不明确的技术偏差等等。（如上图），我们基因差异分析的目的找到的差异是我们实验组和对照组之间真正有意义的生物学差异。因此，足够的生物学重复（3组以上）是必要滴

举个栗子：geneA平均表达在处理后是处理前的两倍，但是这么大的组内离散真的会有显著的组间差异吗？当然不行。而且这种离散度的趋势只有3个以上的重复才能充分体现，所以不能偷懒哦。

还不信？上图！

https://doi.org/10.1093/bioinformatics/btt688 这篇2014发表在bioinformatics的文章表明，生物学重复比增加测序深度更能提高差异分析的精确度和可信度。

二. DEseq2 测序流程详解 讲了半天铺垫，想必大家已经迫不及待了，来，上正菜！：

DEseq2测序的流程

1. Quality control：质控分两种：标准化基因的表达和排除低质量样本

（1）DEseq2的基因标准化原理：DEseq有自己的一套count标准化程序：其实TPM之类的标准化方法虽然解决了基因长度和测序深度的影响的问题，但还是不能解决一个问题：那就是测序文库组成不同造成的差异 这种差异的来源是一个基因被敲减了，完全没表达了，因而影响到了其他基因。

DEseq2使用提高中位表达基因的办法解决了这个问题。想知道具体如何解决的？生信菜鸟团的这篇文章可以完美满足你的需求（向大佬膜拜）：https://mp.weixin.qq.com/s/fw0muJwF-cz2ki2LkflpVw （2） 初始步骤：评估样本之间整体的相似性 请回答下列问题：哪些样本与别的样本相似，哪些不同？样本间的区别来自哪里？这种差异与我们需要找的是否一致？

使用聚类和PCA来得到答案 至于PCA和聚类的原理？你同样可以在生信菜鸟团的专栏里得到答案：https://mp.weixin.qq.com/s/h0pZ0_ZK9BK-K1XI6ZKDzg --PCA的https://mp.weixin.qq.com/s/1rOVTuAxwjgTDThTThtViw --聚类的

1. 差异分析流程：好啦好啦，我知道你已经看不下去了，说了这么多居然还没说到差异分析的重点。行，我保证在500字内完篇。

差异分析分为多个部分：

1.计算离散度

2.拟合并压缩基因的分布使之更适合建模

3.建模并进行统计学检验

1 计算size factors 使用size factors对reads进行标准化（就是我上面说的那个原理，刚刚只是说了原理，这个是在软件中的运行步骤）

1. 计算基因层面的离散度 怎么理解基因表达的离散度？为什么要计算基因的离散度？我们知道：我们需要通过计算基因平均表达的差异找出差异基因，同时考虑组内方差（这个很重要！记不记得生物学重复那里那个方差很大的假差异基因） 我们之前也说过count数方差随着表达量的增加而增大。

DEseq2需要处理这种表达量和方差的关系。我们不希望差异基因都是方差大的高表达基因。所以，DEseq2使用dispersion代替variation, 来表达基因的表达和方差，公式为 Dispersion is calculated by Var = μ + α*μ^2, where α = dispersion, Var = variance, and μ = mean, giving the relationship 因此dispersion与variation还有mean expression关系为：

离散度（dispersion）与表达量成负相关，与方差成正相关 Dispersion可以反应同一平均表达量下基因的离散程度。

如何计算呢？只有3-6个重复难以计算出可信的基因dispersion，所以DEseq2假设有着相同表达的基因有大致相同的统计分布。所以share the information across gene ，用多个大致相同表达的基因的方差，一起算出几个基因的dispersion

1. 拟合基因dispersion曲线 上文我们说过，我们用dispersion合并了方差和表达量，后面我们的统计学检验都是基于dispersion 虽然所有基因都有不同的dispersion，但是所有基因会形成一个规律的分布，而这个分布就能使用负二项分布模型进行线性拟合

红色就是拟合曲线，拟合基因表达与dispersion的关系 黑色是基因点

1. shrink dispersion 压缩离散度 拟合的效果很大程度影响了差异基因的鉴定。为了完成更棒的拟合，我们还需要把dispersion 压缩一下：

压缩的程度由：

（1） gene离线的距离

（2） 样本量 决定 Shrink方法对降低差异分析中的假阳性率至关重要。

因为只有更靠近曲线才能更好的拟合模型。但是过高的gene不会被压缩，因为他们可能是因为技术或生物学原因产生的离群点。（图上被圈起来的小点） 这是一个表示你的数据有没有很好的拟合模型的方法。如果点在曲线周围，说明拟合度较高，如果分散，说明拟合度较低。

plotDispEsts()

这是差的

如果像这样下降又增加。高表达对应的应该是高dispersion这说明高表达基因对应的variation降低了，提示可能有离群样本或者污染。

1. 统计学检验 DEseq2使用Wald test进行统计学检验。Wald也可以用在连续分组上

三．解读结果

1.p-value Adj-p是经过统计学检验后经过FDR调整过的结果 啥？什么是FDR，那你应该看看这篇文章：https://mp.weixin.qq.com/s/dDi7I8LcWSl40JmbkSqCfQ

为啥P-value会出现NA呢？

可能因为在差异分析之前被筛掉了，这样搞更能提高差异分析的效能，DEseq2不会物理移走gene，但是会出现NA，可能出现NA的情况有：

（1） gene在所有样本中都是0

（2） gene中有一个样本出现离群—cook检验

（3） gene有着极低的表达—对应着high dispersion，只要你提高p值阈值，应该不会出现这种情况。也有可能被independent filting给干掉。但只有adj p会受到影响

2.log2Foldchange 超过P值的样本很多，有没有更严格的指标呢？

Log2FC：（需要注意，这里差异的倍数取Log2了哦） log2 (normalized_counts_group1 / normalized_counts_group2) 但是log2FC也有问题，它只考虑了表达值，没考虑组内差异，这该怎么办呢？

DEseq软件已经帮你想好了：

为了更精确地计算log2FC，DEseq2对

（1）low counts

（2）high dispersion 的基因进行了零压缩，在全基因表达的范围内。

例子：绿色和紫色是两个在C57小鼠中表达的gene，紫色gene方差更大。在压缩之前是实线数据，在压缩后是虚线数据。明显，经过压缩，方差的数据LFC减少。因此，虽然两个数据标准化后的counts相同，但LFC不同。

注意！压缩不会改变显著差异的基因数，只会调整logFC的大小。

结束啦，肝完了！皆大欢喜，yes！