在 R 中按范围合并 - 应用循环答案

【问题标题】：Merge by Range in R - Applying Loops在 R 中按范围合并 - 应用循环
【发布时间】：2012-08-07 04:47:17
【问题描述】：

我在这里发布了一个问题：Matched Range Merge in R 关于根据一个文件中的数字合并两个文件，该数字落入第二个文件的范围内。到目前为止，我一直未能成功拼凑代码来完成此任务。我遇到的问题是我使用的代码逐行比较文件。这是一个问题，因为 1.) 一个文件比另一个文件长得多，并且 2.) 我需要通过较长文件中的每个范围对扫描较短文件中的行 - 而不仅仅是同一行中的范围.

我一直在使用原始问题中发布的函数，我觉得应该有一种方法可以将其应用于更通用的循环，将第一个文件中的每一行与第二个文件中的每一行进行比较，但是我还没弄明白。如果有人有任何建议，我将不胜感激。

**** 已编辑。

数据的性质是这样的：每个范围不一定是唯一的，尽管大多数是唯一的。它们的大小也不相同，有些完全属于其他范围。 findInterval 因此会产生错误，因为无法对范围进行排序以符合“非降序”顺序。

这里是每个数据框的前 6 行：

file1test <- data.frame(SNP=c("rs2343", "rs211", "rs754", "rs854", "rs343", "rs626"), BP=c(860269, 369640, 861822, 367934, 706940, 717244))


file2 <- data.frame(Gene=c("E613", "E92", "E49", "E3543", "E11", "E233"), BP_start=c(367640, 621059, 721320, 860260, 861322, 879584), BP_end = c(368634, 622053, 722513, 879955, 879533, 894689))

所以，如您所见，第 5 行的范围在第 4 行的范围内，第一个文件中的两个 SNP 在第 4 行的范围内，但只有一个在第 4 行的范围内第二行。

包含 SNP 的第一个文件只有约 400 行。然而，包含范围的第二个文件大约有 20K。我想作为输出生成一个数据框，其中包含来自第一个文件（SNP）的行，其中 BP 属于第二个文件中的 BP 范围。如果一个 SNP 属于两个范围，那么它会出现两次，等等。

【问题讨论】：

this question 根本没有帮助吗？当您搜索解决方案时，您肯定看到了，不是吗？
谢谢，我确实在最初的搜索中看到了这个问题，但我的范围不一定是唯一的，这带来了一个额外的问题。
好的。然后，您必须提供一个可重现的示例，说明您的数据的该功能，以便很明显 findInterval 不起作用。（当然还要解释您希望合并如何工作，因为如果范围不是唯一的，则合并的定义不明确。

标签： r loops merge genetics

【解决方案1】：

Bioconductor 中的GenomicRanges 包专为此类操作而设计。使用例如 read.delim 读取您的数据，以便

