ChIP-Seq分析和作用

1：ChIP-Seq数据是基因组特异性富集的序列的测序结果，包括组蛋白修饰ChIP-Seq(H3K4me3/启动子相关/narrowpeak、H3K4me1/增强子相关/narrowpeak、H3K27ac/增强子相关/broadpeak)、转录因子ChIP-Seq(CTCF/绝缘子相关/narrowpeak、pol II/转录起始/narrowpeak)、DNA富集序列(DNase-Seq/弱DNA酶消化/活性区域、MNase-Seq/强DNA酶消化/核小体不活跃区域、ATAC-Seq//前两者的结果的集合)。

　　通过互补染色质分析实验分析的基因组位点揭示了染色质结构的不同方面：ChIP-seq显示特异性转录因子（TF）的结合位点; DNase-seq，ATAC-seq和FAIRE-seq显示开放染色质的区域;和MNase-seq鉴定良好定位的核小体。在ChIP-seq中，特异性抗体用于直接或通过包含靶因子的复合物中的其他蛋白质提取结合至靶蛋白的DNA片段。在DNase-seq中，染色质被DNA酶I内切核酸酶轻微消化。大小选择用于富集在DNA对DNA酶I攻击高度敏感的染色质区域产生的片段（在初期会生成包含各种长度的DNA小片段，但是一般来书保留100~300bp长度的小片段建库测序）。 ATAC-seq是DNase-seq的替代方法，其使用工程改造的Tn5转座酶来切割DNA并将引物DNA序列整合到切割的基因组DNA中（即，标记）。微球菌核酸酶（MNase）是内切核酸外切酶，其连续地消化DNA直到达到阻塞（和DNA酶相比（DNase-seq），属于强切，开放的区域全部都被消化），例如核小体。在FAIRE-seq中，甲醛用于交联染色质，并且苯酚 - 氯仿用于分离剪切的DNA。

　　详细介绍可参考(http://www.nature.com/nrg/journal/v15/n11/fig_tab/nrg3788_F1.html)

　　

 

 

2：ChIP-Seq数据的作用：a:构建物种的epigenome，利用chromHMM将基因组分成一个一个的区域；b:与交互数据(HiC/chia-pet)联合分析；c:和RNA-Seq联合分析(chirp-seq)。

 

3：ChIP-Seq数据的分析流程：

　　1：预处理，步骤与RNA-Seq的一致，详情见RNA-Seq分析的1、2两步。

　　2：比对：DNA数据比对软件用的比较多的是bwa和bowtie，bwa比对结果相较于bowtie的比对结果更加准确，但是跑得慢，bowtie正好相反。我一般是用bwa比对，很少使用bowtie比对，慢不了多少。而且方便下游分析(比如使用GATK  call snv/indel，bwa的结果明显比bowtie的要好)。

　　bwa比对提供了两种主流的比对aln和mem，一般序列长度小于70bp选用aln比对，大于70bp的时候选用mem。

　　首先第一步都是根据参考基因组建索引：bwa index genome.fa

　　当长度小于70bp时：

　　针对于单端数据：bwa aln -t 4 -f file.sai genome.fa file.fastq

　　　　　　　　　　bwa samse -n 1 -f file.sam genome.fa file.sai file.fastq

　　针对于双端数据：bwa aln -t 4 -f file1.sai genome.fa file1.fastq

　　　　　　　　　　bwa aln -t 4 -f file2.sai genome.fa file2.fastq

　　　　　　　　　　bwa sampe -n 1 -f file.sam genome.fa file1.sai file1.fastq file2.sai file2.fastq

　　　　　　　　　　-t 线程数目　　-f 输出文件名　　-n 最大输出条目(说白了就是多位点比对，为了省事可以考虑设置为1，就只输出一个比对的位点，一般不超过3)

　　　　　　　　　　可加参数 -r ，头文件是为样本加入分组信息。为了下游分析的正常进行，我一般不在此处加，在后续需要的时候用picard软件加上头文件信息来分组

　　当长度大于70bp时：

　　针对于单端数据：bwa mem -t 4 genome.fa file.fastq > file.sam

　　针对于双端数据：bwa mem -t 4 genome.fa file1.fastq file2.fastq > file.sam

　　　　　　　　　　针对于mem比对一般不考虑加入其他的参数。

　　　　　　　　　　可加参数 -r ，头文件是为样本加入分组信息。为了下游分析的正常进行，我一般不在此处加，在后续需要的时候用picard软件加上头文件信息来分组

 

　　3：sam2bam：相对来说bam格式作为sam文件二进制的形式，存在占用内存更小、软件读取效率更高等优点，所以在此处将sam文件转成bam文件

　　sam文件转bam文件并排序：samtools view -Su file.sam | samtools sort -o file.bam -(samtools view -q 20 -Su file.sam | samtools sort -o file.bam - 或者可以直接在该步骤完成质量过滤但不建议)。

　　去除duplication：java -jar ~/Tools/picard-tools-1.119/MarkDuplicates.jar I=file.bam O=file.rmd.bam  CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT M=file.metrics

　　可以加入的参数：REMOVE_DUPLICATES=true。删除duplication的reads而不是修改该reads的flag(第二列的值)。但是为了保证数据的完整性不建议这么做。

　　提取高质量的比对结果：samtools view -F 1024 -q 20 -o file.flt.bam file.bam

　　此时的file.flt.bam文件就可以导入到可视化软件(IGV)中查看了，而不应该导入file.bam，因为file.bam中包含了很多的duplication和低比对的reads，导致和最终的结果有差异，但是查看file.flt.bam就不会。此外，file.bam主要是为了后续出什么问题方便来检查，而且可以作为其他分析的输入文件，最重要的是bam文件的大小是小于原始数据(fastq)的，可以用来保存原始数据，在需要的时候可以通过bedtools将ban文件转成fastq文件。

　　额外的分析：在该步骤的过程中可以加入样本重复性的分析，deeptools中的multiBamSummary和plotCorrelation可以用来检查样本的重复性

 

　　4：peak calling：该步骤是ChIP-Seq中最重要的一个步骤，做的好坏直接影响最终的结果，在peak calling时候重要有两个软件：macs2和F-seq。一般来说macs2用的比较多，针对于各种不同数据之间的参数选择：a:转录因子和组蛋白修饰中的H3K4me3、H3K4me3直接用默认的参数运行，生成narrowpeak文件；b:除了H3K4me3和H3K4me1以外的组蛋白修饰的ChIP-Seq 加上 -broad 参数即可；c:针对DNase-Seq数据加上参数 --nomodel --shift 100 --extsize 200；针对MNase-Seq数据加上参数 --nomodel --shift 37 --extsize 73。结果文件说明见官网。

 

　　5：peak annotation：获得peak文件后能做的分析就比较多了，其中的一个就是对peak文件的注释。注释软件一般用homers和chipseeker，不过后来仔细比较了两个软件注释的结果，发现homer注释的结果不是很正确，相反chipseeker要好一些。chipseeker是一个R包，需要R3.3.0以上的版本，需要安装clusterProfiler和GenomicFeatures包。

注释的代码见下(chipseeker还可以做很多事，可以去官方网站查看):

library(ChIPseeker)
library(clusterProfiler)
library(GenomicFeatures)
files <- list(randA=c("/home/xyhuang/program/randomA20.peak"),randB=c("/home/xyhuang/program/randomB20.peak"))
txdb <- makeTxDbFromGFF("hg19.gtf",format="gtf")
peakAnno <- annotatePeak(files[[1]], tssRegion=c(-3000, 3000), TxDb=txdb)
pdf("randA.pdf")
plotAnnoPie(peakAnno)
dev.off()
data <- as.data.frame(peakAnno)
aa <- c(7,8,9,10,11,12)
data <- data[-aa]
write.table(data,"randA.txt",row.names=F,quote = F,sep="\t")

peakAnno <- annotatePeak(files[[2]], tssRegion=c(-3000, 3000), TxDb=txdb)
pdf("randB.pdf")
plotAnnoPie(peakAnno)
dev.off()
data <- as.data.frame(peakAnno)
aa <- c(7,8,9,10,11,12)
data <- data[-aa]
write.table(data,"randB.txt",row.names=F,quote = F,sep="\t")

 

　　6：检查peak是否富集：使用bedtools random 随机生成一组peak，然后重复步骤五，比较两个结果是否有差异，差异越大表示富集的程度越高越好。

　　samtools fadix genome.fa

　　awk '{print $1"\t"$2;}' genome.fa.fai > genome.size

　　bedtools random -l 100 -n 100000 -g genome.size

　　-l peak的宽度　　-n 生成的peak数目

 

 　　7：motif分析：针对与转录因子和组蛋白修饰，目前认为是蛋白结构特异性的识别基因组上的一段特异的序列，而这段序列是具有保守性的。可以用来鉴定基因组上的调控元件

　　motif分析主要也是两个软件homer和meme，根据注释结果，我认为meme寻找的结果更加可靠，而且meme专门为ChIP-Seq数据开发了一个meme-chip软件包。但是考虑到需要将寻找到的motif与motif库进行比较，所以一般推荐网页上做。

　　本地化的后续补充：

 

　　8：未完。。。。。。