专栏名称: 23Plus
23Plus是首个专注于表观遗传学领域的网络社区平台，汇聚全球表观遗传领域专家、学者以及医疗实践者，致力于打造兼专业与科普为一体的的表观遗传互动阵地。
目录
相关文章推荐
生信宝典  ·  全国十大院士之乡及其知名院士一览！一起来看看 ... ·  昨天  
生物制品圈  ·  肽类药物：现状、进展与未来展望 ·  2 天前  
生物探索  ·  Reprod Biol ... ·  4 天前  
生物学霸  ·  DeepSeek 联手 ... ·  2 天前  
BioArt  ·  Nature | GZMK通过激活补体系统加剧炎症 ·  2 天前  
51好读  ›  专栏  ›  23Plus

教你学会ChIP-seq分析 | 第六讲

23Plus  · 公众号  · 生物  · 2017-07-14 07:00

正文

写在前面

本次系列文章为大家带来的是生信菜鸟图案的经典文章合辑： 《教你学会ChIP-seq分析》 ， 共九讲内容 。 带领你从相关文献解读、资料收集和公共数据下载开始，通过软件安装、数据比对、寻找并注释peak、寻找motif等ChIP-seq分析主要步骤入手学习，最后还会介绍相关可视化工具。


第六讲：寻找peaks

ChIP-seq测序的本质还是目标片段捕获测序，但跟WES不同的是， 你选择的IP不同，细胞或者机体状态不同，捕获到的序列差异很大。 而我们研究的重点就是捕获到的差异。我们对ChIP-seq测序数据寻找peaks的本质就是得到所有测序数据比对在全基因组之后在基因组上测序深度里面寻找比较突出的部分。比如对WES数据来说，各个外显子，或者外显子的5端到3端，理论上测序深度应该是一致的，都是50X~200X，画一个测序深度曲线，应该是近似于一条直线。但对我们的ChIP-seq测序数据来说， 在所捕获的区域上面，理论上测序深度是绝对不一样的，应该是近似于一个山峰。


WGS,WES,RNA-seq组与ChIP-seq之间的异同


而那些覆盖度高的地方，山顶，就是我们的IP所结合的热点，也就是我们想要找的peaks，在IGV里面看到大致是下面这样：


可以看到测序的reads是分布不均匀的，我们通常说的 ChIP-seq测序的IP，可以是各个组蛋白上各修饰位点对应的抗体，或者是各种转录因子的抗体等等。


如何定义热点呢？通俗地讲，热点是这样一些位置，这些位置多次被测得的read所覆盖（我们测的是一个细胞群体，read出现次数多，说明该位置被TF结合的几率大）。那么，read数达到多少才叫多？这就要用到统计检验来分析。假设TF在基因组上的分布是没有任何规律的，那么，测序得到的read在基因组上的分布也必然是随机的，某个碱基上覆盖的read的数目应该服从二项分布。


具体统计学原理可以看这篇博客文章： http://www.plob.org/2014/05/08/7227.html


为了达到作者文献里面的结果，我换了 8个软件 ： MACS2/HOMER/SICERpy/PePr/SWEMBL/SISSRs/BayesPeak/PeakRanger


