code_perl/Maasha/NClist.pm

   1 package Maasha::NClist;
   2
   3 # Copyright (C) 2010 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Nested Containment List (NCList): A new algorithm for accelerating
  26 # interval query of genome alignment and interval databases
  27 # http://bioinformatics.oxfordjournals.org/cgi/content/abstract/btl647v1
  28
  29 # Nested lists are composed of intervals defined by begin and end positions,
  30 # and any interval that is contained within another interval is rooted to this.
  31 # Thus, fast interval lookups can be performed using binary search.
  32
  33
  34 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  35
  36
  37 use warnings;
  38 use strict;
  39
  40 use Maasha::Filesys;
  41 use Data::Dumper;
  42 use Time::HiRes;
  43 use JSON::XS;
  44
  45 use vars qw( @ISA @EXPORT );
  46
  47 @ISA = qw( Exporter );
  48
  49
  50 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  51
  52
  53 sub nc_list_create
  54 {
  55     # Martin A. Hansen, February 2010.
  56
  57     # Creates a Nested Containment (NC) list from a stack of features.
  58     # The features consits of an AoA where beg and end specifies the
  59     # elements containing the begin and end position of a feature, and
  60     # index specifies a element used for nesting lists.
  61
  62     my ( $features,   # list of features AoA
  63          $beg,        # feature element with begin position
  64          $end,        # feature element with end position
  65          $index,      # feature element with index position
  66        ) = @_;
  67
  68     # Returns a list.
  69
  70     my ( $nc_list );
  71
  72     @{ $features } = sort { $a->[ $beg ] <=> $b->[ $beg ] or $b->[ $end ] <=> $a->[ $end ] } @{ $features };
  73
  74     $nc_list = [ shift @{ $features } ];
  75
  76     map { nc_list_add( $nc_list, $_, $end, $index ) } @{ $features };
  77
  78     return wantarray ? @{ $nc_list } : $nc_list;
  79 }
  80
  81
  82 sub nc_list_add
  83 {
  84     # Martin A. Hansen, February 2010.
  85
  86     # Recursively construct a Nested Containment (NC) list by adding
  87     # a given feature to an existing NC list.
  88
  89     my ( $nc_list,   # NC list
  90          $feat,      # feature (AoA)
  91          $end,       # feature element with end position
  92          $index      # feature element with index position
  93        ) = @_;
  94
  95     # Returns nothing.
  96
  97     if ( $feat->[ $end ] <= $nc_list->[ -1 ]->[ $end ] ) # feature is nested.
  98     {
  99         if ( defined $nc_list->[ -1 ]->[ $index ] ) {   # sublist exists so recurse to this.
 100             nc_list_add( $nc_list->[ -1 ]->[ $index ], $feat, $end, $index );
 101         } else {
 102             $nc_list->[ -1 ]->[ $index ] = [ $feat ];   # creating a new sublist.
 103         }
 104     }
 105     else
 106     {
 107         push @{ $nc_list }, $feat;
 108     }
 109 }
 110
 111
 112 sub nc_list_lookup
 113 {
 114     # Martin A. Hansen, February 2010.
 115
 116     # Given a Nested Containment (NC) list use binary search to locate
 117     # the NC list containing a given search position. The index of the NC
 118     # list containing the search position is returned. Extremely fast.
 119
 120     my ( $nc_list,   # NC list
 121          $pos,       # search position
 122          $beg,       # feature element with begin position
 123          $end,       # feature element with end position
 124          $index,     # feature element with index position
 125        ) = @_;
 126
 127     # Returns an integer.
 128
 129     my ( $low, $high, $try );
 130
 131     $low  = 0;
 132     $high = scalar @{ $nc_list };
 133
 134     while ( $low < $high )
 135     {
 136         $try = int( ( $high + $low ) / 2 );
 137
 138         if ( $pos < $nc_list->[ $try ]->[ $beg ] ) {
 139             $high = $try;
 140         } elsif ( $nc_list->[ $try ]->[ $end ] < $pos ) {
 141             $low = $try + 1;
 142         } else {
 143             last;
 144         }
 145     }
 146
 147     return $try;
 148 }
 149
 150
 151 sub nc_list_count
 152 {
 153     # Martin A. Hansen, February 2010.
 154
 155     # Traverses a Nested Containment (NC) list recursively from a
 156     # given index begin to a given index end and counts all
 157     # features. The count is returned.
 158
 159     my ( $nc_list,     # NC list
 160          $index_beg,   # index begin
 161          $index_end,   # index end
 162          $index,       # feature element with index position
 163        ) = @_;
 164
 165     # Returns an integer.
 166
 167     my ( $i, $count );
 168
 169     for ( $i = $index_beg; $i <= $index_end; $i++ )
 170     {
 171         $count++;
 172
 173         if ( defined $nc_list->[ $i ]->[ $index ] ) {   # sublist exists so recurse to this.
 174             $count += nc_list_count( $nc_list->[ $i ]->[ $index ], 0, scalar @{ $nc_list->[ $i ]->[ $index ] } - 1, $index );
 175         }
 176     }
 177
 178     return $count;
 179 }
 180
 181
 182 sub nc_list_count_interval
 183 {
 184     # Martin A. Hansen, February 2010.
 185
 186     # Counts all features in a Nested Containment (NC) list within a
 187     # specified interval. The count is returned.
 188
 189     my ( $nc_list,   # NC list
 190          $int_beg,   # interval begin
 191          $int_end,   # interval end
 192          $beg,       # feature element with begin position
 193          $end,       # feature element with end position
 194          $index,     # feature element with index position
 195        ) = @_;
 196
 197     # Returns an integer.
 198
 199     my ( $index_beg, $index_end, $count );
 200
 201     $index_beg = nc_list_lookup( $nc_list, $int_beg, $beg, $end, $index );
 202     $index_end = nc_list_lookup( $nc_list, $int_end, $beg, $end, $index );
 203
 204     $count = nc_list_count( $nc_list, $index_beg, $index_end, $index );
 205
 206     return $count;
 207 }
 208
 209
 210 sub nc_list_get
 211 {
 212     # Martin A. Hansen, February 2010.
 213
 214     # Recursively retrieve all features from a Nested Containment (NC) list
 215     # from a specified index begin to a specified index end.
 216
 217     # WARNING: The NC list is distroyed because the sublists are stripped.
 218
 219     my ( $nc_list,     # NC list
 220          $index_beg,   # index begin
 221          $index_end,   # index end
 222          $index,       # feature element with index position
 223        ) = @_;
 224
 225     # Returns a list.
 226
 227     my ( $i, @features );
 228
 229     for ( $i = $index_beg; $i <= $index_end; $i++ )
 230     {
 231         push @features, $nc_list->[ $i ];
 232
 233         if ( defined $nc_list->[ $i ]->[ $index ] ) # sublist exists so recurse to this.
 234         {
 235             push @features, nc_list_get( $nc_list->[ $i ]->[ $index ], 0, scalar @{ $nc_list->[ $i ]->[ $index ] } - 1, $index );
 236
 237             delete $nc_list->[ $i ]->[ $index ];
 238         }
 239     }
 240
 241     return wantarray ? @features : \@features;
 242 }
 243
 244
 245 sub nc_list_get_interval
 246 {
 247     # Martin A. Hansen, February 2010.
 248
 249     # Retrieve all features from a Nested Containment (NC) list from within
 250     # a specified interval.
 251
 252     my ( $nc_list,   # NC list
 253          $int_beg,   # interval begin
 254          $int_end,   # interval end
 255          $beg,       # feature element with begin position
 256          $end,       # feature element with end position
 257          $index,     # feature element with index position
 258        ) = @_;
 259
 260     # Returns a list.
 261
 262     my ( $index_beg, $index_end, $features );
 263
 264     $index_beg = nc_list_lookup( $nc_list, $int_beg, $beg, $end, $index );
 265     $index_end = nc_list_lookup( $nc_list, $int_end, $beg, $end, $index );
 266
 267     $features = nc_list_get( $nc_list, $index_beg, $index_end, $index );
 268
 269     return wantarray ? @{ $features } : $features;
 270 }
 271
 272
 273 sub nc_list_search
 274 {
 275     # Martin A. Hansen, February 2010.
 276
 277     # Recursively search a Nested Containment (NC) list for features matching
 278     # a given REGEX.
 279
 280     my ( $nc_list,   # NC list
 281          $regex,     # regex to search for
 282          $index,     # feature element with index position
 283        ) = @_;
 284
 285     # Returns a list.
 286
 287     my ( $feature, @features );
 288
 289     foreach $feature ( @{ $nc_list } )
 290     {
 291         push @features, $feature if grep { $_ =~ /$regex/i if defined $_ } @{ $feature };
 292
 293         if ( defined $feature->[ $index ] ) # sublist exists so recurse to this.
 294         {
 295             push @features, nc_list_search( $feature->[ $index ], $regex, $index );
 296
 297             delete $feature->[ $index ];
 298         }
 299     }
 300
 301     return wantarray ? @features : \@features;
 302 }
 303
 304
 305 sub nc_list_store
 306 {
 307     # Martin A. Hansen, February 2010.
 308
 309     # Store a Nested Containment (NC) list to a
 310     # given file.
 311
 312     my ( $nc_list,   # NC list
 313          $file,      # path to file
 314        ) = @_;
 315
 316     # Returns nothing.
 317
 318     my ( $fh, $json );
 319
 320     $json = JSON::XS::encode_json( $nc_list );
 321
 322     $fh = Maasha::Filesys::file_write_open( $file );
 323
 324     print $fh $json;
 325
 326     close $fh;
 327 }
 328
 329
 330 sub nc_list_retrieve
 331 {
 332     # Martin A. Hansen, February 2010.
 333
 334     # Retrieves a Nested Containment (NC) list from a
 335     # JSON file.
 336
 337     my ( $file,   # path to JSON file
 338        ) = @_;
 339
 340     # Returns NC list.
 341
 342     my ( $fh, $json, $nc_list );
 343
 344     local $/ = undef;
 345
 346     $fh = Maasha::Filesys::file_read_open( $file );
 347
 348     $json = <$fh>;
 349
 350     close $fh;
 351
 352     $nc_list = JSON::XS::decode_json( $json );
 353
 354     return wantarray ? @{ $nc_list } : $nc_list;
 355 }
 356
 357
 358 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 359
 360
 361 1;
 362
 363 __END__
 364
 365     my $t0 = Time::HiRes::gettimeofday();
 366     my $t1 = Time::HiRes::gettimeofday(); print STDERR "Time: " . ( $t1 - $t0 ) . "\n";