CS0: ChIPseq从入门到放弃

接下来要出一个ChIPseq系列，讲一讲ChIPseq和我的ChIPseeker包，从入门到放弃是我自己的个人写照。我做ChIPseq总共也就3个月的时间，做的事情并不多，在一知半解的情况下写下了ChIPseeker包。

正如我在《话题投票》里说的，我当时被要求做ChIPseq分析是为他人做嫁衣，而且是完全白干那种，但做为学生，白干也得干。

当时一开始使用ChIPpeakAnno做注释，但用UCSC genome browser检验结果的时候，发现对不上。在对ChIPpeakAnno包不满意的情况下，开始着手写ChIPseeker，其实在使用ChIPpeakAnno的时候，我就有写代码对结果做一些可视化，所以未有ChIPseeker先有ChIPseeker的部分可视化功能。当时写了篇博客文说ChIPpeakAnno的问题，一个月后就在Bioconductor上发表了ChIPseeker，这包完全是我半夜在宿舍里写出来的。

当时还在生物系，被我炒掉的前老板每天要求必须起码在实验室待够12小时，我每天都待到10点半左右才回宿舍，日常在实验室里啥都干不了，白天各种瞎折腾，晚上还要陪他聊天，但说来说去，每天几乎都差不多，无非是他很牛逼，我们这帮人读他phd实在太幸运，日复一日传销式洗脑。而我因为结婚了，家又离得近，周末回家，白天经常多一段单独对我的洗脑，做为一个PhD学生，在发表文章之前是不能够有周末的。每天10半从实验室里出来，回到宿舍11点，跟老婆打电话再洗澡，12点。然后从12点开始写代码到2点睡觉，才有了这个包。

虽然是一知半解的时候开发的，但还是受到大家的欢迎，半年前Matt邀请我去人大做报告时，也专门提到了ChIPseeker。

也有美国的助理教授，跟我要paper，说是上课的时候，要给学生读的，这广告效果我给满分。

文章发表了一年，已经被33篇文章引用，其中不乏有影响因子比较高的杂志：

下面是其中一些引用文章的图：

虽然ChIPseeker是我写给自己做ChIPseq注释的，但Ming Tang (https://github.com/crazyhottommy/ChIP-seq-analysis）用它去做DNA breakdown注释，当然像lincRNA注释也是有人做并且完全是支持的。有一些我以前从没在文档里提到的东西，也应该会在这个系列里写出来。

这个系列基本上是围绕着ChIPseeker的功能而来，名副其实从入门到放弃，因为我自己也是入了门然后放弃，如果想看从入门到精通的，这显然不适合你。

然而今天只是个剧透，敬请期待。

CS1: ChIPseq简介

ChIP是指染色质免疫沉淀，它通特异结合抗体将DNA结合蛋白免疫沉淀，可以用于捕获蛋白质（如转录因子，组蛋白修饰）的DNA靶点。这技术存在非常久了，在二代测序之前，结合microarray，它的名字叫ChIP-on-chip，二代测序出来之后，显而易见的，免疫沉淀拉下来的DNA拿去NGS测序，这必然是下一代的ChIP技术，优点也是显而易见的，不再需要设计探针（往往存在着一定的偏向性）。所以NGS出来以后，不差钱的牛逼实验室显然占据上风，谁先做出来，谁就定义了新技术。这是有钱人的竞赛，没钱的只能等着技术烂大街的时候跟风做。 这是显而易见的下一代技术，外加技术上完全是可行的，所以这是一场单纯的时间竞赛，于是几乎同时出来CNS文章，基本上谁也不比谁差地同时扔出来。

• Johnson DS, Mortazavi A et al. (2007) Genome-wide mapping of in vivo protein–DNA interactions. Science 316: 1497–1502
• Robertson G et al.(2007) Genome-wide profiles of STAT1 DNA association using chromatin immunoprecipitation and massively parallel sequencing. Nature Methods 4: 651–657
• Schmid et al. (2007) ChIP-Seq Data reveal nucleosome architecture of human promoters. Cell 131: 831–832 2007年来自三个不同的实验室，几乎是同时间出来（最长差不了3个月），分别发CNS，一起定义了这个ChIPseq技术。

这个技术分为4步：

• Cross-linking
• Sonication
• IP
• Sequencing DNA和蛋白质交联(cross-linking)，超声(sonication)将染色体随机切割，利用抗原抗体的特异性识别(IP)，把目标蛋白相结合的DNA片段沉淀下来，反交联释放DNA片段，最后是测序(sequencing)。 一个典型的分析流程如下：

测序之后，我们当然首先需要做质量控制，然后就是做mapping，拿到这些DNA片段在染色体上的位置信息，ChIPseq的数据我们还需要做peak calling，把背景噪声去掉，比如上图中使用MACS做peak calling，这样我们就得到了protein binding site (peak)，就可以做下游的分析，比如可视化、相关的基因（比如最近的基因、宿主基因）、Motif分析等等。

Peak annotation做的就是binding site的相关基因注释，在讲解ChIPseeker的注释功能之前，下次先讲解一下peak calling的输出，BED文件。

CS2: BED文件

文件格式

BED的全称是Browser Extensible Data，顾名思义是为genome browser设计的，大名鼎鼎的bedtools可不是什么「床上用品」。

BED包含有3个必须的字段和9个可选字段。

三个字段包括：

• 1 chrom - 染色体名字
• 2 chromStart - 染色体起始位点
• 3 chromEnd - 染色体终止位点

这里必须指出的是chromStart是起始于0，而不是1。很多分析软件都忽略了这一点，会有一个碱基的位移，据我说知Homer和ChIPseeker没有这个问题，而像peakAnalyzer, ChIPpeakAnno等都有位移的问题。

可选的9个字段包括：

• 4 name - 名字
• 5 score - 分值(0-1000), 用于genome browser展示时上色。
• 6 strand - 正负链，对于ChIPseq数据来说，一般没有正负链信息。
• 7 thickStart - 画矩形的起点
• 8 thickEnd - 画矩形的终点
• 9 itemRgb - RGB值
• 10 blockCount - 子元件（比如外显子）的数目
• 11 blockSizes - 子元件的大小
• 12 blockStarts - 子元件的起始位点

一般情况下，我们只用到前面5个字段，这也是做peak calling的MACS输出的字段。

其中第5个字段，MACS的解释是这样子的：

The 5th column in this file is the summit height of fragment pileup.

是片段堆积的峰高，这也不难理解，为什么我在ChIPseeker是画peak coverage的函数covplot要有个weightCol的参数了。

数据可视化

从名字上看，它是为genome browser而生，相应的，ChIPseeker实现了covplot来可视化BED数据。

covplot支持直接读文件出图：

library(ChIPseeker)
library(ggplot2)

files <- getSampleFiles()
covplot(files[[4]])

支持GRanges对象，同时可以多个文件或者GRangesList

peak=GenomicRanges::GRangesList(CBX6=readPeakFile(files[[4]]),CBX7=readPeakFile(files[[5]]))
covplot(peak, weightCol="V5") + facet_grid(chr ~ .id)

支持可视化某个窗口

covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7)) + facet_grid(chr ~ .id)

