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_sort
 158 {
 159     # Martin A. Hansen, November 2009.
 160
 161     # Sorts a KISS file on S_BEG and S_END
 162
 163     my ( $file,   # KISS file
 164        ) = @_;
 165
 166     # Returns nothing.
 167
 168     `sort -k 2,2n -k 3,3n $file > $file.sort`;
 169
 170     rename "$file.sort", $file;
 171 }
 172
 173
 174 sub kiss_index
 175 {
 176     # Martin A, Hansen, November 2009.
 177
 178     # Creates an index of a sorted KISS file.
 179
 180     my ( $file,   # KISS file
 181        ) = @_;
 182
 183     # Returns a hashref.
 184
 185     my ( $tree, $offset, $fh, $line, $beg );
 186
 187     $tree   = {};
 188     $offset = 0;
 189
 190     $fh = Maasha::Filesys::file_read_open( $file );
 191
 192     while ( $line = <$fh> )
 193     {
 194         ( undef, $beg ) = split "\t", $line, 3;
 195
 196         kiss_index_node_add( $tree, INDEX_LEVEL, INDEX_FACTOR, $beg, $offset );
 197
 198         $offset += length $line;
 199     }
 200
 201     close $fh;
 202
 203     kiss_index_store( "$file.index", $tree );
 204 }
 205
 206
 207 sub kiss_index_node_add
 208 {
 209     # Martin A, Hansen, November 2009.
 210
 211     # Recursive routine to add nodes to a tree.
 212
 213     my ( $node,
 214          $level,
 215          $factor,
 216          $beg,
 217          $offset,
 218          $sum,
 219        ) = @_;
 220
 221     my ( $bucket );
 222
 223     $sum  ||= 0;
 224     $bucket = int( $beg / $level );
 225
 226     if ( $level >= $factor )
 227     {
 228         $sum += $bucket * $level;
 229         $beg -= $bucket * $level;
 230
 231         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'COUNT' }++;
 232         # $node->{ 'CHILDREN' }->[ $bucket ]->{ 'LEVEL' }  = $level;
 233         # $node->{ 'CHILDREN' }->[ $bucket ]->{ 'BUCKET' } = $bucket;
 234         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'BEG' }    = $sum;
 235         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'END' }    = $sum + $level - 1;
 236         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'OFFSET' } = $offset if not defined $node->{ 'CHILDREN' }->[ $bucket ]->{ 'OFFSET' };
 237
 238         kiss_index_node_add( $node->{ 'CHILDREN' }->[ $bucket ], $level / $factor, $factor, $beg, $offset, $sum );
 239     }
 240 }
 241
 242
 243 sub kiss_index_offset
 244 {
 245     # Martin A. Hansen, November 2009.
 246
 247     # Given a KISS index and a begin position,
 248     # locate the offset closest to the begin position,
 249     # and return this.
 250
 251     my ( $index,    # KISS index
 252          $beg,      # begin position
 253          $level,    # index level  - OPTIONAL
 254          $factor,   # index factor - OPTIONAL
 255        ) = @_;
 256
 257     # Returns a number.
 258
 259     my ( $child, $offset );
 260
 261     $level  ||= INDEX_LEVEL;
 262     $factor ||= INDEX_FACTOR;
 263
 264     foreach $child ( @{ $index->{ 'CHILDREN' } } )
 265     {
 266         next if not defined $child;
 267
 268         if ( $child->{ 'BEG' } <= $beg and $beg <= $child->{ 'END' } )
 269         {
 270             if ( $level == $factor ) {
 271                 $offset = $child->{ 'OFFSET' };
 272             } else {
 273                 $offset = kiss_index_offset( $child, $beg, $level / $factor, $factor );
 274             }
 275         }
 276     }
 277
 278     return $offset;
 279 }
 280
 281
 282 sub kiss_index_count
 283 {
 284     # Martin A. Hansen, November 2009.
 285
 286     # Given a KISS index and a begin/end interval
 287     # sum the number of counts in that interval,
 288     # and return this.
 289
 290     my ( $index,   # KISS index
 291          $beg,     # Begin position
 292          $end,     # End position
 293          $level,   # index level  - OPTIONAL
 294          $factor,  # index factor - OPTIONAL
 295        ) = @_;
 296
 297     # Returns a number.
 298
 299     my ( $count, $child );
 300
 301     $level  ||= INDEX_LEVEL;
 302     $factor ||= INDEX_FACTOR;
 303     $count  ||= 0;
 304
 305     foreach $child ( @{ $index->{ 'CHILDREN' } } )
 306     {
 307         next if not defined $child;
 308
 309         if ( $level >= $factor )
 310         {
 311             if ( Maasha::Calc::overlap( $beg, $end, $child->{ 'BEG' }, $child->{ 'END' } ) )
 312             {
 313                 $count += $child->{ 'COUNT' } if $level == $factor;
 314                 $count += kiss_index_count( $child, $beg, $end, $level / $factor, $factor );
 315             }
 316         }
 317     }
 318
 319     return $count;
 320 }
 321
 322
 323 sub kiss_index_store
 324 {
 325     my ( $path,
 326          $index,
 327        ) = @_;
 328
 329     Maasha::Filesys::file_store( $path, $index );
 330 }
 331
 332
 333 sub kiss_index_retrieve
 334 {
 335     my ( $path,
 336        ) = @_;
 337
 338     my $index;
 339
 340     $index = Maasha::Filesys::file_retrieve( $path );
 341
 342     return wantarray ? @{ $index } : $index;
 343 }
 344
 345
 346 sub kiss_index_get_entries
 347 {
 348     my ( $file,
 349          $index,
 350          $beg,
 351          $end,
 352        ) = @_;
 353
 354     my ( $offset, $fh, $entry, @entries );
 355
 356     $offset = kiss_index_offset( $index, $beg );
 357
 358     $fh = Maasha::Filesys::file_read_open( $file );
 359
 360     sysseek( $fh, $offset, 0 );
 361
 362     while ( $entry = kiss_entry_get( $fh ) )
 363     {
 364         push @entries, $entry if $entry->{ 'S_END' } > $beg;
 365
 366         last if $entry->{ 'S_BEG' } > $end;
 367     }
 368
 369     close $fh;
 370
 371     return wantarray ? @entries : \@entries;
 372 }
 373
 374
 375 sub kiss_index_get_blocks
 376 {
 377     my ( $index,
 378          $beg,
 379          $end,
 380          $level,   # index level  - OPTIONAL
 381          $factor,  # index factor - OPTIONAL
 382          $size,
 383        ) = @_;
 384
 385     # Returns a list.
 386
 387     my ( $len, @blocks, $child );
 388
 389     $level  ||= INDEX_LEVEL;
 390     $factor ||= INDEX_FACTOR;
 391
 392     $size ||= 100;   # TODO: lazy list loading?
 393
 394 #    if ( not defined $size )
 395 #    {
 396 #        $len = $end - $beg + 1;
 397 #
 398 #        if ( $len > 100_000_000 ) {
 399 #            $size = 1_000_000;
 400 #        } elsif ( $len > 1_000_000 ) {
 401 #            $size = 10_000;
 402 #        } else {
 403 #            $size = 100;
 404 #        }
 405 #    }
 406
 407     if ( $level >= $size )
 408     {
 409         foreach $child ( @{ $index->{ 'CHILDREN' } } )
 410         {
 411             next if not defined $child;
 412
 413             if ( Maasha::Calc::overlap( $beg, $end, $child->{ 'BEG' }, $child->{ 'END' } ) )
 414             {
 415                 if ( $level == $size )
 416                 {
 417                     push @blocks, {
 418                         BEG   => $child->{ 'BEG' },
 419                         END   => $child->{ 'END' },
 420                         COUNT => $child->{ 'COUNT' },
 421                     };
 422                 }
 423
 424                 push @blocks, kiss_index_get_blocks( $child, $beg, $end, $level / $factor, $factor, $size );
 425             }
 426         }
 427     }
 428
 429     return wantarray ? @blocks : \@blocks;
 430 }
 431
 432
 433 sub kiss_align
 434 {
 435     # S.aur_COL       41      75      5_vnlh2ywXsN1/1 1       -       .       17:A>T,31:A>N   .       .       .       .
 436
 437     my ( $s_seqref,   # subject sequence reference
 438          $entry,      # KISS entry
 439        ) = @_;
 440
 441     # Returns a string.
 442
 443     my ( $align, $pos, $s_symbol, $q_symbol, $s_seq, $q_seq, $insertions );
 444
 445     $s_seq = substr ${ $s_seqref }, $entry->{ 'S_BEG' }, $entry->{ 'S_END' } - $entry->{ 'S_BEG' } + 1;
 446     $q_seq = $s_seq;
 447
 448     $insertions = 0;
 449
 450     foreach $align ( split /,/, $entry->{ 'ALIGN' } )
 451     {
 452         if ( $align =~ /(\d+):(.)>(.)/ )
 453         {
 454             $pos      = $1;
 455             $s_symbol = $2;
 456             $q_symbol = $3;
 457
 458             if ( $s_symbol eq '-' ) # insertion
 459             {
 460                 substr $s_seq, $pos + $insertions, 0, $s_symbol;
 461                 substr $q_seq, $pos + $insertions, 0, $q_symbol;
 462
 463                 $insertions++;
 464             }
 465             elsif ( $q_symbol eq '-' ) # deletion
 466             {
 467                 substr $q_seq, $pos + $insertions, 1, $q_symbol;
 468             }
 469             else # mismatch
 470             {
 471                 substr $q_seq, $pos + $insertions, 1, $q_symbol;
 472             }
 473         }
 474         else
 475         {
 476             Maasha::Common::error( qq(Bad alignment descriptor: "$align") );
 477         }
 478     }
 479
 480     Maasha::Align::align_print_pairwise( [ $entry->{ 'S_ID' }, $s_seq ], [ $entry->{ 'Q_ID' }, $q_seq ] );
 481 }
 482
 483
 484 sub kiss2biopiece
 485 {
 486     my ( $entry,   # KISS entry
 487        ) = @_;
 488
 489     return wantarray ? %{ $entry } : $entry;
 490 }
 491
 492
 493 sub biopiece2kiss
 494 {
 495     my ( $record,   # Biopiece record
 496        ) = @_;
 497
 498     $record->{ 'SCORE' }       ||= $record->{ 'E_VAL' } || ".";
 499     $record->{ 'HITS' }        ||= ".";
 500     $record->{ 'BLOCK_COUNT' } ||= ".";
 501     $record->{ 'BLOCK_BEGS' }  ||= ".";
 502     $record->{ 'BLOCK_LENS' }  ||= ".";
 503     $record->{ 'BLOCK_TYPE' }  ||= ".";
 504     $record->{ 'ALIGN' }       ||= $record->{ 'DESCRIPTOR' } || ".";
 505
 506     return wantarray ? %{ $record } : $record;
 507 }
 508
 509
 510 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 511
 512 1;
 513
 514 __END__
 515
 516 sub kiss_index
 517 {
 518     # Martin A, Hansen, October 2009.
 519
 520     # Creates an index of a sorted KISS file that
 521     # allowing the location of the byte position
 522     # from where records can be read given a
 523     # specific S_BEG position. The index consists of
 524     # triples: [ beg, end, bytepos ], where beg and
 525     # end denotes the interval where the next KISS
 526     # record begins at bytepos.
 527
 528     my ( $fh,   # filehandle to KISS file
 529        ) = @_;
 530
 531     # Returns a list.
 532
 533     my ( $line, @fields, $beg, $end, $pos, @index );
 534
 535     $beg = 0;
 536     $pos = 0;
 537
 538     while ( $line = <$fh> )
 539     {
 540         chomp $line;
 541
 542         @fields = split /\t/, $line, 3;
 543
 544         $end = $fields[ S_BEG ];
 545
 546         if ( $end == 0 )
 547         {
 548             push @index, [ $beg, $end, $pos ];
 549             $beg = 1;
 550         }
 551         elsif ( $end > $beg )
 552         {
 553             push @index, [ $beg, $end - 1, $pos ];
 554             $beg = $end;
 555         }
 556         elsif ( $end < $beg )
 557         {
 558             Maasha::Common::error( qq(KISS file not sorted: beg > end -> $beg > $end) );
 559         }
 560
 561         $pos += 1 + length $line;
 562     }
 563
 564     return wantarray ? @index : \@index;
 565 }
 566
 567
 568
 569 sub kiss_index_search
 570 {
 571     my ( $index,
 572          $num,
 573        ) = @_;
 574
 575     # Returns a number.
 576
 577     my ( $high, $low, $try );
 578
 579     $low  = 0;
 580     $high = scalar @{ $index };
 581
 582     while ( $low <= $high )
 583     {
 584         $try = int( ( $high + $low ) / 2 );
 585
 586         if ( $num < $index->[ $try ]->[ 0 ] ) {
 587             $high = $try;
 588         } elsif ( $num > $index->[ $try ]->[ 1 ] ) {
 589             $low = $try + 1;
 590         } else {
 591             return $index->[ $try ]->[ 2 ];
 592         }
 593     }
 594
 595     Maasha::Common::error( "Could not find number->$num in index" );
 596 }
 597
 598