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

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

第一讲：文献选择与解读

文献；CARM1 Methylates Chromatin Remodeling Factor BAF155 to Enhance Tumor Progression and Metastasis

我很早以前想自学CHIP-seq的时候关注过这篇文章，那时候懂得还不多，甚至都没有仔细看这篇文章就随便下载了数据进行分析，也只是跑一些软件而已。这次仔细阅读这篇文章才发现里面门道很多，尤其是ChIP-seq的实验基础和表观遗传学的生物学基础知识。

作者首先实验证明了用small haripin RNA来knockout CARM1 只能达到90%的敲除效果，有趣的是，对CARM1的功能影响非常小，说明只需要极少量的CARM1就可以发挥很好的作用，因此作者通过zinc finger nuclease这种基因组编辑技术设计了100%敲除CARM1的实验材料。(当然，现在有更好的基因编辑技术啦)

这样就能比较CARM1有无时各种蛋白被催化状态了，其中SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物的一个亚基 BAF155，非常明显的只有在CARM1这个基因完好无损的细胞系里面才能被正常的甲基化。作者证明了BAF155是CARM1这个基因非常好的一个底物， 而且通过巧妙的实验设计，证明了BAF155这个蛋白的第1064位氨基酸(R) 是 CARM1的作用位点。

因为早就有各种文献说明了SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物在癌症的重要作用， 所以作者也很自然想探究BAF155在癌症的功能详情，这里作者选择的是ChIP-seq技术。BAF155是作为SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物的一个组分，必然neng 直接或者间接的结合DNA咯。而ChIP-seq技术最适合来探究能直接或者间接结合DNA的蛋白的功能，所以作者构造了一种细胞系（MCF7），它的BAF155蛋白的第1064位氨基酸(R) 突变而无法被CARM1这个基因催化而甲基化，然后比较突变的细胞系和野生型细胞系的BAF155的两个ChIP-seq结果，这样就可以研究BAF155是否必须要被CARM1这个基因催化而甲基化后才能行使生物学功能。

作者用me-BAF155特异性抗体+western bloting 证明了正常的野生型MCF7细胞系里面有~74%的BAF155被甲基化。

有一个细胞系SKOV3，可以正常表达除了BAF155之外的其余14种SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物，而不管是把突变的细胞系和野生型细胞系的BAF155混在里面都可以促进染色质重构复合物的组装，所以甲基化与否并不影响这个染色质重构复合物的组装，重点应该研究的是甲基化会影响BAF155在基因组其它地方结合。

结果显示，突变的细胞系和野生型细胞系种BAF155在基因组结合位置(peaks)还是有较大的overlap的，重点是看它们的peaks在各种基因组区域(基因上下游，5,3端UTR，启动子，内含子，外显子，基因间区域，microRNA区域)分布情况的差别，还有它们距离转录起始位点的距离的分布区别，还有它们注释到的基因区别，已经基因富集到什么通路等等。

虽然作者在人的细胞系(MCF7)上面做ChIP-seq，但是在老鼠细胞系(MDA-MB-231)做了mRNA芯片数据分析，BAF155这个蛋白的第1064位氨基酸(R) 突变细胞系和野生型细胞系，用的是Affymetrix HG U133 Plus 2.0这个常用平台。

which was hybridized to Affymetrix HG U133 Plus 2.0 microarrays containing 54,675 probesets for >47,000 transcripts and variants, including 38,500 human genes.

To identify genes differentially expressed between MDA-MB-231-BAF155WT and MDA-MB-231-BAF155R1064K

表达矩阵下载地址：

我简单摘抄作者ChIP-seq数据的生物信息学分析结果

• All samples were mapped from fastq files using BOWTIE [-m 1 -- best] to mm9 [UCSCmouse genome build 9]

• Sequences were mapped to the human genome (hg19) using BOWTIE (--best –m 1) to yield unique alignments

• Peaks were called by using HOMER [] and QuEST [].

用到的软件有

• *QuEST 2.4 *(Valouev et al., 2008) was run using the recommend settings for transcription factor (TF) like binding with the following exceptions:

  kdebandwith=30, regionsize=600, ChIP threshold=35, enrichment fold=3, rescue fold=3.

