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

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

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

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

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

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

我写ChIPseeker的时候，我做的物种是人，ChIPseeker在线一周就有剑桥大学的人写信跟我说在用ChIPseeker做果蝇，在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

Continue reading

这一次讲解非常重要的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对象），你就可以跑注释了，然后就可以出结果了。

Continue reading

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 - 子元件的起始位点

Continue reading

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技术。

Continue reading

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

我当时被要求做ChIPseq分析是为他人做嫁衣，而且是完全白干那种，但做为学生，白干也得干。

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

Continue reading

发现Youtube上有一个视频叫Evolution of clusterProfiler, 是Landon Wilkins用Gource做的。于是我也来玩一下，看一下自己这几年码代码的过程。

Continue reading