专栏名称: 23Plus
23Plus是首个专注于表观遗传学领域的网络社区平台，汇聚全球表观遗传领域专家、学者以及医疗实践者，致力于打造兼专业与科普为一体的的表观遗传互动阵地。
目录
相关文章推荐
生物学霸  ·  Cancer ... ·  2 天前  
BioArt  ·  第十一届国际三维基因组学研讨会第三轮通知 ·  3 天前  
生物学霸  ·  这也行？IF 高达 98 ... ·  4 天前  
生物学霸  ·  5 分毕业神刊出圈，中科院二区，平均 2 ... ·  6 天前  
BioArt  ·  招聘脑科学、心理学、AI等科研人员：深圳理工 ... ·  6 天前  
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-42385599 SRR1042593.sorted.bam

  96. samtools depth -r chr10:42385331-42385599 SRR1042594.sorted.bam

  97. samtools depth -r chr20:45810382-45810662 SRR1042593.sorted.bam

  98. samtools depth -r chr20:45810382-45810662 SRR1042594.sorted.bam

  100. ##我也check了paper里面得到的peak，但是在我的比对文件里面，肉眼看起来根本不像，所以我很纠结

  102. paper results:

  103. chr20 45796362 46384917

  104. chr1 121482722 121485861

  105. samtools depth -r chr1:121482722-121485861 SRR1042593.sorted.bam

  106. samtools depth -r chr1:121482722-121485861 SRR1042594.sorted.bam

  107. samtools depth -r chr20:45796362-46384917 SRR1042593.sorted.bam

  108. samtools depth -r chr20:45796362-46384917 SRR1042594.sorted.bam


很不幸，最后还是没能达到作者的结果，我没搞清楚是为什么，我还用了BayesPeak/PeakRanger这两个软件，结果也不理想。

  • peak finder软件大全： http://wodaklab.org/nextgen/data/peakfinders.html

  • Peak Calling for ChIP-Seq :　http://epigenie.com/guide-peak-calling-for-chip-seq/


本系列历史文章列表

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

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

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

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

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


本文转载自


“生信技能树”公众号

初与大家分享自己的生信学习笔记及心得体会。促进生信的学习和交流，构建出完整的生信技能树。搭建生信技术人员联盟，从入门到进阶帮助到每一位生信人。最期待看到团队成员的成长，以及论坛稳健发展和各版块完善。带领团队和论坛成员完善生信技能树的同时，自己也收获前所未有的锻炼，希望自己不忘初心。


"生信技能树"论坛

生信技能树创建于2016年8月，是中国第一家专注于生信知识体系完善、促进生信学习交流的论坛。我们通过收集国内外生信学习资源，邀请大神分享的领域专业知识，发布菜鸟的真实学习笔记，搭建生信技术人员联盟，从入门到进阶帮助每一位生信人。

欢迎点击文末“阅读原文”了解“生信技能树”论坛，上面有本文作者jimmy原创的一千多篇教程


23Plus欢迎表观遗传领域的同行们投稿，分享学术成果、学术见解和学术故事。

投稿请联系：

[email protected]

微信添加23Plus小秘书：

plus23_sec

拉您入群参加更深入的讨论。


23Plus: 首个专注于表观遗传学领域的网络社区

微信号：epi23plus

推荐文章
生物学霸  ·  Cancer Cell：华山医院毛颖/花玮/迟喻丹团队肿瘤学领域获重磅进展
2 天前
BioArt  ·  第十一届国际三维基因组学研讨会第三轮通知
3 天前
生物学霸  ·  这也行？IF 高达 98 分的顶刊论文，配图竟都是火柴人，网友：换我也能发
4 天前
生物学霸  ·  5 分毕业神刊出圈，中科院二区，平均 2 个月可接收，录用率高达 80%
6 天前
BioArt  ·  招聘脑科学、心理学、AI等科研人员：深圳理工大学-“智能与行为交叉研究中心”
6 天前
风青杨  ·  一位美国妈妈给孩子定的手机使用家规，值得中国爸妈看看！
7 年前
萌宠萌  ·  喵：铲屎的，别动，这是朕的皇冠！
7 年前
中国好文章  ·  中国的南方人和北方人有什么区别？总算说透了！
7 年前
婚姻家庭那些事儿  ·  为什么在爱情里，女人往往比男人勇敢？
7 年前
星座不求人  ·  我是双鱼我骄傲
7 年前
Sov5搜索   ·   小百科   ·   移动版
51好读 - 好文章就要读起来!