博文
神技能！转录调控必读
||
本文转载自嘉因微信公众号，已获得授权。查看最新文章，敬请关注嘉因，微信ID：rainbow-genome
作者：小丫 来源： 嘉因
转录调控的经典问题：哪个蛋白质调控我感兴趣的基因？为了帮助小伙伴儿解决这个问题，小丫写了系列技术贴：
motif系列Plan B：去哪找motif、找启动子区的motif、互补链上的motif有意义吗？怎样展示结果？
Plan C：不需要那么多ChIP-seq
要实现PlanA，常常遇到sample太多的情况，就需要写代码批量处理。物以类聚，我把这部分代码发布在生信技能树上，是因为生信技能树积累了海量的生信原创学习资料和原创代码资源。就像一下子跳进了水果堆，左手🍌，右手🍒，嘴边儿是🍓，眼前是🍍。回复“大树”，获得海量生信学习资源。
本文转载自生信技能树公众号，已获授权
果子导读：
我的导师曾经跟我讲过，10年前，CELL杂志每期一半以上都是在做转录调控。10年后，我们发现，在很多杂志，这个现象依然存在。
如果已知转录因子，找他的靶基因，用ChIP-seq就可以搞定。
但是！如果反过来，团队已经确定所研究的基因，那么找出能够调控他的分子，确实是个难题。
所幸！这篇来自嘉因小丫的帖子把这个事情一次性解决。
从此，国自然申请和课题设计，如虎添翼，一日千里。
以下是正文：
ChIP系列带你实现Plan A。
原理
在线快速查看结果
局限性
速查表
从哪里下载数据
怎样批量处理数据
怎样展示结果
用上篇《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文件。
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-seq、ATAC-seq研究感兴趣的基因？想整合ChIP-seq、ATAC-seq、eCLIP/RIP-seq、RNA-seq数据寻找线索？找嘉因生物吧！从实验、测序，到多种数据整合分析，为您一站式解决。（点击文中蓝字了解详情）
嘉因公众号定位：客户共性问题解答，生信学习资源导航，高通量实验导购 | 为您提供高通量实验-测序-分析-验证一站式解决方案
电话：021-61539657
Email：marketing@rainbow-genome.com
地址：上海市杨浦区国伟路135号10号楼305室
https://blog.sciencenet.cn/blog-3372875-1106115.html
上一篇：画泡泡图，发CNS
下一篇：Y叔又炫技了！有啥用？还真能发CNS！
