一篇文章学会ChIP-seq分析（下）

写在前面：《一篇文章学会ChIP-seq分析（上）》《一篇文章学会ChIP-seq分析（下）》为生信菜鸟团博客相关文章合集，共九讲内容。带领你从相关文献解读、资料收集和公共数据下载开始，通过软件安装、数据比对、寻找并注释peak、寻找motif等ChIP-seq分析主要步骤入手学习，最后还会介绍相关可视化工具。

第五讲：测序数据比对

比对就很简单的了，各种mapping工具层出不穷，我们一般常用的就是BWA和bowtie了，我这里就挑选bowtie2吧，反正别人已经做好了各种工具效果差异的比较，我们直接用就好了，代码如下：

  1. ## step5 : alignment to hg19/ using bowtie2 to do alignment

  2. ## ~/biosoft/bowtie/bowtie2-2.2.9/bowtie2-build ~/biosoft/bowtie/hg19_index /hg19.fa ~/biosoft/bowtie/hg19_index/hg19

  3. ## cat >run_bowtie2.sh

  4. ls *.fastq | while read id ;

  5. do

  6. echo $id

  7. #~/biosoft/bowtie/bowtie2-2.2.9/bowtie2 -p 8 -x ~/biosoft/bowtie/hg19_index/hg19 -U $id -S ${id%%.*}.sam 2>${id%%.*}.align.log;

  8. #samtools view -bhS -q 30 ${id%%.*}.sam > ${id%%.*}.bam ## -F 1548 https://broadinstitute.github.io/picard/explain-flags.html

  9. # -F 0x4 remove the reads that didn't match

  10. samtools sort ${id%%.*}.bam ${id%%.*}.sort ## prefix for the output

  11. # samtools view -bhS a.sam | samtools sort -o - ./ > a.bam

  12. samtools index ${id%%.*}.sorted.bam

  13. done

这个索引~/biosoft/bowtie/hg19_index/hg19需要自己提取建立好，见前文

初步比对的sam文件到底该如何过滤，我查了很多文章都没有给出个子丑寅卯，各执一词，我也没办法给大家一个标准，反正我测试了好几种，看起来call peaks的差异不大，就是得不到文章给出的那些结果！！

一般来说，初步比对的sam文件只能选取unique mapping的结果，所以我用了#samtools view -bhS -q 30，但是结果并没什么改变，有人说是peak caller这些工具本身就会做这件事，所以取决于你下游分析所选择的工具。

给大家看比对的日志吧：

  1. SRR1042593.fastq

  2. 16902907 reads; of these:

  3. 16902907 (100.00%) were unpaired; of these:

  4. 667998 (3.95%) aligned 0 times

  5. 12467095 (73.76%) aligned exactly 1 time

  6. 3767814 (22.29%) aligned >1 times

  7. 96.05% overall alignment rate

  8. ......

  9. SRR1042598.fastq

  10. 58068816 reads; of these:

  11. 58068816 (100.00%) were unpaired; of these:

  12. 8433671 (14.52%) aligned 0 times

  13. 37527468 (64.63%) aligned exactly 1 time

  14. 12107677 (20.85%) aligned >1 times

  15. 85.48% overall alignment rate

  16. [samopen] SAM header is present: 93 sequences.

  17. SRR1042599.fastq

  18. 24019489 reads; of these:

  19. 24019489 (100.00%) were unpaired; of these:

  20. 1411095 (5.87%) aligned 0 times

  21. 17528479 (72.98%) aligned exactly 1 time

  22. 5079915 (21.15%) aligned >1 times

  23. 94.13% overall alignment rate

  24. [samopen] SAM header is present: 93 sequences.

  25. SRR1042600.fastq

  26. 76361026 reads; of these:

  27. 76361026 (100.00%) were unpaired; of these:

  28. 8442054 (11.06%) aligned 0 times

  29. 50918615 (66.68%) aligned exactly 1 time

  30. 17000357 (22.26%) aligned >1 times

  31. 88.94% overall alignment rate

  32. [samopen] SAM header is present: 93 sequences.

可以看到比对非常成功。

我这里就不用表格的形式来展现了，毕竟我又不是给客户写报告，大家就将就着看吧。