拿到数据后，我们首先会可视化看一下数据，接下来就会想知道这些peak都和什么样的基因有关，这将在下次讲解，如何做注释。（昨天biostars上就有人在问用ChIPseeker给schizosaccharomyces pombe裂殖酵母做注释，大家可去围观）

CS3: peak注释

这一次讲解非常重要的peak注释，注释在ChIPseeker里只需要用到一个函数annotatePeak，它可以满足大家各方面的需求。

输入

当然需要我们上次讲到的BED文件，ChIPseeker自带了5个BED文件，用getSampleFiles()可以拿到文件的全路径，它返回的是个named list，我这里取第4个文件来演示。annotatePeak的输入也可以是GRanges对象，你如果用R做peak calling的话,直接就可以衔接上ChIPseeker了。

> require(ChIPseeker)
> f = getSampleFiles()[[4]]

巧妇难为无米之炊，就像拿到fastq要跑BWA，你需要全基因组的序列一样，做注释当然需要注释信息，基因的起始终止，基因有那些内含子，外显子，以及它们的起始终止，非编码区的位置，功能元件的位置等各种信息。 很多软件会针对特定的物种去整理这些信息供软件使用，但这样就限制了软件的物种支持，有些开发者写软件本意也是解决自己的问题，可能对自己的研究无关的物种也没兴趣去支持。 然而ChIPseeker支持所有的物种，你没有看错，ChIPseeker没有物种限制，当然这是有前提的，物种本身起码是有基因的位置这些注释信息，不然就变无米之炊了。 这里我们需要的是一个TxDb对象，这个TxDb就包含了我们需要的各种信息，ChIPseeker会把信息抽取出来，用于注释时使用。

> require(TxDb.Hsapiens.UCSC.hg19.knownGene)
> txdb = TxDb.Hsapiens.UCSC.hg19.knownGene
> x = annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb)
>> loading peak file...                 2017-03-09 11:29:18 PM 
>> preparing features information...         2017-03-09 11:29:18 PM 
>> identifying nearest features...         2017-03-09 11:29:19 PM 
>> calculating distance from peak to TSS...     2017-03-09 11:29:20 PM 
>> assigning genomic annotation...         2017-03-09 11:29:20 PM 
>> assigning chromosome lengths             2017-03-09 11:29:42 PM 
>> done...                     2017-03-09 11:29:42 PM

