差异表达基因分析概念篇

Posted 2023-05-02

参考技术A

Differential gene expression analysis：差异表达基因分析

Differentially expressed gene (DEG)：差异表达基因

Volcano Plot：火山图

fold change翻译过来就是倍数变化，假设A基因表达值为1，B表达值为3，那么B的表达就是A的3倍。一般我们都用count、TPM或FPKM来衡量基因表达水平，所以基因表达值肯定是非负数，那么fold change的取值就是(0, +∞).

为什么我们经常看到差异基因里负数代表下调、正数代表上调？因为我们用了log2 fold change。当expr(A) < expr(B)时，B对A的fold change就大于1，log2 fold change就大于0（见下图），B相对A就是上调；当expr(A) > expr(B)时，B对A的fold change就小于1，log2 fold change就小于0。通常为了防止取log2时产生NA，我们会给表达值加1（或者一个极小的数），也就是log2(B+1) - log2(A+1). 【需要一点对数函数的基础知识】

为什么不直接用表达之差，差直接有正负啊？假设A表达为1，B表达为8，C表达为64；直接用差B相对A就上调了7，C就相对B上调了56；用log2 fold change，B相对A就上调了3，C相对B也只上调了3. 通过测序观察我们发现，不同基因在细胞里的表达差异非常巨大，所以直接用差显然不合适，用log2 fold change更能表示相对的变化趋势。

虽然大家都在用log2 fold change，但显然也是有缺点的：一、到底是5到10的变化大，还是100到120的变化大？二、5到10可能是由于技术误差导致的。所以当基因总的表达值很低时，log2 fold change的可信度就低了，尤其是在接近0的时候。A disadvantage and serious risk of using fold change in this setting is that it is biased[7] and may misclassify differentially expressed genes with large differences (B − A) but small ratios (B/A), leading to poor identification of changes at high expression levels. Furthermore, when the denominator is close to zero, the ratio is not stable, and the fold change value can be disproportionately affected by measurement noise.

这就是统计学的范畴了，显著性就是根据假设检验算出来的。

假设检验首先必须要有假设，我们假设A和B的表达没有差异（H0，零假设），然后基于此假设，通过t test（以RT-PCR为例）算出我们观测到的A和B出现的概率，就得到了P-value，如果P-value<0.05，那么说明小概率事件出现了，我们应该拒绝零假设，即A和B的表达不一样，即有显著差异。

显著性只能说明我们的数据之间具有统计学上的显著性，要看上调下调必须回去看差异倍数。

这里只说了最基本的原理，真正的DESeq2等工具里面的算法肯定要复杂得多。

这张图对q-value（校正了的p-value）取了负log，相当于越显著，负log就越大，所以在火山图里，越外层的岩浆就越显著，差异也就越大。

只需要看懂DEG结果的可以就此止步，想深入了解的可以继续。

另一篇关于建库的文章： RNA-seq建库技术 | RNA sequence library construction

下面会讨论的问题有：

做生物生理生化生信数据分析时，最常听到的肯定是“差异(表达)基因分析”了，从最开始的RT-PCR，到基因芯片microarray，再到RNA-seq，最后到现在的single cell RNA-seq，统统都在围绕着差异表达基因做文章。

（开个脑洞：再下一步应该会测细胞内特定空间内特定基因的动态表达水平了）

表达量 ：我们假设基因转录表达形成的mRNA的数量反映了基因的活性，也会影响下游蛋白和代谢物的变化。我们关注的是基因的表达，不是结构，也是不是isoform。

为什么差异基因分析这么流行？一是中心法则得到了确立，基因表达是核心的一个环节，决定了下游的蛋白组和代谢组；二是建库测序的普及，获取基因的表达水平变得容易。

在生物体内，基因的表达时刻都在动态变化，不一定服从均匀分布，在不同时间、发育程度、组织和环境刺激下，基因的表达肯定会发生变化。

差异基因分析主要应用在：

目前我们对基因和转录组的了解到什么程度了？

基本的建库方法？建库直接决定了我们能测到什么序列，也决定了我们能做什么分析！

基因表达的normalization方法有哪些？

第一类错误、第二类错误是什么？

多重检验的校正？FDR

10x流程解释

The mean UMI counts per cell of this gene in cluster I
The log2 fold-change of this gene\'s expression in cluster i relative to other clusters
The p-value denoting significance of this gene\'s expression in cluster i relative to other clusters, adjusted to account for the number of hypotheses (i.e. genes) being tested.

The differential expression analysis seeks to find, for each cluster, genes that are more highly expressed in that cluster relative to the rest of the sample. Here a differential expression test was performed between each cluster and the rest of the sample for each gene.

The Log2 fold-change (L2FC) is an estimate of the log2 ratio of expression in a cluster to that in all other cells. A value of 1.0 indicates 2-fold greater expression in the cluster of interest.

The p-value is a measure of the statistical significance of the expression difference and is based on a negative binomial test. The p-value reported here has been adjusted for multiple testing via the Benjamini-Hochberg procedure.

In this table you can click on a column to sort by that value. Also, in this table genes were filtered by (Mean UMI counts > 1.0) and the top N genes by L2FC for each cluster were retained. Genes with L2FC < 0 or adjusted p-value >= 0.10 were grayed out. The number of top genes shown per cluster, N, is set to limit the number of table entries shown to 10000; N=10000/K^2 where K is the number of clusters. N can range from 1 to 50. For the full table, please refer to the "differential_expression.csv" files produced by the pipeline.

