code_perl/Maasha/Fasta.pm

   1 package Maasha::Fasta;
   2
   3 # Copyright (C) 2006-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 manipulation of FASTA files and FASTA entries.
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use warnings;
  32 use strict;
  33 use Exporter;
  34 use Data::Dumper;
  35 use Maasha::Common;
  36 use Maasha::Filesys;
  37 use Maasha::Seq;
  38 use vars qw ( @ISA @EXPORT );
  39
  40 @ISA = qw( Exporter );
  41
  42 use constant {
  43     SEQ_NAME => 0,
  44     SEQ  => 1,
  45 };
  46
  47
  48 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  49
  50
  51 sub get_entry
  52 {
  53     # Martin A. Hansen, January 2007.
  54
  55     # Given a filehandle to an FASTA file,
  56     # fetches the next FASTA entry, and returns
  57     # this as a tuple of [ header, sequence ].
  58
  59     my ( $fh,    # filehandle to FASTA file
  60        ) = @_;
  61
  62     # Returns string.
  63
  64     my ( $block, @lines, $seq_name, $seq, $entry );
  65
  66     local $/ = "\n>";
  67
  68     while ( $block = <$fh> )
  69     {
  70         chomp $block;
  71
  72         last if $block !~ /^\s+$/;
  73     }
  74
  75     return if not defined $block;
  76
  77     # $block =~ />?([^\n]+)\n/m;
  78     $block =~ /^>?(.+)\n/m;
  79     $seq_name = $1;
  80     $seq      = $';
  81
  82     local $/ = "\n";
  83
  84     chomp $seq;
  85
  86     $seq_name =~ tr/\r//d;
  87     $seq      =~ tr/ \t\n\r//d;
  88
  89     $entry = [ $seq_name, $seq ];
  90
  91     return wantarray ? @{ $entry } : $entry;
  92 }
  93
  94
  95 sub put_entry
  96 {
  97     # Martin A. Hansen, January 2007.
  98
  99     # Writes a FASTA entry to a file handle.
 100
 101     my ( $entry,     # a FASTA entries
 102          $fh,        # file handle to output file - OPTIONAL
 103          $wrap,      # line width                 - OPTIONAL
 104        ) = @_;
 105
 106     # Returns nothing.
 107
 108     Maasha::Common::error( qq(FASTA entry has no header) )   if not defined $entry->[ SEQ_NAME ];
 109     Maasha::Common::error( qq(FASTA entry has no sequence) ) if not defined $entry->[ SEQ ];
 110
 111     if ( $wrap ) {
 112         Maasha::Fasta::wrap( $entry, $wrap );
 113     }
 114
 115     if ( defined $fh ) {
 116         print $fh ">$entry->[ SEQ_NAME ]\n$entry->[ SEQ ]\n";
 117     } else {
 118         print ">$entry->[ SEQ_NAME ]\n$entry->[ SEQ ]\n";
 119     }
 120 }
 121
 122
 123 sub put_entries
 124 {
 125     # Martin A. Hansen, September 2009
 126
 127     # Write FASTA entries to a file.
 128
 129     my ( $entries,   # list of FASTA entries
 130          $file,      # output file
 131          $wrap,      # line width  -  OPTIONAL
 132        ) = @_;
 133
 134     # Returns nothing
 135
 136     my ( $fh );
 137
 138     $fh = Maasha::Filesys::file_write_open( $file );
 139
 140     map { put_entry( $_, $fh, $wrap ) } @{ $entries };
 141
 142     close $fh;
 143 }
 144
 145
 146 sub wrap
 147 {
 148     # Martin A. Hansen, June 2007
 149
 150     # Wraps the sequence of a given FASTA entry
 151     # to a given length.
 152
 153     my ( $entry,   # FASTA entry
 154          $wrap,    # wrap length
 155        ) = @_;
 156
 157     # Returns nothing.
 158
 159     Maasha::Seq::wrap( \$entry->[ SEQ ], $wrap );
 160 }
 161
 162
 163 sub fasta_index
 164 {
 165     # Martin A. Hansen, July 2008.
 166
 167     # Given a path to a FASTA file create a FASTA index
 168     # and save to the specified path.
 169
 170     my ( $fasta_path,   # path to FASTA file
 171          $index_dir,    # directory
 172          $index_name,   # path to index file
 173        ) = @_;
 174
 175     # Returns nothing.
 176
 177     my ( $index );
 178
 179     $index = index_create( $fasta_path );
 180
 181     Maasha::Filesys::dir_create_if_not_exists( $index_dir );
 182
 183     Maasha::Fasta::index_store( "$index_dir/$index_name", $index );
 184 }
 185
 186
 187 sub index_check
 188 {
 189     # Martin A. Hansen, September 2009.
 190
 191     # Given a path to a FASTA file and a FASTA index
 192     # check if the FILESIZE matches.
 193
 194     my ( $fasta_path,   # path to FASTA file
 195          $index_path,   # path to index file
 196        ) = @_;
 197
 198     # Returns boolean.
 199
 200     my ( $index );
 201
 202     $index = Maasha::Fasta::index_retrieve( $index_path );
 203
 204     if ( Maasha::Filesys::file_size( $fasta_path ) == $index->{ 'FILE_SIZE' } ) {
 205         return 1;
 206     }
 207
 208     return 0;
 209 }
 210
 211
 212 sub index_create
 213 {
 214     # Matin A. Hansen, December 2004.
 215
 216     # Given a FASTA file formatted with one line of header and one line of sequence,
 217     # returns a list of header, seq offset and seq length (first nucleotide is 0). Also,
 218     # the file size of the indexed file is written to the index for checking purposes.
 219
 220     my ( $path,   #  full path to file with multiple FASTA entries
 221        ) = @_;
 222
 223     # Returns a hashref.
 224
 225     my ( $fh, $entry, $offset, $len, $tot, %index );
 226
 227     $index{ 'FILE_SIZE' } = Maasha::Filesys::file_size( $path );
 228
 229     $offset = 0;
 230     $len    = 0;
 231     $tot    = 0;
 232
 233     $fh = Maasha::Filesys::file_read_open( $path );
 234
 235     while ( $entry = get_entry( $fh ) )
 236     {
 237         warn qq(WARNING: header->$entry->[ SEQ_NAME ] alread exists in index) if exists $index{ $entry->[ SEQ_NAME ] };
 238
 239         $offset += $len + 2 + length $entry->[ SEQ_NAME ];
 240         $len  = length $entry->[ SEQ ];
 241
 242         $index{ $entry->[ SEQ_NAME ] } = [ $offset, $len ];
 243
 244         $offset++;
 245
 246         $tot += length( $entry->[ SEQ_NAME ] ) + length( $entry->[ SEQ ] ) + 3;
 247     }
 248
 249     close $fh;
 250
 251     if ( $index{ 'FILE_SIZE' } != $tot ) {
 252         Maasha::Common::error( qq(Filesize "$index{ 'FILE_SIZE' }" does not match content "$tot" -> malformatted FASTA file) );
 253     }
 254
 255     return wantarray ? %index : \%index;
 256 }
 257
 258
 259 sub index_lookup
 260 {
 261     # Martin A. Hansen, July 2007.
 262
 263     # Lookup a list of exact matches in the index and returns these.
 264
 265     my ( $index,     # index list
 266          $headers,   # headers to lookup
 267        ) = @_;
 268
 269     # Returns a list.
 270
 271     my ( $head, @results );
 272
 273     foreach $head ( @{ $headers } )
 274     {
 275         if ( exists $index->{ $head } ) {
 276             push @results, [ $head, $index->{ $head }->[ 0 ], $index->{ $head }->[ 1 ] ];
 277         }
 278     }
 279
 280     return wantarray ? @results : \@results;
 281 }
 282
 283
 284 sub index_store
 285 {
 286     # Martin A. Hansen, May 2007.
 287
 288     # Stores a FASTA index to binary file.
 289
 290     my ( $path,   # full path to file
 291          $index,  # list with index
 292        ) = @_;
 293
 294     # returns nothing
 295
 296     Maasha::Filesys::file_store( $path, $index );
 297 }
 298
 299
 300 sub index_retrieve
 301 {
 302     # Martin A. Hansen, May 2007.
 303
 304     # Retrieves a FASTA index from binary file.
 305
 306     my ( $path,   # full path to file
 307        ) = @_;
 308
 309     # returns list
 310
 311     my $index;
 312
 313     $index = Maasha::Filesys::file_retrieve( $path );
 314
 315     return wantarray ? @{ $index } : $index;
 316 }
 317
 318
 319 sub biopiece2fasta
 320 {
 321     # Martin A. Hansen, July 2008.
 322
 323     # Given a biopiece record converts it to a FASTA record.
 324     # If no generic SEQ or SEQ_NAME is found, the Q_* and S_* are
 325     # tried in that order.
 326
 327     my ( $record,    # record
 328        ) = @_;
 329
 330     # Returns a tuple.
 331
 332     my ( $seq_name, $seq );
 333
 334     $seq_name = $record->{ "SEQ_NAME" } || $record->{ "Q_ID" }  || $record->{ "S_ID" };
 335     $seq      = $record->{ "SEQ" }      || $record->{ "Q_SEQ" } || $record->{ "S_SEQ" };
 336
 337     if ( defined $seq_name and defined $seq ) {
 338         return wantarray ? ( $seq_name, $seq ) : [ $seq_name, $seq ];
 339     } else {
 340         return;
 341     }
 342 }
 343
 344
 345 sub fasta2biopiece
 346 {
 347     # Martin A. Hansen, May 2009
 348
 349     # Convert a FASTA record to a Biopiece record.
 350
 351     my ( $entry,   # FASTA entry
 352        ) = @_;
 353
 354     # Returns a hashref.
 355
 356     my ( $record );
 357
 358     if ( defined $entry->[ SEQ_NAME ] and $entry->[ SEQ ] )
 359     {
 360         $record = {
 361             SEQ_NAME => $entry->[ SEQ_NAME ],
 362             SEQ      => $entry->[ SEQ ],
 363             SEQ_LEN  => length $entry->[ SEQ ],
 364         };
 365     }
 366
 367     return wantarray ? %{ $record } : $record;
 368 }
 369
 370
 371
 372 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<