code_perl/Maasha/KISS.pm

   1 package Maasha::KISS;
   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 vars qw( @ISA @EXPORT );
  40
  41 @ISA = qw( Exporter );
  42
  43 use constant {
  44     S_ID             => 0,
  45     S_BEG            => 1,
  46     S_END            => 2,
  47     Q_ID             => 3,
  48     SCORE            => 4,
  49     STRAND           => 5,
  50     HITS             => 6,
  51     ALIGN            => 7,
  52     BLOCK_COUNT      => 8,
  53     BLOCK_BEGS       => 9,
  54     BLOCK_LENS       => 10,
  55     BLOCK_TYPE       => 11,
  56     INDEX_BLOCK_SIZE => 100,
  57     INDEX_LEVEL      => 100_000_000,
  58     INDEX_FACTOR     => 100,
  59 };
  60
  61
  62 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  63
  64
  65 sub kiss_entry_get
  66 {
  67     # Martin A. Hansen, November 2009.
  68
  69     # Gets the next KISS entry from a file handle.
  70
  71     my ( $fh,    # file handle
  72        ) = @_;
  73
  74     # Returns a hashref.
  75
  76     my ( $line, @fields, %entry );
  77
  78     while ( $line = <$fh> )
  79     {
  80         chomp $line;
  81
  82         next if $line =~ /^$|^#/;
  83
  84         @fields = split /\t/, $line;
  85
  86         Maasha::Common::error( qq(BAD kiss entry: $line) ) if not @fields == 12;
  87
  88         $entry{ 'S_ID' }        = $fields[ S_ID ];
  89         $entry{ 'S_BEG' }       = $fields[ S_BEG ];
  90         $entry{ 'S_END' }       = $fields[ S_END ];
  91         $entry{ 'Q_ID' }        = $fields[ Q_ID ];
  92         $entry{ 'SCORE' }       = $fields[ SCORE ];
  93         $entry{ 'STRAND' }      = $fields[ STRAND ];
  94         $entry{ 'HITS' }        = $fields[ HITS ];
  95         $entry{ 'ALIGN' }       = $fields[ ALIGN ];
  96         $entry{ 'BLOCK_COUNT' } = $fields[ BLOCK_COUNT ];
  97         $entry{ 'BLOCK_BEGS' }  = $fields[ BLOCK_BEGS ];
  98         $entry{ 'BLOCK_LENS' }  = $fields[ BLOCK_LENS ];
  99         $entry{ 'BLOCK_TYPE' }  = $fields[ BLOCK_TYPE ];
 100
 101         return wantarray ? %entry : \%entry;
 102     }
 103 }
 104
 105
 106 sub kiss_entry_parse
 107 {
 108     # Martin A. Hansen, December 2009.
 109
 110     # Parses a line with a KISS entry.
 111
 112     my ( $line,   #  KISS line to parse
 113        ) = @_;
 114
 115     # Returns a hash
 116
 117     my ( @fields, %entry );
 118
 119     @fields = split /\t/, $line;
 120
 121     Maasha::Common::error( qq(BAD kiss entry: $line) ) if not @fields == 12;
 122
 123     $entry{ 'S_ID' }        = $fields[ S_ID ];
 124     $entry{ 'S_BEG' }       = $fields[ S_BEG ];
 125     $entry{ 'S_END' }       = $fields[ S_END ];
 126     $entry{ 'Q_ID' }        = $fields[ Q_ID ];
 127     $entry{ 'SCORE' }       = $fields[ SCORE ];
 128     $entry{ 'STRAND' }      = $fields[ STRAND ];
 129     $entry{ 'HITS' }        = $fields[ HITS ];
 130     $entry{ 'ALIGN' }       = $fields[ ALIGN ];
 131     $entry{ 'BLOCK_COUNT' } = $fields[ BLOCK_COUNT ];
 132     $entry{ 'BLOCK_BEGS' }  = $fields[ BLOCK_BEGS ];
 133     $entry{ 'BLOCK_LENS' }  = $fields[ BLOCK_LENS ];
 134     $entry{ 'BLOCK_TYPE' }  = $fields[ BLOCK_TYPE ];
 135
 136     return wantarray ? %entry : \%entry;
 137 }
 138
 139
 140 sub kiss_entry_put
 141 {
 142     # Martin A. Hansen, November 2009.
 143
 144     # Outputs a KISS record to a given filehandle
 145     # or STDOUT if no filehandle is given.
 146
 147     my ( $entry,   # KISS entry to output
 148          $fh,      # file handle  -  OPTIONAL
 149        ) = @_;
 150
 151     # Returns nothing.
 152
 153     my ( @fields );
 154
 155     if ( defined $entry->{ 'S_ID' }  and
 156          defined $entry->{ 'S_BEG' } and
 157          defined $entry->{ 'S_END' }
 158        )
 159     {
 160         Maasha::Common::error( qq(Bad S_BEG value: $entry->{ 'S_BEG' } < 0 ) ) if $entry->{ 'S_BEG' } < 0;
 161         Maasha::Common::error( qq(Bad S_END value: $entry->{ 'S_END' } < $entry->{ 'S_BEG' }) ) if $entry->{ 'S_END' } < $entry->{ 'S_BEG' };
 162
 163         $fh ||= \*STDOUT;
 164
 165         $fields[ S_ID ]        = $entry->{ 'S_ID' };
 166         $fields[ S_BEG ]       = $entry->{ 'S_BEG' };
 167         $fields[ S_END ]       = $entry->{ 'S_END' };
 168         $fields[ Q_ID ]        = $entry->{ 'Q_ID' }        || ".";
 169         $fields[ SCORE ]       = $entry->{ 'SCORE' }       || ".";
 170         $fields[ STRAND ]      = $entry->{ 'STRAND' }      || ".";
 171         $fields[ HITS ]        = $entry->{ 'HITS' }        || ".";
 172         $fields[ ALIGN ]       = $entry->{ 'ALIGN' }       || ".";
 173         $fields[ BLOCK_COUNT ] = $entry->{ 'BLOCK_COUNT' } || ".";
 174         $fields[ BLOCK_BEGS ]  = $entry->{ 'BLOCK_BEGS' }  || ".";
 175         $fields[ BLOCK_LENS ]  = $entry->{ 'BLOCK_LENS' }  || ".";
 176         $fields[ BLOCK_TYPE ]  = $entry->{ 'BLOCK_TYPE' }  || ".";
 177
 178         print $fh join( "\t", @fields ), "\n";
 179     }
 180 }
 181
 182
 183 sub kiss_sort
 184 {
 185     # Martin A. Hansen, November 2009.
 186
 187     # Sorts a KISS file on S_BEG and S_END
 188
 189     my ( $file,   # KISS file
 190        ) = @_;
 191
 192     # Returns nothing.
 193
 194     `sort -k 2,2n -k 3,3n $file > $file.sort`;
 195
 196     rename "$file.sort", $file;
 197 }
 198
 199
 200 sub kiss_intersect
 201 {
 202     # Martin A. Hansen, December 2009.
 203
 204     # Given filehandles to two different unsorted KISS files
 205     # intersect the entries so that all entries from file1 that
 206     # overlap entries in file2 are returned - unless the inverse flag
 207     # is given in which case entreis from file1 that does not
 208     # overlap any entries in file2 are returned.
 209
 210     my ( $fh1,       # filehandle to file1
 211          $fh2,       # filehandle to file2
 212          $inverse,   # flag indicating inverse matching - OPTIONAL
 213        ) = @_;
 214
 215     # Returns a list
 216
 217     my ( $entry, %lookup, $pos, $overlap, @entries );
 218
 219     while ( $entry = kiss_entry_get( $fh2 ) ) {
 220         map { $lookup{ $_ } = 1 } ( $entry->{ 'S_BEG' } .. $entry->{ 'S_END' } );
 221     }
 222
 223     while ( $entry = kiss_entry_get( $fh1 ) )
 224     {
 225         $overlap = 0;
 226
 227         foreach $pos ( $entry->{ 'S_BEG' } .. $entry->{ 'S_END' } )
 228         {
 229             if ( exists $lookup{ $pos } )
 230             {
 231                 $overlap = 1;
 232
 233                 last;
 234             }
 235         }
 236
 237         if ( $overlap and not $inverse ) {
 238             push @entries, $entry;
 239         } elsif ( not $overlap and $inverse ) {
 240             push @entries, $entry;
 241         }
 242     }
 243
 244     return wantarray ? @entries : \@entries;
 245 }
 246
 247
 248 sub kiss_index
 249 {
 250     # Martin A, Hansen, November 2009.
 251
 252     # Creates an index of a sorted KISS file.
 253
 254     my ( $file,   # KISS file
 255        ) = @_;
 256
 257     # Returns a hashref.
 258
 259     my ( $tree, $offset, $fh, $line, $beg );
 260
 261     $tree   = {};
 262     $offset = 0;
 263
 264     $fh = Maasha::Filesys::file_read_open( $file );
 265
 266     while ( $line = <$fh> )
 267     {
 268         ( undef, $beg ) = split "\t", $line, 3;
 269
 270         kiss_index_node_add( $tree, INDEX_LEVEL, INDEX_FACTOR, $beg, $offset );
 271
 272         $offset += length $line;
 273     }
 274
 275     close $fh;
 276
 277     kiss_index_store( "$file.index", $tree );
 278 }
 279
 280
 281 sub kiss_index_node_add
 282 {
 283     # Martin A, Hansen, November 2009.
 284
 285     # Recursive routine to add nodes to a tree.
 286
 287     my ( $node,
 288          $level,
 289          $factor,
 290          $beg,
 291          $offset,
 292          $sum,
 293        ) = @_;
 294
 295     my ( $bucket );
 296
 297     $sum  ||= 0;
 298     $bucket = int( $beg / $level );
 299
 300     if ( $level >= $factor )
 301     {
 302         $sum += $bucket * $level;
 303         $beg -= $bucket * $level;
 304
 305         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'COUNT' }++;
 306         # $node->{ 'CHILDREN' }->[ $bucket ]->{ 'LEVEL' }  = $level;
 307         # $node->{ 'CHILDREN' }->[ $bucket ]->{ 'BUCKET' } = $bucket;
 308         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'BEG' }    = $sum;
 309         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'END' }    = $sum + $level - 1;
 310         $node->{ 'CHILDREN' }->[ $bucket ]->{ 'OFFSET' } = $offset if not defined $node->{ 'CHILDREN' }->[ $bucket ]->{ 'OFFSET' };
 311
 312         kiss_index_node_add( $node->{ 'CHILDREN' }->[ $bucket ], $level / $factor, $factor, $beg, $offset, $sum );
 313     }
 314 }
 315
 316
 317 sub kiss_index_offset
 318 {
 319     # Martin A. Hansen, November 2009.
 320
 321     # Given a KISS index and a begin position,
 322     # locate the offset closest to the begin position,
 323     # and return this.
 324
 325     my ( $index,    # KISS index
 326          $beg,      # begin position
 327          $level,    # index level  - OPTIONAL
 328          $factor,   # index factor - OPTIONAL
 329        ) = @_;
 330
 331     # Returns a number.
 332
 333     my ( $child, $offset );
 334
 335     $level  ||= INDEX_LEVEL;
 336     $factor ||= INDEX_FACTOR;
 337
 338     foreach $child ( @{ $index->{ 'CHILDREN' } } )
 339     {
 340         next if not defined $child;
 341
 342         if ( $child->{ 'BEG' } <= $beg and $beg <= $child->{ 'END' } )
 343         {
 344             if ( $level == $factor ) {
 345                 $offset = $child->{ 'OFFSET' };
 346             } else {
 347                 $offset = kiss_index_offset( $child, $beg, $level / $factor, $factor );
 348             }
 349         }
 350     }
 351
 352     # Maasha::Common::error( "No offset" ) if not defined $offset;   # FIXME
 353
 354     return $offset;
 355 }
 356
 357
 358 sub kiss_index_count
 359 {
 360     # Martin A. Hansen, November 2009.
 361
 362     # Given a KISS index and a begin/end interval
 363     # sum the number of counts in that interval,
 364     # and return this.
 365
 366     my ( $index,   # KISS index
 367          $beg,     # Begin position
 368          $end,     # End position
 369          $level,   # index level  - OPTIONAL
 370          $factor,  # index factor - OPTIONAL
 371        ) = @_;
 372
 373     # Returns a number.
 374
 375     my ( $count, $child );
 376
 377     $level  ||= INDEX_LEVEL;
 378     $factor ||= INDEX_FACTOR;
 379     $count  ||= 0;
 380
 381     foreach $child ( @{ $index->{ 'CHILDREN' } } )
 382     {
 383         next if not defined $child;
 384
 385         if ( $level >= $factor )
 386         {
 387             if ( Maasha::Calc::overlap( $beg, $end, $child->{ 'BEG' }, $child->{ 'END' } ) )
 388             {
 389                 $count += $child->{ 'COUNT' } if $level == $factor;
 390                 $count += kiss_index_count( $child, $beg, $end, $level / $factor, $factor );
 391             }
 392         }
 393     }
 394
 395     return $count;
 396 }
 397
 398
 399 sub kiss_index_store
 400 {
 401     # Martin A. Hansen, November 2009.
 402
 403     # Stores a KISS index to a file.
 404
 405     my ( $path,    # path to KISS index
 406          $index,   # KISS index
 407        ) = @_;
 408
 409     # Returns nothing.
 410
 411     Maasha::Filesys::file_store( $path, $index );
 412 }
 413
 414
 415 sub kiss_index_retrieve
 416 {
 417     # Martin A. Hansen, November 2009.
 418
 419     # Retrieves a KISS index from a file.
 420
 421     my ( $path,   # Path to KISS index
 422        ) = @_;
 423
 424     # Returns a data structure.
 425
 426     my ( $index );
 427
 428     $index = Maasha::Filesys::file_retrieve( $path );
 429
 430     return wantarray ? @{ $index } : $index;
 431 }
 432
 433
 434 sub kiss_index_get_entries
 435 {
 436     # Martin A. Hansen, November 2009.
 437
 438     # Given a path to a KISS file and a KISS index
 439     # along with a beg/end interval, locate all entries
 440     # in that interval and return those.
 441
 442     my ( $file,    # path to KISS file
 443          $index,   # KISS index
 444          $beg,     # interval begin
 445          $end,     # interval end
 446        ) = @_;
 447
 448     # Returns a list.
 449
 450     my ( $offset, $fh, $entry, @entries );
 451
 452     $offset = kiss_index_offset( $index, $beg );
 453
 454     $fh = Maasha::Filesys::file_read_open( $file );
 455
 456     sysseek( $fh, $offset, 0 );
 457
 458     while ( $entry = kiss_entry_get( $fh ) )
 459     {
 460         push @entries, $entry if $entry->{ 'S_END' } > $beg;
 461
 462         last if $entry->{ 'S_BEG' } > $end;
 463     }
 464
 465     close $fh;
 466
 467     return wantarray ? @entries : \@entries;
 468 }
 469
 470
 471 sub kiss_index_get_blocks
 472 {
 473     # Martin A. Hansen, November 2009.
 474
 475     # Given a KISS index recursively traverse
 476     # this into the appropriate node size determined
 477     # by the size of the given beg/end interval.
 478     # Blocks consisting of hashes of BEG/END/COUNT
 479     # are returned from the requested node size.
 480
 481     my ( $index,   # KISS index node
 482          $beg,     # interval begin
 483          $end,     # interval end
 484          $level,   # index level  - OPTIONAL
 485          $factor,  # index factor - OPTIONAL
 486          $size,    # requested node size
 487        ) = @_;
 488
 489     # Returns a list.
 490
 491     my ( $len, @blocks, $child );
 492
 493     $level  ||= INDEX_LEVEL;
 494     $factor ||= INDEX_FACTOR;
 495
 496     $size ||= 100;   # TODO: lazy list loading?
 497
 498 #    if ( not defined $size )
 499 #    {
 500 #        $len = $end - $beg + 1;
 501 #
 502 #        if ( $len > 100_000_000 ) {
 503 #            $size = 1_000_000;
 504 #        } elsif ( $len > 1_000_000 ) {
 505 #            $size = 10_000;
 506 #        } else {
 507 #            $size = 100;
 508 #        }
 509 #    }
 510
 511     if ( $level >= $size )
 512     {
 513         foreach $child ( @{ $index->{ 'CHILDREN' } } )
 514         {
 515             next if not defined $child;
 516
 517             if ( Maasha::Calc::overlap( $beg, $end, $child->{ 'BEG' }, $child->{ 'END' } ) )
 518             {
 519                 if ( $level == $size )
 520                 {
 521                     push @blocks, {
 522                         BEG   => $child->{ 'BEG' },
 523                         END   => $child->{ 'END' },
 524                         COUNT => $child->{ 'COUNT' },
 525                     };
 526                 }
 527
 528                 push @blocks, kiss_index_get_blocks( $child, $beg, $end, $level / $factor, $factor, $size );
 529             }
 530         }
 531     }
 532
 533     return wantarray ? @blocks : \@blocks;
 534 }
 535
 536
 537 sub kiss_align_enc
 538 {
 539     # Martin A. Hansen, November 2009.
 540
 541     # Encodes alignment descriptors for two
 542     # aligned sequences.
 543
 544     my ( $s_seq,    # Subject sequence reference
 545          $q_seq,    # Query sequence reference
 546          $offset,   # alternate offset - OPTIONAL
 547        ) = @_;
 548
 549     # Returns a list
 550
 551     my ( $i, $s, $q, @align );
 552
 553     Maasha::Common::error( "Sequence lengths don't match" ) if length ${ $s_seq } != length ${ $q_seq };
 554
 555     $offset ||= 0;
 556
 557     for ( $i = 0; $i < length ${ $s_seq }; $i++ )
 558     {
 559         $s = uc substr ${ $s_seq }, $i, 1;
 560         $q = uc substr ${ $q_seq }, $i, 1;
 561
 562         if ( $s eq '-' and $q eq '-' )
 563         {
 564             # do nothing
 565         }
 566         elsif ( $s eq $q )
 567         {
 568             $offset++;
 569         }
 570         else
 571         {
 572             push @align, "$offset:$s>$q";
 573
 574             $offset++ if not $s eq '-';
 575         }
 576     }
 577
 578     return wantarray ? @align : \@align;
 579 }
 580
 581
 582 sub kiss_align_dec
 583 {
 584     # Martin A. Hansen, November 2009.
 585
 586     # Test routine of correct resolve of ALIGN descriptors.
 587
 588     # S.aur_COL       41      75      5_vnlh2ywXsN1/1 1       -       .       17:A>T,31:A>N   .       .       .       .
 589
 590     my ( $s_seqref,   # subject sequence reference
 591          $entry,      # KISS entry
 592        ) = @_;
 593
 594     # Returns a string.
 595
 596     my ( $align, $pos, $s_symbol, $q_symbol, $s_seq, $q_seq, $insertions );
 597
 598     $s_seq = substr ${ $s_seqref }, $entry->{ 'S_BEG' }, $entry->{ 'S_END' } - $entry->{ 'S_BEG' } + 1;
 599     $q_seq = $s_seq;
 600
 601     $insertions = 0;
 602
 603     foreach $align ( split /,/, $entry->{ 'ALIGN' } )
 604     {
 605         if ( $align =~ /(\d+):(.)>(.)/ )
 606         {
 607             $pos      = $1;
 608             $s_symbol = $2;
 609             $q_symbol = $3;
 610
 611             if ( $s_symbol eq '-' ) # insertion
 612             {
 613                 substr $s_seq, $pos + $insertions, 0, $s_symbol;
 614                 substr $q_seq, $pos + $insertions, 0, $q_symbol;
 615
 616                 $insertions++;
 617             }
 618             elsif ( $q_symbol eq '-' ) # deletion
 619             {
 620                 substr $q_seq, $pos + $insertions, 1, $q_symbol;
 621             }
 622             else # mismatch
 623             {
 624                 substr $q_seq, $pos + $insertions, 1, $q_symbol;
 625             }
 626         }
 627         else
 628         {
 629             Maasha::Common::error( qq(Bad alignment descriptor: "$align") );
 630         }
 631     }
 632
 633     Maasha::Align::align_print_pairwise( [ $entry->{ 'S_ID' }, $s_seq ], [ $entry->{ 'Q_ID' }, $q_seq ] );
 634 }
 635
 636
 637 sub kiss2biopiece
 638 {
 639     # Martin A. Hansen, November 2009.
 640
 641     # Converts a KISS entry to a Biopiece record.
 642     # TODO: Consistency checking
 643
 644     my ( $entry,   # KISS entry
 645        ) = @_;
 646
 647     # Returns a hashref
 648
 649     return wantarray ? %{ $entry } : $entry;
 650 }
 651
 652
 653 sub biopiece2kiss
 654 {
 655     # Martin A. Hansen, November 2009.
 656
 657     # Converts a Biopiece record to a KISS entry.
 658
 659     my ( $record,   # Biopiece record
 660        ) = @_;
 661
 662     # Returns a hashref
 663
 664     $record->{ 'SCORE' }       ||= $record->{ 'E_VAL' } || ".";
 665     $record->{ 'HITS' }        ||= ".";
 666     $record->{ 'BLOCK_COUNT' } ||= ".";
 667     $record->{ 'BLOCK_BEGS' }  ||= ".";
 668     $record->{ 'BLOCK_LENS' }  ||= ".";
 669     $record->{ 'BLOCK_TYPE' }  ||= ".";
 670     $record->{ 'ALIGN' }       ||= $record->{ 'DESCRIPTOR' } || ".";
 671
 672     return wantarray ? %{ $record } : $record;
 673 }
 674
 675
 676 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 677
 678 1;
 679