从聚类和共现因素列表生成文氏图

Question

我有一个输入文件，其中包含大约50000个聚类以及每个聚类中存在的若干因素（总共约1000万条记录），以下是一个更小的示例:

set.seed(1)
x = paste("cluster-",sample(c(1:100),500,replace=TRUE),sep="")
y = c(
  paste("factor-",sample(c(letters[1:3]),300, replace=TRUE),sep=""),
  paste("factor-",sample(c(letters[1]),100, replace=TRUE),sep=""),
  paste("factor-",sample(c(letters[2]),50, replace=TRUE),sep=""),
  paste("factor-",sample(c(letters[3]),50, replace=TRUE),sep="")
)
data = data.frame(cluster=x,factor=y)

在另一个问题的帮助下，我设法生成了如下所示的因素共现饼图：

counts = with(data, table(tapply(factor, cluster, function(x) paste(as.character(sort(unique(x))), collapse='+'))))
pie(counts[counts>1])

但现在我想要一个显示因子共现的维恩图。最好可以设置每个因子的最小计数阈值。例如，对于不同因子的维恩图，每个因子在每个群集中都必须出现n>10次才能被考虑。

我尝试使用aggregate生成计数表，但无法使其工作。

- avilella

---

你看过任何R语言的Venn图包吗？可以看看G. Jay Kerns最近使用venneuler库的这个例子，或者是在《统计软件杂志》上使用venn库的简短文章(Murdoch, 2004)。如果这只是关于R编程，那么应该迁移到SO。 - Andy W

Avilella，这个问题可能没有得到任何答案，因为它略微偏离了主题。您可以在SO上找到更好的解决方法，那里有一个活跃的R用户社区。但请不要交叉发布：如果您希望将其迁移，请仅标记该问题以引起管理员的注意。 - whuber

我已经标记了它，但我还没有看到它被移动到SO上... - 719016

嗯...没有引起警觉。我来这里只是因为我记得要来。无论如何，我们开始吧。我相信你会在SO上得到一些好的回复。 - whuber

请关注此线程：https://dev59.com/HGoy5IYBdhLWcg3wMLOj#22988825 - ktyagi

Answer 1

我提供了两种解决方案，使用两个具有 Venn 图功能的不同软件包。正如您所预期的那样，这两种解决方案都涉及使用 aggregate() 函数的初始步骤。

我倾向于使用 venneuler 软件包生成的结果。它的默认标签位置不理想，但您可以通过查看相关的 plot 方法（可能使用 locator() 选择坐标）来进行调整。

第一种解决方案：

有一种可能性是使用 venneuler 软件包中的 venneuler() 画出你的 Venn 图。

library(venneuler)

## Modify the "factor" column, by renaming it and converting
## it to a character vector.
levels(data$factor) <- c("a", "b", "c")
data$factor <- as.character(data$factor)

## FUN is an anonymous function that determines which letters are present
## 2 or more times in the cluster and then pastes them together into 
## strings of a form that venneuler() expects.
##
inter <- aggregate(factor ~ cluster, data=data,
                   FUN = function(X) {
                       tab <- table(X)
                       names <- names(tab[tab>=2])
                       paste(sort(names), collapse="&")
                   })            
## Count how many clusters contain each combination of letters
counts <- table(inter$factor)
counts <- counts[names(counts)!=""]  # To remove groups with <2 of any letter
#  a   a&b a&b&c   a&c     b   b&c     c 
# 19    13    12    14    13     9    12 

## Convert to proportions for venneuler()
ps <- counts/sum(counts)

## Calculate the Venn diagram
vd <- venneuler(c(a=ps[["a"]], b = ps[["b"]], c = ps[["c"]],
                  "a&b" = ps[["a&b"]],
                  "a&c" = ps[["a&c"]],
                  "b&c" = ps[["b&c"]],
                  "a&b&c" = ps[["a&b&c"]]))
## Plot it!
plot(vd)

在编写这段代码时，我做了一些选择:

• 我将因子的名称从"factor-a"更改为"a"。你可以随意改回去。

• 我只需要每个因子在每个聚类中至少出现>=2次（而不是>10）即可被计算。 （这是为了使用您数据的这个小子集演示代码。）

• 如果您查看中间对象counts，您会发现它包含一个初始的未命名元素。该元素是包含小于2个任何字母的聚类数。您可以比我更好地决定是否要在计算随后的ps（“比例”）对象时包括它们。

解决方案第二种:

另一种可能性是利用Bioconductor包limma中的vennCounts()和vennDiagram()。要下载该软件包，请按照此处的说明进行操作。与上面的venneuler解决方案不同，结果图中的重叠不是与实际交集的比例成比例的。相反，它用实际频率注释图。 （请注意，此解决方案不涉及对data$factor列的任何编辑。）

library(limma)

out <- aggregate(factor ~ cluster, data=data, FUN=table)
out <- cbind(out[1], data.frame(out[2][[1]]))

counts <- vennCounts(out[, -1] >= 2)
vennDiagram(counts, names = c("Factor A", "Factor B", "Factor C"),
            cex = 1, counts.col = "red")

- Josh O'Brien