r - 使用igraph采样不同大小的子图

Question

我有一个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)
  }

Answer 1

使用Rcpp软件包可以比在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)$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次复制的时间非常快(请注意，对于大多数基准测试，我测试了1次复制的原始 josilber和 random_network函数，而不是1000次复制，因为测试了1000将花费很长的时间):

大小= 10:0.06秒(比我先前建议的josilber函数高1200倍，比原始random_network函数高4000倍)

Size = 100:0.08秒(比我先前建议的josilber函数高8700倍，比原始random_network函数高162000x)

大小= 1000:0.13秒(比我以前建议的josilber函数高32000倍，而比原始random_network函数高2040万倍)

大小= 5000:0.32秒(比我以前建议的josilber函数高68000倍，而比原始random_network函数高2.9亿倍)

简而言之，Rcpp可能使计算从5到500的每个大小的1000个重复项变得可行(此操作可能总共花费大约1分钟)，而使用以下命令为每个大小计算1000个重复项将太慢了。到目前为止已经提出的纯R代码。