块 bootstrap 用于基因组数据

Question

我正在尝试实施块 bootstrap 过程，但我还没有想出有效的方法。

我的 data.frame 具有以下结构：

CHR POS var_A var_B
1 192 0.9 0.7
1 2000  0.8 0.3
2 3 0.21  0.76 
2 30009 0.36  0.15
...

第一列是染色体标识，第二列是位置，最后两列是我要计算相关性的变量。问题是每一行并不完全相互独立，这取决于它们之间的距离（越近越依赖），所以我不能简单地做 cor(df$var_A, df$var_B).

通常用于此类数据的解决此问题的方法是执行块 bootstrap。也就是说，我需要将我的数据分成长度为 X 的块，在该块内随机 select 一行，然后计算我感兴趣的统计数据。但是请注意，这些块需要根据列 POS 定义，而不是根据行号定义。此外，需要对每个染色体执行此过程。

我试图实现它，但我想出了最慢的代码（它甚至没有完成运行）而且我不能 100% 确定它是否有效。

x = 1000
cors = numeric()
iter = 1000
for(j in 1:iter) {
  df=freq[0,]
  for (i in unique(freq$CHR)) {
    t = freq[freq$CHR==i,]
    fim = t[nrow(t),2]
    i = t[1,2]
    f = i + x
    while(f < fim) {
      rows = which(t$POS>=i & t$POS<f)
      s = sample(rows)
      df = rbind(df,t[s,])
      i = f
      f = f + x
    }
  }
  cors = c(cors, cor(df$var_A, df$var_B))
}

有人能帮帮我吗？我相信有一种更有效的方法可以做到这一点。

提前致谢。

Answer 1

希望我没听错：

# needed for round_any()
library(plyr)

res <- lapply(unique(freq$CHR),function(x){
  
  freq_sel <- freq[freq$CHR==x,]
  blocks <- lapply(seq(1,round_any(max(freq_sel$POS),1000,ceiling),1000), function(ix) freq_sel[freq_sel$POS > ix & freq_sel$POS <= ix+999,])
  do.call(rbind,lapply(blocks,function(x) if (nrow(x) > 1) x[sample(1:nrow(x),1),] else x))
  
})

这应该 return 一个包含每条染色体条目的列表。在每个条目中，如果存在的话，每 1kb 块都有一个观察值。块数由最大POS值决定。

编辑：

library(doParallel)
library(foreach)
library(plyr)

cl <-  makeCluster(detectCores())
registerDoParallel(cl)


res <- foreach(x=unique(freq$CHR),.packages = 'plyr') %dopar% {
  
  freq_sel <- freq[freq$CHR==x,]
  blocks <- lapply(seq(1,round_any(max(freq_sel$POS),1000,ceiling),1000), function(ix) freq_sel[freq_sel$POS > ix & freq_sel$POS <= ix+999,])
  do.call(rbind,lapply(blocks,function(x) if (nrow(x) > 1) x[sample(1:nrow(x),1),] else x))
  
}

stopCluster(cl)

这是每个染色体上与 foreach 的简单并行化。重构函数并将并行处理基于另一个级别（例如 1000 次迭代或块）可能会更好。在任何情况下，我都可以再次强调我在评论中所说的话：在您开始并行化代码之前，您应该确保它尽可能高效。这意味着您可能需要查看 boot 包或类似包以提高效率。也就是说，根据您计划的迭代次数，一旦您对函数感到满意，并行处理可能会有用。

Answer 2

一种有效的尝试方法是使用 'boot' 包，其中的功能包括并行处理能力。

特别是 'tsboot' 或时间序列引导函数，将 select 有序数据块。如果您的 POS 变量是某种有序观察，这可能会起作用。

引导包功能很棒，但首先需要一点帮助。要在启动包中使用 bootstrap 函数，必须首先将感兴趣的统计信息包装在一个包含 index 参数的函数中。这是 bootstrap 生成的索引将用于将采样数据传递给您的统计信息的设备。

cor_hat <- function(data, index) cor(y = data[index,]$var_A, x = data[index,]$var_B)

在下面的参数中注意 cor_hat。 sim = "fixed", l = 1000 参数，表示您需要 fixed 个长度块 (l) 1000。但是，如果您试图捕捉随时间移动的最近邻动态，您可以使用任何大小的块，5 个或 10 个。 multicore 参数不言而喻，但如果您使用 windows，它可能 "snow"。

library(boot)
tsboot(data, cor_hat, R = 1000, sim = "fixed", l = 1000, parallel = "multicore", ncpus = 4)

此外，Elements of Statistical Learning 的第 194 页提供了使用传统 boot 函数的框架的一个很好的示例，所有这些都与 tsboot.

希望对你有所帮助，祝你好运。

贾斯汀

r

Answer 3

所以，过了一会儿，我想出了一个问题的答案。开始了。

您需要 dplyr.

包

l = 1000
teste = freq %>%
  mutate(w = ceiling(POS/l)) %>%
  group_by(CHR, w) %>%
  sample_n(1)

此代码根据基因组 (POS) 中的位置创建一个名为 w 的新变量。这个变量 w 是每一行分配给的 window，它取决于 l，也就是你的 window.

的长度

您可以多次重复此代码，每次对每个 window/CHR（使用 sample_n(1)）采样一行并应用您想要的任何感兴趣的统计数据。

块 bootstrap 用于基因组数据

Block bootstrap for genomic data

r

bioinformatics

genome

编辑：