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

Annotation of /Sprout/NeighboringFeatures.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

1 : parrello 1.1 #!/usr/bin/perl -w
2 :     use strict;
3 :    
4 :     use ERDB;
5 :     use BasicLocation;
6 :    
7 :     if (! @ARGV) {
8 :     print "usage: NeighboringFeatures <alias or id> <size>\n";
9 :     } else {
10 :     my $erdb = ERDB::GetDatabase('Sprout');
11 :     my $chunkSize = $erdb->MaxSegment();
12 :     my ($peg_alias_or_id, $size) = @ARGV;
13 :     my $feature;
14 :     if ($peg_alias_or_id =~ /^fig\|/) {
15 :     # Look for a feature with the given ID.
16 :     $feature = $erdb->GetEntity(Feature => $peg_alias_or_id);
17 :     } else {
18 :     my ($peg_id) = $erdb->GetFlat('FeatureAlias IsAliasOf Feature', "FeatureAlias(id) = ?", [$peg_alias_or_id],
19 :     'Feature(id)');
20 :     if (defined $peg_id) {
21 :     # Look for a feature with the computed ID.
22 :     $feature = $erdb->GetEntity(Feature => $peg_id);
23 :     }
24 :     }
25 :     if (! defined $feature) {
26 :     print "No feature found named \"$peg_alias_or_id\".\n";
27 :     } else {
28 :     # Keep the features we've found in here.
29 :     my %found;
30 :     # Get the location.
31 :     my @locations = split /,/, $feature->PrimaryValue('location-string');
32 :     for my $location (@locations) {
33 :     my $loc = BasicLocation->new($location);
34 :     # Now we need to come up with a way to ask if a location overlaps the
35 :     # BasicLocation we just created. The obvious filter is
36 :     #
37 :     # IsLocatedIn(beg) + IsLocatedIn(len) >= $left AND
38 :     # IsLocatedIn(beg) <= $right
39 :     #
40 :     # This, however, will be slow. We know the maximum size of an
41 :     # IsLocatedIn region is $chunksize, so we can add
42 :     #
43 :     # IsLocatedIn(beg) + $chunkSize >= $left
44 :     #
45 :     # to the query. This enables the SQL processor to grab a section
46 :     # of the index and look inside that.
47 :     # So, first compute the region of interest.
48 :     my $left = $loc->Left - $size;
49 :     my $right = $loc->Right + $size;
50 :     # Ask for the features.
51 :     my $qh = $erdb->Get('Feature IsLocatedIn',
52 :     'IsLocatedIn(beg) + ? >= ? AND ' .
53 :     'IsLocatedIn(beg) <= ? AND ' .
54 :     'IsLocatedIn(beg) + IsLocatedIn(len) >= ?',
55 :     [$chunkSize, $left, $right, $left]);
56 :     while (my $feature = $qh->Fetch()) {
57 :     # Get this feature's data.
58 :     my ($id, $location) = $feature->Values(['id', 'location-string']);
59 :     print "$id\t$location\n";
60 :     }
61 :     }
62 :     }
63 :     }
64 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3