Find nearest range higher and lower inside another range

I have rewritten this question many times, as while writing it out, I have come up with some ideas and have managed to get a very slow, but working solution. It seems writing out what one needs to do tends to help with the problem. However I am still having some issues with speed and would very much like some insight/input/comments/help. Anyway, here is what I am trying to do:

I have two files of genome coordinate data. File1 contains non-coding coordinates and file2 contains coding coordinates. Both contain five fields. The first field contains a uniqID. The second and third fields contain the start and stop coordinates respectively. The fourth field contains the chromosome ID. There can sometimes be many lines with the same chromosome ID (as there are multiple non-coding and coding ranges per chromosome). The fifth field contains the length/span of the coordinates.

For each set of coordinates in file2 I wish to find the closest (below and above, but currently just working on the "below" first and will reuse code to deal with the "above") set of coordinates in file1 matching the same ID that also has a length/span large enough to hold the coordinates from file2.

I then want to offset the coordinates by a random number that can have a min value of 1 through to a max value of the difference in size of the two coordinate span lengths.

To explain some more, here are some lines for each file.

File1 - non-coding

chr_17_0_3614        0      3614   chr_17 3614
chr_17_16177_17009   16177  17009  chr_17 832
chr_20_0_103982      0      103982 chr_20 103982
chr_20_105978_107383 105978 107383 chr_20 1405
chr_20_107502_108744 107502 108744 chr_20 1242
chr_20_108897_110439 108897 110439 chr_20 1542

File2 - coding

chr_17_3614_4570     3614   4570   chr_17 956
chr_17_7901_8585     7901   8585   chr_17 684
chr_17_10351_10698   10351  10698  chr_17 347
chr_20_104882_104971 104882 104971 chr_20 89
chr_20_105763_105978 105763 105978 chr_20 215
chr_20_107383_107502 107383 107502 chr_20 119

What I think I need to do is:

1. Split non-coding file into multiple files based on chromosome ID.
2. Run a loop that takes in each line of coding.txt and assigns variables to each field.
3. Check if coding 'start' value is higher or equal to non-coding 'end' and the 'span' of coding coordinates is less than the 'span' of the non-coding.
4. If true, then set variables for holding needed values for printing.
5. Once the coding 'start' value is less than the non-coding 'end' value we can stop looking and print out the necessary values (after some modification for randomness) and then move onto the next line.

I have written some awk that looks like it works, but it is in need of some speed increases. Here is some output from the above I get with this script:

OffsetBelow_chr_17_3614_4570     1621   2577   chr_17
OffsetBelow_chr_17_7901_8585     293    977    chr_17
OffsetBelow_chr_17_10351_10698   3010   3357   chr_17
OffsetBelow_chr_20_104882_104971 77516  77605  chr_20
OffsetBelow_chr_20_105763_105978 11841  12056  chr_20
OffsetBelow_chr_20_107383_107502 106549 106668 chr_20

And here is the script I wrote:

#!/usr/bin/bash

# split non-coding into individual chromosome files
awk '{print > "split/"$4".txt"}' non_coding.txt

# find nearest coords
while read -r uniqId start stop chrID span; do
 awk -v rseed="$RANDOM" -v qStart="$start" -v qStop="$stop" -v uId="$uniqId" -v span="$span" 'BEGIN {
  belowEnd = 0; belowStart = 0; availSpan = 0; belowUniqId = "none"; srand(rseed)
  }
  {if (qStart < $3) {next} else if (qStart >= $3 && span <= $5)
  {belowEnd = $3; belowStart = $3-span; availSpan = $5; belowUniqId = $1;
  min=1; max=$5-span; offset=int(min+rand()*(max-min+1))}} END {
  print "OffsetBelow_"uId, belowStart-offset, belowEnd-offset, $4}' split/"$chrID".txt
done < coding.txt

