R中的范围合并 - 应用循环

Question

我在这里发布了一个问题: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行的范围内,但只有一个属于第二行的范围.

第一个包含SNP的文件只有大约400行.但是,包含范围的第二个文件大约有20K.我想要作为输出产生的是一个数据框,其中包含来自第一个文件(SNP)的行,其中BP属于第二个文件中的BP范围.如果SNP落入两个范围,那么它将出现两次,等等.

Answer 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期望开始和结束都包含在范围内;同时,end> =当end = start时开始期望0宽度范围 - 1).然后找到与基因('主题')重叠的SNP(IRanges术语中的'查询')

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), "+"))

Answer 2

我相信你所要求的是一个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所属的组.