第六讲: 寻找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. ## 我先看了看说明书：

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

  5. We used the following options:

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

  7. -c: The control or mock data file

  8. -n: The name string of the experiment

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

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

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

  12. GSM1278641 Xu_MUT_rep1_BAF155_MUT        SRR1042593

  13. GSM1278642 Xu_MUT_rep1_Input        SRR1042594

  14. GSM1278643 Xu_MUT_rep2_BAF155_MUT        SRR1042595

  15. GSM1278644 Xu_MUT_rep2_Input        SRR1042596

  16. GSM1278645 Xu_WT_rep1_BAF155        SRR1042597

  17. GSM1278646 Xu_WT_rep1_Input        SRR1042598

  18. GSM1278647 Xu_WT_rep2_BAF155        SRR1042599

  19. GSM1278648 Xu_WT_rep2_Input         SRR1042600

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

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

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

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

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

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

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

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

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

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

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

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

  32. 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 &

  33. 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 &

  34. 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 &

  35. 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 &

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

  37. ## forget to sort the bam files:

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

  39. 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 &

  40. 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 &

  41. 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 &

  42. 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 &

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

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

  45. ## Then change the control and treatment

  46. 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 &

  47. 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 &

  48. 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 &

  49. 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 &

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

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

  52. ##

  53. ## then just use the sam files.

  54. 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 &

  55. 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 &

  56. 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 &

  57. 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 &

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

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

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

  61. 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 &

  62. 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 &

  63. 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 &

  64. 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 &

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

  66. 18919 Xu_MUT_rep1_peaks.xls

  67. 36277 Xu_MUT_rep2_peaks.xls

  68. 32494 Xu_WT_rep1_peaks.xls

  69. 56080 Xu_WT_rep2_peaks.xls

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

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

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

  73. ## masc results :

  74. samtools depth -r chr10:42385331-42385599 SRR1042593.sorted.bam

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

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

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

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

  79. paper results:

  80. chr20 45796362 46384917

  81. chr1 121482722 121485861

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

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

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

  85. 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/

第七讲: peaks注释

经过前面的ChIP-seq测序数据处理的常规分析，我们已经成功的把测序仪下机数据变成了BED格式的peaks记录文件，我选取的这篇文章里面做了4次ChIP-seq实验，分别是两个重复的野生型MCF7细胞系的 BAF155 immun =oprecipitates和两个重复的突变型MCF7细胞系的 BAF155 immunoprecipitates，这样通过比较野生型和突变型MCF7细胞系的 BAF155 immunoprecipitates的不同结果就知道该细胞系的BAF155 突变，对它在全基因组的结合功能的影响。

  1. #我这里直接从GEO里面下载了peaks结果，它们详情如下：wc -l *bed

  2. 6768 GSM1278641_Xu_MUT_rep1_BAF155_MUT.peaks.bed

  3. 3660 GSM1278643_Xu_MUT_rep2_BAF155_MUT.peaks.bed

  4. 11022 GSM1278645_Xu_WT_rep1_BAF155.peaks.bed

  5. 5260 GSM1278647_Xu_WT_rep2_BAF155.peaks.bed

  6. 49458 GSM601398_Ini1HeLa-peaks.bed

  7. 24477 GSM601398_Ini1HeLa-peaks-stringent.bed

  8. 12725 GSM601399_Brg1HeLa-peaks.bed

  9. 12316 GSM601399_Brg1HeLa-peaks-stringent.bed

  10. 46412 GSM601400_BAF155HeLa-peaks.bed

  11. 37920 GSM601400_BAF155HeLa-peaks-stringent.bed

  12. 30136 GSM601401_BAF170HeLa-peaks.bed

  13. 25432 GSM601401_BAF170HeLa-peaks-stringent.bed

每个BED的peaks记录，本质是就3列是需要我们注意的，就是染色体，以及在该染色体上面的起始和终止坐标，如下

  1. #PeakID chr start end strand Normalized Tag Count region size findPeaks Score Clonal Fold Change

  2. chr20 52221388 52856380 chr20-8088 41141 +

  3. chr20 45796362 46384917 chr20-5152 31612 +

  4. chr17 59287502 59741943 chr17-2332 29994 +

  5. chr17 59755459 59989069 chr17-667 19943 +

  6. chr20 52993293 53369574 chr20-7059 12642 +

  7. chr1 121482722 121485861 chr1-995 9070 +

  8. chr20 55675229 55855175 chr20-6524 7592 +

  9. chr3 64531319 64762040 chr3-4022 7213 +

  10. chr20 49286444 49384563 chr20-4482 6165 +

我们所谓的peaks注释，就是想看看该peaks在基因组的哪一个区段，看看它们在各种基因组区域(基因上下游，5,3端UTR，启动子，内含子，外显子，基因间区域，microRNA区域)分布情况，但是一般的peaks都有近万个，所以需要批量注释，如果脚本学的好，自己下载参考基因组的GFF注释文件，完全可以自己写一个，我这里会介绍一个R的bioconductor包ChIPpeakAnno来做CHIP-seq的peaks注释，下面的包自带的示例

  1. library(ChIPpeakAnno)

  2. bed <- system.file("extdata", "MACS_output.bed", package="ChIPpeakAnno")

  3. gr1 <- toGRanges(bed, format="BED", header=FALSE)

  4. ## one can also try import from rtracklayer

  5. library(rtracklayer)

  6. gr1.import <- import(bed, format="BED")

  7. identical(start(gr1), start(gr1.import))

  8. gr1[1:2]

  9. gr1.import[1:2] #note the name slot is different from gr1

  10. gff <- system.file("extdata", "GFF_peaks.gff", package="ChIPpeakAnno")

  11. gr2 <- toGRanges(gff, format="GFF", header=FALSE, skip=3)

  12. ol <- findOverlapsOfPeaks(gr1, gr2)

  13. makeVennDiagram(ol)

  14. ##还可以用binOverFeature来根据特定的GRanges对象(通常是TSS)来画分布图

  15. ## Distribution of aggregated peak scores or peak numbers around transcript start sites.

