日常瞎掰

  对于ChIP-seq、ATAC-seq等这样捕获基因组富集区域(即分析结果中peak)的技术，大家多多少少应该有所耳闻。在分析这类测序数据的时候，必不可少的步骤就是将peak注释到基因组上，以便了解peak出现在哪些基因的周边区域，从而发现生物学上的意义。目前，注释peak的软件不在少数，如ChIPseeker、homer2等。今天我们主要来说说如何利用ChIPseeker绘制peak的分布饼图和条形图。

注释

  ChIPseeker的用法相当简单，这里就顺便简单介绍一下，下面使用包里面的数据演示一下：

1. 使用R里面的数据库

library(ChIPseeker) library(TxDb.Hsapiens.UCSC.hg19.knownGene) library(org.Hs.eg.db) files <- getSampleFiles() peakanno <- annotatePeak(files[[1]], TxDb=txdb, annoDb = "org.Hs.eg.db") peakanno Annotated peaks generated by ChIPseeker 1331/1331 peaks were annotated Genomic Annotation Summary: Feature Frequency 9 Promoter (1-2kb) 6.0856499 10 Promoter (<=1kb) 48.1592787 11 Promoter (2-3kb) 4.5078888 4 5' UTR 0.3005259 3 3' UTR 2.1036814 1 1st Exon 0.6010518 7 Other Exon 2.4042074 2 1st Intron 3.6814425 8 Other Intron 4.2824944 6 Downstream (<=300) 0.9767092 5 Distal Intergenic 26.8970699

  如上面所示，调用一下annotatePeak函数注释就完成了。其中，TxDb.Hsapiens.UCSC.hg19.knownGene数据库是必须的，org.Hs.eg.db数据库为可选项，如果提供了注释结果里面会多出ENSEMBL、SYMBOL、GENENAME三列信息。当然，annotatePeak还有其他参数可以根据需要设置，这里就不展示了。

2. 使用gtf文件

library(GenomicFeatures) txdb <- makeTxDbFromGFF("hg19.gtf") # 或者在线获取 txdb <- makeTxDbFromUCSC(genome="hg19", table="refGene") peakanno <- annotatePeak(files[[1]], TxDb=txdb)

  这种方式注释起来也非常方便。不论使用哪种注释方法，最终生成的注释结果都是以csAnno类对象的方式存储，后续操作就可以基于这个对象了。

过滤

  注释后，我们可以根据需求对其进行过滤，筛选想关注的区域来看看分布情况，比如，这里只保留Promoter、Exon、Intron三个大类的结果，可以看出过滤后每个Feature的比例会重新计算：

annofilter <- subset(peakanno ,annotation %in% grep('Promoter|Exon|Intron',annotation,value=T)) annofilter Annotated peaks generated by ChIPseeker 928/928 peaks were annotated Genomic Annotation Summary: Feature Frequency 5 Promoter (1-2kb) 8.728448 6 Promoter (<=1kb) 69.073276 7 Promoter (2-3kb) 6.465517 1 1st Exon 0.862069 3 Other Exon 3.448276 2 1st Intron 5.280172 4 Other Intron 6.142241

  当然，在过滤之前，我们要先知道都有哪些列可以用来过滤，可以使用下面的命令来查看：

head(peakanno@anno) GRanges object with 6 ranges and 11 metadata columns: seqnames ranges strand | V4 V5 annotation <Rle> <IRanges> <Rle> | <character> <numeric> <character> [1] chr1 815093-817883 * | MACS_peak_1 295.76 Promoter (2-3kb) [2] chr1 1243288-1244338 * | MACS_peak_2 63.19 Promoter (<=1kb) [3] chr1 2979977-2981228 * | MACS_peak_3 100.16 Promoter (<=1kb) [4] chr1 3566182-3567876 * | MACS_peak_4 558.89 Promoter (<=1kb) [5] chr1 3816546-3818111 * | MACS_peak_5 57.57 Promoter (<=1kb) [6] chr1 6304865-6305704 * | MACS_peak_6 54.57 Promoter (<=1kb) geneChr geneStart geneEnd geneLength geneStrand geneId <integer> <integer> <integer> <integer> <integer> <character> [1] 1 803451 812182 8732 2 284593 [2] 1 1243994 1247057 3064 1 126789 [3] 1 2976181 2980350 4170 2 440556 [4] 1 3547331 3566671 19341 2 49856 [5] 1 3816968 3832011 15044 1 100133612 [6] 1 6304252 6305638 1387 1 390992 transcriptId distanceToTSS <character> <numeric> [1] uc001abt.4 -2911 [2] uc001aed.3 0 [3] uc001aka.3 0 [4] uc001ako.3 0 [5] uc001alg.3 0 [6] uc009vly.2 613 ------- seqinfo: 24 sequences from hg19 genome

