我想根据第二个数据框中存在的值过滤数据框。
例如,匹配第一个数据框中"BP"列中高于"start_pos"列的第一个值、小于"end_pos"列或仅小于第二个数据框中"end_pos"的行。
我需要对第二个数据框中的所有值重复此过程。目前,我正在使用 for 循环执行这些操作。但是,我想在单个命令中执行此操作。
数据框 1
CHR BP
29 836019
29 4417047
29 7589996
29 11052921
29 14009294
29 33174196
数据框 2
start_pos end_pos gene_id
19774 19899 ENSBTAG00000046619
34627 35558 ENSBTAG00000006858
69695 71121 ENSBTAG00000039257
83323 84281 ENSBTAG00000035349
124849 179713 ENSBTAG00000001753
264298 264843 ENSBTAG00000005540
for(j in 1:nrow(tmp_markers)){
temp_out_markers<- tmp_markers[j,]
tmp_search<-tmp_gene[which((tmp_markers[j,"BP"]>=tmp_gene[,"start_pos"] & tmp_markers[j,"BP"]<= tmp_gene[,"end_pos"]) | (tmp_markers[j,"BP"]+interval>=tmp_gene[,"start_pos"] & tmp_markers[j,"BP"]+interval <=tmp_gene[,"end_pos"]) | (tmp_markers[j,"BP"]+interval>=tmp_gene[,"start_pos"] & tmp_markers[j,"BP"]+interval <=tmp_gene[,"end_pos"]) | (tmp_markers[j,"BP"]+interval>=tmp_gene[,"start_pos"] & tmp_markers[j,"BP"]+interval >=tmp_gene[,"end_pos"]& tmp_markers[j,"BP"]<=tmp_gene[,"start_pos"])| (tmp_markers[j,"BP"]-interval<=tmp_gene[,"end_pos"] & tmp_markers[j,"BP"]-interval >=tmp_gene[,"start_pos"])|(tmp_markers[j,"BP"]-interval<=tmp_gene[,"end_pos"] & tmp_markers[j,"BP"]-interval<=tmp_gene[,"start_pos"] & tmp_markers[j,"BP"]>=tmp_gene[,"end_pos"])),]
if(nrow(tmp_search)>0){
temp_out<-cbind(temp_out_markers[rep(seq_len(nrow(tmp_search))),],tmp_search)
temp_out[,"Distance_from_gene_start"]<-temp_out[,"BP"]-temp_out[,"start_pos"]
temp_out[,"Distance_from_gene_end"]<-temp_out[,"BP"]-temp_out[,"end_pos"]
output_genes<-rbind(temp_out,output_genes)
}
}
最后,我想要一个数据框,其中包含测试间隔内的所有行。
正如我在评论中所说,您的模拟数据不会导致匹配,因为最小的BP
值 (836019( 大于最大的end_pos
(264843(。
也可能是我完全误解了你的问题!
我知道您想将df1
中的行与BP >= start_pos
和BP <= end_pos
df2
中的行相匹配.如果是这样,我们可以使用包data.table
提供的非 equi 连接来实现这一点。
library(data.table)
result <- dt1[dt2,
.(BP, CHR, gene_id),
on = .(BP >= start_pos, BP <= end_pos),
nomatch = NULL,
by = .EACHI]
setnames(result, 1:2, names(dt2)[1:2])
result
start_pos end_pos BP CHR gene_id
1: 0.000000 2.000000 0 29 ABCD01
2: 4.571429 6.571429 6 30 ABCD03
3: 11.428571 13.428571 12 31 ABCD06
4: 16.000000 18.000000 18 32 ABCD08
5: 22.857143 24.857143 24 33 ABCD011
6: 29.714286 31.714286 30 34 ABCD014
如果您需要完整的 15 行dt2
,只需省略nomatch = NULL
部分即可。
使用的数据:
dt1 <- data.table(CHR = 29:34,
BP = seq(0, 30, length.out = 6),
key = "BP")
dt2 <- data.table(start_pos = seq(0, 32, length.out = 15),
gene_id = paste0("ABCD", rep(0, 3), 1:15))
dt2[, end_pos := start_pos + 2]
setcolorder(dt2, c(1, 3, 2))
foverlaps
的替代方案
正如@r2evans评论中提到的,data.table
还有另一个功能,foverlaps
在这里可能很有用。它检查一个范围是否与另一个表中的范围重叠,因此我们需要做一个小技巧来创建一个 0 宽度的范围dt1
:
dt1[, BP2 := BP]
我们还需要有键控数据表:
setkey(dt1, "BP", "BP2")
setkey(dt2, "start_pos", "end_pos")
然后将其传递给工作马:
foverlaps(dt1, dt2)
start_pos end_pos gene_id CHR BP BP2
1: 0.000000 2.000000 ABCD01 29 0 0
2: 4.571429 6.571429 ABCD03 30 6 6
3: 11.428571 13.428571 ABCD06 31 12 12
4: 16.000000 18.000000 ABCD08 32 18 18
5: 22.857143 24.857143 ABCD011 33 24 24
6: 29.714286 31.714286 ABCD014 34 30 30
当然,我们以后可以通过BP2 := NULL
摆脱BP2
。
如果我们想要完整的 15 行dt2
,它只是颠倒了调用中对象的顺序:
foverlaps(dt2, dt1)
非常感谢!
我以这个解决方案结束,它运行良好。
foverlaps(tmp_gene, tmp_markers, by.x = c("start_pos","end_pos"), by.y =
key(tmp_markers),nomatch = 0)
干杯。