前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >保守性分析并作图

保守性分析并作图

作者头像
生信编程日常
发布2020-07-16 17:19:03
1.8K0
发布2020-07-16 17:19:03
举报
文章被收录于专栏:生物信息学、python、R、linux

基因组的保守性常常与功能有关,保守性强的序列可能在细胞的发育和调控等方面发挥着至关重要的作用。

首先要下载保守性分值的文件,可以在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发表的的一个工具。

安装:

代码语言:javascript
复制
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安装:

代码语言:javascript
复制
conda install -c pwwang bwtool

接下来可以用用bwtool agg提出来+/-2kb的平均分数:

代码语言:javascript
复制
bwtool agg 2000:2000 gencode_tss.bed mm10.60way.phyloP60way.bw  gencode_score.txt

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

最后我们可以得到一个每个位置的分数的文件,ggplot2作图即可:

代码语言:javascript
复制
ggplot(data = gencode_pc, aes(x = V1, y = V2)) + geom_line(color = 'steelblue') + theme_classic()

第一列是位置,第二列是分数。

本文参与 腾讯云自媒体同步曝光计划,分享自作者个人站点/博客。
如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 作者个人站点/博客 前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档