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         # print "low: $low   high: $high   try: $try\n";
 260
 261         if ( $num < $index->[ $try ]->[ 0 ] ) {
 262             $high = $try;
 263         } elsif ( $num > $index->[ $try ]->[ 1 ] ) {
 264             $low = $try + 1;
 265         } else {
 266             return $index->[ $try ]->[ 2 ];
 267         }
 268     }
 269
 270     Maasha::Common::error( "Could not find number->$num in index" );
 271 }
 272
 273
 274 sub kiss_index_get
 275 {
 276     my ( $file,
 277          $beg,
 278          $end,
 279        ) = @_;
 280
 281     my ( $index, $offset, $fh, $entry, @entries );
 282
 283     $index = Maasha::KISS::IO::kiss_index_retrieve( "$file.index" );
 284
 285     $offset = Maasha::KISS::IO::kiss_index_search( $index, $beg );
 286
 287     $fh = Maasha::Filesys::file_read_open( $file );
 288
 289     sysseek( $fh, $offset, 0 );
 290
 291     while ( $entry = kiss_entry_get( $fh ) )
 292     {
 293         push @entries, $entry;
 294
 295         last if $entry->{ 'S_END' } > $end;
 296     }
 297
 298     close $fh;
 299
 300     return wantarray ? @entries : \@entries;
 301 }
 302
 303
 304 sub kiss2biopiece
 305 {
 306     my ( $entry,   # KISS entry
 307        ) = @_;
 308
 309     return wantarray ? %{ $entry } : $entry;
 310 }
 311
 312
 313 sub biopiece2kiss
 314 {
 315     my ( $record,   # Biopiece record
 316        ) = @_;
 317
 318     $record->{ 'HITS' }        ||= ".";
 319     $record->{ 'BLOCK_COUNT' } ||= ".";
 320     $record->{ 'BLOCK_BEGS' }  ||= ".";
 321     $record->{ 'BLOCK_LENS' }  ||= ".";
 322     $record->{ 'BLOCK_TYPE' }  ||= ".";
 323     $record->{ 'ALIGN' }       ||= $record->{ 'DESCRIPTOR' } || ".";
 324
 325     return wantarray ? %{ $record } : $record;
 326 }
 327
 328
 329 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 330
 331 1;
 332
 333