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