主页 / unix / 问题 / 510325
Accepted
Anna1364
Asked: 2019-04-04 08:58:23 +0800 CST

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

  • 772

我有两个大文件，第一个文件包含一些 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 脚本完成此任务？

text-processing command-line

3 个回答

  1. Best Answer

    您可以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。否则"NA"通过。
    • 2

  2. 对于这种情况，您可能会考虑使用小型数据库——例如使用csvsqlcsvkit （它也提供了一个方便的csvformat实用程序）。

    例如，假设您的数据位于名为intervalsand的制表符分隔文件中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
    • 1

  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.

    • -1

相关问题