• *HOMER *(Heinz et al., 2010) analysis was run using the default settings for peak finding.

  False Discovery Rate (FDR) cut off was *0.001 (0.1%) for all peaks. *

  The tag density for each factor was normalized to 1x107 tags and displayed using the UCSC genome browser.

• Motif analysis (de novo and known), was performed using the* HOMER software and Genomatix. *

• *Peak overlaps *were processed with HOMER and Galaxy (Giardine et al., 2005).

• *Peak comparisons *between replicates were processed with EdgeR statistical package in R

以上就是我们接下来需要学习的流程化分析步骤，下面我给一个主要流程的截图，但主要是实验是如何设计

这里有一个文章发表了关于CHIP-seq的流程的：

同时我还推荐大家看几篇相关文献

• Genome-wide maps of chromatin state in pluripotent and lineage-committed cells. 

• Mapping and analysis of chromatin state dynamics in nine human cell types(GSE26386): 

• Promiscuous RNA binding by Polycomb Repressive Complex 2 

---

第二讲：资料收集

CHIP-seq的确是非常完善的NGS流程，各种资料层出不穷。

大家首先可以看下面几个完整流程的PPT来对CHIP-seq流程有个大致的印象，我对前面提到的文献数据处理的几个要点，就跟下面这个图片类似。

• QuEST is a statistical software for analysis of ChIP-Seq data with data and analysis results visualization through UCSC Genome Browser. ~valouev/QuEST/QuEST.html

• peak calling 阈值的选择： 

• MeDIP-seq and histone modification ChIP-seq analysis 

• 2011-review-CHIP-seq-high-quaility-data: 

• 不同处理条件的CHIP-seq的差异peaks分析： 

• 一个实际的CHIP-seq数据分析例子： ~mthomas/other/chip-seq-training/

然后下面的各种资料，是针对CHIP-seq流程的各个环境的，还有一些是针对于表观遗传学知识

• ppt : 

• best practise: 

• pipeline : 

•  ## samtools view -b -F 1548 -q 30 chipSampleRep1.bam

• pipeline : 

• pipeline : 

• Hands-on introduction to ChIP-seq analysis - VIB Training ~mthomas/other/chip-seq-training/

• video(A Step-by-Step Guide to ChIP-Seq Data Analysis Webinar) : 

• Using ChIP-Seq to identify and/or quantify bound regions (peaks)

• 公开课： 

• EBI的教程：

• 台湾教程： 徐唯哲 Paul Wei-Che HSU

• peak finder软件大全： 

• paper： Large-Scale Quality Analysis of Published ChIP-seq Data 

• paper： Chip-seq data analysis: from quality check to motif discovery and more 

• Workshop hands on session(RNA-Seq / ChIP-Seq ) : 

• paper supplement : 

•  "Analyzing ChIP-seq data: preprocessing, normalization, differential identification, and binding pattern characterization."

•  "Normalization, bias correction, and peak calling for ChIP-seq." (stat heavy)

•  "Practical guidelines for the comprehensive analysis of ChIP-seq data."

•  "Identifying and mitigating bias in next-generation sequencing methods for chromatin biology."

•  "Impact of sequencing depth in ChIP-seq experiments."

• figures: 

可视化工具

bioconductor系列工具和教程 :

• ~tgirke/HTMLPresentations/Manuals/WorkshopDec610_2012/Rchipseq/Rchipseq.pdf

公司教程

---

第三讲：公共数据下载

这一步跟自学其它高通量测序数据处理一样，就是仔细研读paper，在里面找到作者把原始测序数据放在了哪个公共数据库里面，一般是NCBI的GEO，SRA，本文也不例外，然后解析样本数，找到下载链接规律。

1. ## step1 : download raw data

2. > cd ~

3. > mkdir CHIPseq_test && cd CHIPseq_test

4. > mkdir rawData && cd rawData

5. > ## batch download the raw data by shell script :

6. > for ((i=593;i<601;i++)) ;do wget [ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033492/SRR1042](ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033492/SRR1042)$i/SRR1042$i.sra;done

很容易就下载了8个测序文件，每个样本的数据大小，测序量如下

