r-统计数据帧中给定的基因组区域周围的出现次数

Question

我需要计算基因组中发生在某些点或范围内的突变。突变具有基因组位置（染色体和碱基对，例如Chr110658324）。范围或斑点分别定义为基因组中给定位置的上下游（+-）10000个碱基对。突变的位置和"斑点"的位置都存储在数据帧中。

示例：

set.seed(1)
Chr <- 1
Pos <- as.integer(runif(5000 , 0, 1e8))
mutations <- data.frame(Pos, Chr)
Chr <- 1
Pos <- as.integer(runif(50 , 0, 1e8))
spots <- data.frame(Pos, Chr)

所以我要问的问题是：在"斑点"中给出的位置周围有多少个+/-10k碱基对的突变。（例如，如果光斑为100k，则范围为90k-110k）真实的数据当然包含所有24条染色体，但为了简单起见，我们现在可以专注于一条染色体。最终数据应该包含"点"及其附近的突变数量，最好是在数据帧或矩阵中。

非常感谢您的任何建议或帮助！

---

这是第一次尝试，但我很高兴有一种更优雅的方法。

w <- 10000  #setting range to 10k basepairs
loop <- spots$Pos  #creating vector of positions to loop through
out <- data.frame(0,0)
colnames(out) <- c("Pos", "Count")
for (l in loop) {
  temp <- nrow(filter(mutations, Pos>=l-w, Pos<=l+w))
  temp2 <- cbind(l,temp)
  colnames(temp2) <- c("Pos", "Count")
  out <- rbind(out, temp2)
}
 out <- out[-1,]

Answer 1

使用data.table中心点，然后聚合：

library(data.table)
#set the flank
myFlank <- 100000
#convert to ranges with flank
spotsRange <- data.table(
  chr = spots$Chr,
  start = spots$Pos - myFlank,
  end = spots$Pos + myFlank,
  posSpot = spots$Pos,
  key = c("chr", "start", "end"))
#convert to ranges start end same as pos
mutationsRange <- data.table(
  chr = mutations$Chr,
  start = mutations$Pos,
  end = mutations$Pos,
  key = c("chr", "start", "end"))
#merge by overlap
res <- foverlaps(mutationsRange, spotsRange, nomatch = 0)
#count mutations
resCnt <- data.frame(table(res$posSpot))
colnames(resCnt) <- c("Pos", "MutationCount")
merge(spots, resCnt, by = "Pos")
#         Pos Chr MutationCount
# 1   3439618   1            10
# 2   3549952   1            15
# 3   4375314   1            11
# 4   7337370   1            13
# ...

Answer 2

我不熟悉R中的床操作，所以我将用bedtools提出一个答案，这里的人可以尝试转换到GRanges或其他R生物信息学库。

本质上，你有两个床档案，一个是你的斑点，另一个是突变（我假设后者中每个都有1bp的坐标）。在这种情况下，你应该使用closestBed来获得最近的斑点和每个突变的距离（以bp为单位），然后过滤那些距离斑点10KB的斑点。UNIX环境中的代码如下所示：

# Assuming 4-column file structure (chr start end name)
closestBed -d -a mutations.bed -b spots.bed | awk '$9 <= 10000 {print}'

其中，第9列（$9）将是距离最近点的距离（以bp为单位）。根据您想要的具体程度，您可以在http://bedtools.readthedocs.io/en/latest/content/tools/closest.html.我确信R中至少有一个类似bedtools的包。如果功能相似，您可以应用完全相同的解决方案。

希望能有所帮助！