这里需要注意的是，启动子区域是没有明确的定义的，所以你可能需要指定tssRegion，把基因起始转录位点的上下游区域来做为启动子区域。 有了这两个输入（BED文件和TxDb对象），你就可以跑注释了，然后就可以出结果了。

输出

如果在R里打输出的对象，它会告诉我们ChIPseq的位点落在基因组上什么样的区域，分布情况如何。

> x
Annotated peaks generated by ChIPseeker
1331/1331  peaks were annotated
Genomic Annotation Summary:
             Feature  Frequency
9           Promoter 48.1592787
4             5' UTR  0.7513148
3             3' UTR  4.2073629
1           1st Exon  0.7513148
7         Other Exon  3.9068370
2         1st Intron  3.6814425
8       Other Intron  7.7385424
6 Downstream (<=3kb)  1.1269722
5  Distal Intergenic 29.6769346

如果我想看具体的信息呢？你可以用as.GRanges方法，这里我只打印前三行：

Bioconductor里有很多包是针对GRanges对象的，这样方便你在R里做后续的处理，如果你说你不懂这些，只想输出个Excel表格。那么也很容易，用as.data.frame就可以转成data.frame，然后你就可以用write.table输出表格了。

两种不同的注释

这里注释有两种，一种是genomic annotation (也就是annotation这一列）还有就是nearest gene annotation（也就是多出来的其它列）。

经常有人问我问题，把这两种搞混。genomic annotation注释的是peak的位置，它落在什么地方了，可以是UTR，可以是内含子或者外显子。

而最近基因是peak相对于转录起始位点的距离，不管这个peak是落在内含子或者别的什么位置上，即使它落在基因间区上，我都能够找到一个离它最近的基因（即使它可能非常远）。 这两种注释的策略是不一样的。针对不同的问题。第一种策略peak所在位置，可能就是调控的根本，比如你要做可变剪切的，最近基因的注释显然不是你关注的点。

而做基因表达调控的，当然promoter区域是重点，离结合位点最近的基因更有可能被调控。

最近基因的注释信息虽然是以基因为单位给出，但我们针对的是转录起始位点来计算距离，针对于不同的转录本，一个基因可能有多个转录起始位点，所以注释是在转录本的水平上进行的，我们可以看到输出有一列是transcriptId.

第三种注释

两种注释有时候还不够，我想看peak上下游某个范围内（比如说-5k到5k的距离）都有什么基因，annotatePeak也可以做到。 你只要传个参数说你要这个信息，还有什么区间内，就可以了。

x = annotatePeak(f[[4]], tssRegion=c(-1000, 1000), TxDb=txdb, addFlankGeneInfo=TRUE, flankDistance=5000)

输出中多三列: flank_txIds, flank_geneIds和flank_gene_distances，在指定范围内所有的基因都被列出。

基因注释

对于通常情况找最近基因的策略，最近基因给出来了，但都是各种人类不友好的ID，我们不能把一切都交给计算机，输出的结果我们还是要看一看的，能不能把基因的ID换成对人类友好的基因名，并给出描述性的全称，这个必然可以有。

只需要给annotatePeak传入annoDb参数就行了。 如果你的TxDb的基因ID是Entrez，它会转成ENSEMBL，反之亦然，当然不管是那一种，都会给出SYMBOL，还有描述性的GENENAME.

CS4：关于ChIPseq注释的几个问题

为什么我要用某个基因组版本？

在上一篇文章中，我用了TxDb.Hsapiens.UCSC.hg19.knownGene。 hg19的TxDb, 或者有人就要问了，为什么不用hg38？

这个问题，不是说要用那一个，不能用那一个。而是你必须得用某一个，这取决于你最初fastq用BWA/Bowtie2比对于某个版本的基因组，你最初用了某个版本，后面就得用相应的版本，不能混，因为不同版本的位置信息有所不同。

当然如果要（贵圈喜欢的）强搞，也不是不可以，你得有chain file，先跑个liftOver，实际上就是在两个基因组版本之间做了位置转换。

为什么说ChIPseeker支持所有物种？

背景注释信息用了TxDb就能保证所有物种都支持了？我去哪里找我要的TxDb?

我写ChIPseeker的时候，我做的物种是人，ChIPseeker在线一周就有剑桥大学的人写信跟我说在用ChIPseeker做果蝇，在《CS2: BED文件》一文中，也提到了最近有人在Biostars上问用ChIPseeker做裂殖酵母。

首先Bioconductor提供了30个TxDb包，可以供我们使用，这当然只能覆盖到一小部分物种，我们的物种基因组信息，多半要从UCSC或者Ensembl获得，我敢说支持所有物种，就是因为UCSC和ensembl上所有的基因组都可以被ChIPseeker支持。

因为我们可以使用GenomicFeatures包函数来制作TxDb对象：

• makeTxDbFromUCSC： 通过UCSC在线制作TxDb
• makeTxDbFromBiomart: 通过ensembl在线制作TxDb
• makeTxDbFromGRanges：通过GRanges对象制作TxDb
• makeTxDbFromGFF：通过解析GFF文件制作TxDb

比如我想用人的参考基因信息来做注释，我们可以直接在线从UCSC生成TxDb:

require(GenomicFeatures)
hg19.refseq.db <- makeTxDbFromUCSC(genome="hg19", table="refGene")

比如最近在biostars上有用户问到的，做裂殖酵母的注释，我们可以下载相应的GFF文件，然后通过makeTxDbFromGFF函数生成TxDb对象，像下面的命令所演示，spombe就是生成的TxDb，就可以拿来做裂殖酵母的ChIPseq注释。

download.file("ftp://ftp.ebi.ac.uk/pub/databases/pombase/pombe/Chromosome_Dumps/gff3/schizosaccharomyces_pombe.chr.gff3", "schizosaccharomyces_pombe.chr.gff3")require(GenomicFeatures)
spombe <- makeTxDbFromGFF("schizosaccharomyces_pombe.chr.gff3")

所以我敢说，所有物种都支持。像Johns Hopkins出品的CisGenome就只支持到12个物种而已，极大地限制了它的应用。

ChIPseeker有什么不能注释的吗？

这个我还没想到，像CpG是不支持的，但也有人「黑」出来了：

当然他的做法是把CpG也整合进去，如果你单纯只想看那些peak落在CpG上，或者说离CpG最近，不需要「黑」也能做到的，因为annotatePeak的背景注释信息除了TxDb之外，其实它还可以是自定义的GRanges对象，这保证了用户各种各样的需求，因为TxDb也不是万能的，如果能自定义，比如说我就只想看蛋白的结合位点会不会在内含子和外显子的交界处，再比如说我做的并不是编码蛋白的基因表达调控，而是非编码RNA，那么我想要用lncRNA的位置信息来做注释。像这样的需求，ChIPseeker都是可以满足的。

可以按正负链分开注释吗？

上一篇文章《CS3: peak注释》中就有人问了，能否同时给出正负链上最近的基因。首先ChIPseq数据通常情况下是没有正负链信息的（有特殊的实验可以有），annotatePeak函数有参数是sameStrand，默认是FALSE，你可以给你的peak分别赋正负链，然后传入sameStrand=TRUE，分开做两次，你就可以分开拿到正链和负链的最近基因。

最近基因位置是相对于TSS的，能否相对于整个基因？

答案也是可以的！

首先如果peak和TSS有overlap，genomic annotation就是promoter，而最近基因也是同一个，所以在这种情况下，两种注释都指向同一个基因，可以说信息是冗余的，能不能不要冗余信息？这个是可以的，你可以传入参数ignoreOverlap=TRUE，那么最近基因就会去找不overlap的。

最近基因是相对于TSS，如果和TSS有overlap，距离是0, 必须是最近。回到标题的问题，如果我想说只要和基因有overlap就是最近基因，这种情况其实你的最近基因就是host gene，也就是annotation这个column给出来的是相对应的，我们就想找peak所在位置的基因信息，那么这当然也是可以的，默认参数overlap=”TSS”, 如果改为overlap=”all”，它看的就是整个基因而不是TSS，当然distanceToTSS也还是会计算，如果overlap的不是TSS，而是基因体里，并不会因为而设为0。

如果我只想注释上游或者是下游的基因呢？

当然也可以，我们有ignoreUpstream和ignoreDownstream参数，默认都是FALSE。随便你想看上游还是下游，都可以。

为什么要有这么多参数？

我在前一篇，只讲了输入，输出，你知道两个输入，会看输出，你就可以做ChIPseq注释了，非常简单。但是我不能把annotatePeak能做的全列出来，会让大家觉得复杂。（而且简单的情况是最常见的行为）

在大家知道输入输出，觉得简单之后，再讲一讲，它有一些参数，可以应对别的情况，这些情况可能并不是我们做ChIPseq所需要的，但不同参数的灵活组合，是可以解决和应对不同的需求的。

比如说, DNA是可以断的，如果我们要对（DNA breakpoints）断点位置做注释，优先是overlap基因，再者是上游，你会发现很多软件都不能做了，但ChIPseeker可以。

做软件关键是要注重细节，不单注重自己需求的细节，更主要是要注重别人的细节