|
0
|
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
|
