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

在上一篇文章中,我用了TxDb.Hsapiens.UCSC.hg19.knownGenehg19TxDb, 或者有人就要问了,为什么不用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注释

这一次讲解非常重要的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

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


Continue reading

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

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

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


Continue reading

Author's picture

Guangchuang Yu

a senior-in-age-but-not-senior-in-knowledge bioinformatician

PhD student

Hong Kong