提高慢 R 代码的效率

Question

我有一个 I/GRanges Views 对象

** 数据简化版，实际数据量巨大

Views on a 10000000-length Rle subject

 views:
      start      end   width
 [1]      1     1000    1000 [100 100 100 100 100 100 100 100 100 100 ...]
 [2]   1001     2000    1000 [190 190 190 190 190 190 190 190 190 190 ...]
 [3]   2001     3000    1000 [280 280 280 280 280 280 280 280 280 280 ...]
 [4]   3001     4000    1000 [370 370 370 370 370 370 370 370 370 370 ...]
 [5]   4001     5000    1000 [460 460 460 460 460 460 460 460 460 460 ...]
 ...    ...      ...     ... ...
 [9996] 995001  9996000 9001000 [89650 89650 89650 89650 89650 89650 ...]
 [9997] 996001  9997000 9001000 [89740 89740 89740 89740 89740 89740 ...]
 [9998] 997001  9998000 9001000 [89830 89830 89830 89830 89830 89830 ...]
 [9999] 998001  9999000 9001000 [89920 89920 89920 89920 89920 89920 ...]
[10000] 999001 10000000 9001000 [90010 90010 90010 90010 90010 90010 ...]

每个视图（线）的宽度为 1000，这意味着 1000 个数据点，每个数据点 100 个。 现在，我想将一组数据点分成 20 个 bin（在这种情况下，每个 bin 50 个），然后取平均值，因此输出将是一个包含 20 个数字的向量，每个数字都是该 bin 的平均值。

输出：

[1] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100

现在，在实际情况中，我有 20 多个这样的视图，每行的宽度不同，有些行 > 5K。我的代码工作正常，但速度很慢，对于我的数据，每行返回一个包含 20 个 bin 的向量，大约需要 1.5 秒，我有 > 30K 行，大约需要 12.5 小时。

我敢肯定，有一些方法可以加快这些计算，如果没有的话，我可以以某种方式使用集群的并行节点。你有什么建议。

生成数据的测试代码：

library('GenomicRanges')
# generating data frame 
df=data.frame(chrom=rep('Chr1',100000),start=seq(1,1000000,by=1000),end=seq(1000,10000000,by=1000),strand=rep("+",100000))

# making GRanges object
gr=GRanges(seqnames=as.vector(df[,1]),IRanges(start=df[,2],end=df[,3]),strand=df[,4])
# obtaining coverage using function coverage in the form of RLE object
gr.cov=coverage(gr)
# generating views for specific start and end
gr.views=Views(gr.cov[[1]],start=seq(1,1000000,by=1000),end=seq(1000,10000000,by=1000))
# putting in temp variable
d=gr.views

# this following code calculates the matrix (where each line is 20 points) for 10 lines
# reduce or increase the number in the outermost sapply loop to increase/decrease the lines to be calculated

sapply(1:10,function(j)
   sapply(1:20,
   function(i)as.numeric(
     format(
       mean(
         as(d[[j]][(
           seq(0,length(d[[j]]),floor(length(d[[j]])/20))+1)[i]:
             c((seq(0,length(d[[j]]),floor(length(d[[j]])/20)))[
               -length((seq(0,length(d[[j]]),floor(length(d[[j]])/20))))
               ],length(d[[j]]))[i+1]],
            "RangedData")$score),
       digits=2)
     )
   )
)

Answer 1

与其根据基因创建视图，为什么不根据要进行计算的窗口创建视图，然后使用viewSums 或viewMaxs 之类的东西来计算视图的统计信息？假设您有一个 GRanges 描述“基因”的开始和结束（成绩单？）

genes <- GRanges(seqnames, IRanges(geneStarts, geneEnds))

你可以计算窗口的开始和结束

n <- 50L
starts0 <- Map(function(...) floor(seq(...)), start(genes), end(genes),
               MoreArgs=list(length.out=n + 1L))
ends <- lapply(starts0, function(x) floor(x)[-1])
starts <- lapply(starts0, function(x) x[-length(x)])

然后创建您的视图

v <- Views(gr.cov[[1]], start=unlist(starts), end=unlist(ends))

（参见?RleViews“Views,RleList-method”）计算统计数据并按基因分割

split(viewMeans(v), rep(seq_along(genes), each=n))

在Bioconductor mailing list 上提问可能会带来许多巧妙的解决方案。

v 是“RleViews”类的一个实例； v[[1]] 是 Rle 的一个实例。您可以计算mean(v[[1]]) 作为viewMeans 的确认，或者更进一步，将v[[1]] 强制转换为一个普通的旧向量并计算mean(as.vector(v[[1]])))。 runValue(v[[1]])（与 v[[1]]@values 相同，但使用适当的访问器，而不是在后台查看）返回 Rle 中的值，例如，

> (x <- Rle(c(rep(100, 10), rep(200, 10))))
numeric-Rle of length 20 with 2 runs
  Lengths:  10  10
  Values : 100 200
> runValue(x)
[1] 100 200
> runLength(x)
[1] 10 10

显然是mean(runValue(x)) != mean(x)。