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     $block =~ /^>?(.+)\n/m;
 179     $seq_name = $1;
 180     $seq      = $';
 181
 182     local $/ = "\n";
 183
 184     chomp $seq;
 185
 186     $seq_name =~ tr/\r//d;
 187     $seq      =~ tr/ \t\n\r//d;
 188
 189     $entry = [ $seq_name, $seq ];
 190
 191     return wantarray ? @{ $entry } : $entry;
 192 }
 193
 194
 195 sub put_entry
 196 {
 197     # Martin A. Hansen, January 2007.
 198
 199     # Writes FASTA entries to STDOUT or file.
 200
 201     my ( $entry,     # a FASTA entries
 202          $fh,        # file handle to output file - OPTIONAL
 203          $wrap,      # line width                 - OPTIONAL
 204        ) = @_;
 205
 206     # Returns nothing.
 207
 208     Maasha::Common::error( qq(FASTA entry has no header) )   if not defined $entry->[ HEAD ];
 209     Maasha::Common::error( qq(FASTA entry has no sequence) ) if not defined $entry->[ SEQ ];
 210
 211     if ( $wrap ) {
 212         Maasha::Fasta::wrap( $entry, $wrap );
 213     }
 214
 215     if ( defined $fh ) {
 216         print $fh ">$entry->[ HEAD ]\n$entry->[ SEQ ]\n";
 217     } else {
 218         print ">$entry->[ HEAD ]\n$entry->[ SEQ ]\n";
 219     }
 220 }
 221
 222
 223 sub find_shortest
 224 {
 225     # Martin A. Hansen, June 2007.
 226
 227     # Given a stack of FASTA entries, find and return
 228     # the shortest entry.
 229
 230     my ( $entries,   # list of FASTA entries
 231        ) = @_;
 232
 233     # returns tuple
 234
 235     my ( $min, $entry, $min_entry );
 236
 237     $min = 99999999999;
 238
 239     foreach $entry ( @{ $entries } )
 240     {
 241         if ( length( $entry->[ SEQ ] ) < $min )
 242         {
 243             $min_entry = $entry;
 244             $min = length $entry->[ SEQ ];
 245         }
 246     }
 247
 248     return wantarray ? @{ $min_entry } : $min_entry;
 249 }
 250
 251
 252 sub find_longest
 253 {
 254     # Martin A. Hansen, June 2007.
 255
 256     # Given a stack of FASTA entries, find and return
 257     # the longest entry.
 258
 259     my ( $entries,   # list of FASTA entries
 260        ) = @_;
 261
 262     # returns tuple
 263
 264     my ( $max, $entry, $max_entry );
 265
 266     $max = 0;
 267
 268     foreach $entry ( @{ $entries } )
 269     {
 270         if ( length( $entry->[ SEQ ] ) > $max )
 271         {
 272             $max_entry = $entry;
 273             $max = length $entry->[ SEQ ];
 274         }
 275     }
 276
 277     return wantarray ? @{ $max_entry } : $max_entry;
 278 }
 279
 280
 281 sub fasta_get_headers
 282 {
 283     # Martin A. Hansen, May 2007.
 284
 285     # Gets the header names of a FASTA file,
 286     # and returns these in a list.
 287
 288     my ( $path,   # full path to FASTA file
 289        ) = @_;
 290
 291     # returns list
 292
 293     my ( $fh, $entry, @list );
 294
 295     $fh = Maasha::Common::read_open( $path );
 296
 297     while ( $entry = get_entry( $fh ) ) {
 298         push @list, $entry->[ HEAD ];
 299     }
 300
 301     close $fh;
 302
 303     return wantarray ? @list : \@list;
 304 }
 305
 306
 307 sub fasta_reformat
 308 {
 309     # Martin A. Hansen, December 2004.
 310
 311     # Given a file of one or more FASTA entries, reformats these so
 312     # each entry consits of one line with header and one line with sequence.
 313
 314     my ( $path,   # full path to file with multiple FASTA entries
 315        ) = @_;
 316
 317     my ( $fh_in, $fh_out, $entry );
 318
 319     $fh_in  = Maasha::Common::read_open( $path );
 320     $fh_out = Maasha::Common::write_open( "$path.temp" );
 321
 322     while ( $entry = get_entry( $fh_in ) ) {
 323         put_entry( $entry, $fh_out );
 324     }
 325
 326     close $fh_in;
 327     close $fh_out;
 328
 329     rename( "$path.temp", $path );
 330 }
 331
 332
 333 sub fasta_index
 334 {
 335     # Martin A. Hansen, July 2008.
 336
 337     # Given a path to a FASTA file create a FASTA index
 338     # and save to the specified path.
 339
 340     my ( $fasta_path,   # path to FASTA file
 341          $index_path,   # path to index file
 342        ) = @_;
 343
 344     # Returns nothing.
 345
 346     my ( $index );
 347
 348     $index = index_create( $fasta_path );
 349
 350     Maasha::Fasta::index_store( $index_path, $index );
 351 }
 352
 353 sub index_create
 354 {
 355     # Matin A. Hansen, December 2004.
 356
 357     # Given a FASTA file formatted with one line of header and one line of sequence,
 358     # returns a list of header, seq beg and seq length (first nucleotide is 0). Also,
 359     # the file size of the indexed file is written to the index for checking purposes.
 360
 361     my ( $path,   #  full path to file with multiple FASTA entries
 362        ) = @_;
 363
 364     # returns a hashref
 365
 366     my ( $file_size, $fh, $entry, $beg, $len, %hash, @index );
 367
 368     $file_size = Maasha::Common::file_size( $path );
 369
 370     push @index, "FILE_SIZE=$file_size";
 371
 372     $fh = Maasha::Common::read_open( $path );
 373
 374     $beg = 0;
 375     $len = 0;
 376
 377     while ( $entry = get_entry( $fh ) )
 378     {
 379         warn qq(WARNING: header->$entry->[ HEAD ] alread exists in index) if exists $hash{ $entry->[ HEAD ] };
 380
 381         $beg += $len + 2 + length $entry->[ HEAD ];
 382         $len = length $entry->[ SEQ ];
 383
 384         push @index, [ $entry->[ HEAD ], $beg, $len ];
 385
 386         $hash{ $entry->[ HEAD ] } = 1;
 387
 388         $beg++;
 389     }
 390
 391     close $fh;
 392
 393     return wantarray ? @index : \@index;
 394 }
 395
 396
 397 sub index_search
 398 {
 399     # Martin A. Hansen, December 2004.
 400
 401     # Searches the index for matching entries.
 402
 403     my ( $index,   # index list
 404          $regex,   # regex to match FASTA headers [OPTIONAL]
 405          $invert,  # invert matching
 406        ) = @_;
 407
 408     # returns list
 409
 410     my ( @results );
 411
 412     if ( not $regex )
 413     {
 414         @results = @{ $index };
 415     }
 416     else
 417     {
 418         if ( $invert ) {
 419             @results = grep { $_->[ 0 ] !~ /$regex/ } @{ $index };
 420         } else {
 421             @results = grep { $_->[ 0 ] =~ /$regex/ } @{ $index };
 422         }
 423     }
 424
 425     return wantarray ? @results : \@results;
 426 }
 427
 428
 429 sub index_lookup
 430 {
 431     # Martin A. Hansen, July 2007.
 432
 433     # Lookup a list of exact matches in the index and returns these
 434
 435     my ( $index,     # index list
 436          $headers,   # headers to lookup
 437        ) = @_;
 438
 439     # returns a list
 440
 441     my ( %hash, $head, @results );
 442
 443     map { $hash{ $_->[ 0 ] } = [ $_->[ 1 ], $_->[ 2 ] ] } @{ $index };
 444
 445     foreach $head ( @{ $headers } )
 446     {
 447         if ( exists $hash{ $head } ) {
 448             push @results, [ $head, $hash{ $head }->[ 0 ], $hash{ $head }->[ 1 ] ];
 449         }
 450     }
 451
 452     return wantarray ? @results : \@results;
 453 }
 454
 455
 456 sub index_store
 457 {
 458     # Martin A. Hansen, May 2007.
 459
 460     # Stores a FASTA index to binary file.
 461
 462     my ( $path,   # full path to file
 463          $index,  # list with index
 464        ) = @_;
 465
 466     # returns nothing
 467
 468     Maasha::Common::file_store( $path, $index );
 469 }
 470
 471
 472 sub index_retrieve
 473 {
 474     # Martin A. Hansen, May 2007.
 475
 476     # Retrieves a FASTA index from binary file.
 477
 478     my ( $path,   # full path to file
 479        ) = @_;
 480
 481     # returns list
 482
 483     my $index;
 484
 485     $index = Maasha::Common::file_retrieve( $path );
 486
 487     return wantarray ? @{ $index } : $index;
 488 }
 489
 490
 491 sub biopiece2fasta
 492 {
 493     # Martin A. Hansen, July 2008.
 494
 495     # Given a biopiece record converts it to a FASTA record.
 496     # If no generic SEQ or SEQ_NAME is found, the Q_* and S_* are
 497     # tried in that order.
 498
 499     my ( $record,    # record
 500        ) = @_;
 501
 502     # Returns a tuple.
 503
 504     my ( $seq_name, $seq );
 505
 506     $seq_name = $record->{ "SEQ_NAME" } || $record->{ "Q_ID" }  || $record->{ "S_ID" };
 507     $seq      = $record->{ "SEQ" }      || $record->{ "Q_SEQ" } || $record->{ "S_SEQ" };
 508
 509     if ( defined $seq_name and defined $seq ) {
 510         return wantarray ? ( $seq_name, $seq ) : [ $seq_name, $seq ];
 511     } else {
 512         return;
 513     }
 514 }
 515
 516
 517 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<