![如何为具有重叠区域的间隔赋值?](https://linux22.com/image/151155/%E5%A6%82%E4%BD%95%E4%B8%BA%E5%85%B7%E6%9C%89%E9%87%8D%E5%8F%A0%E5%8C%BA%E5%9F%9F%E7%9A%84%E9%97%B4%E9%9A%94%E8%B5%8B%E5%80%BC%EF%BC%9F.png)
我有两个大文件,第一个文件包含一些 85K 行的间隔:
head data.intervals
id id_uniq numberA numberB
1 g1 5 20
1 g2 6 29
1 g3 17 35
1 g4 37 46
1 g5 50 63
1 g6 70 95
1 g7 87 93
2 g8 3 15
2 g9 10 33
2 g10 60 77
2 g11 90 132
第二个文件包含一些超过 200 万行的位置:
head data.posiitons
id number
1 4
1 19
1 36
1 49
1 90
2 1
2 20
2 89
2 93
2 120
我想要做的是:对于位置文件“number”列中的每个值,搜索它是否等于或介于 data.intervals 文件的“numberA”和“numberB”对值中的任何一个之间。
此外,对于“numberA”和“numberB”对值,其各自的“id”必须与 data.position 中的“id”匹配。如果这一切都是真的,那么我想将 data.intervals 中的相应“id.uniq”插入到 data.posiitons 文件中相应行的列中。
这里还存在另一个问题:这些间隔中的一些间隔彼此重叠,并且一个位置可能落入2个或2个以上间隔的范围内。我想将它们分别分配给每个间隔。
这是我希望获得的最终输出(NA 意味着该位置不属于任何间隔的范围内):
id number assigned1
1 4 NA
1 19 g1,g2,g3
1 36 NA
1 49 NA
1 90 g6,g7
2 1 NA
2 20 g9
2 89 NA
2 93 g11
2 120 g11
有没有解决方案可以使用 bash 或 perl 脚本完成此任务?
答案1
您可以Perl
使用以下方法来完成此操作:
$ perl -lane '
my($id, $uniq_id, $lower, $upper) = @F;
$h{$id}{$uniq_id}{MIN} = $lower;
$h{$id}{$uniq_id}{MAX} = $upper;
push @{$order{$id}}, $uniq_id;
}{
while(<STDIN>) {
chomp;
my($id, $number) = split;
print join "\t", $id, $number,
join(",", grep { $h{$id}{$_}{MIN} < $number and $h{$id}{$_}{MAX} > $number } @{$order{$id}})
|| qw/NA/;;
}
' data.intervals < data.posiitons
输出:
1 4 NA
1 19 g1,g2,g3
1 36 NA
1 49 NA
1 90 g6,g7
2 1 NA
2 20 g9
2 89 NA
2 93 g11
2 120 g11
作品:
- 首先读取间隔文件,并构建以 ID、唯一 ID 为键并包含范围端点的哈希数据结构。
- 散列
%order
存储遇到唯一 ID 的顺序,以便以相同的顺序进行播放。 OTW,哈希顺序是随机的。 - 接下来读取位置文件并首先解压每个记录(或行)并将它们放入 $id 和 $number 标量中。
grep
应选择满足范围内的数字约束的唯一 ID。否则 a"NA"
被传递。
答案2
对于这种情况,您可能会考虑使用小型数据库 - 例如使用csvsql
来自csvkit(这也提供了一个方便的csvformat
实用程序)。
例如,假设您的数据位于名为intervals
和的制表符分隔文件中positions
,并使用默认sqlite
方言:
csvsql --tabs --query '
SELECT id,number,group_concat(id_uniq) AS "assigned1"
FROM positions JOIN intervals USING(id)
WHERE number >= numberA AND number <= numberB
GROUP BY id,number ORDER BY id,number
' positions intervals | csvformat --out-tabs
id number assigned1
1 19 g1,g2,g3
1 90 g6,g7
2 20 g9
2 93 g11
2 120 g11
获取条目也有些复杂N/A
:为此,您可以将原始positions
表与结果左连接并查找字段NULL
的值assigned1
:
csvsql --tabs --query '
SELECT id,number,IFNULL(assigned1,"NA") assigned1 FROM positions
LEFT JOIN (
SELECT id,number,group_concat(id_uniq) AS "assigned1"
FROM positions JOIN intervals USING(id)
WHERE number >= numberA AND number <= numberB
GROUP BY id,number
) USING(id,number) ORDER BY id,number
' positions intervals | csvformat --out-tabs
id number assigned1
1 4 NA
1 19 g1,g2,g3
1 36 NA
1 49 NA
1 90 g6,g7
2 1 NA
2 20 g9
2 89 NA
2 93 g11
2 120 g11
答案3
假设两个文件都已使用 进行排序sort -b
,则可以使用以下命令组合两个文件中具有相同 ID 的每一行的所有可能组合
join ranges.txt positions.txt
为了简单起见,我还忽略了文件具有标题的事实(请考虑删除它们)。
对于给定的数据,这将产生:
1 g1 5 20 4
1 g1 5 20 19
1 g1 5 20 36
1 g1 5 20 49
1 g1 5 20 90
1 g2 6 29 4
1 g2 6 29 19
1 g2 6 29 36
[...] (in total 55 lines)
这里有 ID、“唯一 ID”、范围的两个值以及您要测试的位置。
这可以由程序解析awk
:
join ranges.txt positions.txt |
awk '!($1 SUBSEP $5 in count) { count[$1,$5]=0 }
$5 >= $3 && $5 <= $4 && ++count[$1,$5]
END {
for (i in count)
if (count[i] == 0) {
split(i,s,SUBSEP)
print s[1], s[2], "NA"
}
}'
这将跟踪看到的 ID 和位置作为数组中的键count
。该值将保存一个位置被放置在一个范围内的次数。我们需要这个能够说“这个位置已经不是在任何范围内都被发现”。
如果输出中的当前行join
包含第 5 个字段中位于字段 3 和 4 范围内的位置,则该计数将递增(并且输出该行)。
这会输出输出中join
与范围内的位置相对应的所有行。该END
块通过循环遍历数组并以您在值为零时count
请求的格式打印出您在问题中请求的信息来处理不匹配的位置。count
对于给定的数据,这会产生
1 g1 5 20 19
1 g2 6 29 19
1 g3 17 35 19
1 g6 70 95 90
1 g7 87 93 90
2 g9 10 33 20
2 g11 90 132 93
2 g11 90 132 120
2 89 NA
1 36 NA
1 4 NA
2 1 NA
1 49 NA
到坍塌该数据基于“唯一ID”,我们可以将其传递给另一个awk
程序。 (我避免将其全部保存在数组中相同 awk
程序,因为它可能包含大量数据)。
awk '$NF == "NA" { print; next }
{ key = $1 SUBSEP $NF }
key == prev { group = group "," $2; next }
{ if (group != "") print id, pos, group; id = $1; pos = $NF; group = $2; }
END { if (group != "") print id, pos, group }'
这会穿过最后一列为的任何行,NA
因为它们已经采用正确的格式。然后它构造一个 ID 和位置的“键”。如果此键与上一行相同,则“唯一 ID”将添加到以group
逗号作为分隔符调用的字符串中。
如果关键是不是与之前相同,那么我们找到了一组新的 ID-位置匹配,并且必须输出我们刚刚解析的组的数据。在块中再次执行此操作END
以输出最后一组的数据。
将所有这些放在一起并记住两个输入文件都需要排序,我们最终得到
join ranges.txt positions.txt |
awk '!($1 SUBSEP $5 in count) { count[$1,$5]=0 }
$5 >= $3 && $5 <= $4 && ++count[$1,$5]
END {
for (i in count)
if (count[i] == 0) {
split(i,s,SUBSEP)
print s[1], s[2], "NA"
}
}' |
awk '$NF == "NA" { print; next }
{ key = $1 SUBSEP $NF }
key == prev { group = group "," $2; next }
{ if (group != "") print id, pos, group
prev = key; id = $1; pos = $NF; group = $2; }
END { if (group != "") print id, pos, group }'
这个的输出是
1 19 g1,g2,g3
1 90 g6,g7
2 20 g9
2 93 g11
2 89 NA
1 36 NA
1 4 NA
2 1 NA
1 49 NA
2 120 g11
除了顺序之外,它与您的预期相同。要修复顺序,请将其传递给sort -k1,1n -k2,2n
。