ChIPseeker的upsetplot是怎么写的

  1. 云栖社区>
  2. 博客>
  3. 正文

ChIPseeker的upsetplot是怎么写的

徐洲更 2019-05-23 21:42:25 浏览334
展开阅读全文

在距今2年前的一个日子里,我写了一篇关于upsetR的简单教程。两年过去了,这篇教程是百度搜索的第一位。

在这篇教程中,我定了一个计划

下一期写一篇Y叔的upsetplot是如何写的。

然而两年过去了,我说的那个教程都没有出现,今天有人提问upsetR的时候,我又看到了这个计划,羞愧难当,遂有了这篇教程。

通常而言,看一个函数的源代码的方式就是直接在R中输入这个函数名

> upsetplot
standardGeneric for "upsetplot" defined from package "enrichplot"

function (x, ...) 
standardGeneric("upsetplot")
<bytecode: 0x0000000047c188c8>
<environment: 0x0000000047c10318>
Methods may be defined for arguments: x
Use  showMethods("upsetplot")  for currently available ones.

但是对于upsetplot这个函数而言,并不奏效,因为它是一个泛型函数。所谓的泛型函数,指的是对于同一个函数而言,它会根据不同数据对象调取实际的函数。例如print函数,既可以用来输出字符串,还可以输出图片。

那么如何查看它的实际函数呢?我们可以根据提示,showMethods("upsetplot")去找它到底有哪些可用函数

> showMethods("upsetplot")
Function: upsetplot (package enrichplot)
x="csAnno"
x="enrichResult"

不难发现,它提供了两种数据类型的展示方法,一种是csAnno,另一种是enrichResult.对于csAnno,我们可以用包名:::不可见函数的形式查看到实际代码

ChIPseeker:::upsetplot.csAnno
function (x, sets = NULL, order.by = "freq", sets.bar.color = NULL, 
    vennpie = FALSE, vp = viewport(x = 0.6, y = 0.7, width = 0.6, 
        height = 0.6), ...) 
{
    y <- x@detailGenomicAnnotation
    y <- as.matrix(y)
    y[y] <- 1
    y <- as.data.frame(y)
    if (is.null(sets)) {
        sets <- c("distal_intergenic", "downstream", 
            "threeUTR", "fiveUTR", "Intron", 
            "Exon", "Promoter")
        if (vennpie && is.null(sets.bar.color)) {
            sets.bar.color <- c("#d95f0e", "#fee0d2", 
                "#98D277", "#6F9E4C", "#fc9272", 
                "#9ecae1", "#ffeda0")
        }
    }
    if (is.null(sets.bar.color)) {
        sets.bar.color <- "black"
    }
    if (vennpie) {
        plot.new()
        upset(y, sets = sets, sets.bar.color = sets.bar.color, 
            order.by = order.by, ...)
        pushViewport(vp)
        par(plt = gridPLT(), new = TRUE)
        vennpie(x)
        popViewport()
    }
    else {
        upset(y, sets = sets, sets.bar.color = sets.bar.color, 
            order.by = order.by, ...)
    }
}

为了能够解释代码,我们需要先运行测试数据,获取能够进行演示的数据

require(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
peakfile <- system.file("extdata", "sample_peaks.txt", package="ChIPseeker")
peakAnno <- annotatePeak(peakfile, tssRegion=c(-3000, 3000), TxDb=txdb)
x <- peakAnno # x是csAnno类

y <- x@detailGenomicAnnotation用来提取详细的注释信息

注释信息

y <- as.matrix(y)转换成矩阵,方便后面的y[y]<-1将所有为TRUE的结果转成1,在y<- as.data.frame(y)的时候让FALSE转成了0。

调整数据结构

如果是我的话,我应该会用apply

y <- x@detailGenomicAnnotation
y2 <- apply(y, 2, as.integer)
y2 <- as.data.frame(y2)

这里的y的数据结构刚好是的upset需要的输入,所以假如Y叔没有提供upsetplot函数的话,那么代码你可以直接用upset(y)作图

upset(y)

直接作图吧

下面代码定义需要提供给upset函数进行集合运算的属性,默认就是ChIPseeker注释的8个属性

设置集合

之后的代码是根据是否绘制饼图走向两条分支。

绘图

在绘制饼图的循环里,Y叔先用graphics的plot.new新建了一个绘图框架,方便在这个绘图框架中作图。

plot.new()

空空如也

再用upset把图画好

# 绘图的准备工作
sets <- c("distal_intergenic", "downstream", 
              "threeUTR", "fiveUTR", "Intron", 
              "Exon", "Promoter")
sets.bar.color <- c("#d95f0e", "#fee0d2", 
                          "#98D277", "#6F9E4C", "#fc9272", 
                          "#9ecae1", "#ffeda0")
sets.bar.color <- "black"
order.by <- "freq"
library(UpSetR)
upset(y, sets = sets, sets.bar.color = sets.bar.color, 
       order.by = order.by)

最初的结果

之后的pushViewport(vp)会调出一个视图对象,也就是vp = viewport(x = 0.6, y = 0.7, width = 0.6, height = 0.6), 指定了后续要画图的位置和大小.

par(plt = gridPLT(), new = TRUE)设置了绘图的基本参数。gridPLT函数来自于gridBase,他的作用就是返回用于设置plt的值。

最终用vennpie(x)绘制饼图,然后用popViewport()回到原来视图,结束画图。

网友评论

登录后评论
0/500
评论
徐洲更
+ 关注