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