快速版本是我想在特定值的特定距离内找到所有匹配项。示例:
library(data.table)
reads <- structure(list(seqnames = c(1, 1, 1, 1, 1, 1, 1), start = c(100,
100, 100, 130, 130, 132, 133), end = c(130, 130, 131, 160, 160,
155, 160), strand = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "+", class = "factor")), row.names = c(NA,
-7L), class = c("data.table", "data.frame"), .internal.selfref = <pointer: 0x15af570>, sorted = c("start",
"end", "strand"))
# take only the end position
p3 <- reads[, list(N3 = .N), by = c("seqnames", "end", "strand")]
setnames(p3, "end", "loc.3p")
# take only the start position
p5 <- reads[, list(N5 = .N), by = c("seqnames", "start", "strand")]
setnames(p5, "start", "loc.5p")
# find the start positions close to the end positions
p3[, loc := loc.3p]
p5[, loc := loc.5p]
p3[p5, roll=50, rollends=c(FALSE, TRUE), nomatch=0, on = c("seqnames", "strand", "loc")]
seqnames loc.3p strand N3 loc loc.5p N5
1: 1 130 + 2 130 130 2
2: 1 131 + 1 132 132 1
3: 1 131 + 1 133 133 1
这是有效的,但只返回第一个匹配项。我现在的目标是找到所有匹配项(在指定的50范围内(。预期输出:
seqnames loc.3p strand N3 loc loc.5p N5
1: 1 130 + 2 130 130 2
2: 1 130 + 2 130 132 1
3: 1 130 + 2 133 133 1
4: 1 131 + 1 132 132 1
4: 1 131 + 1 133 133 1
注意CCD_ 1 130应该如何也与CCD_ 2 132和132匹配。
如何做到这一点?我需要为下一组操作保留分数。
生物信息学版本,我正试图找到所有下游5'读数到一定距离(在同一条链中(。对于这个例子,我只使用"+"链,但我也必须使用"-"链。由于这将用于数百万次读取,因此data.table
似乎是合适的。我研究了GenomicRanges
,但保留两组数据(5'和3'位置(的元数据有点复杂。
当你说时
在50 的指定范围内
这让我觉得不平等。在未来,如果你展示你想要的结果,这可能会有所帮助:(
对p3
稍作修改,注意我跳过了创建您所拥有的loc
列。
p3[, "upper" := 50 + loc.3p]
非平等加入:
p3[p5, .(seqnames, x.loc.3p, strand, N3, loc.5p, N5),
on = .(seqnames, strand, loc.3p <= loc.5p, upper >= loc.5p), nomatch = 0
][order(x.loc.3p, loc.5p),
.(seqnames, loc.3p = x.loc.3p, strand, N3, loc.5p, N5)
]
seqnames loc.3p strand N3 loc.5p N5
1: 1 130 + 2 130 2
2: 1 130 + 2 132 1
3: 1 130 + 2 133 1
4: 1 131 + 1 132 1
5: 1 131 + 1 133 1
最后一个括号用于排序和重命名,我不完全确定为什么我必须在第一个括号中指定x.loc.3p
,但我在out it…中出错