ChIP-seq分析的一般流程方法

现在ChIP-seq的数据基本是最常见的测序数据类型之一，主要有Transcription factor ChIP-seq和Histone ChIP-seq。前者是看转录因子的结合位置，后者是组蛋白修饰发生的位置。下面分享一下一般流程。

ChIP-seq

1.质控 (quality control)

首先要看一下ChIP-seq数据的质量，数据的信号最好比background要强很很多。一般要有control，这样call peaks更准确可信， control主要有Input DNA 和 IgG两种，前一种更常用。

检测质量的一些方式： 1). peaks中reads的数量，如果peaks的reads普遍较少，则质量一般。 2). peaks信号高，背景低。 3). 测序深度深 。 4). Diverse library (与重复duplications有关，如下图)

library

4). 有重复并且与重复之间相似性较高… ……

做质控的软件方法： 1). ChIPQC (T Carroll, Front Genet, 2014.) 2). SPP package - Unix/Linux (PV Karchenko, Nature Biotechnol, 2008.) 3). ENCODE中的标准流程

2.序列比对 (mapping of fastq)

序列比对一般用BWA或者Bowtie2，两者效果差不多。BWA的bwa samse（单端数据）和bwa sampe (双端数据) 跑的速度比较慢，但是效果很不错，用法如下：

bwa index reference.fa   # 建立索引 -p可设置前缀，不设置前缀就是reference.fa。

# 单端数据
bwa aln -t 8 reference.fa  test.fq.gz > test.sai
bwa samse -n 10 reference.fa  test.sai test.fq.gz > test_se.sam 

# 双端数据：
bwa aln reference.fa test_reads1.fq > test1.sai
bwa aln reference.fa test_reads2.fq > test2.sai
bwa sampe reference.fa test1.sai test2.sai test_reads1.fq test_reads2.fq >  test_pe.sam

BWA的mem，速度很快：

bwa mem reference.fa reads.fq > test_se.sam # 单端
bwa mem reference.fa read1.fq read2.fq > test_pe.sam # 双端

bowtie2的用法：

bowtie2-build reference.fa index # 创建index

bowtie2  -p 8 -x index -1 test_read1.fq -2 test_read2.fq -S test_se.sam # 单端比对
bowtie2  -p 8 -x index -1 test_read1.fq -2 test_read2.fq -S test_pe.sam # 双端比对

效果个人感觉差不多。

3.去除重复 (remove duplicates)

由于PCR实验存在不可避免的实验误差，所以会存在重复 (duplicates)。比如两条不同的reads，起止位置完全一致。比如：

example_dup

其中第二条已经被picard标注出来了。被标注的第二列flag会加1024。

去重的软件中samtools rmdup （基本已不用），samtools markdup（更新后的）和picard最常用。rmdup效果不怎么好，而且如果有遇到相同位置的reads， 会优先选择质量高的那一条read。picard与samtools markdup效果相似（仿佛调用的同一个？并不确定）。都可以标记重复，也可以选择直接去掉。以下是用法：

samtools markdup -@ 8 -r test.bam filter_test.bam # -r是直接去掉重复，不加是直接标记

picard去重有三种方式可选，在DUPLICATE_SCORING_STRATEGY参数中，分别是SUM_OF_BASE_QUALITIES, TOTAL_MAPPED_REFERENCE_LENGTH和RANDOM。即当有重复时分别选择留下总碱基质量最高的、匹配上参考基因组最长的和随机。

picard MarkDuplicates I=test.bam  O= filter_test.bam M=dup_metrics.txt REMOVE_DUPLICATES=true

4.peak calling peaks是reads信号比较强的区域，也就是我们找到的转录因子或者组蛋白修饰最有可能结合的地方。call peaks仍然有不少软件，比较常用的是MACS2和Hotspot2。 示例：

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

针对不同的数据考虑用不同的参数。

5.下游分析 (downstream analysis)

分析完之后下游可以做的事情很多，视情况而定。可以同时分析DNase-seq或者ATAC-seq的数据，看转录因子与染色质开放区的关系；或者Homer等工具注释peaks，看不同转录因子/组蛋白修饰之间的关系，或者分析TF的target gene。也可以用MEME等做motif分析。

欢迎关注公众号！

生信编程日常