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 );
  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 );
 175
 176     if ( $args )
 177     {
 178         while ( $lines->[ $i ] =~ /^(\w{2})\s+(.*)/ and $1 ne 'XX' and exists $args->{ $1 } )
 179         {
 180             if ( exists $ref{ $1 } ) {
 181                 $ref{ $1 } .= " " . $2;
 182             } else {
 183                 $ref{ $1 } = $2;
 184             }
 185
 186             $i++;
 187         }
 188     }
 189     else
 190     {
 191         while ( $lines->[ $i ] =~ /^(\w{2})\s+(.*)/ and $1 ne 'XX' )
 192         {
 193             if ( exists $ref{ $1 } ) {
 194                 $ref{ $1 } .= " " . $2;
 195             } else {
 196                 $ref{ $1 } = $2;
 197             }
 198
 199             $i++;
 200         }
 201     }
 202
 203     push @{ $hash->{ 'REF' } }, \%ref;
 204 }
 205
 206
 207 sub parse_feature_table
 208 {
 209     # Martin A. Hansen, June 2006.
 210
 211     # parses the feature table of a EMBL/GenBank/DDBJ entry.
 212     # parsing takes place in 4 steps. 1) the feature key is
 213     # located. 2) the locator is located taking into # consideration
 214     # that it may be split over multiple lines, which is dealt with
 215     # by counting the params that always are present in multiline
 216     # locators. 3) the locator is used to fetch the corresponding
 217     # sequence. 4) qualifier key/value pars are located again taking
 218     # into consideration multiline values, which are dealt with by
 219     # keeping track of the "-count (value-less qualifers are also
 220     # included). only feature keys and qualifers defined in the
 221     # argument hash are returned.
 222
 223     my ( $ft,     # feature table
 224          $seq,    # entry sequence
 225          $args,   # argument hash
 226        ) = @_;
 227
 228     # returns data structure
 229
 230     my ( @lines, $key_regex, $i, $p, $q, %key_hash, $key, $locator, %qual_hash, $qual_name, $qual_val, $subseq );
 231
 232     @lines = split "\n", $ft;
 233
 234     $key_regex = "[A-Za-z0-9_']+";   # this regex should match every possible feature key (gene, misc_feature, 5'UTR ...)
 235
 236     $i = 0;
 237
 238     while ( $lines[ $i ] )
 239     {
 240         if ( $lines[ $i ] =~ /^($key_regex)\s+(.+)/ )
 241         {
 242             $key     = $1;
 243             $locator = $2;
 244
 245             undef %qual_hash;
 246
 247             # ---- getting locator
 248
 249             $p = 1;
 250
 251             if ( not balance_params( $locator ) )
 252             {
 253                 while ( not balance_params( $locator ) )
 254                 {
 255                     $locator .= $lines[ $i + $p ];
 256                     $p++;
 257                 }
 258             }
 259
 260             push @{ $qual_hash{ "_locator" } }, $locator;
 261
 262             if ( $seq )
 263             {
 264                 # ---- getting subsequence
 265
 266                 $subseq = parse_locator( $locator, $seq );
 267
 268                 push @{ $qual_hash{ "_seq" } }, $subseq;
 269             }
 270
 271             # ----- getting qualifiers
 272
 273             while ( defined( $lines[ $i + $p ] ) and not $lines[ $i + $p ] =~ /^$key_regex/ )
 274             {
 275                 if ( $lines[ $i + $p ] =~ /^\// )
 276                 {
 277                     if ( $lines[ $i + $p ] =~ /^\/([^=]+)=(.*)$/ )
 278                     {
 279                         $qual_name = $1;
 280                         $qual_val  = $2;
 281                     }
 282                     elsif ( $lines[ $i + $p ] =~ /^\/(.*)$/ )
 283                     {
 284                         $qual_name = $1;
 285                         $qual_val  = "";
 286                     }
 287
 288                     # ----- getting qualifier value
 289
 290                     $q = 1;
 291
 292                     if ( not balance_quotes( $qual_val ) )
 293                     {
 294                         while ( not balance_quotes( $qual_val ) )
 295                         {
 296                             $qual_val .= " " . $lines[ $i + $p + $q ];
 297                             $q++;
 298                         }
 299                     }
 300
 301                     $qual_val =~ s/^"(.*)"$/$1/;
 302                     $qual_val =~ tr/ //d if $qual_name =~ /translation/i;
 303
 304                     if ( exists $args->{ "quals" } ) {
 305                         push @{ $qual_hash{ $qual_name } }, $qual_val if exists $args->{ "quals" }->{ $qual_name };
 306                     } else {
 307                         push @{ $qual_hash{ $qual_name } }, $qual_val;
 308                     }
 309                 }
 310
 311                 $p += $q;
 312             }
 313
 314             if ( scalar keys %qual_hash > 0 )
 315             {
 316                 if ( exists $args->{ "feats" } ) {
 317                     push @{ $key_hash{ $key } }, dclone \%qual_hash if exists $args->{ "feats" }->{ $key };
 318                 } else {
 319                     push @{ $key_hash{ $key } }, dclone \%qual_hash;
 320                 }
 321             }
 322         }
 323
 324         $i += $p;
 325     }
 326
 327     return wantarray ? %key_hash : \%key_hash;
 328 }
 329
 330
 331 sub parse_locator
 332 {
 333     # Martin A. Hansen, June 2006.
 334
 335     # uses recursion to parse a locator string from a feature
 336     # table and fetches the appropriate subsequence. the operators
 337     # join(), complement(), and order() are handled.
 338     # the locator string is broken into a comma separated lists, and
 339     # modified if the params donnot balance. otherwise the comma separated
 340     # list of ranges are stripped from operators, and the subsequence are
 341     # fetched and handled according to the operators.
 342     # SNP locators are also dealt with (single positions).
 343
 344     my ( $locator,   # locator string
 345          $seq,       # nucleotide sequence
 346          $subseq,    # subsequence being joined
 347          $join,      # join sequences
 348          $comp,      # complement sequence or not
 349          $order,     # order sequences
 350        ) = @_;
 351
 352     # returns string
 353
 354     my ( @intervals, $interval, $beg, $end, $newseq );
 355
 356     @intervals = split ",", $locator;
 357
 358     if ( not balance_params( $intervals[ 0 ] ) )                          # locator includes a join/comp/order of several ranges
 359     {
 360         if ( $locator =~ /^join\((.*)\)$/ )
 361         {
 362             $join = 1;
 363             $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 364         }
 365         elsif ( $locator =~ /^complement\((.*)\)$/ )
 366         {
 367             $comp = 1;
 368             $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 369
 370         }
 371         elsif ( $locator =~ /^order\((.*)\)$/ )
 372         {
 373             $order = 1;
 374             $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 375         }
 376     }
 377     else
 378     {
 379         foreach $interval ( @intervals )
 380         {
 381             if ( $interval =~ /^join\((.*)\)$/ )
 382             {
 383                 $join = 1;
 384                 $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 385             }
 386             elsif ( $interval =~ /^complement\((.*)\)$/ )
 387             {
 388                 $comp = 1;
 389                 $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 390
 391             }
 392             elsif ( $interval =~ /^order\((.*)\)$/ )
 393             {
 394                 $order = 1;
 395                 $subseq = parse_locator( $1, $seq, $subseq, $join, $comp, $order );
 396             }
 397             elsif ( $interval =~ /^[<>]?(\d+)[^\d]+(\d+)$/ )
 398             {
 399                 $beg = $1;
 400                 $end = $2;
 401
 402                 $newseq = substr $seq, $beg - 1, $end - $beg + 1;
 403
 404                 $newseq = Maasha::Seq::dna_revcomp( $newseq ) if $comp;
 405
 406                 if ( $order ) {
 407                     $subseq .= " " . $newseq;
 408                 } else {
 409                     $subseq .= $newseq;
 410                 }
 411             }
 412             elsif ( $interval =~ /^(\d+)$/ )
 413             {
 414                 $beg = $1;
 415
 416                 $newseq = substr $seq, $beg - 1, 1 ;
 417
 418                 $newseq = Maasha::Seq::dna_revcomp( $newseq ) if $comp;
 419
 420                 if ( $order ) {
 421                     $subseq .= " " . $newseq;
 422                 } else {
 423                     $subseq .= $newseq;
 424                 }
 425             }
 426             else
 427             {
 428                 warn qq(WARNING: Could not match locator -> $locator\n);
 429                 # die qq(ERROR: Could not match locator -> $locator\n);
 430                 $subseq .= "";
 431             }
 432         }
 433     }
 434
 435     return $subseq;
 436 }
 437
 438
 439 sub balance_params
 440 {
 441     # Martin A. Hansen, June 2006.
 442
 443     # given a string checks if left and right params
 444     # balances. returns 1 if balanced, else 0.
 445
 446     my ( $string,   # string to check
 447        ) = @_;
 448
 449     # returns boolean
 450
 451     my ( $param_count );
 452
 453     $param_count  = 0;
 454     $param_count += $string =~ tr/(//;
 455     $param_count -= $string =~ tr/)//;
 456
 457     if ( $param_count == 0 ) {
 458         return 1;
 459     } else {
 460         return 0;
 461     }
 462 }
 463
 464
 465 sub balance_quotes
 466 {
 467     # Martin A. Hansen, June 2006.
 468
 469     # given a string checks if the number of double quotes
 470     # balances. returns 1 if balanced, else 0.
 471
 472     my ( $string,   # string to check
 473        ) = @_;
 474
 475     # returns boolean
 476
 477     my ( $quote_count );
 478
 479     $quote_count = $string =~ tr/"//;
 480
 481     if ( $quote_count % 2 == 0 ) {
 482         return 1;
 483     } else {
 484         return 0;
 485     }
 486 }
 487
 488
 489 sub embl2biopieces
 490 {
 491     # Martin A. Hansen, July 2009.
 492
 493     # Given a complete EMBL entry and an option hash,
 494     # configure the arguments for the EMBL parser so
 495     # only wanted information is retrieved and returned
 496     # as a Biopiece record.
 497
 498     my ( $entry,     # EMBL entry to parse
 499          $options,   # Biopiece options
 500        ) = @_;
 501
 502     # Returns a hashref.
 503
 504     my ( %args, $data, $record, $key_record, $key, $feat, $feat2, $qual, $qual_val, @records, $no_seq );
 505
 506     local $Data::Dumper::Terse  = 1;
 507     local $Data::Dumper::Indent = 0;
 508
 509     map { $args{ 'keys' }{ $_ } = 1 }  @{ $options->{ 'keys' } };
 510     map { $args{ 'feats' }{ $_ } = 1 } @{ $options->{ 'features' } };
 511     map { $args{ 'quals' }{ $_ } = 1 } @{ $options->{ 'qualifiers' } };
 512
 513     if ( @{ $options->{ 'features' } } > 0 or @{ $options->{ 'qualifiers' } } > 0 ) {
 514         $args{ 'keys' }{ 'FT' } = 1;
 515     }
 516
 517     if ( ( $args{ 'feats' } and $args{ 'feats' }{ 'SEQ' } ) or ( $args{ 'quals' } and $args{ 'quals' }{ 'SEQ' } ) )
 518     {
 519         $no_seq = 1 if not $args{ 'keys' }{ 'SEQ' };
 520
 521         $args{ 'keys' }{ 'SEQ' } = 1;
 522         delete $args{ 'feats' }{ 'SEQ' } if $args{ 'feats' };
 523         delete $args{ 'quals' }{ 'SEQ' } if $args{ 'quals' };
 524     }
 525
 526     $data = parse_embl_entry( $entry, \%args );
 527
 528     # print Dumper( $data );
 529
 530     foreach $key ( keys %{ $data } )
 531     {
 532         if ( $key eq 'SEQ' and $no_seq ) {
 533             next;
 534         }
 535
 536         if ( $key eq 'REF' ) {
 537             $key_record->{ $key } = Dumper( $data->{ $key } );
 538         }
 539
 540         if ( $key ne 'FT' and $key ne 'REF' ) {
 541             $key_record->{ $key } = $data->{ $key };
 542         }
 543     }
 544
 545     if ( exists $data->{ 'FT' } )
 546     {
 547         $record->{ 'FT' } = Dumper( $data->{ 'FT' } );
 548
 549         map { $record->{ $_ } = $key_record->{ $_ } } keys %{ $key_record };
 550
 551         push @records, $record;
 552     }
 553     else
 554     {
 555         push @records, $key_record;
 556     }
 557
 558     return wantarray ? @records : \@records;
 559 }
 560
 561
 562 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<