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