code_perl/Maasha/KISS/IO.pm

   1 package Maasha::KISS::IO;
   2
   3 # Copyright (C) 2009 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 # Routines for parsing and emitting KISS records.
  26
  27 # http://code.google.com/p/biopieces/wiki/KissFormat
  28
  29
  30 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  31
  32
  33 use warnings;
  34 use strict;
  35 use Data::Dumper;
  36 use Maasha::Common;
  37 use Maasha::Filesys;
  38 use Maasha::Align;
  39 use Maasha::SQL;
  40 use vars qw( @ISA @EXPORT );
  41
  42 @ISA = qw( Exporter );
  43
  44 use constant {
  45     S_ID             => 0,
  46     S_BEG            => 1,
  47     S_END            => 2,
  48     Q_ID             => 3,
  49     SCORE            => 4,
  50     STRAND           => 5,
  51     HITS             => 6,
  52     ALIGN            => 7,
  53     BLOCK_COUNT      => 8,
  54     BLOCK_BEGS       => 9,
  55     BLOCK_LENS       => 10,
  56     BLOCK_TYPE       => 11,
  57     INDEX_BLOCK_SIZE => 100,
  58     INDEX_LEVEL      => 100_000_000,
  59     INDEX_FACTOR     => 100,
  60 };
  61
  62
  63 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  64
  65
  66 sub kiss_entry_get
  67 {
  68     my ( $fh,    # file handle
  69        ) = @_;
  70
  71     # Returns a hashref.
  72
  73     my ( $line, @fields, %entry );
  74
  75     while ( $line = <$fh> )
  76     {
  77         chomp $line;
  78
  79         next if $line =~ /^$|^#/;
  80
  81         @fields = split /\t/, $line;
  82
  83         Maasha::Common::error( qq(BAD kiss entry: $line) ) if not @fields == 12;
  84
  85         $entry{ 'S_ID' }        = $fields[ S_ID ];
  86         $entry{ 'S_BEG' }       = $fields[ S_BEG ];
  87         $entry{ 'S_END' }       = $fields[ S_END ];
  88         $entry{ 'Q_ID' }        = $fields[ Q_ID ];
  89         $entry{ 'SCORE' }       = $fields[ SCORE ];
  90         $entry{ 'STRAND' }      = $fields[ STRAND ];
  91         $entry{ 'HITS' }        = $fields[ HITS ];
  92         $entry{ 'ALIGN' }       = $fields[ ALIGN ];
  93         $entry{ 'BLOCK_COUNT' } = $fields[ BLOCK_COUNT ];
  94         $entry{ 'BLOCK_BEGS' }  = $fields[ BLOCK_BEGS ];
  95         $entry{ 'BLOCK_LENS' }  = $fields[ BLOCK_LENS ];
  96         $entry{ 'BLOCK_TYPE' }  = $fields[ BLOCK_TYPE ];
  97
  98         return wantarray ? %entry : \%entry;
  99     }
 100 }
 101
 102
 103 sub kiss_entry_put
 104 {
 105     my ( $entry,   # KISS entry to output
 106          $fh,      # file handle  -  OPTIONAL
 107        ) = @_;
 108
 109     # Returns nothing.
 110
 111     my ( @fields );
 112
 113     if ( defined $entry->{ 'S_ID' }  and
 114          defined $entry->{ 'S_BEG' } and
 115          defined $entry->{ 'S_END' }
 116        )
 117     {
 118         $fh ||= \*STDOUT;
 119
 120         $fields[ S_ID ]        = $entry->{ 'S_ID' };
 121         $fields[ S_BEG ]       = $entry->{ 'S_BEG' };
 122         $fields[ S_END ]       = $entry->{ 'S_END' };
 123         $fields[ Q_ID ]        = $entry->{ 'Q_ID' };
 124         $fields[ SCORE ]       = $entry->{ 'SCORE' };
 125         $fields[ STRAND ]      = $entry->{ 'STRAND' };
 126         $fields[ HITS ]        = $entry->{ 'HITS' };
 127         $fields[ ALIGN ]       = $entry->{ 'ALIGN' };
 128         $fields[ BLOCK_COUNT ] = $entry->{ 'BLOCK_COUNT' };
 129         $fields[ BLOCK_BEGS ]  = $entry->{ 'BLOCK_BEGS' };
 130         $fields[ BLOCK_LENS ]  = $entry->{ 'BLOCK_LENS' };
 131         $fields[ BLOCK_TYPE ]  = $entry->{ 'BLOCK_TYPE' };
 132
 133         print $fh join( "\t", @fields ), "\n";
 134     }
 135 }
 136
 137
 138 sub kiss_sql_get
 139 {
 140     my ( $dbh,     # Database handle
 141          $table,   # Table name
 142          $s_beg,   # Subject begin
 143          $s_end,   # Subject end
 144        ) = @_;
 145
 146     my ( $sql, $entries );
 147
 148     # $sql = "SELECT * FROM $table WHERE S_BEG >= $s_beg AND S_END <= $s_end ORDER BY S_BEG,S_END";
 149     $sql = "SELECT S_BEG,S_END,Q_ID,ALIGN FROM $table WHERE S_BEG >= $s_beg AND S_END <= $s_end";
 150
 151     $entries = Maasha::SQL::query_hashref_list( $dbh, $sql );
 152
 153     return wantarray ? @{ $entries } : $entries;
 154 }
 155
 156
 157 sub kiss_index
 158 {
 159     # Martin A, Hansen, November 2009.
 160
 161     # Creates an index of a sorted KISS file.
 162
 163     my ( $file,   # KISS file to index
 164        ) = @_;
 165
 166     # Returns a hashref.
 167
 168     my ( $tree, $offset, $fh, $line, $beg );
 169
 170     $tree   = {};
 171     $offset = 0;
 172
 173     $fh = Maasha::Filesys::file_read_open( $file );
 174
 175     while ( $line = <$fh> )
 176     {
 177         ( undef, $beg ) = split "\t", $line, 3;
 178
 179         kiss_index_node_add( $tree, INDEX_LEVEL, INDEX_FACTOR, $beg, $offset );
 180
 181         $offset += length $line;
 182     }
 183
 184     close $fh;
 185
 186     kiss_index_store( "$file.index", $tree );
 187 }
 188
 189
 190 sub kiss_index_node_add
 191 {
 192     # Martin A, Hansen, November 2009.
 193
 194     # Recursive routine to add nodes to a tree.
 195
 196     my ( $node,
 197          $level,
 198          $factor,
 199          $beg,
 200          $offset,
 201          $sum,
 202        ) = @_;
 203
 204     my ( $bucket );
 205
 206     $sum  ||= 0;
 207     $bucket = int( $beg / $level );
 208
 209     if ( $level >= $factor )
 210     {
 211         $sum += $bucket * $level;
 212         $beg -= $bucket * $level;
 213
 214         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'COUNT' }++;
 215         # $node->{ 'CHILDREN' }->[ $bucket ]->{ 'LEVEL' }  = $level;
 216         # $node->{ 'CHILDREN' }->[ $bucket ]->{ 'BUCKET' } = $bucket;
 217         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'BEG' }    = $sum;
 218         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'END' }    = $sum + $level - 1;
 219         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'OFFSET' } = $offset if not defined $node->{ 'CHILDREN' }->[ $bucket ]->{ 'OFFSET' };
 220
 221         kiss_index_node_add( $node->{ 'CHILDREN' }->[ $bucket ], $level / $factor, $factor, $beg, $offset, $sum );
 222     }
 223 }
 224
 225
 226 sub kiss_index_offset
 227 {
 228     # Martin A. Hansen, November 2009.
 229
 230     # Given a KISS index and a begin position,
 231     # locate the offset closest to the begin position,
 232     # and return this.
 233
 234     my ( $index,    # KISS index
 235          $beg,      # begin position
 236          $level,    # index level  - OPTIONAL
 237          $factor,   # index factor - OPTIONAL
 238        ) = @_;
 239
 240     # Returns a number.
 241
 242     my ( $child, $offset );
 243
 244     $level  ||= INDEX_LEVEL;
 245     $factor ||= INDEX_FACTOR;
 246
 247     foreach $child ( @{ $index->{ 'CHILDREN' } } )
 248     {
 249         next if not defined $child;
 250
 251         if ( $child->{ 'BEG' } <= $beg and $beg <= $child->{ 'END' } )
 252         {
 253             if ( $level == $factor ) {
 254                 $offset = $child->{ 'OFFSET' };
 255             } else {
 256                 $offset = kiss_index_offset( $child, $beg, $level / $factor, $factor );
 257             }
 258         }
 259     }
 260
 261     return $offset;
 262 }
 263
 264
 265 sub kiss_index_count
 266 {
 267     # Martin A. Hansen, November 2009.
 268
 269     # Given a KISS index and a begin/end interval
 270     # sum the number of counts in that interval,
 271     # and return this.
 272
 273     my ( $index,   # KISS index
 274          $beg,     # Begin position
 275          $end,     # End position
 276          $level,   # index level  - OPTIONAL
 277          $factor,  # index factor - OPTIONAL
 278        ) = @_;
 279
 280     # Returns a number.
 281
 282     my ( $count, $child );
 283
 284     $level  ||= INDEX_LEVEL;
 285     $factor ||= INDEX_FACTOR;
 286     $count  ||= 0;
 287
 288     foreach $child ( @{ $index->{ 'CHILDREN' } } )
 289     {
 290         next if not defined $child;
 291
 292         if ( $level >= $factor )
 293         {
 294             if ( Maasha::Calc::overlap( $beg, $end, $child->{ 'BEG' }, $child->{ 'END' } ) )
 295             {
 296                 $count += $child->{ 'COUNT' } if $level == $factor;
 297                 $count += kiss_index_count( $child, $beg, $end, $level / $factor, $factor );
 298             }
 299         }
 300     }
 301
 302     return $count;
 303 }
 304
 305
 306 sub kiss_index_store
 307 {
 308     my ( $path,
 309          $index,
 310        ) = @_;
 311
 312     Maasha::Filesys::file_store( $path, $index );
 313 }
 314
 315
 316 sub kiss_index_retrieve
 317 {
 318     my ( $path,
 319        ) = @_;
 320
 321     my $index;
 322
 323     $index = Maasha::Filesys::file_retrieve( $path );
 324
 325     return wantarray ? @{ $index } : $index;
 326 }
 327
 328
 329 sub kiss_index_get_entries
 330 {
 331     my ( $file,
 332          $index,
 333          $beg,
 334          $end,
 335        ) = @_;
 336
 337     my ( $offset, $fh, $entry, @entries );
 338
 339     $offset = kiss_index_offset( $index, $beg );
 340
 341     $fh = Maasha::Filesys::file_read_open( $file );
 342
 343     sysseek( $fh, $offset, 0 );
 344
 345     while ( $entry = kiss_entry_get( $fh ) )
 346     {
 347         push @entries, $entry if $entry->{ 'S_END' } > $beg;
 348
 349         last if $entry->{ 'S_BEG' } > $end;
 350     }
 351
 352     close $fh;
 353
 354     return wantarray ? @entries : \@entries;
 355 }
 356
 357
 358 sub kiss_index_get_blocks
 359 {
 360     my ( $index,
 361          $beg,
 362          $end,
 363          $level,   # index level  - OPTIONAL
 364          $factor,  # index factor - OPTIONAL
 365          $size,
 366        ) = @_;
 367
 368     # Returns a list.
 369
 370     my ( $len, @blocks, $child );
 371
 372     $level  ||= INDEX_LEVEL;
 373     $factor ||= INDEX_FACTOR;
 374
 375     $size ||= 100;   # TODO: lazy list loading?
 376
 377 #    if ( not defined $size )
 378 #    {
 379 #        $len = $end - $beg + 1;
 380 #
 381 #        if ( $len > 100_000_000 ) {
 382 #            $size = 1_000_000;
 383 #        } elsif ( $len > 1_000_000 ) {
 384 #            $size = 10_000;
 385 #        } else {
 386 #            $size = 100;
 387 #        }
 388 #    }
 389
 390     if ( $level >= $size )
 391     {
 392         foreach $child ( @{ $index->{ 'CHILDREN' } } )
 393         {
 394             next if not defined $child;
 395
 396             if ( Maasha::Calc::overlap( $beg, $end, $child->{ 'BEG' }, $child->{ 'END' } ) )
 397             {
 398                 if ( $level == $size )
 399                 {
 400                     push @blocks, {
 401                         BEG   => $child->{ 'BEG' },
 402                         END   => $child->{ 'END' },
 403                         COUNT => $child->{ 'COUNT' },
 404                     };
 405                 }
 406
 407                 push @blocks, kiss_index_get_blocks( $child, $beg, $end, $level / $factor, $factor, $size );
 408             }
 409         }
 410     }
 411
 412     return wantarray ? @blocks : \@blocks;
 413 }
 414
 415
 416 sub kiss_align
 417 {
 418     # S.aur_COL       41      75      5_vnlh2ywXsN1/1 1       -       .       17:A>T,31:A>N   .       .       .       .
 419
 420     my ( $s_seqref,   # subject sequence reference
 421          $entry,      # KISS entry
 422        ) = @_;
 423
 424     # Returns a string.
 425
 426     my ( $align, $pos, $s_symbol, $q_symbol, $s_seq, $q_seq, $insertions );
 427
 428     $s_seq = substr ${ $s_seqref }, $entry->{ 'S_BEG' }, $entry->{ 'S_END' } - $entry->{ 'S_BEG' } + 1;
 429     $q_seq = $s_seq;
 430
 431     $insertions = 0;
 432
 433     foreach $align ( split /,/, $entry->{ 'ALIGN' } )
 434     {
 435         if ( $align =~ /(\d+):(.)>(.)/ )
 436         {
 437             $pos      = $1;
 438             $s_symbol = $2;
 439             $q_symbol = $3;
 440
 441             if ( $s_symbol eq '-' ) # insertion
 442             {
 443                 substr $s_seq, $pos + $insertions, 0, $s_symbol;
 444                 substr $q_seq, $pos + $insertions, 0, $q_symbol;
 445
 446                 $insertions++;
 447             }
 448             elsif ( $q_symbol eq '-' ) # deletion
 449             {
 450                 substr $q_seq, $pos + $insertions, 1, $q_symbol;
 451             }
 452             else # mismatch
 453             {
 454                 substr $q_seq, $pos + $insertions, 1, $q_symbol;
 455             }
 456         }
 457         else
 458         {
 459             Maasha::Common::error( qq(Bad alignment descriptor: "$align") );
 460         }
 461     }
 462
 463     Maasha::Align::align_print_pairwise( [ $entry->{ 'S_ID' }, $s_seq ], [ $entry->{ 'Q_ID' }, $q_seq ] );
 464 }
 465
 466
 467 sub kiss2biopiece
 468 {
 469     my ( $entry,   # KISS entry
 470        ) = @_;
 471
 472     return wantarray ? %{ $entry } : $entry;
 473 }
 474
 475
 476 sub biopiece2kiss
 477 {
 478     my ( $record,   # Biopiece record
 479        ) = @_;
 480
 481     $record->{ 'SCORE' }       ||= $record->{ 'E_VAL' } || ".";
 482     $record->{ 'HITS' }        ||= ".";
 483     $record->{ 'BLOCK_COUNT' } ||= ".";
 484     $record->{ 'BLOCK_BEGS' }  ||= ".";
 485     $record->{ 'BLOCK_LENS' }  ||= ".";
 486     $record->{ 'BLOCK_TYPE' }  ||= ".";
 487     $record->{ 'ALIGN' }       ||= $record->{ 'DESCRIPTOR' } || ".";
 488
 489     return wantarray ? %{ $record } : $record;
 490 }
 491
 492
 493 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 494
 495 1;
 496
 497 __END__
 498
 499 sub kiss_index
 500 {
 501     # Martin A, Hansen, October 2009.
 502
 503     # Creates an index of a sorted KISS file that
 504     # allowing the location of the byte position
 505     # from where records can be read given a
 506     # specific S_BEG position. The index consists of
 507     # triples: [ beg, end, bytepos ], where beg and
 508     # end denotes the interval where the next KISS
 509     # record begins at bytepos.
 510
 511     my ( $fh,   # filehandle to KISS file
 512        ) = @_;
 513
 514     # Returns a list.
 515
 516     my ( $line, @fields, $beg, $end, $pos, @index );
 517
 518     $beg = 0;
 519     $pos = 0;
 520
 521     while ( $line = <$fh> )
 522     {
 523         chomp $line;
 524
 525         @fields = split /\t/, $line, 3;
 526
 527         $end = $fields[ S_BEG ];
 528
 529         if ( $end == 0 )
 530         {
 531             push @index, [ $beg, $end, $pos ];
 532             $beg = 1;
 533         }
 534         elsif ( $end > $beg )
 535         {
 536             push @index, [ $beg, $end - 1, $pos ];
 537             $beg = $end;
 538         }
 539         elsif ( $end < $beg )
 540         {
 541             Maasha::Common::error( qq(KISS file not sorted: beg > end -> $beg > $end) );
 542         }
 543
 544         $pos += 1 + length $line;
 545     }
 546
 547     return wantarray ? @index : \@index;
 548 }
 549
 550
 551
 552 sub kiss_index_search
 553 {
 554     my ( $index,
 555          $num,
 556        ) = @_;
 557
 558     # Returns a number.
 559
 560     my ( $high, $low, $try );
 561
 562     $low  = 0;
 563     $high = scalar @{ $index };
 564
 565     while ( $low <= $high )
 566     {
 567         $try = int( ( $high + $low ) / 2 );
 568
 569         if ( $num < $index->[ $try ]->[ 0 ] ) {
 570             $high = $try;
 571         } elsif ( $num > $index->[ $try ]->[ 1 ] ) {
 572             $low = $try + 1;
 573         } else {
 574             return $index->[ $try ]->[ 2 ];
 575         }
 576     }
 577
 578     Maasha::Common::error( "Could not find number->$num in index" );
 579 }
 580
 581