ChIPseeker的upsetplot是怎么写的

2023-10-13 15:33:22

在距今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()回到原来视图，结束画图。

码农公寓

相关文章