2200 lines
63 KiB
Fortran
2200 lines
63 KiB
Fortran
C DIASUB SUBROUTINES FOR DIAGON
|
|
C AUTHOR: RODGER STADEN
|
|
C 16-7-92 Now passing isame to qicks so that main diagonal can be ignored
|
|
C 16-7-92 Now passing isame to cfsq so that main diagonal can be ignored
|
|
C 7-2-90 mhist changed to using reals (rm and rmsq) for sd calc
|
|
C 8-2-90 switched idim1,idim2 in 2 places in mhist, once in qicks
|
|
C 14-2-90 removed call getreg from actout for compatibility with other
|
|
C programs
|
|
C 12-6-90 Changed all occurrences of lh to lf
|
|
C 9-7-90 removed menu routine
|
|
C 13-11-90 replaced all radio by radion
|
|
C 12-1-91 ALIGNM changed call to alignd to use nmax instead of maxseq
|
|
C Changed alignd to check for pout exceeding array bounds
|
|
C 25-4-91 Fixed bug in cfsq that allowed cfseq to extend off end of arrays
|
|
C 6-6-91 Added a check for pout over end of array to alignd
|
|
C added ctonum to pcon to allow for difference case letters
|
|
C 2-3-92 added filnam = ' ' for som ecalls to openf1
|
|
C ROUTINES IN THIS LIBRARY:
|
|
C DMENU
|
|
C DIAPRW
|
|
C SHOBOX
|
|
C WRITAL
|
|
C ACTOUT
|
|
C DIAPER
|
|
C DIAPRO
|
|
C DIABOX
|
|
C DIAEXP
|
|
C DIAOBS
|
|
C DIALIN
|
|
C ALIGN
|
|
C FMT2
|
|
C FRAME
|
|
C MOVEI
|
|
C HELPD
|
|
C FILEDG
|
|
C QVORH
|
|
C CFSQ
|
|
C ENCONB
|
|
C CFSEQ
|
|
C NCODE
|
|
C SETCON
|
|
C EXTNDM
|
|
C DSTAT
|
|
C REDEFD
|
|
C PRINTD
|
|
C MXPTHD
|
|
C SCORED
|
|
C APBIAS
|
|
C ACALCD
|
|
C GAPRM
|
|
C ALIGND
|
|
SUBROUTINE SWTCHD(SAME,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH)
|
|
INTEGER SAME,OPT
|
|
CHARACTER HELPF*(*)
|
|
CALL YESNO(OPT,'Show main diagonal',
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(OPT.LT.0) RETURN
|
|
SAME = 0
|
|
IF(OPT.EQ.1) SAME = 1
|
|
END
|
|
SUBROUTINE SWTCHI(MARKI,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH)
|
|
INTEGER OPT
|
|
CHARACTER HELPF*(*)
|
|
CALL YESNO(OPT,'Plot identities in matching spans',
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(OPT.LT.0) RETURN
|
|
MARKI = 1
|
|
IF(OPT.EQ.1) MARKI = 0
|
|
END
|
|
C DIAPRW
|
|
SUBROUTINE DIAPRW(S1M,S1P,IDIM1P,S2M,S2P,IDIM2P,
|
|
+M,IDM,LINEB,
|
|
+LINEC,LINEE,LENGTH,MINS,SEQVC,IDIM1,SEQHC,IDIM2,
|
|
+KBOUT,IDEV1,IV1,IH1,ISTARH,ISTARV,SAME)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER S1M(IDIM1P),S1P(IDIM1P)
|
|
INTEGER S2M(IDIM2P),S2P(IDIM2P)
|
|
INTEGER LINEB(IDIM2P),LINEC(IDIM2P),M(IDM,IDM),LINEE(IDIM1P)
|
|
INTEGER TEMPB,TEMPP,SAME
|
|
CHARACTER SEQVC(IDIM1),SEQHC(IDIM2)
|
|
1 CONTINUE
|
|
C WRITE(KBOUT,1000)
|
|
C1000 FORMAT(' List matching spans')
|
|
CALL BUSY(KBOUT)
|
|
LB=LENGTH/2
|
|
IDIM1T=IDIM1P-LENGTH
|
|
IDIM2T=IDIM2P-LENGTH
|
|
C NEED TO SET LINEB TO INITAL VALUE FOR FIRST LENGTH POSITIONS
|
|
C AND LINEE TO LEFT EDGE VALUES
|
|
C FIX LEFT EDGE BY SUMMING ALL THE VALUES FOR A LENGTH LENGTH/2
|
|
C EITHER SIDE OF THE REAL SEQUENCE EDGE
|
|
DO 10 I=1,IDIM1T
|
|
LINEE(I)=0
|
|
IM1 = I - 1
|
|
DO 10 J=1,LENGTH
|
|
K = IM1 + J
|
|
LINEE(I)=LINEE(I)+M(S1M(K),S2M(J))
|
|
10 CONTINUE
|
|
C NOW DO TOP EDGE
|
|
DO 20 I=1,IDIM2T
|
|
C SET UP POINTERS
|
|
LINEC(I)=0
|
|
IM1 = I - 1
|
|
DO 15 J=1,LENGTH
|
|
K = IM1 + J
|
|
LINEC(I)=LINEC(I)+M(S1M(J),S2M(K))
|
|
15 CONTINUE
|
|
20 CONTINUE
|
|
DO 21 I=1,IDIM2T
|
|
LINEC(I)=LINEC(I+1)
|
|
21 CONTINUE
|
|
LINEC(IDIM2T)=0
|
|
C MAIN LOOPS NOW
|
|
C LOOP FOR EACH ROW
|
|
C
|
|
ITV1P=LB+IV1-1-ISTARV+1
|
|
ITV1M=IV1-LB-1-ISTARV+1
|
|
ITH1P=LB+IH1-1-ISTARH+1
|
|
ITH1M=IH1-LB-1-ISTARH+1
|
|
ITV1MM=ITV1M+ISTARV-1
|
|
ITH1MM=ITH1M+ISTARH-1
|
|
DO 200 I=1,IDIM1T
|
|
C SET LINEB TO LINEC, THEN ZERO LINEC
|
|
DO 110 J=1,IDIM2T
|
|
LINEB(J+1)=LINEC(J)
|
|
LINEC(J)=0
|
|
110 CONTINUE
|
|
C SET LINEB(1) TO EDGE VALUE AS ITS OFF PAGE
|
|
LINEB(1)=LINEE(I)
|
|
C
|
|
C NOW COMPARE THIS CHAR OF SEQ1 WITH WHOLE OF SEQ2
|
|
TEMPP=S1P(I)
|
|
TEMPB=S1M(I)
|
|
DO 150 J=1,IDIM2T
|
|
LINEC(J)=LINEB(J)+M(TEMPP,S2P(J))-M(TEMPB,S2M(J))
|
|
IF(LINEC(J).LT.MINS)GO TO 150
|
|
C ABOVE CUTOFF SO WRITE
|
|
IF(I.EQ.J)THEN
|
|
IF(SAME.NE.1)THEN
|
|
WRITE(IDEV1,1006)I+ITV1MM
|
|
WRITE(IDEV1,1005)(SEQVC(K),K=MAX(1,I+ITV1M),MIN(IDIM1,I+ITV1P))
|
|
WRITE(IDEV1,1005)(SEQHC(K),K=MAX(1,J+ITH1M),MIN(IDIM2,J+ITH1P))
|
|
WRITE(IDEV1,1006)J+ITH1MM
|
|
END IF
|
|
ELSE
|
|
WRITE(IDEV1,1006)I+ITV1MM
|
|
WRITE(IDEV1,1005)(SEQVC(K),K=MAX(1,I+ITV1M),MIN(IDIM1,I+ITV1P))
|
|
WRITE(IDEV1,1005)(SEQHC(K),K=MAX(1,J+ITH1M),MIN(IDIM2,J+ITH1P))
|
|
WRITE(IDEV1,1006)J+ITH1MM
|
|
END IF
|
|
1006 FORMAT(' ',I7)
|
|
1005 FORMAT(' ',200A1)
|
|
150 CONTINUE
|
|
200 CONTINUE
|
|
END
|
|
SUBROUTINE SHOBOX(SEQVC,IDIMV,SEQHC,IDIMH,IWX1,IWX2,IWX3,IWX4,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT,KSTARH,KSTARV,IENDH,IENDV)
|
|
C AUTHOR: RODGER STADEN
|
|
CHARACTER HELPF*(*)
|
|
CHARACTER SEQVC(IDIMV),SEQHC(IDIMH)
|
|
IF(0.EQ.0) GO TO 671
|
|
1000 FORMAT(' Zoom-in to matrix',/)
|
|
WRITE(KBOUT,1000)
|
|
1002 FORMAT(' Horizontal sequence')
|
|
601 CONTINUE
|
|
WRITE(KBOUT,1002)
|
|
WRITE(KBOUT,10091)
|
|
10091 FORMAT(' MAX SIZE=36,DEFAULT=XHAIR POSITION - SPAN/2 TO',
|
|
+' XHAIR POSITION + SPAN/2')
|
|
C CALL FSTLST(IWB3,IWB4,KBIN,KBOUT)
|
|
IF(IWB3.LT.0)RETURN
|
|
IF(IWB4.LT.0)RETURN
|
|
IF((IWB3.EQ.0).AND.(IWB4.EQ.0))GO TO 651
|
|
IF((IWB4-IWB3).GT.35)GO TO 601
|
|
IF(IWB3.LT.KSTARH)GO TO 601
|
|
IF(IWB3.GT.IENDH)GO TO 601
|
|
IF(IWB4.GT.IENDH)GO TO 601
|
|
IF(IWB3.GE.IWB4)GO TO 601
|
|
C NE 0 SO USE
|
|
IWX3=IWB3
|
|
IWX4=IWB4
|
|
651 CONTINUE
|
|
IF((IWX4-IWX3).GT.35)GO TO 601
|
|
IF(IWX3.LT.1)GO TO 601
|
|
IF(IWX4.LT.IWX3)GO TO 601
|
|
661 CONTINUE
|
|
WRITE(KBOUT,1001)
|
|
1001 FORMAT(' Vertical sequence')
|
|
WRITE(KBOUT,10091)
|
|
C CALL FSTLST(IWB1,IWB2,KBIN,KBOUT)
|
|
IF((IWB1.EQ.-99).OR.(IWB2.EQ.-99))THEN
|
|
CALL HELP2(IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
GO TO 661
|
|
END IF
|
|
IF(IWB1.LT.0)RETURN
|
|
IF(IWB2.LT.0)RETURN
|
|
IF((IWB1.EQ.0).AND.(IWB2.EQ.0))GO TO 671
|
|
IF((IWB2-IWB1).GT.35)GO TO 661
|
|
IF(IWB1.LT.KSTARV)GO TO 661
|
|
IF(IWB1.GT.IENDV)GO TO 661
|
|
IF(IWB2.GT.IENDV)GO TO 661
|
|
IF(IWB1.GE.IWB2)GO TO 661
|
|
C NE 0 SO USE
|
|
IWX1=IWB1
|
|
IWX2=IWB2
|
|
671 CONTINUE
|
|
C IF((IWX2-IWX1).GT.35)GO TO 661
|
|
C IF(IWX1.LT.1)GOTO 661
|
|
C IF(IWX2.LT.IWX1)GO TO 661
|
|
IF(IWX1.LT.1) RETURN
|
|
IF(IWX2.LT.1) RETURN
|
|
IF(IWX3.LT.1) RETURN
|
|
IF(IWX4.LT.1) RETURN
|
|
IF(IWX1.GT.IDIMV) RETURN
|
|
IF(IWX2.GT.IDIMV) RETURN
|
|
IF(IWX3.GT.IDIMH) RETURN
|
|
IF(IWX4.GT.IDIMH) RETURN
|
|
IF(IWX2-IWX1.GT.35) RETURN
|
|
IF(IWX4-IWX3.GT.35) RETURN
|
|
CALL DIABOX(SEQVC(KSTARV),IDIMV,SEQHC(KSTARH),
|
|
+IDIMH,IWX1,IWX2,IWX3,IWX4,KBOUT)
|
|
RETURN
|
|
END
|
|
SUBROUTINE WRITAL(SEQHC,IDIMH,ISTARH,IENDH,
|
|
+SEQVC,IDIMV,ISTARV,IENDV,IDEVOT,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
C AUTHOR: RODGER STADEN
|
|
CHARACTER HELPF*(*),SEQHC(IDIMH),SEQVC(IDIMV)
|
|
C WRITE(KBOUT,*)' Write out aligned sequences'
|
|
IW3 = ISTARH
|
|
IW4 = IENDH
|
|
CALL GTREG(KBIN,KBOUT,ISTARH,IENDH,IW3,IW4,
|
|
+'Region of horizontal sequence',
|
|
+IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
IW1 = ISTARV
|
|
IW2 = IENDV
|
|
CALL GTREG(KBIN,KBOUT,ISTARV,IENDV,IW1,IW2,
|
|
+'Region of vertical sequence',
|
|
+IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
IDIML=MIN(IW2-IW1,IW4-IW3) + 1
|
|
CALL FMT2(IDEVOT,SEQVC(IW1-ISTARV+1),SEQHC(IW3-ISTARH+1),
|
|
+IDIML,IW1,IW3)
|
|
RETURN
|
|
END
|
|
SUBROUTINE ACTOUT(SEQ,IDIM,ISTART,IEND,IDEVOT,FILOUT,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
CHARACTER FILOUT*(*),HELPF*(*),SEQ(IDIM)
|
|
C AUTHOR: RODGER STADEN
|
|
C FILE SEQUENCE TO DISK
|
|
FILOUT = ' '
|
|
CALL OPENF1(IDEVOT,FILOUT,1,IOK,KBIN,KBOUT,
|
|
+'File name',IHELPS,IHELPE,HELPF,IDEVH)
|
|
IF(IOK.NE.0)RETURN
|
|
IW1 = ISTART
|
|
IW2 = IEND
|
|
C CALL GTREG(KBIN,KBOUT,ISTART,IEND,IW1,IW2,
|
|
C +'Region of sequence',
|
|
C +IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
C IF(IOK.NE.0) RETURN
|
|
IDIML=IW2-IW1+1
|
|
CALL TITOUT(IDEVOT,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH)
|
|
CALL FMTDK(IDEVOT,SEQ(IW1-ISTART+1),IDIML)
|
|
CLOSE(UNIT=IDEVOT)
|
|
RETURN
|
|
END
|
|
C
|
|
C DIAPER
|
|
SUBROUTINE DIAPER(SEQ1,IDIM1,SEQ2,IDIM2,LINEB,LINEC,IDL,MINS,
|
|
+MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX,KBOUT)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER SEQ1(IDIM1),SEQ2(IDIM2)
|
|
INTEGER LINEB(IDL),LINEC(IDL)
|
|
INTEGER TEMPB
|
|
CALL BUSY(KBOUT)
|
|
C WRITE(KBOUT,*)'Working'
|
|
XMAX=IDIM2
|
|
XMIN=1.
|
|
YMAX=IDIM1
|
|
YMIN=1.
|
|
DO 10 I=1,IDIM2+1
|
|
LINEB(I)=0
|
|
LINEC(I)=0
|
|
10 CONTINUE
|
|
C LOOP FOR EACH ROW
|
|
DO 200 I=1,IDIM1
|
|
C ZERO LINE C AND SET LINEB TO LINEC
|
|
DO 110 J=1,IDIM2
|
|
LINEB(J+1)=LINEC(J)
|
|
110 LINEC(J)=0
|
|
C SET LINEB(1) TO ZERO AS ITS OFF PAGE
|
|
LINEB(1)=0
|
|
C
|
|
C NOW COMPARE THIS CHAR OF SEQ1 WITH WHOLE OF SEQ2
|
|
TEMPB=SEQ1(I)
|
|
DO 150 J=1,IDIM2
|
|
IF(TEMPB.NE.SEQ2(J))GO TO 150
|
|
C MATCH SO COUNT LENGTH BY ADDING CURRENT LENGTH OF THIS DIAGONAL
|
|
LINEC(J)=LINEB(J)+1
|
|
C IS THIS HIGH ENOUGH?
|
|
IF(LINEC(J).LT.MINS)GO TO 150
|
|
X=J
|
|
Y=I
|
|
CALL POINT(X,Y,XMAX,XMIN,YMAX,YMIN,
|
|
1MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX)
|
|
150 CONTINUE
|
|
200 CONTINUE
|
|
RETURN
|
|
END
|
|
C
|
|
C DIAPRO
|
|
SUBROUTINE DIAPRO(S1M,S1P,IDIM1P,S2M,S2P,IDIM2P,
|
|
+M,IDM,LINEB,
|
|
+LINEC,LINEE,LENGTH,MINS,
|
|
+MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX,KBOUT,SAME)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER S1M(IDIM1P),S1P(IDIM1P)
|
|
INTEGER S2M(IDIM2P),S2P(IDIM2P)
|
|
INTEGER LINEB(IDIM2P),LINEC(IDIM2P),M(IDM,IDM),LINEE(IDIM1P)
|
|
INTEGER TEMPB,TEMPP,SAME
|
|
CALL BUSY(KBOUT)
|
|
C WRITE(KBOUT,*)'Working'
|
|
IDIM1=IDIM1P-LENGTH
|
|
IDIM2=IDIM2P-LENGTH
|
|
XMAX=IDIM2
|
|
XMIN=1.
|
|
YMAX=IDIM1
|
|
YMIN=1.
|
|
C NEED TO SET LINEB TO INITAL VALUE FOR FIRST LENGTH POSITIONS
|
|
C AND LINEE TO LEFT EDGE VALUES
|
|
C FIX LEFT EDGE BY SUMMING ALL THE VALUES FOR A LENGTH LENGTH/2
|
|
C EITHER SIDE OF THE REAL SEQUENCE EDGE
|
|
DO 10 I=1,IDIM1
|
|
LINEE(I)=0
|
|
IM1 = I - 1
|
|
DO 9 J=1,LENGTH
|
|
K = IM1 + J
|
|
LINEE(I)=LINEE(I)+M(S1M(K),S2M(J))
|
|
9 CONTINUE
|
|
10 CONTINUE
|
|
C NOW DO TOP EDGE
|
|
DO 20 I=1,IDIM2
|
|
C SET UP POINTERS
|
|
LINEC(I)=0
|
|
IM1 = I - 1
|
|
DO 15 J=1,LENGTH
|
|
K = IM1 + J
|
|
LINEC(I)=LINEC(I)+M(S1M(J),S2M(K))
|
|
15 CONTINUE
|
|
20 CONTINUE
|
|
DO 21 I=1,IDIM2
|
|
LINEC(I)=LINEC(I+1)
|
|
21 CONTINUE
|
|
LINEC(IDIM2)=0
|
|
C MAIN LOOPS NOW
|
|
C LOOP FOR EACH ROW
|
|
C
|
|
DO 200 I=1,IDIM1
|
|
C SET LINEB TO LINEC, THEN ZERO LINEC
|
|
DO 110 J=1,IDIM2
|
|
LINEB(J+1)=LINEC(J)
|
|
LINEC(J)=0
|
|
110 CONTINUE
|
|
C SET LINEB(1) TO EDGE VALUE AS ITS OFF PAGE
|
|
LINEB(1)=LINEE(I)
|
|
C
|
|
C NOW COMPARE THIS CHAR OF SEQ1 WITH WHOLE OF SEQ2
|
|
TEMPP=S1P(I)
|
|
TEMPB=S1M(I)
|
|
DO 150 J=1,IDIM2
|
|
LINEC(J)=LINEB(J)+M(TEMPP,S2P(J))-M(TEMPB,S2M(J))
|
|
IF(LINEC(J).LT.MINS)GO TO 150
|
|
C ABOVE CUTOFF SO PLOT
|
|
IF(I.EQ.J)THEN
|
|
IF(SAME.NE.1)THEN
|
|
X=J
|
|
Y=I
|
|
CALL POINT(X,Y,XMAX,XMIN,YMAX,YMIN,
|
|
1 MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX)
|
|
END IF
|
|
ELSE
|
|
X=J
|
|
Y=I
|
|
CALL POINT(X,Y,XMAX,XMIN,YMAX,YMIN,
|
|
1 MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX)
|
|
END IF
|
|
150 CONTINUE
|
|
200 CONTINUE
|
|
END
|
|
C DIAPRI
|
|
SUBROUTINE DIAPRI(S1M,S1P,IDIM1P,S2M,S2P,IDIM2P,
|
|
+M,IDM,LINEB,
|
|
+LINEC,LINEE,LENGTH,MINS,
|
|
+MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX,KBOUT,SAME)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER S1M(IDIM1P),S1P(IDIM1P)
|
|
INTEGER S2M(IDIM2P),S2P(IDIM2P)
|
|
INTEGER LINEB(IDIM2P),LINEC(IDIM2P),M(IDM,IDM),LINEE(IDIM1P)
|
|
INTEGER TEMPB,TEMPP,SAME
|
|
C VERSION TO MARK ALL IDENTITIES
|
|
CALL BUSY(KBOUT)
|
|
IDIM1=IDIM1P-LENGTH
|
|
IDIM2=IDIM2P-LENGTH
|
|
XMAX=IDIM2
|
|
XMIN=1.
|
|
YMAX=IDIM1
|
|
YMIN=1.
|
|
C NEED TO SET LINEB TO INITAL VALUE FOR FIRST LENGTH POSITIONS
|
|
C AND LINEE TO LEFT EDGE VALUES
|
|
C FIX LEFT EDGE BY SUMMING ALL THE VALUES FOR A LENGTH LENGTH/2
|
|
C EITHER SIDE OF THE REAL SEQUENCE EDGE
|
|
DO 10 I=1,IDIM1
|
|
LINEE(I)=0
|
|
IM1 = I - 1
|
|
DO 9 J=1,LENGTH
|
|
K = IM1 + J
|
|
LINEE(I)=LINEE(I)+M(S1M(K),S2M(J))
|
|
9 CONTINUE
|
|
10 CONTINUE
|
|
C NOW DO TOP EDGE
|
|
DO 20 I=1,IDIM2
|
|
C SET UP POINTERS
|
|
LINEC(I)=0
|
|
IM1 = I - 1
|
|
DO 15 J=1,LENGTH
|
|
K = IM1 + J
|
|
LINEC(I)=LINEC(I)+M(S1M(J),S2M(K))
|
|
15 CONTINUE
|
|
20 CONTINUE
|
|
DO 21 I=1,IDIM2
|
|
LINEC(I)=LINEC(I+1)
|
|
21 CONTINUE
|
|
LINEC(IDIM2)=0
|
|
C MAIN LOOPS NOW
|
|
C LOOP FOR EACH ROW
|
|
C
|
|
DO 200 I=1,IDIM1
|
|
C SET LINEB TO LINEC, THEN ZERO LINEC
|
|
DO 110 J=1,IDIM2
|
|
LINEB(J+1)=LINEC(J)
|
|
LINEC(J)=0
|
|
110 CONTINUE
|
|
C SET LINEB(1) TO EDGE VALUE AS ITS OFF PAGE
|
|
LINEB(1)=LINEE(I)
|
|
C
|
|
C NOW COMPARE THIS CHAR OF SEQ1 WITH WHOLE OF SEQ2
|
|
TEMPP=S1P(I)
|
|
TEMPB=S1M(I)
|
|
DO 150 J=1,IDIM2
|
|
LINEC(J)=LINEB(J)+M(TEMPP,S2P(J))-M(TEMPB,S2M(J))
|
|
IF(LINEC(J).LT.MINS)GO TO 150
|
|
C ABOVE CUTOFF SO PLOT
|
|
IF(I.EQ.J)THEN
|
|
IF(SAME.NE.1)THEN
|
|
K1 = I - 1
|
|
K2 = J - 1
|
|
DO 130 K = 1,LENGTH
|
|
K1 = K1 + 1
|
|
K2 = K2 + 1
|
|
IF(S1M(K1).EQ.S2M(K2))THEN
|
|
X = K2
|
|
Y = K1
|
|
CALL POINT(X,Y,XMAX,XMIN,YMAX,YMIN,
|
|
1 MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX)
|
|
END IF
|
|
130 CONTINUE
|
|
END IF
|
|
ELSE
|
|
K1 = I - 1
|
|
K2 = J - 1
|
|
DO 140 K = 1,LENGTH
|
|
K1 = K1 + 1
|
|
K2 = K2 + 1
|
|
IF(S1M(K1).EQ.S2M(K2))THEN
|
|
X = K2
|
|
Y = K1
|
|
CALL POINT(X,Y,XMAX,XMIN,YMAX,YMIN,
|
|
1 MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX)
|
|
END IF
|
|
140 CONTINUE
|
|
END IF
|
|
150 CONTINUE
|
|
200 CONTINUE
|
|
END
|
|
C ALIGN
|
|
C ROUTINE TO INSERT PADDING AS DASHES INTO DIAGON SEQS
|
|
SUBROUTINE ALIGN(SEQVV,IDIMV,SEQHH,IDIMH,
|
|
+MAXSEQ,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,ISTARV,IENDV,ISTARH,IENDH,
|
|
+FILEH,FILEV)
|
|
CHARACTER HELPF*(*),FILEH*(*),FILEV*(*)
|
|
C AUTHOR: RODGER STADEN
|
|
CHARACTER SEQVV(MAXSEQ),SEQHH(MAXSEQ)
|
|
CHARACTER VORH,DASH,IORD
|
|
PARAMETER (MAXPRM = 6)
|
|
CHARACTER PROMPT(2)*(MAXPRM)
|
|
SAVE DASH
|
|
DATA DASH/'-'/
|
|
C ONLY ALLOW EDITING TO SEQUENCES THAT CAN BE TOTALLY CONTAINED IN THE
|
|
C RAM BUFFER, AND WHICH START AT 1 (WE DONT PLAN TO WRITE OUT TO THE
|
|
C DISK BUFFER)
|
|
C IE ISTARH=1, IDIMH<MAXSEQ
|
|
C IE ISTARV=1, IDIMV<MAXSEQ
|
|
C
|
|
1 CONTINUE
|
|
CALL QVORH(VORH,IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT,
|
|
+FILEH,FILEV)
|
|
1004 FORMAT(A)
|
|
IF(VORH.EQ.' ')RETURN
|
|
C CHECK IF OK
|
|
IF(VORH.EQ.'H')THEN
|
|
IF(ISTARH.NE.1)GO TO 100
|
|
IF(IDIMH.GT.MAXSEQ)GO TO 100
|
|
END IF
|
|
IF(VORH.EQ.'V')THEN
|
|
IF(ISTARV.NE.1)GO TO 100
|
|
IF(IDIMV.GT.MAXSEQ)GO TO 100
|
|
END IF
|
|
IF(VORH.EQ.'H') THEN
|
|
MN = ISTARH
|
|
MX = IENDH
|
|
ELSE IF(VORH.EQ.'V') THEN
|
|
MN = ISTARV
|
|
MX = IENDV
|
|
END IF
|
|
MPOS = 1
|
|
CALL GETINT(MN,MX,MPOS,'Position to edit',
|
|
+IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
IPOS = IVAL
|
|
MN = 0
|
|
MX = MIN(IENDH,IENDV)
|
|
MPOS = 0
|
|
CALL GETINT(MN,MX,MPOS,'Number of characters',
|
|
+IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
IF(IVAL.LT.1) RETURN
|
|
NCHAR = IVAL
|
|
I = 1
|
|
PROMPT(1) = 'Insert'
|
|
PROMPT(2) = 'Delete'
|
|
CALL RADION('Select edit command',PROMPT,2,I,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(I.LT.1) RETURN
|
|
IF(I.EQ.2) IORD = 'D'
|
|
IF(I.EQ.1) IORD = 'I'
|
|
IF(IORD.EQ.'D')NCHAR=SIGN(NCHAR,-1)
|
|
C
|
|
C V
|
|
IF(VORH.EQ.'V')THEN
|
|
IF((IDIMV+NCHAR).GT.MAXSEQ)THEN
|
|
WRITE(KBOUT,1008)MAXSEQ
|
|
1008 FORMAT(' This would make sequence longer than ',I6,
|
|
+ ' so not done.')
|
|
GO TO 1
|
|
END IF
|
|
CALL MOVEC(SEQVV,MAXSEQ,IDIMV,IPOS,NCHAR)
|
|
IF(IORD(1:1).EQ.'I')CALL FILLC(SEQVV(IPOS),NCHAR,DASH)
|
|
IDIMV=IDIMV+NCHAR
|
|
IENDV = IDIMV
|
|
WRITE(KBOUT,1010)
|
|
GO TO 1
|
|
END IF
|
|
C H
|
|
IF(VORH.EQ.'H')THEN
|
|
IF((IDIMH+NCHAR).GT.MAXSEQ)THEN
|
|
WRITE(KBOUT,1008)MAXSEQ
|
|
GO TO 1
|
|
END IF
|
|
CALL MOVEC(SEQHH,MAXSEQ,IDIMH,IPOS,NCHAR)
|
|
IF(IORD(1:1).EQ.'I')CALL FILLC(SEQHH(IPOS),NCHAR,DASH)
|
|
IDIMH=IDIMH+NCHAR
|
|
IENDH = IDIMH
|
|
WRITE(KBOUT,1010)
|
|
END IF
|
|
1010 FORMAT(' Edit done')
|
|
GO TO 1
|
|
100 CONTINUE
|
|
C COME HERE IF DISK BUFFERING SEQUENCE AND ITS TOO LONG
|
|
WRITE(KBOUT,1011)
|
|
1011 FORMAT(' Sorry but sequences of this length cannot be edited',
|
|
+' on this machine')
|
|
GO TO 1
|
|
END
|
|
C DIABOX
|
|
C WRITES BOX OF PERFECT MATCHES
|
|
SUBROUTINE DIABOX(SEQ1,IDIM1,SEQ2,IDIM2,I11,I12,I21,I22,
|
|
+KBOUT)
|
|
C AUTHOR: RODGER STADEN
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2),LINE(36),SPACE
|
|
DATA SPACE/' '/
|
|
1000 FORMAT(4X,36A1)
|
|
WRITE(KBOUT,*)
|
|
DO 100 I=I12,I11,-1
|
|
DO 10 J=1,36
|
|
10 LINE(J)=SPACE
|
|
L=0
|
|
DO 50 J=I21,I22
|
|
L=L+1
|
|
50 IF(SEQ1(I).EQ.SEQ2(J))LINE(L)=SEQ1(I)
|
|
WRITE(KBOUT,1001)SEQ1(I),(LINE(K),K=1,L)
|
|
1001 FORMAT(3X,A1,36A1)
|
|
100 CONTINUE
|
|
WRITE(KBOUT,1000)(SEQ2(K),K=I21,I22)
|
|
RETURN
|
|
END
|
|
C DIAEXP COPIED ORIGINALLY FROM ANDREW MC LACHLAN
|
|
C NA=IDIM1
|
|
C NB=IDIM2
|
|
C MATS=MATRIX THE SCORE MATRIX
|
|
C NACID=IDM THE DIMENSION OF THE SCORE MATRIX
|
|
C IS=LENGTH THE SPAN LENGTH
|
|
C LIMIT=HIGHEST INDIVIDUAL SCORE IN SCORE MATRIX MATS
|
|
SUBROUTINE DIAEXP(SEQ1,NA,SEQ2,NB,MATS,NACID,IS1,
|
|
+LIMIT,POLYA,POLYB,IDPOLY,KBIN,KBOUT,IDEV,CHRSET,NORO,
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
CHARACTER HELPF*(*)
|
|
C AUTHOR: RODGER STADEN
|
|
REAL POLYG(50),POLYA(IDPOLY),POLYB(IDPOLY)
|
|
INTEGER JCOMPA(26),JCOMPB(26)
|
|
REAL FCOMPA(26),FCOMPB(26),PR(10)
|
|
INTEGER MATS(NACID,NACID)
|
|
INTEGER NUM(10),OPT
|
|
CHARACTER SEQ1(NA),SEQ2(NB),CHRSET(NACID),NORO,CHART
|
|
INTEGER CTONUM
|
|
EXTERNAL CTONUM
|
|
PARAMETER (SMALL = 1.0E-30)
|
|
PARAMETER (MAXPRM = 29)
|
|
CHARACTER PROMPT(3)*(MAXPRM)
|
|
C DATA SMALL/0.00000000000000000000000000001/
|
|
C IF HAVE ALREADY DONE CALC FOR THIS SPAN SKIP
|
|
IF(NORO.EQ.'O')GO TO 533
|
|
C SET NORO TO OLD IE BEEN DONE BEFORE
|
|
NORO='O'
|
|
CALL BUSY(KBOUT)
|
|
C SET SPAN
|
|
IS=IS1
|
|
C SEQ'S ENTER ALREADY TRANSLATED TO AMINOACID POINTERS 1-26
|
|
C CALC COMPOSITION
|
|
DO 300 IC=1,NACID
|
|
JCOMPA(IC)=0
|
|
JCOMPB(IC)=0
|
|
300 CONTINUE
|
|
DO 301 JA=1,NA
|
|
CHART = SEQ1(JA)
|
|
ITEMP = CTONUM(CHART)
|
|
JCOMPA(ITEMP) = JCOMPA(ITEMP) + 1
|
|
301 CONTINUE
|
|
DO 302 JB=1,NB
|
|
CHART = SEQ2(JB)
|
|
ITEMP = CTONUM(CHART)
|
|
JCOMPB(ITEMP) = JCOMPB(ITEMP) + 1
|
|
302 CONTINUE
|
|
ANA=NA
|
|
ANB=NB
|
|
DO 303 IC=1,NACID
|
|
FCOMPA(IC)=JCOMPA(IC)
|
|
FCOMPA(IC)=FCOMPA(IC)/(ANA+0.0000001)
|
|
FCOMPB(IC)=JCOMPB(IC)
|
|
FCOMPB(IC)=FCOMPB(IC)/(ANB+0.0000001)
|
|
303 CONTINUE
|
|
C WRITE(KBOUT,1004)JCOMPA
|
|
C WRITE(KBOUT,1004)JCOMPB
|
|
1004 FORMAT(2((1X,13I4)/))
|
|
C CALC EXPECTED AV SCORE AND MEAN SQUARE
|
|
S1=0.0
|
|
S2=0.0
|
|
DO 310 L=1,NACID
|
|
DO 310 M=1,NACID
|
|
AM=MATS(L,M)
|
|
S1=S1+FCOMPA(L)*AM*FCOMPB(M)
|
|
S2=S2+FCOMPA(L)*(AM**2)*FCOMPB(M)
|
|
310 CONTINUE
|
|
C RMS SCORE PER POSITION
|
|
IF(S3.GT.0.0)S3=SQRT(S3)
|
|
S3=SQRT(S2-(S1**2))
|
|
S4=IS
|
|
S3=S3*SQRT(S4)
|
|
S1=S1*IS
|
|
C
|
|
WRITE(KBOUT,1000)S1
|
|
1000 FORMAT(' Average score=',F12.5)
|
|
WRITE(KBOUT,1001)S3
|
|
1001 FORMAT(' RMS deviation=',F12.5)
|
|
C
|
|
C SET UP GENERATING POLYNOMIAL WHICH IS THE PROBABILITY OF SELECTION
|
|
C FOR EACH OF THE LEVELS OF SCORE IN THE SCORE MATRIX MATS AND SO
|
|
C DEPENDS ONLY ON THE COMPOSITION OF THE TWO SEQUENCES AND THE
|
|
C DISTRIBUTION OF SCORES IN THE SCORE MATRIX
|
|
DO 320 J=1,50
|
|
POLYG(J)=0.
|
|
320 CONTINUE
|
|
C SUM THE PROBABILITY OF SELECTION FOR EACH LEVEL OF SCORE
|
|
DO 322 L=1,NACID
|
|
DO 321 M=1,NACID
|
|
C SCORE?
|
|
J=MATS(L,M)+1
|
|
POLYG(J)=POLYG(J)+FCOMPA(L)*FCOMPB(M)
|
|
321 CONTINUE
|
|
322 CONTINUE
|
|
C SET THOSE LESS THAN SMALL TO ZERO
|
|
DO 325 J=1,50
|
|
IF(ABS(POLYG(J)).LT.SMALL)POLYG(J)=0.0
|
|
325 CONTINUE
|
|
C
|
|
C COMPUTE PROBABILITIES
|
|
C ZERO POLYA,POLYB AND SET POLYA(1) TO 1.
|
|
DO 400 I=1,IDPOLY
|
|
POLYA(I)=0.
|
|
POLYB(I)=0.
|
|
400 CONTINUE
|
|
POLYA(1)=1.0
|
|
C SET UP
|
|
LORDER=LIMIT
|
|
MORDER=1
|
|
C LOOP SPAN TIMES TO RAISE THE POLYNOMIAL TO THE POWER SPAN
|
|
DO 500 JS=1,IS
|
|
LR=LORDER+1
|
|
MR=MORDER+1
|
|
DO 450 JG=1,LR
|
|
DO 451 JA=1,MR
|
|
C CAN SUM THE JB TERMS BECAUSE THEY GIVE THE SAME SCORE WHEN COMBINED
|
|
JB=JA+JG-1
|
|
POLYB(JB)=POLYB(JB)+POLYA(JA)*POLYG(JG)
|
|
451 CONTINUE
|
|
450 CONTINUE
|
|
MORDER=MORDER+LORDER
|
|
MR=MORDER+1
|
|
C SET VERY SMALL VALUES TO ZERO AND PREPARE FOR NEXT MULTIPLICATION STEP
|
|
DO 453 JA=1,MR
|
|
P=POLYB(JA)
|
|
IF(ABS(P).LT.SMALL)P=0.0
|
|
POLYA(JA)=P
|
|
POLYB(JA)=0.
|
|
453 CONTINUE
|
|
500 CONTINUE
|
|
C CALC CUMULATIVE PROBABILITY
|
|
MR=MORDER+1
|
|
C THIS LOOP LEAVES THE CUMULATIVE PROBABILITY FOR SCORE I IN POLYA(I)
|
|
C IE ALL THE PROBABILITY OF ACHEIVING AT LEAST SCORE I IS STORED IN POLYA(I)
|
|
DO 510 JB=1,MR
|
|
JA=MR+1-JB
|
|
POLYA(JA)=POLYA(JA)+POLYA(JA+1)
|
|
510 CONTINUE
|
|
MORDER=LIMIT*IS
|
|
XMN = 0.0000000001
|
|
MXS = MORDER
|
|
DO 520 J=1,MORDER
|
|
IF(POLYA(J).GT.XMN) GO TO 520
|
|
MXS = J
|
|
GO TO 530
|
|
520 CONTINUE
|
|
530 CONTINUE
|
|
C
|
|
533 CONTINUE
|
|
OPT = 1
|
|
590 CONTINUE
|
|
PROMPT(1) = 'Show probability for a score'
|
|
PROMPT(2) = 'Show score for a probability'
|
|
PROMPT(3) = 'List scores and probabilities'
|
|
CALL RADION('Select probability display mode',PROMPT,3,OPT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(OPT.LT.1) RETURN
|
|
IF(OPT.EQ.2)GO TO 700
|
|
IF(OPT.EQ.3)GO TO 800
|
|
IF(OPT.NE.1)GO TO 590
|
|
MN = 1
|
|
ISCORE = S1 + (MXS - S1) / 2
|
|
CALL GETINT(MN,MXS,ISCORE,'Show probability for score',
|
|
+IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
ISCORE = IVAL
|
|
WRITE(KBOUT,1008)ISCORE,POLYA(ISCORE+1)
|
|
1008 FORMAT(' Probability of score ',I6,' is ',F12.10)
|
|
GO TO 590
|
|
700 CONTINUE
|
|
XMN = 0.0000000001
|
|
XMX = 1.
|
|
PROB = 0.00001
|
|
CALL GETRLS(XMN,XMX,PROB,'Show score for probability',
|
|
+VAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
PROB = VAL
|
|
C LOOK FOR FIRST SCORE GIVING THIS PROBABILITY
|
|
DO 710 J=1,MORDER
|
|
IF(POLYA(J).GT.PROB)GO TO 710
|
|
WRITE(KBOUT,1012)PROB,J
|
|
1012 FORMAT(' Score for probability ',F12.10,' is ',I5)
|
|
GO TO 590
|
|
710 CONTINUE
|
|
GO TO 590
|
|
800 CONTINUE
|
|
MN = 1
|
|
MX = 10
|
|
JSTEP = 5
|
|
CALL GETINT(MN,MX,JSTEP,'Number of steps between scores',
|
|
+IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
JSTEP = IVAL
|
|
MMAX=MORDER+1
|
|
LROW=MMAX/(3*JSTEP)
|
|
DO 541 JROW=1,LROW
|
|
DO 540 J=1,3
|
|
NUM(J)=(JROW-1+(J-1)*LROW)*JSTEP
|
|
PR(J)=POLYA(NUM(J)+1)
|
|
540 CONTINUE
|
|
WRITE(IDEV,1002)(NUM(J),PR(J),J=1,3)
|
|
541 CONTINUE
|
|
1002 FORMAT(3(2X,I5,1X,E12.5))
|
|
GO TO 590
|
|
END
|
|
C
|
|
C DIAOBS
|
|
SUBROUTINE DIAOBS(S1M,S1P,IDIM1P,S2M,S2P,
|
|
+IDIM2P,M,IDM,LINEB,
|
|
+LINEC,LINEE,LENGTH,NORO,PERC,SCORE,IDSCOR,KBIN,KBOUT,
|
|
+IDEV,IHELPS,IHELPE,HELPF,IDEVH,SAME)
|
|
CHARACTER HELPF*(*)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER S1M(IDIM1P),S1P(IDIM1P)
|
|
INTEGER S2M(IDIM2P),S2P(IDIM2P)
|
|
INTEGER LINEB(IDIM2P),LINEC(IDIM2P),M(IDM,IDM),LINEE(IDIM1P)
|
|
INTEGER SCORE(IDSCOR)
|
|
REAL PERC(IDSCOR)
|
|
INTEGER TEMPB,TEMPP,OPT,SAME
|
|
PARAMETER (MAXPRM = 32)
|
|
CHARACTER PROMPT(3)*(MAXPRM)
|
|
CHARACTER NORO
|
|
C ALREADY DONE CALC FOR THIS LENGTH AND REGION?
|
|
IF(NORO.EQ.'O')GO TO 533
|
|
C SET NORO TO OLD IE BEEN DONE BEFORE
|
|
NORO='O'
|
|
CALL BUSY(KBOUT)
|
|
DO 5 I=1,IDSCOR
|
|
SCORE(I)=0
|
|
5 CONTINUE
|
|
IDIM1=IDIM1P-LENGTH
|
|
IDIM2=IDIM2P-LENGTH
|
|
C NEED TO SET LINEB TO INITAL VALUE FOR FIRST LENGTH POSITIONS
|
|
C AND LINEE TO LEFT EDGE VALUES
|
|
C FIX LEFT EDGE BY SUMMING ALL THE VALUES FOR A LENGTH LENGTH/2
|
|
C EITHER SIDE OF THE REAL SEQUENCE EDGE
|
|
DO 10 I=1,IDIM1
|
|
LINEE(I)=0
|
|
IM1 = I - 1
|
|
DO 9 J=1,LENGTH
|
|
K = IM1 + J
|
|
LINEE(I)=LINEE(I)+M(S1M(K),S2M(J))
|
|
9 CONTINUE
|
|
10 CONTINUE
|
|
C NOW DO TOP EDGE
|
|
DO 20 I=1,IDIM2
|
|
C SET UP POINTERS
|
|
LINEC(I)=0
|
|
IM1 = I - 1
|
|
DO 15 J=1,LENGTH
|
|
K = IM1 + J
|
|
LINEC(I)=LINEC(I)+M(S1M(J),S2M(K))
|
|
15 CONTINUE
|
|
20 CONTINUE
|
|
DO 21 I=1,IDIM2
|
|
LINEC(I)=LINEC(I+1)
|
|
21 CONTINUE
|
|
LINEC(IDIM2)=0
|
|
C MAIN LOOPS NOW
|
|
C LOOP FOR EACH ROW
|
|
C
|
|
DO 200 I=1,IDIM1
|
|
C SET LINEB TO LINEC, THEN ZERO LINEC
|
|
DO 110 J=1,IDIM2
|
|
LINEB(J+1)=LINEC(J)
|
|
LINEC(J)=0
|
|
110 CONTINUE
|
|
C SET LINEB(1) TO EDGE VALUE AS ITS OFF PAGE
|
|
LINEB(1)=LINEE(I)
|
|
C
|
|
C NOW COMPARE THIS CHAR OF SEQ1 WITH WHOLE OF SEQ2
|
|
TEMPP=S1P(I)
|
|
TEMPB=S1M(I)
|
|
DO 150 J=1,IDIM2
|
|
LINEC(J)=LINEB(J)+M(TEMPP,S2P(J))-M(TEMPB,S2M(J))
|
|
IF(LINEC(J).NE.0)THEN
|
|
IF(I.EQ.J)THEN
|
|
IF(SAME.NE.1)THEN
|
|
SCORE(LINEC(J)) = SCORE(LINEC(J)) + 1
|
|
END IF
|
|
ELSE
|
|
SCORE(LINEC(J)) = SCORE(LINEC(J)) + 1
|
|
END IF
|
|
END IF
|
|
150 CONTINUE
|
|
200 CONTINUE
|
|
C FIND MAX SCORE
|
|
IDSCOS=IDSCOR+1
|
|
DO 300 I=1,IDSCOR
|
|
J=IDSCOS-I
|
|
IF(SCORE(J).NE.0)THEN
|
|
C FOUND A NON ZERO SCORE. IT MUST BE LARGEST
|
|
MAX=J
|
|
GO TO 301
|
|
END IF
|
|
300 CONTINUE
|
|
301 CONTINUE
|
|
WRITE(KBOUT,1000)MAX
|
|
1000 FORMAT(' Maximum observed score is ',I6)
|
|
C CALC CUMULATIVE SCORE
|
|
DO 400 I=1,MAX-1
|
|
J=MAX-I
|
|
SCORE(J)=SCORE(J)+SCORE(J+1)
|
|
400 CONTINUE
|
|
ISQ=IDIM1*IDIM2
|
|
TOTPER=100./ISQ
|
|
C CALC PERCENT SCORES
|
|
DO 410 I=1,MAX
|
|
PERC(I)=SCORE(I)*TOTPER
|
|
410 CONTINUE
|
|
401 CONTINUE
|
|
533 CONTINUE
|
|
OPT = 1
|
|
590 CONTINUE
|
|
PROMPT(1) = 'Show percentage reaching a score'
|
|
PROMPT(2) = 'Show score for a percentage'
|
|
PROMPT(3) = 'List scores and percentages'
|
|
CALL RADION('Select score display mode',PROMPT,3,OPT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(OPT.LT.1) RETURN
|
|
IF(OPT.EQ.2)GO TO 700
|
|
IF(OPT.EQ.3)GO TO 800
|
|
IF(OPT.NE.1)GO TO 590
|
|
MN = 1
|
|
ISCORE = 0.75 * MAX
|
|
CALL GETINT(MN,MAX,ISCORE,'Show percentage for score',
|
|
+IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
ISCORE = IVAL
|
|
WRITE(KBOUT,1008)ISCORE,PERC(ISCORE)
|
|
1008 FORMAT(' Percentage of points with score ',I6,' is ',F13.9)
|
|
GO TO 590
|
|
700 CONTINUE
|
|
XMN = 0.00001
|
|
XMX = 1.
|
|
PROB = 0.001
|
|
CALL GETRLS(XMN,XMX,PROB,'Show score for percentage',
|
|
+VAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
PROB = VAL
|
|
C LOOK FOR FIRST SCORE GIVING THIS PROBABILITY
|
|
DO 710 J=1,MAX+1
|
|
IF(PERC(J).GT.PROB)GO TO 710
|
|
WRITE(KBOUT,1012)PROB,J
|
|
1012 FORMAT(' Score for percentage ',F13.9,' is ',I5)
|
|
GO TO 590
|
|
710 CONTINUE
|
|
GO TO 590
|
|
800 CONTINUE
|
|
MN = 1
|
|
MX = 10
|
|
JSTEP = 5
|
|
CALL GETINT(MN,MX,JSTEP,'Number of steps between scores',
|
|
+IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
JSTEP = IVAL
|
|
C LIST
|
|
DO 420 I=1,MAX,JSTEP
|
|
IF(PERC(I).NE.100.)WRITE(IDEV,1005)I,SCORE(I),PERC(I)
|
|
1005 FORMAT(1X,I5,1X,I8,1X,E12.5)
|
|
420 CONTINUE
|
|
GO TO 590
|
|
END
|
|
SUBROUTINE DIALIN(ISXMAX,ISYMAX,MARGL,MARGR,MARGB,MARGT)
|
|
C AUTHOR: RODGER STADEN
|
|
XF=MARGL
|
|
XT=XF+MARGR
|
|
YF=MARGB
|
|
YT=YF+MARGT
|
|
XMIN=XF
|
|
XMAX=XT
|
|
YMIN=YF
|
|
YMAX=YT
|
|
CALL LINE(XF,XT,YF,YT,XMAX,XMIN,YMAX,YMIN,
|
|
1MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
CALL VT100M
|
|
RETURN
|
|
END
|
|
C SUBROUTINE TO WRITE 3 LINES OF SEQUENCES
|
|
C
|
|
SUBROUTINE FMT2(IDEV,SEQ1,SEQ2,IDIM,ISW,ISX)
|
|
C AUTHOR: RODGER STADEN
|
|
CHARACTER SEQ1(IDIM),SEQ2(IDIM),MATCH(60)
|
|
INTEGER KL(6)
|
|
C
|
|
C
|
|
ISXX=ISX
|
|
ISWW=ISW
|
|
IE=0
|
|
1 CONTINUE
|
|
C DO 100 I=1,100
|
|
C SET UP DECIMAL COUNTERS
|
|
DO 50 J=1,6
|
|
KL(J)=ISWW
|
|
ISWW=ISWW+10
|
|
50 CONTINUE
|
|
IS=IE+1
|
|
IE=IE+60
|
|
IF(IE.GT.IDIM)IE=IDIM
|
|
C COMPARE THE TWO SECTIONS OF SEQUENCE
|
|
IL=IE-IS+1
|
|
WRITE(IDEV,1001)'V',(KL(K),K=1,1+(IL-1)/10)
|
|
CALL SQMTCH(SEQ1(IS),SEQ2(IS),MATCH,IL)
|
|
WRITE(IDEV,1002)(SEQ1(K),K=IS,IE)
|
|
WRITE(IDEV,1002)(MATCH(K),K=1,IL)
|
|
WRITE(IDEV,1002)(SEQ2(K),K=IS,IE)
|
|
1002 FORMAT( 10X,6(10A1,1X))
|
|
C SET UP DECIMAL COUNTERS
|
|
DO 60 J=1,6
|
|
KL(J)=ISXX
|
|
ISXX=ISXX+10
|
|
60 CONTINUE
|
|
WRITE(IDEV,1001)'H',(KL(K),K=1,1+(IL-1)/10)
|
|
1001 FORMAT( ' ',A,' ',6(I5,6X))
|
|
C1001 FORMAT( 6X,6(I5,6X))
|
|
IF(IE.EQ.IDIM)RETURN
|
|
C100 CONTINUE
|
|
GO TO 1
|
|
END
|
|
C MOVEI
|
|
C MOVES INTEGER ARRAY CONTENTS LEFT OR RIGHT
|
|
SUBROUTINE MOVEI(SEQ,IDIMX,IDIM,IPOS,NCHAR1)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER SEQ(IDIMX)
|
|
INTEGER TO,FROM
|
|
NCHAR=ABS(NCHAR1)
|
|
C LEFT OR RIGHT?
|
|
IF(NCHAR1.LT.0)GO TO 20
|
|
C RIGHT
|
|
FROM=IDIM
|
|
TO=IDIM+NCHAR
|
|
C NUMBER TO MOVE?
|
|
NUM=IDIM-IPOS+1
|
|
DO 10 I=1,NUM
|
|
SEQ(TO)=SEQ(FROM)
|
|
TO=TO-1
|
|
FROM=FROM-1
|
|
10 CONTINUE
|
|
RETURN
|
|
20 CONTINUE
|
|
C LEFT
|
|
FROM=IPOS+NCHAR
|
|
TO=IPOS
|
|
C NUMBER TO MOVE?
|
|
NUM=IDIM-FROM+1
|
|
DO 30 I=1,NUM
|
|
SEQ(TO)=SEQ(FROM)
|
|
TO=TO+1
|
|
FROM=FROM+1
|
|
30 CONTINUE
|
|
RETURN
|
|
END
|
|
SUBROUTINE FILEDG(SEQC,IDIMC,I1,I2,
|
|
+SEQN,IDIMN,SPAN,CHRSET,IDCHR,MSPO2,ISTART,IEND)
|
|
C AUTHOR: RODGER STADEN
|
|
C ROUTINE TO SORT OUT EDGES OF ACTIVE SEQUENCE:
|
|
C WE NOW COPY THE ACTIVE SEQUENCE INTO A TEMPORARY ARRAY
|
|
C FOR PROCESSING (ACTUALLY ALSO CONVERTING CHARS TO NUMBERS)
|
|
C AND WE NEED TO EITHER FILL UP THE EDGES WITH NULL CHARS
|
|
C OR, IF THERE IS REAL DATA AVAILABLE, WITH THAT
|
|
CHARACTER SEQC(IDIMC),CHRSET(IDCHR)
|
|
INTEGER SEQN(IDIMN),SPAN,SPANO2
|
|
INTEGER CTONUM
|
|
EXTERNAL CTONUM
|
|
C DO LEFT EDGE
|
|
C FILL WITH VALUES IDCHR FOR SPAN/2
|
|
C CALC ACTIVE SEQUENCE LENGTH
|
|
IDIMCA=I2-I1+1
|
|
SPANO2=SPAN/2
|
|
DO 1 I=MSPO2-SPANO2,MSPO2
|
|
1 SEQN(I)=IDCHR
|
|
C NOW PUT ANY IN REAL SEQ CHARACTERS AVAILABLE (IF I1>ISTART)
|
|
IF(I1.GT.ISTART)THEN
|
|
C NEED TO INSERT AT MOST SPAN/2 ELEMENTS
|
|
N=MIN(SPANO2+1,I1-ISTART)
|
|
C FIRST TO COPY?
|
|
J1=I1-1-ISTART+1
|
|
K1=MSPO2
|
|
DO 2 I=1,N
|
|
SEQN(K1)=CTONUM(SEQC(J1))
|
|
J1=J1-1
|
|
K1=K1-1
|
|
2 CONTINUE
|
|
END IF
|
|
C DO RIGHT EDGE
|
|
DO 3 I=MSPO2+IDIMCA+1,MSPO2+IDIMCA+SPANO2
|
|
3 SEQN(I)=IDCHR
|
|
IF(I2.LT.IEND)THEN
|
|
N=MIN(SPANO2,(IEND-I2))
|
|
J1=I2+1-ISTART+1
|
|
K1=MSPO2+1+IDIMCA
|
|
DO 4 I=1,N
|
|
SEQN(K1)=CTONUM(SEQC(J1))
|
|
J1=J1+1
|
|
K1=K1+1
|
|
4 CONTINUE
|
|
END IF
|
|
RETURN
|
|
END
|
|
SUBROUTINE QVORH(VHOUT,IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT,
|
|
+FILEH,FILEV)
|
|
C AUTHOR: RODGER STADEN
|
|
CHARACTER HELPF*(*),VHOUT,FILEH*(*),FILEV*(*)
|
|
PARAMETER (MAXPRM = 19)
|
|
CHARACTER PROMPT(2)*(MAXPRM)
|
|
PROMPT(1) = 'Horizontal sequence'
|
|
PROMPT(2) = 'Vertical sequence'
|
|
WRITE(KBOUT,1001)FILEH
|
|
1001 FORMAT(' Horizontal sequence is ',A)
|
|
WRITE(KBOUT,1002)FILEV
|
|
1002 FORMAT(' Vertical sequence is ',A)
|
|
VHOUT = ' '
|
|
I = 1
|
|
CALL RADION('Select sequence',PROMPT,2,I,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(I.LT.1) RETURN
|
|
IF(I.EQ.2) VHOUT = 'V'
|
|
IF(I.EQ.1) VHOUT = 'H'
|
|
END
|
|
SUBROUTINE CFSQ(SEQ1,IDIM1,SEQ2,IDIM2,POSN,WORDP,IDE,IDCHAR,
|
|
+CONSTS,LCONST,LENGTH,MINMAT,
|
|
+MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX,KBOUT,ISAME)
|
|
C 25-4-91 Fixed bug that allowed cfseq to extend off end of arrays
|
|
C by changing call to enconc from idim1 to idim1-minmat+1
|
|
C
|
|
C NOTE !!!!!!! this "fix" does NOT work well for short sequences:
|
|
C enconc only process 1 to idim-length+1 so if idim = idim1-minmat+1
|
|
C <= length-1 nothing gets hashed !!!!!!!!!! FIX IT sometime.
|
|
C should it be idim1-length+1 that gets sent to enconc ??
|
|
C
|
|
INTEGER SEQ1(IDIM1),SEQ2(IDIM2)
|
|
INTEGER POSN(IDIM1),WORDP(IDE),CONSTS(0:LCONST)
|
|
CALL BUSY(KBOUT)
|
|
CALL SETCN(CONSTS,LENGTH,IDCHAR,LCONST)
|
|
CALL ENCONC(SEQ1,IDIM1-MINMAT+1,POSN,WORDP,IDE,IDCHAR,
|
|
+CONSTS,LENGTH,LCONST)
|
|
CALL CFSEQ(SEQ1,IDIM1,POSN,WORDP,IDE,SEQ2,IDIM2,CONSTS,LCONST,
|
|
+LENGTH,IDCHAR,MINMAT,
|
|
+MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX,ISAME)
|
|
END
|
|
SUBROUTINE CFSEQ(SEQ1,IDIM1,POSN,WORDP,IDE,SEQ2,IDIM2,CONSTS,
|
|
+LCONST,
|
|
+LENGTH,IDCHAR,MINMAT,
|
|
+MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX,ISAME)
|
|
INTEGER SEQ1(IDIM1),SEQ2(IDIM2)
|
|
INTEGER POSN(IDIM1),WORDP(IDE),CONSTS(0:LCONST)
|
|
INTEGER NCODEA,EXTNDM
|
|
EXTERNAL NCODEA,EXTNDM
|
|
LEX = MINMAT - LENGTH
|
|
XMAX=IDIM2
|
|
XMIN=1.
|
|
YMAX=IDIM1
|
|
YMIN=1.
|
|
DO 20 I = 1,IDIM2-MINMAT+1
|
|
J = NCODEA(SEQ2(I),LENGTH,CONSTS,IDCHAR,LCONST)
|
|
IF(J.NE.0)THEN
|
|
J1 = WORDP(J)
|
|
IF(J1.NE.0)THEN
|
|
IF (((ISAME.EQ.1).AND.(I.NE.J1)).OR.(ISAME.EQ.0)) THEN
|
|
IMATCH = 0
|
|
IF(MINMAT.GT.LENGTH) IMATCH = EXTNDM(SEQ1(J1+LENGTH),
|
|
+ SEQ2(I+LENGTH),LEX)
|
|
IF(IMATCH.EQ.0)THEN
|
|
X=I
|
|
Y=J1
|
|
CALL POINT(X,Y,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX)
|
|
END IF
|
|
END IF
|
|
10 CONTINUE
|
|
J2 = J1
|
|
J1 = POSN(J2)
|
|
IF(J1.NE.0)THEN
|
|
IF (((ISAME.EQ.1).AND.(I.NE.J1)).OR.(ISAME.EQ.0)) THEN
|
|
IMATCH = 0
|
|
IF(MINMAT.GT.LENGTH) IMATCH = EXTNDM(SEQ1(J1+LENGTH),
|
|
+ SEQ2(I+LENGTH),LEX)
|
|
IF(IMATCH.EQ.0)THEN
|
|
X=I
|
|
Y=J1
|
|
CALL POINT(X,Y,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGXL,MARGXR,MARGYL,MARGYR,ISXMAX,ISYMAX)
|
|
END IF
|
|
END IF
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
END IF
|
|
20 CONTINUE
|
|
END
|
|
INTEGER FUNCTION EXTNDM(SEQ1,SEQ2,LENGTH)
|
|
INTEGER SEQ1(LENGTH),SEQ2(LENGTH)
|
|
EXTNDM = 1
|
|
DO 10 I = 1,LENGTH
|
|
IF(SEQ1(I).NE.SEQ2(I))RETURN
|
|
10 CONTINUE
|
|
EXTNDM = 0
|
|
END
|
|
SUBROUTINE DSTAT(FILEH,ISH1,ISH2,FILEV,ISV1,ISV2,
|
|
+LENGTH,MINPRO,MINPER,KBOUT,ISAME,IMARK)
|
|
CHARACTER FILEH*(*),FILEV*(*)
|
|
WRITE(KBOUT,*)'Horizontal sequence'
|
|
WRITE(KBOUT,1023)FILEH
|
|
WRITE(KBOUT,10010)ISH1,ISH2
|
|
WRITE(KBOUT,*)'Vertical sequence'
|
|
WRITE(KBOUT,1023)FILEV
|
|
1023 FORMAT(' ',A)
|
|
WRITE(KBOUT,10010)ISV1,ISV2
|
|
10010 FORMAT(' Positions',/,' ',I6,' TO ',I6)
|
|
WRITE(KBOUT,10011)LENGTH
|
|
10011 FORMAT(' Span length=',I6)
|
|
WRITE(KBOUT,10014)
|
|
10014 FORMAT(' Scores')
|
|
WRITE(KBOUT,10012)MINPRO
|
|
10012 FORMAT(' Proportional=',I6)
|
|
WRITE(KBOUT,10013)MINPER
|
|
10013 FORMAT(' Identities=',I6)
|
|
IF(IMARK.EQ.1)THEN
|
|
WRITE(KBOUT,1001)
|
|
1001 FORMAT(' Identities on')
|
|
ELSE
|
|
WRITE(KBOUT,1002)
|
|
1002 FORMAT(' Identites off')
|
|
END IF
|
|
IF(ISAME.EQ.1)THEN
|
|
WRITE(KBOUT,1004)
|
|
1004 FORMAT(' Main diagonal blank')
|
|
ELSE
|
|
WRITE(KBOUT,1005)
|
|
1005 FORMAT(' Main diagonal shown')
|
|
END IF
|
|
END
|
|
C REDEFD
|
|
SUBROUTINE REDEFD(IDIMT,J1,J2,MAXSEQ,IDIMA,ISTART,IEND,IDIMB)
|
|
C AUTHOR RODGER STADEN
|
|
C DIAGON IS COMPLICATED BY HAVING SEVERAL ARRARY LIMITS
|
|
C THE SMALL MACHINE VERSION USES A DISK BUFFER
|
|
C ALL VERSIONS HAVE A CHUNK OF SEQUENCE IN RAM (ON THE LARGE MACHINE
|
|
C VERSION THIS WILL BE THE WHOLE SEQUENCE IF < MAXSEQ)
|
|
C ALL VERSIONS HAVE AN ACTIVE REGION FOR COMPARISON
|
|
C FOR ALL PROGRAMS OTHER THAN DIAGON THE ACTIVE REGION IS THE CHUNK IN THE
|
|
C RAM BUFFER, BUT FOR DIAGON IT MAY BE ONLY A PART OF THIS UP TO MAXCOM
|
|
C THIS CHUNK IS ALSO KEPT AS INTEGERS.
|
|
C THIS ROUTINE IS TO INITIALIZE THE ARRAY POINTERS AND SIZES
|
|
C MAXSEQ = THE DIMENSION OF THE RAM BUFFER SEQ
|
|
C IDIMT = THE ACTUAL SEQUENCE LENGTH (AND THEREFORE THE NUMBER OF ELEMENTS
|
|
C IN THE DISK BUFFER)
|
|
C ISTART = THE SEQUENCE NUMBER OF THE CHARACTER OCCUPYING SEQ(1)
|
|
C J1 = THE SEQUENCE NUMBER OF THE FIRST CHARACTER IN THE ACTIVE REGION
|
|
C J2 = THE SEQUENCE NUMBER OF THE LAST CHARACTER IN THE ACTIVE REGION
|
|
C IDIMA = J2-J1+1 I.E. THE NUMBER OF ELEMENTS IN THE ACTIVE REGION
|
|
C IEND = THE SEQUENCE NUMBER OF THE LAST ELEMENT OF SEQ
|
|
C IDIMB = IEND-ISTART+1 I.E. THE NUMBER OF ELEMENTS IN THE RAM BUFFER
|
|
K1 = 0
|
|
K2 = 0
|
|
IF(K1.EQ.0)K1=J1
|
|
K2MAX=MIN((K1+MAXSEQ-1),IDIMT)
|
|
IF(K2.GT.K2MAX)K2=K2MAX
|
|
IF(K2.EQ.0)K2=MIN(J2,K2MAX)
|
|
J1=K1
|
|
J2=K2
|
|
IDIMA=J2-J1+1
|
|
IDIMB=IEND-ISTART+1
|
|
END
|
|
SUBROUTINE MATTIN(MATRIX,IDM,FILNAM,FILEP,CHRSET,
|
|
+KBIN,KBOUT,IDEV,MATMAX,IHELPS,IHELPE,HELPF,IDEVH)
|
|
INTEGER MATRIX(IDM,IDM)
|
|
CHARACTER FILEP*(*),FILNAM*(*),CHRSET(IDM),HELPF*(*)
|
|
PARAMETER (MAXPRM = 20)
|
|
CHARACTER PROMPT(3)*(MAXPRM)
|
|
IN = 1
|
|
PROMPT(1) = 'Identity matrix'
|
|
PROMPT(2) = 'MDM78 matrix'
|
|
PROMPT(3) = 'Personal matrix file'
|
|
10 CONTINUE
|
|
CALL RADION('Select score matrix',PROMPT,3,IN,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(IN.LT.1) RETURN
|
|
IF(IN.EQ.1)THEN
|
|
C SET IDENTITIES
|
|
DO 15 I = 1,IDM
|
|
DO 15 J = 1,IDM
|
|
IF(I.EQ.J)THEN
|
|
MATRIX(I,J) = 1
|
|
ELSE
|
|
MATRIX(I,J) = 0
|
|
END IF
|
|
15 CONTINUE
|
|
ELSE IF(IN.EQ.2)THEN
|
|
C READ IN MDM78
|
|
CALL GETMAT(IDEV,FILEP,MATRIX,IDM,CHRSET,KBOUT,IOK)
|
|
IF(IOK.NE.0)GO TO 100
|
|
ELSE IF(IN.EQ.3)THEN
|
|
FILNAM = ' '
|
|
CALL OPENF1(IDEV,FILNAM,0,IOK,KBIN,KBOUT,
|
|
+ 'Matrix file name',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH)
|
|
IF(IOK.NE.0)GO TO 10
|
|
CALL GETMAT(IDEV,FILNAM,MATRIX,IDM,CHRSET,KBOUT,IOK)
|
|
IF(IOK.NE.0)GO TO 100
|
|
ELSE
|
|
GO TO 10
|
|
END IF
|
|
C NEED TO KNOW LARGEST SCORE IN MATRIX FOR STATS
|
|
MATMAX=0
|
|
DO 20 I=1,IDM
|
|
DO 20 J=1,IDM
|
|
IF(MATMAX.LT.MATRIX(I,J))MATMAX=MATRIX(I,J)
|
|
20 CONTINUE
|
|
MATMAX=MATMAX+1
|
|
WRITE(KBOUT,1004)
|
|
1004 FORMAT(' Remember to reset the score for',/,
|
|
+ ' the proportional algorithm')
|
|
RETURN
|
|
100 CONTINUE
|
|
CALL ERROM(KBOUT,'Error in score matrix file')
|
|
GO TO 10
|
|
END
|
|
SUBROUTINE ALIGNM(SEQ1N,SEQ2N,IDIM1,IDIM2,CC,DD,RR,SS,SOP,
|
|
+NMAX,IV1,IH1,IDIMV,IDIMH,MAXSEQ,WTS,IDM,ISCORE,IG,IH,KBIN,
|
|
+KBOUT,IDEV,SEQ1,SEQ2,SEQ1A,SEQ2A,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KEEP,STACK,MAXSTK,STKREC,IOK)
|
|
INTEGER CC(0:NMAX+1),DD(0:NMAX+1),RR(0:NMAX+1),SS(0:NMAX+1)
|
|
INTEGER WTS(0:IDM,0:IDM),SOP(0:2*NMAX),R1,R2
|
|
INTEGER SEQ1N(NMAX),SEQ2N(NMAX),STKREC,STACK(0:MAXSTK)
|
|
CHARACTER SEQ1(MAXSEQ),SEQ2(MAXSEQ),SEQ1A(NMAX),SEQ2A(NMAX)
|
|
CHARACTER HELPF*(*),PAD
|
|
SAVE PAD
|
|
DATA PAD/','/
|
|
C
|
|
C 12-1-91 modified call to alignd to send nmax instead of maxseq
|
|
C
|
|
C MN = 1
|
|
C MX = NMAX
|
|
C WINDOW = NMAX
|
|
C CALL GETINT(MN,MX,WINDOW,'Window size',
|
|
C +IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
C IF(IOK.NE.0) RETURN
|
|
C WINDOW = IVAL
|
|
C10 CONTINUE
|
|
CALL BUSY(KBOUT)
|
|
LW1 = MIN(IDIM1,NMAX)
|
|
LW2 = MIN(IDIM2,NMAX)
|
|
N = MAX(LW1,LW2)
|
|
CALL DIFF(SEQ1N,SEQ2N,LW1,LW2,
|
|
+CC,DD,RR,SS,SOP,N,
|
|
+WTS,IDM,ISCORE,IG,IH,KBOUT,STACK,MAXSTK,STKREC,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
CALL ALIGND(SEQ1(IV1),SEQ2(IH1),SOP,
|
|
+SEQ1A,SEQ2A,N,NMAX,R1,R2,PAD,NP1,NP2,IDIM1,IDIM2)
|
|
KPOUT = MAX(R1,R2)
|
|
PC = PCON(SEQ1A,SEQ2A,KPOUT,PAD)
|
|
CALL FMT2(IDEV,SEQ1A,SEQ2A,KPOUT,IV1,IH1)
|
|
WRITE(IDEV,1001)PC
|
|
1001 FORMAT(' Conservation ',F5.1,'%')
|
|
WRITE(IDEV,1002)NP1,NP2
|
|
1002 FORMAT(' Number of padding characters inserted',I6,' and',I6)
|
|
CALL YESNO(KEEP,'Keep alignment',
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(KEEP.LT.0) RETURN
|
|
IF(KEEP.EQ.1) RETURN
|
|
NMOV = KPOUT - IDIM1
|
|
IV2 = IV1 + IDIM1
|
|
CALL MOVEC(SEQ1,MAXSEQ,IDIMV,IV2,NMOV)
|
|
IDIMV = IDIMV + NMOV
|
|
CALL SQCOPY(SEQ1A,SEQ1(IV1),KPOUT)
|
|
NMOV = KPOUT - IDIM2
|
|
IV2 = IH1 + IDIM2
|
|
CALL MOVEC(SEQ2,MAXSEQ,IDIMH,IV2,NMOV)
|
|
IDIMH = IDIMH + NMOV
|
|
CALL SQCOPY(SEQ2A,SEQ2(IH1),KPOUT)
|
|
END
|
|
SUBROUTINE ALIGND(SEQ1,SEQ2,S,SEQ1A,SEQ2A,NMAX,MAXSEQ,
|
|
+R1,R2,PAD,NP1,NP2,IDIM1,IDIM2)
|
|
CHARACTER SEQ1(MAXSEQ),SEQ2(MAXSEQ),SEQ1A(MAXSEQ),SEQ2A(MAXSEQ)
|
|
CHARACTER PAD
|
|
INTEGER S(0:2*NMAX),P1,P2,POUT,R1,R2
|
|
C 12-1-91 Added checks for pout going off end of array
|
|
C 6-6-91 Added another check!
|
|
P1 = 1
|
|
P2 = 1
|
|
NP1 = 0
|
|
NP2 = 0
|
|
POUT = 1
|
|
I = -1
|
|
10 CONTINUE
|
|
IF((P1.LE.IDIM1).AND.(P2.LE.IDIM2).AND.(POUT.LE.MAXSEQ)) THEN
|
|
I = I + 1
|
|
IF(S(I).EQ.0) THEN
|
|
C WRITE(*,*)P1,'=',P2,SEQ1(P1),SEQ2(P2)
|
|
SEQ1A(POUT) = SEQ1(P1)
|
|
SEQ2A(POUT) = SEQ2(P2)
|
|
P1 = P1 + 1
|
|
P2 = P2 + 1
|
|
POUT = POUT + 1
|
|
ELSE
|
|
IF(S(I).LT.0) THEN
|
|
K = ABS(S(I))
|
|
C WRITE(*,*)'INSERT ',K,' AT',P2
|
|
CALL SQCOPY(SEQ1(P1),SEQ1A(POUT),K)
|
|
CALL FILLC(SEQ2A(POUT),K,PAD)
|
|
NP2 = NP2 + K
|
|
P1 = P1 + K
|
|
POUT = POUT + K
|
|
ELSE
|
|
K = S(I)
|
|
C WRITE(*,*)'INSERT ',K,' AT',P1
|
|
CALL SQCOPY(SEQ2(P2),SEQ2A(POUT),K)
|
|
CALL FILLC(SEQ1A(POUT),K,PAD)
|
|
NP1 = NP1 + K
|
|
P2 = P2 + K
|
|
POUT = POUT + K
|
|
END IF
|
|
END IF
|
|
GO TO 10
|
|
END IF
|
|
R1 = POUT - 1
|
|
R2 = POUT - 1
|
|
J = 0
|
|
K = 0
|
|
IF((P1.LE.IDIM1).AND.(POUT.LE.MAXSEQ)) THEN
|
|
J = IDIM1 - P1 + 1
|
|
J = MIN(J,MAXSEQ-POUT+1)
|
|
CALL SQCOPY(SEQ1(P1),SEQ1A(POUT),J)
|
|
R1 = R1 + J
|
|
END IF
|
|
IF((P2.LE.IDIM2).AND.(POUT.LE.MAXSEQ)) THEN
|
|
K = IDIM2 - P2 + 1
|
|
K = MIN(K,MAXSEQ-POUT+1)
|
|
CALL SQCOPY(SEQ2(P2),SEQ2A(POUT),K)
|
|
R2 = R2 + K
|
|
END IF
|
|
I = R1 - R2
|
|
IF(I.GT.0) THEN
|
|
CALL FILLC(SEQ2A(POUT+K),I,PAD)
|
|
ELSE IF(I.LT.0) THEN
|
|
I = ABS(I)
|
|
CALL FILLC(SEQ1A(POUT+J),I,PAD)
|
|
END IF
|
|
END
|
|
REAL FUNCTION PCON(SEQ1,SEQ2,L,PAD)
|
|
CHARACTER SEQ1(L),SEQ2(L),PAD
|
|
INTEGER CTONUM
|
|
EXTERNAL CTONUM
|
|
N = 0
|
|
DO 10 I = 1,L
|
|
IF(SEQ1(I).EQ.PAD) SEQ1(I) = '-'
|
|
IF(SEQ2(I).EQ.PAD) SEQ2(I) = '-'
|
|
IF(CTONUM(SEQ1(I)).EQ.CTONUM(SEQ2(I))) N = N + 1
|
|
10 CONTINUE
|
|
PCON = 100. * REAL(N)/REAL(L)
|
|
END
|
|
SUBROUTINE DIFF(SA,SB,IDIM1,IDIM2,CC,DD,RR,SS,SOP,NMAX,
|
|
+W,IDM,SCORE,G,H,KBOUT,STACK,MAXSTK,STKREC,IOK)
|
|
IMPLICIT INTEGER(A-Z)
|
|
INTEGER CC(0:NMAX+1),DD(0:NMAX+1),RR(0:NMAX+1),SS(0:NMAX+1)
|
|
INTEGER W(0:IDM,0:IDM),SOP(0:2*NMAX),STACK(0:MAXSTK)
|
|
INTEGER SA(NMAX),SB(NMAX)
|
|
EXTERNAL GAP
|
|
IOK = 0
|
|
CALL FILLI(SOP,2*NMAX+1,0)
|
|
M = IDIM1
|
|
N = IDIM2
|
|
LM = G + H
|
|
LAST = 0
|
|
SAPP = 0
|
|
A = 1
|
|
B = 1
|
|
TB = G
|
|
TE = H
|
|
MIDI = 0
|
|
MIDJ = 0
|
|
MIDC = 0
|
|
TYPE = 0
|
|
C INITIALISE THE STACK WITH A PUSH
|
|
CALL STACKH(1,A,B,M,N,TB,TE,MIDI,MIDJ,TYPE,MIDC,990,
|
|
+STACK,MAXSTK,STKREC,IOK)
|
|
IF(IOK.NE.0) GO TO 999
|
|
90 CONTINUE
|
|
IF(N.LE.0)THEN
|
|
IF(M.GT.0) CALL DEL(M,LAST,SOP,NMAX,SAPP)
|
|
ANS = GAP(M,G,H)
|
|
GO TO 980
|
|
END IF
|
|
IF(M.LE.1) THEN
|
|
IF(M.LE.0) THEN
|
|
CALL INS(N,LAST,SOP,NMAX,SAPP)
|
|
ANS = GAP(N,G,H)
|
|
GO TO 980
|
|
END IF
|
|
IF(TB.GT.TE) TB = TE
|
|
MIDC = TB + H + GAP(N,G,H)
|
|
MIDJ = 0
|
|
ITP = B - 1
|
|
DO 100 J = 1,N
|
|
C = GAP(J-1,G,H) + W(SA(A),SB(J+ITP)) + GAP(N-J,G,H)
|
|
IF (C.LT.MIDC) THEN
|
|
MIDC = C
|
|
MIDJ = J
|
|
END IF
|
|
100 CONTINUE
|
|
IF(MIDJ.EQ.0)THEN
|
|
CALL INS(N,LAST,SOP,NMAX,SAPP)
|
|
CALL DEL(1,LAST,SOP,NMAX,SAPP)
|
|
ELSE
|
|
IF(MIDJ.GT.1) CALL INS(MIDJ-1,LAST,SOP,NMAX,SAPP)
|
|
CALL REP(LAST,SOP,NMAX,SAPP)
|
|
IF(MIDJ.LT.N) CALL INS(N-MIDJ,LAST,SOP,NMAX,SAPP)
|
|
END IF
|
|
ANS = MIDC
|
|
GO TO 980
|
|
END IF
|
|
MIDI = M/2
|
|
C FORWARD PHASE
|
|
CC(0) = 0
|
|
T = G
|
|
DO 200 J=1,N
|
|
T = T + H
|
|
CC(J) = T
|
|
DD(J) = T + G
|
|
200 CONTINUE
|
|
T = TB
|
|
JTP = B - 1
|
|
DO 400 I = 1,MIDI
|
|
S = CC(0)
|
|
T = T + H
|
|
C = T
|
|
CC(0) = C
|
|
E = T + G
|
|
ITP = I + A - 1
|
|
DO 300 J = 1,N
|
|
C = C + LM
|
|
E = E + H
|
|
IF(C.LT.E) E = C
|
|
C = CC(J) + LM
|
|
D = DD(J) + H
|
|
IF(C.LT.D) D = C
|
|
C = S + W(SA(ITP),SB(J+JTP))
|
|
IF(E.LT.C) C = E
|
|
IF(D.LT.C) C = D
|
|
S = CC(J)
|
|
CC(J) = C
|
|
DD(J) = D
|
|
300 CONTINUE
|
|
400 CONTINUE
|
|
DD(0) = CC(0)
|
|
C REVERSE PHASE
|
|
RR(N) = 0
|
|
T = G
|
|
DO 500 J =N-1,0,-1
|
|
T = T + H
|
|
RR(J) = T
|
|
SS(J) = T + G
|
|
500 CONTINUE
|
|
T = TE
|
|
DO 700 I=M-1,MIDI,-1
|
|
S = RR(N)
|
|
T = T + H
|
|
C = T
|
|
RR(N) = C
|
|
E = T + G
|
|
ITP = I + A
|
|
DO 600 J = N-1,0,-1
|
|
C = C + LM
|
|
E = E + H
|
|
IF(C.LT.E) E = C
|
|
C = RR(J) + LM
|
|
D = SS(J) + H
|
|
IF(C.LT.D) D = C
|
|
C = S + W(SA(ITP),SB(J+B))
|
|
IF(E.LT.C) C = E
|
|
IF(D.LT.C) C = D
|
|
S = RR(J)
|
|
RR(J) = C
|
|
SS(J) = D
|
|
600 CONTINUE
|
|
700 CONTINUE
|
|
SS(N) = RR(N)
|
|
C FIND OPTIMAL MIDPOINT
|
|
MIDC = CC(0) + RR(0)
|
|
MIDJ = 0
|
|
TYPE = 1
|
|
DO 800 J = 0,N
|
|
C = CC(J) + RR(J)
|
|
IF(C.LE.MIDC) THEN
|
|
IF((C.LT.MIDC).OR.(CC(J).NE.DD(J)).AND.(RR(J).EQ.SS(J)))THEN
|
|
MIDC = C
|
|
MIDJ = J
|
|
END IF
|
|
END IF
|
|
800 CONTINUE
|
|
DO 900 J = N,0,-1
|
|
C = DD(J) + SS(J) - G
|
|
IF(C.LT.MIDC) THEN
|
|
MIDC = C
|
|
MIDJ = J
|
|
TYPE = 2
|
|
END IF
|
|
900 CONTINUE
|
|
C CONQUER RECURSIVELY AROUND MIDPOINT
|
|
IF(TYPE.NE.1) GO TO 960
|
|
IF(TYPE.EQ.1) THEN
|
|
CALL STACKH(1,A,B,M,N,TB,TE,MIDI,MIDJ,TYPE,MIDC,950,
|
|
+STACK,MAXSTK,STKREC,IOK)
|
|
IF(IOK.NE.0) GO TO 999
|
|
M = MIDI
|
|
N = MIDJ
|
|
TE = G
|
|
GO TO 90
|
|
END IF
|
|
950 CONTINUE
|
|
CALL STACKH(1,A,B,M,N,TB,TE,MIDI,MIDJ,TYPE,MIDC,980,
|
|
+STACK,MAXSTK,STKREC,IOK)
|
|
IF(IOK.NE.0) GO TO 999
|
|
A = A + MIDI
|
|
B = B + MIDJ
|
|
M = M - MIDI
|
|
N = N - MIDJ
|
|
TB = G
|
|
GO TO 90
|
|
960 CONTINUE
|
|
CALL STACKH(1,A,B,M,N,TB,TE,MIDI,MIDJ,TYPE,MIDC,970,
|
|
+STACK,MAXSTK,STKREC,IOK)
|
|
IF(IOK.NE.0) GO TO 999
|
|
M = MIDI - 1
|
|
N = MIDJ
|
|
TE = 0
|
|
GO TO 90
|
|
970 CONTINUE
|
|
CALL DEL(2,LAST,SOP,NMAX,SAPP)
|
|
CALL STACKH(1,A,B,M,N,TB,TE,MIDI,MIDJ,TYPE,MIDC,980,
|
|
+STACK,MAXSTK,STKREC,IOK)
|
|
IF(IOK.NE.0) GO TO 999
|
|
A = A + MIDI + 1
|
|
B = B + MIDJ
|
|
M = M - MIDI - 1
|
|
N = N - MIDJ
|
|
TB = 0
|
|
GO TO 90
|
|
980 CONTINUE
|
|
SCORE = MIDC
|
|
CALL STACKH(2,A,B,M,N,TB,TE,MIDI,MIDJ,TYPE,MIDC,ADDR,
|
|
+STACK,MAXSTK,STKREC,IOK)
|
|
IF(IOK.NE.0) GO TO 999
|
|
IF(ADDR.EQ.950) GO TO 950
|
|
IF(ADDR.EQ.970) GO TO 970
|
|
IF(ADDR.EQ.980) GO TO 980
|
|
IF(ADDR.EQ.990) GO TO 990
|
|
WRITE(KBOUT,*)'Unexpected address in align'
|
|
IOK = 4
|
|
RETURN
|
|
999 CONTINUE
|
|
IF(IOK.EQ.1) THEN
|
|
WRITE(KBOUT,*)'Stack overflow'
|
|
ELSE IF(IOK.EQ.2) THEN
|
|
WRITE(KBOUT,*)'Stack underflow'
|
|
ELSE IF (IOK.EQ.3) THEN
|
|
WRITE(KBOUT,*)'Unexpected stack task'
|
|
END IF
|
|
990 CONTINUE
|
|
END
|
|
SUBROUTINE GETGAP(KBIN,KBOUT,IG,IH,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,IDM,IOK)
|
|
CHARACTER HELPF*(*)
|
|
C FOR GAP OF LENGTH K, COST IS G + H*K
|
|
C WHERE G IS COST OF STARTING GAP, AND H IS ADDED FOR EACH ELEMENT IN GAP
|
|
C WRITE(KBOUT,*)'THE COST OF A GAP OF LENGTH K = G + H*K'
|
|
MN = 1
|
|
MX = 100
|
|
IG = 10
|
|
IF(IDM.EQ.5) IG = 20
|
|
CALL GETINT(MN,MX,IG,'Penalty for starting a gap',
|
|
+IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
IG = IVAL
|
|
MN = 1
|
|
MX = 100
|
|
IH = 10
|
|
IF(IDM.EQ.5) IH = 5
|
|
CALL GETINT(MN,MX,IH,'Penalty for each residue in gap',
|
|
+IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
IH = IVAL
|
|
END
|
|
INTEGER FUNCTION GAP(K,G,H)
|
|
INTEGER K,G,H
|
|
C NOTE FOR DES: HE HAS EQ HERE
|
|
IF(K.LE.0) THEN
|
|
GAP = 0
|
|
ELSE
|
|
GAP = G + H * K
|
|
END IF
|
|
END
|
|
SUBROUTINE DEL(K,LAST,S,NMAX,SAPP)
|
|
INTEGER S(0:NMAX*2),SAPP
|
|
IF(LAST.LT.0) THEN
|
|
S(SAPP-1) = S(SAPP-1) - K
|
|
LAST = -K
|
|
ELSE
|
|
S(SAPP) = -K
|
|
SAPP = SAPP + 1
|
|
LAST = -K
|
|
END IF
|
|
END
|
|
SUBROUTINE INS(K,LAST,S,NMAX,SAPP)
|
|
INTEGER S(0:NMAX*2),SAPP
|
|
IF (LAST.LT.0) THEN
|
|
S(SAPP-1) = K
|
|
S(SAPP) = LAST
|
|
C DES HAS NEXT LINE. I REMOVED IT 16-5-89
|
|
C LAST = K
|
|
SAPP = SAPP + 1
|
|
ELSE
|
|
S(SAPP) = K
|
|
SAPP = SAPP + 1
|
|
LAST = K
|
|
END IF
|
|
END
|
|
SUBROUTINE REP(LAST,S,NMAX,SAPP)
|
|
INTEGER S(0:NMAX*2),SAPP
|
|
S(SAPP) = 0
|
|
SAPP = SAPP + 1
|
|
LAST = 0
|
|
END
|
|
SUBROUTINE STACKH(JOB,A,B,M,N,TB,TE,MIDI,MIDJ,TYPE,MIDC,ADDR,
|
|
+STACK,MAXSTK,STKREC,IOK)
|
|
IMPLICIT INTEGER (A-Z)
|
|
INTEGER STACK(0:MAXSTK)
|
|
SAVE
|
|
DATA SP/-1/
|
|
C HANDLE STACK: 1 = PUSH, 2 = POP
|
|
IF(JOB.EQ.1) THEN
|
|
SP = SP + STKREC
|
|
IF(SP.GT.MAXSTK) THEN
|
|
C WRITE(*,*)'HELP, STACK OVERFLOW'
|
|
IOK = 1
|
|
RETURN
|
|
END IF
|
|
STACK(SP-10) = A
|
|
STACK(SP-9) = B
|
|
STACK(SP-8) = M
|
|
STACK(SP-7) = N
|
|
STACK(SP-6) = TB
|
|
STACK(SP-5) = TE
|
|
STACK(SP-4) = MIDI
|
|
STACK(SP-3) = MIDJ
|
|
STACK(SP-2) = TYPE
|
|
STACK(SP-1) = MIDC
|
|
STACK(SP) = ADDR
|
|
ELSE IF(JOB.EQ.2) THEN
|
|
SP = SP - STKREC
|
|
IF(SP.LT.-1) THEN
|
|
C WRITE(*,*)'HELP, STACK UNDERFLOW'
|
|
IOK = 2
|
|
RETURN
|
|
END IF
|
|
A = STACK(SP+1)
|
|
B = STACK(SP+2)
|
|
M = STACK(SP+3)
|
|
N = STACK(SP+4)
|
|
TB = STACK(SP+5)
|
|
TE = STACK(SP+6)
|
|
MIDI = STACK(SP+7)
|
|
MIDJ = STACK(SP+8)
|
|
TYPE = STACK(SP+9)
|
|
MIDC = STACK(SP+10)
|
|
ADDR = STACK(SP+11)
|
|
ELSE
|
|
C WRITE(*,*)'HELP, STACK COCKUP'
|
|
IOK = 3
|
|
C STOP
|
|
END IF
|
|
END
|
|
SUBROUTINE PAMDIS(MATRIX,WTS,IDM)
|
|
INTEGER MATRIX(IDM,IDM),WTS(0:IDM,0:IDM)
|
|
IF(IDM.EQ.26) THEN
|
|
N = 8
|
|
M = 0
|
|
DO 10 I = 1,IDM
|
|
DO 5 J = 1,IDM
|
|
K = MATRIX(I,J)
|
|
M = MAX(M,K)
|
|
5 CONTINUE
|
|
10 CONTINUE
|
|
DO 20 I = 1,IDM
|
|
DO 15 J = 1,IDM
|
|
K = M - MATRIX(I,J)
|
|
WTS(I,J) = K
|
|
15 CONTINUE
|
|
20 CONTINUE
|
|
DO 30 I = 0,IDM
|
|
WTS(I,0) = N
|
|
WTS(0,I) = N
|
|
30 CONTINUE
|
|
RETURN
|
|
END IF
|
|
DO 40 I = 1,IDM
|
|
DO 40 J = 1,IDM
|
|
IF(I.EQ.J) THEN
|
|
WTS(I,J) = 0
|
|
ELSE
|
|
WTS(I,J) = 20
|
|
END IF
|
|
40 CONTINUE
|
|
DO 50 I = 0,IDM
|
|
WTS(I,0) = 10
|
|
WTS(0,I) = 10
|
|
50 CONTINUE
|
|
END
|
|
SUBROUTINE ENCONC(SEQ,IDIM,POSN,WORDP,IDE,IDCHAR,CONSTS,LENGTH,
|
|
+LCONST)
|
|
C AUTHOR RODGER STADEN
|
|
INTEGER SEQ(IDIM)
|
|
INTEGER POSN(IDIM),WORDP(IDE),CONSTS(0:LCONST)
|
|
INTEGER NCODEA
|
|
EXTERNAL NCODEA
|
|
C ENCODES A SEQUENCE OF LENGTH IDIM AND CHARACTERSET SIZE IDCHAR
|
|
C INTO TWO ARRAYS: WORDP(I) CONTAINS THE POSITION OF THE FIRST OCCURRENCE
|
|
C OF WORD(I), POSN(I) CONTAINS A LINKED LIST OF SECOND, THIRD,... OCCURENCES
|
|
C OF WORD
|
|
CALL FILLI(POSN,IDIM,0)
|
|
CALL FILLI(WORDP,IDE,0)
|
|
DO 20 I = 1, IDIM-LENGTH+1
|
|
J = NCODEA(SEQ(I),LENGTH,CONSTS,IDCHAR,LCONST)
|
|
IF(J.NE.0)THEN
|
|
J1 = WORDP(J)
|
|
IF(J1.EQ.0)THEN
|
|
WORDP(J) = I
|
|
ELSE
|
|
10 CONTINUE
|
|
J2 = J1
|
|
J1 = POSN(J2)
|
|
IF(J1.NE.0) GO TO 10
|
|
POSN(J2) = I
|
|
END IF
|
|
END IF
|
|
20 CONTINUE
|
|
END
|
|
SUBROUTINE QICKS(SEQ1,IDIM1,POSN,WORDP,IDE,SEQ2,IDIM2,CONSTS,
|
|
+LENGTH,IDM,LCONST,HIST,MAXSEQ,MATRIX,SEQS,MAIND,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX,SPAN,MINPRO,KBOUT,RNSD,
|
|
+ISAME)
|
|
INTEGER SEQ1(IDIM1),SEQ2(IDIM2),SEQS(IDIM1),SPAN
|
|
INTEGER POSN(IDIM1),WORDP(IDE),CONSTS(0:LCONST)
|
|
INTEGER MAIND(IDM)
|
|
INTEGER HIST(-MAXSEQ:MAXSEQ)
|
|
PARAMETER (MAXDIA = 20)
|
|
INTEGER TOPD(MAXDIA),TOPI(MAXDIA),TOPJ(MAXDIA),MATRIX(IDM,IDM)
|
|
EXTERNAL NCODEA
|
|
NDIAG = MAXDIA
|
|
CALL BUSY(KBOUT)
|
|
CALL FILLI(HIST(-IDIM1),IDIM2+IDIM1+1,0)
|
|
CALL SETCN(CONSTS,LENGTH,IDM,LCONST)
|
|
C WRITE(*,*)'CONSTS'
|
|
C WRITE(*,*)CONSTS
|
|
CALL ENCONC(SEQ1,IDIM1,POSN,WORDP,IDE,IDM,CONSTS,LENGTH,
|
|
+LCONST)
|
|
DO 4 I = 1,IDM
|
|
MAIND(I) = MATRIX(I,I)
|
|
4 CONTINUE
|
|
C WRITE(*,*)'MAIND'
|
|
C WRITE(*,*)MAIND
|
|
CALL WDSCR(SEQ1,SEQS,IDIM1,LENGTH,MAIND,IDM)
|
|
C WRITE(*,*)'SEQS'
|
|
C WRITE(*,*)(SEQS(K),K=1,30)
|
|
C WRITE(*,*)'SEQ1'
|
|
C WRITE(*,*)(SEQ1(K),K=1,30)
|
|
C WRITE(*,*)IDIM1,IDIM2,LENGTH
|
|
DO 20 I = 1,IDIM2-LENGTH+1
|
|
J = NCODEA(SEQ2(I),LENGTH,CONSTS,IDM,LCONST)
|
|
IF(J.NE.0)THEN
|
|
J1 = WORDP(J)
|
|
IF(J1.NE.0)THEN
|
|
K = I - J1
|
|
C FOR IDENTITIES ADD 1 ON NEXT LINE (NOT SEQS)
|
|
HIST(K) = HIST(K) + SEQS(J1)
|
|
10 CONTINUE
|
|
J2 = J1
|
|
J1 = POSN(J2)
|
|
IF(J1.NE.0)THEN
|
|
K = I - J1
|
|
C FOR IDENTITIES ADD 1 ON NEXT LINE (NOT SEQS)
|
|
HIST(K) = HIST(K) + SEQS(J1)
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
END IF
|
|
20 CONTINUE
|
|
IF (ISAME.EQ.1) HIST(0) = 0
|
|
CALL MHIST(HIST,IDIM1,IDIM2,TOPD,TOPI,TOPJ,
|
|
+NDIAG,MAXSEQ,RNSD)
|
|
IF(NDIAG.EQ.0) THEN
|
|
WRITE(KBOUT,*)' No diagonals found scoring',RNSD,
|
|
+ ' sd above mean'
|
|
RETURN
|
|
END IF
|
|
ISPO2 = SPAN/2
|
|
XMIN = 1.
|
|
XMAX = IDIM2
|
|
YMIN = 1.
|
|
YMAX = IDIM1
|
|
CALL VECTOM
|
|
DO 40 I = 1,NDIAG
|
|
IF(TOPD(I).NE.0) THEN
|
|
L = MIN(IDIM1-TOPJ(I),IDIM2-TOPI(I)) + 1
|
|
CALL DSCORP(SEQ1(TOPJ(I)),SEQ2(TOPI(I)),L,MATRIX,IDM,
|
|
+ SPAN,MINPRO,ISPO2,TOPJ(I),TOPI(I),
|
|
+XMAX,XMIN,YMAX,YMIN,MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
END IF
|
|
40 CONTINUE
|
|
CALL VT100M
|
|
END
|
|
SUBROUTINE DSCORP(SEQ1,SEQ2,L,MATRIX,IDM,SPAN,MINSCR,SPO2,JS,IS,
|
|
+XMAX,XMIN,YMAX,YMIN,MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
INTEGER SEQ1(L),SEQ2(L),MATRIX(IDM,IDM),SPAN,FRONT,BACK,SPO2
|
|
C 8-6-91 Fixed bug that allowed span>L
|
|
M = 0
|
|
FRONT = SPAN
|
|
BACK = 0
|
|
DO 10 I = 1,MIN(SPAN,L)
|
|
M = M + MATRIX(SEQ1(I),SEQ2(I))
|
|
10 CONTINUE
|
|
IF(M.GE.MINSCR) THEN
|
|
Y = JS + SPO2
|
|
X = IS + SPO2
|
|
C WRITE(*,*)IX,IY
|
|
CALL POINT(X,Y,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
END IF
|
|
DO 20 I = 2,L-SPAN+1
|
|
FRONT = FRONT + 1
|
|
BACK = BACK + 1
|
|
MM = MATRIX(SEQ1(BACK),SEQ2(BACK))
|
|
MP = MATRIX(SEQ1(FRONT),SEQ2(FRONT))
|
|
M = M - MM + MP
|
|
IF(M.GE.MINSCR) THEN
|
|
Y = JS + BACK + SPO2
|
|
X = IS + BACK + SPO2
|
|
C WRITE(*,*)IX,IY
|
|
CALL POINT(X,Y,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
END IF
|
|
20 CONTINUE
|
|
END
|
|
SUBROUTINE MHIST(HIST,IDIM1,IDIM2,TOPD,TOPI,TOPJ,NDIAG,
|
|
+ MAXSEQ,RNSD)
|
|
INTEGER HIST(-MAXSEQ:MAXSEQ),TOPD(NDIAG),TOPI(NDIAG),TOPJ(NDIAG)
|
|
C ANALYSE HISTOGRAM TO FIND BEST NDIAG DIAGONALS
|
|
C LOOK AT THE TOP TEN SCORES, TOP DOWN
|
|
C IF ELEMENT I > TOP(J) THEN MOVE J+1 TO NDIAG-1 DOWN
|
|
C PUT I AT TOP(J)
|
|
C DIAGO FINDS THE INTERSECTION OF DIAGONAL I ON THE 2 AXES
|
|
C FIRST FIND THE TOP SCORES AND THERE HIST ELEMENT NO
|
|
C THEN GET THEIR AXES INTERSECTIONS
|
|
IDIAG = 0
|
|
CALL FILLI(TOPD,NDIAG,0)
|
|
CALL FILLI(TOPI,NDIAG,0)
|
|
CALL FILLI(TOPJ,NDIAG,0)
|
|
RMSQ = 0.
|
|
RM = 0.
|
|
DO 20 I = -IDIM1,IDIM2
|
|
RJ = HIST(I)
|
|
RM = RM + RJ
|
|
RMSQ = RMSQ + RJ * RJ
|
|
20 CONTINUE
|
|
N = IDIM1 + IDIM2
|
|
RM = RM / N
|
|
RMSQ = RMSQ / N
|
|
RM2 = RM * RM
|
|
SD = 0.
|
|
T = RMSQ - RM2
|
|
IF(T.GT.0.) SD = SQRT(T)
|
|
MINS = NINT(RM + RNSD * SD)
|
|
C TRY TO FIND HIGHEST SCORES FIRST
|
|
DO 100 I = 0,IDIM2
|
|
M = HIST(I)
|
|
IF(M.GT.MINS) THEN
|
|
IDIAG = IDIAG + 1
|
|
IF(M.GT.TOPD(NDIAG)) THEN
|
|
DO 50 J = 1,NDIAG
|
|
IF(M.GT.TOPD(J)) THEN
|
|
DO 40 K = NDIAG-1,J,-1
|
|
TOPD(K+1) = TOPD(K)
|
|
TOPI(K+1) = TOPI(K)
|
|
40 CONTINUE
|
|
TOPD(J) = M
|
|
TOPI(J) = I
|
|
GO TO 60
|
|
END IF
|
|
50 CONTINUE
|
|
END IF
|
|
END IF
|
|
60 CONTINUE
|
|
100 CONTINUE
|
|
DO 200 I = -1,-IDIM1,-1
|
|
M = HIST(I)
|
|
IF(M.GT.MINS) THEN
|
|
IDIAG = IDIAG + 1
|
|
IF(M.GT.TOPD(NDIAG)) THEN
|
|
DO 150 J = 1,NDIAG
|
|
IF(M.GT.TOPD(J)) THEN
|
|
DO 140 K = NDIAG-1,J,-1
|
|
TOPD(K+1) = TOPD(K)
|
|
TOPI(K+1) = TOPI(K)
|
|
140 CONTINUE
|
|
TOPD(J) = M
|
|
TOPI(J) = I
|
|
GO TO 160
|
|
END IF
|
|
150 CONTINUE
|
|
END IF
|
|
END IF
|
|
160 CONTINUE
|
|
200 CONTINUE
|
|
NDIAG = MIN(IDIAG,NDIAG)
|
|
DO 300 I = 1,NDIAG
|
|
J = TOPI(I)
|
|
CALL DIAGO(J,TOPI(I),TOPJ(I))
|
|
300 CONTINUE
|
|
END
|
|
SUBROUTINE DIAGO(I,II,JJ)
|
|
IF(I.GE.0) THEN
|
|
JJ = 1
|
|
II = I + 1
|
|
ELSE
|
|
II = 1
|
|
JJ = ABS(I) + 1
|
|
END IF
|
|
END
|
|
SUBROUTINE WDSCR(SEQ1N,SEQS,IDIM1,KTUP,MAIND,IDM)
|
|
INTEGER SEQ1N(IDIM1),SEQS(IDIM1),MAIND(IDM)
|
|
DO 10 I =1,IDIM1-KTUP+1
|
|
K = 0
|
|
DO 5 J = I,I+KTUP-1
|
|
K = K + MAIND(SEQ1N(J))
|
|
5 CONTINUE
|
|
SEQS(I) = K
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE DP21(MATMAX,AVSCOR,LENGTH,MINPRO,
|
|
+KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
CHARACTER HELPF*(*)
|
|
MN = 1
|
|
MX = MAX(1,MATMAX-1) * LENGTH
|
|
MINPRO = INT(AVSCOR*LENGTH)
|
|
CALL GETINT(MN,MX,MINPRO,'Proportional score',
|
|
+ IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.EQ.0) MINPRO = IVAL
|
|
END
|
|
SUBROUTINE DP22(IDM,MINPER,
|
|
+KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
CHARACTER HELPF*(*)
|
|
MN = 1
|
|
IF(IDM.EQ.5) THEN
|
|
MX = 100
|
|
ELSE
|
|
MX = 20
|
|
END IF
|
|
CALL GETINT(MN,MX,MINPER,'Identity score',
|
|
+ IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.EQ.0) MINPER = IVAL
|
|
END
|
|
SUBROUTINE DP20(LENGTH,MXSPAN,IDIMVA,IDIMHA,IDIMVP,IDIMHP,
|
|
+NOROO,NOROE,MAXSEQ,
|
|
+SEQVC,IDIMBV,ISV1,ISV2,SEQV,MXCOMP,
|
|
+CHRSET,IDM,MSPO2,ISTARV,IENDV,
|
|
+SEQHC,IDIMBH,ISH1,ISH2,SEQH,
|
|
+ISTARH,IENDH,LB,LF,
|
|
+KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
CHARACTER HELPF*(*),NOROO,NOROE,CHRSET(IDM)
|
|
CHARACTER SEQHC(MAXSEQ),SEQVC(MAXSEQ)
|
|
INTEGER SEQH(MXCOMP),SEQV(MXCOMP)
|
|
LTEMP=LENGTH
|
|
5 CONTINUE
|
|
MN = 1
|
|
MX = MIN(MXSPAN,MIN(IDIMVA,IDIMHA))
|
|
CALL GETINT(MN,MX,LENGTH,'Odd span length',
|
|
+ IVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
IF(MOD(IVAL,2).NE.1)GO TO 5
|
|
LENGTH = IVAL
|
|
LB=(LENGTH+1)/2
|
|
LF=LENGTH/2
|
|
IDIMVP=IDIMVA+LENGTH
|
|
IDIMHP=IDIMHA+LENGTH
|
|
NOROO='N'
|
|
NOROE='N'
|
|
C IF LENGTH LONGER THAN BEFORE EXTEND INTEGER BUFFERS
|
|
IF(LTEMP.LT.LENGTH)THEN
|
|
CALL FILEDG(SEQVC,IDIMBV,ISV1,ISV2,SEQV,MXCOMP,LENGTH,
|
|
+ CHRSET,IDM,MSPO2,ISTARV,IENDV)
|
|
CALL FILEDG(SEQHC,IDIMBH,ISH1,ISH2,SEQH,MXCOMP,LENGTH,
|
|
+ CHRSET,IDM,MSPO2,ISTARH,IENDH)
|
|
END IF
|
|
END
|
|
SUBROUTINE DP33(RNSD,
|
|
+KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
CHARACTER HELPF*(*)
|
|
RMN = 0.
|
|
RMX = 10.
|
|
CALL GETRL(RMN,RMX,RNSD,'Number of sd above mean',
|
|
+ RVAL,KBIN,KBOUT,IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.EQ.0) RNSD = RVAL
|
|
END
|
|
SUBROUTINE CFSQL(SEQ1,IDIM1,SEQ2,IDIM2,POSN,WORDP,IDE,IDCHAR,
|
|
+CONSTS,LCONST,LENGTH,MINMAT,IDEV,KBOUT,ISAME,SEQ,HIST,MAXSEQ)
|
|
INTEGER SEQ1(IDIM1),SEQ2(IDIM2),HIST(-MAXSEQ:MAXSEQ)
|
|
INTEGER POSN(IDIM1),WORDP(IDE),CONSTS(0:LCONST)
|
|
CHARACTER SEQ(IDIM1)
|
|
CALL BUSY(KBOUT)
|
|
CALL SETCN(CONSTS,LENGTH,IDCHAR,LCONST)
|
|
CALL ENCONC(SEQ1,IDIM1-MINMAT+1,POSN,WORDP,IDE,IDCHAR,
|
|
+CONSTS,LENGTH,LCONST)
|
|
CALL CFSEQL(SEQ1,IDIM1,POSN,WORDP,IDE,SEQ2,IDIM2,CONSTS,LCONST,
|
|
+LENGTH,IDCHAR,MINMAT,IDEV,ISAME,SEQ,HIST,MAXSEQ)
|
|
END
|
|
SUBROUTINE CFSEQL(SEQ1,IDIM1,POSN,WORDP,IDE,SEQ2,IDIM2,CONSTS,
|
|
+LCONST,
|
|
+LENGTH,IDCHAR,MINMAT,IDEV,ISAME,SEQ,HIST,MAXSEQ)
|
|
INTEGER SEQ1(IDIM1),SEQ2(IDIM2)
|
|
INTEGER POSN(IDIM1),WORDP(IDE),CONSTS(0:LCONST)
|
|
INTEGER NCODEA,EXTNDN
|
|
INTEGER HIST(-MAXSEQ:MAXSEQ)
|
|
CHARACTER SEQ(IDIM1)
|
|
EXTERNAL NCODEA,EXTNDN
|
|
LEX = MINMAT - LENGTH - 1
|
|
C
|
|
C if we are looking for internal repeats we switch off half the matrix
|
|
C by making the saved score high
|
|
C
|
|
IF (ISAME.EQ.1) THEN
|
|
CALL FILLI(HIST(-IDIM1),IDIM1,IDIM1+IDIM2)
|
|
CALL FILLI(HIST(0),IDIM2,0)
|
|
ELSE
|
|
CALL FILLI(HIST(-IDIM1),IDIM1+IDIM2+1,0)
|
|
END IF
|
|
DO 20 I = 1,IDIM2-MINMAT+1
|
|
J = NCODEA(SEQ2(I),LENGTH,CONSTS,IDCHAR,LCONST)
|
|
IF(J.NE.0)THEN
|
|
J1 = WORDP(J)
|
|
IF(J1.NE.0)THEN
|
|
IF (((ISAME.EQ.1).AND.(I.NE.J1)).OR.(ISAME.EQ.0)) THEN
|
|
LT = 1 + MIN(IDIM1-(J1+LENGTH),IDIM2-(I+LENGTH))
|
|
IMATCH = EXTNDN(SEQ1(J1+LENGTH),SEQ2(I+LENGTH),LT)
|
|
IF(IMATCH.GT.LEX)THEN
|
|
L = I + IMATCH + LENGTH
|
|
K = I - J1
|
|
C
|
|
C if this match ends furthest away we display it and save it
|
|
C
|
|
IF (HIST(K).LT.L) THEN
|
|
HIST(K) = L
|
|
WRITE(IDEV,1000)I,J1,IMATCH+LENGTH
|
|
1000 FORMAT(' Positions',I7,'h',I7,'v and length',I7)
|
|
WRITE(IDEV,1001)(SEQ(K),K=J1,J1+IMATCH+LENGTH-1)
|
|
1001 FORMAT(' ',50A1)
|
|
END IF
|
|
END IF
|
|
END IF
|
|
10 CONTINUE
|
|
J2 = J1
|
|
J1 = POSN(J2)
|
|
IF(J1.NE.0)THEN
|
|
IF (((ISAME.EQ.1).AND.(I.NE.J1)).OR.(ISAME.EQ.0)) THEN
|
|
LT = 1 + MIN(IDIM1-(J1+LENGTH),IDIM2-(I+LENGTH))
|
|
IMATCH = EXTNDN(SEQ1(J1+LENGTH),SEQ2(I+LENGTH),LT)
|
|
IF(IMATCH.GT.LEX)THEN
|
|
L = I + IMATCH + LENGTH
|
|
K = I - J1
|
|
IF (HIST(K).LT.L) THEN
|
|
HIST(K) = L
|
|
WRITE(IDEV,1000)I,J1,IMATCH+LENGTH
|
|
WRITE(IDEV,1001)(SEQ(K),K=J1,J1+IMATCH+LENGTH-1)
|
|
END IF
|
|
END IF
|
|
END IF
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
END IF
|
|
20 CONTINUE
|
|
END
|
|
INTEGER FUNCTION EXTNDN(SEQ1,SEQ2,LENGTH)
|
|
INTEGER SEQ1(LENGTH),SEQ2(LENGTH)
|
|
DO 10 I = 1,LENGTH
|
|
IF(SEQ1(I).NE.SEQ2(I)) THEN
|
|
EXTNDN = I - 1
|
|
RETURN
|
|
END IF
|
|
10 CONTINUE
|
|
EXTNDN = LENGTH
|
|
END
|