你考虑过computeMatrix的工作原理么?

写在前面的话

你既然点进来了，说明你多半用过这个软件。但我猜，大多数人早期使用这个软件的时候，一定是赶鸭子上架：因为要画图，所以必须会。

但你真的懂了么？

太长不看系列

对于TSS到TES之间的区域，computeMatrix将不同基因的长度均一化到同样的长度。而对于TSS上游或者下游，则无需均一化。

那么如何均一化到同样的长度呢？请细细往下看

废话超多系列

我们在这里以scale-region为例进行分析

computeMatrix scale-regions -S test.bw -R test.bed --regionBodyLength 5000 -b 2000 -a 2000 -bs 50 -p 24 --skipZeros -o ./test.matrix.gz

经常处理ChIP-seq数据的同学，应该很快就能明白这条命令的具体含义。当然不懂得的话我这里也简单介绍下命令的意思。

大概就是计算绘制ChIP-seq的信号趋势图，或者说是metaplot所需要的矩阵文件。test.bw文件中存储的是ChIP-seq的信号，而test.bed存储的基因的坐标。信号起始位点取转录前起始位点2000bp，信号终止位点取转录起始位点之后2000bp，基因body区域设置为5000bp。如果我们使用的是H3K4me3的ChIP-seq数据，那么我们使用plotprofile绘图得到的结果将如下：

image.png

到这里大家也许会产生一个疑惑，起码我产生了一个疑惑：不同的基因具有不同的长度，computeMatrix是如何做到把不同长度的基因均一化到相同长度呢？

首先，我们假设有存在一个geneA，他的长度是2000bp。按照刚刚上面的命令，我们将binsize设置成了50bp，也就说把基因分成若干块，每一块的长度均是50bp。

这样的话，长度为2000bp的geneA被切成了40块。随着这基因的长度分bin，每一个bin的取值也分别对应着这个bin所在区域，所有碱基信号加和然后取平均的值。如果这个bin所在的碱基信号是1，2，3，4，……50，那么最后得到的该bin所对应的信号值则是(1+2+3+4+……50)/50=25.5

这个地方bin对应的信号是所有碱基加和的平均值，最早我以为是信号的总和(sum)……虽然看起来似乎没有什么不同

那么如果此时存在一个长度为4000bp的geneB，该怎么办呢？刚刚bin设置为50bp后，geneA分成了40块，那么为了对应，我们就应该也将geneB分成40块，于是geneB的binsize就成了4000/40=100bp。而每一个bin对应的score也可以按照刚刚对geneA的处理进行计算。

这样处理之后，我们就可以将不同长度的基因分成相同份数的bin，从而将他们统统均一化到相同的长度。至于我们如何确定应该分为多少个bin呢？

在使用computeMatrix时，我们定义了要normalized后的基因长度(--regionBodyLength)。上面使用的是5000，那么我们就将所有的基因都分成5000/50=100块。

根据在上述原理，我们就可以自己去计算得分矩阵，也就是computeMatrix所得到的matrix.gz文件。当然了，我不建议大家自己去计算，因为我自己做过类似的事情，费了老大力气改bug，后来发现速度还不如computeMatrix，何必呢？

image.png

虽然不建议仿照computeMatrix造轮子，但是我们可以利用得到的打分矩阵自己画图啊。我是实在无法忍受，画图的各种参数不能受我支配的感觉。

那么简单介绍一下最后生成的得分矩阵的格式，也方便大家以后自己使用ggplot2绘制metaplot。信息记录大概是下图这个样子：

image.png

第一行呢，主要是一些注释信息，不太重要，一般读取数据的时候可以直接忽略。从第二行开始就是每个基因的位置信息，和其所有bin值所对应的得分信息。从第一列到第六列，对应的信息分别是基因的染色体信息，起始位置，终止位置，基因名，打分(均为0)，和链方向。

从第七列开始就是打分信息，如果我们只有一个bw文件，并切成了100块，那么这个得分矩阵从第七列到第106列对应着每个基因不同bin的得分。如果有两个，则第二个bw文件的信息是从第107列开始，到206列结束。

利用上述内容，我相信你多半可以自己使用ggplot2绘制metaplot了。

顺便说一句，无论你是用正链还是负链计算的得分矩阵，最后绘制的图都是从左往右方向，即起始位点在左，终止位点在右

写在后面的话

如果你还是不知道怎么画图，那你可以联系我，我考虑发一下代码……

说实话，代码其实还是手写的踏实