这里就不一一介绍peaks caller软件的安装以及使用了，因为MACS2是最常用的，所以简单贴一下我关于MACS2的学习代码：

  1. ## step6 : peak calling

  2. ### step6.1: with MACS2

  3. ## 我先看了看说明书：

  5. macs2 callpeak -t TF_1.bam -c Input .bam -n mypeaks

  6. We used the following options:

  7. - t: This is the only required parameter for MACS, refers to the name of the file with the ChIP -seq data

  8. - c: The control or mock data file

  9. - n: The name string of the experiment

  10. MAC2 creates 4 files (mypeaks peaks.narrowPeak, mypeaks summits.bed, mypeaks peaks.xls and mypeaks model.r)

  12. # MACS首先的工作是要确定一个模型，这个模型最关键的参数就是峰宽d。这个d就是bw(band width)，而它的一半就是shiftsize。

  13. ### 然后根据文章确定了下载的测序数据的分类

  15. GSM1278641 Xu_MUT_rep1_BAF155_MUT SRR1042593

  16. GSM1278642 Xu_MUT_rep1_Input SRR1042594

  17. GSM1278643 Xu_MUT_rep2_BAF155_MUT SRR1042595

  18. GSM1278644 Xu_MUT_rep2_Input SRR1042596

  19. GSM1278645 Xu_WT_rep1_BAF155 SRR1042597

  20. GSM1278646 Xu_WT_rep1_Input SRR1042598

  21. GSM1278647 Xu_WT_rep2_BAF155 SRR1042599

  22. GSM1278648 Xu_WT_rep2_Input SRR1042600

  24. ## 这里有个很奇怪的问题，input的测序数据居然比IP的测序数据多？？？

  26. 848M Jun 28 14 : 31 SRR1042593.bam

  27. 2.7G Jun 28 14 : 52 SRR1042594.bam

  28. 716M Jun 28 14 : 58 SRR1042595.bam

  29. 2.9G Jun 28 15 : 20 SRR1042596.bam

  30. 1.1G Jun 28 15 : 28 SRR1042597.bam

  31. 2.6G Jun 28 15 : 48 SRR1042598.bam

  32. 1.2G Jun 28 15 : 58 SRR1042599.bam

  33. 3.5G Jun 28 16 : 26 SRR1042600.bam

  35. ## 我没有想明白为什么

  37. ## http://www2.uef.fi/documents/1698400/2466431/Macs2/f4d12870-34f9-43ef-bf0d-f5d087267602

  38. ## http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3120977/我首先用的是下面这些代码

  39. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042594.bam -t SRR1042593.bam -f BAM -B -g hs -n Xu_MUT_rep1 2 > Xu_MUT_rep1 .masc2.log &

  40. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042596.bam -t SRR1042595.bam -f BAM -B -g hs -n Xu_MUT_rep2 2 > Xu_MUT_rep2 .masc2.log &

  41. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042598.bam -t SRR1042597.bam -f BAM -B -g hs -n Xu_WT_rep1 2 > Xu_WT_rep1 .masc2.log &

  42. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042600.bam -t SRR1042599.bam -f BAM -B -g hs -n Xu_WT_rep2 2 > Xu_WT_rep2 .masc2.log &

  44. 得到的 peaks少的可怜，我第一次检查，以为是因为自己没有sort 比对的bam文件导致

  46. ## forget to sort the bam files:

  47. ## 首先把bam文件sort好，构建了inde，然后继续运行！

  49. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042594.sorted.bam -t SRR1042593.sorted.bam -f BAM -B -g hs -n Xu_MUT_rep1 2 > Xu_MUT_rep1 .masc2.log &

  50. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042596.sorted.bam -t SRR1042595.sorted.bam -f BAM -B -g hs -n Xu_MUT_rep2 2 > Xu_MUT_rep2 .masc2.log &

  51. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042598.sorted.bam -t SRR1042597.sorted.bam -f BAM -B -g hs -n Xu_WT_rep1 2 > Xu_WT_rep1 .masc2.log &

  52. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042600.sorted.bam -t SRR1042599.sorted.bam -f BAM -B -g hs -n Xu_WT_rep2 2 > Xu_WT_rep2 .masc2.log &

  54. ##此时得到peaks跟上面为sort的bam文件得到的peaks一模一样，看来不是这个原因

  55. ##然后我怀疑是不是作者上传数据的时候把input和IP标记反了，所以我认为的调整过来

  57. ## Then change the control and treatment

  58. nohup time ~ /.local/ bin/macs2 callpeak -t SRR1042594.sorted.bam -c SRR1042593.sorted.bam -f BAM -B -g hs -n Xu_MUT_rep1 2 > Xu_MUT_rep1 .masc2.log &

  59. nohup time ~ /.local/ bin/macs2 callpeak -t SRR1042596.sorted.bam -c SRR1042595.sorted.bam -f BAM -B -g hs -n Xu_MUT_rep2 2 > Xu_MUT_rep2 .masc2.log &

  60. nohup time ~ /.local/ bin/macs2 callpeak -t SRR1042598.sorted.bam -c SRR1042597.sorted.bam -f BAM -B -g hs -n Xu_WT_rep1 2 > Xu_WT_rep1 .masc2.log &

  61. nohup time ~ /.local/ bin/macs2 callpeak -t SRR1042600.sorted.bam -c SRR1042599.sorted.bam -f BAM -B -g hs -n Xu_WT_rep2 2 > Xu_WT_rep2 .masc2.log &

  63. ##结果，压根就没有peaks了！！！！看了作者并没有搞错

  64. ##接下来我怀疑是自己用samtools view -bhS -q 30   处理了sam文件，这个标准太严格了！！

  65. ##

  66. ## then just use the sam files.

  68. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042594.sam -t SRR1042593.sam -f SAM -B -g hs -n Xu_MUT_rep1 2 > Xu_MUT_rep1 .masc2.log &

  69. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042596.sam -t SRR1042595.sam -f SAM -B -g hs -n Xu_MUT_rep2 2 > Xu_MUT_rep2 .masc2.log &

  70. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042598.sam -t SRR1042597.sam -f SAM -B -g hs -n Xu_WT_rep1 2 > Xu_WT_rep1 .masc2.log &

  71. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042600.sam -t SRR1042599.sam -f SAM -B -g hs -n Xu_WT_rep2 2 > Xu_WT_rep2 .masc2.log &

  73. ## 也没有多几个peaks，最后我只能想到是我的p值太严格了

  74. ## then chang the criteria for p values :

  76. https : //github.com/taoliu/MACS/

  78. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042594.sam -t SRR1042593.sam -f SAM -p 0.01 -g hs -n Xu_MUT_rep1 2 > Xu_MUT_rep1 .masc2.log &

  79. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042596.sam -t SRR1042595.sam -f SAM -p 0.01 -g hs -n Xu_MUT_rep2 2 > Xu_MUT_rep2 .masc2.log &

  80. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042598.sam -t SRR1042597.sam -f SAM -p 0.01 -g hs -n Xu_WT_rep1 2 > Xu_WT_rep1 .masc2.log &

  81. nohup time ~ /.local/ bin/macs2 callpeak -c SRR1042600.sam -t SRR1042599.sam -f SAM -p 0.01 -g hs -n Xu_WT_rep2 2 > Xu_WT_rep2 .masc2.log &

  83. ##我大大减小了P值的标准，结果是输出一大堆的peaks

  85. 18919 Xu_MUT_rep1_peaks .xls

  86. 36277 Xu_MUT_rep2_peaks .xls

  87. 32494 Xu_WT_rep1_peaks .xls

  88. 56080 Xu_WT_rep2_peaks .xls

  90. 问题是这些 peaks根本就都是假阳性！！！

  91. 我手动的 check了几个之前严格过滤条件下的peaks，的确可以看到测序深度是两个山峰形状的曲线

  93. ## check some peaks 手动的 ## chr1 121484235 121485608

  94. ## masc results :

  95. samtools depth -r chr10: 42385331







请到「今天看啥」查看全文


推荐文章
生信宝典  ·  全国十大院士之乡及其知名院士一览！一起来看看你的家乡上榜了吗？
昨天
生物制品圈  ·  肽类药物：现状、进展与未来展望
2 天前
生物探索  ·  Reprod Biol Endocrinol | ViLR技术：无家族史染色体异常检测实现精准胚胎筛选
4 天前
生物学霸  ·  DeepSeek 联手 Kimi，助力研究生轻松搞定 PPT
2 天前
BioArt  ·  Nature | GZMK通过激活补体系统加剧炎症
2 天前
APPSO  ·  买 5000 送 2000？苹果降价大促 1 月 6 日开抢
8 年前
读书小分队  ·  我哪里都可以湿
8 年前
扬子晚报  ·  请放过我们的五星红旗！戛纳红毯上的人间幺蛾，祖国人民快被你气哭了
7 年前
书法在线  ·  喝水的最佳时间：多喝不如会喝
7 年前
华泰睿思  ·  【有色李斌团队】盛和资源(600392,买入): 稀土系列报告一：稀土画龙，钛锆点睛
7 年前
Sov5搜索   ·   小百科   ·   今天看啥   ·   移动版
51好读 - 好文章就要读起来!