Re: What is the best approach to check if a position falls within a target range?

[umasuresh]

Hi Monks, I am posting an updated version of the code which significantly reduced the time for assigning target status (less than 1/2 the time compared to the original version).

use strict;
use Data::Dumper;
use List::Util qw/max/;

# target region
my @a = ("100_200","210_310","400_450","475_600", "680_900");
my $query_snp = 205;

my $line = quicksort(\@a, \$query_snp);

print "$$line";

# loosely based on the quicksort algorithm (Mastering Alorithms with P
+erl)
sub quicksort
{
    my $array     = shift;
    my $snp       = shift;
    my $start_idx = 0;
    my $end_idx = scalar @$array - 1; 
    my $mid_point    = int( ( $end_idx - $start_idx)/2 );
    my $pivot = $start_idx + $mid_point;
    
    my ($current_start , $current_end) = split /\_/, $array->[$pivot];
    my $current_start_idx = $$snp < $current_start ? $start_idx: $pivo
+t; 
    my $current_end_idx   = $current_start_idx ==$start_idx ? $pivot: 
+$end_idx;
    
    my $new_array = [ @$array[$current_start_idx..$current_end_idx] ];
    my $out_line; 
    
    
    
    # check if the query_snp falls in the pivot region
    if($end_idx > 1 && $$snp >= $current_start && $$snp <= $current_en
+d)
    {
        $out_line = "$$snp\t$current_start\t$current_end";
        return(\$out_line);
    }
    # If there are only two elements in the input array, three cases a
+re possible
    # case1: (100_200, 280_380)
    #    a. query_snp <  100
    #     b. query_snp > 380
    #    c. query_snp > 200 && query_snp < 280
    elsif($end_idx ==1 )
    {
        if( $$snp > ( split /\_/, $array->[0] )[1] && $$snp < ( split 
+/\_/, $array->[1] )[0] ||
            $$snp > ( split /\_/, $array->[1] )[1] ||
            $$snp < ( split /\_/, $array->[0] )[0] )
        {
            $out_line = "$$snp\tNOT_IN_TARGET";
            return(\$out_line);
        }
        
        # case 2: query _snp lies in one of the target regions
        #     a. query_snp >=100 && <=200
        #    b. query_snp >=280 && <=380
        elsif( $$snp >= ( split /\_/, $array->[0] )[0] && $$snp <= ( s
+plit /\_/, $array->[0] )[1]
            || $$snp >= ( split /\_/, $array->[1] )[0] && $$snp <= ( s
+plit /\_/, $array->[1] )[1]
            )
        {
            $out_line = "$$snp\tyes";
            return(\$out_line);
        }
        # case3: Doesn't satisfy any of the above and hence is an erro
+r
        else
        {
            $out_line = "$$snp\terror\n";
            return(\$out_line);
        }


    }
    # if it doesn't lie in the pivot region, do a recursive check of h
+alf the array elements each time it enters this loop
    # this significantly reduces the search time
    else
    {
        quicksort($new_array, $snp);
    }

Any input is greatly appreciated!

[umasuresh]

Hi Monks,
I would very much appreciate any suggestions on optimizing this algorithm. Typically, the query file contains some 2 million records and the target region has about 200,000 records. In these cases, it takes several hours to assign the target status even with the quicksort algorithm. I am thinking of first dividing the target array into pivot regions of 12.5, 25, 50 75, 87.5 and 100, then checking to see if a query snp lies in a bigger chunk before passing the subset array into the quicksort function. It appears that the recursion in the function is a big drag!
Is there a quicker way of achieving this?
Thanks Much,
Uma

[BrowserUk]

the query file contains some 2 million records and the target region has about 200,000 records.

Could you post 3 query file records and 3 target region records?

[umasuresh]

Hi BrowserUk,

query
=====
100
200
300

target_region
=============
10_50
60_150
180_250

expected output
===============
100 60_150
200 180_250

Thanks Much, Uma

[BrowserUk]

FWIW. I've a solution that matches 2e6 random integers (0 .. 1000) against 200,000 randomly generated ranges (0..700, 1..300) in 45 minutes using 3GB of ram.

Of course, of those 2e6 queries only 1000 are unique so it's doing 2000 more work than it needs to.

[umasuresh]

Hi BrowserUk,
Thanks much for all your replies. I did some major refactoring of the code these past few days and was able to achieve significant improvement in speed.
Major changes are:
1. split the target region into 24 chunks for 24 chromosomes (chr) and only loaded the chr of interest in memory.
2. Converted the target hash to a target array. This caused a big gain in speed.
UPDATE
3.Divided target region into 8 chunks 10-12.5-25-50 and so on. Checked the if the query snp is in which chunk before assigning target status. Uma