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