code_perl/Maasha/EMBL.pm

   1 package Maasha::EMBL;
   2
   3 # Copyright (C) 2007 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 to parse EMBL records.
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use warnings;
  32 use strict;
  33 use Data::Dumper;
  34 use Storable qw( dclone );
  35 use Maasha::Common;
  36 use Maasha::Fasta;
  37 use Maasha::Calc;
  38 use Maasha::Seq;
  39 use vars qw ( @ISA @EXPORT );
  40
  41
  42 @ISA = qw( Exporter );
  43
  44
  45 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  46
  47
  48 sub get_embl_entry
  49 {
  50     # Martin A. Hansen, June 2006.
  51
  52     # Given a filehandle to an embl file,
  53     # fetches the next embl entry, and returns
  54     # this as a string with multiple lines.
  55
  56     my ( $fh,    # filehandle to embl file
  57        ) = @_;
  58
  59     # returns string
  60
  61     my ( $entry );
  62
  63     local $/ = "//\n";
  64
  65     $entry = <$fh>;
  66
  67     return $entry;
  68 }
  69
  70
  71 sub parse_embl_entry
  72 {
  73     # Martin A. Hansen, June 2006.
  74
  75     # given an embl entry extracts the keys
  76     # given in an argument hash. Special care
  77     # is taken to parse the feature table if
  78     # requested.
  79
  80     my ( $entry,   # embl entry
  81          $args,    # argument hash
  82        ) = @_;
  83
  84     # returns data structure
  85
  86     my ( @lines, $i, %hash, $ft, $seq, $key, $val, $ref );
  87
  88     @lines = split "\n", $entry;
  89
  90     $i = 0;
  91
  92     while ( $i < @lines )
  93     {
  94         if ( exists $args->{ "keys" } )
  95         {
  96             $args->{ "keys" }->{ "RN" } = 1 if grep { $_ =~ /^R/ } keys %{ $args->{ "keys" } };
  97
  98             if ( $lines[ $i ] =~ /^(\w{2})\s+(.*)/ and exists $args->{ "keys" }->{ $1 } )
  99             {
 100                 $key = $1;
 101                 $val = $2;
 102
 103                 if ( $key =~ /RN|RX|RP|RG|RA|RT|RL/ ) {
 104                      add_ref( \%hash, \@lines, $i, $args->{ "keys" } ) if $key eq "RN";
 105                 } elsif ( exists $hash{ $key } and $key eq "FT" ) {
 106                     $hash{ $key } .= "\n" . $val;
 107                 } elsif ( exists $hash{ $1 } ) {
 108                     $hash{ $key } .= " " . $val;
 109                 } else {
 110                     $hash{ $key } = $val;
 111                 }
 112             }
 113             elsif ( $lines[ $i ] =~ /^\s+(.*)\s+\d+$/ and exists $args->{ "keys" }->{ "SEQ" } )
 114             {
 115                 $seq .= $1;
 116             }
 117         }
 118         else
 119         {
 120             if ( $lines[ $i ] =~ /^(\w{2})\s+(.*)/ )
 121             {
 122                 $key = $1;
 123                 $val = $2;
 124
 125                 if ( $key =~ /RN|RX|RP|RG|RA|RT|RL/ ) {
 126                     add_ref( \%hash, \@lines, $i ) if $key eq "RN";
 127                 } elsif ( exists $hash{ $1 } and $key eq "FT" ) {
 128                     $hash{ $key } .= "\n" . $val;
 129                 } elsif ( exists $hash{ $key } ) {
 130                     $hash{ $key } .= " " . $val;
 131                 } else {
 132                     $hash{ $key } = $val;
 133                 }
 134             }
 135             elsif ( $lines[ $i ] =~ /^\s+(.*)\s+\d+$/ )
 136             {
 137                 $seq .= $1;
 138             }
 139         }
 140
 141         $i++;
 142     }
 143
 144     if ( $seq )
 145     {
 146         $seq =~ tr/ //d;
 147         $hash{ "SEQ" } = $seq;
 148     }
 149
 150     if ( exists $hash{ "FT" } )
 151     {
 152         $ft = parse_feature_table( $hash{ "FT" }, $seq, $args );
 153         $hash{ "FT" } = $ft;
 154     }
 155
 156     return wantarray ? %hash : \%hash;
 157 }
 158
 159
 160 sub add_ref
 161 {
 162     # Martin A. Hansen, August 2009.
 163
 164     # Add a EMBL reference.
 165
 166     my ( $hash,    # Parsed EMBL data
 167          $lines,   # EMBL entry lines
 168          $i,       # line index
 169          $args,    # hashref with keys to save - OPTIONAL
 170        ) = @_;
 171
 172     # Returns nothing.
 173
 174     my ( %ref, $key, $val );
 175
 176     if ( $args )
 177     {
 178         while ( $lines->[ $i ] =~ /^(\w{2})\s+(.*)/ )
 179         {
 180             $key = $1;
 181             $val = $2;
 182
 183             last if $key eq "XX";
 184
 185             if ( exists $args->{ $key } )
 186             {
 187                 if ( exists $ref{ $key } ) {
 188                     $ref{ $key } .= " " . $val;
 189                 } else {
 190                     $ref{ $key } = $val;
 191                 }
 192             }
 193
 194             $i++;
 195         }
 196     }
 197     else
 198     {
 199         while ( $lines->[ $i ] =~ /^(\w{2})\s+(.*)/ and $1 ne 'XX' )
 200         {
 201             if ( exists $ref{ $1 } ) {
 202                 $ref{ $1 } .= " " . $2;
 203             } else {
 204                 $ref{ $1 } = $2;
 205             }
 206
 207             $i++;
 208         }
 209     }
 210
 211     push @{ $hash->{ 'REF' } }, \%ref;
 212 }
 213
 214
 215 sub parse_feature_table
 216 {
 217     # Martin A. Hansen, June 2006.
 218
 219     # parses the feature table of a EMBL/GenBank/DDBJ entry.
 220     # parsing takes place in 4 steps. 1) the feature key is
 221     # located. 2) the locator is located taking into # consideration
 222     # that it may be split over multiple lines, which is dealt with
 223     # by counting the params that always are present in multiline
 224     # locators. 3) the locator is used to fetch the corresponding
 225     # sequence. 4) qualifier key/value pars are located again taking
 226     # into consideration multiline values, which are dealt with by
 227     # keeping track of the "-count (value-less qualifers are also
 228     # included). only feature keys and qualifers defined in the
 229     # argument hash are returned.
 230
 231     my ( $ft,     # feature table
 232          $seq,    # entry sequence
 233          $args,   # argument hash
 234        ) = @_;
 235
 236     # returns data structure
 237
 238     my ( @lines, $key_regex, $i, $p, $q, %key_hash, $key, $locator, %qual_hash, $qual_name, $qual_val, $subseq );
 239
 240     @lines = split "\n", $ft;
 241
 242     $key_regex = "[A-Za-z0-9_']+";   # this regex should match every possible feature key (gene, misc_feature, 5'UTR ...)
 243
 244     $i = 0;
 245
 246     while ( $lines[ $i ] )
 247     {
 248         if ( $lines[ $i ] =~ /^($key_regex)\s+(.+)/ )
 249         {
 250             $key     = $1;
 251             $locator = $2;
 252
 253             undef %qual_hash;
 254
 255             # ---- getting locator
 256
 257             $p = 1;
 258
 259             if ( not balance_params( $locator ) )
 260             {
 261                 while ( not balance_params( $locator ) )
 262                 {
 263                     $locator .= $lines[ $i + $p ];
 264                     $p++;
 265                 }
 266             }
 267
 268             push @{ $qual_hash{ "_locator" } }, $locator;
 269
 270             if ( $seq )
 271             {
 272                 # ---- getting subsequence
 273
 274                 $subseq = parse_locator( $locator, $seq );
 275
 276                 push @{ $qual_hash{ "_seq" } }, $subseq;
 277             }
 278
 279             # ----- getting qualifiers
 280
 281             while ( defined( $lines[ $i + $p ] ) and not $lines[ $i + $p ] =~ /^$key_regex/ )
 282             {
 283                 if ( $lines[ $i + $p ] =~ /^\// )
 284                 {
 285                     if ( $lines[ $i + $p ] =~ /^\/([^=]+)=(.*)$/ )
 286                     {
 287                         $qual_name = $1;
 288                         $qual_val  = $2;
 289                     }
 290                     elsif ( $lines[ $i + $p ] =~ /^\/(.*)$/ )
 291                     {
 292                         $qual_name = $1;
 293                         $qual_val  = "";
 294                     }
 295
 296                     # ----- getting qualifier value
 297
 298                     $q = 1;
 299
 300                     if ( not balance_quotes( $qual_val ) )
 301                     {
 302                         while ( not balance_quotes( $qual_val ) )
 303                         {
 304                             $qual_val .= " " . $lines[ $i + $p + $q ];
 305                             $q++;
 306                         }
 307                     }
 308
 309                     $qual_val =~ s/^"(.*)"$/$1/;
 310                     $qual_val =~ tr/ //d if $qual_name =~ /translation/i;
 311
 312                     if ( exists $args->{ "quals" } ) {
 313                         push @{ $qual_hash{ $qual_name } }, $qual_val if exists $args->{ "quals" }->{ $qual_name };
 314                     } else {
 315                         push @{ $qual_hash{ $qual_name } }, $qual_val;
 316                     }
 317                 }
 318
 319                 $p += $q;
 320             }
 321
 322             if ( scalar keys %qual_hash > 0 )
 323             {
 324                 if ( exists $args->{ "feats" } ) {
 325                     push @{ $key_hash{ $key } }, dclone \%qual_hash if exists $args->{ "feats" }->{ $key };
 326                 } else {
 327                     push @{ $key_hash{ $key } }, dclone \%qual_hash;
 328                 }
 329             }
 330         }
 331
 332         $i += $p;
 333     }
 334
 335     return wantarray ? %key_hash : \%key_hash;
 336 }
 337
 338
 339 sub parse_locator
 340 {
 341     # Martin A. Hansen, June 2006.
 342
 343     # uses recursion to parse a locator string from a feature
 344     # table and fetches the appropriate subsequence. the operators
 345     # join(), complement(), and order() are handled.
 346     # the locator string is broken into a comma separated lists, and
 347     # modified if the params donnot balance. otherwise the comma separated
 348     # list of ranges are stripped from operators, and the subsequence are
 349     # fetched and handled according to the operators.
 350     # SNP locators are also dealt with (single positions).
 351
 352     my ( $locator,   # locator string
 353          $seq,       # nucleotide sequence
 354          $subseq,    # subsequence being joined
 355          $join,      # join sequences
 356          $comp,      # complement sequence or not
 357          $order,     # order sequences
 358        ) = @_;
 359
 360     # returns string
 361
 362     my ( @intervals, $interval, $beg, $end, $newseq );
 363
 364     @intervals = split ",", $locator;
 365
 366     if ( not balance_params( $intervals[ 0 ] ) )                          # locator includes a join/comp/order of several ranges
 367     {
 368         if ( $locator =~ /^join\((.*)\)$/ )
 369         {
 370             $join = 1;
 371             $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 372         }
 373         elsif ( $locator =~ /^complement\((.*)\)$/ )
 374         {
 375             $comp = 1;
 376             $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 377
 378         }
 379         elsif ( $locator =~ /^order\((.*)\)$/ )
 380         {
 381             $order = 1;
 382             $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 383         }
 384     }
 385     else
 386     {
 387         foreach $interval ( @intervals )
 388         {
 389             if ( $interval =~ /^join\((.*)\)$/ )
 390             {
 391                 $join = 1;
 392                 $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 393             }
 394             elsif ( $interval =~ /^complement\((.*)\)$/ )
 395             {
 396                 $comp = 1;
 397                 $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 398
 399             }
 400             elsif ( $interval =~ /^order\((.*)\)$/ )
 401             {
 402                 $order = 1;
 403                 $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 404             }
 405             elsif ( $interval =~ /^[<>]?(\d+)[^\d]+(\d+)$/ )
 406             {
 407                 $beg = $1;
 408                 $end = $2;
 409
 410                 $newseq = substr $seq, $beg - 1, $end - $beg + 1;
 411
 412                 $newseq = Maasha::Seq::dna_revcomp( $newseq ) if $comp;
 413
 414                 if ( $order ) {
 415                     $subseq .= " " . $newseq;
 416                 } else {
 417                     $subseq .= $newseq;
 418                 }
 419             }
 420             elsif ( $interval =~ /^(\d+)$/ )
 421             {
 422                 $beg = $1;
 423
 424                 $newseq = substr $seq, $beg - 1, 1 ;
 425
 426                 $newseq = Maasha::Seq::dna_revcomp( $newseq ) if $comp;
 427
 428                 if ( $order ) {
 429                     $subseq .= " " . $newseq;
 430                 } else {
 431                     $subseq .= $newseq;
 432                 }
 433             }
 434             else
 435             {
 436                 warn qq(WARNING: Could not match locator -> $locator\n);
 437                 # die qq(ERROR: Could not match locator -> $locator\n);
 438                 $subseq .= "";
 439             }
 440         }
 441     }
 442
 443     return $subseq;
 444 }
 445
 446
 447 sub balance_params
 448 {
 449     # Martin A. Hansen, June 2006.
 450
 451     # given a string checks if left and right params
 452     # balances. returns 1 if balanced, else 0.
 453
 454     my ( $string,   # string to check
 455        ) = @_;
 456
 457     # returns boolean
 458
 459     my ( $param_count );
 460
 461     $param_count  = 0;
 462     $param_count += $string =~ tr/(//;
 463     $param_count -= $string =~ tr/)//;
 464
 465     if ( $param_count == 0 ) {
 466         return 1;
 467     } else {
 468         return 0;
 469     }
 470 }
 471
 472
 473 sub balance_quotes
 474 {
 475     # Martin A. Hansen, June 2006.
 476
 477     # given a string checks if the number of double quotes
 478     # balances. returns 1 if balanced, else 0.
 479
 480     my ( $string,   # string to check
 481        ) = @_;
 482
 483     # returns boolean
 484
 485     my ( $quote_count );
 486
 487     $quote_count = $string =~ tr/"//;
 488
 489     if ( $quote_count % 2 == 0 ) {
 490         return 1;
 491     } else {
 492         return 0;
 493     }
 494 }
 495
 496
 497 sub embl2biopieces
 498 {
 499     # Martin A. Hansen, July 2009.
 500
 501     # Given a complete EMBL entry and an option hash,
 502     # configure the arguments for the EMBL parser so
 503     # only wanted information is retrieved and returned
 504     # as a Biopiece record.
 505
 506     my ( $entry,     # EMBL entry to parse
 507          $options,   # Biopiece options
 508        ) = @_;
 509
 510     # Returns a hashref.
 511
 512     my ( %args, $data, $record, $key_record, $key, $feat, $feat2, $qual, $qual_val, @records, $no_seq );
 513
 514     local $Data::Dumper::Terse  = 1;
 515     local $Data::Dumper::Indent = 0;
 516
 517     map { $args{ 'keys' }{ $_ } = 1 }  @{ $options->{ 'keys' } };
 518     map { $args{ 'feats' }{ $_ } = 1 } @{ $options->{ 'features' } };
 519     map { $args{ 'quals' }{ $_ } = 1 } @{ $options->{ 'qualifiers' } };
 520
 521     if ( @{ $options->{ 'features' } } > 0 or @{ $options->{ 'qualifiers' } } > 0 ) {
 522         $args{ 'keys' }{ 'FT' } = 1;
 523     }
 524
 525     if ( ( $args{ 'feats' } and $args{ 'feats' }{ 'SEQ' } ) or ( $args{ 'quals' } and $args{ 'quals' }{ 'SEQ' } ) )
 526     {
 527         $no_seq = 1 if not $args{ 'keys' }{ 'SEQ' };
 528
 529         $args{ 'keys' }{ 'SEQ' } = 1;
 530         delete $args{ 'feats' }{ 'SEQ' } if $args{ 'feats' };
 531         delete $args{ 'quals' }{ 'SEQ' } if $args{ 'quals' };
 532     }
 533
 534     $data = parse_embl_entry( $entry, \%args );
 535
 536     # print Dumper( $data );
 537
 538     foreach $key ( keys %{ $data } )
 539     {
 540         if ( $key eq 'SEQ' and $no_seq ) {
 541             next;
 542         }
 543
 544         if ( $key eq 'REF' ) {
 545             $key_record->{ $key } = Dumper( $data->{ $key } );
 546         }
 547
 548         if ( $key ne 'FT' and $key ne 'REF' ) {
 549             $key_record->{ $key } = $data->{ $key };
 550         }
 551     }
 552
 553     if ( exists $data->{ 'FT' } )
 554     {
 555         $record->{ 'FT' } = Dumper( $data->{ 'FT' } );
 556
 557         map { $record->{ $_ } = $key_record->{ $_ } } keys %{ $key_record };
 558
 559         push @records, $record;
 560     }
 561     else
 562     {
 563         push @records, $key_record;
 564     }
 565
 566     return wantarray ? @records : \@records;
 567 }
 568
 569
 570 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<