con <- textConnection("SNP     BP
rs064   12292
rs319   345367
rs285   700042")
snps <- read.delim(con, head=TRUE, sep="")

con <- textConnection("Gene    BP_start    BP_end
E613    345344      363401
E92     694501      705408
E49     362370      368340") ## missing trailing digit on BP_end??
genes <- read.delim(con, head=TRUE, sep="")

然后从每个中创建“IRanges”

library(IRanges)
isnps <- with(snps, IRanges(BP, width=1, names=SNP))
igenes <- with(genes, IRanges(BP_start, BP_end, names=Gene))

（注意坐标系，IRanges 期望 start 和 end 包含在范围内；另外，当 end = start - 1 时，end >= start 期望 0 宽度范围）。然后找到与基因（'subject'）重叠的 SNP（IRanges 术语中的'query'）

olaps <- findOverlaps(isnps, igenes)

两个 SNP 重叠

> queryHits(olaps)
[1] 2 3

它们与第一个和第二个基因重叠

> subjectHits(olaps)
[1] 1 2

如果一个查询与多个基因重叠，它会在 queryHits 中重复（反之亦然）。然后，您可以将数据框合并为

> cbind(snps[queryHits(olaps),], genes[subjectHits(olaps),])
    SNP     BP Gene BP_start BP_end
2 rs319 345367 E613   345344 363401
3 rs285 700042  E92   694501 705408

通常基因和 SNP 具有染色体和链（“+”、“-”或“*”表示该链不重要）信息，并且您希望在这些上下文中进行重叠；与其创建“IRanges”实例，不如创建“GRanges”（基因组范围），随后的簿记将为您处理

library(GenomicRanges)
isnps <-
    with(snps, GRanges("chrA", IRanges(BP, width=1, names=SNP), "*")
igenes <-
    with(genes, GRanges("chrA", IRanges(BP_start, BP_end, names=Gene), "+"))

【讨论】：

博士。 Morgan - 感谢您向我介绍 GenomicRanges 和 Bioconductor。很酷。随着我的工作向前发展，这些工具将非常方便。我很感激！

【解决方案2】：

您可以使用apply 和which 在base 中按范围合并。只返回那些匹配的（内连接）：

do.call(rbind, apply(file1test, 1, function(x) {
  if(length(idx <- which(x["BP"] >= file2$BP_start & x["BP"] <= file2$BP_end)) > 0) {
    cbind(rbind(x), file2[idx,])
  }
}))
#      SNP     BP  Gene BP_start BP_end
#x  rs2343 860269 E3543   860260 879955
#1   rs754 861822 E3543   860260 879955
#2   rs754 861822   E11   861322 879533
#x1  rs854 367934  E613   367640 368634

或全部来自 file1test（左连接）：

do.call(rbind, apply(file1test, 1, function(x) {
  if(length(idx <- which(x["BP"] >= file2$BP_start & x["BP"] <= file2$BP_end)) > 0) {
    cbind(rbind(x), file2[idx,])
  } else {
    cbind(rbind(x), file2[1,][NA,])
  }
}))
#      SNP     BP  Gene BP_start BP_end
#x  rs2343 860269 E3543   860260 879955
#x1  rs211 369640  <NA>       NA     NA
#1   rs754 861822 E3543   860260 879955
#2   rs754 861822   E11   861322 879533
#x2  rs854 367934  E613   367640 368634
#x3  rs343 706940  <NA>       NA     NA
#x4  rs626 717244  <NA>       NA     NA

【讨论】：

【解决方案3】：

我相信您要的是conditional join。它们在 SQL 中很容易，sqldf 包使使用 SQL 在 R 中查询数据帧变得容易。

只需根据您希望如何处理不匹配的 SNP 来选择一个版本。

内联版：

> sqldf("select * from file1test f1 inner join file2 f2 
+       on (f1.BP > f2.BP_start and f1.BP<= f2.BP_end) ")

输出：

     SNP     BP  Gene BP_start BP_end
1 rs2343 860269 E3543   860260 879955
2  rs754 861822 E3543   860260 879955
3  rs754 861822   E11   861322 879533
4  rs854 367934  E613   367640 368634
>

左加入版本：

> sqldf("select * from file1test f1 left join file2 f2 
+       on (f1.BP > f2.BP_start and f1.BP<= f2.BP_end) ")

输出：

     SNP     BP  Gene BP_start BP_end
1 rs2343 860269 E3543   860260 879955
2  rs211 369640  <NA>       NA     NA
3  rs754 861822 E3543   860260 879955
4  rs754 861822   E11   861322 879533
5  rs854 367934  E613   367640 368634
6  rs343 706940  <NA>       NA     NA
7  rs626 717244  <NA>       NA     NA
>

请注意，如果在 BP 与 BP_start 或 BP_end 完全匹配的情况下，BP 属于哪个组很重要，那么您可能需要小心放置 = 的位置。

【讨论】：