c c Los Alamos National Laboratory c c Copyright, 1989. The Regents of the University of California. c This software was produced under a U.S. Government contract c (W-7405-ENG-36) by Los Alamos National Laboratory, which is c operated by the University of California for the U.S. Department c of Energy. The U.S. Government is licensed to use, reproduce, c and distribute this software. Permission is granted to the c public to copy and use this software without charge, provided c that this Notice and any statement of authorship are reproduced c on all copies. Neither the Government nor the University makes c any warranty, express or implied, or assumes any liability or c responsibility for the use of this software. c c c seqf1, CFT77 Fortran version, based on VAXTRAN version c c c CFT77: VAX, SUN, Multiflow dbd c c modified (11/08/83) by M.Kanehisa c modified (12/30/83) by M.Kanehisa c modified (07/11/84) by C.Burks for longer query sequences. c note that, while increasing query length (ndim), c the following relations must hold -- c ldim=2*ndim c lwork=.5*ndim c modified (04/05/87) by C.Burks, added some header comments for clarity c c very modified by Dan Davison 10/16/87 for VMS; also handles up to c 80,000 character sequences. This should be modified to use the VAX c LIB$GET_VM routine (the moral equivalent of the CALL MEMADJ() routine c in the CFTLIB library). My changes are erratically annotated; for best c results, use the VMS "diff" program to compare this program with the Cray c program of the same name, SEQF1.C (<- this "C" refers to Cray, not the c language!) c c modified (4/27/88) by D. Davison for Sun, Multiflow c modified 5/3/89 by dbd for Y/MP & generic ftn 77 c c a fortran 77 program to search the dna database for good local c homologies. first the program locates separate diagonal positions c containing candidates of good matches by the k-tuple method (wilbur c and lipman, proc. natl. acad. sci. usa 80, 726-730, 1973). then c the needleman-wunsch-sellers algorithm (goad and kanehisa, nucleic c acids res. 10, 247-263, 1982) is performed to find the best local c homology along each of these diagonals. an option also exists to c use the vectorizable nws algorithm for the initial scan of finding c the position of the best local homology. c c note that all GenBank-like files used as input for this program c must be in all-lower-case. c c input parameters c --------------- c sequence data c iseq(1)...iseq(n), sequence to be compared (length n) c imin, start value in real numbering scheme c iname, sequence identifier c c measure of similarity and output control c dmat, weight for a match (default, -1) c drep, weight for a mismatch (default, 2) c del, weight for a deletion (default, 3) c maxd, maximum homology score for each alignment to be printed c ktpl, k-tuple size for optional fast scan (default, 4) c nw, window size for searching along the diagonal (default, 20) c lwid, number of characters per line at the terminal (60 to 120) c c - if ktpl is zero, the nws algorithm is used to scan the sequence c for finding the location of the best local homology. c c machine dependence c ------------------ c bit and character manipulations (read instructions in program seqh) c lbit, number of bits per word c c storage requirements c ndim, maximum length of the query sequence c mdim, maximum length of the library sequence c ldim, size of work areas c c - the total storage area required is the sum of the following: c iseq ... ndim c kseq ... ndim c jseq ... mdim c dist ... max( 2*(ndim+1), ldim ) c path ... max( 2*ldim, ((2*nwmax)/lpat)*ndim ) c c - note the following storage areas are equivalenced. c nwint <- next -> c ktuple <- dist -> <- diag -> <- last -> c nwacr <-- dist --> c nwvect <---- zist ----> c nwband <-- dist --> <------ path -------> c c - the size of mdim can be dynamically adjusted in cdc and cray. c c files used c ---------- c unit 1, input file of query sequences c unit 2, input file of library sequences (database to be searched) c unit 6, output file of results (default, user's terminal) c c----------------------------------------------------------------------- c c database search program seqf(tty,input=tty,output=tty) c program seqf(tty,input=tty,output=tty,tape7=tty,tape6=tty) c program seqf implicit integer(a-y) parameter (ndim=2000,mdim=8000,ldim=ndim*2,ktmax=4,nwmax=100) character iname*10,jname*70,defin*80,ifile*10,jfile*10,ans*1 character list(100)*10 character kfile*80 character nucl(5), cmpl(5) integer dis(100),loc(100) logical logvar * level 2,// common jseq(mdim) common /metric/dmat,drep,del,imin1,lwid common /bits/lpat,lres common /base/nucl,cmpl common /info/jname,defin data dmat,drep,del/-1,2,3/,nw/20/ data nucl/'a','c','g','t','x'/,cmpl/'t','g','c','a','x'/ *** * vax, multiflow, sun * data lbit/32/ * cdc7600 * data lbit/60/ * cray-1 data lbit/64/ *** lsize=mdim lpat=lbit/3 lres=lbit-lpat*3 print *,' ' print *,'los alamos sequence analysis system -- dna data bank sear .ch' print *,' ' 1 print 101 101 format(' input file or help ') read 102,ifile 102 format(a) if(ifile.eq.'help') go to 91 c dd 10/15/87 inquire(file=ifile,exist=logvar) if(.NOT.logvar) then print 220, ifile goto 92 else open(unit=1,file=ifile,status='UNKNOWN') endif read(1,102) iname if(iname.eq.'locus'.or.iname.eq.'LOCUS') then lf=1 rewind 1 else lf=2 rewind 1 endif 2 print 103 103 format(' library file or help ') read 102,jfile if(jfile.eq.'help') go to 93 c if queryfile ne jfile, open the library file if(jfile.eq.ifile) then goto 3 c call addlun(2,1) else logvar = .TRUE. inquire(file=jfile,exist=logvar) if(.NOT.logvar) then goto 94 else open(unit=2,file=jfile,status='UNKNOWN') endif endif 3 print *,' output file ([cr] for your terminal) ' kfile='' read 102,kfile if(kfile.ne.'') then c dd 9/29/87 c call addlun(6,7) open(unit=7,file=kfile, status='new') c else c call assign(6,kfile,-1) endif maxd=-15 ktpl=4 lwid=70 c added by dan davison, at the request of Christian Burks print 106,maxd,ktpl,nw,lwid,dmat,drep,del c end addition. just prints out the values of the alignment parameters c before asking for changes! 995 print 105 105 format(' change parameters? (y/n/h) ') read 102,ans if(ans.eq.'y') go to 4 if(ans.eq.'n') go to 5 if(ans.ne.'h') go to 995 print 106,maxd,ktpl,nw,lwid,dmat,drep,del 106 format(/' maxd, alignments with more negative scores are printed ( .current value:',i4,')'/ .' ktpl, k-tuple size for initial rapid scan (current valuue:',i2, .')'/ .' nw, window size of diagonal band (current value:',i3,')'/ .' lwid, terminal line width (current value:',i4,')'/ .' dmat, negative weight for a match (current value:',i3,')'/ .' drep, positive weight for a mismatch (current value:',i2,')'/ .' del, positive weight for a deletion (current value:',i2,')'/) go to 3 4 print 107 107 format(' enter maxd,ktpl,nw,lwid,dmat,drep,del ') read *,maxd,ktpl,nw,lwid,dmat,drep,del if(ktpl.gt.ktmax) ktpl=ktmax if(nw.gt.nwmax) nw=nwmax lwid=(lwid/10)*10 if(lwid.lt.60) lwid=60 if(lwid.gt.120) lwid=120 5 print 109,maxd,ktpl,nw,lwid,dmat,drep,del 109 format(/' parameters are:'/3x,'maxd =',i4, . 4x,'ktpl =',i2,4x,'nw =',i3,4x,'lwid =',i5/ . 3x,'dmat =',i3,4x,'drep =',i3,4x,'del =',i3/) 11 print 111 111 format(' query sequence ') iname='end' read 102,iname if(iname.eq.'end'.or.iname.eq.'END') goto 99 call seqin(lf,iname,imin,n,jseq,*11) if(n.gt.ndim) go to 95 n1=n+1 imin1=imin-1 imax=imin1+n print 110,iname,imin,imax,n if(kfile.ne.'') write(7,110) iname,imin,imax,n rewind 2 kount=0 length=0 best=0 call nwint(n,jseq,ktpl) if(ktpl.ne.0) go to 31 21 call genin(lsize,m,jseq,jname,defin,*41) kount=kount+1 length=length+m *** * non-vectorized nws algorithm VAX, multiflow * call nwacr(m,jseq,dmin,ii,jj) * vectorized nws algorithm call nwvect(m,jseq,dmin,ii,jj) *** if(dmin.lt.best) then best=dmin l=1 list(1)=jname(13:22) else if(dmin.eq.best.and.l.lt.100) then l=l+1 list(l)=jname(13:22) endif if(dmin.gt.maxd) go to 21 call nwband(ii,jj,jseq,dmin,dmin,maxd,nw) go to 21 31 call genin(lsize,m,jseq,jname,defin,*41) kount=kount+1 length=length+m call ktuple(m,jseq,dmin,maxd,nt,dis,loc,nw) if(dmin.lt.best) then best=dmin l=1 list(1)=jname(13:22) else if(dmin.eq.best.and.l.lt.100) then l=l+1 list(l)=jname(13:22) endif if(nt.le.0) go to 31 do 32 i=1,nt ij=loc(i) if(ij.le.m) then ii=n jj=ij else ii=n+m-ij jj=m endif 32 call nwband(ii,jj,jseq,-999999,dis(i),maxd,nw) go to 31 41 print 120,length,kount if(ktpl.eq.0) then print 130,best,(list(i),i=1,l) if(kfile.ne.'') write(7,130) best,(list(i),i=1,l) else print 140,best,(list(i),i=1,l) if(kfile.ne.'') write(7,140) best,(list(i),i=1,l) endif 43 print 112 112 format(/' continue? (y/n)') read 102,ans if(ans.eq.'y') go to 11 if(ans.eq.'n') go to 99 go to 43 91 print 200 go to 1 92 print 220,ifile go to 1 93 print 210 go to 2 94 print 230,jfile go to 2 95 print *,'*** query sequence too long.' go to 11 99 stop 110 format(' query sequence',a12,' from',i8,' to',i8,' total',i7) 120 format(/' total of ',i8,' bases in ',i5,' sequences searched') 130 format(/' the best homology score ',i6, .' was found in:'/5(a13)) 140 format(/' the best homology score (uncorrected) ',i6, .' was found in:'/5(a13)) 200 format(/' input file of query sequences must be either one of the .following:'/ .' 1. genbank format (if the first line is a locus line)'/ .' 2. stanford seq format (otherwise)'/) 210 format(/' choose one of the following library files:'/ .' dna entire genbank database'/ .' dnam mammalian sequences'/ .' dnan non-mammalian vertebrate sequences'/ .' dnai invertebrate sequences'/ .' dnap plant sequences'/ .' dnao organelle sequences'/ .' dnab bacteria sequences'/ .' dnas structural rna sequences'/ .' dnav virus sequences'/ .' dnaf phage sequences'/ .' dnar recombined or synthetic sequences'/ .' or any file in the genbank format.'/) 220 format(' *** file',a12,' does not exist. try again.') 230 format(' *** file',a12,' does not exist as library. try again.') end c initialization for nws algorithm subroutine nwint(m,jseq,ktpl) implicit integer(a-y) parameter (ndim=2000,mdim=8000,ldim=ndim*2,ktmax=4,nwmax=100) parameter (lwork=250) * level 2,jseq character*1 blank, hyphen, colon character kfile*80 character nucl(5), cmpl(5) dimension iseq(ndim),jseq(mdim),dist(ndim*2+2),diff(5,5) dimension zist(ndim*3+3) dimension kseq(ndim),next(ndim),mark(4**ktmax+1) dimension diag(ldim),last(ldim),work(2,lwork),dis(1),loc(1) c dimension pat(nwmax*2),path(ldim*2) integer pat(nwmax*2),path(ldim*2) dimension ishf(13),jshf(13),inum(13),jnum(13),midl(120) equivalence (dist(1),next(1),zist(1)) equivalence (path(1),diag(1)),(path(ldim+1),last(1)) common /metric/ddmat,ddrep,ddel,imin1,lwid common /bits/lpat,lres common /base/nucl,cmpl common /info/jname,defin character jname*70,defin*80,fm*17,lchr*2 logical final,first save n,n1,iseq,diff,dmat,drep,del,dmin,dmkt save kseq,len,kt,kt1,first data blank,hyphen,colon/' ','-',':'/ iand(i,j)=i.and.j ior(i,j)=i.or.j c shift(i,j)= ishft(i,j) ishft(i,j)=shift(i,j) c print *,' statement functions made' lpat3=lpat*3 lc=lwid-1 dmat=ddmat drep=ddrep del=ddel zmat=dmat zrep=drep zel=del c print *,' constants set in nwint' do 12 i=1,4 do 10 j=1,4 10 diff(i,j)=drep diff(i,5)=0 diff(5,i)=0 12 diff(i,i)=dmat diff(5,5)=0 n=m do 14 i=1,n 14 iseq(i)=jseq(i) n1=n+1 if(ktpl.eq.0) return c initialization for k-tuple search c print *,' initializing k-tuple search' do 20 i=1,n 20 kseq(i)=0 kt=ktpl kt1=kt-1 dmkt=dmat*kt len=4**kt do 22 k=1,len+1 22 mark(k)=0 do 28 i=1,n-kt1 i1=i-1 k=0 do 24 j=1,kt if(iseq(i1+j).eq.5) go to 25 24 k=k*4+iseq(i1+j)-1 go to 26 25 k=len 26 k=k+1 if(mark(k).eq.0) then mark(k)=i else kseq(next(k))=i endif 28 next(k)=i c print *,' exiting ktuple init' return c k-tuple search entry ktuple(m,jseq,ddmin,maxd,nt,dis,loc,nw) c print *,' entry into ktuple' first=.true. msafe=-6 limit=ldim-n+1 lwork1=lwork/(m/limit+1) nm1=n+m-1 do 30 i=1,min0(ldim,nm1) dist(i)=0 diag(i)=0 30 last(i)=0 mk=m-kt1 mt=0 ja=1 31 ja1=ja-1 jb=ja+limit-1 if(jb.gt.mk) then jb=mk limit=nm1-ja1 endif c print *,' beginning loop 34 in nwint' do 34 j=ja,jb j1=j-1 jkt=j-kt k=0 do 32 i=1,kt 32 k=k*4+jseq(j1+i)-1 if(k.gt.len) go to 34 k=k+1 if(mark(k).eq.0) go to 34 i=mark(k) 33 ij=j-i+n-ja1 if(last(ij).eq.0) then dist(ij)=dmkt c print *,' calc distance A' else if(jkt-last(ij).le.0) then dist(ij)=dist(ij)+dmat c print *,' calc distance b' else dist(ij)=min0(dist(ij)+(jkt-last(ij))*drep,0)+dmkt c print *,' calc distance c' endif endif if(dist(ij).lt.diag(ij)) diag(ij)=dist(ij) last(ij)=j if(kseq(i).eq.0) go to 34 i=kseq(i) go to 33 34 continue 35 l=0 c print *,' starting limit check' do 36 i=1,limit if(diag(i).lt.msafe) then l=l+1 if(l.gt.lwork1) go to 39 work(1,l+mt)=diag(i) work(2,l+mt)=i+ja1 endif 36 continue mt=mt+l if(jb.ge.mk) go to 41 ll=ldim-limit do 37 i=1,ll dist(i)=dist(i+limit) diag(i)=diag(i+limit) 37 last(i)=last(i+limit) do 38 i=ll+1,ldim dist(i)=0 diag(i)=0 38 last(i)=0 ja=ja+limit go to 31 39 msafe=msafe-1 * print *,'*** work area overflow. recovering...' go to 35 c locate good alignments on separate diagonals c print *,' checking diagonals' 41 if(mt-1) 42,43,44 42 nt=0 ddmin=0 return 43 nt=1 dis(1)=work(1,1) loc(1)=work(2,1) ddmin=dis(1) return 44 nt=0 nw2=nw/2 ddmin=0 do 48 j=1,100 k=0 dmin=0 do 45 i=1,mt if(work(1,i).ne.0.and.work(1,i).lt.dmin) then c print *,' checking work arrays a' dmin=work(1,i) k=i endif 45 continue if(k.eq.0) return work(1,k)=0 ij=work(2,k) do 46 i=1,mt if(work(1,i).ne.0) then if(iabs(work(2,i)-ij).lt.nw2) then c print *,' checking work arrays b' dmin=dmin+min0(work(1,i)+iabs(work(2,i)-ij)*del,0) work(1,i)=0 endif endif 46 continue if(dmin.lt.ddmin) ddmin=dmin if(dmin.le.maxd) then nt=nt+1 dis(nt)=dmin loc(nt)=ij endif 48 continue return c row-wise nws algorithm (non-vectorizable) entry nwacr(m,jseq,ddmin,ii2,jj2) do 50 i=1,n1+1 50 dist(i)=0 jn1=n1 jn2=0 dmin=0 do 54 j=1,m js=jseq(j) do 52 i=1,n dist(i+1+jn1)=min0(dist(i+jn2)+diff(iseq(i),js), . dist(i+1+jn2)+del,dist(i+jn1)+del,0) if(dist(i+1+jn1).lt.dmin) then dmin=dist(i+1+jn1) i2=i j2=j end if 52 continue jj=jn1 jn1=jn2 jn2=jj 54 continue ddmin=dmin ii2=i2 jj2=j2 first=.true. return c anti-diagonal nws algorithm (vectorizable) entry nwvect(m,jseq,ddmin,ii2,jj2) c print *,' entry into nwvec' n2=n1+n1 n3=n2+n1 m1=m-1 nm1=n+m1 do 60 i=1,n3 60 zist(i)=0. jn=0 jn1=n2 jn2=n1 zmin=0 do 66 k=1,nm1 k1=k+1 ia=max0(1,k-m1) ib=min0(k,n) c print *, ' only important vector loop in the pgm' cdir$ ivdep do 62 i=ia,ib z=cvmgt(zmat,zrep,iseq(i).eq.jseq(k1-i)) zist(i+1+jn)=amin1(zist(i+jn2)+z, . zist(i+1+jn1)+zel,zist(i+jn1)+zel,0.) 62 continue i=ismin(ib-ia+1,zist(ia+1+jn),1)+ia if(zist(i+jn).lt.zmin) then zmin=zist(i+jn) i2=i j2=k1 end if jn=mod(jn+n1,n3) jn1=mod(jn1+n1,n3) jn2=mod(jn2+n1,n3) 66 continue i2=i2-1 j2=j2-i2 ddmin=zmin ii2=i2 jj2=j2 first=.true. return c nws algorithm along the diagonal band entry nwband(nn,mm,jseq,ddmin,dd,maxd,nw) final=.false. c print *, ' entry into nwband' nn1=nn+1 mm1=mm+1 nw1=nw-1 nw2=nw+nw dist(1)=0 dist(1+nn1)=0 do 70 i=1,nw dist(i+1)=0 70 pat(i)=0 do 72 i=nw+2,nn1 dist(i)=1000 72 dist(i+nn1)=1000 k=1 jn=-nw jn1=nn1 jn2=0 dmin=0 do 78 j=1,mm if(j.eq.nn) final=.true. js=jseq(mm1-j) ia=max0(j-nw,1) if(ia.gt.1) dist(ia+jn1)=1000 ib=min0(j+nw1,nn) do 74 i=ia,ib c print *, ' distance check in nwband' d1=dist(i+1+jn2)+del d2=dist(i+jn2)+diff(iseq(nn1-i),js) d4=dist(i+jn1)+del d=min0(d1,d2,d4,0) dist(i+1+jn1)=d pat(i-jn)=0 if(d.eq.0) go to 74 if(d1.eq.d) pat(i-jn)=ior(pat(i-jn),1) if(d2.eq.d) pat(i-jn)=ior(pat(i-jn),2) if(d4.eq.d) pat(i-jn)=ior(pat(i-jn),4) if(d.lt.dmin) then dmin=d i2=i j2=j if(dmin.le.ddmin) final=.true. end if 74 continue i=0 75 patwd=0 do 76 l=1,lpat i=i+1 if (pat(i).lt.0.or.pat(i).gt.7) go to 76 c patwd=(patwd.or.ishft(pat(i),lpat3-l*3)) c print *,' piddling on patwd' c patwd=ior(patwd,ishift(pat(i),lpat3-l*3)) patwd=ior(patwd,ishft(pat(i),lpat3-l*3)) 76 if(i.eq.nw2) go to 77 77 path(k)=patwd k=k+1 if(i.ne.nw2) go to 75 if(final) go to 81 jn=jn+1 jj=jn1 jn1=jn2 jn2=jj 78 continue c trace back 81 if(dmin.gt.maxd) return c print *,' beeginning traceback' if(first) then write(*,110) jname,defin if (kfile.ne.'') write(7,110) jname,defin first=.false. mt=0 endif imin=nn1-i2 jmin=mm1-j2 do 82 i=1,mt 82 if(work(1,i).eq.imin.and.work(2,i).eq.jmin) return mt=mt+1 work(1,mt)=imin work(2,mt)=jmin l=0 i=i2 j=j2 kn=(nw2-1)/lpat+1 k=(j2-1)*kn+(i2-(j2-1-nw)-1)/lpat+1 ks=(mod(i2-(j2-1-nw)-1,lpat)+1)*3 83 if(i.eq.0.or.j.eq.0) go to 87 *** * rightward (negative) shift VAX * jump=iand(ishft(path(k),ks-lpat3),7) * jump=and(lshift(path(k),ks-lpat3),7) * leftward circular shift (cray) jump=iand(ishft(path(k),ks+lres),7) *** if(iand(jump,2).eq.0) go to 84 c if(and(jump,2).eq.0) go to 84 l=l+1 dist(l)=ichar(nucl(iseq(nn1-i))) dist(l+ndim)=ichar(nucl(jseq(mm1-j))) i=i-1 j=j-1 k=k-kn go to 83 84 if(iand(jump,1).eq.0) go to 85 c 84 if(and(jump,1).eq.0) go to 85 l=l+1 dist(l)=ichar(hyphen) dist(l+ndim)=ichar(nucl(jseq(mm1-j))) j=j-1 k=k-kn ks=ks+3 if(ks.le.lpat3) go to 83 ks=3 k=k+1 go to 83 85 if(iand(jump,4).eq.0) go to 87 c 85 if(and(jump,4).eq.0) go to 87 l=l+1 dist(l)=ichar(nucl(iseq(nn1-i))) dist(l+ndim)=ichar(hyphen) i=i-1 ks=ks-3 if(ks.ge.3) go to 83 ks=lpat3 k=k-1 go to 83 87 imin=nn1-i2+imin1 jmin=mm1-j2 imax=nn1-i+imin1 jmax=mm1-j c alignment output nl=l match=0 mismat=0 idel=0 jdel=0 if(imin.lt.0) then inum(1)=(imin/10)*10 ishf(1)=inum(1)-imin+1 if(iabs(inum(1)).ge.10**(min0(5,ishf(1))-1)) then inum(1)=inum(1)+10 ishf(1)=ishf(1)+10 endif else inum(1)=(imin/10+1)*10 ishf(1)=inum(1)-imin+1 if(inum(1).ge.10**min0(5,ishf(1))) then inum(1)=inum(1)+10 ishf(1)=ishf(1)+10 endif endif jnum(1)=(jmin/10+1)*10 jshf(1)=jnum(1)-jmin+1 if(jnum(1).ge.10**min0(5,jshf(1))) then jnum(1)=jnum(1)+10 jshf(1)=jshf(1)+10 endif i=imin j=jmin l1=1 91 l2=l1+lc if(l2.gt.nl) l2=nl k=1 ki=1 kj=1 do 92 l=l1,l2 if(dist(l).eq.dist(l+ndim)) then midl(k)=ichar(colon) match=match+1 else midl(k)=ichar(blank) endif k=k+1 c dd 4/28/88 for f77 if(char(dist(l)).eq.hyphen) then ishf(ki)=ishf(ki)+1 idel=idel+1 else i=i+1 if(i.gt.inum(ki)) then ki=ki+1 inum(ki)=inum(ki-1)+10 ishf(ki)=10 endif endif c dd 4/28/88 for f77 if(char(dist(l+ndim)).eq.hyphen) then jshf(kj)=jshf(kj)+1 jdel=jdel+1 else j=j+1 if(j.gt.jnum(kj)) then kj=kj+1 jnum(kj)=jnum(kj-1)+10 jshf(kj)=10 endif endif 92 continue k1=k-1 ki1=max0(1,ki-1) kj1=max0(1,kj-1) call numbrp(ishf,inum,ki1) write(*,100) (char(dist(l)),l=l1,l2) if (kfile.ne.'') write(7,100) (char(dist(l)),l=l1,l2) write(*,100) (char(midl(k)),k=1,k1) if (kfile.ne.'') write(7,100) (char(midl(k)),k=1,k1) write(*,100) (char(dist(l+ndim)),l=l1,l2) if (kfile.ne.'') write(7,100) (char(dist(l+ndim)),l=l1,l2) call numbrp(jshf,jnum,kj1) if(l2.eq.nl) go to 93 l1=l2+1 ishf(1)=inum(ki)-i+1 jshf(1)=jnum(kj)-j+1 inum(1)=inum(ki) jnum(1)=jnum(kj) if(imin.lt.0) then if(iabs(inum(1)).ge.10**(min0(5,ishf(1))-1)) then inum(1)=inum(1)+10 ishf(1)=ishf(1)+10 endif else if(inum(1).ge.10**min0(5,ishf(1))) then inum(1)=inum(1)+10 ishf(1)=ishf(1)+10 endif endif if(jnum(1).ge.10**min0(5,jshf(1))) then jnum(1)=jnum(1)+10 jshf(1)=jshf(1)+10 endif go to 91 93 mismat=nl-(match+idel+jdel) z=real(match)/real(nl)*100. write(*,120) dmin,dd,z,match,nl if (kfile.ne.' ') write(7,120) dmin,dd,z,match,nl return 100 format(120a1) 110 format(/a70/a80/) 120 format(/6x,'homology score =',i6,4x,'(',i5,' )'/ . 6x,'percent match =',f6.1,'% (',i4,' /',i4,' )'/) end c variable format subroutine numbrp(ishf,inum,ki) integer ishf(1),inum(1) character line(1) character fm*52,lchr*2 fm(1:4)='(' l=2 do 10 k=1,ki write(unit=lchr,fmt='(i2)') ishf(k) fm(l:l+3)='i'//lchr//',' 10 l=l+4 fm(l-1:l-1)=')' write(*,fm) (inum(k),k=1,ki) return end c genbank input subroutine genin(lsize,n,seq,locus,defin,*) parameter (ndim=2000,mdim=8000,ldim=ndim*2,ktmax=4,nwmax=100) character text*80,locus*70,defin*80, standrd*80 * level 2,seq integer seq(mdim) character kfile*80 character nucl(5), cmpl(5) character alphseq(ndim) common /base/nucl, cmpl do 10 i=1,1000 read(2,100,end=21) text 10 if(text(1:10).eq.'locus'.or.text(1:10).eq.'LOCUS') goto 11 go to 21 11 locus=text if(locus(13:22).eq.'end'.or.locus(13:22).eq.'END') goto 21 read(unit=text,fmt='(22x,i7)') n if(n.gt.lsize) then c call memadj(n-lsize,ierr) c lsize=n write(*,120) locus, n, lsize n=lsize else if(lsize.gt.mdim) then c call memadj(mdim-lsize,ierr) lsize=mdim endif read(2,100) standrd read(2,100) defin do 12 i=1,1000 read(2,100,end=21) text 12 if(text(1:10).eq.'origin'.or.text(1:10).eq.'ORIGIN') goto 13 go to 21 13 read(2,110) (alphseq(i),i=1,n) do 16 i=1,n c print *, ' letter read in: ', alphseq(i) if(alphseq(i).eq.'u') alphseq(i)='t' do 14 j=1,4 14 if(alphseq(i).eq.nucl(j)) go to 16 j=5 16 seq(i)=j return 21 return 1 100 format(a) 110 format((9x,6(1x,10a1))) 120 format(//,4x,' Truncated ',a22,', length ', i6,' avail space= ', . i6,'.') end c sequence input subroutine seqin(lf,locus,min,n,seq,*) parameter (ndim=2000, mdim=8000) implicit integer(a-z) * level 2,seq dimension seq(mdim) c line(80) character*10 locus,name,head,ans*1 character nucl(5), cmpl(5) character alphseq(ndim), line(80) common /base/ nucl, cmpl rewind 1 if(lf.eq.1) go to 13 11 read(1,100,end=51) name if(name.ne.locus) go to 11 go to 17 13 read(1,100,end=51) head,name if(head.eq.'locus') then if (name.eq.locus) then goto 15 else goto 13 endif else if(head.eq.'LOCUS') then if (name.eq.locus) then goto 15 else goto 13 endif else goto 13 endif c if(name.ne.locus) go to 13 15 read(1,100) head if(head.eq.'origin') then goto 17 elseif (head.eq.'ORIGIN') then goto 17 else goto 15 endif 17 print 101 101 format(' start,end (1,0 for the entire sequence) ') min=1 max=0 read *,min,max min1=min-1 i=0 if(lf.eq.1) go to 23 21 read(1,110,end=25) (line(j),j=1,80) do 22 j=1,80 if(line(j).eq.'1') go to 25 if(line(j).eq.' ') go to 22 i=i+1 if(i.lt.min) go to 22 alphseq(i-min1)=line(j) if(i.eq.max) go to 25 22 continue go to 21 23 read(1,120) head,(line(j),j=1,60) if(head.ne.' ') go to 25 do 24 j=1,60 if(line(j).eq.' ') go to 25 i=i+1 if(i.lt.min) go to 24 alphseq(i-min1)=line(j) if(i.eq.max) go to 25 24 continue go to 23 25 max=i n=max-min+1 do 32 i=1,n if(alphseq(i).eq.'u') alphseq(i)='t' do 30 j=1,4 30 if(alphseq(i).eq.nucl(j)) go to 32 j=5 32 seq(i)=j 33 print 102 102 format(' complement? (y/n) ') read 110,ans if(ans.eq.'n') return if(ans.ne.'y') go to 33 41 min=-max do 42 i=1,n 42 seq(i)=ichar(cmpl(seq(i))) n1=n+1 n2=n/2 do 44 i=1,n2 s=seq(i) seq(i)=seq(n1-i) 44 seq(n1-i)=s c print *,seq return 51 print 200,locus return 1 100 format(a10,2x,a10) 110 format(80a1) 120 format(a1,8x,6(1x,10a1)) 200 format(' *** ',a10,' was not found *** try again ***') end