code_perl/Maasha/KISS/IO.pm

   1 package Maasha::KISS::IO;
   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::SQL;
  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 };
  57
  58
  59 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  60
  61
  62 sub kiss_entry_get
  63 {
  64     my ( $fh,    # file handle
  65        ) = @_;
  66
  67     # Returns a hashref.
  68
  69     my ( $line, @fields, %entry );
  70
  71     while ( $line = <$fh> )
  72     {
  73         chomp $line;
  74
  75         next if $line =~ /^$|^#/;
  76
  77         @fields = split /\t/, $line;
  78
  79         Maasha::Common::error( qq( BAD kiss entry: $line) ) if not @fields == 12;
  80
  81         $entry{ 'S_ID' }        = $fields[ S_ID ];
  82         $entry{ 'S_BEG' }       = $fields[ S_BEG ];
  83         $entry{ 'S_END' }       = $fields[ S_END ];
  84         $entry{ 'Q_ID' }        = $fields[ Q_ID ];
  85         $entry{ 'SCORE' }       = $fields[ SCORE ];
  86         $entry{ 'STRAND' }      = $fields[ STRAND ];
  87         $entry{ 'HITS' }        = $fields[ HITS ];
  88         $entry{ 'ALIGN' }       = $fields[ ALIGN ];
  89         $entry{ 'BLOCK_COUNT' } = $fields[ BLOCK_COUNT ];
  90         $entry{ 'BLOCK_BEGS' }  = $fields[ BLOCK_BEGS ];
  91         $entry{ 'BLOCK_LENS' }  = $fields[ BLOCK_LENS ];
  92         $entry{ 'BLOCK_TYPE' }  = $fields[ BLOCK_TYPE ];
  93
  94         return wantarray ? %entry : \%entry;
  95     }
  96 }
  97
  98
  99 sub kiss_entry_put
 100 {
 101     my ( $entry,   # KISS entry to output
 102          $fh,      # file handle  -  OPTIONAL
 103        ) = @_;
 104
 105     # Returns nothing.
 106
 107     my ( @fields );
 108
 109     $fh ||= \*STDOUT;
 110
 111     $fields[ S_ID ]        = $entry->{ 'S_ID' };
 112     $fields[ S_BEG ]       = $entry->{ 'S_BEG' };
 113     $fields[ S_END ]       = $entry->{ 'S_END' };
 114     $fields[ Q_ID ]        = $entry->{ 'Q_ID' };
 115     $fields[ SCORE ]       = $entry->{ 'SCORE' };
 116     $fields[ STRAND ]      = $entry->{ 'STRAND' };
 117     $fields[ HITS ]        = $entry->{ 'HITS' };
 118     $fields[ ALIGN ]       = $entry->{ 'ALIGN' };
 119     $fields[ BLOCK_COUNT ] = $entry->{ 'BLOCK_COUNT' };
 120     $fields[ BLOCK_BEGS ]  = $entry->{ 'BLOCK_BEGS' };
 121     $fields[ BLOCK_LENS ]  = $entry->{ 'BLOCK_LENS' };
 122     $fields[ BLOCK_TYPE ]  = $entry->{ 'BLOCK_TYPE' };
 123
 124     print $fh join( "\t", @fields ), "\n";
 125 }
 126
 127
 128 sub kiss_sql_get
 129 {
 130     my ( $dbh,     # Database handle
 131          $table,   # Table name
 132          $s_beg,   # Subject begin
 133          $s_end,   # Subject end
 134        ) = @_;
 135
 136     my ( $sql, $entries );
 137
 138     # $sql = "SELECT * FROM $table WHERE S_BEG >= $s_beg AND S_END <= $s_end ORDER BY S_BEG,S_END";
 139     $sql = "SELECT S_BEG,S_END,Q_ID,ALIGN FROM $table WHERE S_BEG >= $s_beg AND S_END <= $s_end";
 140
 141     $entries = Maasha::SQL::query_hashref_list( $dbh, $sql );
 142
 143     return wantarray ? @{ $entries } : $entries;
 144 }
 145
 146
 147 sub kiss_index
 148 {
 149     # Martin A, Hansen, October 2009.
 150
 151     # Creates an index of a sorted KISS file that
 152     # allowing the location of the byte position
 153     # from where records can be read given a
 154     # specific S_BEG position. The index consists of
 155     # triples: [ beg, end, bytepos ], where beg and
 156     # end denotes the interval where the next KISS
 157     # record begins at bytepos.
 158
 159     my ( $fh,   # filehandle to KISS file
 160        ) = @_;
 161
 162     # Returns a list.
 163
 164     my ( $line, @fields, $beg, $end, $pos, @index );
 165
 166     $beg = 0;
 167     $pos = 0;
 168
 169     while ( $line = <$fh> )
 170     {
 171         chomp $line;
 172
 173         @fields = split /\t/, $line, 3;
 174
 175         $end = $fields[ S_BEG ];
 176
 177         if ( $end == 0 )
 178         {
 179             push @index, [ $beg, $end, $pos ];
 180             $beg = 1;
 181         }
 182         elsif ( $end > $beg )
 183         {
 184             push @index, [ $beg, $end - 1, $pos ];
 185             $beg = $end;
 186         }
 187         elsif( $end < $beg )
 188         {
 189             Maasha::Common::error( qq(KISS file not sorted: $end < $beg) );
 190         }
 191
 192         $pos += 1 + length $line;
 193     }
 194
 195     return wantarray ? @index : \@index;
 196 }
 197
 198
 199 sub kiss_index_store
 200 {
 201     my ( $path,
 202          $index,
 203        ) = @_;
 204
 205     Maasha::Filesys::file_store( $path, $index );
 206 }
 207
 208
 209 sub kiss_index_retrieve
 210 {
 211     my ( $path,
 212        ) = @_;
 213
 214     my $index;
 215
 216     $index = Maasha::Filesys::file_retrieve( $path );
 217
 218     return wantarray ? @{ $index } : $index;
 219 }
 220
 221
 222 sub kiss_index_search
 223 {
 224     my ( $index,
 225          $num,
 226        ) = @_;
 227
 228     # Returns a number.
 229
 230     my ( $high, $low, $try );
 231
 232     $low  = 0;
 233     $high = scalar @{ $index };
 234
 235     while ( $low <= $high )
 236     {
 237         $try = int( ( $high + $low ) / 2 );
 238
 239         if ( $num < $index->[ $try ]->[ 0 ] ) {
 240             $high = $try;
 241         } elsif ( $num > $index->[ $try ]->[ 1 ] ) {
 242             $low = $try + 1;
 243         } else {
 244             return $index->[ $try ]->[ 2 ];
 245         }
 246     }
 247
 248     Maasha::Common::error( "Could not find number->$num in index" );
 249 }
 250
 251
 252 sub kiss_index_get
 253 {
 254     my ( $file,
 255          $beg,
 256          $end,
 257        ) = @_;
 258
 259     my ( $index, $offset, $fh, $entry, @entries );
 260
 261     $index = Maasha::KISS::IO::kiss_index_retrieve( "$file.index" );
 262
 263     $offset = Maasha::KISS::IO::kiss_index_search( $index, $beg );
 264
 265     $fh = Maasha::Filesys::file_read_open( $file );
 266
 267     sysseek( $fh, $offset, 0 );
 268
 269     while ( $entry = kiss_entry_get( $fh ) )
 270     {
 271         push @entries, $entry;
 272
 273         last if $entry->{ 'S_END' } > $end;
 274     }
 275
 276     close $fh;
 277
 278     return wantarray ? @entries : \@entries;
 279 }
 280
 281
 282 sub kiss_align
 283 {
 284     # S.aur_COL       41      75      5_vnlh2ywXsN1/1 1       -       .       17:A>T,31:A>N   .       .       .       .
 285
 286     my ( $s_seqref,   # subject sequence reference
 287          $entry,      # KISS entry
 288        ) = @_;
 289
 290     # Returns a string.
 291
 292     my ( $align, $pos, $s_symbol, $q_symbol, $s_seq, $q_seq, $insertions );
 293
 294     $s_seq = substr ${ $s_seqref }, $entry->{ 'S_BEG' }, $entry->{ 'S_END' } - $entry->{ 'S_BEG' } + 1;
 295     $q_seq = $s_seq;
 296
 297     $insertions = 0;
 298
 299     foreach $align ( split /,/, $entry->{ 'ALIGN' } )
 300     {
 301         if ( $align =~ /(\d+):(.)>(.)/ )
 302         {
 303             $pos      = $1;
 304             $s_symbol = $2;
 305             $q_symbol = $3;
 306
 307             if ( $s_symbol eq '-' ) # insertion
 308             {
 309                 substr $s_seq, $pos + $insertions, 0, $s_symbol;
 310                 substr $q_seq, $pos + $insertions, 0, $q_symbol;
 311
 312                 $insertions++;
 313             }
 314             elsif ( $q_symbol eq '-' ) # deletion
 315             {
 316                 substr $q_seq, $pos + $insertions, 1, $q_symbol;
 317             }
 318             else # mismatch
 319             {
 320                 substr $q_seq, $pos + $insertions, 1, $q_symbol;
 321             }
 322         }
 323         else
 324         {
 325             Maasha::Common::error( qq(Bad alignment descriptor: "$align") );
 326         }
 327     }
 328
 329     print Dumper( "s_seq: $s_seq", "q_seq: $q_seq" ) and exit;
 330 }
 331
 332
 333 sub sam2kiss
 334 {
 335     my ( $sam_entry,   # SAM entry
 336        ) = @_;
 337
 338     # Returns a hashref
 339
 340     my ( $cigar, $offset, $op, $len, @descriptors, $kiss_entry );
 341
 342     $cigar = $sam_entry->{ 'CIGAR' };
 343
 344     $cigar =~ tr/\*//d;
 345
 346     $offset = 0;
 347
 348     while ( length $cigar > 0 )
 349     {
 350         if ( $cigar =~ s/^([MINDSHP])(\d+)// )
 351         {
 352             $op  = $1;
 353             $len = $2;
 354
 355             print "CIGAR: $cigar   OP: $op   LEN: $len\n";
 356
 357             if ( $op eq 'I' )
 358             {
 359
 360             }
 361             elsif ( $op eq 'D' )
 362             {
 363
 364             }
 365
 366             $offset += $len;
 367         }
 368         else
 369         {
 370             Maasha::Common::error( qq(Bad CIGAR format: "$cigar") );
 371         }
 372     }
 373
 374     return wantarray ? %{ $kiss_entry } : $kiss_entry;
 375 }
 376
 377
 378 sub kiss2biopiece
 379 {
 380     my ( $entry,   # KISS entry
 381        ) = @_;
 382
 383     return wantarray ? %{ $entry } : $entry;
 384 }
 385
 386
 387 sub biopiece2kiss
 388 {
 389     my ( $record,   # Biopiece record
 390        ) = @_;
 391
 392     $record->{ 'HITS' }        ||= ".";
 393     $record->{ 'BLOCK_COUNT' } ||= ".";
 394     $record->{ 'BLOCK_BEGS' }  ||= ".";
 395     $record->{ 'BLOCK_LENS' }  ||= ".";
 396     $record->{ 'BLOCK_TYPE' }  ||= ".";
 397     $record->{ 'ALIGN' }       ||= $record->{ 'DESCRIPTOR' } || ".";
 398
 399     return wantarray ? %{ $record } : $record;
 400 }
 401
 402
 403 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 404
 405 1;
 406
 407