how to assign values to intervals with overlapping regions?

Question

I have two big files, the first file contain some intervals with 85K rows:

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

the second file contains some positions with over 2 million rows:

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

What I want to do is this: For each value in the position file "number" column, search if it is equal to or between ANY of the "numberA" and "numberB" pair values of the data.intervals file.

Additionally, for this "numberA" and "numberB" pair values, its respective "id" must match the "id" in data.position. If this is all true, then I want to insert the respective "id.uniq" in the data.intervals to column of the respective row in data.posiitons file.

There is another problem here as well: some of these intervals overlap with each other and a position may fall within the range of 2 or more than 2 intervals. I want to assign them to each interval separately.

here is the final output that I desire to get (NA means that position does not fall within the range of any intervals):

   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

is there any solution to do this task with a bash or perl script?

Answer 1

You could do this with Perl using the following method:

$ 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

Output:

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

Works:

• Read the intervals file first and build the data structure of a hash keyed on the ID, unique ID, and containing the range endpoints.
• The %order hash stores the order in which the unique IDs were encountered for playback purposes in the same order. OTW, the hash ordering is random.
• Next read the positions file and first unpack each record(or, line) and put them up in $id and $number scalars.
• grep shall select the unique IDs that satisfy the constraint for the number to be in range. Else a "NA" is passed.

Answer 2

This is the kind of thing for which you might consider using a small database - using for example csvsql from csvkit (which also provides a handy csvformat utility as well).

For example, supposing your data is in tab-separated files named intervals and positions, and using the default sqlite dialect:

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

It's somewhat more complicated to get the N/A entries as well: to do so, you can left join the original positions table against the results and look for NULL values of the assigned1 field:

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

Answer 3

Assuming that both files have been sorted using sort -b, you may compose all possible combinations of every line in the two files that have the same ID with

join ranges.txt positions.txt

For simplicity, I'm also ignoring the fact that the files have headers (do consider removing them).

This would produce, for the given data,

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)

Here you have the ID, the "unique ID", the two values of the range, and the position that you'd like to test.

This could be parsed by an awk program:

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"
             }
     }'

This would keep track of the seen IDs and position as keys in a count array. The value would hold the number of times a position has been placed in a range. We need this to be able to say "this position has not been found in any range".

If the current line in the output of join contains a position in the 5th field that is within the range of field 3 and 4, this count is incremented (and the line is outputted).

This outputs all lines in the output of join that corresponds to positions within ranges. The END block handles unmatched positions by looping over the count array and printing out the information that you requested in the question, in the format you requested there whenever the count value is zero.

For the given data, this produces

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

To collapse this data based on the "unique ID", we could pass it through another awk program. (I'm avoiding saving it all in an array in the same awk program as it could be quite a lot of data).

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 }'

This passes through any line whose last column is NA since these are already on the correct format. It then constructs a "key" of the ID and the position. If this key is the same as for the previous line, the "unique ID" is added onto a string called group with a comma as delimiter.

If the key is not the same as the previous, then we've found a new set of ID-position matches and must output the data for the group we've just parsed. This is done again in the END block to output the data for the last group.

Taking all of this together and remembering that both input files need to be sorted, we end up with

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 }'

The output of this is

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

which, apart form the ordering, is identical to what you expected. To fix the ordering, pass this through sort -k1,1n -k2,2n.