不同单细胞DEG鉴定工具的比较

Comparative analysis of differential gene expression analysis tools for single-cell RNA sequencing data

For data with a high level of multimodality, methods that consider the behavior of each individual gene, such as DESeq2, EMDomics, Monocle2, DEsingle, and SigEMD, show better TPRs. 这些工具敏感性高，就是说不会漏掉很多真的DEG，但是会包含很多假的DEG。
If the level of multimodality is low, however, SCDE, MAST, and edgeR can provide higher precision. 这些工具精准性很高，意味着得到的DEG里假的很少，所以会漏掉很多真的DEG，不会引入假的DEG。

time-course DEG analysis

Comparative analysis of differential gene expression tools for RNA sequencing time course data

参考：

Question: How to calculate "fold changes" in gene expression?

Exact Negative Binomial Test with edgeR

Differential gene expression analysis

相关文章：

ggplot的boxplot添加显著性 | Add P-values and Significance Levels to ggplots | 方差分析

差异表达edgeR，limma（下）

参考技术A 在我们看来，用于检查基因表达分析的最重要的探索性图表之一便是MDS图或其余类似的图。这种图表使用无监督聚类方法展示出了样品间的相似性和不相似性，能让我们在进行正式的检验之前对于能检测到多少差异表达基因有个大致概念。理想情况下，样本会在不同的实验组内很好的聚类，且可以鉴别出远离所属组的样本，并追踪误差或额外方差的来源。如果存在技术重复，它们应当互相非常接近。

这样的图可以用 limma 中的 plotMDS 函数绘制。第一个维度表示能够最好地分离样品且解释最大比例的方差的引导性的倍数变化（leading-fold-change），而后续的维度的影响更小，并与之前的维度正交。当实验设计涉及到多个因子时，建议在多个维度上检查每个因子。如果在其中一些维度上样本可按照某因子聚类，这说明该因子对于表达差异有影响，在线性模型中应当将其包括进去。反之，没有或者仅有微小影响的因子在下游分析时应当被剔除。

主要说明一下 edgeR 中的 glmQLFTest , exactTest 以及 limma 中的 voom 这几种获取差异基因的不同方式

Limma包基于线性模型建模。 它最初设计用于分析微阵列数据，但最近已扩展到RNA-seq数据。 根据limma用户指南的当前建议是使用edgeR包的TMM标准化和“voom”转换，其本质上将标准化数据取对数并估计它们的均值 - 方差关系以确定在线性建模之前每次观察的权重。 默认情况下，Benjamini-Hochberg程序用于估计FDR 。

首先先建立design矩阵,在此研究中，我们想知道哪些基因在我们研究的两组之间以不同水平表达。在我们的分析中，假设基础数据是正态分布的，为其拟合一个线性模型。在此之前，需要创建一个包含细胞类型design矩阵。

据显示对于RNA-seq计数数据而言，当使用原始计数或当其被转换为log-CPM值时，方差并不独立于均值。使用负二项分布来模拟计数的方法假设均值与方差间具有二次的关系。在limma中，假设log-CPM值符合正态分布，并使用由 voom 函数计算得到的精确权重来调整均值与方差的关系，从而对log-CPM值进行线性建模。

当操作DGEList对象时，voom从x中自动提取文库大小和归一化因子，以此将原始计数转换为log-CPM值。在 voom 中，对于log-CPM值额外的归一化可以通过设定normalize.method参数来进行。

下图左侧展示了这个数据集log-CPM值的均值-方差关系。通常而言，方差是测序实验中的技术差异和不同细胞类型的重复样本之间的生物学差异的结合，而voom图会显示出一个在均值与方差之间递减的趋势。 生物学差异高的实验通常会有更平坦的趋势，其方差值在高表达处稳定。 生物学差异低的实验更倾向于急剧下降的趋势。

不仅如此，voom图也提供了对于上游所进行的过滤水平的可视化检测。如果对于低表达基因的过滤不够充分，在图上表达低的一端，受到非常低的表达计数的影响，可以观察到方差水平的下降。如果观察到了这种情况，应当回到最初的过滤步骤并提高用于该数据集的表达阈值。

edgeR使用经验贝叶斯估计和基于负二项模型的精确检验来确定差异基因。 特别地，经验贝叶斯用于通过在基因之间来调节跨基因的过度离散程度。 使用类似于Fisher精确检验但适应过度分散数据的精确检验用于评估每个基因的差异表达。edgeR 在默认情况下，执行TMM标准化程序以考虑样本之间的不同测序深度，通过Benjamini-Hochberg用于控制FDR 。

精确检验适用于多组实验的精确统计法，使用函数 exactTest 估计差异基因，人们将其称为classic edgeR。

estimateDisp 函数在一组离散网格点上为每个标签(基因)计算一个似然矩阵，然后应用加权似然经验贝叶斯方法获得后验离散度估计。如果没有设计矩阵，它计算每个标签的分位数的条件似然，然后将其最大化。在这种情况下，它类似于函数 estimateCommonDisp 和 estimateTagwiseDisp 。

同样首先利用 calcNormFactors 对因子进行矫正

似然比检验是用广义线性模型（glms）的统计学方法，适用于不同复杂程度的多因素实验，而 QLFTest 则是首选，因为它反映了估计每个基因离散度的不确定性。当重复次数较少时，它可以提供更可靠的错误率控制。