基因組的保守性常常與功能有關(guān)，保守性強(qiáng)的序列可能在細(xì)胞的發(fā)育和調(diào)控等方面發(fā)揮著至關(guān)重要的作用。分享一下一般的分析保守型的方法：
首先要下載保守性分值的文件，可以在UCSC中選擇對應(yīng)的物種的文件，比如小鼠mm10的phyloP conservation文件(http://hgdownload.cse.ucsc.edu/goldenpath/mm10/phyloP60way/)或者phastCons文件(http://hgdownload.cse.ucsc.edu/goldenpath/mm10/phastCons60way/mm10.60way.phastCons.bw)。
接下來我們用bwtool提取出來區(qū)域的分值均值就可以了。bwtool是2014年發(fā)表在Bioinformatics發(fā)表的的一個工具。
安裝：
git clone https://github.com/CRG-Barcelona/libbeato.git
git clone https://github.com/CRG-Barcelona/bwtool.git
cd libbeato/
./configure --prefix=$HOME CFLAGS="-g -O0 -I${HOME}/include" LDFLAGS=-L${HOME}/lib
make
make install
cd ../bwtool/
./configure --prefix=$HOME CFLAGS="-g -O0 -I${HOME}/include" LDFLAGS=-L${HOME}/lib
make
make install
這是Manuel中的安裝方法，不過我安裝了好久一直失敗。實在不行可以用conda安裝：
conda install -c pwwang bwtool
接下來可以用用bwtool agg提出來+/-2kb的平均分?jǐn)?shù)：
bwtool agg 2000:2000 gencode_tss.bed mm10.60way.phyloP60way.bw gencode_score.txt
注意這個bed文件，當(dāng)提供的是BED6格式時，第六列是正負(fù)鏈的信息，當(dāng)是負(fù)鏈時，數(shù)據(jù)會被反轉(zhuǎn)（與正鏈上下游是相反的）。我們一般是需要這么處理的。如果只提供三列，染色體加位置的話，則不能判斷正負(fù)鏈，會導(dǎo)致錯誤。
最后我們可以得到一個每個位置的分?jǐn)?shù)的文件，ggplot2作圖即可：
ggplot(data = new_merged_score, aes(x = distance, y = score, color = group)) + \
geom_line(size = 2) + \
theme_classic() + \
theme(legend.position = 'top') + \
guides(color = guide_legend(title = NULL, override.aes = list(size = 1))) + \
labs(x = '', y = 'mm10 60-way phyloP score')
歡迎關(guān)注～
