得到一个物种所有基因的TSS(转录起始位点)区域的bed文件。

首先在UCSC的table browser 里面下载下面这个文件:

可以看到我这里选择的mm10的refseq系统的所有基因,共有29037个不同的tss,36872个转录本,只有24540个基因,说明有部分基因有多个tss,这个其实挺麻烦的。

#bin    name    chrom   strand  txStart txEnd   cdsStart    cdsEnd  exonCount   exonStarts  exonEnds    score   name2   cdsStartStat    cdsEndStat  exonFrames0    NM_001282945    chr1    -   134199214   134235457   134202950   134234355   3   134199214,134234014,134235227,  134203590,134234446,134235457,  0   Adora1  cmpl    cmpl    2,0,-1,0    NM_001039510    chr1    -   134199214   134235457   134202950   134234355   3   134199214,134234014,134235227,  134203590,134234412,134235457,  0   Adora1  cmpl    cmpl    2,0,-1,0    NM_001291930    chr1    -   134199214   134235457   134202950   134203505   2   134199214,134235227,    134203590,134235457,    0   Adora1  cmpl    cmpl    0,-1,0    NM_001291928    chr1    -   134199214   134234856   134202950   134234733   2   134199214,134234662,    134203590,134234856,    0   Adora1  cmpl    cmpl    2,0,0    NM_001008533    chr1    -   134199214   134235457   134202950   134234355   2   134199214,134234014,    134203590,134235457,    0   Adora1  cmpl    cmpl    2,0,

其实里面可以设置直接下载所有基因的TSS区域的bed文件,可是我不会设置各种参数,也懒得去摸索,直接对上面的文件我可以写脚本处理得到需要的数据形式。 需要输出的是bed格式文件,如下: chrom / chromStart /chromEnd /name /score /strand 我这里定义的TSS(转录起始位点)区域上下游2.5kb,所以代码如下:

perl -alne '{next if /^#/;if($F[3] eq "+"){$start=$F[4]-2500;$end=$F[4]+2500}else{$start=$F[5]-2500;$end=$F[5]+2500}print join("\t",$F[2],$start,$end,$F[12],0,$F[3])}' ucsc.refseq.txt |sort -u >ucsc.refseq.tss.bed

最后得到的文件如下:

tail ucsc.refseq.tss.bed chrY    816212  821212  Uba1y   0   +chrY    81798997    81803997    Gm20747 0   -chrY    82222714    82227714    Gm20736 0   +chrY    83925411    83930411    Gm20854 0   -chrY    85527019    85532019    Gm20854 0   -chrY    8832669 8837669 Gm20815 0   -chrY    895287  900287  Kdm5d   0   +chrY    90752550    90757550    G530011O06Rik   0   -chrY    90782941    90787941    Erdr1   0   +chrY    90836906    90841906    G530011O06Rik   0   

这里面会有一个问题,对于部分基因在非正常染色体的,会出现如下诡异的结果,建议干脆删除掉。

chr4_GL456216_random    13380   18380   Dhrsx   0   +chr4_GL456350_random    -1369   3631    Ccl21c  0   -chr4_GL456350_random    -1369   3631    Gm10591 0   -chr4_GL456350_random    -1369   3631    Gm13304 0   -

记住,这个时候,部分基因还有多个tss哦,反正取决于你的下游分析流程啦。

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-07-26

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏IT派

这一篇就够了 python语音识别指南终极版

【导读】亚马逊的 Alexa 的巨大成功已经证明:在不远的将来,实现一定程度上的语音支持将成为日常科技的基本要求。整合了语音识别的 Python 程序提供了其他...

1100
来自专栏杨熹的专栏

Python 爬虫 2 爬取多页网页

参考资料:极客学院: Python单线程爬虫 代码:2.Single-thread-crawler.ipynb 本文内容: Requests.get 爬取多个页...

3355
来自专栏AI科技大本营的专栏

python语音识别终极指南

译者 | 廉洁 编辑 | 明明 【AI科技大本营导读】亚马逊的 Alexa 的巨大成功已经证明:在不远的将来,实现一定程度上的语音支持将成为日常科技的基本要求。...

4357
来自专栏极客猴

学爬虫之道

Django 已经算是入门,所以自己把学习目标转到爬虫。自己接下来会利用三个月的时间来专攻 Python 爬虫。这几天,我使用“主题阅读方法”阅读 Python...

712
来自专栏肖洒的博客

基于Python的微信好友分析

“如果我比别人看得远,那是因为我站在巨人的肩膀上”–不知道牛顿说了没 本文利用Python3的itchat包简单的分析了一下自己的微信好友。

542
来自专栏玉树芝兰

如何用VOSviewer分析CNKI关键词共现?

用VOSviewer尝试CNKI中文文献关键词共现(keyword co-occurence)分析时,你可能会踩到一个大坑。本文帮助你绕开这个坑,或是从坑里爬出...

512
来自专栏13blog.site

微信公众号问题:{"errcode":40125,"errmsg":"invalid appsecret, view more at http:\/\/t.cn\/LOEdzVq, hints: [

在调试微信公众号授权登录时遇到了这个错误,着实是心烦了半天,公众号相关开发以前是经常做的,很久没有接触了,而且遇到了这么个以前没遇到的问题。 {"errcode...

3886
来自专栏PPV课数据科学社区

大规模爬虫流程总结

爬虫是一个比较容易上手的技术,也许花5分钟看一篇文档就能爬取单个网页上的数据。但对于大规模爬虫,完全就是另一回事,并不是1*n这么简单,还会衍生出许多别的问题。...

27211
来自专栏影子

一张图解析 编译器编译流程

35415
来自专栏华章科技

7大笔记应用,让你的代码效率翻7倍

但是大多数笔记应用的设计并不是以程序员作为目标受众,这些程序可能会让使用者用起来很难受,甚至完全放弃这些工具。这就是为什么我们为你找来了这些最好的笔记工具。快来...

712

扫描关注云+社区