1. > 621M Jun 27 14:03 SRR1042593.sra (16.9M reads)

2. > 2.2G Jun 27 15:58 SRR1042594.sra (60.6M reads)

3. > 541M Jun 27 16:26 SRR1042595.sra (14.6M reads)

4. > 2.4G Jun 27 18:24 SRR1042596.sra (65.9M reads)

5. > 814M Jun 27 18:59 SRR1042597.sra (22.2M reads)

6. > 2.1G Jun 27 20:30 SRR1042598.sra (58.1M reads)

7. > 883M Jun 27 21:08 SRR1042599.sra (24.0M reads)

8. > 2.8G Jun 28 11:53 SRR1042600.sra (76.4M reads)

虽然下载的SRA格式数据也是一个很流行的标准，但它只是数据压缩的标准，几乎没有软件能直接跟SRA的格式的测序数据来进行分析，我们需要转成fastq格式，代码如下：

1. > ## step2 :  change sra data to fastq files.

2. > ## cell line: MCF7 //  Illumina HiSeq 2000 //  50bp // Single ends // phred+33

3. > ## [http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52964](http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52964)

4. > ## [ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033492](ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033492)

5. > ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump $id;done

6. > rm *sra

解压的详情如下，可以看到SRA格式有6~9倍的压缩了，比zip格式压缩的2~3倍高多了

1. ##  621M --> 3.9G

2. ##  2.2G --> 14G

3. ##  541M --> 3.3G

4. ##  2.4G --> 15G

---

第四讲：必要软件安装及结果下载

博文的顺序有点乱，因为怕读到前面的公共测序数据下载这篇文章的朋友搞不清楚，我如何调用各种软件的，所以我这里强势插入一篇博客来描述这件事，当然也只是略过，我所有的软件理论上都是安装在我的home目录下的biosoft文件夹，所以你看到我一般安装程序都是:

1. cd ~/biosoft

2. mkdir macs2 && cd macs2 ##指定的软件安装在指定文件夹里面

这只是我个人的安装习惯，因为我不是root，所以不能在linux系统下做太多事，我这里贴出我所有的软件安装代码：

1. ## pre-step: download sratoolkit /fastx_toolkit_0.0.13/fastqc/bowtie2/bwa/MACS2/HOMER/QuEST/mm9/hg19/bedtools

2. ## http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software

3. ## http://www.ncbi.nlm.nih.gov/books/NBK158900/

4. ## Download and install sratoolkit

5. cd ~/biosoft

6. mkdir sratoolkit && cd sratoolkit

7. wget http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.6.3/sratoolkit.2.6.3-centos_linux64.tar.gz

8. ##

9. ## Length: 63453761 (61M) [application/x-gzip]

10. ## Saving to: "sratoolkit.2.6.3-centos_linux64.tar.gz"

11. tar zxvf sratoolkit.2.6.3-centos_linux64.tar.gz

12. ## Download and install bedtools

13. cd ~/biosoft

14. mkdir bedtools && cd bedtools

15. wget https://github.com/arq5x/bedtools2/releases/download/v2.25.0/bedtools-2.25.0.tar.gz

16. ## Length: 19581105 (19M) [application/octet-stream]

17. tar -zxvf bedtools-2.25.0.tar.gz

18. cd bedtools2

19. make

20. ## Download and install PeakRanger

21. cd ~/biosoft

22. mkdir PeakRanger && cd PeakRanger

23. wget https://sourceforge.net/projects/ranger/files/PeakRanger-1.18-Linux-x86_64.zip/

24. ## Length: 1517587 (1.4M) [application/octet-stream]

25. unzip PeakRanger-1.18-Linux-x86_64.zip

26. ~/biosoft/PeakRanger/bin/peakranger -h

27. ## Download and install bowtie

28. cd ~/biosoft

29. mkdir bowtie && cd bowtie

30. wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.9/bowtie2-2.2.9-linux-x86_64.zip/download

31. #Length: 27073243 (26M) [application/octet-stream]

32. #Saving to: "download" ## I made a mistake here for downloading the bowtie2

33. mv download bowtie2-2.2.9-linux-x86_64.zip

34. unzip bowtie2-2.2.9-linux-x86_64.zip

