Выборка подграфов разных размеров с помощью igraph
У меня есть объект igraph mygraph с ~10000 узлов и ~145000 ребер, и мне нужно создать несколько подграфов из этого графа, но с разными размерами. Мне нужно создать подграфы определенного размера (от 5 до 500 узлов), где все узлы соединены в каждом подграфе. Мне нужно создать ~1000 подграфов для каждого размера (т. Е. 1000 подграфов для размера 5, 1000 для размера 6 и т. Д.), А затем вычислить некоторые значения для каждого графа в соответствии с различными атрибутами узла. У меня есть некоторый код, но все расчеты занимают много времени. Я думал, используя graphlets Функция для получения разных размеров, но каждый раз, когда я запускаю ее на своем компьютере, происходит сбой из-за проблем с памятью.
Вот код, который я использую:
Первым шагом было создание функции для создания подграфов разных размеров и выполнения необходимых расчетов.
random_network<-function(size,G){
score_fun<-function(g){
subsum <- sum(V(g)$weight*V(g)$RWRNodeweight)/sqrt(sum(V(g)$RWRNodeweight^2))
subsum
}
genes.idx <- V(G)$name
perm <- c()
while(length(perm)<1000){
seed<-sample(genes.idx,1)
while( length(seed)<size ){
tmp.neigh <- V(G)[unlist(neighborhood(G,1,seed))]$name
tmp.neigh <- setdiff(tmp.neigh, seed)
if( length(tmp.neigh)>0 )
seed<-c(seed,sample(tmp.neigh,1)) else break
}
if( length(seed)==size )
perm <- c(perm,score_fun(induced.subgraph(G,seed)))
}
perm
}
Вторым шагом было применение функции к реальному графику.
### generate some example data
library(igraph)
my_graph <- erdos.renyi.game(10000, 0.0003)
V(my_graph)$name <- 1:vcount(my_graph)
V(my_graph)$weight <- rnorm(10000)
V(my_graph)$RWRNodeweight <- runif(10000, min=0, max=0.05)
### Run the code to get the subgraphs from different size and do calculations based on nodes
genesets.length<- seq(5:500)
genesets.length.null.dis <- list()
for(k in 5:max(genesets.length){
genesets.length.null.dis[[as.character(k)]] <- random_network(size=k,G=my_graph)
}
5 ответов
Один из способов ускорить ваш код дальше, чем это возможно в base R, - использовать пакет Rcpp. Рассмотрим следующую реализацию Rcpp полной операции. Он принимает в качестве входных данных следующее:
valid: Индексы всех узлов, которые находятся в достаточно большом компонентеel,deg,firstPos: Представление списка ребер графа (узелiсоседиel[firstPos[i]],el[firstPos[i]+1],...,el[firstPos[i]+deg[i]-1]).size: Размер подграфа для выборкиnrep: Количество повторенийweights: Вес ребер хранится вV(G)$weightRWRNodeweight: Вес ребер хранится вV(G)$RWRNodeweight
library(Rcpp)
cppFunction(
"NumericVector scores(IntegerVector valid, IntegerVector el, IntegerVector deg,
IntegerVector firstPos, const int size, const int nrep,
NumericVector weights, NumericVector RWRNodeweight) {
const int n = deg.size();
std::vector<bool> used(n, false);
std::vector<bool> neigh(n, false);
std::vector<int> neighList;
std::vector<double> scores(nrep);
for (int outerIter=0; outerIter < nrep; ++outerIter) {
// Initialize variables
std::fill(used.begin(), used.end(), false);
std::fill(neigh.begin(), neigh.end(), false);
neighList.clear();
// Random first node
int recent = valid[rand() % valid.size()];
used[recent] = true;
double wrSum = weights[recent] * RWRNodeweight[recent];
double rrSum = RWRNodeweight[recent] * RWRNodeweight[recent];
// Each additional node
for (int idx=1; idx < size; ++idx) {
// Add neighbors of recent
for (int p=firstPos[recent]; p < firstPos[recent] + deg[recent]; ++p) {
if (!neigh[el[p]] && !used[el[p]]) {
neigh[el[p]] = true;
neighList.push_back(el[p]);
}
}
// Compute new node to add from all neighbors
int newPos = rand() % neighList.size();
recent = neighList[newPos];
used[recent] = true;
wrSum += weights[recent] * RWRNodeweight[recent];
rrSum += RWRNodeweight[recent] * RWRNodeweight[recent];
// Remove from neighList
neighList[newPos] = neighList[neighList.size() - 1];
neighList.pop_back();
}
// Compute score from wrSum and rrSum
scores[outerIter] = wrSum / sqrt(rrSum);
}
return NumericVector(scores.begin(), scores.end());
}
")
Теперь все, что нам нужно сделать в базе R, это сгенерировать аргументы для scores, что можно сделать довольно легко:
josilber.rcpp <- function(size, num.rep, G) {
n <- length(V(G)$name)
# Determine which nodes fall in sufficiently large connected components
comp <- components(G)
valid <- which(comp$csize[comp$membership] >= size)
# Construct an edge list representation for use in the Rcpp code
el <- get.edgelist(G, names=FALSE) - 1
el <- rbind(el, el[,2:1])
el <- el[order(el[,1]),]
deg <- degree(G)
first.pos <- c(0, cumsum(head(deg, -1)))
# Run the proper number of replications
scores(valid-1, el[,2], deg, first.pos, size, num.rep,
as.numeric(V(G)$weight), as.numeric(V(G)$RWRNodeweight))
}
Время выполнения 1000 копий молниеносно по сравнению с исходным кодом и всеми igraph решения, которые мы видели до сих пор (обратите внимание, что для большей части этого бенчмаркинга я проверил оригинал josilber а также random_network функции для 1 репликации вместо 1000, потому что тестирование на 1000 заняло бы слишком много времени):
- Размер =10: 0,06 секунды (ускорение в 1200 раз по сравнению с предложенным ранее
josilberфункция и 4000-кратное ускорение по сравнению с оригиналомrandom_networkфункция) - Размер =100: 0,08 секунды (ускорение в 8700 раз по сравнению с ранее предложенным
josilberфункция и ускорение 162000x по сравнению с оригиналомrandom_networkфункция) - Размер =1000: 0,13 секунды (ускорение в 32000 раз по сравнению с ранее предложенным
josilberфункция и ускорение в 20,4 миллиона раз по сравнению с оригиналомrandom_networkфункция) - Размер =5000: 0,32 секунды (ускорение в 68000 раз по сравнению с предложенным ранее)
josilberфункция и ускорение в 290 миллионов раз по сравнению с оригиналомrandom_networkфункция)
Короче говоря, Rcpp, вероятно, делает возможным вычисление 1000 повторов для каждого размера от 5 до 500 (всего эта операция, вероятно, займет около 1 минуты), в то время как вычисление 1000 повторений для каждого из этих размеров будет чрезмерно медленным с использованием чистый код R, который был предложен до сих пор.
По сути, ваш алгоритм выборки графа можно описать как инициализацию набора узлов в качестве случайно выбранного узла, а затем итеративное добавление соседа из вашего текущего набора до тех пор, пока не останется больше соседей или у вас не будет желаемого размера подмножества.
Дорогая повторная операция здесь определяет соседей текущего набора, что вы делаете со следующим:
tmp.neigh <- V(G)[unlist(neighborhood(G,1,seed))]$name
tmp.neigh <- setdiff(tmp.neigh, seed)
Короче говоря, вы смотрите на соседей каждого узла в выбранном вами подмножестве на каждой итерации. Более эффективный подход - хранить вектор соседей и обновлять его каждый раз, когда вы добавляете новый узел; это будет более эффективным, потому что вам нужно только рассмотреть соседей нового узла.
josilber <- function(size, num.rep, G) {
score_fun <- function(vert) sum(vert$weight*vert$RWRNodeweight)/sqrt(sum(vert$RWRNodeweight^2))
n <- length(V(G)$name)
# Determine which nodes fall in sufficiently large connected components
comp <- components(G)
valid <- which(comp$csize[comp$membership] >= size)
perm <- replicate(num.rep, {
first.node <- sample(valid, 1)
used <- (1:n) == first.node # Is this node selected?
neigh <- (1:n) %in% neighbors(G, first.node) # Does each node neighbor our selections?
for (iter in 2:size) {
new.node <- sample(which(neigh & !used), 1)
used[new.node] <- TRUE
neigh[neighbors(G, new.node)] <- TRUE
}
score_fun(V(G)[used])
})
perm
}
Для одного экземпляра это приводит к значительному ускорению по сравнению с одним повторением кода в вопросе:
- Для размера =50 один экземпляр занимает 0,3 секунды для этого кода и 3,8 секунды для отправленного кода
- Для размера =100 один экземпляр занимает 0,6 секунды для этого кода и 15,2 секунды для отправленного кода
- Для размера =200 один экземпляр занимает 1,5 секунды для этого кода и 69,4 секунды для отправленного кода
- Для размера =500 одна реплика для этого кода занимает 2,7 секунды (таким образом, 1000 копий должны занимать около 45 минут); Я не тестировал ни одной копии опубликованного кода.
Как упоминалось в других ответах, распараллеливание может еще больше повысить производительность выполнения 1000 повторов для заданного размера графа. Следующее использует doParallel пакет для распараллеливания повторяющегося шага (корректировка практически идентична корректировке, сделанной @Chris в его ответе):
library(doParallel)
cl <- makeCluster(4)
registerDoParallel(cl)
josilber2 <- function(size, num.rep, G) {
score_fun <- function(vert) sum(vert$weight*vert$RWRNodeweight)/sqrt(sum(vert$RWRNodeweight^2))
n <- length(V(G)$name)
# Determine which nodes fall in sufficiently large connected components
comp <- components(G)
valid <- which(comp$csize[comp$membership] >= size)
perm <- foreach (i=1:num.rep, .combine='c') %dopar% {
library(igraph)
first.node <- sample(valid, 1)
used <- (1:n) == first.node # Is this node selected?
neigh <- (1:n) %in% neighbors(G, first.node) # Does each node neighbor our selections?
for (iter in 2:size) {
new.node <- sample(which(neigh & !used), 1)
used[new.node] <- TRUE
neigh[neighbors(G, new.node)] <- TRUE
}
score_fun(V(G)[used])
}
perm
}
На моем MacBook Air, josilber(100, 1000, my_graph) занимает 670 секунд (это не параллельная версия), в то время как josilber2(100, 1000, my_graph) запуск занимает 239 секунд (это параллельная версия с 4 рабочими). Для size=100 В этом случае мы получили 20-кратное ускорение от улучшений кода и дополнительное 3-кратное ускорение от распараллеливания для общего ускорения в 60 раз.
Я думаю, что было бы гораздо эффективнее использовать функцию igraph в igraph, так как клика - это подграф полностью связанных узлов. Просто установите min и max равными размеру подграфа, который вы ищете, и он вернет все клики размера 5. Вы можете взять любое подмножество из них, которое соответствует вашим потребностям. К сожалению, на примере графа Эрдоша-Рени, который вы генерировали часто, наибольшая клика меньше 5, поэтому в этом примере это не сработает. Тем не менее, он должен работать просто отлично для реальной сети, которая демонстрирует больше кластеризации, чем график Эрдоша-Рени, как, скорее всего, у вас.
library(igraph)
##Should be 0.003, not 0.0003 (145000/choose(10000,2))
my_graph <- erdos.renyi.game(10000, 0.003)
cliques(my_graph,min=5,max=5)
У вас есть ряд проблем с вашим кодом (вы не выделяете векторы и т. Д.). Пожалуйста, посмотрите код, который я придумал ниже. Я только проверил это до подграфа размера 100, все же. Тем не менее, при увеличении размера подграфа, экономия скорости возрастает по сравнению с вашим кодом. Вы должны установить foreach пакет, а также. Я запускал это на ноутбуке с 4 ядрами, 2,1 ГГц.
random_network_new <- function(gsize, G) {
score_fun <- function(g) {
subsum <- sum(V(g)$weight * V(g)$RWRNodeweight) / sqrt(sum(V(g)$RWRNodeweight^2))
}
genes.idx <- V(G)$name
perm <- foreach (i=seq_len(1e3), .combine='c') %dopar% {
seed <- rep(0, length=gsize)
seed[1] <- sample(genes.idx, 1)
for (j in 2:gsize) {
tmp.neigh <- neighbors(G, as.numeric(seed[j-1]))
tmp.neigh <- setdiff(tmp.neigh, seed)
if (length(tmp.neigh) > 0) {
seed[j] <- sample(tmp.neigh, 1)
} else {
break
}
}
score_fun(induced.subgraph(G, seed))
}
perm
}
Обратите внимание, что я переименовал функцию в random_network_new и аргумент gsize,
system.time(genesets <- random_network_new(gsize=100, G=my_graph))
user system elapsed
1011.157 2.974 360.925
system.time(genesets <- random_network_new(gsize=50, G=my_graph))
user system elapsed
822.087 3.119 180.358
system.time(genesets <- random_network_new(gsize=25, G=my_graph))
user system elapsed
379.423 1.130 74.596
system.time(genesets <- random_network_new(gsize=10, G=my_graph))
user system elapsed
144.458 0.677 26.508
Один пример использования вашего кода (у меня более чем в 10 раз быстрее для размера подграфа 10; с большими подграфами это было бы намного быстрее):
system.time(genesets_slow <- random_network(10, my_graph))
user system elapsed
350.112 0.038 350.492
У меня нет полного ответа, но вот некоторые вещи, которые следует учитывать, чтобы ускорить его (при условии, что не существует более быстрого подхода с использованием другого метода).
Удалите из своего графика любые узлы, которые не являются частью такого большого размера, как вы ищете. Это действительно будет зависеть от структуры вашей сети, но похоже, что ваши сети - это гены, поэтому, вероятно, существует много генов с очень низкой степенью, и вы можете получить некоторые ускорения, удалив их. Что-то вроде этого кода:
cgraph <- clusters(G) tooSmall <- which(cgraph$csize < size) toKeep <- setdiff(1:length(V(G)), which(cgraph$membership %in% tooSmall)) graph <- induced.subgraph(G, vids=toKeep)Попробуйте запустить это параллельно, чтобы использовать преимущества нескольких ядер. Например, используя
parallelпакет иmclapply,library(parallel) genesets.length<- seq(5, 500) names(genesets.length) <- genesets.length genesets.length.null.dis <- mclapply(genesets.length, mc.cores=7, function(length) { random_network(size=length, G=my_graph) })