code_perl/Maasha/Fasta.pm

   1 package Maasha::Fasta;
   2
   3 # Copyright (C) 2006 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 strict;
  32 use Data::Dumper;
  33 use Maasha::Common;
  34 use Maasha::Seq;
  35 use vars qw ( @ISA @EXPORT );
  36
  37 @ISA = qw( Exporter );
  38
  39 use constant {
  40     HEAD => 0,
  41     SEQ  => 1,
  42 };
  43
  44
  45 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  46
  47
  48 sub fasta_format_ok
  49 {
  50     # Martin A. Hansen, March 2007.
  51
  52     # Checks if a given FASTA file is formatted with
  53     # one header per line and one sequence per line.
  54     # returns 1 if so, else 0.
  55
  56     my ( $path,   # full path to FASTA file
  57        ) = @_;
  58
  59     # Returns boolean
  60
  61     my ( $fh, $line, $count );
  62
  63     $fh = &Maasha::Common::read_open( $path );
  64
  65     $count = 0;
  66
  67     while ( $line = <$fh> )
  68     {
  69         if ( not $count % 2 and substr( $line, 0, 1 ) ne ">" ) {
  70             return 0;
  71         }
  72
  73         $count++;
  74     }
  75
  76     close $fh;
  77
  78     return 1;
  79 }
  80
  81
  82 sub get_entries
  83 {
  84     # Martin A. Hansen, December 2006.
  85
  86     # Parses a fasta file and returns a list of headers and sequence tuples.
  87
  88     my ( $path,   # full path to FASTA file
  89          $count,  # number of sequences to read - OPTIONAL
  90        ) = @_;
  91
  92     # returns list of tuples
  93
  94     my ( $fh, $entry, @entries );
  95
  96     $fh = &Maasha::Common::read_open( $path );
  97
  98     while ( $entry = &get_entry( $fh ) )
  99     {
 100         push @entries, $entry;
 101
 102         if ( $count and $count == @entries ) {
 103             last;
 104         }
 105     }
 106
 107     close $fh;
 108
 109     return wantarray ? @entries : \@entries;
 110 }
 111
 112
 113 sub put_entries
 114 {
 115     # Martin A. Hansen, March 2004.
 116
 117     # writes fasta sequences to STDOUT or file
 118
 119     my ( $entries,   # list of fasta entries
 120          $path,      # full path to file - OPTIONAL
 121          $wrap,      # line width        - OPTIONAL
 122        ) = @_;
 123
 124     my ( $fh );
 125
 126     $fh = &Maasha::Common::write_open( $path ) if $path;
 127
 128     map { &put_entry( $_, $fh, $wrap ) } @{ $entries };
 129
 130     close $fh if defined;
 131 }
 132
 133
 134 sub wrap
 135 {
 136     # Martin A. Hansen, June 2007
 137
 138     # Wraps the sequence of a given FASTA entry
 139     # to a given length.
 140
 141     my ( $entry,   # FASTA entry
 142          $wrap,    # wrap length
 143        ) = @_;
 144
 145     # Returns nothing.
 146
 147     &Maasha::Seq::wrap( \$entry->[ SEQ ], $wrap );
 148 }
 149
 150
 151 sub get_entry
 152 {
 153     # Martin A. Hansen, January 2007.
 154
 155     # Given a filehandle to an FASTA file,
 156     # fetches the next FASTA entry, and returns
 157     # this as a tuple of [ header, sequence ].
 158
 159     my ( $fh,    # filehandle to FASTA file
 160        ) = @_;
 161
 162     # Returns string.
 163
 164     my ( $block, @lines, $seq_name, $seq, $entry );
 165
 166     local $/ = "\n>";
 167
 168     while ( $block = <$fh> )
 169     {
 170         chomp $block;
 171
 172         last if $block !~ /^\s+$/;
 173     }
 174
 175     return if not defined $block;
 176
 177     $block =~ />?([^\n]+)\n/m;
 178     $seq_name = $1;
 179     $seq      = $';
 180
 181     local $/ = "\n";
 182
 183     chomp $seq;
 184
 185     $seq =~ tr/ \t\n//d;
 186
 187     $entry = [ $seq_name, $seq ];
 188
 189     return wantarray ? @{ $entry } : $entry;
 190 }
 191
 192
 193 sub put_entry
 194 {
 195     # Martin A. Hansen, January 2007.
 196
 197     # Writes FASTA entries to STDOUT or file.
 198
 199     my ( $entry,     # a FASTA entries
 200          $fh,        # file handle to output file - OPTIONAL
 201          $wrap,      # line width                 - OPTIONAL
 202        ) = @_;
 203
 204     # Returns nothing.
 205
 206     &Maasha::Common::error( qq(FASTA entry has no header) )   if not defined $entry->[ HEAD ];
 207     &Maasha::Common::error( qq(FASTA entry has no sequence) ) if not defined $entry->[ SEQ ];
 208
 209     if ( $wrap ) {
 210         &Maasha::Fasta::wrap( $entry, $wrap );
 211     }
 212
 213     if ( defined $fh ) {
 214         print $fh ">$entry->[ HEAD ]\n$entry->[ SEQ ]\n";
 215     } else {
 216         print ">$entry->[ HEAD ]\n$entry->[ SEQ ]\n";
 217     }
 218 }
 219
 220
 221 sub find_shortest
 222 {
 223     # Martin A. Hansen, June 2007.
 224
 225     # Given a stack of FASTA entries, find and return
 226     # the shortest entry.
 227
 228     my ( $entries,   # list of FASTA entries
 229        ) = @_;
 230
 231     # returns tuple
 232
 233     my ( $min, $entry, $min_entry );
 234
 235     $min = 99999999999;
 236
 237     foreach $entry ( @{ $entries } )
 238     {
 239         if ( length( $entry->[ SEQ ] ) < $min )
 240         {
 241             $min_entry = $entry;
 242             $min = length $entry->[ SEQ ];
 243         }
 244     }
 245
 246     return wantarray ? @{ $min_entry } : $min_entry;
 247 }
 248
 249
 250 sub find_longest
 251 {
 252     # Martin A. Hansen, June 2007.
 253
 254     # Given a stack of FASTA entries, find and return
 255     # the longest entry.
 256
 257     my ( $entries,   # list of FASTA entries
 258        ) = @_;
 259
 260     # returns tuple
 261
 262     my ( $max, $entry, $max_entry );
 263
 264     $max = 0;
 265
 266     foreach $entry ( @{ $entries } )
 267     {
 268         if ( length( $entry->[ SEQ ] ) > $max )
 269         {
 270             $max_entry = $entry;
 271             $max = length $entry->[ SEQ ];
 272         }
 273     }
 274
 275     return wantarray ? @{ $max_entry } : $max_entry;
 276 }
 277
 278
 279 sub fasta_get_headers
 280 {
 281     # Martin A. Hansen, May 2007.
 282
 283     # Gets the header names of a FASTA file,
 284     # and returns these in a list.
 285
 286     my ( $path,   # full path to FASTA file
 287        ) = @_;
 288
 289     # returns list
 290
 291     my ( $fh, $entry, @list );
 292
 293     $fh = &Maasha::Common::read_open( $path );
 294
 295     while ( $entry = &get_entry( $fh ) ) {
 296         push @list, $entry->[ HEAD ];
 297     }
 298
 299     close $fh;
 300
 301     return wantarray ? @list : \@list;
 302 }
 303
 304
 305 sub fasta_reformat
 306 {
 307     # Martin A. Hansen, December 2004.
 308
 309     # Given a file of one or more FASTA entries, reformats these so
 310     # each entry consits of one line with header and one line with sequence.
 311
 312     my ( $path,   # full path to file with multiple FASTA entries
 313        ) = @_;
 314
 315     my ( $fh_in, $fh_out, $entry );
 316
 317     $fh_in  = &Maasha::Common::read_open( $path );
 318     $fh_out = &Maasha::Common::write_open( "$path.temp" );
 319
 320     while ( $entry = &get_entry( $fh_in ) ) {
 321         &put_entry( $entry, $fh_out );
 322     }
 323
 324     close $fh_in;
 325     close $fh_out;
 326
 327     rename( "$path.temp", $path );
 328 }
 329
 330
 331 sub index_create
 332 {
 333     # Matin A. Hansen, December 2004.
 334
 335     # Given a FASTA file formatted with one line of header and one line of sequence,
 336     # returns a list of header, seq beg and seq length (first nucleotide is 0). Also,
 337     # the file size of the indexed file is written to the index for checking purposes.
 338
 339     my ( $path,   #  full path to file with multiple FASTA entries
 340        ) = @_;
 341
 342     # returns a hashref
 343
 344     my ( $file_size, $fh, $entry, $beg, $len, %hash, @index );
 345
 346     $file_size = &Maasha::Common::file_size( $path );
 347
 348     push @index, "FILE_SIZE=$file_size";
 349
 350     $fh = &Maasha::Common::read_open( $path );
 351
 352     $beg = 0;
 353     $len = 0;
 354
 355     while ( $entry = &get_entry( $fh ) )
 356     {
 357         warn qq(WARNING: header->$entry->[ HEAD ] alread exists in index) if exists $hash{ $entry->[ HEAD ] };
 358
 359         $beg += $len + 2 + length $entry->[ HEAD ];
 360         $len = length $entry->[ SEQ ];
 361
 362         push @index, [ $entry->[ HEAD ], $beg, $len ];
 363
 364         $hash{ $entry->[ HEAD ] } = 1;
 365
 366         $beg++;
 367     }
 368
 369     close $fh;
 370
 371     return wantarray ? @index : \@index;
 372 }
 373
 374
 375 sub index_search
 376 {
 377     # Martin A. Hansen, December 2004.
 378
 379     # Searches the index for matching entries.
 380
 381     my ( $index,   # index list
 382          $regex,   # regex to match FASTA headers [OPTIONAL]
 383          $invert,  # invert matching
 384        ) = @_;
 385
 386     # returns list
 387
 388     my ( @results );
 389
 390     if ( not $regex )
 391     {
 392         @results = @{ $index };
 393     }
 394     else
 395     {
 396         if ( $invert ) {
 397             @results = grep { $_->[ 0 ] !~ /$regex/ } @{ $index };
 398         } else {
 399             @results = grep { $_->[ 0 ] =~ /$regex/ } @{ $index };
 400         }
 401     }
 402
 403     return wantarray ? @results : \@results;
 404 }
 405
 406
 407 sub index_lookup
 408 {
 409     # Martin A. Hansen, July 2007.
 410
 411     # Lookup a list of exact matches in the index and returns these
 412
 413     my ( $index,     # index list
 414          $headers,   # headers to lookup
 415        ) = @_;
 416
 417     # returns a list
 418
 419     my ( %hash, $head, @results );
 420
 421     map { $hash{ $_->[ 0 ] } = [ $_->[ 1 ], $_->[ 2 ] ] } @{ $index };
 422
 423     foreach $head ( @{ $headers } )
 424     {
 425         if ( exists $hash{ $head } ) {
 426             push @results, [ $head, $hash{ $head }->[ 0 ], $hash{ $head }->[ 1 ] ];
 427         }
 428     }
 429
 430     return wantarray ? @results : \@results;
 431 }
 432
 433
 434 sub index_store
 435 {
 436     # Martin A. Hansen, May 2007.
 437
 438     # Stores a FASTA index to binary file.
 439
 440     my ( $path,   # full path to file
 441          $index,  # list with index
 442        ) = @_;
 443
 444     # returns nothing
 445
 446     &Maasha::Common::file_store( $path, $index );
 447 }
 448
 449
 450 sub index_retrieve
 451 {
 452     # Martin A. Hansen, May 2007.
 453
 454     # Retrieves a FASTA index from binary file.
 455
 456     my ( $path,   # full path to file
 457        ) = @_;
 458
 459     # returns list
 460
 461     my $index;
 462
 463     $index = &Maasha::Common::file_retrieve( $path );
 464
 465     return wantarray ? @{ $index } : $index;
 466 }
 467
 468
 469 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<