35. mkdir -p ~/biosoft/bowtie/hg19_index

36. cd ~/biosoft/bowtie/hg19_index

37. # download hg19 chromosome fasta files

38. wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz

39. # unzip and concatenate chromosome and contig fasta files

40. tar zvfx chromFa.tar.gz

41. cat *.fa > hg19.fa

42. rm chr*.fa

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

44. ## Download and install BWA

45. cd ~/biosoft

46. mkdir bwa && cd bwa

47. http://sourceforge.net/projects/bio-bwa/files/

48. tar xvfj bwa-0.7.12.tar.bz2 # x extracts, v is verbose (details of what it is doing), f skips prompting for each individual file, and j tells it to unzip .bz2 files

49. cd bwa-0.7.12

50. make

51. export PATH=$PATH:/path/to/bwa-0.7.12 # Add bwa to your PATH by editing ~/.bashrc file (or .bash_profile or .profile file)

52. # /path/to/ is an placeholder. Replace with real path to BWA on your machine

53. source ~/.bashrc

54. # bwa index [-a bwtsw|is] index_prefix reference.fasta

55. bwa index -p hg19bwaidx -a bwtsw ~/biosoft/bowtie/hg19_index/hg19.fa

56. # -p index name (change this to whatever you want)

57. # -a index algorithm (bwtsw for long genomes and is for short genomes)

58. ## Download and install macs2

59. ## // https://pypi.python.org/pypi/MACS2/

60. cd ~/biosoft

61. mkdir macs2 && cd macs2

62. wget ~~~~~~~~~~~~~~~~~~~~~~MACS2-2.1.1.20160309.tar.gz

63. tar zxvf MACS2-2.1.1.20160309.tar.gz

64. cd MACS2-2.1.1.20160309

65. python setup.py install --user

66. #################### The log for installing MACS2:

67. Creating ~/.local/lib/python2.7/site-packages/site.py

68. Processing MACS2-2.1.1.20160309-py2.7-linux-x86_64.egg

69. Copying MACS2-2.1.1.20160309-py2.7-linux-x86_64.egg to ~/.local/lib/python2.7/site-packages

70. Adding MACS2 2.1.1.20160309 to easy-install.pth file

71. Installing macs2 script to ~/.local/bin

72. Finished processing dependencies for MACS2==2.1.1.20160309

73. ############################################################

74. ~/.local/bin/macs2 --help

75. Example for regular peak calling:

76. macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01

77. Example for broad peak calling:

78. macs2 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1

79. ## Download and install homer (Hypergeometric Optimization of Motif EnRichment)

80. ## // http://homer.salk.edu/homer/

81. ## // http://blog.qiubio.com:8080/archives/3024

82. ## pre-install: Ghostscript，seqlogo,blat

83. cd ~/biosoft

84. mkdir homer && cd homer

85. wget http://homer.salk.edu/homer/configureHomer.pl

86. perl configureHomer.pl -install

87. perl configureHomer.pl -install hg19

一般来说，对我这样水平的人来说，软件安装就跟家常便饭一样，没有什么问题了，但如果你是初学者呢，肯定没那么轻松，所以请加强学习，我无法在这里讲解太具体的知识了。

所有软件安装完毕后就可以下载文章对这些ChIP-seq的处理结果了，这个很重要，检验我们是否重复了人家的数据分析过程。

1. ## step3 : download the results from paper

2. ## http://www.bio-info-trainee.com/1571.html

3. mkdir paper_results && cd paper_results

4. wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE52nnn/GSE52964/suppl/GSE52964_RAW.tar

5. tar xvf GSE52964_RAW.tar

6. ls *gz |xargs gunzip

7. ## step4 : run FastQC to check the sequencing quality.

8. ##这里可以看到我们下载的原始数据已经被作者处理好了，去了接头，去了低质量序列

9. ls *.fastq | while read id ; do ~/biosoft/fastqc/FastQC/fastqc $id;done

10. ## Sequence length 51

11. ## %GC 39

12. ## Adapter Content passed

13. The quality of the reads is pretty good, we don't need to do any filter or trim

14. mkdir QC_results

15. mv *zip *html QC_results/

---

编辑校对：思考问题的熊