使用 igraph 对不同大小的子图进行采样

sampling subgraphs from different sizes using igraph

我有一个 igraph 对象 mygraph,它有大约 10,000 个节点和大约 145,000 个边,我需要从这个图创建一些不同大小的子图。 我需要的是根据确定的大小(从 5 个节点到 500 个节点)创建子图,其中所有节点都连接在每个子图中。我需要为每个尺寸创建约 1,000 个子图(即 5 号 1000 个子图,6 号 1000 个,依此类推),然后根据不同的节点属性为每个图计算一些值。 我有一些代码,但完成所有计算需要很长时间。我想使用 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)
  }

我没有完整的答案,但这里有一些可以帮助加快速度的注意事项(假设没有更快的方法使用不同的方法)。

  1. 从图表中删除不属于您要查找的大组件的任何节点。这实际上取决于您的网络结构,但看起来您的网络是基因,因此可能有许多基因的度数非常低,您可以通过删除它们来获得一些加速。像这样的代码:

    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)
    
  2. 并行考虑 运行 以利用多核。例如,使用 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)
                                         })
    

我认为在 igraph 中使用 cliques 函数会更有效,因为 clique 是完全连接的节点的子图。只需将 min 和 max 设置为等于您正在搜索的子图的大小,它将 return 所有大小为 5 的派系。您可以采用满足您需要的任何子集。不幸的是,对于您生成的示例 Erdos-Renyi 图,最大的集团通常小于 5,因此这不适用于该示例。但是,它应该适用于真实网络,它比 Erdos-Renyi 图表现出更多的聚类,就像您最有可能的那样。

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 GHz。

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 

基本上,您用于对图进行采样的算法可以描述为将节点集初始化为随机选择的节点,然后迭代地添加当前集的邻居,直到没有更多的邻居或您具有所需的子集大小。

此处昂贵的重复操作是确定当前集合的邻居,您可以使用以下方法执行此操作:

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 个 worker 的并行版本)。因此,对于 size=100 案例,我们通过代码改进获得了 20 倍的加速,并通过并行化获得了额外的 3 倍加速,总加速为 60 倍。

一种比 base R 中可能的方法进一步加速代码的方法是使用 Rcpp 包。考虑以下完整操作的 Rcpp 实现。它将以下内容作为输入:

  • valid: 一个足够大的组件中所有节点的索引
  • eldegfirstPos:表示图的边列表(节点i的邻居是el[firstPos[i]]el[firstPos[i]+1], ..., el[firstPos[i]+deg[i]-1]).
  • size: 要采样的子图大小
  • nrep:重复次数
  • weightsV(G)$weight
  • 中存储的边权重
  • RWRNodeweightV(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());
}
")

现在我们需要在 base 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))
}

与原始代码和我们迄今为止看到的所有 igraph 解决方案相比,执行 1000 次复制的时间快得惊人(请注意,对于大部分基准测试,我测试了原始 josilber并且 random_network 函数用于 1 个复制而不是 1000 个,因为测试 1000 个将花费非常长的时间):

  • Size=10:0.06 秒(比我之前提出的 josilber 函数提速 1200 倍,比原来的 random_network 函数提速 4000 倍)
  • Size=100:0.08 秒(比我之前提出的 josilber 函数提速 8700 倍,比原来的 random_network 函数提速 162000 倍)
  • Size=1000:0.13 秒(比我之前提出的 josilber 函数提速 32000 倍,比原始 random_network 函数提速 2040 万倍 )
  • Size=5000:0.32 秒(比我之前提出的 josilber 函数提速 68000 倍,比原来的 random_network 函数提速 2.9 亿倍 )

简而言之,Rcpp 可能使计算从 5 到 500 的每个大小的 1000 个重复成为可能(此操作总共可能需要大约 1 分钟),而计算每个大小的 1000 个重复将是非常慢的这些大小使用迄今为止提出的纯 R 代码。