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