☰

如何为具有重叠区域的间隔赋值？

tag-icon • tag-icon text-processing command-line perl
如何为具有重叠区域的间隔赋值？

我有两个大文件，第一个文件包含一些 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。

相关内容