Диаграмма Венна из списка кластеров и сопутствующих факторов
У меня есть входной файл со списком ~50000 кластеров и наличием ряда факторов в каждом из них (всего ~10 миллионов записей), см. Меньший пример ниже:
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 в каждом кластере, который должен быть принят во внимание.
Я пытался найти способ произвести подсчет таблиц с помощью агрегата, но не смог заставить его работать.
1 ответ
Я предоставил два решения, используя два разных пакета с возможностями диаграммы Венна. Как вы и ожидали, оба включают начальный шаг с использованием aggregate()
функция.
Я склонен предпочесть результаты из venneuler
пакет. Это положение меток по умолчанию не идеально, но вы можете изменить их, посмотрев на plot
метод (возможно, с использованием locator()
выбрать координаты).
Решение 1-е:
Одна из возможностей заключается в использовании venneuler()
в venneuler
пакет, чтобы нарисовать вашу диаграмму Венна.
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
("пропорции") объект.
Решение 2-е:
Другая возможность состоит в том, чтобы нанять vennCounts()
а также vennDiagram()
в пакете Биокондуктор limma
, Чтобы скачать пакет, следуйте инструкциям здесь. в отличие от 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")