Mercurial > repos > cbib > fibronectin

comparison fibronectin/NWalign_PAM30.f @ 0:0c6cfb9906f3 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author cbib
date Wed, 10 Nov 2021 15:15:50 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:0c6cfb9906f3
1 *************************************************************************
2 * This is a program for protein sequence alignment using the standard
3 * Needleman-Wunsch dynamic programming. The mutation matrix is from
4 * PAM30 with gap openning penaly=-11 and gap extension panalty=-1.
5 * The program can be freely copied and modified provided the notices
6 * on the head are retained. Comments and bug report should be addressed
7 * to Yang Zhang (Email: zhng@umich.edu). Last update is in 2010/08/03.
8 *
9 * Instructions:
10 * 1, the program can be compiled by
11 * >gfortran -static -O3 -ffast-math -lm -o align align.f
12 * 2, simply running the program will give a brief note on how to use it
13 * 3, You can run the program in following convenient ways:
14 * >align F1.fasta F2.fasta (align two sequences in fasta file)
15 * >align F1.pdb F2.pdb 1 (align two sequences in PDB file)
16 * >align F1.fasta F2.pdb 2 (align Sequence 1 in fasta and 2 in pdb)
17 * >align GKDGL EVADELVSE 3 (align sequences typed by keyboard)
18 * >align GKDGL F.fasta 4 (align Seq-1 by keyboard and 2 in fasta)
19 * >align GKDGL F.pdb 5 (align Seq-1 by keyboard and 2 in pdb)
20 *************************************************************************
21
22 program compares
23 PARAMETER(ndim=6000)
24 parameter(naa=24) !number of amino acid
25 common/dpc/score(ndim,ndim),gap_open,gap_extn,j2i(ndim)
26 & ,nseq1,nseq2
27 common/matra/imut(naa,naa) !b,z,x are additional
28
29 integer seq1(ndim),seq2(ndim)
30 character*10000 fnam1,fnam2,fnam3,fnam4
31 character*10000 s
32 character*3 aa(naa),aanam
33 character seqw(naa),upper
34 character*100 du,ad
35 character sequenceA(ndim),sequenceB(ndim),sequenceM(ndim)
36
37 *---------------------- 24 amino acids ---------------------
38 data aa/'ALA','ARG','ASN','ASP','CYS','GLN','GLU',
39 & 'GLY','HIS','ILE','LEU','LYS','MET','PHE','PRO','SER',
40 & 'THR','TRP','TYR','VAL','ASX','GLX','UNK','STOP'/
41 data seqw/'A','R','N','D','C','Q','E','G','H','I','L','K',
42 & 'M','F','P','S','T','W','Y','V','B','Z','X','*'/
43
44 call getarg(1,fnam1)
45 call getarg(2,fnam2)
46 call getarg(3,fnam3)
47 call getarg(4,fnam4)
48
49 if(fnam1.eq.' ')then
50 write(*,*)'align F1.fasta F2.fasta ',
51 & '(align two sequences in fasta file)'
52 write(*,*)'align F1.pdb F2.pdb 1 ',
53 & '(align two sequences in PDB file)'
54 write(*,*)'align F1.fasta F2.pdb 2 ',
55 & '(align Sequence 1 in fasta and 2 in pdb)'
56 write(*,*)'align GKDGL EVADELVSE 3 ',
57 & '(align two sequences typed by keyboard)'
58 write(*,*)'align GKDGL F.fasta 4 ',
59 & '(align Sequence 1 by keyboard and 2 in fasta)'
60 write(*,*)'align GKDGL F.pdb 5 ',
61 & '(align Sequence 1 by keyboard and 2 in pdb)'
62 goto 999
63 endif
64
65 *1******* read sequences ------------------------->
66 if(fnam3.eq.'5')then !direct, 555555555555555555
67 *** read sequence1:
68 i=0
69 do k=1,10000
70 fnam1(k:k)=upper(fnam1(k:k))
71 do j=1,naa
72 if(fnam1(k:k).eq.seqw(j))then
73 i=i+1
74 seq1(i)=j
75 goto 5
76 endif
77 enddo
78 if(fnam1(k:k).ne.'-')goto 55 !same time
79 5 continue
80 if(i.ge.ndim)goto 55
81 enddo
82 55 continue
83 nseq1=i
84 *** read sequence2:
85 open(unit=10,file=fnam2,status='old')
86 i=0
87 do while (.true.)
88 read(10,1,end=551) s
89 if(i.gt.0.and.s(1:3).eq.'TER')goto 551
90 if(s(1:3).eq.'ATO')then
91 if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.
92 & or.s(13:16).eq.' CA')then
93 i=i+1
94 read(s,111)du,aanam
95 do j=1,naa
96 if(aanam.eq.aa(j))seq2(i)=j
97 enddo
98 endif
99 endif
100 if(i.ge.ndim)goto 551
101 enddo
102 551 continue
103 close(10)
104 nseq2=i
105 elseif(fnam3.eq.'4')then !direct, 444444444444444444444444444
106 *** read sequence1:
107 i=0
108 do k=1,10000
109 fnam1(k:k)=upper(fnam1(k:k))
110 do j=1,naa
111 if(fnam1(k:k).eq.seqw(j))then
112 i=i+1
113 seq1(i)=j
114 goto 4
115 endif
116 enddo
117 if(fnam1(k:k).ne.'-')goto 44
118 4 continue
119 if(i.ge.ndim)goto 44
120 enddo
121 44 continue
122 nseq1=i
123 *** read sequence2:
124 open(unit=10,file=fnam2,status='old')
125 i=0
126 do while(.true.)
127 read(10,1,end=443)s
128 if(s(1:1).eq.'>')goto 442
129 do k=1,10000
130 s(k:k)=upper(s(k:k))
131 do j=1,naa
132 if(s(k:k).eq.seqw(j))then
133 i=i+1
134 seq2(i)=j
135 goto 441
136 endif
137 enddo
138 if(s(k:k).ne.'-')goto 442 !same time
139 441 continue
140 enddo
141 442 continue
142 if(i.ge.ndim)goto 443
143 enddo
144 443 continue
145 close(10)
146 nseq2=i
147
148
149 elseif(fnam3.eq.'3')then !direct, 33333333333333333333333333333333333
150 *** read sequence1:
151 i=0
152 do k=1,10000
153 fnam1(k:k)=upper(fnam1(k:k))
154 do j=1,naa
155 if(fnam1(k:k).eq.seqw(j))then
156 i=i+1
157 seq1(i)=j
158 goto 3
159 endif
160 enddo
161 if(fnam1(k:k).ne.'-')goto 33
162 3 continue
163 if(i.ge.ndim)goto 33
164 enddo
165 33 continue
166 nseq1=i
167 *** read sequence2:
168 i=0
169 do k=1,10000
170 fnam2(k:k)=upper(fnam2(k:k))
171 do j=1,naa
172 if(fnam2(k:k).eq.seqw(j))then
173 i=i+1
174 seq2(i)=j
175 goto 331
176 endif
177 enddo
178 if(fnam2(k:k).ne.'-')goto 332
179 331 continue
180 if(i.ge.ndim)goto 332
181 enddo
182 332 continue
183 nseq2=i
184 elseif(fnam3.eq.'1')then !pdb,pdb, 11111111111111111111111111111
185 *** read sequence1:
186 open(unit=10,file=fnam1,status='old')
187 i=0
188 do while (.true.)
189 read(10,1,end=11) s
190 if(i.gt.0.and.s(1:3).eq.'TER')goto 11
191 if(s(1:3).eq.'ATO')then
192 if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.
193 & or.s(13:16).eq.' CA')then
194 i=i+1
195 read(s,111)du,aanam
196 do j=1,naa
197 if(aanam.eq.aa(j))seq1(i)=j
198 enddo
199 endif
200 endif
201 if(i.ge.ndim)goto 11
202 enddo
203 1 format(A10000)
204 11 continue
205 111 format(A17,A3)
206 close(10)
207 nseq1=i
208 *** read sequence2:
209 open(unit=10,file=fnam2,status='old')
210 i=0
211 do while (.true.)
212 read(10,1,end=112) s
213 if(i.gt.0.and.s(1:3).eq.'TER')goto 112
214 if(s(1:3).eq.'ATO')then
215 if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.
216 & or.s(13:16).eq.' CA')then
217 i=i+1
218 read(s,111)du,aanam
219 do j=1,naa
220 if(aanam.eq.aa(j))seq2(i)=j
221 enddo
222 endif
223 endif
224 if(i.ge.ndim)goto 112
225 enddo
226 112 continue
227 close(10)
228 nseq2=i
229 elseif(fnam3.eq.'2')then !seq,pdb 2222222222222222222222222222222
230 *** read sequence1:
231 open(unit=10,file=fnam1,status='old')
232 i=0
233 do while(.true.)
234 read(10,1,end=221)s
235 if(s(1:1).eq.'>')goto 22
236 do k=1,10000
237 s(k:k)=upper(s(k:k))
238 do j=1,naa
239 if(s(k:k).eq.seqw(j))then
240 i=i+1
241 seq1(i)=j
242 goto 2
243 endif
244 enddo
245 if(s(k:k).ne.'-')goto 22
246 2 continue
247 enddo
248 22 continue
249 if(i.ge.ndim)goto 221
250 enddo
251 221 continue
252 close(10)
253 nseq1=i
254 *** read sequence2:
255 open(unit=10,file=fnam2,status='old')
256 i=0
257 do while (.true.)
258 read(10,1,end=222) s
259 if(i.gt.0.and.s(1:3).eq.'TER')goto 222
260 if(s(1:3).eq.'ATO')then
261 if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.
262 & or.s(13:16).eq.' CA')then
263 i=i+1
264 read(s,111)du,aanam
265 do j=1,naa
266 if(aanam.eq.aa(j))seq2(i)=j
267 enddo
268 endif
269 endif
270 if(i.ge.ndim)goto 222
271 enddo
272 222 continue
273 close(10)
274 nseq2=i
275 else !seq,seq 00000000000000000000000000000000
276 *** read sequence1:
277 open(unit=10,file=fnam1,status='old')
278 i=0
279 do while(.true.)
280 read(10,1,end=881)s
281 if(s(1:1).eq.'>')goto 88
282 do k=1,10000
283 s(k:k)=upper(s(k:k))
284 do j=1,naa
285 if(s(k:k).eq.seqw(j))then
286 i=i+1
287 seq1(i)=j
288 goto 8
289 endif
290 enddo
291 if(s(k:k).ne.'-')goto 88
292 8 continue
293 enddo
294 88 continue
295 if(i.ge.ndim)goto 881
296 enddo
297 881 continue
298 close(10)
299 nseq1=i
300 *** read sequence2:
301 open(unit=10,file=fnam2,status='old')
302 i=0
303 do while(.true.)
304 read(10,1,end=884)s
305 if(s(1:1).eq.'>')goto 883
306 do k=1,10000
307 s(k:k)=upper(s(k:k))
308 do j=1,naa
309 if(s(k:k).eq.seqw(j))then
310 i=i+1
311 seq2(i)=j
312 goto 882
313 endif
314 enddo
315 if(s(k:k).ne.'-')goto 883
316 882 continue
317 enddo
318 883 continue
319 if(i.ge.ndim)goto 884
320 enddo
321 884 continue
322 close(10)
323 nseq2=i
324 endif
325
326 *2** read mutation matrix ---------->
327 call matrix !take pam
328 *** set unit mutation matrix ---------->
329 c do i=1,naa
330 c do j=1,naa
331 c imut(i,j)=0
332 c enddo
333 c enddo
334 c do i=1,naa
335 c imut(i,i)=1
336 c enddo
337
338 *3** score------------------>
339 do i=1,nseq1
340 do j=1,nseq2
341 score(i,j)=imut(seq1(i),seq2(j))
342 enddo
343 enddo
344
345 *4*****************************************************************
346 * dynamatic program:
347 ******************************************************************
348 gap_open=-11
349 gap_extn=-1
350 call DP(score0) !W(k)=Go+Ge*k1+Go+Ge*k2, standard NW
351 c call DPalt(score0) !W(k)=Go+Ge*k1+Ge*k2, alternative NW
352
353 *5** calculate sequence identity---------------------------->
354 L_id=0
355 L_ali=0
356 do j=1,nseq2
357 if(j2i(j).gt.0)then
358 i=j2i(j)
359 L_ali=L_ali+1
360 if(seq1(i).eq.seq2(j))L_id=L_id+1
361 endif
362 enddo
363
364 write(*,*)
365 write(*,101)nseq1,fnam1
366 101 format('Length of sequence 1: ',I4,' ->',A10)
367 write(*,102)nseq2,fnam2
368 102 format('Length of sequence 2: ',I4,' ->',A10)
369 write(*,103)L_ali
370 103 format('Aligned length: ',I4)
371 write(*,104)L_id
372 104 format('Identical length: ',I4)
373 write(*,105)float(L_id)/(nseq2+0.00000001),L_id,nseq2
374 105 format('Sequence identity: ',F8.3,' (=',I4,'/',I4,')')
375 write(*,*)
376
377 *6******************************************************************
378 *** output aligned sequences
379 k=0 !final aligned order
380 i=1 !on sequence 1
381 j=1 !on sequence 2
382 800 continue
383 if(i.gt.nseq1.and.j.gt.nseq2)goto 802
384 if(i.gt.nseq1.and.j.le.nseq2)then !unaligned C on 1
385 k=k+1
386 sequenceA(k)='-'
387 sequenceB(k)=seqw(seq2(j))
388 sequenceM(k)=' '
389 j=j+1
390 goto 800
391 endif
392 if(i.le.nseq1.and.j.gt.nseq2)then !unaligned C on 2
393 k=k+1
394 sequenceA(k)=seqw(seq1(i))
395 sequenceB(k)='-'
396 sequenceM(k)=' '
397 i=i+1
398 goto 800
399 endif
400 if(i.eq.j2i(j))then !if aligned
401 k=k+1
402 sequenceA(k)=seqw(seq1(i))
403 sequenceB(k)=seqw(seq2(j))
404 if(seq1(i).eq.seq2(j))then !identical
405 sequenceM(k)=':'
406 else
407 sequenceM(k)=' '
408 endif
409 i=i+1
410 j=j+1
411 goto 800
412 elseif(j2i(j).lt.0)then !if gap on 1
413 k=k+1
414 sequenceA(k)='-'
415 sequenceB(k)=seqw(seq2(j))
416 sequenceM(k)=' '
417 j=j+1
418 goto 800
419 elseif(j2i(j).gt.0)then !if gap on 2
420 k=k+1
421 sequenceA(k)=seqW(seq1(i))
422 sequenceB(k)='-'
423 sequenceM(k)=' '
424 i=i+1
425 goto 800
426 endif
427 802 continue
428
429 write(*,601)(sequenceA(i),i=1,k)
430 write(*,601)(sequenceM(i),i=1,k)
431 write(*,601)(sequenceB(i),i=1,k)
432 write(*,602)(mod(i,10),i=1,k)
433 601 format(2000A1)
434 602 format(2000I1)
435 write(*,*)
436
437 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
438 c STOP
439 999 END
440
441 ********************************************************************
442 * This is a standard Needleman-Wunsch dynamic program (by Y. Zhang 2005)
443 * 1. Count multiple-gap.
444 * 2. The gap penality W(k)=Go+Ge*k1+Go+Ge*k2 if gap open on both sequences
445 *
446 * Input: score(i,j), gap_open, gap_extn
447 * Output: j2i(j)
448 * idir(i,j)=1,2,3, from diagonal, horizontal, vertical
449 * val(i,j) is the cumulative score of (i,j)
450 ********************************************************************
451 subroutine DP(score0)
452 PARAMETER(ndim=6000)
453 common/dpc/score(ndim,ndim),gap_open,gap_extn,j2i(ndim)
454 & ,nseq1,nseq2
455
456 dimension val(0:ndim,0:ndim),idir(0:ndim,0:ndim)
457 dimension jpV(0:ndim,0:ndim),jpH(0:ndim,0:ndim)
458 dimension preV(0:ndim,0:ndim),preH(0:ndim,0:ndim)
459 real D,V,H
460
461 ccc initializations --------------->
462 val(0,0)=0.0
463 do i=1,nseq1
464 val(i,0)=gap_extn*i
465 preV(i,0)=val(i,0) !not use preV at the beginning
466 idir(i,0)=0 !useless
467 jpV(i,0)=1 !useless
468 jpH(i,0)=i !useless
469 enddo
470 do j=1,nseq2
471 val(0,j)=gap_extn*j
472 preH(0,j)=val(0,j)
473 idir(0,j)=0
474 jpV(0,j)=j
475 jpH(0,j)=1
476 enddo
477
478 ccc DP ------------------------------>
479 do 111 j=1,nseq2
480 do 222 i=1,nseq1
481 ccc D=VAL(i-1,j-1)+SCORE(i,j)--------------->
482 D=val(i-1,j-1)+score(i,j) !from diagonal, val(i,j) is val(i-1,j-1)
483 ccc H=H+gap_open ------->
484 jpH(i,j)=1
485 val1=val(i-1,j)+gap_open !gap_open from both D and V
486 val2=preH(i-1,j)+gap_extn !gap_extn from horizontal
487 if(val1.gt.val2) then !last step from D or V
488 H=val1
489 else !last step from H
490 H=val2
491 if(i.gt.1)jpH(i,j)=jpH(i-1,j)+1 !record long-gap
492 endif
493 ccc V=V+gap_open --------->
494 jpV(i,j)=1
495 val1=val(i,j-1)+gap_open
496 val2=preV(i,j-1)+gap_extn
497 if(val1.gt.val2) then
498 V=val1
499 else
500 V=val2
501 if(j.gt.1)jpV(i,j)=jpV(i,j-1)+1
502 endif
503 preH(i,j)=H !unaccepted H
504 preV(i,j)=V !unaccepted V
505
506 if(D.gt.H.and.D.gt.V)then
507 idir(i,j)=1
508 val(i,j)=D
509 elseif(H.gt.V)then
510 idir(i,j)=2
511 val(i,j)=H
512 else
513 idir(i,j)=3
514 val(i,j)=V
515 endif
516 222 continue
517 111 continue
518 score0=val(nseq1,nseq2) !alignment score
519
520 c tracing back the pathway:
521 do j=1,nseq2
522 j2i(j)=-1 !all are not aligned
523 enddo
524 i=nseq1
525 j=nseq2
526 do while(i.gt.0.and.j.gt.0)
527 if(idir(i,j).eq.1)then !from diagonal
528 j2i(j)=i
529 i=i-1
530 j=j-1
531 elseif(idir(i,j).eq.2)then !from horizonal
532 it=jpH(i,j)
533 do me=1,it
534 if(i.gt.0) then
535 i=i-1
536 endif
537 enddo
538 else
539 it=jpV(i,j)
540 do me=1,it
541 if(j.gt.0) then
542 j=j-1
543 endif
544 enddo
545 endif
546 enddo
547
548 *^^^^^^^^^^^DP finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
549 return
550 end
551
552 ********************************************************************
553 * This is an alternative implementation of Needleman-Wunsch dynamic program
554 * (by Y. Zhang 2005)
555 * 1. Count two-layer iteration and multiple-gaps
556 * 2. The gap penality W(k)=Go+Ge*k1+Ge*k2 if gap open on both sequences
557 *
558 * Input: score(i,j), gap_open, gap_extn
559 * Output: j2i(j)
560 * idir(i,j)=1,2,3, from diagonal, horizontal, vertical
561 * val(i,j) is the cumulative score of (i,j)
562 ********************************************************************
563 subroutine DPalt(score0)
564 PARAMETER(ndim=6000)
565 common/dpc/score(ndim,ndim),gap_open,gap_extn,j2i(ndim)
566 & ,nseq1,nseq2
567
568 dimension val(0:ndim,0:ndim),idir(0:ndim,0:ndim)
569 dimension preV(0:ndim,0:ndim),preH(0:ndim,0:ndim),
570 & preD(0:ndim,0:ndim)
571 dimension idirH(0:ndim,0:ndim),idirV(0:ndim,0:ndim)
572
573 ccc initializations --------------->
574 val(0,0)=0.0
575 do i=1,nseq1
576 val(i,0)=0
577 idir(i,0)=0
578 preD(i,0)=0.0
579 preH(i,0)=-1000.0
580 preV(i,0)=-1000.0
581 enddo
582 do j=1,nseq2
583 val(0,j)=0
584 idir(0,j)=0
585 preD(0,j)=0.0
586 preH(0,j)=-1000.0
587 preV(0,j)=-1000.0
588 enddo
589
590 ccc DP ------------------------------>
591 do 111 j=1,nseq2
592 do 222 i=1,nseq1
593 ccc preD=VAL(i-1,j-1)+SCORE(i,j)--------------->
594 preD(i,j)=val(i-1,j-1)+score(i,j)
595 ccc preH: pre-accepted H----------------------->
596 D=preD(i-1,j)+gap_open
597 H=preH(i-1,j)+gap_extn
598 V=preV(i-1,j)+gap_extn
599 if(D.gt.H.and.D.gt.V)then
600 preH(i,j)=D
601 idirH(i-1,j)=1
602 elseif(H.gt.V)then
603 preH(i,j)=H
604 idirH(i-1,j)=2
605 else
606 preH(i,j)=V
607 idirH(i-1,j)=3
608 endif
609 ccc preV: pre-accepted V----------------------->
610 D=preD(i,j-1)+gap_open
611 H=preH(i,j-1)+gap_extn
612 V=preV(i,j-1)+gap_extn
613 if(D.gt.H.and.D.gt.V)then
614 preV(i,j)=D
615 idirV(i,j-1)=1
616 elseif(H.gt.V)then
617 preV(i,j)=H
618 idirV(i,j-1)=2
619 else
620 preV(i,j)=V
621 idirV(i,j-1)=3
622 endif
623
624 ccc decide idir(i,j)----------->
625 if(preD(i,j).gt.preH(i,j).and.preD(i,j).gt.preV(i,j))then
626 idir(i,j)=1
627 val(i,j)=preD(i,j)
628 elseif(preH(i,j).gt.preV(i,j))then
629 idir(i,j)=2
630 val(i,j)=preH(i,j)
631 else
632 idir(i,j)=3
633 val(i,j)=preV(i,j)
634 endif
635 222 continue
636 111 continue
637 score0=val(nseq1,nseq2) !alignment score
638
639 ccc tracing back the pathway:
640 do j=1,nseq2
641 j2i(j)=-1 !all are not aligned
642 enddo
643 i=nseq1
644 j=nseq2
645 do while(i.gt.0.and.j.gt.0)
646 if(idir(i,j).eq.1)then !from diagonal
647 j2i(j)=i
648 i=i-1
649 j=j-1
650 elseif(idir(i,j).eq.2)then
651 i=i-1
652 idir(i,j)=idirH(i,j)
653 else
654 j=j-1
655 idir(i,j)=idirV(i,j)
656 endif
657 enddo
658
659 *^^^^^^^^^^^DP finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
660 return
661 end
662
663 ********************************************************************
664 * read matrix
665 *
666 subroutine matrix
667 parameter(naa=24) !number of amino acid
668 common/matra/imut(naa,naa) !b,z,x are additional
669
670 * following from PAM30:
671 imut(1,1)=6
672 imut(1,2)=-7
673 imut(1,3)=-4
674 imut(1,4)=-3
675 imut(1,5)=-6
676 imut(1,6)=-4
677 imut(1,7)=-2
678 imut(1,8)=-2
679 imut(1,9)=-7
680 imut(1,10)=-5
681 imut(1,11)=-6
682 imut(1,12)=-7
683 imut(1,13)=-5
684 imut(1,14)=-8
685 imut(1,15)=-2
686 imut(1,16)=0
687 imut(1,17)=-1
688 imut(1,18)=-13
689 imut(1,19)=-8
690 imut(1,20)=-2
691 imut(1,21)=-3
692 imut(1,22)=-3
693 imut(1,23)=-1
694 imut(1,24)=-17
695 imut(2,1)=-7
696 imut(2,2)=8
697 imut(2,3)=-6
698 imut(2,4)=-10
699 imut(2,5)=-8
700 imut(2,6)=-2
701 imut(2,7)=-9
702 imut(2,8)=-9
703 imut(2,9)=-2
704 imut(2,10)=-5
705 imut(2,11)=-8
706 imut(2,12)=0
707 imut(2,13)=-4
708 imut(2,14)=-9
709 imut(2,15)=-4
710 imut(2,16)=-3
711 imut(2,17)=-6
712 imut(2,18)=-2
713 imut(2,19)=-10
714 imut(2,20)=-8
715 imut(2,21)=-7
716 imut(2,22)=-4
717 imut(2,23)=-1
718 imut(2,24)=-17
719 imut(3,1)=-4
720 imut(3,2)=-6
721 imut(3,3)=8
722 imut(3,4)=2
723 imut(3,5)=-11
724 imut(3,6)=-3
725 imut(3,7)=-2
726 imut(3,8)=-3
727 imut(3,9)=0
728 imut(3,10)=-5
729 imut(3,11)=-7
730 imut(3,12)=-1
731 imut(3,13)=-9
732 imut(3,14)=-9
733 imut(3,15)=-6
734 imut(3,16)=0
735 imut(3,17)=-2
736 imut(3,18)=-8
737 imut(3,19)=-4
738 imut(3,20)=-8
739 imut(3,21)=6
740 imut(3,22)=-3
741 imut(3,23)=-1
742 imut(3,24)=-17
743 imut(4,1)=-3
744 imut(4,2)=-10
745 imut(4,3)=2
746 imut(4,4)=8
747 imut(4,5)=-14
748 imut(4,6)=-2
749 imut(4,7)=2
750 imut(4,8)=-3
751 imut(4,9)=-4
752 imut(4,10)=-7
753 imut(4,11)=-12
754 imut(4,12)=-4
755 imut(4,13)=-11
756 imut(4,14)=-15
757 imut(4,15)=-8
758 imut(4,16)=-4
759 imut(4,17)=-5
760 imut(4,18)=-15
761 imut(4,19)=-11
762 imut(4,20)=-8
763 imut(4,21)=6
764 imut(4,22)=1
765 imut(4,23)=-1
766 imut(4,24)=-17
767 imut(5,1)=-6
768 imut(5,2)=-8
769 imut(5,3)=-11
770 imut(5,4)=-14
771 imut(5,5)=10
772 imut(5,6)=-14
773 imut(5,7)=-14
774 imut(5,8)=-9
775 imut(5,9)=-7
776 imut(5,10)=-6
777 imut(5,11)=-15
778 imut(5,12)=-14
779 imut(5,13)=-13
780 imut(5,14)=-13
781 imut(5,15)=-8
782 imut(5,16)=-3
783 imut(5,17)=-8
784 imut(5,18)=-15
785 imut(5,19)=-4
786 imut(5,20)=-6
787 imut(5,21)=-12
788 imut(5,22)=-14
789 imut(5,23)=-1
790 imut(5,24)=-17
791 imut(6,1)=-4
792 imut(6,2)=-2
793 imut(6,3)=-3
794 imut(6,4)=-2
795 imut(6,5)=-14
796 imut(6,6)=8
797 imut(6,7)=1
798 imut(6,8)=-7
799 imut(6,9)=1
800 imut(6,10)=-8
801 imut(6,11)=-5
802 imut(6,12)=-3
803 imut(6,13)=-4
804 imut(6,14)=-13
805 imut(6,15)=-3
806 imut(6,16)=-5
807 imut(6,17)=-5
808 imut(6,18)=-13
809 imut(6,19)=-12
810 imut(6,20)=-7
811 imut(6,21)=-3
812 imut(6,22)=6
813 imut(6,23)=-1
814 imut(6,24)=-17
815 imut(7,1)=-2
816 imut(7,2)=-9
817 imut(7,3)=-2
818 imut(7,4)=2
819 imut(7,5)=-14
820 imut(7,6)=1
821 imut(7,7)=8
822 imut(7,8)=-4
823 imut(7,9)=-5
824 imut(7,10)=-5
825 imut(7,11)=-9
826 imut(7,12)=-4
827 imut(7,13)=-7
828 imut(7,14)=-14
829 imut(7,15)=-5
830 imut(7,16)=-4
831 imut(7,17)=-6
832 imut(7,18)=-17
833 imut(7,19)=-8
834 imut(7,20)=-6
835 imut(7,21)=1
836 imut(7,22)=6
837 imut(7,23)=-1
838 imut(7,24)=-17
839 imut(8,1)=-2
840 imut(8,2)=-9
841 imut(8,3)=-3
842 imut(8,4)=-3
843 imut(8,5)=-9
844 imut(8,6)=-7
845 imut(8,7)=-4
846 imut(8,8)=6
847 imut(8,9)=-9
848 imut(8,10)=-11
849 imut(8,11)=-10
850 imut(8,12)=-7
851 imut(8,13)=-8
852 imut(8,14)=-9
853 imut(8,15)=-6
854 imut(8,16)=-2
855 imut(8,17)=-6
856 imut(8,18)=-15
857 imut(8,19)=-14
858 imut(8,20)=-5
859 imut(8,21)=-3
860 imut(8,22)=-5
861 imut(8,23)=-1
862 imut(8,24)=-17
863 imut(9,1)=-7
864 imut(9,2)=-2
865 imut(9,3)=0
866 imut(9,4)=-4
867 imut(9,5)=-7
868 imut(9,6)=1
869 imut(9,7)=-5
870 imut(9,8)=-9
871 imut(9,9)=9
872 imut(9,10)=-9
873 imut(9,11)=-6
874 imut(9,12)=-6
875 imut(9,13)=-10
876 imut(9,14)=-6
877 imut(9,15)=-4
878 imut(9,16)=-6
879 imut(9,17)=-7
880 imut(9,18)=-7
881 imut(9,19)=-3
882 imut(9,20)=-6
883 imut(9,21)=-1
884 imut(9,22)=-1
885 imut(9,23)=-1
886 imut(9,24)=-17
887 imut(10,1)=-5
888 imut(10,2)=-5
889 imut(10,3)=-5
890 imut(10,4)=-7
891 imut(10,5)=-6
892 imut(10,6)=-8
893 imut(10,7)=-5
894 imut(10,8)=-11
895 imut(10,9)=-9
896 imut(10,10)=8
897 imut(10,11)=-1
898 imut(10,12)=-6
899 imut(10,13)=-1
900 imut(10,14)=-2
901 imut(10,15)=-8
902 imut(10,16)=-7
903 imut(10,17)=-2
904 imut(10,18)=-14
905 imut(10,19)=-6
906 imut(10,20)=2
907 imut(10,21)=-6
908 imut(10,22)=-6
909 imut(10,23)=-1
910 imut(10,24)=-17
911 imut(11,1)=-6
912 imut(11,2)=-8
913 imut(11,3)=-7
914 imut(11,4)=-12
915 imut(11,5)=-15
916 imut(11,6)=-5
917 imut(11,7)=-9
918 imut(11,8)=-10
919 imut(11,9)=-6
920 imut(11,10)=-1
921 imut(11,11)=7
922 imut(11,12)=-8
923 imut(11,13)=1
924 imut(11,14)=-3
925 imut(11,15)=-7
926 imut(11,16)=-8
927 imut(11,17)=-7
928 imut(11,18)=-6
929 imut(11,19)=-7
930 imut(11,20)=-2
931 imut(11,21)=-9
932 imut(11,22)=-7
933 imut(11,23)=-1
934 imut(11,24)=-17
935 imut(12,1)=-7
936 imut(12,2)=0
937 imut(12,3)=-1
938 imut(12,4)=-4
939 imut(12,5)=-14
940 imut(12,6)=-3
941 imut(12,7)=-4
942 imut(12,8)=-7
943 imut(12,9)=-6
944 imut(12,10)=-6
945 imut(12,11)=-8
946 imut(12,12)=7
947 imut(12,13)=-2
948 imut(12,14)=-14
949 imut(12,15)=-6
950 imut(12,16)=-4
951 imut(12,17)=-3
952 imut(12,18)=-12
953 imut(12,19)=-9
954 imut(12,20)=-9
955 imut(12,21)=-2
956 imut(12,22)=-4
957 imut(12,23)=-1
958 imut(12,24)=-17
959 imut(13,1)=-5
960 imut(13,2)=-4
961 imut(13,3)=-9
962 imut(13,4)=-11
963 imut(13,5)=-13
964 imut(13,6)=-4
965 imut(13,7)=-7
966 imut(13,8)=-8
967 imut(13,9)=-10
968 imut(13,10)=-1
969 imut(13,11)=1
970 imut(13,12)=-2
971 imut(13,13)=11
972 imut(13,14)=-4
973 imut(13,15)=-8
974 imut(13,16)=-5
975 imut(13,17)=-4
976 imut(13,18)=-13
977 imut(13,19)=-11
978 imut(13,20)=-1
979 imut(13,21)=-10
980 imut(13,22)=-5
981 imut(13,23)=-1
982 imut(13,24)=-17
983 imut(14,1)=-8
984 imut(14,2)=-9
985 imut(14,3)=-9
986 imut(14,4)=-15
987 imut(14,5)=-13
988 imut(14,6)=-13
989 imut(14,7)=-14
990 imut(14,8)=-9
991 imut(14,9)=-6
992 imut(14,10)=-2
993 imut(14,11)=-3
994 imut(14,12)=-14
995 imut(14,13)=-4
996 imut(14,14)=9
997 imut(14,15)=-10
998 imut(14,16)=-6
999 imut(14,17)=-9
1000 imut(14,18)=-4
1001 imut(14,19)=2
1002 imut(14,20)=-8
1003 imut(14,21)=-10
1004 imut(14,22)=-13
1005 imut(14,23)=-1
1006 imut(14,24)=-17
1007 imut(15,1)=-2
1008 imut(15,2)=-4
1009 imut(15,3)=-6
1010 imut(15,4)=-8
1011 imut(15,5)=-8
1012 imut(15,6)=-3
1013 imut(15,7)=-5
1014 imut(15,8)=-6
1015 imut(15,9)=-4
1016 imut(15,10)=-8
1017 imut(15,11)=-7
1018 imut(15,12)=-6
1019 imut(15,13)=-8
1020 imut(15,14)=-10
1021 imut(15,15)=8
1022 imut(15,16)=-2
1023 imut(15,17)=-4
1024 imut(15,18)=-14
1025 imut(15,19)=-13
1026 imut(15,20)=-6
1027 imut(15,21)=-7
1028 imut(15,22)=-4
1029 imut(15,23)=-1
1030 imut(15,24)=-17
1031 imut(16,1)=0
1032 imut(16,2)=-3
1033 imut(16,3)=0
1034 imut(16,4)=-4
1035 imut(16,5)=-3
1036 imut(16,6)=-5
1037 imut(16,7)=-4
1038 imut(16,8)=-2
1039 imut(16,9)=-6
1040 imut(16,10)=-7
1041 imut(16,11)=-8
1042 imut(16,12)=-4
1043 imut(16,13)=-5
1044 imut(16,14)=-6
1045 imut(16,15)=-2
1046 imut(16,16)=6
1047 imut(16,17)=0
1048 imut(16,18)=-5
1049 imut(16,19)=-7
1050 imut(16,20)=-6
1051 imut(16,21)=-1
1052 imut(16,22)=-5
1053 imut(16,23)=-1
1054 imut(16,24)=-17
1055 imut(17,1)=-1
1056 imut(17,2)=-6
1057 imut(17,3)=-2
1058 imut(17,4)=-5
1059 imut(17,5)=-8
1060 imut(17,6)=-5
1061 imut(17,7)=-6
1062 imut(17,8)=-6
1063 imut(17,9)=-7
1064 imut(17,10)=-2
1065 imut(17,11)=-7
1066 imut(17,12)=-3
1067 imut(17,13)=-4
1068 imut(17,14)=-9
1069 imut(17,15)=-4
1070 imut(17,16)=0
1071 imut(17,17)=7
1072 imut(17,18)=-13
1073 imut(17,19)=-6
1074 imut(17,20)=-3
1075 imut(17,21)=-3
1076 imut(17,22)=-6
1077 imut(17,23)=-1
1078 imut(17,24)=-17
1079 imut(18,1)=-13
1080 imut(18,2)=-2
1081 imut(18,3)=-8
1082 imut(18,4)=-15
1083 imut(18,5)=-15
1084 imut(18,6)=-13
1085 imut(18,7)=-17
1086 imut(18,8)=-15
1087 imut(18,9)=-7
1088 imut(18,10)=-14
1089 imut(18,11)=-6
1090 imut(18,12)=-12
1091 imut(18,13)=-13
1092 imut(18,14)=-4
1093 imut(18,15)=-14
1094 imut(18,16)=-5
1095 imut(18,17)=-13
1096 imut(18,18)=13
1097 imut(18,19)=-5
1098 imut(18,20)=-15
1099 imut(18,21)=-10
1100 imut(18,22)=-14
1101 imut(18,23)=-1
1102 imut(18,24)=-17
1103 imut(19,1)=-8
1104 imut(19,2)=-10
1105 imut(19,3)=-4
1106 imut(19,4)=-11
1107 imut(19,5)=-4
1108 imut(19,6)=-12
1109 imut(19,7)=-8
1110 imut(19,8)=-14
1111 imut(19,9)=-3
1112 imut(19,10)=-6
1113 imut(19,11)=-7
1114 imut(19,12)=-9
1115 imut(19,13)=-11
1116 imut(19,14)=2
1117 imut(19,15)=-13
1118 imut(19,16)=-7
1119 imut(19,17)=-6
1120 imut(19,18)=-5
1121 imut(19,19)=10
1122 imut(19,20)=-7
1123 imut(19,21)=-6
1124 imut(19,22)=-9
1125 imut(19,23)=-1
1126 imut(19,24)=-17
1127 imut(20,1)=-2
1128 imut(20,2)=-8
1129 imut(20,3)=-8
1130 imut(20,4)=-8
1131 imut(20,5)=-6
1132 imut(20,6)=-7
1133 imut(20,7)=-6
1134 imut(20,8)=-5
1135 imut(20,9)=-6
1136 imut(20,10)=2
1137 imut(20,11)=-2
1138 imut(20,12)=-9
1139 imut(20,13)=-1
1140 imut(20,14)=-8
1141 imut(20,15)=-6
1142 imut(20,16)=-6
1143 imut(20,17)=-3
1144 imut(20,18)=-15
1145 imut(20,19)=-7
1146 imut(20,20)=7
1147 imut(20,21)=-8
1148 imut(20,22)=-6
1149 imut(20,23)=-1
1150 imut(20,24)=-17
1151 imut(21,1)=-3
1152 imut(21,2)=-7
1153 imut(21,3)=6
1154 imut(21,4)=6
1155 imut(21,5)=-12
1156 imut(21,6)=-3
1157 imut(21,7)=1
1158 imut(21,8)=-3
1159 imut(21,9)=-1
1160 imut(21,10)=-6
1161 imut(21,11)=-9
1162 imut(21,12)=-2
1163 imut(21,13)=-10
1164 imut(21,14)=-10
1165 imut(21,15)=-7
1166 imut(21,16)=-1
1167 imut(21,17)=-3
1168 imut(21,18)=-10
1169 imut(21,19)=-6
1170 imut(21,20)=-8
1171 imut(21,21)=6
1172 imut(21,22)=0
1173 imut(21,23)=-1
1174 imut(21,24)=-17
1175 imut(22,1)=-3
1176 imut(22,2)=-4
1177 imut(22,3)=-3
1178 imut(22,4)=1
1179 imut(22,5)=-14
1180 imut(22,6)=6
1181 imut(22,7)=6
1182 imut(22,8)=-5
1183 imut(22,9)=-1
1184 imut(22,10)=-6
1185 imut(22,11)=-7
1186 imut(22,12)=-4
1187 imut(22,13)=-5
1188 imut(22,14)=-13
1189 imut(22,15)=-4
1190 imut(22,16)=-5
1191 imut(22,17)=-6
1192 imut(22,18)=-14
1193 imut(22,19)=-9
1194 imut(22,20)=-6
1195 imut(22,21)=0
1196 imut(22,22)=6
1197 imut(22,23)=-1
1198 imut(22,24)=-17
1199 imut(23,1)=-1
1200 imut(23,2)=-1
1201 imut(23,3)=-1
1202 imut(23,4)=-1
1203 imut(23,5)=-1
1204 imut(23,6)=-1
1205 imut(23,7)=-1
1206 imut(23,8)=-1
1207 imut(23,9)=-1
1208 imut(23,10)=-1
1209 imut(23,11)=-1
1210 imut(23,12)=-1
1211 imut(23,13)=-1
1212 imut(23,14)=-1
1213 imut(23,15)=-1
1214 imut(23,16)=-1
1215 imut(23,17)=-1
1216 imut(23,18)=-1
1217 imut(23,19)=-1
1218 imut(23,20)=-1
1219 imut(23,21)=-1
1220 imut(23,22)=-1
1221 imut(23,23)=-1
1222 imut(23,24)=-17
1223 imut(24,1)=-17
1224 imut(24,2)=-17
1225 imut(24,3)=-17
1226 imut(24,4)=-17
1227 imut(24,5)=-17
1228 imut(24,6)=-17
1229 imut(24,7)=-17
1230 imut(24,8)=-17
1231 imut(24,9)=-17
1232 imut(24,10)=-17
1233 imut(24,11)=-17
1234 imut(24,12)=-17
1235 imut(24,13)=-17
1236 imut(24,14)=-17
1237 imut(24,15)=-17
1238 imut(24,16)=-17
1239 imut(24,17)=-17
1240 imut(24,18)=-17
1241 imut(24,19)=-17
1242 imut(24,20)=-17
1243 imut(24,21)=-17
1244 imut(24,22)=-17
1245 imut(24,23)=-17
1246 imut(24,24)=1
1247
1248 return
1249 end
1250
1251 function upper(A)
1252 CHARACTER A,upper
1253 IF(A.LE.'z'.and.A.GE.'a')then
1254 A=CHAR(ICHAR(A)-32)
1255 endif
1256 upper=A
1257 RETURN
1258 END