数据normalize和批次的思考

数据normalize和批次的思考

内容目录

前言数据normalizeFPKM和TPMnormalizeTCGA和GEO/GTEx数据一些思考批次效应关于批次效应处理批次效应的工具如何去批次的差异小结后记

前言

最近想给pubmed做一次浇灌，所以开展了一次利用TCGA数据和GEO数据进行的灌溉计划。之前总觉得这种事情应该比较简单吧，的确来说，看那些文章非常的简单（可能这些人把其中的艰苦给省略了，也可能是他们选择忽略了其中难以处理的问题），但是自己做的时候遇到了很多问题。今天这里只是想简单聊一聊关于normalize和批次效应的问题。

数据normalize

FPKM和TPM

什么是normalize呢？我查了百度，感觉没有什么文章写的比较清楚，或者说没有什么说明适合我们这种不是很懂生信的人来看，为了简单理解什么是normalize，我这里就写上我对normalize过程的理解吧。

首先推荐一个视频教程：

https://www.bilibili.com/video/BV1T4411T7EL

这里的FPKM/RPKM和TPM实际上就是一种normalize，为了让数据之间具有可比性。为了更简单的去理解这个概念，让我们回到生活中的例子：

我们知道考试中，语文数学的总分可能是150，而体育可能只有60分总分。那么，这时语文考了60分和体育考了60分，这两个60分具有可比性吗？

当然是没有可比性的！那么我们这么来比较这个个体是语文更好还是体育更好呢？现实中我们可能会比较百分比：

语文——60/150=40%

体育——60/60=100%

所以我们认为这个同学的体育成绩更好！

类比到测序中，A基因可能长10kb，而B基因可能长100bp，那么A基因的raw reads有50和B基因的raw reads有50，这两者之间也是没有可比性的，因为A基因更长，所以肯定来说A测到的概率更高！

同样的，我们知道在测序过程中，即使是同一个文库，我一个测了100x，一个测了10x，那么这两个数据之间肯定也是相差较大的！

下面用一个具体例子来说明吧：

上面这个例子我写的非常理化，我假设有3个基因，长度之比为3:2:1，那么自然而然，对于个体来说每个基因的读数可能就是3:2:1。然后对于A个体来说，A的测序深度可能比较深，一共测了60个reads，而B的测序深度可能比较浅，一共只测了6个reads。

根据上面的表格，如果我问：A和B之间基因表达有差异吗？

结果是没有的！！

虽然30和3相差很大，30和10也相差很大，但是这些差别是来自基因长度和测序深度的，而非真正的生物学差异！

为了可以比较，就引入了FPKM和TPM：

• FPKM是矫正了测序深度和基因长度的值
• TPM是FPKM的百分比

具体公式我就不写了，大家可以百度到。还是以上述表格为例，我们计算下相关的FPKM值（我没有计算真正的具体值，我只是矫正了长度和深度，具体值可能还需在数据后加若干个0）：

可以看到所有的值都是0.167，所以这些个体之间基因没有差异！！

如果我们转为TPM的话：

可以看到所有的值都是0.33，所以这些个体之间基因没有差异！！

normalizeTCGA和GEO/GTEx数据

因为我现在想要做的事是对LAML进行分析，我找了找Pubmed文章，发现LAML数据集被人挖掘的不是特别多吧（相对于其他数据集，如肝癌、胃肠道肿瘤以及乳腺癌之类的），我猜原因之一是——LAML没有正常对照样本！！

关键是吧，其他project中，如果TCGA里没有正常对照，可以去GTEx里去找，因为GTEx官网里的数据都是正常人的数据，放上它的官网地址：

https://www.gtexportal.org/home/samplingSitePage

然后我就去TGTEx里找了找，居然没有！！因为LAML里用的样本来自骨髓，所以作为对照，我也应该找GTEx里的骨髓样本才对，但是，它没有！！！这就呵呵了！

果然是，天将降大灾大难于你，就会一直劳你筋苦！于是，只能去GEO里找了，然后找到了几个数据集可用！！哈哈哈哈哈哈哈~

那么接下来问题就来了！怎么进行合并分析呢？？

当然是要normalize啦！！最简单的当然是直接算TPM/FPKM啦！不过，我没有直接算TPM/FPKM，我用了TCGAbiolinks这个R包来进行的normalize。