可以看到这个包使用起来非常简单，只需要把我们做好的peaks文件(GSM1278641XuMUTrep1BAF155_MUT.peaks.bed等等)用 toGRanges或者 import读进去，成一个GRanges对象即可，上面的代码是比较两个peaks文件的overlap。然后还可以根据R很多包都自带的数据来注释基因组特征

  1. data(TSS.human.GRCh37) ## 主要是借助于这个GRanges对象来做注释，也可以用getAnnotation来获取其它GRanges对象来做注释

  2. ## featureType ： TSS, miRNA, Exon, 5'UTR, 3'UTR, transcript or Exon plus UTR

  3. peaks=MUT_rep1_peaks

  4. macs.anno <- annotatePeakInBatch(peaks, AnnotationData=TSS.human.GRCh37,

  5. output="overlapping", maxgap=5000L)

  6. ## 得到的macs.anno对象就是已经注释好了的，每个peaks是否在基因上，或者距离基因多远，都是写的清清楚楚

  7. if(require(TxDb.Hsapiens.UCSC.hg19.knownGene)){

  8. aCR<-assignChromosomeRegion(peaks, nucleotideLevel=FALSE,

  9. precedence=c("Promoters", "immediateDownstream",

  10. "fiveUTRs", "threeUTRs",

  11. "Exons", "Introns"),

  12. TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene)

  13. barplot(aCR$percentage)

  14. }

得到的条形图如下，虽然很丑，但这就是peaks注释的精髓，搞清楚每个peaks在基因组的位置特征：

同理，对每个peaks文件，都可以做类似的分析！

但是对多个peaks文件，比如本文中的，想比较野生型和突变型MCF7细胞系的 BAF155 immunoprecipitates的结果的不同，就需要做peaks之间的差异分析，已经后续的差异基因注释啦

当然，值得一提的是peaks注释我更喜欢网页版工具，反正peaks文件非常小，直接上传到别人做好的web tools，就可立即出一大堆可视化图表分析结果啦，大家可以去试试看：

http://chipseek.cgu.edu.tw/

http://bejerano.stanford.edu/great/public/html/

http://liulab.dfci.harvard.edu/CEAS/

虽然我花费了大部分篇幅来描述ChIPpeakAnno这个包的用法，但是真正的重点是你得明白peaks记录了什么，要注释什么，以及把这3个网页工具的可视化图表分析结果全部看懂，这网页版工具才是重点！

第八讲：寻找motif

motif是比较有特征的短序列，会多次出现的，一般认为它的生物学意义重大，做完CHIP-seq分析之后，一般都会寻找motif 。

查找有两种，一种是de novo的，要求的输入文件的fasta序列，一般是根据peak的区域的坐标提取好序列；另一种是依赖于数据库的搜寻匹配，很多课题组会将现有的ChIP-seq数据进行整合，提供更全面，更准确的motif数据库。

motif的定义如下：

motif: recurring pattern. eg, sequence motif, structure motif or network motif

DNA sequence motif: short, recurring patterns in DNA that are presumed to have a biological function.

从上边的定义可以看出，其实motif这个单词就是形容一种反复出现的模式，而序列motif往往是DNA上的反复出现的模式，并被假设拥有生物学功能。而且，经常是一些具有序列特异性的蛋白的结合位点（如，转录因子）或者是涉及到重要生物过程的（如，RNA 起始，RNA 终止， RNA 剪切等等）。

摘抄自：http://blog.163.com/zju_whw/blog/static/225753129201532104815301/

motif最先是通过实验的方法发现的，换句话说，不是说有了ChIP-seq才有了motif分析，起始很早人们就开始研究motif了！例如，'TATAAT’ box在1975年就被pribnow发现了，它与上游的'TTGACA’motif是RNA聚合酶结合位点的特异性序列。而且，当时的人们就知道，不是所有的结合位点都一定完美地与motif匹配，大部分都只匹配了12个碱基中的7-9个。结合位点与motif的匹配程度往往也与蛋白质与DNA的结合强弱有关。

目前被人们识别出来的motif也越来越多，如TRANSFAC和JASPAR数据库都有着大量转录因子的motif。而随着ChIP-seq数据的大量产出，motif的研究会进一步深入，有一些课题组会将现有的ChIP-seq数据进行整合，提供更全面，更准确的motif数据库。

