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