寫在前面的話

你既然點(diǎn)進(jìn)來(lái)了，說(shuō)明你多半用過(guò)這個(gè)軟件。但我猜，大多數(shù)人早期使用這個(gè)軟件的時(shí)候，一定是趕鴨子上架：因?yàn)橐媹D，所以必須會(huì)。

但你真的懂了么？

太長(zhǎng)不看系列

對(duì)于TSS到TES之間的區(qū)域，computeMatrix將不同基因的長(zhǎng)度均一化到同樣的長(zhǎng)度。而對(duì)于TSS上游或者下游，則無(wú)需均一化。

那么如何均一化到同樣的長(zhǎng)度呢？請(qǐng)細(xì)細(xì)往下看

廢話超多系列

我們?cè)谶@里以scale-region為例進(jìn)行分析

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

經(jīng)常處理ChIP-seq數(shù)據(jù)的同學(xué)，應(yīng)該很快就能明白這條命令的具體含義。當(dāng)然不懂得的話我這里也簡(jiǎn)單介紹下命令的意思。

大概就是計(jì)算繪制ChIP-seq的信號(hào)趨勢(shì)圖，或者說(shuō)是metaplot所需要的矩陣文件。test.bw文件中存儲(chǔ)的是ChIP-seq的信號(hào)，而test.bed存儲(chǔ)的基因的坐標(biāo)。信號(hào)起始位點(diǎn)取轉(zhuǎn)錄前起始位點(diǎn)2000bp，信號(hào)終止位點(diǎn)取轉(zhuǎn)錄起始位點(diǎn)之后2000bp，基因body區(qū)域設(shè)置為5000bp。如果我們使用的是H3K4me3的ChIP-seq數(shù)據(jù)，那么我們使用plotprofile繪圖得到的結(jié)果將如下：

image.png

到這里大家也許會(huì)產(chǎn)生一個(gè)疑惑，起碼我產(chǎn)生了一個(gè)疑惑：不同的基因具有不同的長(zhǎng)度，computeMatrix是如何做到把不同長(zhǎng)度的基因均一化到相同長(zhǎng)度呢？

首先，我們假設(shè)有存在一個(gè)geneA，他的長(zhǎng)度是2000bp。按照剛剛上面的命令，我們將binsize設(shè)置成了50bp，也就說(shuō)把基因分成若干塊，每一塊的長(zhǎng)度均是50bp。

這樣的話，長(zhǎng)度為2000bp的geneA被切成了40塊。隨著這基因的長(zhǎng)度分bin，每一個(gè)bin的取值也分別對(duì)應(yīng)著這個(gè)bin所在區(qū)域，所有堿基信號(hào)加和然后取平均的值。如果這個(gè)bin所在的堿基信號(hào)是1，2，3，4，……50，那么最后得到的該bin所對(duì)應(yīng)的信號(hào)值則是(1+2+3+4+……50)/50=25.5

這個(gè)地方bin對(duì)應(yīng)的信號(hào)是所有堿基加和的平均值，最早我以為是信號(hào)的總和(sum)……雖然看起來(lái)似乎沒(méi)有什么不同

那么如果此時(shí)存在一個(gè)長(zhǎng)度為4000bp的geneB，該怎么辦呢？剛剛bin設(shè)置為50bp后，geneA分成了40塊，那么為了對(duì)應(yīng)，我們就應(yīng)該也將geneB分成40塊，于是geneB的binsize就成了4000/40=100bp。而每一個(gè)bin對(duì)應(yīng)的score也可以按照剛剛對(duì)geneA的處理進(jìn)行計(jì)算。

這樣處理之后，我們就可以將不同長(zhǎng)度的基因分成相同份數(shù)的bin，從而將他們統(tǒng)統(tǒng)均一化到相同的長(zhǎng)度。至于我們?nèi)绾未_定應(yīng)該分為多少個(gè)bin呢？

在使用computeMatrix時(shí)，我們定義了要normalized后的基因長(zhǎng)度(--regionBodyLength)。上面使用的是5000，那么我們就將所有的基因都分成5000/50=100塊。

根據(jù)在上述原理，我們就可以自己去計(jì)算得分矩陣，也就是computeMatrix所得到的matrix.gz文件。當(dāng)然了，我不建議大家自己去計(jì)算，因?yàn)槲易约鹤鲞^(guò)類似的事情，費(fèi)了老大力氣改bug，后來(lái)發(fā)現(xiàn)速度還不如computeMatrix，何必呢？

image.png

雖然不建議仿照computeMatrix造輪子，但是我們可以利用得到的打分矩陣自己畫圖啊。我是實(shí)在無(wú)法忍受，畫圖的各種參數(shù)不能受我支配的感覺(jué)。

那么簡(jiǎn)單介紹一下最后生成的得分矩陣的格式，也方便大家以后自己使用ggplot2繪制metaplot。信息記錄大概是下圖這個(gè)樣子：

image.png

第一行呢，主要是一些注釋信息，不太重要，一般讀取數(shù)據(jù)的時(shí)候可以直接忽略。從第二行開(kāi)始就是每個(gè)基因的位置信息，和其所有bin值所對(duì)應(yīng)的得分信息。從第一列到第六列，對(duì)應(yīng)的信息分別是基因的染色體信息，起始位置，終止位置，基因名，打分(均為0)，和鏈方向。

從第七列開(kāi)始就是打分信息，如果我們只有一個(gè)bw文件，并切成了100塊，那么這個(gè)得分矩陣從第七列到第106列對(duì)應(yīng)著每個(gè)基因不同bin的得分。如果有兩個(gè)，則第二個(gè)bw文件的信息是從第107列開(kāi)始，到206列結(jié)束。

利用上述內(nèi)容，我相信你多半可以自己使用ggplot2繪制metaplot了。

順便說(shuō)一句，無(wú)論你是用正鏈還是負(fù)鏈計(jì)算的得分矩陣，最后繪制的圖都是從左往右方向，即起始位點(diǎn)在左，終止位點(diǎn)在右

寫在后面的話

如果你還是不知道怎么畫圖，那你可以聯(lián)系我，我考慮發(fā)一下代碼……

說(shuō)實(shí)話，代碼其實(shí)還是手寫的踏實(shí)