从间隔列表中模拟随机位置

Question

我正在尝试在 R 中开发一个函数来输出给定间隔列表中的随机位置。

我的间隔文件（14,600 行）是一个制表符分隔的 bed 文件（chromosome start end name），如下所示：

1      4953    16204   1
1      16284   16612   1
1      16805   17086   1
1      18561   18757   1
1      18758   19040   1
1      19120   19445   1

目前我的函数将在这些间隔内生成 N 个随机位置。

sim_dat <- bpSim(N=10)
head(sim_dat)

  seqnames    start      end width strand
1       1 22686939 22686939     1      *
2       1 14467770 14467770     1      *
3       2 10955472 10955472     1      *
4        X   823201   823201     1      *
5        6 10421738 10421738     1      *
6       17 21827745 21827745     1      *

library(GenomicRanges)
library(rtracklayer)

bpSim <- function(intervals="intervals.bed", N=100, write=F) {
  intFile <- import.bed(intervals)
  space <- sum(width(intFile))
  positions <- sample(c(1:space), N)
  cat("Simulating", N, "breakpoints", sep = " ", "\n")
  new_b <- GRanges(
    seqnames = as.character(rep(seqnames(intFile), width(intFile))),
    ranges = IRanges(start = unlist(mapply(seq, from = start(intFile), to = end(intFile))), width = 1)
  )
  bedOut <- new_b[positions]
  if (write) {
    export.bed(new_b[positions], "simulatedBPs.bed")
  }
  remove(new_b)
  return(data.frame(bedOut))
}

这个有效，但是由于我对 GenomicRanges 包不是特别熟悉，所以我宁愿把它拼凑在一起。我更希望能够使用基础 R 或来自 tidyverse 的包重写它，以便我可以调整它，例如，允许用户指定染色体。

它也需要很长时间 - 即使是 N=10:

system.time(sim_dat <- bpSim(N=10))
Simulating 10 breakpoints 
   user  system elapsed 
 10.689   3.267  13.970

最终，我试图模拟基因组中的随机位置，因此需要为每个 N 模拟数百次数据。

如果有任何建议，我将不胜感激：

减少运行时间
删除 GenomicRanges

此外 - 如果有人知道任何已经这样做的包，我宁愿使用现有包而不是重新发明轮子。

Answer 1

由于范围长度不同，我假设您希望这些随机选择的位置与线段的长度成比例。换句话说，基于范围内的实际碱基对，选择是统一的。否则你将 over-representing 小范围（较高标记密度）和 under-representing 大范围（较低标记密度）。

这是一个 data.table 解决方案，可以在我的机器上几乎立即创建一千个站点，并在大约 10 秒内创建一百万个随机站点。它随机抽样您想要的网站数量，首先通过抽样行（按每行的范围大小加权），然后在该范围内均匀抽样。

library(data.table)

nSites <- 1e4

bed <- data.table(chromosome=1, start=c(100,1050,3600,4000,9050), end=c(1000,3000,3700,8000,20000))

# calculate size of range
bed[, size := 1 + end-start]

# Randomly sample bed file rows, proportional to the length of each range
simulated.sites <- bed[sample(.N, size=nSites, replace=TRUE, prob=bed$size)]

# Randomly sample uniformly within each chosen range
simulated.sites[, position := sample(start:end, size=1), by=1:dim(simulated.sites)[1]]

# Remove extra columns and format as needed
simulated.sites[, start  := position]
simulated.sites[, end := position]
simulated.sites[, c("size", "position") := NULL]

以 table 开头，例如：

 chromosome start   end  size
          1   100  1000   901
          1  1050  3000  1951
          1  3600  3700   101
          1  4000  8000  4001
          1  9050 20000 10951

输出如下：

       chromosome start   end
    1:          1 10309 10309
    2:          1  4578  4578
    3:          1  1984  1984
    4:          1 14703 14703
    5:          1 10090 10090
   ---
 9996:          1  1601  1601
 9997:          1  5317  5317
 9998:          1 18918 18918
 9999:          1  1154  1154
10000:          1  7343  7343

从间隔列表中模拟随机位置

Simulate random positions from a list of intervals

r

simulation

bioinformatics

bioconductor

genomicranges