基因组数据的块引导程序答案

【问题标题】：Block bootstrap for genomic data基因组数据的块引导程序
【发布时间】：2017-06-24 15:00:22
【问题描述】：

我正在尝试实现一个块引导程序，但我还没有找到有效的方法。

我的 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)。

解决此类数据常用的问题的方法是执行块引导。也就是说，我需要将我的数据分成长度为 X 的块，在该块内随机选择一行，然后计算我感兴趣的统计量。但是请注意，这些块需要基于列 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))
}

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

提前谢谢你。

【问题讨论】：

如何在POS上定义块？
因为我需要根据基因组中的位置每 1kb 块采样一行。在某些情况下，我在这些 1kb 块内有不止一行，但在某些情况下不会发生这种情况。
你的代码有点混乱。基于上面的小数据示例，基于POS 的条件是什么才能选择一行？
你必须沿着这些位置“走”，每 1kb 块采样一行。如果您的 POS 从 0 开始，那么您必须查找 [0,1000[ 范围内的行并采样一行。如果此块内没有行，则继续。如果只有一行，请保留。

标签： r bioinformatics genome

【解决方案1】：

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

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

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

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，它可能会“下雪”。

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

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

希望有帮助，祝你好运。

贾斯汀

【讨论】：

感谢您的回复。我只有几个问题：(i) cor_hat 函数不应该如下：cor_hat <- function(data, index) cor(y = data$var_A, x = data$var_B)?。 (ii) 我不必告诉 tsboot 列 POS 应该用于块？ (iii) 不同的染色体呢？如果我只考虑 POS，那么不同染色体中的东西就会混合在一起。
另外，我可以将我的 data.frame 传递给 tsboot 函数还是必须将其更改为时间序列对象？
听起来这不是一个时间序列，观察结果也不是完全有序的，而且我对该领域还不够熟悉，不知道 POS 和染色体是如何相关的。对于问题 (i)，您需要索引，它是原始包装函数的点。抱歉，我的语法来自另一个模型，我在上面和这里更正了 cor_hat <- function(data, index) cor(y = data[index,]$var_A, x = data[index,]$var_B) (ii) 如果它是一个时间序列并沿 POS 索引，tsboot 将自动生成它。但情况似乎并非如此。
我认为时间序列和基因组数据之间存在相似之处。 POS是基因组中观察的位置。基因组被分成相互独立的染色体。也许我们应该将每条染色体视为不同的时间序列？如何将我的数据转换为通过 POS 对其进行索引的时间序列？
它们是有序的，但也有一个结构，因为在每条染色体内部，POS 反映了观察的顺序。

【解决方案2】：

希望我的理解是正确的：

# 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))
  
})

这应该返回一个包含每个染色体条目的列表。在每个条目中，每 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 包或类似包以提高效率。也就是说，根据您计划的迭代次数，一旦您对自己的函数感到满意，并行处理可能会很有用。

【讨论】：

有效！谢谢你。我只是担心，因为它需要很长时间才能运行，并且要进行引导，我必须运行它可能 1k 次或更多。但无论如何，很高兴看到你的代码。我从中学到了很多。
你认为我们可以做些什么来加快代码速度吗？
如果我像这样运行代码，需要 52 小时才能完成。
你能分享一些关于你的数据大小/结构的信息吗？你需要迭代多少次？
如果您需要从 1 到 3000 万次遍历所有 1kb 块 1000 次，那将是很多迭代，我猜需要一些时间。我不确定我的方法是否最有效。使用简单的并行化，我可以在大约 1 分钟内完成一次迭代，如果你需要做 1000 次，这仍然是很多时间。我不熟悉boot 包，但也许有更有效的使用方法。你见过this question吗？也许这对你来说可能很有趣

【解决方案3】：

所以，过了一会儿，我想出了一个答案来解决我的问题。就这样吧。

你需要 dplyr 包。

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

此代码根据基因组中的位置 (POS) 创建一个名为 w 的新变量。这个变量w是每行被分配到的窗口，它取决于l，这是你的窗口的长度。

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

【讨论】：