[Bio] / Sprout / NeighboringFeatures.pl Repository:
ViewVC logotype

View of /Sprout/NeighboringFeatures.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (download) (as text) (annotate)
Mon Jan 19 21:49:02 2009 UTC (10 years, 1 month ago) by parrello
Branch: MAIN
CVS Tags: mgrast_dev_08112011, mgrast_dev_08022011, rast_rel_2014_0912, mgrast_dev_04082011, rast_rel_2010_0928, mgrast_version_3_2, mgrast_dev_12152011, mgrast_dev_06072011, rast_rel_2009_0925, rast_rel_2010_0526, rast_rel_2014_0729, rast_rel_2009_05_18, rast_rel_2010_1206, mgrast_release_3_0, mgrast_dev_03252011, rast_rel_2010_0118, rast_rel_2009_02_05, rast_rel_2011_0119, mgrast_release_3_0_4, mgrast_release_3_0_2, mgrast_release_3_0_3, mgrast_release_3_0_1, mgrast_dev_03312011, mgrast_release_3_1_2, mgrast_release_3_1_1, mgrast_release_3_1_0, mgrast_dev_04132011, mgrast_dev_04012011, rast_rel_2009_07_09, rast_rel_2010_0827, myrast_33, rast_rel_2011_0928, mgrast_dev_04052011, rast_rel_2009_03_26, mgrast_dev_10262011, HEAD
Tutorial examples.

#!/usr/bin/perl -w
use strict;

use ERDB;
use BasicLocation;

if (! @ARGV) {
    print "usage: NeighboringFeatures <alias or id> <size>\n";
} else {
    my $erdb = ERDB::GetDatabase('Sprout');
    my $chunkSize = $erdb->MaxSegment(); 
    my ($peg_alias_or_id, $size) = @ARGV;
    my $feature;
    if ($peg_alias_or_id =~ /^fig\|/) {
        # Look for a feature with the given ID.
        $feature = $erdb->GetEntity(Feature => $peg_alias_or_id);
    } else {
        my ($peg_id) = $erdb->GetFlat('FeatureAlias IsAliasOf Feature', "FeatureAlias(id) = ?", [$peg_alias_or_id],
                                      'Feature(id)');
        if (defined $peg_id) {
            # Look for a feature with the computed ID.
            $feature = $erdb->GetEntity(Feature => $peg_id);
        }
    }
    if (! defined $feature) {
        print "No feature found named \"$peg_alias_or_id\".\n";
    } else {
        # Keep the features we've found in here.
        my %found;
        # Get the location.
        my @locations = split /,/, $feature->PrimaryValue('location-string');
        for my $location (@locations) {
            my $loc = BasicLocation->new($location);
            # Now we need to come up with a way to ask if a location overlaps the
            # BasicLocation we just created. The obvious filter is
            #
            # IsLocatedIn(beg) + IsLocatedIn(len) >= $left AND
            # IsLocatedIn(beg) <= $right
            #
            # This, however, will be slow. We know the maximum size of an
            # IsLocatedIn region is $chunksize, so we can add
            #
            # IsLocatedIn(beg) + $chunkSize >= $left
            #
            # to the query. This enables the SQL processor to grab a section
            # of the index and look inside that.
            # So, first compute the region of interest.
            my $left = $loc->Left - $size;
            my $right = $loc->Right + $size;
            # Ask for the features.
            my $qh = $erdb->Get('Feature IsLocatedIn',
                                'IsLocatedIn(beg) + ? >= ? AND ' .
                                  'IsLocatedIn(beg) <= ? AND ' .
                                  'IsLocatedIn(beg) + IsLocatedIn(len) >= ?',
                                [$chunkSize, $left, $right, $left]);
            while (my $feature = $qh->Fetch()) {
                # Get this feature's data.
                my ($id, $location) = $feature->Values(['id', 'location-string']);
                print "$id\t$location\n";
            }
        }
    }
}


MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3