制作了一个banner,用来插入到文后,方便阅读到最后的时候,直接扫描关注。

这图当然用PS一下就可以得到,无非是拼图和加点文字。但做为天天写代码画图的人来说,必然是要纯代码来产生的,而且做科学的人,讲究自动化、可重复性。

就像有些人不理解data scentist为什么讨厌excel一样,觉得无非是工具,没什么好搞阵营的。试想一下,一个分析流程中间有一步要用到excel,需要人工去点点鼠标,这对于讲究自动化、可重复性的data scientist来说是不可接受的。

加载中文字体

为了打几个中文字,需要先加载字体,这里我用showtext。

require(showtext)
font.add("heiti", "/Library/Fonts/华文黑体.ttf")
showtext.auto()


Continue reading

比如下面的代码:

require(ggplot2)
d <- data.frame(x=c(0, 0.002, 0.00575), y = 1:3)
p <- ggplot(d, aes(x, y)) + geom_point() + xlab(NULL) + ylab(NULL)
print(p)

上面图中x轴的文本0.006,这个数字中的6几乎看不到,因为一半过界了。


Continue reading

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

在上一篇文章中,我用了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

News in emojifont

面向对象有多种实现方式,R里面就有3种,class-based, method-based, object-based,R6与C++/JAVA一样是class-based的,S3/S4是method-based的,还有一种是object-based的,这在proto包中实现,很多人可能没听说过,但是ggplot2你们一定听过,ggplot2就是object-based的实现,它现在是自己的定制实现,称之为ggproto。

emojifont就是用proto实现的,属于我的练手之作,很高兴深受大家的喜欢。


Continue reading

how to bug author

http://blog.devtang.com/2017/03/05/how-to-get-help/ 这篇文章讲问问题的礼仪,会问问题的人容易得到别人的帮助,并不是作者拽,而是对着一群不会问问题,且一副理所当然的用户,早就有了深深的无力感,也看看我写的这篇吧,https://guangchuangyu.github.io/2016/07/how-to-bug-author/.

这是前几天的文字推送,不知道大家有没有阅读这两篇博客文?如何提问,这是一项重要的技能,很遗憾很多人并没有这项技能!

就在上周,我正好就在微信群里跟某位老师说不要给我个截屏,说能不能在ggtree里实现某功能。懂的人知道,我只是希望大家用正确的方式来做事情,不懂的人,当然觉得我拽,JJYY什么的。


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

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


Continue reading

Author's picture

Guangchuang Yu

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

Postdoc researcher

Hong Kong