ChIPseeker绘图函数借用

Zad • • 杂文

日常瞎掰

  对于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还有其他参数可以根据需要设置，这里就不展示了。

  1. 使用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
                   |        
  [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
             
  [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
            
  [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

结果如下：

  1. 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。

  1. 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')

结果如下：

  1. 自定义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

版权声明：
作者：Zad
链接：https://www.techfm.club/p/42898.html
来源：TechFM
文章版权归作者所有，未经允许请勿转载。

THE END
分享
二维码
< <上一篇
下一篇>>