从算法上来讲，这是很复杂的，我就不多说了，我这里主要讲best practice：

一篇文献列出了2014年以前的近乎所有知名的A survey of motif finding Web tools for detecting binding site motifs in ChIP-Seq data 链接见：https://biologydirect.biomedcentral.com/articles/10.1186/1745-6150-9-4

最常用的是 MEME工具套件 ：

http://meme-suite.org/ 输入文件是fasta序列，需要对peaks进行转换，根据bed的基因坐标从基因组里面提取对应的序列咯： http://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html

它里面集成了4个寻找motif 的工具，每个工具都是一篇文章，里面有详细描述具体原理，但是整个网页给人的感觉是too busy，让初学者无从下手。

把自己的fasta序列上传上去即可，还是选取我们本次系列教程的数据

  1. $ ls -lh  *fasta

  2. -rw-r--r-- 1 Jimmy 197121  18M Jul  7 19:40 GSM1278641_Xu_MUT_rep1_BAF155_MUT_sequence.fasta

  3. -rw-r--r-- 1 Jimmy 197121 9.9M Jul  7 19:38 GSM1278643_Xu_MUT_rep2_BAF155_MUT_sequence.fasta

  4. -rw-r--r-- 1 Jimmy 197121  26M Jul  7 19:41 GSM1278645_Xu_WT_rep1_BAF155_sequence.fasta

  5. -rw-r--r-- 1 Jimmy 197121  14M Jul  7 19:41 GSM1278647_Xu_WT_rep2_BAF155_sequence.fasta

然后就可以看到所有结果啦，大家可以试试看。

最后值得一提的是现在流行的R的bioconductor系列包，也可以寻找motif：

一般的R包都可以直接从BED文件里面记录的基因坐标来找motif，有点需要输入fasta序列，就需要自己根据bed的基因坐标从基因组里面提取对应的序列咯：

rGADEM (motif discovery): http://bioconductor.org/packages/devel/bioc/html/rGADEM.html

MotIV (motif validation): http://bioconductor.org/packages/devel/bioc/html/MotIV.html

http://lgsun.grc.nia.nih.gov/CisFinder/

http://bioinfo.cs.technion.ac.il/drim/

http://www.ncbi.nlm.nih.gov/pubmed/20736340

还有一个PICS (ChIP-seq): 虽然不是bioconductor的包 http://www.rglab.org/pics-probabilistic-inference-for-chip-seq/ 貌似国内被墙了，无法打开。

第九讲：ChIP-seq可视化大全

讲到这里，我们的自学ChIP-seq分析系列教程就告一段落了，当然，我会随时查漏补缺，根据读者的反馈来更新着系列教程。

其实可视化这已经是一个比较复杂的方向了，不仅仅是针对于ChIP-seq数据。可视化本身是发文章的先决条件，而让人一目了然图片也说明了数据分析人员对数据本身的理解。我这里就列出一些目录和一些工具和ppt。这个主要靠大家自学，而且我博客空间有限，就不上传一大堆图片了，大家随便找一些经典的paper里面都会有很多可视化分析。

首先强烈推荐两个网页版工具，针对找到的peaks可视化:

http://chipseek.cgu.edu.tw/

http://bejerano.stanford.edu/great/public/html/

然后再推荐一个哈佛刘小乐实验室出品的软件，也是专门为了作图 http://liulab.dfci.harvard.edu/CEAS/usermanual.html

还有一个java工具：也可以可视化CHIP-seq的peaks结果EXPANDER (EXpression Analyzer and DisplayER) is a java-based tool for analysis of gene expression data.http://acgt.cs.tau.ac.il/expander/help/ver7.0Help/html/InputData.htm

如图所示

然后我所了解的图片大概有下面这些，都是有专门的软件，甚至自己写脚本也可以做的

  • peaks长度分布柱状图

  • 每个peak的测序情况可视化(IGV,sushi)

  • 测序reads在全基因组各个染色体的分布(Chromosome ideograms)

  • reads相对基因位置分布统计

  • peaks相对基因位置分布统计

  • reads在基因组位置分布统计（染色体分开作图）

  • peaks在基因组位置分布统计（染色体分开作图）

  • 统计peaks在各种基因组区域(基因上下游，5,3端UTR，启动子，内含子，外显子，基因间区域，- microRNA区域)分布情况，条形图和饼图均可

  • peak与转录起始位点距离的分析（曲线图和热图）

最后总结一下

CHIP-seq pipeline :　http://www.slideshare.net/COST-events/chipseq-data-analysis

大家一定要看这个ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia.http://www.ncbi.nlm.nih.gov/pubmed/22955991

编辑校对：思考问题的熊

赞 (0)

相关推荐