So, it works, but is super slow as I'm running over the same file many, many times (coding.txt is about 100,000 lines long), and I'm just not sure where to go from here to make it work better. Also, I need to now make this work for the nearest span 'above' (which I'm pretty sure should be mostly a jostling of variables and operators), and maybe somehow that can work in with making it quicker... Anyway, some comments/thoughts/ideas to make this process run quicker would be much appreciated. Thanks for your time and brain power!

> So, it works, but is super slow as I'm running over the same file many, many times (coding.txt is about 100,000 lines long), and I'm just not sure where to go from here to make it work better.

You do save yourself some time by splitting the non-coding regions into per-chromosome files, but the fundamental problem here is that you are performing a separate pass for each coding region. It would be cleaner and much faster to make a single, joint pass through both of the original files to extract the wanted information for all the coding regions in one go.

There are several ways to arrange the details, but I'd probably go with something along these lines:

• Add a field to each line of each of the files to indicate whether it is a coding or non-coding region, then
• sort the two files together into a single one, ordered by (chromosome id, start coordinate)
• As you read through that file line by line, track the last-seen non-coding region
• When you see a coding region, you know that either
  • there is no preceding non-coding region (if there is not yet a non-coding region tracked, or if the tracked one belongs to a different chromosome), OR
  • the non-coding region immediately preceding this coding region is the one you currently have tracked

With a bit more attention you could also extract each coding region's nearest subsequent non-coding region from the same pass, but you could also just get the subsequent regions by reversing the sort order and doing exactly the same thing as before.

Partial example:

#!/bin/bash

sort -k4,4 -k2,2n \
  <(sed '/^\s*$/d; s/$/ N/' non_coding.txt) \
  <(sed '/^\s*$/d; s/$/ C/' coding.txt) |
  while read -r uid start stop chr length kind; do
    case $kind in
      C) # handle a coding region
        ;;
      N) # handle a non-coding region
        ;;
      *) # report an error
        ;;
    esac
  done

The sort -k4,4 -k2,2n expresses a line-by-line sort (of the command's standard input), using a primary key spanning fields 4 through 4, and a secondary, numeric key spanning fields 2 through 2. Results are written to the command's standard output.

The <(sed '/^\s*$/d; s/$/ N/' non_coding.txt) is a bash-ism that expands to the name of a file from which the output of sed '/^\s*$/d; s/$/ N/' non_coding.txt can be read. The sed command reads input from non_coding.txt, drops empty / blank lines (/^\s*$/d), and appends a space and an 'N' to all other lines (s/$/ N/), which are written to stdout.

The <(sed '/^\s*$/d; s/$/ C/' coding.txt) is analogous to the preceding.

The read -r uid start stop chr length kind reads one line, splits it into shell words, and assigns those, in order, to the specified shell variables. In particular, this puts the region-type indicator into $kind, which is used inside the loop to choose the appropriate behavior.

... and that's where I leave the rest of the details to you.

Taking into account John and Ed's help, I think I've managed to come up with a solution.

First, add the two files together into one, add an identifier of coding/non-coding as a sixth field and sort on field 4 then 2. Then run an awk script. This is super fast. So thanks again for your help in direction John and Ed.

It seems to work as I want, but I'll need to some more testing to be fully sure. Anyway, for those who are interested and anyone else who stumbles upon this, this is what I came up with. If anyone sees any issues with this, let me know.

#!/bin/bash
awk '
BEGIN {OFS="\t"}
{
    # Store each line in an array indexed by record id
    lines[$1] = $0;
    
    # Store record ids, start, ref, length, and type in separate arrays
    ids[NR] = $1;
    starts[$1] = $2;
    refs[$1] = $4;
    lengths[$1] = $5;
    types[$1] = $6;
}

END {
    srand();
    for(i = 1; i <= NR; i++) { 
        # Only consider coding regions
        if(types[ids[i]] == "C") {
            
            # Find nearest below non-coding region
            for(j = i - 1; j > 0; j--) {
                if(types[ids[j]] == "N" && refs[ids[i]] == refs[ids[j]]) {
                    OffNeg = ids[j];
                    break;
                }
            }

            # Find nearest above non-coding region
            for(j = i + 1; j <= NR; j++) {
                if(types[ids[j]] == "N" && refs[ids[i]] == refs[ids[j]]) {
                    OffPos = ids[j];
                    break;
                }
            }

            # Check that we found both above and below regions
            if(OffPos && OffNeg) {
                # Determine length of coding region
                len = lengths[ids[i]];

                # Make sure both non-coding regions are long enough
                if(lengths[OffPos] >= len && lengths[OffNeg] >= len) {

                    # Calculate random start positions within each non-coding region
                    # that allow for a segment of equal length to the coding region
                    OffPosStart = starts[OffPos] + int(rand() * (lengths[OffPos] - len + 1));
                    OffNegStart = starts[OffNeg] + int(rand() * (lengths[OffNeg] - len + 1));

                    # Output the new record ids with the calculated start and end positions
                    print "OffNeg_" ids[i], OffNegStart, OffNegStart+len, refs[ids[i]]
                    print "OffPos_" ids[i], OffPosStart, OffPosStart+len, refs[ids[i]]
                }
            }

            # Reset OffPos and OffNeg for the next coding region
            OffPos = "";
            OffNeg = "";
        }
    }
}' file.txt