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