2111 lines
59 KiB
Fortran
2111 lines
59 KiB
Fortran
C NIPF
|
|
SUBROUTINE FMAIN()
|
|
PARAMETER (MAXSEQ=400000,MAXFIL=200,MAXCOM=20000)
|
|
CHARACTER SEQ(MAXSEQ),PAA(5,5,5)
|
|
INTEGER BOTOPT,TOPOPT
|
|
PARAMETER (NAMLEN = 60)
|
|
CHARACTER*(NAMLEN) HELPF,POINTF,FILMAR,FILNAM,NAMES(MAXFIL)
|
|
INTEGER STARTS(MAXFIL),ST(MAXFIL),LE(MAXFIL)
|
|
INTEGER SUMCOD(4096),MTT(64,3)
|
|
REAL MT(4,4),ET(4,4)
|
|
INTEGER COUNTS(MAXCOM),POINTS(MAXCOM),TOTCOU,SUMCS,DIFFER(MAXSEQ)
|
|
REAL WORKR(MAXSEQ)
|
|
PARAMETER (BOTOPT=-10,TOPOPT=70,
|
|
+ MAXMEN=-3,
|
|
+ MAXOPT=39,
|
|
+ MAXDEV=9,
|
|
+ IDM = 5)
|
|
PARAMETER (
|
|
+ HELPF='NIPFHELP',
|
|
+ POINTF='NIPFHPNT',
|
|
+ FILMAR='NIPFMARG')
|
|
INTEGER HELPS(BOTOPT:TOPOPT),HELPE(BOTOPT:TOPOPT),DEVNOS(MAXDEV)
|
|
INTEGER OPT,MARGB(MAXOPT),MARGT(MAXOPT)
|
|
INTEGER BATCH,AVALL,FVREST,ON,OFF,SORT,GRAPH,BCHANG(5,5)
|
|
PARAMETER (AVALL = 1, FVREST = 2, OFF = 1, ON = 2)
|
|
CHARACTER CHRSET(IDM),ACID,CANACD,TRANF
|
|
PARAMETER (MAXPRM = 28)
|
|
CHARACTER PROMPT(2)*(MAXPRM)
|
|
INTEGER AMINOP,FACID(22),FBASE(5),CTONUM,BCHAMG(5,5,3)
|
|
EXTERNAL IMAXA,IMINA,TRANF,VARIM,AMINOP,FINF,CTONUM,RMAXA,RMINA
|
|
EQUIVALENCE (WORKR,DIFFER)
|
|
DATA CHRSET/'T','C','A','G','-'/
|
|
DATA PAA/'F','F','L','L','-','S','S','S','S','S',
|
|
1'Y','Y','*','*','-','C','C','*','W','-',
|
|
1'-','-','-','-','-','L','L','L','L','L',
|
|
1'P','P','P','P','P','H','H','Q','Q','-',
|
|
1'R','R','R','R','R','-','-','-','-','-','I','I','I','M','-',
|
|
1'T','T','T','T','T',
|
|
1'N','N','K','K','-','S','S','R','R','-','-','-','-','-','-',
|
|
1'V','V','V','V','V','A','A','A','A','A','D','D','E','E','-',
|
|
1'G','G','G','G','G',
|
|
1'-','-','-','-','-','-','-','-','-','-',
|
|
1'-','-','-','-','-','-','-','-','-','-',
|
|
1'-','-','-','-','-','-','-','-','-','-'/
|
|
C Initialise help
|
|
CALL INTHLP('nipf', TOPOPT)
|
|
C GET DEVICE NUMBERS
|
|
CALL UNITNO(KBIN,KBOUT,DEVNOS,MAXDEV)
|
|
CALL OPENGR(DEVNOS(3))
|
|
IGORT = 0
|
|
MTCALL = 0
|
|
WRITE(KBOUT,1000)
|
|
1000 FORMAT(
|
|
+' NIPF v1.0 Copyright Rodger Staden'/)
|
|
C READ IN THE POINTERS TO THE HELP FILE
|
|
IDEVH = DEVNOS(4)
|
|
CALL SETHLP(HELPS,HELPE,BOTOPT,TOPOPT,POINTF,IDEVH,KBOUT)
|
|
CALL INITGR(KBIN,KBOUT,HELPS(0),HELPE(0),HELPF,IDEVH)
|
|
IOK=0
|
|
CALL INITLU(IDM)
|
|
C GET SCREEN AND MARGIN SIZES
|
|
CALL GETMRG(ISXMAX,ISYMAX,MARGL,MARGR,MARGB,MARGT,
|
|
+MAXOPT,DEVNOS(1),FILMAR)
|
|
IDEV=KBOUT
|
|
MOPT=0
|
|
IDIMT = 0
|
|
CALL EMT(MTT)
|
|
C
|
|
C
|
|
1 CONTINUE
|
|
C
|
|
C
|
|
C IF(IOK.NE.0)GO TO 9999
|
|
C give menu, get option
|
|
C
|
|
IDEV=KBOUT
|
|
LINLEN=60
|
|
PROMPT(1) = 'Use file of file names'
|
|
PROMPT(2) = 'Type in a pair of file names'
|
|
IB = 1
|
|
CALL RADION('Select input mode',PROMPT,2,IB,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(IB.LT.1) GO TO 999
|
|
IF(IB.EQ.1) THEN
|
|
CALL GETSQS(DEVNOS(1),DEVNOS(2),SEQ,MAXSEQ,
|
|
+ STARTS,NAMES,
|
|
+ FILNAM,NFILE,MAXFIL,MAXNAM,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH)
|
|
SORT = ON
|
|
BATCH = FVREST
|
|
ELSE
|
|
C
|
|
C Read in the two sequences
|
|
C
|
|
CALL TGETF(DEVNOS(1),FILNAM,KBIN,KBOUT,
|
|
+ SEQ,MAXSEQ,STARTS,NAMES,MAXFIL,MAXNAM,NFILE,IOK,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH)
|
|
SORT = OFF
|
|
BATCH = FVREST
|
|
END IF
|
|
WRITE(KBOUT,*)'Working with ',NFILE,' files'
|
|
C
|
|
C SET DEFAULTS FOR ACTIVE REGION
|
|
C
|
|
C ST(I) contains the active region start for sequence I relative to its start
|
|
C LE(I) contains the active region end of sequence I
|
|
C STARTS(I) contains the start for sequence I in seq
|
|
C
|
|
DO 2 I=1,NFILE
|
|
ST(I) = 1
|
|
LE(I) = STARTS(I+1) - STARTS(I)
|
|
2 CONTINUE
|
|
LONGS = IMAXA(ST,NFILE)
|
|
LONGE = IMINA(LE,NFILE)
|
|
I1=1
|
|
N1=IDIM1
|
|
I2=1
|
|
N2=IDIM2
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
IDIML=MIN(IDIM3,IDIM4)
|
|
BATCH = FVREST
|
|
100 CONTINUE
|
|
CALL BPAUSE(KBIN,KBOUT,IOK)
|
|
CALL MENU(OPT,KOPT,MOPT,MAXOPT,MAXMEN,KBIN,KBOUT,
|
|
+HELPS(0),HELPE(0),HELPF,DEVNOS(4))
|
|
C IF(0.EQ.0) GO TO 99
|
|
C
|
|
C stop
|
|
C
|
|
IF(OPT.EQ.2) GO TO 999
|
|
C
|
|
C type text
|
|
C
|
|
IF(OPT.EQ.6) THEN
|
|
CALL TTEXT(DEVNOS(1),FILNAM,KBIN,KBOUT,
|
|
+ HELPS(OPT),HELPE(OPT),HELPF,DEVNOS(4))
|
|
GO TO 100
|
|
END IF
|
|
C
|
|
C help
|
|
C
|
|
IF(OPT.EQ.1) THEN
|
|
CALL HELP(HELPS,HELPE,BOTOPT,TOPOPT,HELPF,DEVNOS(4),
|
|
+ KBIN,KBOUT)
|
|
GO TO 100
|
|
END IF
|
|
C
|
|
C clear all
|
|
C
|
|
IF(OPT.EQ.10) THEN
|
|
CALL CLEARG
|
|
GO TO 100
|
|
END IF
|
|
C
|
|
C ruler
|
|
C
|
|
IF(OPT.EQ.12)THEN
|
|
CALL RULER(LONGS,LONGE,MARGL,MARGR,
|
|
+ MARGB(OPT),MARGT(OPT),ISXMAX,ISYMAX,KBIN,KBOUT,1,
|
|
+ HELPS(OPT),HELPE(OPT),HELPF,DEVNOS(4))
|
|
GO TO 100
|
|
END IF
|
|
C
|
|
C xhairs
|
|
C
|
|
IF(OPT.EQ.13) THEN
|
|
XMAX = LONGE
|
|
XMIN = LONGS
|
|
YMAX=ISYMAX
|
|
YMIN=0.
|
|
IIIIX=0
|
|
IIIIY=0
|
|
CALL CLEARV
|
|
CALL XHAIRN(XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB(OPT),MARGT(OPT),
|
|
+ ISXMAX,ISYMAX,IIIIX,IIIIY,N,KBOUT,
|
|
+ SEQ(LONGS),LONGS,LONGE-LONGS+1,
|
|
+ SEQ(LONGS),LONGS,LONGE-LONGS+1,1)
|
|
GO TO 100
|
|
END IF
|
|
C
|
|
C clear vt100
|
|
C
|
|
IF(OPT.EQ.11)THEN
|
|
CALL CLEARV
|
|
GO TO 100
|
|
END IF
|
|
C
|
|
C CHANGE MARGINS
|
|
C
|
|
IF(OPT.EQ.14)THEN
|
|
CALL MARGC(ISXMAX,ISYMAX,MARGL,MARGR,MARGB,MARGT,
|
|
+ HELPS(OPT),HELPE(OPT),MAXOPT,HELPF,IDEVH,KBIN,KBOUT)
|
|
GO TO 100
|
|
END IF
|
|
C
|
|
C WRITE LABELS
|
|
C
|
|
IF(OPT.EQ.15) THEN
|
|
CALL LABLER(KBIN,KBOUT,ISXMAX,ISYMAX,
|
|
+ HELPS(OPT),HELPE(OPT),HELPF,DEVNOS(4))
|
|
GO TO 100
|
|
END IF
|
|
C
|
|
C GET DISK OUTPUT FILE ON UNIT DEVNOS(2) IF REQUIRED
|
|
C
|
|
IF(OPT.EQ.7)THEN
|
|
CALL REDIR(IDEV,DEVNOS(2),DEVNOS(3),IGORT,FILNAM,KBIN,KBOUT,
|
|
+ HELPS(OPT),HELPE(OPT),HELPF,IDEVH,KOPT)
|
|
GO TO 100
|
|
END IF
|
|
IF(OPT.EQ.5)THEN
|
|
LINLEN=60
|
|
CALL LISTN(SEQ,MAXSEQ,STARTS,NFILE,LINLEN,LONGS,LONGE,
|
|
+ IDEV,KBOUT,NAMES,MAXFIL)
|
|
GO TO 100
|
|
END IF
|
|
IF(OPT.EQ.32)THEN
|
|
LINLEN=60
|
|
CALL LISTND(SEQ,MAXSEQ,STARTS,NFILE,LINLEN,LONGS,LONGE,
|
|
+ IDEV,KBOUT,NAMES,MAXFIL)
|
|
GO TO 100
|
|
END IF
|
|
IF(OPT.EQ.9)THEN
|
|
LINLEN=60
|
|
CALL LISTNT(SEQ,MAXSEQ,STARTS,NFILE,LINLEN,LONGS,LONGE,
|
|
+ IDEV,KBOUT,NAMES,MAXFIL,PAA)
|
|
GO TO 100
|
|
END IF
|
|
IF(OPT.EQ.33)THEN
|
|
LINLEN=60
|
|
CALL LISTNV(SEQ,MAXSEQ,STARTS,NFILE,LINLEN,LONGS,LONGE,
|
|
+ IDEV,KBOUT,NAMES,MAXFIL,PAA)
|
|
GO TO 100
|
|
END IF
|
|
IF ((OPT.EQ.37).AND.(NFILE.EQ.2)) THEN
|
|
LINLEN=60
|
|
CALL TRAN(SEQ(STARTS(1)),MAXSEQ,ST(1),LE(1),
|
|
+ SEQ(STARTS(2)),MAXSEQ,ST(2),LE(2),
|
|
+ KBIN,KBOUT,IDEV,LINLEN,PAA)
|
|
GO TO 100
|
|
END IF
|
|
IF(OPT.EQ.4)THEN
|
|
ISTART = 1
|
|
IEND = STARTS(2) - 1
|
|
DO 7 I=1,NFILE
|
|
CALL GTREG(KBIN,KBOUT,ISTART,IEND,ST(I),LE(I),
|
|
+ 'Define active region',
|
|
+ HELPS(OPT),HELPE(OPT),HELPF,IDEVH,IOK)
|
|
ST(I+1) = ST(I)
|
|
LE(I+1) = LE(I)
|
|
IF(IOK.NE.0) GO TO 100
|
|
7 CONTINUE
|
|
C
|
|
C Get max values
|
|
C
|
|
LONGS = IMAXA(ST,NFILE)
|
|
LONGE = IMINA(LE,NFILE)
|
|
GO TO 100
|
|
END IF
|
|
IF(OPT.EQ.17) THEN
|
|
PROMPT(1) = 'All versus all'
|
|
PROMPT(2) = 'First versus rest'
|
|
IB = BATCH
|
|
CALL RADION('Select comparison mode',PROMPT,2,IB,
|
|
+ HELPS(OPT),HELPE(OPT),HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(IB.LT.1) GO TO 100
|
|
BATCH = IB
|
|
GO TO 100
|
|
END IF
|
|
IF(OPT.EQ.18) THEN
|
|
PROMPT(ON) = 'Sort on'
|
|
PROMPT(OFF) = 'Sort off'
|
|
IB = ON
|
|
CALL RADION('Select sort mode',PROMPT,2,IB,
|
|
+ HELPS(OPT),HELPE(OPT),HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(IB.LT.1) GO TO 100
|
|
SORT = IB
|
|
GO TO 100
|
|
END IF
|
|
C
|
|
C start of batch processing
|
|
C
|
|
ICOUNT = 0
|
|
TOTCOU = 0
|
|
SUMCS = 0
|
|
C
|
|
C need to allow different batch comparison modes
|
|
C
|
|
IF(BATCH.EQ.AVALL) THEN
|
|
KFILES = 1
|
|
KFILEE = NFILE-1
|
|
NEEDS = (NFILES-1)*NFILES/2
|
|
IF (NEEDS.GT.MAXCOM) THEN
|
|
WRITE(*,*)'Sorry unsufficient space for this comparison'
|
|
GO TO 100
|
|
END IF
|
|
ELSE IF(BATCH.EQ.FVREST) THEN
|
|
KFILES = 1
|
|
KFILEE = KFILES
|
|
ELSE
|
|
WRITE(*,*)'Batch mode not set!!'
|
|
GO TO 100
|
|
END IF
|
|
IF((OPT.EQ.24).OR.(OPT.EQ.35).OR.(OPT.EQ.36).OR.(OPT.EQ.39)) THEN
|
|
DO 111 I=1,5
|
|
DO 111 J=1,5
|
|
BCHANG(I,J) = 0
|
|
DO 111 K = 1,3
|
|
BCHAMG(I,J,K) = 0
|
|
111 CONTINUE
|
|
END IF
|
|
IF(OPT.EQ.38) THEN
|
|
CALL INITMT(MT,SUMCOD,MTCALL)
|
|
MTCALL = MTCALL + 1
|
|
END IF
|
|
GRAPH = -1
|
|
IF((OPT.EQ.25).OR.(OPT.EQ.26).OR.(OPT.EQ.27).OR.(OPT.EQ.31)
|
|
+.OR.(OPT.EQ.34)) THEN
|
|
PROMPT(1) = 'List results'
|
|
PROMPT(2) = 'Plot results'
|
|
GRAPH = ON
|
|
CALL RADION('Select output mode',PROMPT,2,GRAPH,
|
|
+ HELPS(OPT),HELPE(OPT),HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(GRAPH.LT.1) GO TO 100
|
|
CALL FILLI(DIFFER,LONGE-LONGS+1,0)
|
|
END IF
|
|
IF(OPT.EQ.31) THEN
|
|
IPOSM = (LONGE-LONGS+1)/3
|
|
CALL FILLR(WORKR,IPOSM,0.0)
|
|
IW = 0
|
|
DO 600 IPOS = LONGS,LONGS-1+IPOSM*3,3
|
|
CALL FILLI(FACID,22,0)
|
|
DO 500 IFILE = KFILES,KFILEE
|
|
CANACD = TRANF(SEQ(STARTS(IFILE)+IPOS-1),PAA)
|
|
DO 450 JFILE = IFILE+1,NFILE
|
|
ACID = TRANF(SEQ(STARTS(JFILE)+IPOS-1),PAA)
|
|
IF(CANACD.NE.ACID)
|
|
+ FACID(AMINOP(ACID)) = FACID(AMINOP(ACID)) + 1
|
|
450 CONTINUE
|
|
500 CONTINUE
|
|
IW = IW + 1
|
|
WORKR(IW) = VARIM(FACID)
|
|
600 CONTINUE
|
|
IF(GRAPH.EQ.ON) THEN
|
|
XMIN = 1.
|
|
XMAX = REAL(LONGE - LONGS + 1)/3.
|
|
YMIN = 0.
|
|
YMAX = RMAXA(WORKR,IPOSM)
|
|
CALL GETRL(0.,9999.,YMAX,'Maximum Y value for plotting',
|
|
+ VAL,KBIN,KBOUT,HELPS(OPT),HELPE(OPT),HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 100
|
|
YMAX = VAL
|
|
CALL PLTR(WORKR,IPOSM,
|
|
+ XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB(OPT),MARGT(OPT),ISXMAX,ISYMAX)
|
|
GO TO 100
|
|
END IF
|
|
IF(GRAPH.EQ.OFF) THEN
|
|
WRITE(IDEV,6667)
|
|
+ (LONGS+3*(KKKK-1),WORKR(KKKK),KKKK=1,IPOSM)
|
|
6667 FORMAT(I6,F10.3)
|
|
GO TO 100
|
|
END IF
|
|
END IF
|
|
IF(OPT.EQ.34) THEN
|
|
IPOSM = LONGE-LONGS+1
|
|
CALL FILLR(WORKR,IPOSM,0.0)
|
|
IW = 0
|
|
DO 700 IPOS = LONGS,LONGE
|
|
CALL FILLI(FBASE,5,0)
|
|
DO 650 JFILE = 1,NFILE
|
|
C WRITE(*,*)JFILE
|
|
I = CTONUM(SEQ(STARTS(JFILE)+IPOS-1))
|
|
FBASE(I) = FBASE(I) + 1
|
|
650 CONTINUE
|
|
IW = IW + 1
|
|
WORKR(IW) = FINF(FBASE)
|
|
700 CONTINUE
|
|
IF(GRAPH.EQ.ON) THEN
|
|
XMIN = 1.
|
|
XMAX = REAL(LONGE - LONGS + 1)
|
|
YMIN = RMINA(WORKR,IPOSM)
|
|
YMAX = RMAXA(WORKR,IPOSM)
|
|
CALL GETRL(-9999.,9999.,YMAX,'Maximum Y value for plotting',
|
|
+ VAL,KBIN,KBOUT,HELPS(OPT),HELPE(OPT),HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 100
|
|
YMAX = VAL
|
|
CALL PLTR(WORKR,IPOSM,
|
|
+ XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB(OPT),MARGT(OPT),ISXMAX,ISYMAX)
|
|
GO TO 100
|
|
END IF
|
|
IF(GRAPH.EQ.OFF) THEN
|
|
WRITE(IDEV,6667)
|
|
+ (LONGS+KKKK-1,WORKR(KKKK),KKKK=1,IPOSM)
|
|
GO TO 100
|
|
END IF
|
|
END IF
|
|
DO 300 IFILE = KFILES,KFILEE
|
|
DO 200 JFILE = IFILE+1,NFILE
|
|
IF(OPT.EQ.21)THEN
|
|
C
|
|
C MISMAT NUMBER OF BASE DIFFERENCES
|
|
C
|
|
CALL BSCHN1(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ MISMAT,IDEV)
|
|
END IF
|
|
IF(OPT.EQ.22)THEN
|
|
C
|
|
C CODON CHANGES
|
|
C
|
|
CALL COCHAN(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ ICONSA,ISAME,ICONSC,MISMAT,IDEV,SUMCOD)
|
|
C CONSERVED ACIDS ICONSA=KSUM
|
|
C CONSERVATIVE CHANGES ICONSC=SUMC
|
|
C UNCHANGED CODONS ISAME=ISUM
|
|
SUMCS = SUMCS + ICONSC
|
|
END IF
|
|
IF(OPT.EQ.23)THEN
|
|
C NUMBER OF GENETIC EVENTS
|
|
CALL BSCHN2(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ MISMAT,IDEV)
|
|
END IF
|
|
IF(OPT.EQ.24)THEN
|
|
C MISMAT NUMBER OF BASE DIFFERENCES
|
|
CALL BSCHN3(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),BCHANG)
|
|
END IF
|
|
IF(OPT.EQ.35)THEN
|
|
C MISMAT NUMBER OF BASE DIFFERENCES
|
|
CALL BSCHM3(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),BCHAMG)
|
|
END IF
|
|
IF(OPT.EQ.3000)THEN
|
|
C MISMAT NUMBER OF BASE DIFFERENCES
|
|
CALL BSCHN5(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ MISMAT,IDEV,PAA)
|
|
END IF
|
|
IF(OPT.EQ.28)THEN
|
|
C MISMAT NUMBER OF expressed BASE DIFFERENCES
|
|
CALL BSCHN7(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ MISMAT,IDEV,PAA)
|
|
END IF
|
|
IF(OPT.EQ.29)THEN
|
|
C MISMAT NUMBER OF silent BASE DIFFERENCES
|
|
CALL BSCHN8(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ MISMAT,IDEV,PAA)
|
|
END IF
|
|
IF(OPT.EQ.30)THEN
|
|
C MISMAT NUMBER OF Amino acid DIFFERENCES
|
|
CALL BSCHN9(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ MISMAT,IDEV,PAA)
|
|
END IF
|
|
C
|
|
C GET DIFFERENCES IN DIFFER
|
|
C
|
|
IF(OPT.EQ.25)THEN
|
|
CALL BSCHN4(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ DIFFER)
|
|
END IF
|
|
C
|
|
C GET DIFFERENCES IN DIFFER (accumulated counts of bases that change the
|
|
C amino acid)
|
|
C
|
|
IF(OPT.EQ.26)THEN
|
|
CALL COCHN1(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ DIFFER,PAA)
|
|
END IF
|
|
C
|
|
C estimate mutation rate from observed silent changes
|
|
C
|
|
IF(OPT.EQ.38)THEN
|
|
CALL COCHN4(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ SUMCOD)
|
|
END IF
|
|
C
|
|
C GET DIFFERENCES IN DIFFER (accumulated counts of bases that dont change the
|
|
C amino acid)
|
|
C
|
|
IF(OPT.EQ.27)THEN
|
|
CALL COCHN2(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ DIFFER,PAA)
|
|
END IF
|
|
IF(OPT.EQ.36)THEN
|
|
CALL COCHN3(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ BCHANG,PAA)
|
|
END IF
|
|
IF(OPT.EQ.39)THEN
|
|
CALL COCHN5(
|
|
+ SEQ(STARTS(IFILE)),STARTS(IFILE+1)-STARTS(IFILE),
|
|
+ ST(IFILE),LE(IFILE),
|
|
+ SEQ(STARTS(JFILE)),STARTS(JFILE+1)-STARTS(JFILE),
|
|
+ ST(JFILE),LE(JFILE),
|
|
+ BCHANG,PAA)
|
|
END IF
|
|
ICOUNT = ICOUNT + 1
|
|
COUNTS(ICOUNT) = MISMAT
|
|
POINTS(ICOUNT) = ICOUNT
|
|
TOTCOU = TOTCOU + MISMAT
|
|
200 CONTINUE
|
|
300 CONTINUE
|
|
IF(OPT.EQ.38) THEN
|
|
C CALL UPDET(SUMCOD,PAA,MTT)
|
|
CALL UPDE(MT,ET,MTT,SUMCOD,PAA)
|
|
WRITE(IDEV,1010)(CHRSET(K),K=1,4)
|
|
DO 323 I=1,4
|
|
WRITE(IDEV,1012)CHRSET(I),(MT(I,K),K=1,4)
|
|
323 CONTINUE
|
|
1012 FORMAT(' ',A1,' ',4(F5.4,1X))
|
|
GO TO 100
|
|
END IF
|
|
IF(GRAPH.EQ.ON) THEN
|
|
XMIN = 1.
|
|
XMAX = LONGE - LONGS + 1
|
|
YMIN = 0.
|
|
YMAX = IMAXA(DIFFER,LONGE-LONGS+1)
|
|
IYMAX = YMAX
|
|
CALL GETINT(0,9999,IYMAX,'Maximum Y value for plotting',
|
|
+ IVAL,KBIN,KBOUT,HELPS(OPT),HELPE(OPT),HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 100
|
|
YMAX = IVAL
|
|
CALL PLTI(DIFFER,LONGE-LONGS+1,
|
|
+ XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB(OPT),MARGT(OPT),ISXMAX,ISYMAX)
|
|
GO TO 100
|
|
END IF
|
|
IF(GRAPH.EQ.OFF) THEN
|
|
WRITE(IDEV,6666)
|
|
+ (KKKK+LONGS-1,DIFFER(KKKK),KKKK=1,LONGE-LONGS+1)
|
|
6666 FORMAT(2I6)
|
|
C CLOSE(UNIT=DEVNOS(1))
|
|
GO TO 100
|
|
END IF
|
|
IF(GRAPH.EQ.ON) THEN
|
|
XMIN = 1.
|
|
XMAX = LONGE - LONGS + 1
|
|
YMIN = 0.
|
|
YMAX = IMAXA(DIFFER,LONGE-LONGS+1)
|
|
IYMAX = YMAX
|
|
CALL GETINT(0,9999,IYMAX,'Maximum Y value for plotting',
|
|
+ IVAL,KBIN,KBOUT,HELPS(OPT),HELPE(OPT),HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 100
|
|
YMAX = IVAL
|
|
CALL PLTI(DIFFER,LONGE-LONGS+1,
|
|
+ XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB(OPT),MARGT(OPT),ISXMAX,ISYMAX)
|
|
GO TO 100
|
|
END IF
|
|
IF(GRAPH.EQ.OFF) THEN
|
|
WRITE(IDEV,6666)
|
|
+ (KKKK+LONGS-1,DIFFER(KKKK),KKKK=1,LONGE-LONGS+1)
|
|
GO TO 100
|
|
END IF
|
|
IF(OPT.EQ.24) THEN
|
|
WRITE(IDEV,1010)CHRSET
|
|
1010 FORMAT(' ',5(2X,A1,3X))
|
|
DO 320 I=1,5
|
|
WRITE(IDEV,1011)CHRSET(I),(BCHANG(I,K),K=1,5)
|
|
1011 FORMAT(' ',A1,5I6)
|
|
320 CONTINUE
|
|
CALL KIMURA(BCHANG,DIST)
|
|
WRITE(*,*)'Kimura evolutionary distance=',DIST
|
|
GO TO 100
|
|
END IF
|
|
IF(OPT.EQ.35) THEN
|
|
DO 333 IP=1,3
|
|
WRITE(IDEV,1010)CHRSET
|
|
DO 321 I=1,5
|
|
WRITE(IDEV,1011)CHRSET(I),(BCHAMG(I,K,IP),K=1,5)
|
|
321 CONTINUE
|
|
333 CONTINUE
|
|
GO TO 100
|
|
END IF
|
|
IF(OPT.EQ.36) THEN
|
|
WRITE(IDEV,1010)CHRSET
|
|
DO 322 I=1,5
|
|
WRITE(IDEV,1011)CHRSET(I),(BCHANG(I,K),K=1,5)
|
|
322 CONTINUE
|
|
GO TO 100
|
|
END IF
|
|
IF(OPT.EQ.39) THEN
|
|
WRITE(IDEV,1010)CHRSET
|
|
DO 324 I=1,5
|
|
WRITE(IDEV,1011)CHRSET(I),(BCHANG(I,K),K=1,5)
|
|
324 CONTINUE
|
|
GO TO 100
|
|
END IF
|
|
IF(SORT.EQ.OFF) GO TO 100
|
|
CALL BUB2AS(COUNTS,POINTS,ICOUNT)
|
|
WRITE(IDEV,*)' Sorted list'
|
|
DO 400 II=1,ICOUNT
|
|
N = II
|
|
CALL IANDJ(NFILE,POINTS(N),I,J)
|
|
WRITE(IDEV,1015)NAMES(I)(1:20),NAMES(J)(1:20),COUNTS(II)
|
|
1015 FORMAT(' ',A,' ',A,' distance',I6)
|
|
400 CONTINUE
|
|
WRITE(IDEV,1020)TOTCOU
|
|
1020 FORMAT(' Total distance ',I6)
|
|
IF(OPT.EQ.22)WRITE(IDEV,1021)SUMCS
|
|
1021 FORMAT(' Total conservative changes',I6)
|
|
GO TO 100
|
|
999 CONTINUE
|
|
CALL SHUTD
|
|
END
|
|
REAL FUNCTION VARIM(FACID)
|
|
INTEGER FACID(22),IMAXA
|
|
EXTERNAL IMAXA
|
|
C Variability defined (as kabat):
|
|
C
|
|
C Number of different acids / Freq of most common acid
|
|
C
|
|
NACIDS = 0
|
|
TACIDS = 0.
|
|
DO 10 K = 1,22
|
|
TACIDS = TACIDS + FACID(K)
|
|
IF(FACID(K).NE.0) NACIDS = NACIDS + 1
|
|
10 CONTINUE
|
|
VARIM = 0.
|
|
IF(TACIDS.GT.0.) THEN
|
|
FREQC = REAL(IMAXA(FACID,22)) / TACIDS
|
|
VARIM = REAL(NACIDS) / FREQC
|
|
END IF
|
|
END
|
|
INTEGER FUNCTION AMINOP(ACID)
|
|
CHARACTER ALLONE*22,ACID
|
|
SAVE ALLONE
|
|
DATA ALLONE/'ACDEFGHIKLMNPQRSTVWY*-'/
|
|
AMINOP = 22
|
|
DO 10 I=1,21
|
|
IF(ACID.EQ.ALLONE(I:I)) THEN
|
|
AMINOP = I
|
|
RETURN
|
|
END IF
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE TGETF(IDEV,FILNAM,KBIN,KBOUT,
|
|
+SEQ1,MAXSEQ,STARTS,NAMES,MAXFIL,MAXNAM,NFILE,IOK,
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
CHARACTER SEQ1(MAXSEQ),FILNAM*(*),HELPF*(*)
|
|
CHARACTER*(*) NAMES(MAXFIL)
|
|
INTEGER STARTS(MAXFIL)
|
|
C
|
|
C Routine to read in two sequences
|
|
C
|
|
FILNAM = ' '
|
|
CALL OPENF1(IDEV,FILNAM,0,IOK,KBIN,KBOUT,
|
|
+'File name of first sequence',
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
IF(IOK.NE.0)RETURN
|
|
IDIM1 = MAXSEQ
|
|
CALL ARRFIL(IDEV,SEQ1,IDIM1,KBOUT)
|
|
CLOSE(UNIT=IDEV)
|
|
NAMES(1) = FILNAM
|
|
STARTS(1) = 1
|
|
FILNAM = ' '
|
|
CALL OPENF1(IDEV,FILNAM,0,IOK,KBIN,KBOUT,
|
|
+'File name of second sequence',
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
IF(IOK.NE.0)RETURN
|
|
STARTS(2) = IDIM1 + 1
|
|
IDIM1 = MAXSEQ - IDIM1
|
|
CALL ARRFIL(IDEV,SEQ1(STARTS(2)),IDIM1,KBOUT)
|
|
CLOSE(UNIT=IDEV)
|
|
NAMES(2) = FILNAM
|
|
STARTS(3) = STARTS(2) + IDIM1
|
|
NFILE = 2
|
|
END
|
|
SUBROUTINE GETSQS(IDEV1,IDEV2,SEQ,MAXSEQ,STARTS,NAMES,
|
|
+FILNAM,NFILE,MAXFIL,MAXNAM,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
CHARACTER SEQ(MAXSEQ),HELPF*(*)
|
|
CHARACTER*(*) FILNAM,NAMES(MAXFIL)
|
|
INTEGER STARTS(MAXFIL)
|
|
NFILE = 0
|
|
IEND = 0
|
|
FILNAM = ' '
|
|
CALL OPENF1(IDEV1,FILNAM,0,IOK,KBIN,KBOUT,
|
|
+'File of file names',
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
2 CONTINUE
|
|
IF (NFILE.LT.MAXFIL) THEN
|
|
READ(IDEV1,1000,END=3)FILNAM
|
|
1000 FORMAT(A)
|
|
NFILE = NFILE + 1
|
|
STARTS(NFILE) = IEND + 1
|
|
NAMES(NFILE) = FILNAM
|
|
IDIM2 = MAXSEQ - IEND
|
|
IF (IDIM2.GT.0) THEN
|
|
CALL OPENRS(IDEV2,FILNAM,IOK,LRECL,2)
|
|
CALL ARRFIL(IDEV2,SEQ(STARTS(NFILE)),IDIM2,KBOUT)
|
|
CLOSE(UNIT=IDEV2)
|
|
IEND = IEND + IDIM2
|
|
GO TO 2
|
|
END IF
|
|
C
|
|
C if we get here weve run out of array space
|
|
C
|
|
NFILE = NFILE - 1
|
|
WRITE(KBOUT,*)'Only ',NFILE,' files read: no more memory'
|
|
CLOSE(UNIT=IDEV1)
|
|
RETURN
|
|
ELSE
|
|
C
|
|
C if we get here weve got too many files
|
|
C
|
|
NFILE = NFILE - 1
|
|
CLOSE(UNIT=IDEV1)
|
|
WRITE(KBOUT,*)'Maximum files ',NFILE,' read'
|
|
RETURN
|
|
END IF
|
|
3 CONTINUE
|
|
STARTS(NFILE+1) = IEND + 1
|
|
CLOSE(UNIT=IDEV1)
|
|
END
|
|
SUBROUTINE LISTN(SEQ,MAXCHR,LENSEQ,NFILE,LINLEN,I1,I2,
|
|
+IDEV,KBOUT,NAMES,MAXFIL)
|
|
CHARACTER*(*) NAMES(MAXFIL)
|
|
PARAMETER (MAXLIN = 120)
|
|
CHARACTER SEQ(MAXCHR),LINE*(MAXLIN)
|
|
INTEGER LENSEQ(MAXFIL)
|
|
EXTERNAL NOTIRL
|
|
C SET WIDTH FOR LAST PAGE
|
|
LPAGE=MOD(I2-I1+1,LINLEN)
|
|
C HOW MANY PAGE WIDTHS?
|
|
NPAGE=1+(I2-I1+1)/LINLEN
|
|
IF(MOD(I2-I1+1,LINLEN).EQ.0)THEN
|
|
NPAGE=NPAGE-1
|
|
LPAGE=LINLEN
|
|
END IF
|
|
ISTART=I1-LINLEN
|
|
DO 50 I=1,NPAGE
|
|
ISTART=ISTART+LINLEN
|
|
IF(I.EQ.NPAGE)LINLEN=LPAGE
|
|
WRITE(IDEV,1006)(K,K=ISTART+9,ISTART+LINLEN-1,10)
|
|
DO 40 J=1,NFILE
|
|
KF = ISTART + LENSEQ(J) - 1
|
|
KT = MIN(KF+LINLEN,LENSEQ(J+1)) - 1
|
|
LINE(1:) = ' '
|
|
WRITE(LINE,1003,ERR=60)(SEQ(K),K=KF,KT)
|
|
1003 FORMAT(' ',60A1)
|
|
LINE(KT-KF+4:) = NAMES(J)(1:15)
|
|
WRITE(IDEV,1004,ERR=60)LINE(1:NOTIRL(LINE,MAXLIN,' '))
|
|
1004 FORMAT(A)
|
|
1006 FORMAT(' ',10I10)
|
|
40 CONTINUE
|
|
WRITE(IDEV,1008)
|
|
1008 FORMAT( )
|
|
50 CONTINUE
|
|
RETURN
|
|
60 CONTINUE
|
|
WRITE(KBOUT,*)' Error writing file'
|
|
END
|
|
SUBROUTINE LISTND(SEQ,MAXCHR,LENSEQ,NFILE,LINLEN,I1,I2,
|
|
+IDEV,KBOUT,NAMES,MAXFIL)
|
|
CHARACTER*(*) NAMES(MAXFIL)
|
|
PARAMETER (MAXLIN = 120)
|
|
CHARACTER SEQ(MAXCHR),LINE*(MAXLIN)
|
|
INTEGER LENSEQ(MAXFIL),CTONUM
|
|
EXTERNAL NOTIRL,CTONUM
|
|
C SET WIDTH FOR LAST PAGE
|
|
LPAGE=MOD(I2-I1+1,LINLEN)
|
|
C HOW MANY PAGE WIDTHS?
|
|
NPAGE=1+(I2-I1+1)/LINLEN
|
|
IF(MOD(I2-I1+1,LINLEN).EQ.0)THEN
|
|
NPAGE=NPAGE-1
|
|
LPAGE=LINLEN
|
|
END IF
|
|
ISTART=I1-LINLEN
|
|
DO 50 I=1,NPAGE
|
|
ISTART=ISTART+LINLEN
|
|
IF(I.EQ.NPAGE)LINLEN=LPAGE
|
|
WRITE(IDEV,1006)(K,K=ISTART+9,ISTART+LINLEN-1,10)
|
|
KF = ISTART + LENSEQ(1) - 1
|
|
KT = MIN(KF+LINLEN,LENSEQ(2)) - 1
|
|
WRITE(LINE,1003,ERR=60)(SEQ(K),K=KF,KT)
|
|
1003 FORMAT(' ',60A1)
|
|
LINE(KT-KF+4:) = NAMES(1)(1:15)
|
|
WRITE(IDEV,1004,ERR=60)LINE(1:NOTIRL(LINE,MAXLIN,' '))
|
|
DO 40 J=2,NFILE
|
|
KF = ISTART + LENSEQ(J) - 1
|
|
KT = MIN(KF+LINLEN,LENSEQ(J+1)) - 1
|
|
KSF = ISTART + LENSEQ(1) - 2
|
|
LINE(1:) = ' '
|
|
KKL = 1
|
|
DO 35 KK = KF,KT
|
|
KKL = KKL + 1
|
|
KSF = KSF + 1
|
|
IF(CTONUM(SEQ(KK)).EQ.CTONUM(SEQ(KSF))) THEN
|
|
LINE(KKL:KKL) = '-'
|
|
ELSE
|
|
LINE(KKL:KKL) = SEQ(KK)
|
|
END IF
|
|
35 CONTINUE
|
|
LINE(KT-KF+4:) = NAMES(J)(1:15)
|
|
WRITE(IDEV,1004,ERR=60)LINE(1:NOTIRL(LINE,MAXLIN,' '))
|
|
1004 FORMAT(A)
|
|
1006 FORMAT(' ',10I10)
|
|
40 CONTINUE
|
|
WRITE(IDEV,1008)
|
|
1008 FORMAT( )
|
|
50 CONTINUE
|
|
RETURN
|
|
60 CONTINUE
|
|
WRITE(KBOUT,*)' Error writing file'
|
|
END
|
|
SUBROUTINE LISTNT(SEQ,MAXCHR,LENSEQ,NFILE,LINLIN,I1,I2,
|
|
+IDEV,KBOUT,NAMES,MAXFIL,PAA)
|
|
CHARACTER*(*) NAMES(MAXFIL)
|
|
PARAMETER (MAXLIN = 120)
|
|
CHARACTER SEQ(MAXCHR),LINE*(MAXLIN),PAA*(*),TRANF
|
|
INTEGER LENSEQ(MAXFIL)
|
|
EXTERNAL NOTIRL,TRANF
|
|
C
|
|
C line length is linlen acids which needs 3*linlen bases
|
|
C
|
|
LINLEN = LINLIN
|
|
I2MI1 = I2 - I1 + 1
|
|
C number of whole codons
|
|
I2MI1 = I2MI1/3
|
|
C SET WIDTH FOR LAST PAGE
|
|
LPAGE=MOD(I2MI1,LINLEN)
|
|
C HOW MANY PAGE WIDTHS?
|
|
NPAGE=1+I2MI1/LINLEN
|
|
IF(MOD(I2MI1,LINLEN).EQ.0)THEN
|
|
NPAGE=NPAGE-1
|
|
LPAGE=LINLEN
|
|
END IF
|
|
ISTART=I1-3*LINLEN
|
|
DO 50 I=1,NPAGE
|
|
ISTART=ISTART+3*LINLEN
|
|
IF(I.EQ.NPAGE)LINLEN=LPAGE
|
|
WRITE(IDEV,1006)(K,K=ISTART+29,ISTART+3*LINLEN-1,30)
|
|
DO 40 J=1,NFILE
|
|
KF = ISTART + LENSEQ(J) - 1
|
|
KT = MIN(KF+3*LINLEN,LENSEQ(J+1)) - 1
|
|
LINE(1:) = ' '
|
|
LL = KF - 3
|
|
DO 30 L = 1,LINLEN
|
|
LL = LL + 3
|
|
LINE(L:L) = TRANF(SEQ(LL),PAA)
|
|
30 CONTINUE
|
|
LINE(1+(KT-KF+4)/3:) = NAMES(J)(1:15)
|
|
WRITE(IDEV,1004,ERR=60)LINE(1:NOTIRL(LINE,MAXLIN,' '))
|
|
1004 FORMAT(' ',A)
|
|
1006 FORMAT(' ',10I10)
|
|
40 CONTINUE
|
|
WRITE(IDEV,1008)
|
|
1008 FORMAT( )
|
|
50 CONTINUE
|
|
RETURN
|
|
60 CONTINUE
|
|
WRITE(KBOUT,*)' Error writing file'
|
|
END
|
|
SUBROUTINE LISTNV(SEQ,MAXCHR,LENSEQ,NFILE,LINLIN,I1,I2,
|
|
+IDEV,KBOUT,NAMES,MAXFIL,PAA)
|
|
CHARACTER*(*) NAMES(MAXFIL)
|
|
PARAMETER (MAXLIN = 120)
|
|
CHARACTER SEQ(MAXCHR),LINE*(MAXLIN),PAA*(*),TRANF
|
|
INTEGER LENSEQ(MAXFIL)
|
|
EXTERNAL NOTIRL,TRANF
|
|
C
|
|
C line length is linlen acids which needs 3*linlen bases
|
|
C
|
|
LINLEN = LINLIN
|
|
I2MI1 = I2 - I1 + 1
|
|
C number of whole codons
|
|
I2MI1 = I2MI1/3
|
|
C SET WIDTH FOR LAST PAGE
|
|
LPAGE=MOD(I2MI1,LINLEN)
|
|
C HOW MANY PAGE WIDTHS?
|
|
NPAGE=1+I2MI1/LINLEN
|
|
IF(MOD(I2MI1,LINLEN).EQ.0)THEN
|
|
NPAGE=NPAGE-1
|
|
LPAGE=LINLEN
|
|
END IF
|
|
ISTART=I1-3*LINLEN
|
|
DO 50 I=1,NPAGE
|
|
ISTART=ISTART+3*LINLEN
|
|
IF(I.EQ.NPAGE)LINLEN=LPAGE
|
|
WRITE(IDEV,1006)(K,K=ISTART+29,ISTART+3*LINLEN-1,30)
|
|
KF = ISTART + LENSEQ(1) - 1
|
|
KT = MIN(KF+3*LINLEN,LENSEQ(2)) - 1
|
|
LINE(1:) = ' '
|
|
LL = KF - 3
|
|
DO 20 L = 1,LINLEN
|
|
LL = LL + 3
|
|
LINE(L:L) = TRANF(SEQ(LL),PAA)
|
|
20 CONTINUE
|
|
LINE(1+(KT-KF+4)/3:) = NAMES(1)(1:15)
|
|
WRITE(IDEV,1004,ERR=60)LINE(1:NOTIRL(LINE,MAXLIN,' '))
|
|
DO 40 J=2,NFILE
|
|
KF = ISTART + LENSEQ(J) - 1
|
|
KT = MIN(KF+3*LINLEN,LENSEQ(J+1)) - 1
|
|
KFS = ISTART + LENSEQ(1) - 1
|
|
LINE(1:) = ' '
|
|
LL = KF - 3
|
|
LLS = KFS - 3
|
|
DO 30 L = 1,LINLEN
|
|
LL = LL + 3
|
|
LLS = LLS + 3
|
|
IF(TRANF(SEQ(LLS),PAA).EQ.TRANF(SEQ(LL),PAA)) THEN
|
|
LINE(L:L) = '-'
|
|
ELSE
|
|
LINE(L:L) = TRANF(SEQ(LL),PAA)
|
|
END IF
|
|
30 CONTINUE
|
|
LINE(1+(KT-KF+4)/3:) = NAMES(J)(1:15)
|
|
WRITE(IDEV,1004,ERR=60)LINE(1:NOTIRL(LINE,MAXLIN,' '))
|
|
1004 FORMAT(' ',A)
|
|
1006 FORMAT(' ',10I10)
|
|
40 CONTINUE
|
|
WRITE(IDEV,1008)
|
|
1008 FORMAT( )
|
|
50 CONTINUE
|
|
RETURN
|
|
60 CONTINUE
|
|
WRITE(KBOUT,*)' Error writing file'
|
|
END
|
|
SUBROUTINE IANDJ(N,L,I,J)
|
|
C ROUTINE, GIVEN N FILES AND ELEMENT L FROM LIST
|
|
C CALC ROW I AND COL J NOS ASSUMING L LIES IN TOP
|
|
C RIGHT CORNER OF MATRIX
|
|
I=1
|
|
J1=0
|
|
K=0
|
|
10 CONTINUE
|
|
J1=J1+1
|
|
J=J1
|
|
11 CONTINUE
|
|
C REACHED ELEMENT L?
|
|
IF(K.EQ.L)RETURN
|
|
J=J+1
|
|
K=K+1
|
|
C REACHED END OF ROW?
|
|
IF(J.LE.N)GO TO 11
|
|
C YES, POINT TO NEX ROW
|
|
I=I+1
|
|
K=K-1
|
|
GO TO 10
|
|
END
|
|
INTEGER FUNCTION IMAXA(I,N)
|
|
INTEGER I(N)
|
|
IMAXA = I(1)
|
|
DO 10 J=2,N
|
|
IMAXA = MAX(IMAXA,I(J))
|
|
10 CONTINUE
|
|
END
|
|
INTEGER FUNCTION IMINA(I,N)
|
|
INTEGER I(N)
|
|
IMINA = I(1)
|
|
DO 10 J=2,N
|
|
IMINA = MIN(IMINA,I(J))
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE PLTI(IY,NY,
|
|
+XMAX,XMIN,YMAX,YMIN,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
INTEGER IY(NY)
|
|
CALL VECTOM
|
|
CALL FRAME(MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
YF = YMIN
|
|
DO 100 I=1,NY
|
|
XF = I
|
|
YT = IY(I)
|
|
CALL LINE(XF,XF,YF,YT,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
100 CONTINUE
|
|
CALL VT100M
|
|
END
|
|
REAL FUNCTION RMAXA(I,N)
|
|
REAL I(N)
|
|
RMAXA = I(1)
|
|
DO 10 J=2,N
|
|
RMAXA = MAX(RMAXA,I(J))
|
|
10 CONTINUE
|
|
END
|
|
REAL FUNCTION RMINA(I,N)
|
|
REAL I(N)
|
|
RMINA = I(1)
|
|
DO 10 J=2,N
|
|
RMINA = MIN(RMINA,I(J))
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE PLTR(IY,NY,
|
|
+XMAX,XMIN,YMAX,YMIN,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
REAL IY(NY)
|
|
CALL VECTOM
|
|
CALL FRAME(MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
YF = YMIN
|
|
DO 100 I=1,NY
|
|
XF = I
|
|
YT = IY(I)
|
|
CALL LINE(XF,XF,YF,YT,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
100 CONTINUE
|
|
CALL VT100M
|
|
END
|
|
SUBROUTINE COCHAN(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,
|
|
+ICONSA,ISAME,ICONSC,MISMAT,IDEV,SUM)
|
|
INTEGER SJ(6)
|
|
INTEGER SUM(4096),SAME(77),SUMC
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM1)
|
|
CHARACTER AA(18)
|
|
EXTERNAL JCODNO
|
|
DATA AA/'F','L','I','V','S','P','T','A','Y','H','Q','N',
|
|
+'K','D','E','C','R','G'/
|
|
DATA SAME/2,1,2,6,3,4,17,18,19,20,3,33,34,35,4,49,50,51,
|
|
+52,6,5,6,7,8,45,46,4,21,22,23,24,4,37,38,39,40,4,53,54,
|
|
+55,56,2,9,10,2,25,26,2,27,28,2,41,42,2,43,44,2,57,58,2,
|
|
+59,60,2,13,14,6,29,30,31,32,47,48,4,61,62,63,64/
|
|
C CONSERVED ACIDS ICONSA=KSUM
|
|
C CONSERVATIVE CHANGES ICONSC=SUMC
|
|
C UNCHANGED CODONS ISAME=ISUM
|
|
DO 55 I=1,4096
|
|
SUM(I)=0
|
|
55 CONTINUE
|
|
C
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
C SET MIN TO NUMBER OF CODONS
|
|
MINL=MINL/3
|
|
MAXL=MAX(IDIM3,IDIM4)
|
|
MAXL=MAXL/3
|
|
J=I1
|
|
K=I2
|
|
C
|
|
C MAIN LOOP
|
|
C
|
|
DO 100 I=1,MINL
|
|
IP1=JCODNO(SEQ1,IDIM1,J)
|
|
IF(IP1.EQ.0)GO TO 199
|
|
IP2=JCODNO(SEQ2,IDIM2,K)
|
|
IF(IP2.EQ.0)GO TO 199
|
|
IP=IP1*64+IP2-64
|
|
SUM(IP)=SUM(IP)+1
|
|
199 CONTINUE
|
|
J=J+3
|
|
K=K+3
|
|
100 CONTINUE
|
|
C
|
|
C SUM THOSE CONSERVED
|
|
C
|
|
ISUM=0
|
|
SUMC=0
|
|
IK1=2
|
|
C
|
|
C LOOP FOR EACH AMINO ACID WITH MORE THAN 1 CODON
|
|
C
|
|
DO 120 K=1,18
|
|
NUM=SAME(IK1-1)
|
|
C
|
|
C FIRST ELEMENT IN SAME IS THE NUMBER OF CODONS FOR THIS AMINO ACID
|
|
C SUCCEEDING ELEMENTS ARE THE CODON NUMBERS
|
|
C LOOP FOR EACH CODON
|
|
C
|
|
DO 110 I=1,NUM
|
|
IP=SAME(IK1+I-1)
|
|
NP=(IP-1)*64
|
|
C
|
|
C LOOP FOR EACH CODON CALCULATING POINTERS TO EACH CHANGE NOT
|
|
C AFFECTING AMINO ACID
|
|
C
|
|
DO 105 J=1,NUM
|
|
MP=SAME(IK1+J-1)
|
|
JP=NP+MP
|
|
SJ(J)=SUM(JP)
|
|
105 CONTINUE
|
|
DO 104 N=1,NUM
|
|
IF(N.EQ.I)ISUM=ISUM+SJ(N)
|
|
SUMC=SUMC+SJ(N)
|
|
SJ(N)=0
|
|
104 CONTINUE
|
|
110 CONTINUE
|
|
IK1=IK1+NUM+1
|
|
120 CONTINUE
|
|
C ADD MET AND TRP CODONS TO THOSE CONSERVING ACID
|
|
KSUM=SUMC+SUM(2276)+SUM(976)
|
|
FKSUM=100.*FLOAT(KSUM)/FLOAT(MAXL)
|
|
WRITE(IDEV,1005)KSUM,FKSUM
|
|
1005 FORMAT(' Number of conserved amino acids ',I6,' = ',F6.2,'%')
|
|
C1005 FORMAT(' NUMBER OF CONSERVED AMINO ACIDS=',I6,', AS PERCENT='
|
|
C 1,F6.2)
|
|
C
|
|
C SUBTRACT CONSERVED CODONS FROM THOSE CONSERVING AMINO ACID
|
|
SUMC=SUMC-ISUM
|
|
FSUMC=100.*FLOAT(SUMC)/FLOAT(MAXL)
|
|
WRITE(IDEV,1001)SUMC,FSUMC
|
|
1001 FORMAT(' Number of conservative codon changes ',I6,' = ',F6.2,'%')
|
|
C1001 FORMAT(' NUMBER OF CONSERVATIVE CODON CHANGES=',I6,', AS',
|
|
C 1' PERCENT=',
|
|
C 1F6.2)
|
|
C ADD MET AND TRP TO UNCHANGED CODONS
|
|
ISUM=ISUM+SUM(2276)+SUM(976)
|
|
FISUM=100.*FLOAT(ISUM)/FLOAT(MAXL)
|
|
WRITE(IDEV,1004)ISUM,FISUM
|
|
1004 FORMAT(' Number of unchanged codons ',I6,' = ',F6.2,'%')
|
|
C1004 FORMAT(' NUMBER OF UNCHANGED CODONS=',I6,', AS PERCENT=',F6.2)
|
|
C
|
|
C CONSERVED ACIDS
|
|
ICONSA=KSUM
|
|
C CONSERVATIVE CHANGES
|
|
ICONSC=SUMC
|
|
C UNCHANGED CODONS
|
|
ISAME=ISUM
|
|
C CODONS CHANGING AMINO ACID
|
|
MISMAT=MAXL-ICONSA
|
|
RETURN
|
|
END
|
|
SUBROUTINE BSCHN2(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,MISMAT,
|
|
1IDEV)
|
|
C THIS ROUTINE COUNTS GENETIC EVENTS AND 1 BASE INSERTION
|
|
C EQUALS N BASE INSERTION
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2)
|
|
INTEGER CTONUM
|
|
EXTERNAL CTONUM
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
C SET MIN TO NUMBER OF BASES
|
|
MAXL=MAX(IDIM3,IDIM4)
|
|
MISMAT=0
|
|
J=I2-1
|
|
IDASH1=0
|
|
IDASH2=0
|
|
DO 5 I=I1,I1+MINL-1
|
|
J=J+1
|
|
IF(CTONUM(SEQ1(I)).EQ.5)IDASH1=IDASH1+1
|
|
IF(CTONUM(SEQ1(I)).NE.5)IDASH1=0
|
|
IF(CTONUM(SEQ2(J)).EQ.5)IDASH2=IDASH2+1
|
|
IF(CTONUM(SEQ2(J)).NE.5)IDASH2=0
|
|
IF(SEQ1(I).EQ.SEQ2(J))GO TO 5
|
|
MISMAT=MISMAT+1
|
|
IF((IDASH1.GT.1).OR.(IDASH2.GT.1))MISMAT=MISMAT-1
|
|
5 CONTINUE
|
|
C cantor and jukes
|
|
P = REAL(MISMAT)/REAL(MAXL)
|
|
D = -0.75*LOG(1.0 - 4.0*P/3.0)
|
|
V = P*(1.0-P)/(REAL(MAXL)*(1.0-4.0*P/3.0)**2)
|
|
CONS=MINL-MISMAT
|
|
CONS=CONS*100./MAXL
|
|
WRITE(IDEV,1001)MISMAT,CONS
|
|
1001 FORMAT(' Number of events ',I6,' = ',F6.2,'%')
|
|
C1001 FORMAT(' NUMBER OF EVENTS=',I6,' PERCENTAGE CONSERVATION=',
|
|
C 1F6.2)
|
|
WRITE(IDEV,1002)D,V
|
|
1002 FORMAT(' Expected number of substitutions per site and variance '
|
|
+,F10.5,F10.7)
|
|
C1002 FORMAT(
|
|
C +' JC SUBSTITUTIONS PER SITE=',F10.5,' AND VARIANCE=',F10.7)
|
|
END
|
|
C
|
|
SUBROUTINE BSCHN1(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,MISMAT,
|
|
+IDEV)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2)
|
|
INTEGER CTONUM
|
|
EXTERNAL CTONUM
|
|
C COUNT DIFFERENCES (BASE CHANGES)
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
C SET MIN TO NUMBER OF BASES
|
|
MAXL=MAX(IDIM3,IDIM4)
|
|
MISMAT=0
|
|
J=I2-1
|
|
DO 10 I=I1,I1+MINL-1
|
|
J=J+1
|
|
IF(CTONUM(SEQ1(I)).NE.CTONUM(SEQ2(J))) MISMAT=MISMAT+1
|
|
10 CONTINUE
|
|
CONS=MINL-MISMAT
|
|
CONS=CONS*100./MAXL
|
|
WRITE(IDEV,1001)MISMAT,CONS
|
|
1001 FORMAT(' Number of differences ',I6,' = ',F6.2,'%')
|
|
END
|
|
SUBROUTINE BSCHN6(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,MISMAT,
|
|
1IDEV)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2)
|
|
C COUNT DIFFERENCES (BASE CHANGES)
|
|
INTEGER CTONUM
|
|
EXTERNAL CTONUM
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
C SET MIN TO NUMBER OF BASES
|
|
MAXL=MAX(IDIM3,IDIM4)
|
|
MISMAT=0
|
|
J=I2-1
|
|
DO 10 I=I1,I1+MINL-1
|
|
J=J+1
|
|
K1 = CTONUM(SEQ1(I))
|
|
K2 = CTONUM(SEQ2(J))
|
|
IF(K1.EQ.K2) GO TO 10
|
|
IF(K1.EQ.5) GO TO 10
|
|
IF(K2.EQ.5) GO TO 10
|
|
MISMAT=MISMAT+1
|
|
10 CONTINUE
|
|
CONS=MINL-MISMAT
|
|
CONS=CONS*100./MAXL
|
|
WRITE(IDEV,1001)MISMAT,CONS
|
|
1001 FORMAT(' NUMBER OF DIFFERENCES=',I6,' PERCENTAGE CONSERVATION',
|
|
1' = ',F6.2)
|
|
RETURN
|
|
END
|
|
C
|
|
SUBROUTINE BSCHN4(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,DIFFER)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2)
|
|
INTEGER DIFFER(IDIM1)
|
|
INTEGER CTONUM
|
|
EXTERNAL CTONUM
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
J1 = I1 - 1
|
|
J2 = I2 - 1
|
|
DO 10 I=1,MINL
|
|
J1 = J1 + 1
|
|
J2 = J2 + 1
|
|
IF(CTONUM(SEQ1(J1)).NE.CTONUM(SEQ2(J2)))
|
|
+ DIFFER(I)=DIFFER(I)+1
|
|
10 CONTINUE
|
|
END
|
|
C
|
|
SUBROUTINE BSCHN3(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,CHANGE)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2)
|
|
INTEGER CHANGE(5,5)
|
|
INTEGER CTONUM
|
|
EXTERNAL CTONUM
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
C SET MIN TO NUMBER OF BASES
|
|
MAXL=MAX(IDIM3,IDIM4)
|
|
MISMAT=0
|
|
J=I2-1
|
|
DO 10 I=I1,I1+MINL-1
|
|
J=J+1
|
|
II=CTONUM(SEQ1(I))
|
|
JJ=CTONUM(SEQ2(J))
|
|
CHANGE(II,JJ)=CHANGE(II,JJ)+1
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE BSCHM3(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,CHANGE)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2)
|
|
INTEGER CHANGE(5,5,3)
|
|
INTEGER CTONUM
|
|
EXTERNAL CTONUM
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
C SET MIN TO NUMBER OF BASES
|
|
MAXL=MAX(IDIM3,IDIM4)
|
|
MISMAT=0
|
|
J=I2-3
|
|
DO 10 I=I1,I1+MINL-2,3
|
|
J=J+3
|
|
DO 5 K=0,2
|
|
II=CTONUM(SEQ1(I+K))
|
|
JJ=CTONUM(SEQ2(J+K))
|
|
CHANGE(II,JJ,K+1)=CHANGE(II,JJ,K+1)+1
|
|
5 CONTINUE
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE KIMURA(CHANGE,DIST)
|
|
INTEGER CHANGE(5,5)
|
|
EXTERNAL EVDIST
|
|
C
|
|
C Calc P and Q, then calc the evolutionary distance
|
|
C
|
|
P = CHANGE(1,2) + CHANGE(2,1) + CHANGE(3,4) + CHANGE(4,3)
|
|
Q = CHANGE(1,3) + CHANGE(3,1) + CHANGE(1,4) + CHANGE(4,1)
|
|
Q = Q + CHANGE(3,2) + CHANGE(2,3) + CHANGE(2,4) + CHANGE(4,2)
|
|
NSITES = 0
|
|
C
|
|
C dont count changes to dash
|
|
C
|
|
DO 10 I=1,4
|
|
DO 5 J=1,4
|
|
NSITES = NSITES + CHANGE(I,J)
|
|
5 CONTINUE
|
|
10 CONTINUE
|
|
WRITE(*,*)'NP=',P
|
|
WRITE(*,*)'NQ=',Q
|
|
WRITE(*,*)'NSITES=',NSITES
|
|
P = 58
|
|
Q = 63
|
|
NSITES = 438
|
|
WRITE(*,*)'NP=',P
|
|
WRITE(*,*)'NQ=',Q
|
|
WRITE(*,*)'NSITES=',NSITES
|
|
P = P / REAL(NSITES)
|
|
Q = Q / REAL(NSITES)
|
|
DIST = EVDIST(P,Q)
|
|
END
|
|
REAL FUNCTION EVDIST(P,Q)
|
|
C
|
|
C Calc evolutionary distance per site according to kimura
|
|
C
|
|
C P and Q are the fractions if sites with respectively type I and II
|
|
C substitutions (ie transitions and transversions)
|
|
C
|
|
EVDIST = -0.5 * LOG((1.0 - 2.0*P - Q) * SQRT(1.0 - 2.0 * Q))
|
|
C
|
|
C standard error
|
|
C
|
|
A = 1.0 / (1.0 - 2.0 * P - Q)
|
|
B = 0.5 * (1.0 / ( A - 1.0 / (1.0 - 2.0 * Q)))
|
|
C
|
|
C !!!!!!!!!!!!!! NOT SURE I CAN READ THIS IN THE PAPER
|
|
C AND CANT READ THE NEXT BIT
|
|
END
|
|
C
|
|
C SUBROUTINE TO WRITE 3 LINES OF SEQUENCES
|
|
C
|
|
SUBROUTINE FMT2(IDEV,SEQ1,SEQ2,IDIM,ISW,ISX,LINLEN)
|
|
CHARACTER SEQ1(IDIM),SEQ2(IDIM),MATCH(120)
|
|
INTEGER KL(12)
|
|
C
|
|
C
|
|
ISXX=ISX-1
|
|
ISWW=ISW-1
|
|
IE=0
|
|
1 CONTINUE
|
|
WRITE(IDEV,1003)
|
|
1003 FORMAT( )
|
|
DO 50 J=1,LINLEN/10
|
|
ISWW=ISWW+10
|
|
KL(J)=ISWW
|
|
50 CONTINUE
|
|
WRITE(IDEV,1001)(KL(KKK),KKK=1,LINLEN/10)
|
|
IS=IE+1
|
|
IE=IE+LINLEN
|
|
IF(IE.GT.IDIM)IE=IDIM
|
|
C COMPARE THE TWO SECTIONS OF SEQUENCE
|
|
IL=IE-IS+1
|
|
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( ' ',12(10A1,1X))
|
|
C SET UP DECIMAL COUNTERS
|
|
DO 60 J=1,LINLEN/10
|
|
ISXX=ISXX+10
|
|
KL(J)=ISXX
|
|
60 CONTINUE
|
|
WRITE(IDEV,1001)(KL(KKK),KKK=1,LINLEN/10)
|
|
1001 FORMAT( ' ',12(5X,I6))
|
|
IF(IE.EQ.IDIM)RETURN
|
|
GO TO 1
|
|
END
|
|
C
|
|
C
|
|
C
|
|
C RETURNS POINTER TO CODON NUMBER FROM 1 - 64
|
|
INTEGER FUNCTION JCODNO(SEQ,IDIM,I)
|
|
C AUTHOR: RODGER STADEN
|
|
CHARACTER SEQ(IDIM)
|
|
INTEGER CTONUM,JP(3)
|
|
EXTERNAL CTONUM
|
|
C CHECK FOR BAD CHAR IN SEQ
|
|
JCODNO=0
|
|
IP=I-1
|
|
DO 1 J=1,3
|
|
JP(J)=CTONUM(SEQ(IP+J))
|
|
1 IF(JP(J).EQ.5)RETURN
|
|
JCODNO=JP(1)*16+4*JP(2)+JP(3)-20
|
|
RETURN
|
|
END
|
|
C
|
|
SUBROUTINE BSCHN5(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,MISMAT,
|
|
1IDEV,PAA)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2),PAA(5,5,5)
|
|
INTEGER CFCOD1
|
|
EXTERNAL CFCOD1
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
J1=I1 - 3
|
|
J2 = I2 - 3
|
|
ISUM = 0
|
|
JSUM = 0
|
|
DO 10 I=1,MINL-2,3
|
|
J1 = J1 + 3
|
|
J2 = J2 + 3
|
|
ISUM = ISUM + CFCOD1(SEQ1(J1),SEQ2(J2),PAA)
|
|
JSUM = JSUM + CFCOD1(SEQ2(J2),SEQ1(J1),PAA)
|
|
10 CONTINUE
|
|
WRITE(IDEV,1000)ISUM
|
|
WRITE(IDEV,1001)JSUM
|
|
1000 FORMAT(' 1 TO 2 =',I6)
|
|
1001 FORMAT(' 2 TO 1 =',I6)
|
|
MISMAT = ISUM + JSUM
|
|
END
|
|
SUBROUTINE BSCHN7(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,MISMAT,
|
|
+IDEV,PAA)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2),PAA(5,5,5)
|
|
INTEGER CFCOD1
|
|
EXTERNAL CFCOD1
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
J1=I1 - 3
|
|
J2 = I2 - 3
|
|
MISMAT = 0
|
|
DO 10 I=1,MINL-2,3
|
|
J1 = J1 + 3
|
|
J2 = J2 + 3
|
|
MISMAT = MISMAT + CFCOD1(SEQ1(J1),SEQ2(J2),PAA)
|
|
10 CONTINUE
|
|
WRITE(IDEV,1000)MISMAT
|
|
1000 FORMAT(' Expressed base changes ',I6)
|
|
END
|
|
SUBROUTINE BSCHN8(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,MISMAT,
|
|
+IDEV,PAA)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2),PAA(5,5,5)
|
|
INTEGER CFCOD2
|
|
EXTERNAL CFCOD2
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
J1=I1 - 3
|
|
J2 = I2 - 3
|
|
MISMAT = 0
|
|
DO 10 I=1,MINL-2,3
|
|
J1 = J1 + 3
|
|
J2 = J2 + 3
|
|
MISMAT = MISMAT + CFCOD2(SEQ1(J1),SEQ2(J2),PAA)
|
|
10 CONTINUE
|
|
WRITE(IDEV,1000)MISMAT
|
|
1000 FORMAT(' Silent base changes ',I6)
|
|
END
|
|
SUBROUTINE BSCHN9(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,MISMAT,
|
|
+IDEV,PAA)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2),PAA(5,5,5)
|
|
INTEGER CFCOD3
|
|
EXTERNAL CFCOD3
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
J1=I1 - 3
|
|
J2 = I2 - 3
|
|
MISMAT = 0
|
|
DO 10 I=1,MINL-2,3
|
|
J1 = J1 + 3
|
|
J2 = J2 + 3
|
|
MISMAT = MISMAT + CFCOD3(SEQ1(J1),SEQ2(J2),PAA)
|
|
10 CONTINUE
|
|
WRITE(IDEV,1000)MISMAT
|
|
1000 FORMAT(' Number of changed acids ',I6)
|
|
END
|
|
SUBROUTINE COCHN1(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,DIFFER,
|
|
+PAA)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2),PAA(5,5,5)
|
|
INTEGER DIFFER(IDIM1)
|
|
C count number of expressed base changes
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
J1=I1 - 3
|
|
J2 = I2 - 3
|
|
DO 10 I=1,MINL-2,3
|
|
J1 = J1 + 3
|
|
J2 = J2 + 3
|
|
CALL SCFCD1(SEQ1(J1),SEQ2(J2),PAA,DIFFER(I))
|
|
C DIFFER(I) = DIFFER(I) + CFCOD1(SEQ1(J1),SEQ2(J2),PAA)
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE SCFCD1(COD1,COD2,PAA,DIFF)
|
|
CHARACTER COD1(3),COD2(3),COD(3),ACID,TRANF,PAA(5,5,5)
|
|
INTEGER DIFF(3)
|
|
EXTERNAL TRANF
|
|
C
|
|
C returns 1 for each base in a codon cod2 that changes the acid from that
|
|
C encoded by cod1
|
|
C
|
|
ACID = TRANF(COD1,PAA)
|
|
DO 10 I = 1,3
|
|
COD(I) = COD1(I)
|
|
10 CONTINUE
|
|
DO 20 I = 1,3
|
|
COD(I) = COD2(I)
|
|
IF(TRANF(COD,PAA).NE.ACID)DIFF(I) = DIFF(I) + 1
|
|
COD(I) = COD1(I)
|
|
20 CONTINUE
|
|
END
|
|
INTEGER FUNCTION CFCOD1(COD1,COD2,PAA)
|
|
CHARACTER COD1(3),COD2(3),COD(3),ACID,TRANF,PAA(5,5,5)
|
|
EXTERNAL TRANF
|
|
C
|
|
C returns 1 for each base in a codon cod2 that changes the acid from that
|
|
C encoded by cod1
|
|
C
|
|
CFCOD1 = 0
|
|
ACID = TRANF(COD1,PAA)
|
|
DO 10 I = 1,3
|
|
COD(I) = COD1(I)
|
|
10 CONTINUE
|
|
DO 20 I = 1,3
|
|
COD(I) = COD2(I)
|
|
IF(TRANF(COD,PAA).NE.ACID)CFCOD1 = CFCOD1 + 1
|
|
COD(I) = COD1(I)
|
|
20 CONTINUE
|
|
END
|
|
SUBROUTINE COCHN2(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,DIFFER,
|
|
+PAA)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2),PAA(5,5,5)
|
|
INTEGER DIFFER(IDIM1)
|
|
C count number of silent base changes
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
J1=I1 - 3
|
|
J2 = I2 - 3
|
|
DO 10 I=1,MINL-2,3
|
|
J1 = J1 + 3
|
|
J2 = J2 + 3
|
|
CALL SCFCD2(SEQ1(J1),SEQ2(J2),PAA,DIFFER(I))
|
|
C DIFFER(I) = DIFFER(I) + CFCOD2(SEQ1(J1),SEQ2(J2),PAA)
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE SCFCD2(COD1,COD2,PAA,DIFF)
|
|
CHARACTER COD1(3),COD2(3),COD(3),ACID,TRANF,PAA(5,5,5)
|
|
INTEGER DIFF(3)
|
|
EXTERNAL TRANF
|
|
C
|
|
C compare each base, if different does this base alone change the acid
|
|
C if not then add 1 to score
|
|
ACID = TRANF(COD1,PAA)
|
|
DO 10 I = 1,3
|
|
COD(I) = COD1(I)
|
|
10 CONTINUE
|
|
DO 20 I = 1,3
|
|
IF(COD1(I).NE.COD2(I)) THEN
|
|
COD(I) = COD2(I)
|
|
IF(TRANF(COD,PAA).EQ.ACID)DIFF(I) = DIFF(I) + 1
|
|
COD(I) = COD1(I)
|
|
END IF
|
|
20 CONTINUE
|
|
END
|
|
INTEGER FUNCTION CFCOD2(COD1,COD2,PAA)
|
|
CHARACTER COD1(3),COD2(3),COD(3),ACID,TRANF,PAA(5,5,5)
|
|
EXTERNAL TRANF
|
|
C
|
|
C compare each base, if different does this base alone change the acid
|
|
C if not then add 1 to score
|
|
CFCOD2 = 0
|
|
ACID = TRANF(COD1,PAA)
|
|
DO 10 I = 1,3
|
|
COD(I) = COD1(I)
|
|
10 CONTINUE
|
|
DO 20 I = 1,3
|
|
IF(COD1(I).NE.COD2(I)) THEN
|
|
COD(I) = COD2(I)
|
|
IF(TRANF(COD,PAA).EQ.ACID)CFCOD2 = CFCOD2 + 1
|
|
COD(I) = COD1(I)
|
|
END IF
|
|
20 CONTINUE
|
|
END
|
|
INTEGER FUNCTION CFCOD3(COD1,COD2,PAA)
|
|
CHARACTER COD1(3),COD2(3),TRANF,PAA(5,5,5)
|
|
EXTERNAL TRANF
|
|
C
|
|
C compare each acid, count if different
|
|
C
|
|
CFCOD3 = 0
|
|
IF(TRANF(COD1,PAA).NE.TRANF(COD2,PAA)) CFCOD3 = 1
|
|
END
|
|
REAL FUNCTION FINF(COMP)
|
|
INTEGER COMP(5)
|
|
F = 0.
|
|
S = 5.
|
|
DO 10 I=1,5
|
|
S = S + COMP(I)
|
|
10 CONTINUE
|
|
DO 20 I=1,5
|
|
X = (1.0+REAL(COMP(I)))/S
|
|
F = F - X * LOG(X)
|
|
20 CONTINUE
|
|
FINF = F
|
|
END
|
|
SUBROUTINE COCHN3(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,CHANGE,
|
|
+PAA)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2),PAA(5,5,5)
|
|
INTEGER CHANGE(5,5)
|
|
C count number of expressed base changes
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
J1=I1 - 3
|
|
J2 = I2 - 3
|
|
DO 10 I=1,MINL-2,3
|
|
J1 = J1 + 3
|
|
J2 = J2 + 3
|
|
CALL SCFCD3(SEQ1(J1),SEQ2(J2),PAA,CHANGE)
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE SCFCD3(COD1,COD2,PAA,CHANGE)
|
|
CHARACTER COD1(3),COD2(3),COD(3),ACID,TRANF,PAA(5,5,5)
|
|
INTEGER CHANGE(5,5),CTONUM
|
|
EXTERNAL TRANF,CTONUM
|
|
C
|
|
C compare each base, if different does this base alone change the acid
|
|
C if so then add 1 to score
|
|
ACID = TRANF(COD1,PAA)
|
|
DO 10 I = 1,3
|
|
COD(I) = COD1(I)
|
|
10 CONTINUE
|
|
DO 20 I = 1,3
|
|
IF(COD1(I).NE.COD2(I)) THEN
|
|
COD(I) = COD2(I)
|
|
IF(TRANF(COD,PAA).NE.ACID) THEN
|
|
K1 = CTONUM(COD1(I))
|
|
K2 = CTONUM(COD2(I))
|
|
CHANGE(K1,K2) = CHANGE(K1,K2) + 1
|
|
END IF
|
|
COD(I) = COD1(I)
|
|
END IF
|
|
20 CONTINUE
|
|
END
|
|
SUBROUTINE COCHN5(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,CHANGE,
|
|
+PAA)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2),PAA(5,5,5)
|
|
INTEGER CHANGE(5,5)
|
|
C count number of expressed base changes
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
J1=I1 - 3
|
|
J2 = I2 - 3
|
|
DO 10 I=1,MINL-2,3
|
|
J1 = J1 + 3
|
|
J2 = J2 + 3
|
|
CALL SCFCD4(SEQ1(J1),SEQ2(J2),PAA,CHANGE)
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE SCFCD4(COD1,COD2,PAA,CHANGE)
|
|
CHARACTER COD1(3),COD2(3),TRANF,PAA(5,5,5)
|
|
INTEGER CHANGE(5,5),CTONUM
|
|
EXTERNAL TRANF,CTONUM
|
|
C
|
|
C if the acid is the same count the mutations
|
|
C
|
|
IF (TRANF(COD1,PAA).EQ.TRANF(COD2,PAA)) THEN
|
|
DO 20 I = 1,3
|
|
K1 = CTONUM(COD1(I))
|
|
K2 = CTONUM(COD2(I))
|
|
CHANGE(K1,K2) = CHANGE(K1,K2) + 1
|
|
20 CONTINUE
|
|
END IF
|
|
END
|
|
SUBROUTINE TRAN(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,
|
|
+KBIN,KBOUT,IDEV,LINLEN,PAA)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2),SYM1,SYM2,MATCH(120),TRANF
|
|
CHARACTER PAA(5,5,5)
|
|
INTEGER CTONUM
|
|
EXTERNAL CTONUM,TRANF
|
|
SYM2 = ':'
|
|
SYM1 = ' '
|
|
IF(0.EQ.1) THEN
|
|
WRITE(KBOUT,2000)
|
|
2000 FORMAT(' ? IDENTITY SYMBOL=',$)
|
|
READ(KBIN,2001)SYM2
|
|
2001 FORMAT(A1)
|
|
WRITE(KBOUT,2002)
|
|
2002 FORMAT(' ? MISMATCH SYMBOL=',$)
|
|
READ(KBIN,2001)SYM1
|
|
END IF
|
|
J1=I1
|
|
J2=I1+LINLEN-1
|
|
IF(J2.GT.N1)J2=N1
|
|
K1=I2
|
|
K2=I2+LINLEN-1
|
|
IF(K2.GT.N2)K2=N2
|
|
20 CONTINUE
|
|
WRITE(IDEV,1004)
|
|
WRITE(IDEV,1000)(TRANF(SEQ1(K),PAA),K=J1,J2-1,3)
|
|
WRITE(IDEV,1003)(SEQ1(K),K=J1,J2)
|
|
DO 30 JJ=1,LINLEN
|
|
MATCH(JJ)=SYM1
|
|
30 CONTINUE
|
|
DO 40 JJ=1,MIN((J2-J1+1),(K2-K1+1),(N1-J1+1),(N2-K1+1))
|
|
KK=CTONUM(SEQ1(J1-1+JJ))
|
|
LL=CTONUM(SEQ2(K1-1+JJ))
|
|
IF(KK.EQ.LL)MATCH(JJ)=SYM2
|
|
40 CONTINUE
|
|
WRITE(IDEV,1003)(MATCH(K),K=1,MIN((J2-J1+1),(K2-K1+1),
|
|
+(N1-J1+1),(N2-K1+1)))
|
|
WRITE(IDEV,1003)(SEQ2(K),K=K1,K2)
|
|
WRITE(IDEV,1000)(TRANF(SEQ2(K),PAA),K=K1,K2-1,3)
|
|
C
|
|
J1=J2+1
|
|
J2=J2+LINLEN
|
|
IF(J1.GT.N1)RETURN
|
|
IF(J2.GT.N1)J2=N1
|
|
K1=K2+1
|
|
K2=K2+LINLEN
|
|
IF(K1.GT.N2)RETURN
|
|
IF(K2.GT.N2)K2=N2
|
|
GO TO 20
|
|
1000 FORMAT(3X,40(2X,1A1))
|
|
1001 FORMAT(4X,40(2X,1A1))
|
|
1002 FORMAT(5X,40(2X,1A1))
|
|
1003 FORMAT(4X,120A1)
|
|
1004 FORMAT( )
|
|
1005 FORMAT(4X,12I10)
|
|
END
|
|
SUBROUTINE INITMT(M,S,JOB)
|
|
INTEGER S(4096)
|
|
REAL M(16)
|
|
CALL FILLI(S,4096,0)
|
|
C only initialise mutation rates if first pass
|
|
IF(JOB.EQ.0) CALL FILLR(M,16,0.25)
|
|
END
|
|
SUBROUTINE MAKED(M,D)
|
|
REAL M(4,4),D(17)
|
|
INTEGER T,C,A,G
|
|
SAVE T,C,A,G
|
|
DATA T,C,A,G/1,2,3,4/
|
|
C
|
|
C To use these values: 1. add m(x,y) to d where x,y is the selected change
|
|
C 2. correct each of the appropriate m's using the
|
|
C d values and add m(x,y)/d to m(x,y)
|
|
C
|
|
D(1) = M(T,C) + M(T,A) + M(T,G)
|
|
D(2) = M(T,A) + M(T,G)
|
|
D(3) = M(T,G)
|
|
D(4) = 0.
|
|
D(5) = M(C,T) + M(C,A) + M(C,G)
|
|
D(6) = M(C,A) + M(C,G)
|
|
D(7) = 0.
|
|
D(8) = M(A,T) + M(A,C) + M(A,G)
|
|
D(9) = M(A,T) + M(A,C)
|
|
D(10) = M(A,G)
|
|
D(11) = 0.
|
|
D(12) = M(G,T) + M(G,C) + M(G,A)
|
|
D(13) = M(G,T) + M(G,C)
|
|
D(14) = 0.
|
|
D(15) = M(C,T) + M(C,G)
|
|
D(16) = M(C,G)
|
|
D(17) = M(A,T) + M(A,G)
|
|
C DO 10 I=1,17
|
|
C WRITE(*,*)I,D(I)
|
|
C 10 CONTINUE
|
|
END
|
|
SUBROUTINE UPDM(M,D,E,MC,B1,B2,IADD,ND)
|
|
REAL M(4,4),D(17),E(4,4)
|
|
INTEGER T,C,A,G,B1,B2,ND(17)
|
|
SAVE T,C,A,G
|
|
PARAMETER (DELTA = 1E-10)
|
|
DATA T,C,A,G/1,2,3,4/
|
|
IF (IADD.EQ.0) RETURN
|
|
C ND(MC) = ND(MC) + IADD
|
|
ADD = IADD
|
|
C E(B1,B2) = E(B1,B2) + ADD
|
|
C IF(0.EQ.0) RETURN
|
|
C
|
|
C add the rate for the observed mutation to the denominator
|
|
C
|
|
X = D(MC) + M(B1,B2) + DELTA
|
|
C WRITE(*,*)'B1,B2,MC,ADD,D,X',B1,B2,MC,ADD,D(MC),X
|
|
IF (MC.EQ.1) THEN
|
|
C T>T
|
|
IF(0.EQ.0) RETURN
|
|
C WRITE(*,*)'D(1),M(B1,B2)',D(1),M(B1,B2)
|
|
C
|
|
C add in the observed mutation to the appropriate mutation rate
|
|
C which is equivalent to adding (1 - the corrections)
|
|
C
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
C
|
|
C now add in the corrections for the barred bases
|
|
C
|
|
E(B1,C) = E(B1,C) + ADD*M(B1,C)/X
|
|
E(B1,A) = E(B1,A) + ADD*M(B1,A)/X
|
|
E(B1,G) = E(B1,G) + ADD*M(B1,G)/X
|
|
ELSE IF (MC.EQ.2) THEN
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,A) = E(B1,A) + ADD*M(B1,A)/X
|
|
E(B1,G) = E(B1,G) + ADD*M(B1,G)/X
|
|
ELSE IF (MC.EQ.3) THEN
|
|
C WRITE(*,*)'D(3),M(B1,B2)',D(3),M(B1,B2)
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,G) = E(B1,G) + ADD*M(B1,G)/X
|
|
ELSE IF (MC.EQ.4) THEN
|
|
C T>TCAG
|
|
C WRITE(*,*)'D(),M(B1,B2)',D(4),M(B1,B2)
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
ELSE IF (MC.EQ.5) THEN
|
|
C C>C
|
|
IF(0.EQ.0) RETURN
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,T) = E(B1,T) + ADD*M(B1,T)/X
|
|
E(B1,A) = E(B1,A) + ADD*M(B1,A)/X
|
|
E(B1,G) = E(B1,G) + ADD*M(B1,G)/X
|
|
ELSE IF (MC.EQ.6) THEN
|
|
C C>CT
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,A) = E(B1,A) + ADD*M(B1,A)/X
|
|
E(B1,G) = E(B1,G) + ADD*M(B1,G)/X
|
|
ELSE IF (MC.EQ.7) THEN
|
|
C C>TCAG
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
ELSE IF (MC.EQ.8) THEN
|
|
C A>A
|
|
IF(0.EQ.0) RETURN
|
|
C WRITE(*,*)'D(8),M(B1,B2)',D(8),M(B1,B2)
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,T) = E(B1,T) + ADD*M(B1,T)/X
|
|
E(B1,C) = E(B1,C) + ADD*M(B1,C)/X
|
|
E(B1,G) = E(B1,G) + ADD*M(B1,G)/X
|
|
ELSE IF (MC.EQ.9) THEN
|
|
C A>AG
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,T) = E(B1,T) + ADD*M(B1,T)/X
|
|
E(B1,C) = E(B1,C) + ADD*M(B1,C)/X
|
|
ELSE IF (MC.EQ.10) THEN
|
|
C A>TCA
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,G) = E(B1,G) + ADD*M(B1,G)/X
|
|
ELSE IF (MC.EQ.11) THEN
|
|
C A>TCAG
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
ELSE IF (MC.EQ.12) THEN
|
|
C G>G
|
|
IF(0.EQ.0) RETURN
|
|
C WRITE(*,*)'D(12),M(B1,B2)',D(12),M(B1,B2)
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,T) = E(B1,T) + ADD*M(B1,T)/X
|
|
E(B1,C) = E(B1,C) + ADD*M(B1,C)/X
|
|
E(B1,A) = E(B1,A) + ADD*M(B1,A)/X
|
|
ELSE IF (MC.EQ.13) THEN
|
|
C G>AG
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,T) = E(B1,T) + ADD*M(B1,T)/X
|
|
E(B1,C) = E(B1,C) + ADD*M(B1,C)/X
|
|
ELSE IF (MC.EQ.14) THEN
|
|
C G>TCAG
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
ELSE IF (MC.EQ.15) THEN
|
|
C C>CA
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,T) = E(B1,T) + ADD*M(B1,T)/X
|
|
E(B1,G) = E(B1,G) + ADD*M(B1,G)/X
|
|
ELSE IF (MC.EQ.16) THEN
|
|
C C>TCA
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,G) = E(B1,G) + ADD*M(B1,G)/X
|
|
ELSE IF (MC.EQ.17) THEN
|
|
C A>CA
|
|
E(B1,B2) = E(B1,B2) + ADD*M(B1,B2)/X
|
|
E(B1,T) = E(B1,T) + ADD*M(B1,T)/X
|
|
E(B1,G) = E(B1,G) + ADD*M(B1,G)/X
|
|
ELSE
|
|
WRITE(*,*)'SCREAM UPDM',MC
|
|
END IF
|
|
END
|
|
SUBROUTINE UPDE(M,E,MC,SUMS,PAA)
|
|
REAL M(4,4),D(17),E(4,4),T(4),PREV(4,4),DELTA
|
|
INTEGER SUMS(4,4,4,4,4,4),MC(4,4,4,3),CYCLE,ND(17)
|
|
CHARACTER PAA(5,5,5)
|
|
PARAMETER (SMALL = 1E-10)
|
|
SAVE DELTA
|
|
DATA DELTA/0.01/
|
|
DO 1 I=1,4
|
|
DO 1 J=1,4
|
|
PREV(I,J) = M(I,J)
|
|
1 CONTINUE
|
|
DO 3 I=1,17
|
|
ND(I) = 0
|
|
3 CONTINUE
|
|
IR = 0
|
|
5 CONTINUE
|
|
IR = IR + 1
|
|
CALL MAKED(M,D)
|
|
DO 10 I=1,4
|
|
DO 10 J=1,4
|
|
E(I,J) = 0.
|
|
10 CONTINUE
|
|
DO 20 I=1,4
|
|
DO 20 J=1,4
|
|
DO 20 K=1,4
|
|
DO 15 II=1,4
|
|
DO 15 JJ=1,4
|
|
DO 15 KK=1,4
|
|
IF(PAA(K,J,I).EQ.PAA(KK,JJ,II)) THEN
|
|
IF((I.NE.II).OR.(J.NE.JJ).OR.(K.NE.KK)) THEN
|
|
CALL UPDM(M,D,E,MC(K,J,I,1),I,II,SUMS(I,J,K,II,JJ,KK),ND)
|
|
CALL UPDM(M,D,E,MC(K,J,I,2),J,JJ,SUMS(I,J,K,II,JJ,KK),ND)
|
|
CALL UPDM(M,D,E,MC(K,J,I,3),K,KK,SUMS(I,J,K,II,JJ,KK),ND)
|
|
END IF
|
|
END IF
|
|
15 CONTINUE
|
|
20 CONTINUE
|
|
C WRITE(*,*)'NEW E '
|
|
C DO 55 I=1,4
|
|
C WRITE(*,*)(E(I,J),J=1,4)
|
|
C 55 CONTINUE
|
|
CALL FILLR(T,4,SMALL)
|
|
DO 30 I=1,4
|
|
DO 30 J=1,4
|
|
T(I) = T(I) + E(I,J)
|
|
30 CONTINUE
|
|
C WRITE(*,*)'ROW TOTS',T
|
|
DO 40 I=1,4
|
|
DO 40 J=1,4
|
|
M(I,J) = E(I,J)/T(I)
|
|
40 CONTINUE
|
|
C WRITE(*,*)' NEW M'
|
|
C DO 50 I=1,4
|
|
C WRITE(*,1000)(M(I,J),J=1,4)
|
|
C WRITE(*,1000)(PREV(I,J),J=1,4)
|
|
C 50 CONTINUE
|
|
C 1000 FORMAT(4F8.5)
|
|
CYCLE = 0
|
|
DO 60 I=1,4
|
|
DO 60 J=1,4
|
|
IF (ABS(M(I,J)-PREV(I,J)).GT.DELTA) CYCLE = 1
|
|
C PREV(I,J) = M(I,J)
|
|
60 CONTINUE
|
|
C WRITE(*,*)IR
|
|
DO 70 I=1,4
|
|
DO 70 J=1,4
|
|
PREV(I,J) = M(I,J)
|
|
70 CONTINUE
|
|
IF(CYCLE.EQ.1) GO TO 5
|
|
WRITE(*,*)'Number of iterations',IR
|
|
C WRITE(*,*)(K,ND(K),K=1,17)
|
|
END
|
|
SUBROUTINE UPDET(SUMS,PAA,MC)
|
|
INTEGER SUMS(4,4,4,4,4,4),MC(4,4,4,3)
|
|
CHARACTER PAA(5,5,5),BASES(4)
|
|
SAVE BASES
|
|
DATA BASES/'T','C','A','G'/
|
|
DO 20 I=1,4
|
|
DO 20 J=1,4
|
|
DO 20 K=1,4
|
|
DO 15 II=1,4
|
|
DO 15 JJ=1,4
|
|
DO 15 KK=1,4
|
|
IF(PAA(K,J,I).EQ.PAA(KK,JJ,II)) THEN
|
|
C IF(PAA(I,J,K).EQ.PAA(II,JJ,KK)) THEN
|
|
C WRITE(*,*)I,J,K,' to',II,JJ,KK,SUMS(I,J,K,II,JJ,KK)
|
|
WRITE(*,*)BASES(I),BASES(J),BASES(K),
|
|
+BASES(II),BASES(JJ),BASES(KK),MC(K,J,I,1),MC(K,J,I,2),
|
|
+MC(K,J,I,3),SUMS(I,J,K,II,JJ,KK)
|
|
END IF
|
|
15 CONTINUE
|
|
20 CONTINUE
|
|
END
|
|
SUBROUTINE COCHN4(SEQ1,IDIM1,I1,N1,SEQ2,IDIM2,I2,N2,CHANGE)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM2)
|
|
INTEGER CHANGE(4,4,4,4,4,4),CTONUM,K1(3),K2(3)
|
|
EXTERNAL CTONUM
|
|
IDIM3=N1-I1+1
|
|
IDIM4=N2-I2+1
|
|
C GET MINIMUM LENGTH TO SCAN
|
|
MINL=MIN(IDIM3,IDIM4)
|
|
J1=I1 - 3
|
|
J2 = I2 - 3
|
|
DO 10 I=1,MINL-2,3
|
|
J1 = J1 + 3
|
|
J2 = J2 + 3
|
|
DO 5 JJ=1,3
|
|
KK = CTONUM(SEQ1(J1+JJ-1))
|
|
IF(KK.EQ.5) GO TO 10
|
|
K1(JJ) = KK
|
|
KK = CTONUM(SEQ2(J2+JJ-1))
|
|
IF(KK.EQ.5) GO TO 10
|
|
K2(JJ) = KK
|
|
5 CONTINUE
|
|
C WRITE(*,*)I,K1,K2
|
|
CHANGE(K1(1),K1(2),K1(3),K2(1),K2(2),K2(3)) =
|
|
+ CHANGE(K1(1),K1(2),K1(3),K2(1),K2(2),K2(3)) + 1
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE EMT(MTT)
|
|
INTEGER MTT(64,3)
|
|
CHARACTER MTTC(64)*3
|
|
SAVE MTTC
|
|
DATA MTTC/
|
|
+'112','116','219','21D',
|
|
+'154','157','15B','15E',
|
|
+'182','186','199','19D',
|
|
+'1C2','1C6','1D8','1CC',
|
|
+'514','517','61B','61E',
|
|
+'554','557','55B','55E',
|
|
+'582','586','589','58D',
|
|
+'5C4','5C7','FCB','FCE',
|
|
+'813','81G','81A','81C',
|
|
+'854','857','85B','85E',
|
|
+'882','886','889','88D',
|
|
+'8C2','8C6','HC9','HCD',
|
|
+'C14','C17','C1B','C1E',
|
|
+'C54','C57','C5B','C5E',
|
|
+'C82','C86','C89','C8D',
|
|
+'CC4','CC7','CCB','CCE'/
|
|
CALL CTI(MTTC,MTT)
|
|
C CALL WRTCOI(SUM,6,PAA,MTTC)
|
|
C CALL EM(PAA,KBOUT)
|
|
END
|
|
SUBROUTINE CTI(MTTC,MTT)
|
|
CHARACTER MTTC(64)*3
|
|
INTEGER MTT(64,3)
|
|
CHARACTER C*17
|
|
DATA C/'123456789ABCDEFGH'/
|
|
DO 10 J=1,64
|
|
DO 8 M=1,3
|
|
DO 5 K=1,17
|
|
IF(MTTC(J)(M:M).EQ.C(K:K)) THEN
|
|
MTT(J,M) = K
|
|
GO TO 8
|
|
END IF
|
|
5 CONTINUE
|
|
8 CONTINUE
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE WRTCOI(SUM,IDEV,PAA,MTT)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER SUM(4,4,4)
|
|
CHARACTER BASE(4),PAA(5,5,5),MTT(4,4,4)*3
|
|
SAVE BASE
|
|
DATA BASE/'T','C','A','G'/
|
|
C
|
|
WRITE(IDEV,1001)
|
|
1001 FORMAT(6X,'===========================================')
|
|
DO 10 I=1,4
|
|
DO 20 K=1,4
|
|
WRITE(IDEV,1000)(PAA(K,J,I),
|
|
+BASE(I),BASE(J),BASE(K),SUM(I,J,K),J=1,4)
|
|
WRITE(IDEV,1002)(MTT(K,J,I),J=1,4)
|
|
1002 FORMAT(5X,4(1X,1X,A,5X))
|
|
20 CONTINUE
|
|
10 WRITE(IDEV,1001)
|
|
1000 FORMAT(5X,4(1X,A1,1X,3A1,I5))
|
|
END
|