ChengyangWang的个人博客分享 http://blog.sciencenet.cn/u/ChengyangWang

博文

神技能!转录调控必读

已有 5245 次阅读 2018-3-28 15:03 |个人分类:转录调控|系统分类:科普集锦

 本文转载自嘉因微信公众号,已获得授权。查看最新文章,敬请关注嘉因,微信ID:rainbow-genome

作者:小丫  来源: 嘉因

转录调控的经典问题:哪个蛋白质调控我感兴趣的基因?为了帮助小伙伴儿解决这个问题,小丫写了系列技术贴:

 

要实现PlanA,常常遇到sample太多的情况,就需要写代码批量处理。物以类聚,我把这部分代码发布在生信技能树上,是因为生信技能树积累了海量的生信原创学习资料原创代码资源。就像一下子跳进了水果堆,左手🍌,右手🍒,嘴边儿是🍓,眼前是🍍。回复“大树”,获得海量生信学习资源

 


 

本文转载自生信技能树公众号,已获授权

神技能!批量解决哪个转录因子调控你的基因

 

果子导读:

我的导师曾经跟我讲过,10年前,CELL杂志每期一半以上都是在做转录调控。10年后,我们发现,在很多杂志,这个现象依然存在。

 

如果已知转录因子,找他的靶基因,用ChIP-seq就可以搞定。

但是!如果反过来,团队已经确定所研究的基因,那么找出能够调控他的分子,确实是个难题。

所幸!这篇来自嘉因小丫的帖子把这个事情一次性解决。

从此,国自然申请和课题设计,如虎添翼,一千里。

以下是正文:

 

ChIP系列带你实现Plan A

  1. 原理

  2. 在线快速查看结果

  3. 局限性

  4. 速查表

  5. 从哪里下载数据

  6. 怎样批量处理数据

  7. 怎样展示结果

 

用上篇《Plan A详细步骤1234》找到转录因子的小伙伴可以跳过本篇,直接看7,下篇讲7。本文讲56,适用于从几十上百套ChIP-seq中找上游调控因子的情况。如果在嘉因公众号讲这篇,需要铺垫太多基础知识,读者也未必愿意看。思来想去,还是放到生信技能树发布吧。

 


 

书接上文《Plan A详细步骤1234》,如果您关注的细胞类型有几十甚至几百套ChIP-seq,用肉眼挨个看哪个track有peak,就要疯掉了。这时就需要我们懂生信的出手了,批量下载,批量处理。跟2对应,分别讲Cistrome和ENCODE这两个数据库的数据下载和处理。

 

本文要点

  • Cistrome Data Browser,下载peak.bed,用它提供的GEO ID下载原始数据,跑出bigwig。

  • ENCODE,下载bam、bigwig,用bam文件call peak,跑出peak.bed。

  • 用peak.bed找出基因附近有peak的ChIP-seq,后面需要用bigwig画图。

 

5. 从哪里下载数据

Cistrome

Cistrome Data Browser提供Batch download,能批量下载sample信息、peak.bed文件。

http://cistrome.org/db/#/

 

sample信息包括转录因子的名字、细胞类型、器官、细胞系、GEO ID。

peak.bed是ChIP-seq的峰所在的位置,它标志着转录因子结合在这个位置,包括直接结合和间接结合,间接结合例如蛋白A跟蛋白B形成复合物再结合DNA。

 

目前Cistrome Data Browser不提供bigwiggle文件的下载,如果想要本地画图,用batch download提供的sample信息里的GEO ID下载原始数据,参照它的流程自己处理,获得bigwiggle。

 

ENCODE

https://www.encodeproject.org,提供fastq、bam、bigwig、peak.bed文件的下载。点击Download:

 

 

获得这些文件的地址,然后用红色那行命令批量下载需要的文件。

 


ENCODE的peak不理想,所以只在这里下载bam文件,然后用Cistrome Data Browser的参数call peak。

 

6. 数据批量处理

先去UCSC找您感兴趣的基因在基因组上的位置,http://genome.ucsc.edu/,例如TP53,TP53 (uc060aur.1) at chr17:7668402-7687538,整理成这样的bed格式:

chr17    7668402    7687538

中间不是空格,而是用键盘上的tab键输入的制表符,存成文本文件,文件名:TP53.bed


用bedtools的slop算出基因上下游10kb的位置,http://bedtools.readthedocs.io/en/latest/content/tools/slop.html?highlight=slop,用手算也行。

 

用bedtools的intersect找TP53两侧10kb范围跟peak.bed的交集http://bedtools.readthedocs.io/en/latest/content/tools/intersect.html,判断这套ChIP-seq有没有peak落在基因附近。

 

 

如果某个转录因子的ChIP-seq数据在TP53附近有peak,就说明这个转录因子在TP53附近结合,推测它参与了TP53的转录调控。

 

 

bedtools一次只能处理一个ChIP-seq的peak.bed。想批量处理所有的ChIP-seq文件,用shell实现:

for file in ./bed/*.bed;do bedtools intersect -a gene_loci.bed -b $file -c > ./peak/gene_loci_c_${file##*/};done

这样就获得了每个ChIP-seq在基因附近出现peak的数量,0表示没有peak,1或更多表示有1个或多个peak出现在基因附近。先用shell提取peak数量这列:

for file in ./peak/*_c_*.bed;do awk '{print $5}' $file> ./peak/${file##*/}.tmp;done

再把这些文件合并到一起,方便查看,用R:

filename<-dir("path/peak/")filenametmp=read.table(filename[1])colnames(tmp)<-filename[1]tmp_comb=tmpfor (i in 2:length(filename)){  tmp=read.table(filename[i])  colnames(tmp)<-filename[i]  tmp_comb=cbind(tmp_comb,tmp)}iwrite.table(tmp_comb,"combine.tmp",sep="\t",row.names=T,col.names=F,quote=F)

这样就看到哪些ChIP-seq在基因附近有结合信号。

 

然后,就可以用bigwig画自己想看的转录因子了。

 

另外,还可以回到上篇的2,cistrome有选择功能,在有peak的ChIP-seq前面打勾,一起放到WashU Browser查看。

 

7. 怎样展示结果

 

结果图怎样展示?参考这篇介绍的工具《找到了motif,怎样展示结果?》。具体怎么画?且听下回分解。



 



想用ChIP-seqATAC-seq研究感兴趣的基因?想整合ChIP-seq、ATAC-seq、eCLIP/RIP-seq、RNA-seq数据寻找线索?找嘉因生物吧!从实验、测序,到多种数据整合分析,为您一站式解决。(点击文中蓝字了解详情)




嘉因公众号定位:客户共性问题解答,生信学习资源导航,高通量实验导购 | 为您提供高通量实验-测序-分析-验证一站式解决方案


 


电话:021-61539657

Email:marketing@rainbow-genome.com
地址:上海市杨浦区国伟路135号10号楼305室





https://wap.sciencenet.cn/blog-3372875-1106115.html

上一篇:画泡泡图,发CNS
下一篇:Y叔又炫技了!有啥用?还真能发CNS!
收藏 IP: 124.77.56.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-12-3 23:07

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部