logo标签列表

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

rperformanceigraphsubgraph
10

我有一个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)
  }
5个回答
7
使用Rcpp包是加快您的代码速度超出基本R可能性的一种方法。考虑以下Rcpp实现完整操作。它以以下内容作为输入：
- `valid`: 所有位于足够大的组件中的节点的索引 - `el`，`deg`，`firstPos`: 图的边列表的表示法（节点`i`的邻居是 `el [firstPos [i]]`， `el [firstPos [i] +1]`，...， `el [firstPos [i] + deg [i] -1]`） - `size`: 要采样的子图大小 - `nrep`: 重复次数 - `weights`: 存储在`V(G)$weight`中的边权重 - `RWRNodeweight`: 存储在`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))
}

相比于原始代码和我们迄今为止看到的所有igraph解决方案，执行1000个复制的时间非常快（请注意，在大部分基准测试中，我仅对原始的josilber和random_network函数进行了一次复制的测试，而不是1000次，因为测试1000次将耗费极长的时间）：

  • 尺寸= 10：0.06秒（与我之前提出的josilber函数相比加速1200倍，并且比原始的random_network函数快4000倍）
  • 尺寸= 100：0.08秒（与我之前提出的josilber函数相比加速8700倍，并且比原始的random_network函数快162000倍）
  • 尺寸= 1000：0.13秒（与我之前提出的josilber函数相比加速32000倍，并且比原始的random_network函数快2040万倍）
  • 尺寸= 5000：0.32秒（与我之前提出的josilber函数相比加速68000倍，并且比原始的random_network函数快2.9亿倍）

总之，使用Rcpp计算每个大小从5到500的1000个复制品可能是可行的（此操作总共大约需要1分钟），而使用迄今为止提出的纯R代码计算这些大小的每个复制品的1000个副本将速度过慢。

非常感谢@josilber，您的“纯”R答案非常棒，但是这个Rcpp实现似乎非常快！如果图的名称是“字符”而不是数字，您能解释一下如何创建边缘列表吗？ - user2380782
@user2380782 完成了 - 诀窍是将 names=FALSE 传递给 get.edgelist。 - josliber
我有另一个函数想要使用Rcpp实现，你能帮我吗？ - user2380782
@user2380782 当然可以，但请提出一个单独的问题。 - josliber
感谢@josilber，我刚刚发布了一个标题为“基于节点权重构建图形的算法优化”的新问题。如果这个问题需要很多你的时间，请忘记它。 - user2380782
显示剩余2条评论
2
基本上，对于随机抽取图形的算法可以描述为初始化节点集作为随机选择的节点，然后迭代地添加当前集合的邻居，直到没有更多的邻居或者您拥有所需的子集大小。
这里昂贵而重复的操作是确定当前集合的邻居，您可以使用以下方式进行操作：
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
}

对于单次复制，此方法比问题中的代码单次复制速度要快得多：
- 当size=50时，此方法的单次复制仅需0.3秒，而问题中的代码需要3.8秒。 - 当size=100时，此方法的单次复制仅需0.6秒，而问题中的代码需要15.2秒。 - 当size=200时，此方法的单次复制需要1.5秒，而问题中的代码需要69.4秒。 - 当size=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倍。
1

我没有完整的答案，但是以下是一些需要考虑的事项以帮助加快速度（假设没有使用不同方法的更快的方法）。

  1. Remove from your graph any any nodes which are not part of a component as large as you are looking for. It will really depend on your network structure but it looks like your networks are genes so there are likely many genes with very low degree and you could get some speedups by removing them. Something like this code:

    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. Consider running this in parallel to take advantage of multiple cores. For example, using the parallel package and 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)
                                     })
感谢@av1的建议，不过实际上我的基因网络是最大的组件，但并行选项是我将尝试的内容，只是在等待是否有其他更高效的方法来完成它。 - user2380782
这似乎是一个有用且通用的问题，如果它不存在的话，或许可以将其作为一个功能请求添加到 iGraph C 库 中？ - Stefan Avey
1
你的代码存在一些问题（例如，你没有预分配向量等）。请看下面我编写的代码。虽然我只测试了大小为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 <- neighbors(G, as.numeric(seed[j-1]))时，你只获取了最近添加的节点的邻居，而不是你到目前为止选择的所有节点的邻居（这是 OP 所做的）。因此，你的代码将无限循环搜索 graph.data.frame(cbind(c(1, 2, 2, 2), c(2, 3, 4, 5)), directed=F) 的大小为 4 或 5 的子集（实际上，它也会在大小为 3 时崩溃，但可以通过一些错误检查来修复）。 - josliber
另一个问题是，即使在无法获得正确大小的图形的情况下，代码仍会在诱导子图上调用score_fun。请注意，OP的代码将丢弃这些不够大的子图。 - josliber
哦，我还没有测试过小于5个顶点的情况，因为那是OP提到的最低数量。我只在ER图上进行了测试。我没有遇到大小为5的图形问题，但并没有深入研究它。 - Chris Watson
我只是想说，如果你在我在评论中提供的大小为4或5的图上运行你的代码，那么由于你检查邻居的方式，你的代码将会无限循环。 - josliber
1
我认为在igraph中使用cliques函数会更加高效，因为clique是完全连接节点的子图。只需将min和max设置为您要搜索的子图大小，它将返回所有大小为5的clique。然后，您可以选择满足您需求的任何子集。不幸的是，在您生成的Erdos-Renyi图示例中，最大clique的大小通常小于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)
我认为当作者说“在每个子图中，所有节点都连接在一起”时，他们的意思是“节点形成单个连通分量”，而不是“节点形成一个团”，这是你在这里假设的。 - josliber
此外，对于大型图形计算所有5到500个团是非常计算密集的。 - user2380782
啊，我现在明白了。谢谢你提醒我。 - Ryan Haunfelder
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接