【差异表达分析】差异表达分析标准不标准化这是一个问题（含其其它报错问题）

卡西莫多的猫

已于 2023-03-10 19:43:19 修改

阅读量3k

点赞数 6

分类专栏：生信那些坑文章标签：算法

于 2023-03-10 16:13:27 首次发布

本文链接：https://blog.csdn.net/ohmmmmmypyr/article/details/129444059

版权

生信那些坑专栏收录该内容

1 篇文章 1 订阅

订阅专栏

前言

在一开始学习基因差异表达分析时，老师就强调用raw count做差异分析，相关文献和资料我也保存了不少，我之前弄清楚log2/cpm与count fpkm等不是在一个水平上讨论的问题，但是具体用的时候还是要栽个跟头才能印象深刻。

我在复现这篇推文时老文新看，今天来看看两个数据集的整合分析
一般情况下我自以为已经掌握一般差异分析流程，就看看把作者代码跑一边就行了，在这篇推文中作者没有给出差异分析部分代码和数据，但是后面介绍并使用了一种我很感兴趣的RRA算法，所以自己动手分别用DESeq2和limma做了差异分析，结果一动手问题就出来了，两个包的差异分析结果都表明没有显著差异表达的基因！我一开始甚至认为是是数据集的问题，但昨天恰好在讨论群里也有一位uu有相同的问题，我跟他交换了代码才发现了log2这个端倪。

uu的流程

uu用的是是GSE26728芯片数据集，我用我的流程跑一一遍他的数据，起初发现跟他结果一致，没有什么差异表达的基因，但是我的结果数值都比他大一点。查看他自己的代码，发现他在exp表达矩阵range在20以内的情况下仍然进行了log2，所以导致数据比我小（其实不是，是因为直接log2和使用limma::voom算法的不同导致的微小差异，后面会提）。

但是我同时注意到，他在使用limma包进行差异分析时直接使用了exp表达矩阵进行（这里面他的数据框格式下的表达矩阵有一个小坑）

fit=lmFit(exp1,design)
# Error in lmFit(exp1, design) : 
#   Expression object should be numeric, instead it is a data.frame with 11 non-numeric columns
# 我那里这个是DGEList voom后 EList对象
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)

（小坑，他在处理数据框时里面有许多描述信息，属于字符元素，而差异分析需要numeric数值型元素，如果你不把所有字符型元素去除就会发现为了统一所有的元素类型都变成了字符型
在这里插入图片描述）
但是当我将他的log2处代码注释掉后再次进行发现差异分析结果变得“正常了”，而不是预期中与我一致

所以我开始怀疑我的流程中存在我不知的log2或者说标准化操作

我的流程

我查看我的流程发现我在这部分先将exp表达矩阵转换成了一个DEGList对象后使用voom转换成了一个EList对象
在这里插入图片描述

一开始我以为voom函数是线性化矩阵以便后面拟合，而当我再次查看voom函数的功能时发现voom函数同样对表达矩阵进行了log2标准化处理
在这里插入图片描述

真相大白

原来我和他都犯了同一个错误，那就是重复标准化，导致差异也被缩小，换句话讲就是，我和他都使用了标准化后的数据(cpm/log2后数据常用来作图)，而非rawcount进行了一般情况下的差异表达分析流程！
所以我认为，针对这个错误，有两种解决办法，一是将我们一般性流程中的标准化部分去除，如uu那样直接使用表达矩阵，而不通过limma::voom
但这种方法存在局限性:只能适用于limma包，因为limma包voom log2步骤是分开的，有操作空间，而DESeq2 result一体化

在这里插入图片描述
limma这种分布进行的特点在复现这篇推文，差异分析不是这样做的……
时也可以发现，论文作者的失误正是使用limma时，想当然地认为不对表达矩阵取log2就得到FC，实际上

limma包需要接收log2后的表达矩阵(voom/log2)
算法核心是logA-logB= log(A-B) 作者实际操作是rawcount A-B
不能想当然认为不对表达矩阵取log2就得到FC 如果要直接从rawcount入手也要是FC = B/A
如果作者就是想直接拿rawcount走到差异分析这一步，可以使用DESeq2，result后自动得到的虽然也是log2FC，手动2^即可

在实际学习的过程中，
生信菜鸟团：：转录组专题的老师使用DESeq2和limma包时都使用的时rawcount，因为这个专辑涉及到上游分析的学习，limma于是就需要voom
而表达量芯片的老师偏向使用limma包处理log2后的exp矩阵，可能因为芯片的数据大多数GEO上存的都是log2后的，如果是rawcount也还是需要手动log2/voom

其它需要注意的标准化

组间差异

在进行基因差异分析前需要考虑组间差异，如通过limma::normalizeBetweenArrays去除

进行基因差异分析前通常需要考虑组间差异，因为在组间差异显著的情况下，比较基因表达水平的差异才具有生物学意义。

在使用limma包进行基因差异分析时，可以使用其中的normalizeBetweenArrays函数对不同组之间的差异进行标准化处理。该函数可以用于批次效应的去除和样本间的标准化，使得样本之间的比较更加准确和可靠。但需要注意的是，标准化只是一个预处理步骤，对于基因差异分析的结果可能会有一定的影响，因此还需要根据实验的设计和数据的质量进行合理的统计分析。

TODO

待补充

未进行组间差异处理导致的·DESeq2报错

dds2 <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
Error in estimateDispersionsFit(object, fitType = fitType, quiet = quiet) :
all gene-wise dispersion estimates are within 2 orders of magnitude
from the minimum value, and so the standard curve fitting techniques will not work.
One can instead use the gene-wise estimates as final estimates:
dds <- estimateDispersionsGeneEst(dds)
dispersions(dds) <- mcols(dds)$dispGeneEst
…then continue with testing using nbinomWaldTest or nbinomLRT

chatgpt：
这个错误提示表明在使用DESeq对dds对象进行差异分析时，估计基因的分散度时出现了问题。具体来说，这个错误通常是由于基因的分散度估计值非常相近而无法进行标准拟合所导致的。

解决这个问题的方法是使用estimateDispersionsGeneEst函数，该函数可以使用基因特定的分散度估计值来代替标准拟合。具体而言，需要按照错误提示中所示的步骤进行操作：

dds <- estimateDispersionsGeneEst(dds)
dispersions(dds) <- mcols(dds)$dispGeneEst

这将使用基因特定的分散度估计值进行差异分析，可以使用nbinomWaldTest或nbinomLRT函数进行差异分析测试。