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