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