保守性分析并作图

基因组的保守性常常与功能有关，保守性强的序列可能在细胞的发育和调控等方面发挥着至关重要的作用。分享一下一般的分析保守型的方法：

首先要下载保守性分值的文件，可以在UCSC中选择对应的物种的文件，比如小鼠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提取出来区域的分值均值就可以了。bwtool是2014年发表在Bioinformatics发表的的一个工具。

安装：

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的平均分数：

bwtool agg 2000:2000 gencode_tss.bed mm10.60way.phyloP60way.bw  gencode_score.txt

注意这个bed文件，当提供的是BED6格式时，第六列是正负链的信息，当是负链时，数据会被反转（与正链上下游是相反的）。我们一般是需要这么处理的。如果只提供三列，染色体加位置的话，则不能判断正负链，会导致错误。

最后我们可以得到一个每个位置的分数的文件，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')

欢迎关注～