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