使用igraph从不同尺寸中采样子图

Question

我有一个mygraph有10,000个节点和~145,000个边的igraph对象,我需要从这个图中创建一些子图,但是它们的大小不同.我需要的是从确定的大小(从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)
  }

Answer 1

一种加速代码的方法比使用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)$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))
}

执行1000次重复的时间超快相比原来的代码和所有igraph到目前为止,我们已经看到了解决方案(注意,这多少标杆我测试的原始josilber和random_network功能1个复制,而不是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分钟),而使用这些大小计算每个大小的1000个重复将是非常慢的.到目前为止已提出的纯R代码.