可视化

1. barplot
     使用ChIPseeker注释后，我们可以使用其内置的函数plotAnnoBar来展示样本的peak分布情况，适合一个或多个样本。这里我们选取两个样本来演示：

alist <- sapply(files[4:5], function(x){annotatePeak(x, TxDb=txdb)}) p <- plotAnnoBar(alist) p

结果如下：

2. pieplot
     单个样本的时候也可以选择饼图来展示结果，调用函数plotAnnoPie即可绘图：

plotAnnoPie(peakanno)

结果如下：

可视化函数借用

  其实，我们也可以只用ChIPseeker的可视化函数来展示自己的结果而不用其来注释。如果我们已经有了注释的结果，可以将统计结果整理成类似上面所示展示的Genomic Annotation Summary数据框格式，然后借助ChIPseeker包的绘图函数来绘图：

1. barplot

df1 <- data.frame(Feature=factor(c('intron','exon','promoter','utr')),Frequency=prop.table(sample(10:100,4)),sample='sampleA') df2 <- data.frame(Feature=factor(c('intron','exon','promoter','utr')),Frequency=prop.table(sample(10:100,4)),sample='sampleB') df <- rbind(df1,df2) df Feature Frequency sample 1 intron 33.333333 sampleA 2 exon 27.147766 sampleA 3 promoter 30.240550 sampleA 4 utr 9.278351 sampleA 5 intron 14.798206 sampleB 6 exon 33.183857 sampleB 7 promoter 8.968610 sampleB 8 utr 43.049327 sampleB p <- ChIPseeker:::plotAnnoBar.data.frame(df, categoryColumn='sample') p

结果如下：

  如果觉得这个barplot风格不错，想用于绘制自己的数据，借用一下函数是不是很方面快捷。准备数据的时候只需注意两点：一是必须有Feature、Frequency两列，列名不能错，分别用于指定类别及比例；二是Feature的数据类型为factor。第三列做为样本列，名称随意，如果没有categoryColumn设置为1。

2. pieplot
     准备饼图数据时，只需要有Feature、Frequency两列既可，有没有其他列不会影响绘图：

df1 Feature Frequency sample 1 intron 33.333333 sampleA 2 exon 27.147766 sampleA 3 promoter 30.240550 sampleA 4 utr 9.278351 sampleA col <- ChIPseeker:::getCols(nrow(df1)) ChIPseeker:::annoPie(df1,col=col,legend.position='rightside')

结果如下：

3. 自定义pieplot
     其实，ChIPseeker绘制的饼图个人不是很喜欢，有两个小的因素。首先，其绘制饼图时不是基于ggplot2，而是基于基础函数pie，扩展性没有那么方便；再者，绘出的图和图例的间距比较大。虽然这两点可以忽略不计，但是如果能绘出更接近成品图的效果那后续不是更省时省力么。所以，个人借用ChIPseeker的配色方法包装了一个基于ggplot2的绘制饼图的函数，可以接受ChIPseeker的注释对象或者含有Feature、Frequency两列的数据框。下面一起来看看绘图效果：

plotpie <- function(data){ if(class(data)=='csAnno'){ pdata <- data@annoStat }else{ pdata <- data } pdata$ratio <- round(pdata$Frequency, 2) pdata$label <- paste0(pdata$Feature, ' (', pdata$ratio,'%)') col <- ChIPseeker:::getCols(nrow(pdata)) p <- ggplot(pdata,aes(x=1,ratio,fill=factor(label, levels=label))) + geom_bar(stat='identity', color="black",size=0.2, show.legend=T) + labs(fill="") + coord_polar("y", start=0) + scale_y_reverse() + theme_void() + scale_fill_manual(values=col) + theme(legend.spacing.y = unit(0.1, 'cm'),legend.title = element_blank(),legend.box.margin=margin(1,10,1,-12), legend.key.size=unit(0.15,'inches'),legend.text=element_text(size=8)) + guides(fill = guide_legend(byrow = T)) return(p) } p <- plotpie(peakanno) p

结果如下：

结束语

  ChIPseeker做为Y叔的作品，使用体验依然是顺滑流畅没话说，功能也是非常全面强大。网上有很多关于ChIPseeker的帖子，这里列举一个本人觉得写得比较详细的帖子可供参考：https://www.jianshu.com/p/c76e83e6fa57。工欲善其事必先利其器，好的工具可以让我们事半功倍，节省更多的时间去关注数据本身的意义。不过，没事的时候学习一下工具的原理及实现过程还是有好处的，这样用起来才会更加得心应手来去自如。好了，今天到此结束~~~

往期回顾

R语言书籍免费领
可视化：网络图
可视化：Wordcloud
可视化：Dumbbell Chart
可视化：Arc Diagrams