这里我就给出思路吧，虽然我已经做好了，但是因为有些问题，所以我就不给代码了。

但是当我对每个样本单独进行normalize后，我常规对矩阵进行boxplot来看数据情况，发现了一个bug问题：

先看看TCGA内部数据情况：

看上去还不错，内部之间没有什么太大的差异，然后我们看看GEO数据的情况：

看上去也还不错，内部之间没有什么太大的差异，最后我们看看直接cbind后的数据情况：

这就呵呵了。目前我认为这个属于批次效应吧~

一些思考

关于上面合并后的数据问题，我目前可能认为是一种批次效应，打算用combat来矫正下，具体情况等我验证了后再看吧。

同时，在和其他人讨论后，也引出了一些其他思考，例如：

1. normalize数据的顺序问题？大家认为我们是应该先对每个数据单独normalize，然后用cbind合并分析呢？还是应该先把raw counts直接cbind合并后，再进行normalize呢？（如果是利用FFPKM/TPM的话其实无所谓，因为这个是根据每一列来计算的，和行关系不大）。
2. 我们如何判断我们的数据到底有没有批次效应呢？用PCA聚类？画层次聚类图？画boxplot？如果有的话如何矫正呢？

关于上面提到的问题，我近期都会尝试下，如果有精力的话也许会记录下结果，如果没有精力的话，可能有空了会放个结论，让大家知道以后怎么做会比较好吧！

批次效应

放上几个博客：

https://www.plob.org/article/14410.html

https://www.dazhuanlan.com/2019/11/15/5dce4a071589c/

搞定TCGA批次效应！?

关于批次效应

2010年的时候10年nature有一篇综述，专门讲这个问题：

https://www.nature.com/articles/nrg2825.pdf

我只是，先下载了，还没看，这几天有空了我会读一读，但是思想上我大致有个了解。

其实，我们做normalize一定程度上也可以理解为在去批次，去除的是每次测序深度的批次，以及每个基因本身属性的批次。

处理批次效应的工具

批次效应不能被消除，只有尽可能的降低。你可能之前或多或少了解到数据中可能会存在批次效应，或者已经使用 removeBatchEffect /ComBat 等工具做过批次矫正。那么就涉及到一个经典的问题，这两个工具该如何选择？

1. 两者都是针对芯片（强调！！！芯片是芯片，测序是测序）数据。所以在使用前，数据都要经过一定的标准化操作，如log转化。并不能直接应用于 read counts！
2. 经典工具 ComBat 的原作者又开发了ComBat-Seq，针对测序数据，以保证经过批次效应矫正后得到的仍是整数counts，便于后续衔接基于整数进行差异分析的工具，如 edgeR 和 DESeq2 等。具体网址在：https://github.com/zhangyuqing/ComBat-seq
3. DESeq2 包官方推荐使用 limma 包中的 removeBatchEffect() 去批次。

如何去批次的差异

不管是removeBatchEffect还是ComBat都是直接对原数据进行修改，意味着你的后续分析要基于矫正后的数据进行。那么是不是意味着，对批次进行矫正后，就刚好能拿来做差异表达分析了呢？

不是的！！

removeBatchEffect() 函数的文档中有明确的叙述: removeBatchEffect 只能用于衔接聚类、PCA等可视化展示，不能将去批次后的数据用于差异分析！

This function is intended for use with clustering or PCA, not for use prior to linear modelling. If linear modelling is intended, it is better to include the batch effect as part of the linear model.

因为用矫正后的数据进行差异表达分析 (已矫正批次，所以模型只考虑分组因素)，有两个缺陷：

1、批次因素和分组因素可能重叠，所以直接对原数据矫正批次可能会抵消一部分真实生物学因素。

2、低估了误差。

所以，如果想做差异表达分析，但数据中又有已知的批次问题，则应该在构建模型矩阵时加入批次因素！！！

小结

总之，批次效应的处理要依分析目的而定。如果是做差异表达分析，则应该在模型中添加批次因素，而非使用矫正后的数据。如果是做可视化，则考虑先对原数据进行处理矫正，并使用矫正后的数据进行分析。

后记

之前觉得数据挖掘起来应该还比较简单吧，因为代码什么的还是很容易看懂的，但是做了后才发现，难的是细节！！

这次就这样吧，也是随便写写记录下，后面还是要真正运行下，才知道到底怎么做会比较好！！有精力的话再整理推文吧！