*************************************************************************
* This is a program for protein sequence alignment using the standard
* Needleman-Wunsch dynamic programming. The mutation matrix is from
* BLOSUM62 with gap openning penaly=-11 and gap extension panalty=-1.
* The program can be freely copied and modified provided the notices
* on the head are retained. Comments and bug report should be addressed
* to Yang Zhang (Email: zhng@umich.edu).
*
* Instructions:
* 1, the program can be compiled by
* >gfortran -static -O3 -ffast-math -lm -o align align.f
* 2, simply running the program will give a brief note on how to use it
* 3, You can run the program in following convenient ways:
* >align F1.fasta F2.fasta (align two sequences in fasta file)
* >align F1.pdb F2.pdb 1 (align two sequences in PDB file)
* >align F1.fasta F2.pdb 2 (align Sequence 1 in fasta and 2 in pdb)
* >align GKDGL EVADELVSE 3 (align sequences typed by keyboard)
* >align GKDGL F.fasta 4 (align Seq-1 by keyboard and 2 in fasta)
* >align GKDGL F.pdb 5 (align Seq-1 by keyboard and 2 in pdb)
*
* Reference to cite:
* Y. Zhang, http://zhanglab.dcmb.med.umich.edu/NW-align
*
* Permission to use, copy, modify, and distribute this program for
* any purpose, with or without fee, is hereby granted, provided that
* the notices on the head, the reference information, and this
* copyright notice appear in all copies or substantial portions of
* the Software. It is provided "as is" without express or implied
* warranty.
************************ updating history ********************************
* 2006/08/03: first release of NW-align program
* 2016/05/11: A change is made to allow MSE in sequence
*************************************************************************
program compares
PARAMETER(ndim=6000)
parameter(naa=24) !number of amino acid
common/dpc/score(ndim,ndim),gap_open,gap_extn,j2i(ndim)
& ,nseq1,nseq2
common/matra/imut(naa,naa) !b,z,x are additional
integer seq1(ndim),seq2(ndim)
character*10000 fnam1,fnam2,fnam3,fnam4
character*10000 s
character*3 aa(naa),aanam
character seqw(naa),upper
character*100 du,ad
character sequenceA(ndim),sequenceB(ndim),sequenceM(ndim)
*---------------------- 23 amino acids ---------------------
data aa/'ALA','ARG','ASN','ASP','CYS','GLN','GLU',
& 'GLY','HIS','ILE','LEU','LYS','MET','PHE','PRO','SER',
& 'THR','TRP','TYR','VAL','ASX','GLX','UNK','MSE'/
data seqw/'A','R','N','D','C','Q','E','G','H','I','L','K',
& 'M','F','P','S','T','W','Y','V','B','Z','X','M'/
call getarg(1,fnam1)
call getarg(2,fnam2)
call getarg(3,fnam3)
call getarg(4,fnam4)
if(fnam1.eq.' ')then
write(*,*)'align F1.fasta F2.fasta ',
& '(align two sequences in fasta file)'
write(*,*)'align F1.pdb F2.pdb 1 ',
& '(align two sequences in PDB file)'
write(*,*)'align F1.fasta F2.pdb 2 ',
& '(align Sequence 1 in fasta and 2 in pdb)'
write(*,*)'align GKDGL EVADELVSE 3 ',
& '(align two sequences typed by keyboard)'
write(*,*)'align GKDGL F.fasta 4 ',
& '(align Sequence 1 by keyboard and 2 in fasta)'
write(*,*)'align GKDGL F.pdb 5 ',
& '(align Sequence 1 by keyboard and 2 in pdb)'
goto 999
endif
*1******* read sequences ------------------------->
if(fnam3.eq.'5')then !direct, 555555555555555555
*** read sequence1:
i=0
do k=1,10000
fnam1(k:k)=upper(fnam1(k:k))
do j=1,naa
if(fnam1(k:k).eq.seqw(j))then
i=i+1
seq1(i)=j
goto 5
endif
enddo
if(fnam1(k:k).ne.'-')goto 55 !same time
5 continue
if(i.ge.ndim)goto 55
enddo
55 continue
nseq1=i
*** read sequence2:
open(unit=10,file=fnam2,status='old')
i=0
do while (.true.)
read(10,1,end=551) s
if(i.gt.0.and.s(1:3).eq.'TER')goto 551
if(s(1:3).eq.'ATO')then
if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.
& or.s(13:16).eq.' CA')then
i=i+1
read(s,111)du,aanam
do j=1,naa
if(aanam.eq.aa(j))seq2(i)=j
enddo
endif
endif
if(i.ge.ndim)goto 551
enddo
551 continue
close(10)
nseq2=i
elseif(fnam3.eq.'4')then !direct, 444444444444444444444444444
*** read sequence1:
i=0
do k=1,10000
fnam1(k:k)=upper(fnam1(k:k))
do j=1,naa
if(fnam1(k:k).eq.seqw(j))then
i=i+1
seq1(i)=j
goto 4
endif
enddo
if(fnam1(k:k).ne.'-')goto 44
4 continue
if(i.ge.ndim)goto 44
enddo
44 continue
nseq1=i
*** read sequence2:
open(unit=10,file=fnam2,status='old')
i=0
do while(.true.)
read(10,1,end=443)s
if(s(1:1).eq.'>')goto 442
do k=1,10000
s(k:k)=upper(s(k:k))
do j=1,naa
if(s(k:k).eq.seqw(j))then
i=i+1
seq2(i)=j
goto 441
endif
enddo
if(s(k:k).ne.'-')goto 442 !same time
441 continue
enddo
442 continue
if(i.ge.ndim)goto 443
enddo
443 continue
close(10)
nseq2=i
elseif(fnam3.eq.'3')then !direct, 33333333333333333333333333333333333
*** read sequence1:
i=0
do k=1,10000
fnam1(k:k)=upper(fnam1(k:k))
do j=1,naa
if(fnam1(k:k).eq.seqw(j))then
i=i+1
seq1(i)=j
goto 3
endif
enddo
if(fnam1(k:k).ne.'-')goto 33
3 continue
if(i.ge.ndim)goto 33
enddo
33 continue
nseq1=i
*** read sequence2:
i=0
do k=1,10000
fnam2(k:k)=upper(fnam2(k:k))
do j=1,naa
if(fnam2(k:k).eq.seqw(j))then
i=i+1
seq2(i)=j
goto 331
endif
enddo
if(fnam2(k:k).ne.'-')goto 332
331 continue
if(i.ge.ndim)goto 332
enddo
332 continue
nseq2=i
elseif(fnam3.eq.'1')then !pdb,pdb, 11111111111111111111111111111
*** read sequence1:
open(unit=10,file=fnam1,status='old')
i=0
do while (.true.)
read(10,1,end=11) s
if(i.gt.0.and.s(1:3).eq.'TER')goto 11
if(s(1:3).eq.'ATO')then
if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.
& or.s(13:16).eq.' CA')then
i=i+1
read(s,111)du,aanam
do j=1,naa
if(aanam.eq.aa(j))seq1(i)=j
enddo
endif
endif
if(i.ge.ndim)goto 11
enddo
1 format(A10000)
11 continue
111 format(A17,A3)
close(10)
nseq1=i
*** read sequence2:
open(unit=10,file=fnam2,status='old')
i=0
do while (.true.)
read(10,1,end=112) s
if(i.gt.0.and.s(1:3).eq.'TER')goto 112
if(s(1:3).eq.'ATO')then
if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.
& or.s(13:16).eq.' CA')then
i=i+1
read(s,111)du,aanam
do j=1,naa
if(aanam.eq.aa(j))seq2(i)=j
enddo
endif
endif
if(i.ge.ndim)goto 112
enddo
112 continue
close(10)
nseq2=i
elseif(fnam3.eq.'2')then !seq,pdb 2222222222222222222222222222222
*** read sequence1:
open(unit=10,file=fnam1,status='old')
i=0
do while(.true.)
read(10,1,end=221)s
if(s(1:1).eq.'>')goto 22
do k=1,10000
s(k:k)=upper(s(k:k))
do j=1,naa
if(s(k:k).eq.seqw(j))then
i=i+1
seq1(i)=j
goto 2
endif
enddo
if(s(k:k).ne.'-')goto 22
2 continue
enddo
22 continue
if(i.ge.ndim)goto 221
enddo
221 continue
close(10)
nseq1=i
*** read sequence2:
open(unit=10,file=fnam2,status='old')
i=0
do while (.true.)
read(10,1,end=222) s
if(i.gt.0.and.s(1:3).eq.'TER')goto 222
if(s(1:3).eq.'ATO')then
if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.
& or.s(13:16).eq.' CA')then
i=i+1
read(s,111)du,aanam
do j=1,naa
if(aanam.eq.aa(j))seq2(i)=j
enddo
endif
endif
if(i.ge.ndim)goto 222
enddo
222 continue
close(10)
nseq2=i
else !seq,seq 00000000000000000000000000000000
*** read sequence1:
open(unit=10,file=fnam1,status='old')
i=0
do while(.true.)
read(10,1,end=881)s
if(s(1:1).eq.'>')goto 88
do k=1,10000
s(k:k)=upper(s(k:k))
do j=1,naa
if(s(k:k).eq.seqw(j))then
i=i+1
seq1(i)=j
goto 8
endif
enddo
if(s(k:k).ne.'-')goto 88
8 continue
enddo
88 continue
if(i.ge.ndim)goto 881
enddo
881 continue
close(10)
nseq1=i
*** read sequence2:
open(unit=10,file=fnam2,status='old')
i=0
do while(.true.)
read(10,1,end=884)s
if(s(1:1).eq.'>')goto 883
do k=1,10000
s(k:k)=upper(s(k:k))
do j=1,naa
if(s(k:k).eq.seqw(j))then
i=i+1
seq2(i)=j
goto 882
endif
enddo
if(s(k:k).ne.'-')goto 883
882 continue
enddo
883 continue
if(i.ge.ndim)goto 884
enddo
884 continue
close(10)
nseq2=i
endif
*2** read mutation matrix ---------->
call matrix !take blosum
*** set unit mutation matrix ---------->
c do i=1,naa
c do j=1,naa
c imut(i,j)=0
c enddo
c enddo
c do i=1,naa
c imut(i,i)=1
c enddo
*3** score------------------>
do i=1,nseq1
do j=1,nseq2
score(i,j)=imut(seq1(i),seq2(j))
enddo
enddo
*4*****************************************************************
* dynamatic program:
******************************************************************
gap_open=-11
gap_extn=-1
call DP(score0) !W(k)=Go+Ge*k1+Go+Ge*k2, standard NW
c call DPalt(score0) !W(k)=Go+Ge*k1+Ge*k2, alternative NW
*5** calculate sequence identity---------------------------->
L_id=0
L_ali=0
do j=1,nseq2
if(j2i(j).gt.0)then
i=j2i(j)
L_ali=L_ali+1
if(seq1(i).eq.seq2(j))L_id=L_id+1
endif
enddo
write(*,*)
write(*,101)nseq1,fnam1
101 format('Length of sequence 1: ',I4,' ->',A10)
write(*,102)nseq2,fnam2
102 format('Length of sequence 2: ',I4,' ->',A10)
write(*,103)L_ali
103 format('Aligned length: ',I4)
write(*,104)L_id
104 format('Identical length: ',I4)
write(*,105)float(L_id)/(nseq2+0.00000001),L_id,nseq2
105 format('Sequence identity: ',F8.3,' (=',I4,'/',I4,')')
write(*,*)
*6******************************************************************
*** output aligned sequences
k=0 !final aligned order
i=1 !on sequence 1
j=1 !on sequence 2
800 continue
if(i.gt.nseq1.and.j.gt.nseq2)goto 802
if(i.gt.nseq1.and.j.le.nseq2)then !unaligned C on 1
k=k+1
sequenceA(k)='-'
sequenceB(k)=seqw(seq2(j))
sequenceM(k)=' '
j=j+1
goto 800
endif
if(i.le.nseq1.and.j.gt.nseq2)then !unaligned C on 2
k=k+1
sequenceA(k)=seqw(seq1(i))
sequenceB(k)='-'
sequenceM(k)=' '
i=i+1
goto 800
endif
if(i.eq.j2i(j))then !if aligned
k=k+1
sequenceA(k)=seqw(seq1(i))
sequenceB(k)=seqw(seq2(j))
if(seq1(i).eq.seq2(j))then !identical
sequenceM(k)=':'
else
sequenceM(k)=' '
endif
i=i+1
j=j+1
goto 800
elseif(j2i(j).lt.0)then !if gap on 1
k=k+1
sequenceA(k)='-'
sequenceB(k)=seqw(seq2(j))
sequenceM(k)=' '
j=j+1
goto 800
elseif(j2i(j).gt.0)then !if gap on 2
k=k+1
sequenceA(k)=seqW(seq1(i))
sequenceB(k)='-'
sequenceM(k)=' '
i=i+1
goto 800
endif
802 continue
write(*,601)(sequenceA(i),i=1,k)
write(*,601)(sequenceM(i),i=1,k)
write(*,601)(sequenceB(i),i=1,k)
write(*,602)(mod(i,10),i=1,k)
601 format(20000A1)
602 format(20000I1)
write(*,*)
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c STOP
999 END
********************************************************************
* This is a standard Needleman-Wunsch dynamic program (by Y. Zhang 2005)
* 1. Count multiple-gap.
* 2. The gap penality W(k)=Go+Ge*k1+Go+Ge*k2 if gap open on both sequences
*
* Input: score(i,j), gap_open, gap_extn
* Output: j2i(j)
* idir(i,j)=1,2,3, from diagonal, horizontal, vertical
* val(i,j) is the cumulative score of (i,j)
********************************************************************
subroutine DP(score0)
PARAMETER(ndim=6000)
common/dpc/score(ndim,ndim),gap_open,gap_extn,j2i(ndim)
& ,nseq1,nseq2
dimension val(0:ndim,0:ndim),idir(0:ndim,0:ndim)
dimension jpV(0:ndim,0:ndim),jpH(0:ndim,0:ndim)
dimension preV(0:ndim,0:ndim),preH(0:ndim,0:ndim)
real D,V,H
ccc initializations --------------->
val(0,0)=0.0
val(1,0)=gap_open
do i=2,nseq1
val(i,0)=val(i-1,0)+gap_extn
enddo
do i=1,nseq1
preV(i,0)=val(i,0)
idir(i,0)=0
jpV(i,0)=1
jpH(i,0)=i
enddo
val(0,1)=gap_open
do j=2,nseq2
val(0,j)=val(0,j-1)+gap_extn
enddo
do j=1,nseq2
preH(0,j)=val(0,j)
idir(0,j)=0
jpV(0,j)=j
jpH(0,j)=1
enddo
ccc DP ------------------------------>
do 111 j=1,nseq2
do 222 i=1,nseq1
ccc D=VAL(i-1,j-1)+SCORE(i,j)--------------->
D=val(i-1,j-1)+score(i,j) !from diagonal, val(i,j) is val(i-1,j-1)
ccc H=H+gap_open ------->
jpH(i,j)=1
val1=val(i-1,j)+gap_open !gap_open from both D and V
val2=preH(i-1,j)+gap_extn !gap_extn from horizontal
if(val1.gt.val2) then !last step from D or V
H=val1
else !last step from H
H=val2
if(i.gt.1)jpH(i,j)=jpH(i-1,j)+1 !record long-gap
endif
ccc V=V+gap_open --------->
jpV(i,j)=1
val1=val(i,j-1)+gap_open
val2=preV(i,j-1)+gap_extn
if(val1.gt.val2) then
V=val1
else
V=val2
if(j.gt.1)jpV(i,j)=jpV(i,j-1)+1
endif
preH(i,j)=H !unaccepted H
preV(i,j)=V !unaccepted V
if(D.gt.H.and.D.gt.V)then
idir(i,j)=1
val(i,j)=D
elseif(H.gt.V)then
idir(i,j)=2
val(i,j)=H
else
idir(i,j)=3
val(i,j)=V
endif
222 continue
111 continue
score0=val(nseq1,nseq2) !alignment score
c tracing back the pathway:
do j=1,nseq2
j2i(j)=-1 !all are not aligned
enddo
i=nseq1
j=nseq2
do while(i.gt.0.and.j.gt.0)
if(idir(i,j).eq.1)then !from diagonal
j2i(j)=i
i=i-1
j=j-1
elseif(idir(i,j).eq.2)then !from horizonal
it=jpH(i,j)
do me=1,it
if(i.gt.0) then
i=i-1
endif
enddo
else
it=jpV(i,j)
do me=1,it
if(j.gt.0) then
j=j-1
endif
enddo
endif
enddo
*^^^^^^^^^^^DP finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
return
end
********************************************************************
* This is an alternative implementation of Needleman-Wunsch dynamic program
* (by Y. Zhang 2005)
* 1. Count two-layer iteration and multiple-gaps
* 2. The gap penality W(k)=Go+Ge*k1+Ge*k2 if gap open on both sequences
*
* Input: score(i,j), gap_open, gap_extn
* Output: j2i(j)
* idir(i,j)=1,2,3, from diagonal, horizontal, vertical
* val(i,j) is the cumulative score of (i,j)
********************************************************************
subroutine DPalt(score0)
PARAMETER(ndim=6000)
common/dpc/score(ndim,ndim),gap_open,gap_extn,j2i(ndim)
& ,nseq1,nseq2
dimension val(0:ndim,0:ndim),idir(0:ndim,0:ndim)
dimension preV(0:ndim,0:ndim),preH(0:ndim,0:ndim),
& preD(0:ndim,0:ndim)
dimension idirH(0:ndim,0:ndim),idirV(0:ndim,0:ndim)
ccc initializations --------------->
val(0,0)=0.0
val(1,0)=gap_open
do i=2,nseq1
val(i,0)=val(i-1,0)+gap_extn
enddo
do i=1,nseq1
idir(i,0)=0
preD(i,0)=val(i,0)
preH(i,0)=val(i,0)
preV(i,0)=val(i,0)
enddo
val(0,1)=gap_open
do j=2,nseq2
val(0,j)=val(0,j-1)+gap_extn
enddo
do j=1,nseq2
idir(0,j)=0
preD(0,j)=val(0,j)
preH(0,j)=val(0,j)
preV(0,j)=val(0,j)
enddo
ccc DP ------------------------------>
do 111 j=1,nseq2
do 222 i=1,nseq1
ccc preD=VAL(i-1,j-1)+SCORE(i,j)--------------->
preD(i,j)=val(i-1,j-1)+score(i,j)
ccc preH: pre-accepted H----------------------->
D=preD(i-1,j)+gap_open
H=preH(i-1,j)+gap_extn
V=preV(i-1,j)+gap_extn
if(D.gt.H.and.D.gt.V)then
preH(i,j)=D
idirH(i-1,j)=1
elseif(H.gt.V)then
preH(i,j)=H
idirH(i-1,j)=2
else
preH(i,j)=V
idirH(i-1,j)=3
endif
ccc preV: pre-accepted V----------------------->
D=preD(i,j-1)+gap_open
H=preH(i,j-1)+gap_extn
V=preV(i,j-1)+gap_extn
if(D.gt.H.and.D.gt.V)then
preV(i,j)=D
idirV(i,j-1)=1
elseif(H.gt.V)then
preV(i,j)=H
idirV(i,j-1)=2
else
preV(i,j)=V
idirV(i,j-1)=3
endif
ccc decide idir(i,j)----------->
if(preD(i,j).gt.preH(i,j).and.preD(i,j).gt.preV(i,j))then
idir(i,j)=1
val(i,j)=preD(i,j)
elseif(preH(i,j).gt.preV(i,j))then
idir(i,j)=2
val(i,j)=preH(i,j)
else
idir(i,j)=3
val(i,j)=preV(i,j)
endif
222 continue
111 continue
score0=val(nseq1,nseq2) !alignment score
ccc tracing back the pathway:
do j=1,nseq2
j2i(j)=-1 !all are not aligned
enddo
i=nseq1
j=nseq2
do while(i.gt.0.and.j.gt.0)
if(idir(i,j).eq.1)then !from diagonal
j2i(j)=i
i=i-1
j=j-1
elseif(idir(i,j).eq.2)then
i=i-1
idir(i,j)=idirH(i,j)
else
j=j-1
idir(i,j)=idirV(i,j)
endif
enddo
*^^^^^^^^^^^DP finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
return
end
********************************************************************
* read matrix
*
subroutine matrix
parameter(naa=24) !number of amino acid
common/matra/imut(naa,naa) !b,z,x are additional
* folowing from BLOSUM62 used in BLAST:
imut(1,1)=4
imut(1,2)=-1
imut(1,3)=-2
imut(1,4)=-2
imut(1,5)=0
imut(1,6)=-1
imut(1,7)=-1
imut(1,8)=0
imut(1,9)=-2
imut(1,10)=-1
imut(1,11)=-1
imut(1,12)=-1
imut(1,13)=-1
imut(1,14)=-2
imut(1,15)=-1
imut(1,16)=1
imut(1,17)=0
imut(1,18)=-3
imut(1,19)=-2
imut(1,20)=0
imut(1,21)=-2
imut(1,22)=-1
imut(1,23)=0
imut(2,1)=-1
imut(2,2)=5
imut(2,3)=0
imut(2,4)=-2
imut(2,5)=-3
imut(2,6)=1
imut(2,7)=0
imut(2,8)=-2
imut(2,9)=0
imut(2,10)=-3
imut(2,11)=-2
imut(2,12)=2
imut(2,13)=-1
imut(2,14)=-3
imut(2,15)=-2
imut(2,16)=-1
imut(2,17)=-1
imut(2,18)=-3
imut(2,19)=-2
imut(2,20)=-3
imut(2,21)=-1
imut(2,22)=0
imut(2,23)=-1
imut(3,1)=-2
imut(3,2)=0
imut(3,3)=6
imut(3,4)=1
imut(3,5)=-3
imut(3,6)=0
imut(3,7)=0
imut(3,8)=0
imut(3,9)=1
imut(3,10)=-3
imut(3,11)=-3
imut(3,12)=0
imut(3,13)=-2
imut(3,14)=-3
imut(3,15)=-2
imut(3,16)=1
imut(3,17)=0
imut(3,18)=-4
imut(3,19)=-2
imut(3,20)=-3
imut(3,21)=3
imut(3,22)=0
imut(3,23)=-1
imut(4,1)=-2
imut(4,2)=-2
imut(4,3)=1
imut(4,4)=6
imut(4,5)=-3
imut(4,6)=0
imut(4,7)=2
imut(4,8)=-1
imut(4,9)=-1
imut(4,10)=-3
imut(4,11)=-4
imut(4,12)=-1
imut(4,13)=-3
imut(4,14)=-3
imut(4,15)=-1
imut(4,16)=0
imut(4,17)=-1
imut(4,18)=-4
imut(4,19)=-3
imut(4,20)=-3
imut(4,21)=4
imut(4,22)=1
imut(4,23)=-1
imut(5,1)=0
imut(5,2)=-3
imut(5,3)=-3
imut(5,4)=-3
imut(5,5)=9
imut(5,6)=-3
imut(5,7)=-4
imut(5,8)=-3
imut(5,9)=-3
imut(5,10)=-1
imut(5,11)=-1
imut(5,12)=-3
imut(5,13)=-1
imut(5,14)=-2
imut(5,15)=-3
imut(5,16)=-1
imut(5,17)=-1
imut(5,18)=-2
imut(5,19)=-2
imut(5,20)=-1
imut(5,21)=-3
imut(5,22)=-3
imut(5,23)=-2
imut(6,1)=-1
imut(6,2)=1
imut(6,3)=0
imut(6,4)=0
imut(6,5)=-3
imut(6,6)=5
imut(6,7)=2
imut(6,8)=-2
imut(6,9)=0
imut(6,10)=-3
imut(6,11)=-2
imut(6,12)=1
imut(6,13)=0
imut(6,14)=-3
imut(6,15)=-1
imut(6,16)=0
imut(6,17)=-1
imut(6,18)=-2
imut(6,19)=-1
imut(6,20)=-2
imut(6,21)=0
imut(6,22)=3
imut(6,23)=-1
imut(7,1)=-1
imut(7,2)=0
imut(7,3)=0
imut(7,4)=2
imut(7,5)=-4
imut(7,6)=2
imut(7,7)=5
imut(7,8)=-2
imut(7,9)=0
imut(7,10)=-3
imut(7,11)=-3
imut(7,12)=1
imut(7,13)=-2
imut(7,14)=-3
imut(7,15)=-1
imut(7,16)=0
imut(7,17)=-1
imut(7,18)=-3
imut(7,19)=-2
imut(7,20)=-2
imut(7,21)=1
imut(7,22)=4
imut(7,23)=-1
imut(8,1)=0
imut(8,2)=-2
imut(8,3)=0
imut(8,4)=-1
imut(8,5)=-3
imut(8,6)=-2
imut(8,7)=-2
imut(8,8)=6
imut(8,9)=-2
imut(8,10)=-4
imut(8,11)=-4
imut(8,12)=-2
imut(8,13)=-3
imut(8,14)=-3
imut(8,15)=-2
imut(8,16)=0
imut(8,17)=-2
imut(8,18)=-2
imut(8,19)=-3
imut(8,20)=-3
imut(8,21)=-1
imut(8,22)=-2
imut(8,23)=-1
imut(9,1)=-2
imut(9,2)=0
imut(9,3)=1
imut(9,4)=-1
imut(9,5)=-3
imut(9,6)=0
imut(9,7)=0
imut(9,8)=-2
imut(9,9)=8
imut(9,10)=-3
imut(9,11)=-3
imut(9,12)=-1
imut(9,13)=-2
imut(9,14)=-1
imut(9,15)=-2
imut(9,16)=-1
imut(9,17)=-2
imut(9,18)=-2
imut(9,19)=2
imut(9,20)=-3
imut(9,21)=0
imut(9,22)=0
imut(9,23)=-1
imut(10,1)=-1
imut(10,2)=-3
imut(10,3)=-3
imut(10,4)=-3
imut(10,5)=-1
imut(10,6)=-3
imut(10,7)=-3
imut(10,8)=-4
imut(10,9)=-3
imut(10,10)=4
imut(10,11)=2
imut(10,12)=-3
imut(10,13)=1
imut(10,14)=0
imut(10,15)=-3
imut(10,16)=-2
imut(10,17)=-1
imut(10,18)=-3
imut(10,19)=-1
imut(10,20)=3
imut(10,21)=-3
imut(10,22)=-3
imut(10,23)=-1
imut(11,1)=-1
imut(11,2)=-2
imut(11,3)=-3
imut(11,4)=-4
imut(11,5)=-1
imut(11,6)=-2
imut(11,7)=-3
imut(11,8)=-4
imut(11,9)=-3
imut(11,10)=2
imut(11,11)=4
imut(11,12)=-2
imut(11,13)=2
imut(11,14)=0
imut(11,15)=-3
imut(11,16)=-2
imut(11,17)=-1
imut(11,18)=-2
imut(11,19)=-1
imut(11,20)=1
imut(11,21)=-4
imut(11,22)=-3
imut(11,23)=-1
imut(12,1)=-1
imut(12,2)=2
imut(12,3)=0
imut(12,4)=-1
imut(12,5)=-3
imut(12,6)=1
imut(12,7)=1
imut(12,8)=-2
imut(12,9)=-1
imut(12,10)=-3
imut(12,11)=-2
imut(12,12)=5
imut(12,13)=-1
imut(12,14)=-3
imut(12,15)=-1
imut(12,16)=0
imut(12,17)=-1
imut(12,18)=-3
imut(12,19)=-2
imut(12,20)=-2
imut(12,21)=0
imut(12,22)=1
imut(12,23)=-1
imut(13,1)=-1
imut(13,2)=-1
imut(13,3)=-2
imut(13,4)=-3
imut(13,5)=-1
imut(13,6)=0
imut(13,7)=-2
imut(13,8)=-3
imut(13,9)=-2
imut(13,10)=1
imut(13,11)=2
imut(13,12)=-1
imut(13,13)=5
imut(13,14)=0
imut(13,15)=-2
imut(13,16)=-1
imut(13,17)=-1
imut(13,18)=-1
imut(13,19)=-1
imut(13,20)=1
imut(13,21)=-3
imut(13,22)=-1
imut(13,23)=-1
imut(14,1)=-2
imut(14,2)=-3
imut(14,3)=-3
imut(14,4)=-3
imut(14,5)=-2
imut(14,6)=-3
imut(14,7)=-3
imut(14,8)=-3
imut(14,9)=-1
imut(14,10)=0
imut(14,11)=0
imut(14,12)=-3
imut(14,13)=0
imut(14,14)=6
imut(14,15)=-4
imut(14,16)=-2
imut(14,17)=-2
imut(14,18)=1
imut(14,19)=3
imut(14,20)=-1
imut(14,21)=-3
imut(14,22)=-3
imut(14,23)=-1
imut(15,1)=-1
imut(15,2)=-2
imut(15,3)=-2
imut(15,4)=-1
imut(15,5)=-3
imut(15,6)=-1
imut(15,7)=-1
imut(15,8)=-2
imut(15,9)=-2
imut(15,10)=-3
imut(15,11)=-3
imut(15,12)=-1
imut(15,13)=-2
imut(15,14)=-4
imut(15,15)=7
imut(15,16)=-1
imut(15,17)=-1
imut(15,18)=-4
imut(15,19)=-3
imut(15,20)=-2
imut(15,21)=-2
imut(15,22)=-1
imut(15,23)=-2
imut(16,1)=1
imut(16,2)=-1
imut(16,3)=1
imut(16,4)=0
imut(16,5)=-1
imut(16,6)=0
imut(16,7)=0
imut(16,8)=0
imut(16,9)=-1
imut(16,10)=-2
imut(16,11)=-2
imut(16,12)=0
imut(16,13)=-1
imut(16,14)=-2
imut(16,15)=-1
imut(16,16)=4
imut(16,17)=1
imut(16,18)=-3
imut(16,19)=-2
imut(16,20)=-2
imut(16,21)=0
imut(16,22)=0
imut(16,23)=0
imut(17,1)=0
imut(17,2)=-1
imut(17,3)=0
imut(17,4)=-1
imut(17,5)=-1
imut(17,6)=-1
imut(17,7)=-1
imut(17,8)=-2
imut(17,9)=-2
imut(17,10)=-1
imut(17,11)=-1
imut(17,12)=-1
imut(17,13)=-1
imut(17,14)=-2
imut(17,15)=-1
imut(17,16)=1
imut(17,17)=5
imut(17,18)=-2
imut(17,19)=-2
imut(17,20)=0
imut(17,21)=-1
imut(17,22)=-1
imut(17,23)=0
imut(18,1)=-3
imut(18,2)=-3
imut(18,3)=-4
imut(18,4)=-4
imut(18,5)=-2
imut(18,6)=-2
imut(18,7)=-3
imut(18,8)=-2
imut(18,9)=-2
imut(18,10)=-3
imut(18,11)=-2
imut(18,12)=-3
imut(18,13)=-1
imut(18,14)=1
imut(18,15)=-4
imut(18,16)=-3
imut(18,17)=-2
imut(18,18)=11
imut(18,19)=2
imut(18,20)=-3
imut(18,21)=-4
imut(18,22)=-3
imut(18,23)=-2
imut(19,1)=-2
imut(19,2)=-2
imut(19,3)=-2
imut(19,4)=-3
imut(19,5)=-2
imut(19,6)=-1
imut(19,7)=-2
imut(19,8)=-3
imut(19,9)=2
imut(19,10)=-1
imut(19,11)=-1
imut(19,12)=-2
imut(19,13)=-1
imut(19,14)=3
imut(19,15)=-3
imut(19,16)=-2
imut(19,17)=-2
imut(19,18)=2
imut(19,19)=7
imut(19,20)=-1
imut(19,21)=-3
imut(19,22)=-2
imut(19,23)=-1
imut(20,1)=0
imut(20,2)=-3
imut(20,3)=-3
imut(20,4)=-3
imut(20,5)=-1
imut(20,6)=-2
imut(20,7)=-2
imut(20,8)=-3
imut(20,9)=-3
imut(20,10)=3
imut(20,11)=1
imut(20,12)=-2
imut(20,13)=1
imut(20,14)=-1
imut(20,15)=-2
imut(20,16)=-2
imut(20,17)=0
imut(20,18)=-3
imut(20,19)=-1
imut(20,20)=4
imut(20,21)=-3
imut(20,22)=-2
imut(20,23)=-1
imut(21,1)=-2
imut(21,2)=-1
imut(21,3)=3
imut(21,4)=4
imut(21,5)=-3
imut(21,6)=0
imut(21,7)=1
imut(21,8)=-1
imut(21,9)=0
imut(21,10)=-3
imut(21,11)=-4
imut(21,12)=0
imut(21,13)=-3
imut(21,14)=-3
imut(21,15)=-2
imut(21,16)=0
imut(21,17)=-1
imut(21,18)=-4
imut(21,19)=-3
imut(21,20)=-3
imut(21,21)=4
imut(21,22)=1
imut(21,23)=-1
imut(22,1)=-1
imut(22,2)=0
imut(22,3)=0
imut(22,4)=1
imut(22,5)=-3
imut(22,6)=3
imut(22,7)=4
imut(22,8)=-2
imut(22,9)=0
imut(22,10)=-3
imut(22,11)=-3
imut(22,12)=1
imut(22,13)=-1
imut(22,14)=-3
imut(22,15)=-1
imut(22,16)=0
imut(22,17)=-1
imut(22,18)=-3
imut(22,19)=-2
imut(22,20)=-2
imut(22,21)=1
imut(22,22)=4
imut(22,23)=-1
imut(23,1)=0
imut(23,2)=-1
imut(23,3)=-1
imut(23,4)=-1
imut(23,5)=-2
imut(23,6)=-1
imut(23,7)=-1
imut(23,8)=-1
imut(23,9)=-1
imut(23,10)=-1
imut(23,11)=-1
imut(23,12)=-1
imut(23,13)=-1
imut(23,14)=-1
imut(23,15)=-2
imut(23,16)=0
imut(23,17)=0
imut(23,18)=-2
imut(23,19)=-1
imut(23,20)=-1
imut(23,21)=-1
imut(23,22)=-1
imut(23,23)=-1
return
end
function upper(A)
CHARACTER A,upper
IF(A.LE.'z'.and.A.GE.'a')then
A=CHAR(ICHAR(A)-32)
endif
upper=A
RETURN
END