code_perl/Maasha/Align.pm

   1 package Maasha::Align;
   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 perform and print pairwise and multiple alignments
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use strict;
  32 use Data::Dumper;
  33 use IPC::Open2;
  34 use Maasha::Common;
  35 use Maasha::Fasta;
  36 use Maasha::Calc;
  37 use Maasha::Seq;
  38 use vars qw ( @ISA @EXPORT );
  39
  40 use constant {
  41     HEAD  => 0,
  42     SEQ   => 1,
  43 };
  44
  45 @ISA = qw( Exporter );
  46
  47
  48 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  49
  50
  51 sub align
  52 {
  53     # Martin A. Hansen, August 2007.
  54
  55     # Aligns a given list of FASTA entries and returns a
  56     # list of aligned sequences as FASTA entries.
  57     # (currently uses Muscle, but other align engines can
  58     # be used with a bit of tweaking).
  59
  60     my ( $entries,   # Fasta entries
  61          $args,      # additional alignment program specific arguments - OPTIONAL
  62        ) = @_;
  63
  64     # Returns a list.
  65
  66     my ( @aligned_entries, $muscle_args );
  67
  68     $muscle_args  = "-quiet";
  69     $muscle_args .= $args if $args;
  70
  71     @aligned_entries = &align_muscle( $entries, $muscle_args );
  72
  73     return wantarray ? @aligned_entries : \@aligned_entries;
  74 }
  75
  76
  77 sub align_muscle
  78 {
  79     # Martin A. Hansen, June 2007.
  80
  81     # Aligns a given list of FASTA entries using Muscle.
  82     # Returns a list of aligned sequences as FASTA entries.
  83
  84     my ( $entries,   # FASTA entries
  85          $args,      # additional Muscle arguments - OPTIONAL
  86        ) = @_;
  87
  88     # Returns a list.
  89
  90     my ( $pid, $fh_in, $fh_out, $cmd, $entry, @aligned_entries );
  91
  92     $cmd  = "muscle";
  93     $cmd .= " " . $args if $args;
  94
  95     $pid = open2( $fh_out, $fh_in, $cmd );
  96
  97     map { &Maasha::Fasta::put_entry( $_, $fh_in ) } @{ $entries };
  98
  99     close $fh_in;
 100
 101     while ( $entry = &Maasha::Fasta::get_entry( $fh_out ) ) {
 102         push @aligned_entries, $entry;
 103     }
 104
 105     close $fh_out;
 106
 107     waitpid $pid, 0;
 108
 109     return wantarray ? @aligned_entries : \@aligned_entries;
 110 }
 111
 112
 113 sub align_print_pairwise
 114 {
 115     # Martin A. Hansen, June 2007.
 116
 117     # Prints a given pairwise alignment in FASTA format.
 118
 119     my ( $entry1,   # first entry
 120          $entry2,   # second entry
 121          $fh,       # output filehandle - OPTIONAL
 122          $wrap,     # wrap width        - OPTIONAL
 123        ) = @_;
 124
 125     # returns nothing
 126
 127     my ( @entries, $ruler1, $ruler2, $pins );
 128
 129     $ruler1 = &align_ruler( $entry1, 1 );
 130     $ruler2 = &align_ruler( $entry2, 1 );
 131     $pins   = &align_pins( $entry1, $entry2 );
 132
 133     push @entries, $ruler1, $entry1, $pins, $entry2, $ruler2;
 134
 135     &align_print( \@entries, $fh, $wrap );
 136 }
 137
 138
 139 sub align_print_multi
 140 {
 141     # Martin A. Hansen, June 2007.
 142
 143     # Prints a given multiple alignment in FASTA format.
 144
 145     my ( $entries,     # list of aligned FASTA entries
 146          $fh,          # output filehandle    - OPTIONAL
 147          $wrap,        # wrap width           - OPTIONAL
 148          $no_ruler,    # omit ruler flag      - OPTIONAL
 149          $no_cons,     # omit consensus flag  - OPTIONAL
 150        ) = @_;
 151
 152     # returns nothing
 153
 154     my ( @entries, $ruler, $consensus );
 155
 156     $ruler     = &align_ruler( $entries->[ 0 ] );
 157     $consensus = &align_consensus( $entries ) if not $no_cons;
 158
 159     unshift @{ $entries }, $ruler if not $no_ruler;
 160     push    @{ $entries }, $consensus;
 161
 162     &align_print( $entries, $fh, $wrap );
 163 }
 164
 165
 166 sub align_print
 167 {
 168     # Martin A. Hansen, June 2007.
 169
 170     # Prints an alignment.
 171
 172     my ( $entries,   # Alignment as FASTA entries
 173          $fh,        # output filehandle  - OPTIONAL
 174          $wrap,      # wrap alignment     - OPTIONAL
 175        ) = @_;
 176
 177     # returns nothing
 178
 179     my ( $max, $blocks, $block, $entry );
 180
 181     $max = 0;
 182
 183     map { $max = length $_->[ HEAD ] if length $_->[ HEAD ] > $max } @{ $entries };
 184
 185     $blocks = &align_wrap( $entries, $wrap );
 186
 187     foreach $block ( @{ $blocks } )
 188     {
 189         foreach $entry ( @{ $block } )
 190         {
 191             $entry->[ HEAD ] =~ s/stats|ruler|consensus//;
 192
 193             if ( $fh ) {
 194                 print $fh $entry->[ HEAD ], " " x ( $max + 3 - length $entry->[ HEAD ] ), $entry->[ SEQ ], "\n";
 195             } else {
 196                 print $entry->[ HEAD ], " " x ( $max + 3 - length $entry->[ HEAD ] ), $entry->[ SEQ ], "\n";
 197             }
 198         }
 199     }
 200 }
 201
 202
 203 sub align_wrap
 204 {
 205     # Martin A. Hansen, October 2005.
 206
 207     # Given a set of fasta entries wraps these
 208     # according to a given width.
 209
 210     my ( $entries,     # list of fasta_entries
 211          $wrap,        # wrap width - OPTIONAL
 212        ) = @_;
 213
 214     # returns AoA
 215
 216     my ( $ruler, $i, $c, @lines, @blocks );
 217
 218     $wrap ||= 999999999;
 219
 220     $i = 0;
 221
 222     while ( $i < length $entries->[ 0 ]->[ SEQ ] )
 223     {
 224         undef @lines;
 225
 226         for ( $c = 0; $c < @{ $entries }; $c++ )
 227         {
 228             if ( $entries->[ $c ]->[ HEAD ] eq "ruler" )
 229             {
 230                 $ruler = substr $entries->[ $c ]->[ SEQ ], $i, $wrap;
 231
 232                 if ( $ruler =~ /^(\d+)/ ) {
 233                     $ruler =~ s/^($1)/' 'x(length $1)/e;
 234                 }
 235
 236                 if ( $ruler =~ /(\d+)$/ ) {
 237                     $ruler =~ s/($1)$/' 'x(length $1)/e;
 238                 }
 239
 240                 push @lines, [ "ruler", $ruler ];
 241             }
 242             else
 243             {
 244                 push @lines, [ $entries->[ $c ]->[ HEAD ], substr $entries->[ $c ]->[ SEQ ], $i, $wrap ];
 245             }
 246         }
 247
 248         push @blocks, [ @lines ];
 249
 250         $i += $wrap;
 251     }
 252
 253     return wantarray ? @blocks: \@blocks;
 254 }
 255
 256
 257 sub align_pins
 258 {
 259     # Martin A. Hansen, June 2007.
 260
 261     # Given two aligned FASTA entries, generates an entry with pins.
 262
 263     my ( $entry1,   # first entry
 264          $entry2,   # second entry
 265          $type,     # residue type - OPTIONAL
 266        ) = @_;
 267
 268     # returns tuple
 269
 270     my ( $blosum, $i, $char1, $char2, $pins );
 271
 272     $type ||= &Maasha::Seq::seq_guess_type( $entry1->[ SEQ ] );
 273
 274     $blosum = &blosum_read() if $type =~ /protein/;
 275
 276     for ( $i = 0; $i < length $entry1->[ SEQ ]; $i++ )
 277     {
 278         $char1 = substr $entry1->[ SEQ ], $i, 1;
 279         $char2 = substr $entry2->[ SEQ ], $i, 1;
 280
 281         if ( $blosum and $char1 eq $char2 ) {
 282             $pins .= $char1;
 283         } elsif ( $char1 eq $char2 ) {
 284             $pins .= "|";
 285         } elsif ( $blosum and $blosum->{ $char1 }->{ $char2 } > 0 ) {
 286             $pins .= "+";
 287         } else {
 288             $pins .= " ";
 289         }
 290     }
 291
 292     return wantarray ? ( "consensus", $pins ) : [ "consensus", $pins ];
 293 }
 294
 295
 296 sub align_ruler
 297 {
 298     # Martin A. Hansen, February 2007;
 299
 300     # Gererates a ruler for a given FASTA entry (with indels).
 301
 302     my ( $entry,        # FASTA entry
 303          $count_gaps,   # flag for counting indels in pairwise alignments.
 304        ) = @_;
 305
 306     # Returns tuple
 307
 308     my ( $i, $char, $skip, $count, $gap, $tics );
 309
 310     $char = "";
 311     $gap  = 0;
 312     $i    = 1;
 313
 314     while ( $i <= length $entry->[ SEQ ] )
 315     {
 316         $char = substr( $entry->[ SEQ ], $i - 1, 1 ) if $count_gaps;
 317
 318         $gap++ if $char eq "-";
 319
 320         if ( $skip )
 321         {
 322             $skip--;
 323         }
 324         else
 325         {
 326             $count = $i - $gap;
 327             $count = 1 if $char eq "-";
 328
 329             if ( $count % 100 == 0 )
 330             {
 331                 if ( $count + length( $count ) >= length $entry->[ SEQ ] )
 332                 {
 333                     $tics .= "|";
 334                 }
 335                 else
 336                 {
 337                     $tics .= "|" . $count;
 338                     $skip = length $count;
 339                 }
 340             }
 341             elsif ( $count % 50 == 0 ) {
 342                 $tics .= ":";
 343             } elsif ( $count % 10 == 0 ) {
 344                 $tics .= ".";
 345             } else {
 346                 $tics .= " ";
 347             }
 348         }
 349
 350         $i++;
 351     }
 352
 353     return wantarray ? ( "ruler", $tics ) : [ "ruler", $tics ];
 354 }
 355
 356
 357 sub align_consensus
 358 {
 359     # Martin A. Hansen, June 2006.
 360
 361     # Given an alignment as a list of FASTA entries,
 362     # generates a consensus sequences based on the
 363     # entropies for each column similar to the way
 364     # a sequence logo i calculated. Returns the
 365     # consensus sequence as a FASTA entry.
 366
 367     my ( $entries,   # list of aligned FASTA entries
 368          $type,      # residue type       - OPTIONAL
 369          $min_sim,   # minimum similarity - OPTIONAL
 370        ) = @_;
 371
 372     # Returns tuple
 373
 374     my ( $bit_max, $data, $pos, $char, $score, $entry );
 375
 376     $type    ||= &Maasha::Seq::seq_guess_type( $entries->[ 0 ]->[ SEQ ] );
 377     $min_sim ||= 50;
 378
 379     if ( $type =~ /protein/ ) {
 380         $bit_max = 4;
 381     } else {
 382         $bit_max = 2;
 383     }
 384
 385     $data = &Maasha::Seq::seqlogo_calc( $bit_max, $entries );
 386
 387     foreach $pos ( @{ $data } )
 388     {
 389         if ( $pos->[ -1 ] )
 390         {
 391             ( $char, $score ) = @{ $pos->[ -1 ] };
 392
 393             if ( ( $score / $bit_max ) * 100 >= $min_sim ) {
 394                 $entry->[ SEQ ] .= $char;
 395             } else {
 396                 $entry->[ SEQ ] .= "-";
 397             }
 398         }
 399         else
 400         {
 401             $entry->[ SEQ ] .= "-";
 402         }
 403     }
 404
 405     $entry->[ HEAD ] = "Consensus: $min_sim%";
 406
 407     return wantarray ? @{ $entry } : $entry;
 408 }
 409
 410
 411 sub align_sim_global
 412 {
 413     # Martin A. Hansen, June 2007.
 414
 415     # Calculate the global similarity of two aligned entries
 416     # The similarity is calculated as the number of matching
 417     # residues divided by the length of the shortest sequence.
 418
 419     my ( $entry1,   # first  aligned entry
 420          $entry2,   # second aligned entry
 421        ) = @_;
 422
 423     # returns float
 424
 425     my ( $seq1, $seq2, $len1, $len2, $i, $match_tot, $min, $sim );
 426
 427     $seq1 = $entry1->[ SEQ ];
 428     $seq2 = $entry2->[ SEQ ];
 429
 430     # $seq1 =~ tr/-//d;
 431     # $seq2 =~ tr/-//d;
 432
 433     $seq1 =~ s/^-*//;
 434     $seq2 =~ s/^-*//;
 435     $seq1 =~ s/-*$//;
 436     $seq2 =~ s/-*$//;
 437
 438     $len1 = length $seq1;
 439     $len2 = length $seq2;
 440
 441     return 0 if $len1 == 0 or $len2 == 0;
 442
 443     $match_tot = 0;
 444
 445     for ( $i = 0; $i < $len1; $i++ ) {
 446         $match_tot++ if substr( $entry1->[ SEQ ], $i, 1 ) eq substr( $entry2->[ SEQ ], $i, 1 );
 447     }
 448
 449     $min = &Maasha::Calc::min( $len1, $len2 );
 450
 451     $sim = sprintf( "%.2f", ( $match_tot / $min ) * 100 );
 452
 453     return $sim;
 454 }
 455
 456
 457 sub align_tile
 458 {
 459     # Martin A. Hansen, February 2008.
 460
 461     # Tile a list of query sequences agains a reference sequence,
 462     # using pairwise alignments. The result is returned as a list of
 463     # aligned FASTA entries.
 464
 465     my ( $ref_entry,   # reference entry as [ HEAD, SEQ ] tuple
 466          $q_entries,   # list of [ HEAD, SEQ ] tuples
 467          $args,        # argument hash
 468        ) = @_;
 469
 470     # Returns a list.
 471
 472     my ( $entry, $seq1, $seq2, $type, $align1, $align2, $sim1, $sim2, $gaps, @entries );
 473
 474     $args->{ "identity" } ||= 70;
 475
 476     foreach $entry ( @{ $q_entries } )
 477     {
 478         $seq1 = $entry->[ SEQ ];
 479
 480         $type = &Maasha::Seq::seq_guess_type( $seq1 );
 481
 482         if ( $type eq "rna" ) {
 483             $seq2 = &Maasha::Seq::rna_revcomp( $seq1 );
 484         } elsif ( $type eq "dna" ) {
 485             $seq2 = &Maasha::Seq::dna_revcomp( $seq1 );
 486         } else {
 487             &Maasha::Common::error( qq(Bad sequence type->$type) );
 488         }
 489
 490         $align1 = &Maasha::Align::align_muscle( [ $ref_entry, [ $entry->[ HEAD ] . "_+", $seq1 ] ], "-quiet -maxiters 1" );
 491         $align2 = &Maasha::Align::align_muscle( [ $ref_entry, [ $entry->[ HEAD ] . "_-", $seq2 ] ], "-quiet -maxiters 1" );
 492
 493         if ( $args->{ "supress_indels" } )
 494         {
 495             &align_supress_indels( $align1 );
 496             &align_supress_indels( $align2 );
 497         }
 498
 499         $sim1 = &Maasha::Align::align_sim_global( $align1->[ 0 ], $align1->[ 1 ] );
 500         $sim2 = &Maasha::Align::align_sim_global( $align2->[ 0 ], $align2->[ 1 ] );
 501
 502         if ( $sim1 < $args->{ "identity" } and $sim2 < $args->{ "identity" } )
 503         {
 504             # do nothing
 505         }
 506         elsif ( $sim1 > $sim2 )
 507         {
 508             $gaps = $align1->[ 0 ]->[ SEQ ] =~ tr/-//;
 509
 510             $align1->[ 1 ]->[ SEQ ] =~ s/-{$gaps}$// if $gaps;
 511
 512             $entry->[ HEAD ] = "$align1->[ 1 ]->[ HEAD ]_$sim1";
 513             $entry->[ SEQ ]  = $align1->[ 1 ]->[ SEQ ];
 514
 515             push @entries, $entry;
 516         }
 517         else
 518         {
 519             $gaps = $align2->[ 0 ]->[ SEQ ] =~ tr/-//;
 520
 521             $align2->[ 1 ]->[ SEQ ] =~ s/-{$gaps}$// if $gaps;
 522
 523             $entry->[ HEAD ] = "$align2->[ 1 ]->[ HEAD ]_$sim2";
 524             $entry->[ SEQ ]  = $align2->[ 1 ]->[ SEQ ];
 525
 526             push @entries, $entry;
 527         }
 528     }
 529
 530     @entries = sort { $b->[ SEQ ] cmp $a->[ SEQ ] } @entries;
 531
 532     unshift @entries, $ref_entry;
 533
 534     return wantarray ? @entries : \@entries;
 535 }
 536
 537
 538 sub align_supress_indels
 539 {
 540     # Martin A. Hansen, June 2008.
 541
 542     # Given a pairwise alignment, removes
 543     # indels in the first sequence AND corresponding
 544     # sequence in the second.
 545
 546     my ( $align,   # pairwise alignment
 547        ) = @_;
 548
 549     # Returns nothing
 550
 551     my ( $count, $seq, $i );
 552
 553     $count = $align->[ 0 ]->[ SEQ ] =~ tr/-//;
 554
 555     if ( $count > 0 )
 556     {
 557         for ( $i = 0; $i < length $align->[ 0 ]->[ SEQ ]; $i++ )
 558         {
 559             if ( substr( $align->[ 0 ]->[ SEQ ], $i, 1 ) ne '-' ) {
 560                 $seq .= substr( $align->[ 1 ]->[ SEQ ], $i, 1 );
 561             }
 562
 563         }
 564
 565         $align->[ 0 ]->[ SEQ ] =~ tr/-//d;
 566         $align->[ 1 ]->[ SEQ ] = $seq;
 567     }
 568 }
 569
 570
 571 sub align_invert
 572 {
 573     # Martin A. Hansen, February 2008.
 574
 575     # Invert an alignment in such a way that only
 576     # residues differing from the first sequence (the reference sequence)
 577     # are shown. The matching sequence can either be lowercased (soft) or replaced
 578     # with _.
 579
 580     my ( $entries,   # list of FASTA entries.
 581          $soft,
 582        ) = @_;
 583
 584     # Returns nothing.
 585
 586     my ( $i, $c, $char1, $char2 );
 587
 588     map { $_->[ SEQ ] =~ tr/-/_/ } @{ $entries };
 589
 590     for ( $i = 0; $i < length $entries->[ 0 ]->[ SEQ ]; $i++ )
 591     {
 592         $char1 = uc substr $entries->[ 0 ]->[ SEQ ], $i, 1;
 593
 594         for ( $c = 1; $c < @{ $entries }; $c++ )
 595         {
 596             $char2 = uc substr $entries->[ $c ]->[ SEQ ], $i, 1;
 597
 598             if ( $char1 eq $char2 )
 599             {
 600                 if ( $soft ) {
 601                     substr $entries->[ $c ]->[ SEQ ], $i, 1, lc $char2;
 602                 } else {
 603                     substr $entries->[ $c ]->[ SEQ ], $i, 1, "-";
 604                 }
 605             }
 606         }
 607     }
 608 }
 609
 610
 611 sub blosum_read
 612 {
 613     # Martin A. Hansen, January 2006.
 614
 615     # this routine parses the BLOSUM62 matrix,
 616     # which is located in the __DATA__ section
 617
 618     # returns HoH
 619
 620     my ( @lines, @chars, $i, $c, @list, $HoH );
 621
 622     @lines = <DATA>;
 623     @lines = grep { $_ !~ /^$|^#/ } @lines;
 624
 625     @chars = split /\s+/, $lines[ 0 ];
 626
 627     $i = 1;
 628
 629     while( $lines[ $i ] )
 630     {
 631         last if $lines[ $i ] =~ /^__END__/;
 632
 633         @list = split /\s+/, $lines[ $i ];
 634
 635         for ( $c = 1; $c < @list; $c++ ) {
 636             $HoH->{ $list[ 0 ] }->{ $chars[ $c ] } = $list[ $c ];
 637         }
 638
 639         $i++;
 640     }
 641
 642     return wantarray ? %{ $HoH } : $HoH;
 643 }
 644
 645
 646 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 647
 648
 649 __DATA__
 650
 651
 652 #  Matrix made by matblas from blosum62.iij
 653 #  * column uses minimum score
 654 #  BLOSUM Clustered Scoring Matrix in 1/2 Bit Units
 655 #  Blocks Database = /data/blocks_5.0/blocks.dat
 656 #  Cluster Percentage: >= 62
 657 #  Entropy =   0.6979, Expected =  -0.5209
 658    A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X  *
 659 A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0 -2 -1  0 -4
 660 R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3 -1  0 -1 -4
 661 N -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3  3  0 -1 -4
 662 D -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3  4  1 -1 -4
 663 C  0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4
 664 Q -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2  0  3 -1 -4
 665 E -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4
 666 G  0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3 -1 -2 -1 -4
 667 H -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3  0  0 -1 -4
 668 I -1 -3 -3 -3 -1 -3 -3 -4 -3  4  2 -3  1  0 -3 -2 -1 -3 -1  3 -3 -3 -1 -4
 669 L -1 -2 -3 -4 -1 -2 -3 -4 -3  2  4 -2  2  0 -3 -2 -1 -2 -1  1 -4 -3 -1 -4
 670 K -1  2  0 -1 -3  1  1 -2 -1 -3 -2  5 -1 -3 -1  0 -1 -3 -2 -2  0  1 -1 -4
 671 M -1 -1 -2 -3 -1  0 -2 -3 -2  1  2 -1  5  0 -2 -1 -1 -1 -1  1 -3 -1 -1 -4
 672 F -2 -3 -3 -3 -2 -3 -3 -3 -1  0  0 -3  0  6 -4 -2 -2  1  3 -1 -3 -3 -1 -4
 673 P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4  7 -1 -1 -4 -3 -2 -2 -1 -2 -4
 674 S  1 -1  1  0 -1  0  0  0 -1 -2 -2  0 -1 -2 -1  4  1 -3 -2 -2  0  0  0 -4
 675 T  0 -1  0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1  1  5 -2 -2  0 -1 -1  0 -4
 676 W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1  1 -4 -3 -2 11  2 -3 -4 -3 -2 -4
 677 Y -2 -2 -2 -3 -2 -1 -2 -3  2 -1 -1 -2 -1  3 -3 -2 -2  2  7 -1 -3 -2 -1 -4
 678 V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4 -3 -2 -1 -4
 679 B -2 -1  3  4 -3  0  1 -1  0 -3 -4  0 -3 -3 -2  0 -1 -4 -3 -3  4  1 -1 -4
 680 Z -1  0  0  1 -3  3  4 -2  0 -3 -3  1 -1 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4
 681 X  0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2  0  0 -2 -1 -1 -1 -1 -1 -4
 682 * -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4  1
 683
 684
 685 __END__