5041 lines
158 KiB
Fortran
5041 lines
158 KiB
Fortran
C 26-1-93 added new consensus routine
|
|
C added file busy check
|
|
C 9-7-92 Changed remove contig line to remove all readings as well
|
|
C 4-6-92 Added pad shifting routine
|
|
C 6-4-92 fixed a bug in autocn revealed by new update method:
|
|
C for some cases where more than 2 overlaps were being found
|
|
C i was moving the wrong arrays around!
|
|
C 2-4-92 added new dbauto with new update consensus method
|
|
C and new consensus calculation
|
|
C 2-4-92 brought uptodate with dap
|
|
C 29.05.91 IMPLEMENTED REMOVE GEL READING
|
|
C 21-8-91 Added routines to find internal overlaps
|
|
C 2-9-91 Fixed bug in copytg
|
|
C 8-11-91 fixed bugs in "find internal joins"
|
|
C 23-Jun-92 COPYTAG - params to READCC and WRITCC in wrong order
|
|
C
|
|
SUBROUTINE NEWCON(RELPG,LNGTHG,LNBR,RNBR,MAXDB,IDBSIZ,
|
|
+NGELS,NCONTS,MAXGEL,LLINO,LINCON,KOPT,
|
|
+TEMP3,SEQ1,MAXSEQ,SEQ2,SEQ3,
|
|
+MAXGLM,MAXGL2,CHRSIZ,ECHRSZ,
|
|
+ILADD,IRADD,MAXCON,
|
|
+KBIN,KBOUT,IDEV1,IDEV2,IDEV3,IDEV,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,NAMARC,NAMPRO,FILE,
|
|
+PERCD,IDM,IDEVC,IDEVT)
|
|
INTEGER CHRSIZ,ECHRSZ,RREG
|
|
INTEGER RELPG(MAXDB)
|
|
INTEGER LNGTHG(MAXDB),LNBR(MAXDB),RNBR(MAXDB)
|
|
INTEGER TEMP3(ECHRSZ,MAXGL2),ANS
|
|
INTEGER ILADD(MAXCON),IRADD(MAXCON)
|
|
CHARACTER SEQ3(MAXGLM)
|
|
CHARACTER SEQ1(MAXSEQ),SEQ2(MAXGLM)
|
|
CHARACTER NAMARC*(*),NAMPRO*(*),FILE*(*),HELPF*(*)
|
|
PARAMETER (MAXPRM = 6)
|
|
CHARACTER PROMPT(2)*(MAXPRM)
|
|
IF(NGELS.LT.1) RETURN
|
|
CALL DBCHEK(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ TEMP3,IERR,KBOUT)
|
|
IF(IERR.GT.1) RETURN
|
|
PROMPT(1) = 'Staden'
|
|
PROMPT(2) = 'FASTA'
|
|
C
|
|
C idim is current consensus length
|
|
C
|
|
IDIM = 0
|
|
ANS = 0
|
|
IWING = 0
|
|
FILE = ' '
|
|
CALL OPENF1(IDEV,FILE,1,IOK,KBIN,KBOUT,
|
|
+'Name for consensus file',
|
|
+IHELPS,IHELPE,FILEH,IDEVH)
|
|
IF(IOK.NE.0)RETURN
|
|
IF (KOPT.EQ.1) THEN
|
|
CALL YESNO(I,'Use clipped data',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(I.LT.0) GO TO 400
|
|
IF(I.EQ.0) THEN
|
|
MN = 1
|
|
MX = MAXGEL
|
|
IWING = 100
|
|
CALL GETINT(MN,MX,IWING,
|
|
+ 'Window size for good data scan',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 400
|
|
IWING = IVAL
|
|
MN = 1
|
|
MX = MIN(100,IWING)
|
|
NBAD = MIN(IWING,5)
|
|
CALL GETINT(MN,MX,NBAD,
|
|
+ 'Maximum number of dashes in scan window',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 400
|
|
NBAD = IVAL
|
|
END IF
|
|
END IF
|
|
CALL YESNO(ANS,'Make consensus for whole database',
|
|
+IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
IF(ANS.LT.0) GO TO 400
|
|
IF(ANS.EQ.1)GO TO 150
|
|
IDIM2=MAXGEL
|
|
C
|
|
C calc the consensus and add the unused data if required
|
|
C save the lengths of the unused data in iladd, iradd
|
|
C (we dont use them here!)
|
|
C
|
|
CALL FILLI(ILADD,NCONTS,0)
|
|
CALL FILLI(IRADD,NCONTS,0)
|
|
JOB = 0
|
|
CALL JCONS(RELPG,LNGTHG,LNBR,RNBR,NAMPRO,NGELS,NCONTS,
|
|
+SEQ1,MAXSEQ,SEQ2,IDBSIZ,IDIM,JOB,KDUMM,KDUMM,KDUMM,KDUMM,
|
|
+TEMP3,
|
|
+ECHRSZ,MAXGL2,KBOUT,IDEV2,IFAIL,MAXGEL,IDM,PERCD,SEQ3,
|
|
+ILADD,IRADD,MAXCON,KDUMM,IWING,NBAD)
|
|
IF(IFAIL.NE.0) GO TO 400
|
|
ANS = 1
|
|
CALL RADION('Select output format',PROMPT,2,ANS,
|
|
+IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
IF (ANS.LT.1) GO TO 400
|
|
IF(ANS.EQ.1) THEN
|
|
CALL FMTDK(IDEV,SEQ1,IDIM)
|
|
ELSE
|
|
CALL WRITCF(IDEV,SEQ1,IDIM,NAMPRO,KBOUT,IOK)
|
|
END IF
|
|
RETURN
|
|
150 CONTINUE
|
|
CALL GETLN3(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+LINCON,LLINO,NULGEL,IERR,IDBSIZ,KBIN,KBOUT,IDEV3,
|
|
+'Contig identifier',
|
|
+IHELPS,IHELPE,FILEH,IDEVH)
|
|
IF(IERR.NE.0)GO TO 400
|
|
CALL GETREG(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+1,RELPG(LINCON),LREG,RREG,LINCON,LLINO,IDBSIZ,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IERR)
|
|
IF(IERR.NE.0)GO TO 400
|
|
IDIM2=RREG-LREG+1
|
|
IF((IDIM+20+IDIM2).GT.MAXSEQ)THEN
|
|
CALL ERROM(KBOUT,'Maximum consensus length exceeded')
|
|
CLOSE(UNIT=IDEV)
|
|
RETURN
|
|
END IF
|
|
JOB = 1
|
|
JCON = 1
|
|
CALL JCONS(RELPG,LNGTHG,LNBR,RNBR,NAMPRO,NGELS,NCONTS,
|
|
+ SEQ1,MAXSEQ,SEQ2,IDBSIZ,IDIM,JOB,LLINO,LINCON,LREG,RREG,
|
|
+ TEMP3,
|
|
+ ECHRSZ,MAXGL2,KBOUT,IDEV2,IFAIL,MAXGEL,IDM,PERCD,SEQ3,
|
|
+ ILADD,IRADD,MAXCON,JCON,IWING,NBAD)
|
|
IF (IFAIL.NE.0) GO TO 400
|
|
300 CONTINUE
|
|
CALL YESNO(ANS,'Select another contig',
|
|
+IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
IF(ANS.LT.0) GO TO 400
|
|
IF(ANS.EQ.0) GO TO 150
|
|
ANS = 1
|
|
CALL RADION('Select output format',PROMPT,2,ANS,
|
|
+IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
IF (ANS.LT.1) GO TO 400
|
|
IF(ANS.EQ.1) THEN
|
|
CALL FMTDK(IDEV,SEQ1,IDIM)
|
|
ELSE
|
|
CALL WRITCF(IDEV,SEQ1,IDIM,NAMPRO,KBOUT,IOK)
|
|
END IF
|
|
RETURN
|
|
400 CONTINUE
|
|
CLOSE(UNIT=IDEV)
|
|
END
|
|
C SUBROUTINE TO ENTER NEW GEL SEQUENCES INTO DATA BASE.
|
|
C IT READS IN AN ARCHIVE VERSION AND WRITES OUT A WORKING VERSION.
|
|
C IT ALSO SETS UP ANY RELATIONSHIPS WITH OTHER DATA IN THE DATABASE
|
|
C BOTH BY POSITION IN A CONTIG AND POINTERS TO LEFT AND RIGHT
|
|
C NEIGHBOURS.
|
|
SUBROUTINE AENTER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+GEL,NAMARC,X,ITYPE,ISENSE,SEQC2,ITOTPC,
|
|
+IDIM,IDC,NCONTC,LINCON,IFAIL,IDBSIZ,KBOUT,IDEVR,IDEVW,IDEVN,
|
|
+IDEVT,IDEVC,IDEVG,MAXGEL)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),X,Y
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
CHARACTER GEL(MAXGEL),NAMARC*(*)
|
|
CHARACTER SEQC2(IDC)
|
|
CHARACTER NAMARK*16
|
|
C WRITE(*,*)'X,ITYPE,ISENSE,IDIM,IDC'
|
|
C WRITE(*,*)X,ITYPE,ISENSE,IDIM,IDC
|
|
C SET FAIL FLAG
|
|
IFAIL=0
|
|
C WRITE(KBOUT,1000)
|
|
C1000 FORMAT(' TRYING TO ENTER NEW GEL READING INTO DATABASE')
|
|
C IS THERE SPACE?
|
|
IF((IDBSIZ-(NGELS+NCONTS)).GT.2)GO TO 5
|
|
C FULL
|
|
WRITE(KBOUT,1999)IDBSIZ
|
|
1999 FORMAT(' Database full, current size=',I6,' Extend with copy')
|
|
IFAIL=7
|
|
RETURN
|
|
5 CONTINUE
|
|
C NEED TO CHECK TO SEE IF GEL ALREADY IN DB
|
|
C LOOK THRU ARC FILE
|
|
DO 10 J=1,NGELS
|
|
CALL READN(IDEVN,J,NAMARK)
|
|
IF(NAMARK.NE.NAMARC(1:16))GO TO 10
|
|
C FOUND
|
|
WRITE(KBOUT,1013)J
|
|
1013 FORMAT(' New gel already in database with number',I6,
|
|
+' Entry aborted')
|
|
IFAIL=6
|
|
RETURN
|
|
10 CONTINUE
|
|
C INCREMENT NUMBER OF GELS
|
|
NGELS=NGELS+1
|
|
C SET LENGTH THIS GEL
|
|
LNGTHG(NGELS)=IDIM*ISENSE
|
|
C WRITE NAME OF ARCHIVE TO LIST OF ARCHIVES
|
|
C NAMPRO,ARC
|
|
NAMARK=NAMARC(1:16)
|
|
CALL WRITEN(IDEVN,NGELS,NAMARK)
|
|
WRITE(KBOUT,1003)NGELS
|
|
1003 FORMAT(' This gel reading has been given the number ',I6)
|
|
C WRITE GEL TO WORKING VERSION
|
|
CALL WRITEW(IDEVW,NGELS,GEL,MAXGEL)
|
|
IF(IDEVT.GT.0) CALL ENTRD(IDEVG,IDEVT,IDEVC,NAMARC,NGELS,IOK)
|
|
C CREATE TAGS FOR THIS NASTY
|
|
CALL TAGGEL(NGELS,LNGTHG(NGELS),GEL)
|
|
C SET UP RELATIONSHIPS
|
|
C DOES THIS GEL OVERLAP?
|
|
IF(ITYPE.NE.0)GO TO 100
|
|
C
|
|
C DOES NOT OVERLAP SO IT STARTS A CONTIG OF ITS OWN
|
|
C SET LEFT AND RIGHT POINTERS TO ZERO,RELPG TO 1
|
|
LNBR(NGELS)=0
|
|
RNBR(NGELS)=0
|
|
RELPG(NGELS)=1
|
|
C WRITE NEW GEL LINE
|
|
CALL WRITER(IDEVR,NGELS,RELPG(NGELS),LNGTHG(NGELS),
|
|
+LNBR(NGELS),RNBR(NGELS))
|
|
C
|
|
C SET CONTIG POINTERS AND GENERAL VALUES
|
|
C INCREMENT NUMBER OF CONTIGS
|
|
NCONTS=NCONTS+1
|
|
C POINTER TO THIS CONTIG
|
|
N=IDBSIZ-NCONTS
|
|
C POINTER TO LEFT GEL THIS CONTIG
|
|
LNBR(N)=NGELS
|
|
C POINTER TO RIGHT GEL THIS CONTIG
|
|
RNBR(N)=NGELS
|
|
C LENGTH OF CONTIG
|
|
RELPG(N)=IDIM
|
|
C WRITE CONTIG DESCRIPTOR
|
|
CALL WRITER(IDEVR,N,RELPG(N),LNGTHG(N),
|
|
+LNBR(N),RNBR(N))
|
|
C WRITE DB DESCRIPTOR
|
|
CALL WRITER(IDEVR,IDBSIZ,NGELS,NCONTS,NGELS,NCONTS)
|
|
RETURN
|
|
C
|
|
100 CONTINUE
|
|
C
|
|
C
|
|
C DOES OVERLAP
|
|
150 CONTINUE
|
|
C
|
|
C LEFT END OR RIGHT OVERLAP?
|
|
IF(ITYPE.EQ.1)GO TO 400
|
|
C RIGHT END OR INTERNAL OVERLAP
|
|
C
|
|
160 CONTINUE
|
|
C NEED TO SEARCH THRU THIS CONTIG TO FIND LEFT AND RIGHT
|
|
C NEIGHBOURS FOR THIS NEW GEL
|
|
C LINE NUMBER OF LEFT END OF CONTIG
|
|
N=NCONTC
|
|
C LOOK THRU UNTIL CURRENT IS >= THEN IT MUST BE THE PREVIOUS ONE
|
|
200 CONTINUE
|
|
IF(RELPG(N).GT.X)GO TO 250
|
|
C IS THIS THE LAST GEL IN CONTIG?
|
|
IF(RNBR(N).EQ.0)GO TO 350
|
|
C NO SO LOOK AT NEXT
|
|
N=RNBR(N)
|
|
GO TO 200
|
|
250 CONTINUE
|
|
C GEL LIES BETWEEN N AND LNBR(N)
|
|
C NEED TO EDIT DB HERE
|
|
IF(ITOTPC.GT.0)CALL ABEDIN(RELPG,LNGTHG,LNBR,RNBR,
|
|
1NGELS,NCONTS,
|
|
2GEL,LINCON,X,SEQC2,ITOTPC,IDC,IDBSIZ,KBOUT,IDEVR,IDEVW,
|
|
+MAXGEL)
|
|
C
|
|
C
|
|
C SET POINTERS IN NEW GEL
|
|
LNBR(NGELS)=LNBR(N)
|
|
RNBR(NGELS)=N
|
|
RELPG(NGELS)=X
|
|
C WRITE NEW GEL LINE
|
|
CALL WRITER(IDEVR,NGELS,RELPG(NGELS),LNGTHG(NGELS),
|
|
+LNBR(NGELS),RNBR(NGELS))
|
|
C SET POINTERS IN LEFT AND RIGHT NEIGHBOURS
|
|
K=LNBR(N)
|
|
RNBR(K)=NGELS
|
|
C RNBR(LNBR(N))=NGELS
|
|
C WRITE LEFT AND RIGHT NEIGHBOURS
|
|
CALL WRITER(IDEVR,K,RELPG(K),LNGTHG(K),
|
|
+LNBR(K),RNBR(K))
|
|
LNBR(N)=NGELS
|
|
CALL WRITER(IDEVR,N,RELPG(N),LNGTHG(N),
|
|
+LNBR(N),RNBR(N))
|
|
C WRITE NGELS NCONTS
|
|
CALL WRITER(IDEVR,IDBSIZ,NGELS,NCONTS,NGELS,NCONTS)
|
|
C HAVE WE INCREASED LENGTH OF CONTIG?
|
|
C ITS LINE NUMBER IS LINCON
|
|
C NEED TO UPDATE IDIM IN CASE OF EDITS
|
|
IDIM=ABS(LNGTHG(NGELS))
|
|
Y=X+IDIM-1
|
|
IF(Y.LE.RELPG(LINCON))RETURN
|
|
RELPG(LINCON)=Y
|
|
C WRITE NEW CONTIG LINE
|
|
CALL WRITER(IDEVR,LINCON,RELPG(LINCON),LNGTHG(LINCON),
|
|
+LNBR(LINCON),RNBR(LINCON))
|
|
RETURN
|
|
350 CONTINUE
|
|
C MUST BE A RIGHT END OVERLAP
|
|
C NEED TO EDIT DB HERE
|
|
IF(ITOTPC.GT.0)CALL ABEDIN(RELPG,LNGTHG,LNBR,RNBR,
|
|
1NGELS,NCONTS,
|
|
2GEL,LINCON,X,SEQC2,ITOTPC,IDC,IDBSIZ,KBOUT,IDEVR,IDEVW,
|
|
+MAXGEL)
|
|
C
|
|
C
|
|
C SET POINTERS FOR NEW GEL
|
|
LNBR(NGELS)=N
|
|
RNBR(NGELS)=0
|
|
RELPG(NGELS)=X
|
|
C WRITE NEW GEL LINE
|
|
CALL WRITER(IDEVR,NGELS,RELPG(NGELS),LNGTHG(NGELS),
|
|
+LNBR(NGELS),RNBR(NGELS))
|
|
C OLD RIGHT END
|
|
RNBR(N)=NGELS
|
|
C WRITE NEW RIGHT LINE
|
|
CALL WRITER(IDEVR,N,RELPG(N),LNGTHG(N),
|
|
+LNBR(N),RNBR(N))
|
|
C RESET RIGHT NAME IN CONTIG
|
|
C ITS LINE NUMBER IS LINCON
|
|
RNBR(LINCON)=NGELS
|
|
C HAVE WE INCREASED LENGTH OF CONTIG?
|
|
C NEED TO UPDATE LENGTH OF GEL IN CASE OF EDITS
|
|
IDIM=ABS(LNGTHG(NGELS))
|
|
Y=X+IDIM-1
|
|
RELPG(LINCON)=MAX(RELPG(LINCON),Y)
|
|
C WRITE HERE
|
|
C WRITE CONTIG DESCRIPTOR
|
|
CALL WRITER(IDEVR,LINCON,RELPG(LINCON),LNGTHG(LINCON),
|
|
+LNBR(LINCON),RNBR(LINCON))
|
|
CALL WRITER(IDEVR,IDBSIZ,NGELS,NCONTS,NGELS,NCONTS)
|
|
RETURN
|
|
C
|
|
400 CONTINUE
|
|
C
|
|
C ADDING TO LEFT END
|
|
410 CONTINUE
|
|
C NEED TO EDIT DB HERE
|
|
IF(ITOTPC.GT.0)CALL ABEDIN(RELPG,LNGTHG,LNBR,RNBR,
|
|
1NGELS,NCONTS,
|
|
2GEL,LINCON,1,SEQC2,ITOTPC,IDC,IDBSIZ,KBOUT,IDEVR,IDEVW,
|
|
+MAXGEL)
|
|
C
|
|
420 CONTINUE
|
|
C SET POINTERS IN NEW GEL
|
|
RELPG(NGELS)=1
|
|
RNBR(NGELS)=NCONTC
|
|
LNBR(NGELS)=0
|
|
C WRITE NEW GEL LINE
|
|
CALL WRITER(IDEVR,NGELS,RELPG(NGELS),LNGTHG(NGELS),
|
|
+LNBR(NGELS),RNBR(NGELS))
|
|
C SET POINTERS IN OLD LEFT END
|
|
LNBR(NCONTC)=NGELS
|
|
RELPG(NCONTC)=X
|
|
C WRITE NEW LEFT END
|
|
CALL WRITER(IDEVR,NCONTC,RELPG(NCONTC),LNGTHG(NCONTC),
|
|
+LNBR(NCONTC),RNBR(NCONTC))
|
|
C NEW LENGTH OF CONTIG
|
|
RELPG(LINCON)=RELPG(LINCON)+X-1
|
|
C MAY HAVE JUST ADDED A GEL LONGER THAN CONTIG
|
|
IDIM=ABS(LNGTHG(NGELS))
|
|
Y=IDIM
|
|
IF(Y.GT.RELPG(LINCON))RELPG(LINCON)=Y
|
|
C NEW NAME OF LEFT END OF CONTIG
|
|
LNBR(LINCON)=NGELS
|
|
C WRITE CONTIG DESCRIPTOR
|
|
CALL WRITER(IDEVR,LINCON,RELPG(LINCON),LNGTHG(LINCON),
|
|
+LNBR(LINCON),RNBR(LINCON))
|
|
CALL WRITER(IDEVR,IDBSIZ,NGELS,NCONTS,NGELS,NCONTS)
|
|
C NOW GO THRU AND CHANGE ALL RELATIVE POSITIONS
|
|
N=NCONTC
|
|
440 CONTINUE
|
|
IF(RNBR(N).EQ.0)RETURN
|
|
N=RNBR(N)
|
|
RELPG(N)=RELPG(N)+X-1
|
|
C WRITE NEW LINE
|
|
CALL WRITER(IDEVR,N,RELPG(N),LNGTHG(N),
|
|
+LNBR(N),RNBR(N))
|
|
GO TO 440
|
|
END
|
|
SUBROUTINE COPYTG(IDEVT,IDEV,IOK,IDBSIZ,NEWSIZ,NGELS)
|
|
C Read tag details
|
|
IDIFF = NEWSIZ - IDBSIZ
|
|
CALL READTG(IDEVT,IDBSIZ,ICNT,LLEN,LCOM,LTYPE,NEXT)
|
|
IF (NEXT.NE.0) NEXT = NEXT + IDIFF
|
|
CALL WRITTG(IDEV,NEWSIZ,ICNT+IDIFF,LLEN,LCOM,LTYPE,NEXT)
|
|
C Copy headers for each gels
|
|
DO 10 I = 1,NGELS
|
|
CALL READTG(IDEVT,I,LPOS,LLEN,LCOM,LTYPE,NEXT)
|
|
IF (NEXT.NE.0) NEXT = NEXT + IDIFF
|
|
CALL WRITTG(IDEV,I,LPOS,LLEN,LCOM,LTYPE,NEXT)
|
|
10 CONTINUE
|
|
C Copy rest of tags
|
|
DO 20 I = IDBSIZ+1, ICNT
|
|
CALL READTG(IDEVT,I,LPOS,LLEN,LCOM,LTYPE,NEXT)
|
|
IF (NEXT.NE.0) NEXT = NEXT + IDIFF
|
|
CALL WRITTG(IDEV,I+IDIFF,LPOS,LLEN,LCOM,LTYPE,NEXT)
|
|
20 CONTINUE
|
|
IOK = 0
|
|
END
|
|
SUBROUTINE COPYCC(IDEVC,IDEV,IOK)
|
|
C COMMENT_LENGTH
|
|
CHARACTER NOTE*40
|
|
CALL READCC(IDEVC,1,ICNT,NEXT,NOTE)
|
|
CALL WRITCC(IDEV,1,ICNT,NEXT,NOTE)
|
|
DO 10 I = 2,ICNT
|
|
CALL READCC(IDEVC,I,ICNT,NEXT,NOTE)
|
|
CALL WRITCC(IDEV,I,ICNT,NEXT,NOTE)
|
|
10 CONTINUE
|
|
IOK = 0
|
|
END
|
|
SUBROUTINE DBAUTO(RELPG,LNGTHG,LNBR,RNBR,MAXDB,IDBSIZ,
|
|
+NGELS,NCONTS,MAXGEL,
|
|
+TEMP3,WORDP,WORDN,LPOWRC,POSNS,GELN,
|
|
+SEQ1,MAXSEQ,SEQ2,SEQ3,SEQ4,SEQ5,SEQC2,SEQG2,MATCH,
|
|
+MAXGLM,MAXGL2,CHRSIZ,ECHRSZ,LENGTH,
|
|
+SAV1,SAV2,SAV3,MAXSAV,CENDS,NENDS,MAXCON,CONST,
|
|
+KBIN,KBOUT,IDEV1,IDEV2,IDEV3,IDEV4,IDEV7,IDEV8,IDEV,IDEVT,IDEVC,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,NAMARC,NAMPRO,FILE,
|
|
+PERCD,IOPEN,IDM,SEQG3,SEQC3,IOK)
|
|
INTEGER CHRSIZ,ECHRSZ
|
|
INTEGER RELPG(MAXDB),PL(2),PR(2),RMOST
|
|
INTEGER LNGTHG(MAXDB),LNBR(MAXDB),RNBR(MAXDB)
|
|
INTEGER JOINT(2),ITOTPC(2),ITOTPG(2),IDIM22(2),IDOUT(2)
|
|
INTEGER LINCON(2),LLINO(2),ITYPE(2),IFAIL(2)
|
|
INTEGER ILEFTS(2),ILC(2),IPOSC(2),IPOSG(2),ISENSE(2)
|
|
INTEGER LREG,RREG,X,ANSJOK
|
|
INTEGER TEMP3(ECHRSZ,MAXGL2),CONST(LENGTH)
|
|
INTEGER POSNS(MAXSEQ),WORDP(LPOWRC),WORDN(LPOWRC),GELN(MAXGLM)
|
|
INTEGER CENDS(MAXCON),NENDS(MAXCON)
|
|
CHARACTER SEQ3(MAXGLM),SEQC2(MAXGLM,2),SEQG2(MAXGLM,2)
|
|
CHARACTER SEQ1(MAXSEQ),SEQ2(MAXGLM),MATCH(MAXGLM),SEQ4(MAXGLM)
|
|
INTEGER SAV1(MAXSAV),SAV2(MAXSAV),SAV3(MAXSAV),CSTART
|
|
CHARACTER NAMARC*(*),NAMPRO*(*),FILE*(*)
|
|
CHARACTER GET,SEQ5(MAXGLM),HELPF*(*),SEQG3(MAXGLM),SEQC3(MAXGLM)
|
|
PARAMETER (MAXPRM = 32)
|
|
CHARACTER PROMPT(3)*(MAXPRM)
|
|
REAL PERMIS(2)
|
|
INTEGER GNFFOF
|
|
EXTERNAL GNFFOF
|
|
SAVE GET
|
|
DATA GET/'>'/
|
|
WRITE(KBOUT,*)' Automatic sequence assembler'
|
|
C
|
|
C set flag for saving alignment scores
|
|
C
|
|
IREPSC = 1
|
|
ICRAP = 1
|
|
IFAIL(1) = 0
|
|
IEMPTY=0
|
|
IF(NGELS.LT.1)IEMPTY=1
|
|
CALL DBCHEK(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+TEMP3,IERR,KBOUT)
|
|
IF(IERR.GT.1) RETURN
|
|
CALL YESNO(IOKENT,'Permit entry',
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(IOKENT.LT.0) RETURN
|
|
CALL YESNO(ISHOW,'Hide alignments',
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(ISHOW.LT.0) RETURN
|
|
CALL YESNO(INF,'Use file of file names',
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(INF.LT.0) RETURN
|
|
IF(INF.EQ.0) THEN
|
|
FILE = ' '
|
|
CALL OPENF1(IDEV7,FILE,0,IOK,KBIN,KBOUT,
|
|
+ 'File of gel reading names',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH)
|
|
IF(IOK.NE.0) RETURN
|
|
END IF
|
|
IF(IOKENT.EQ.0) THEN
|
|
FILE = ' '
|
|
CALL OPENF1(IDEV8,FILE,1,IOK,KBIN,KBOUT,
|
|
+ 'File for names of failures',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH)
|
|
IF(IOK.NE.0) RETURN
|
|
ELSE
|
|
CALL YESNO(IREPSC,'Save alignment scores in a file',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(IREPSC.LT.0) RETURN
|
|
IF(IREPSC.EQ.0) THEN
|
|
FILE = ' '
|
|
CALL OPENF1(IDEV8,FILE,1,IOK,KBIN,KBOUT,
|
|
+ 'File for names and scores',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH)
|
|
IF(IOK.NE.0) RETURN
|
|
CALL YESNO(ICRAP,'Use poor data',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(ICRAP.LT.0) RETURN
|
|
END IF
|
|
END IF
|
|
PROMPT(1) = 'Perform normal shotgun assembly'
|
|
PROMPT(2) = 'Put all sequences in one contig'
|
|
PROMPT(3) = 'Put all sequences in new contigs'
|
|
IOPT = 1
|
|
CALL RADION('Select entry mode',PROMPT,3,IOPT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(IOPT.LT.1) RETURN
|
|
IF(IOPT.EQ.1) THEN
|
|
C parameters for normal assembly
|
|
ANSJOK = 0
|
|
CALL YESNO(ANSJOK,'Permit joins',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(ANSJOK.LT.0) RETURN
|
|
MN = LENGTH*2
|
|
MX = MAXGLM + 1
|
|
MINMAT = MAX(15,MN)
|
|
CALL GETINT(MN,MX,MINMAT,
|
|
+ 'Minimum initial match',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
MINMAT = IVAL
|
|
MINSLI = 3
|
|
MN = 0
|
|
MX = 25
|
|
MAXPG = 8
|
|
CALL GETINT(MN,MX,MAXPG,
|
|
+ 'Maximum pads per gel',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
MAXPG = IVAL
|
|
MN = 0
|
|
MX = 25
|
|
MAXPC = 8
|
|
CALL GETINT(MN,MX,MAXPC,
|
|
+ 'Maximum pads per gel in contig',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
MAXPC = IVAL
|
|
C IF((IOKENT.EQ.0).OR.(IREPSC.EQ.0)) THEN
|
|
RMN = 0.
|
|
RMX = 100.
|
|
PERMAX = 8.
|
|
CALL GETRL(RMN,RMX,PERMAX,
|
|
+ 'Maximum percent mismatch after alignment',
|
|
+ VAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
PERMAX = VAL
|
|
C END IF
|
|
IDIM1=0
|
|
MAXOVR=MAXGEL-3*MAX(MAXPC,MAXPG)
|
|
JGEL = 0
|
|
JNGEL = 0
|
|
JNJOIN = 0
|
|
JOINF = 0
|
|
IMATC = 0
|
|
IF(IEMPTY.EQ.0) THEN
|
|
JOB = 0
|
|
CALL ACONSN(RELPG,LNGTHG,LNBR,RNBR,NAMPRO,NGELS,NCONTS,
|
|
+ SEQ1,MAXSEQ,SEQ2,IDBSIZ,IDIM1,JOB,KDUMM,KDUMM,KDUMM,TEMP3,
|
|
+ ECHRSZ,MAXGL2,KBOUT,IDEV2,IFAIL(1),MAXGEL,IDM,PERCD)
|
|
IF(IFAIL(1).NE.0) THEN
|
|
CALL ERROM(KBOUT,'Error calculating consensus')
|
|
GO TO 900
|
|
END IF
|
|
END IF
|
|
C
|
|
C init hashing constants
|
|
C
|
|
CALL INITE(CONST,CSTART,LENGTH)
|
|
END IF
|
|
C
|
|
C set intitial values for contig count and number of readings
|
|
C just to get thru the first consensus calc
|
|
NGELSL = NGELS + 2
|
|
NCONTL = NCONTS +1
|
|
C
|
|
C
|
|
C MAIN LOOP
|
|
C
|
|
C
|
|
1 CONTINUE
|
|
C IF(IDIM1.GT.0)CALL FMTDB(SEQ1,IDIM1,1,IDIM1,60,IDEV)
|
|
C
|
|
C
|
|
IDIM2=MAXGEL
|
|
IF(INF.EQ.1) THEN
|
|
3 CONTINUE
|
|
MN = 0
|
|
CALL GTSTR('Gel reading name',' ',NAMARC,MN,KBOUT,KBIN,INFLAG)
|
|
IF(INFLAG.EQ.3) RETURN
|
|
IF(INFLAG.EQ.2) RETURN
|
|
IF(INFLAG.EQ.1) THEN
|
|
CALL HELP2(IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
GO TO 3
|
|
END IF
|
|
ELSE
|
|
IOK = GNFFOF(IDEV7,NAMARC)
|
|
IF(IOK.EQ.1) GO TO 900
|
|
IF(IOK.NE.0) GO TO 1
|
|
END IF
|
|
DO 77 MM=1,80
|
|
MATCH(MM)=GET
|
|
77 CONTINUE
|
|
WRITE(IDEV,1077)(MATCH(KK),KK=1,79)
|
|
1077 FORMAT(' ',79A1)
|
|
JGEL = JGEL + 1
|
|
WRITE(IDEV,*)'Processing',JGEL,' in batch'
|
|
1007 FORMAT(' Gel reading name=',A)
|
|
WRITE(IDEV,1007)NAMARC
|
|
CALL OPENRS(IDEV4,NAMARC,IOK,LRECL,2)
|
|
IF(IOK.NE.0)THEN
|
|
C IF(INF.EQ.1) RETURN
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,0)
|
|
GO TO 1
|
|
END IF
|
|
CALL ARRFIN(IDEV4,SEQ2,IDIM2,KBOUT,ICRAP)
|
|
CLOSE(UNIT=IDEV4)
|
|
WRITE(IDEV,1800)IDIM2
|
|
1800 FORMAT(' Gel reading length=',I6)
|
|
C
|
|
C
|
|
C
|
|
IF(IOPT.NE.1) THEN
|
|
CALL DBAUTP(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ SEQ2,NAMARC,JOINT,ITYPE,ISENSE,SEQC2,ITOTPC,
|
|
+ IDIM2,IDOUT,LLINO,LINCON,IFAIL,IDBSIZ,MAXDB,IDEV,
|
|
+ IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,IDEV4,MAXGEL,IMATC,IEMPTY,IOPT)
|
|
IF(IFAIL(1).NE.0) THEN
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,0)
|
|
ELSE
|
|
JNGEL = JNGEL + 1
|
|
END IF
|
|
GO TO 1
|
|
END IF
|
|
IF(IDIM2.LT.MINMAT)THEN
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,1)
|
|
GO TO 1
|
|
END IF
|
|
CALL SQCOPY(SEQ2,SEQ3,IDIM2)
|
|
IFCOMP=0
|
|
IMATC=0
|
|
JOBC = 2
|
|
IF ((NGELSL.LT.NGELS).AND.(NCONTL.LT.NCONTS)) THEN
|
|
JOBC = 1
|
|
ELSE IF (NGELSL.EQ.NGELS) THEN
|
|
JOBC = 0
|
|
END IF
|
|
NGELSL = NGELS
|
|
NCONTL = NCONTS
|
|
IF(IEMPTY.EQ.0)
|
|
+CALL AUTOCN(SEQ1,IDIM1,SEQ2,IDIM2,ILEFTS,ILC,IPOSC,
|
|
+IPOSG,ISENSE,LLINO,IMATC,IFCOMP,MINMAT,POSNS,WORDP,WORDN,
|
|
+CONST,LENGTH,LPOWRC,IDEV,MATCH,MAXGEL,MAXGLM,SEQ5,GELN,
|
|
+SAV1,SAV2,SAV3,MAXSAV,CENDS,NENDS,MAXCON,
|
|
+SEQG2,SEQC2,SEQ4,IDOUT,IDIM22,ITOTPG,ITOTPC,JOINT,IFAIL,
|
|
+ITYPE,MAXPC,MAXPG,PERMAX,MINSLI,SEQG3,SEQC3,KFAIL,CSTART,
|
|
+JOBC,PERMIS,LENO,ISHOW)
|
|
IF(IREPSC.EQ.0) THEN
|
|
IF(IFCOMP.NE.0) THEN
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,2)
|
|
GO TO 1
|
|
END IF
|
|
IF(IMATC.GT.0) THEN
|
|
WRITE(IDEV8,1022)NAMARC,PERMIS(1),IDIM2,LENO
|
|
ELSE
|
|
PERMIS(1) = 0.
|
|
LENO = 0
|
|
WRITE(IDEV8,1022)NAMARC,PERMIS(1),IDIM2,LENO
|
|
END IF
|
|
END IF
|
|
1022 FORMAT(' ',A,F5.1,2I6)
|
|
IF(IOKENT.NE.0) GO TO 1
|
|
C THIS RETURNS THE FOLLOWING:
|
|
C ILEFTS POSITION IN CONSENSUS OF LEFT END OF MATCHING CONTIGS
|
|
C ILC LENGTHS OF MATCHING CONTIGS
|
|
C IPOSC POSITION OF MATCH RELATIVE TO CONTIG
|
|
C IPOSG POSITION OF MATCH RELATIVE TO NEW GEL
|
|
C ISENSE SENSE OF NEW GEL
|
|
C LLINO LEFT GEL NUMBER IN MATCHING CONTIGS
|
|
C IMATC THE NUMBER OF MATCHING CONTIGS (>2 IS ERROR!)
|
|
C IFCOMP ERROR FLAG FOR COMPARISON (COMPARISON ARRAYS OVERFLOWED)
|
|
IF(IFCOMP.NE.0) THEN
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,2)
|
|
GO TO 1
|
|
END IF
|
|
CALL SQCOPY(SEQ3,SEQ2,IDIM2)
|
|
IF(IMATC.EQ.0) THEN
|
|
C
|
|
C NO OVERLAP NEW CONTIG
|
|
C
|
|
C ITYPE 0 = NO OVERLAP
|
|
C ISENSE 1 = SAME SENSE AS ARCHIVE
|
|
IF(IFAIL(1).NE.0) THEN
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,2)
|
|
GO TO 1
|
|
END IF
|
|
ITYPE(1)=0
|
|
ISENSE(1)=1
|
|
IDOUT(1)=MAXGEL
|
|
WRITE(IDEV,1015)
|
|
1015 FORMAT(' New gel reading does not overlap: start a new contig')
|
|
CALL AENTER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ SEQ2,NAMARC,X,ITYPE,ISENSE,SEQC2(1,1),ITOTPC(1),
|
|
+ IDIM2,IDOUT(1),LLINO,LINCON,IFAIL,IDBSIZ,IDEV,
|
|
+ IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,IDEV4,MAXGEL)
|
|
IF(IFAIL(1).NE.0) THEN
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,3)
|
|
GO TO 1
|
|
END IF
|
|
IEMPTY=0
|
|
IDIM1=IDIM1+1
|
|
IF((IDIM1+19+IDIM2).GT.MAXSEQ)THEN
|
|
WRITE(IDEV,1021)MAXSEQ
|
|
1021 FORMAT(' Database maximum consensus length (',I6,') exceeded')
|
|
GO TO 900
|
|
END IF
|
|
CALL ADDTIT(SEQ1(IDIM1),NAMPRO,NGELS,IDIM1)
|
|
CALL MSTLKL(SEQ2,IDIM2)
|
|
CALL SQCOPY(SEQ2,SEQ1(IDIM1),IDIM2)
|
|
IDIM1=IDIM1+IDIM2-1
|
|
JNGEL = JNGEL + 1
|
|
GO TO 1
|
|
END IF
|
|
C
|
|
C
|
|
C OVERLAP SO TRY TO ALIGN THE SEQUENCES
|
|
C
|
|
C
|
|
DO 100 I=1,IMATC
|
|
N=IDBSIZ-NCONTS
|
|
DO 99 J=N,IDBSIZ-1
|
|
IF(LNBR(J).NE.LLINO(I))GO TO 99
|
|
LINCON(I)=J
|
|
GO TO 100
|
|
99 CONTINUE
|
|
WRITE(IDEV,10077)LLINO(I)
|
|
10077 FORMAT(' Contig line for contig',I6,' not found!')
|
|
GO TO 900
|
|
100 CONTINUE
|
|
C
|
|
IF (IMATC.EQ.1) THEN
|
|
C
|
|
C
|
|
C SINGLE OVERLAP
|
|
C
|
|
C
|
|
C
|
|
WRITE(IDEV,1014)LLINO(1)
|
|
1014 FORMAT(' New gel reading overlaps contig',I6)
|
|
IF(ITOTPG(1).GT.0) CALL CCTA(SEQG2(1,1),IDIM22(1))
|
|
CALL AENTER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ SEQG2(1,1),NAMARC,JOINT(1),ITYPE(1),ISENSE(1),
|
|
+ SEQC2(1,1),
|
|
+ ITOTPC(1),IDIM22(1),IDOUT(1),LLINO(1),LINCON(1),
|
|
+ IFAIL(1),IDBSIZ,IDEV,IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,IDEV4,MAXGEL)
|
|
IF(IFAIL(1).NE.0) THEN
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,3)
|
|
GO TO 1
|
|
END IF
|
|
JNGEL = JNGEL + 1
|
|
CALL UPDCON(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,NGELS,NCONTS,
|
|
+ SEQ1,MAXSEQ,IDIM1,ILEFTS(1),ILC(1),LINCON(1),NAMPRO,SEQ2,TEMP3,
|
|
+ ECHRSZ,MAXGL2,KBOUT,IDEV,IDEV2,IFAIL(1),MAXGEL,IDM,PERCD)
|
|
IF(IFAIL(1).NE.0) THEN
|
|
CALL ERROM(KBOUT,'Error calculating consensus')
|
|
GO TO 900
|
|
END IF
|
|
IF(KFAIL.NE.0) THEN
|
|
C CALL AERROR(IDEV,IDEV8,NAMARC,4)
|
|
END IF
|
|
GO TO 1
|
|
END IF
|
|
C
|
|
C
|
|
C DOUBLE OVERLAP
|
|
C
|
|
C
|
|
WRITE(IDEV,1013)LLINO
|
|
1013 FORMAT(' Overlap between contigs',I6,' and',I6)
|
|
IF(ANSJOK.NE.0) GO TO 1
|
|
IF(LLINO(1).EQ.LLINO(2))THEN
|
|
WRITE(IDEV,*)' Trying to form loop in contig',LLINO(1)
|
|
WRITE(IDEV,*)' Gel not entered'
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,3)
|
|
GO TO 1
|
|
END IF
|
|
CALL AJOIN3(RELPG,IDBSIZ,LINCON,ITYPE,ISENSE,JOINT,
|
|
+IDIM22,KLASS,IOVER,IDEV,PL,PR)
|
|
IF(IOVER.GT.MAXOVR)THEN
|
|
WRITE(IDEV,*)' Overlap too large: entry only'
|
|
C
|
|
C cannot align the two contigs, so try to enter the reading into one of them
|
|
C
|
|
IFAIL(2)=1
|
|
IGOOD=0
|
|
IF(IFAIL(1).EQ.0)IGOOD=1
|
|
IF(IFAIL(2).EQ.0)IGOOD=2
|
|
IF(IGOOD.EQ.0) THEN
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,2)
|
|
JOINF = JOINF + 1
|
|
GO TO 1
|
|
END IF
|
|
IF(ITOTPG(IGOOD).GT.0) CALL CCTA(SEQG2(1,IGOOD),IDIM22(IGOOD))
|
|
WRITE(IDEV,1012)LLINO(IGOOD)
|
|
CALL AENTER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ SEQG2(1,IGOOD),NAMARC,JOINT(IGOOD),ITYPE(IGOOD),
|
|
+ ISENSE(IGOOD),SEQC2(1,IGOOD),ITOTPC(IGOOD),
|
|
+ IDIM22(IGOOD),IDOUT(IGOOD),LLINO(IGOOD),LINCON(IGOOD),
|
|
+ IFAIL(IGOOD),IDBSIZ,IDEV,IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,IDEV4,
|
|
+ MAXGEL)
|
|
IF(IFAIL(IGOOD).NE.0) THEN
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,3)
|
|
JOINF = JOINF + 1
|
|
GO TO 1
|
|
END IF
|
|
JNGEL = JNGEL + 1
|
|
CALL UPDCON(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,NGELS,NCONTS,
|
|
+ SEQ1,MAXSEQ,IDIM1,ILEFTS(IGOOD),ILC(IGOOD),LINCON(IGOOD),
|
|
+ NAMPRO,SEQ2,TEMP3,
|
|
+ ECHRSZ,MAXGL2,KBOUT,IDEV,IDEV2,IFAIL(1),MAXGEL,IDM,PERCD)
|
|
IF(IFAIL(1).NE.0) THEN
|
|
CALL ERROM(KBOUT,'Error calculating consensus')
|
|
GO TO 900
|
|
END IF
|
|
WRITE(IDEV,1020)LLINO,LLINO(IGOOD)
|
|
1020 FORMAT(' Could not join contigs',I4,' and',I4,' but the gel',
|
|
+ ' has been entered into contig',I4,/,
|
|
+ ' If required do the join manually.')
|
|
JOINF = JOINF + 1
|
|
GO TO 1
|
|
END IF
|
|
C WHICH CONTIG IS LEFTMOST?
|
|
LMOST=1
|
|
RMOST=2
|
|
IF(PL(1).GT.PL(2))THEN
|
|
LMOST=2
|
|
RMOST=1
|
|
END IF
|
|
C SAVE LENGTH OF RMOST CONTIG FOR DELETION STEP LATER
|
|
ILCR=ILC(RMOST)
|
|
IF(ITOTPG(LMOST).GT.0) CALL CCTA(SEQG2(1,LMOST),IDIM22(LMOST))
|
|
WRITE(IDEV,1012)LLINO(LMOST)
|
|
1012 FORMAT(' Entering the new gel reading into contig',I6)
|
|
CALL AENTER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+SEQG2(1,LMOST),NAMARC,JOINT(LMOST),ITYPE(LMOST),
|
|
+ISENSE(LMOST),SEQC2(1,LMOST),ITOTPC(LMOST),
|
|
+IDIM22(LMOST),IDOUT(LMOST),LLINO(LMOST),LINCON(LMOST),
|
|
+IFAIL(LMOST),IDBSIZ,IDEV,IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,IDEV4,
|
|
+MAXGEL)
|
|
IF(IFAIL(LMOST).NE.0) THEN
|
|
CALL AERROR(IDEV,IDEV8,NAMARC,3)
|
|
JOINF = JOINF + 1
|
|
GO TO 1
|
|
END IF
|
|
JNGEL = JNGEL + 1
|
|
CALL UPDCON(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,NGELS,NCONTS,
|
|
+ SEQ1,MAXSEQ,IDIM1,ILEFTS(LMOST),ILC(LMOST),LINCON(LMOST),
|
|
+ NAMPRO,SEQ2,TEMP3,
|
|
+ ECHRSZ,MAXGL2,KBOUT,IDEV,IDEV2,IFAIL(1),MAXGEL,IDM,PERCD)
|
|
IF(IFAIL(1).NE.0) THEN
|
|
CALL ERROM(KBOUT,'Error calculating consensus')
|
|
GO TO 900
|
|
END IF
|
|
IF(ITYPE(LMOST).EQ.1)LLINO(LMOST)=NGELS
|
|
IF(ILEFTS(LMOST).LT.ILEFTS(RMOST))THEN
|
|
ILEFTS(RMOST)=ILEFTS(RMOST) +
|
|
+ (RELPG(LINCON(LMOST)) - ILC(LMOST))
|
|
END IF
|
|
ILC(LMOST) = RELPG(LINCON(LMOST))
|
|
DO 500 I=1,2
|
|
IF(ISENSE(I).EQ.-1)THEN
|
|
CALL CMPLMT(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,LINCON(I),
|
|
+ LLINO(I),SEQ2,IDBSIZ,IDEV,IDEV1,IDEV2,MAXGEL)
|
|
CALL SQREV(SEQ1(ILEFTS(I)),ILC(I))
|
|
CALL SQCOM(SEQ1(ILEFTS(I)),ILC(I))
|
|
KT=IDIM1
|
|
CALL ADDTIT(SEQ1((ILEFTS(I)-20)),NAMPRO,LNBR(LINCON(I)),KT)
|
|
END IF
|
|
500 CONTINUE
|
|
C NEED TO KNOW POSITION OF OVERLAP RELATIVE TO CONTIG, TO CONSENSUS
|
|
C WHICH BITS TO SEND TO ALIGNMENT ROUTINES
|
|
C SET UP FOR ALINE (NOTE RMOST IS EQUIVALENT TO THE GEL READING AND
|
|
C SO IS SLID ALONG THE LMOST CONTIG. THE SECTION SENT TO ALINE MUST
|
|
C BE OF LENGTH < MAXGEL-2*MAX(MAXPC,MAXPG)
|
|
C IT MUST START AT POSITION 1 IN THE RMOST CONTIG AND EXTEND
|
|
IPOSC(LMOST)=PL(RMOST)+RELPG(NGELS)-1
|
|
ILCT = RELPG(LINCON(LMOST)) - RELPG(NGELS) - PL(RMOST) + 2
|
|
ILC(RMOST)=MIN(ILCT,ILC(RMOST))
|
|
IPOSC(RMOST)=1
|
|
IDOUT(LMOST)=MAXGEL
|
|
IDOUT(RMOST)=MAXGEL
|
|
IDSAV=MAXSAV
|
|
C ON INPUT TO ALINE ILC(RMOST) CONTAINS THE OVERLAP LENGTH
|
|
C ON OUTPUT IT CONTAINS THE LENGTH OF THE ALIGNED SECTION (IE INCLUDING
|
|
C PADS)
|
|
WRITE(IDEV,1009)
|
|
1009 FORMAT(' Trying to align the two contigs')
|
|
CALL ALINE(SEQ1(ILEFTS(LMOST)),SEQ1(ILEFTS(RMOST)),
|
|
+SEQC2(1,RMOST),SEQC2(1,LMOST),SAV1,SAV2,SAV3,IDSAV,
|
|
+ILC(LMOST),ILC(RMOST),IDOUT(LMOST),IPOSC(LMOST),IPOSC(RMOST),
|
|
+MINSLI,JOINT(LMOST),ITOTPC(LMOST),ITOTPC(RMOST),IFAIL(1),
|
|
+ITYPE(1),MAXPC,MAXPC,PERMAX,IDEV,SEQ4,MAXGEL,Z,LENO,ISHOW)
|
|
C SEQC2(1,LMOST) NOW CONTAINS THE ALIGNED SECTION OF THE LMOST CONTIG
|
|
C SEQC2(1,RMOST) NOW CONTAINS THE ALIGNED SECTION OF THE RMOST CONTIG
|
|
C ILC(RMOST) IS NOW THE LENGTH OF ALIGNED SECTION OF THE RMOST CONTIG
|
|
C IDOUT(LMOST) IS NOW THE LENGTH OF ALIGNED SECTION OF THE LMOST CONTIG
|
|
C JOINT(LMOST) IS THE POSITION OF THE JOIN RLETIVE TO THE LMOST CONTIG
|
|
C ITYPE IS TYPE OF OVERLAP (-1 = RIGHT END OR INTERNAL, 1 = LEFT END)
|
|
C NB SHOULD ALWAYS BE -1
|
|
C IF THIS HAS BEEN DONE OK WE CAN EDIT THE TWO CONTIGS THEN JOIN
|
|
IF(IFAIL(1).NE.0)THEN
|
|
WRITE(IDEV,*)' Failed to align the two overlapping contigs'
|
|
C CALL AERROR(IDEV,IDEV8,NAMARC,4)
|
|
JOINF = JOINF + 1
|
|
GO TO 1
|
|
END IF
|
|
IF(ITOTPC(LMOST).GT.0)THEN
|
|
WRITE(IDEV,1017)LLINO(LMOST)
|
|
1017 FORMAT(' Editing contig',I6)
|
|
CALL ABEDIN(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ NGELS,NCONTS,SEQ3,LINCON(LMOST),JOINT(LMOST),SEQC2(1,LMOST),
|
|
+ ITOTPC(LMOST),IDOUT(LMOST),IDBSIZ,IDEV,IDEV1,IDEV2,
|
|
+ MAXGEL)
|
|
END IF
|
|
JOINT(RMOST)=1
|
|
IDOUT(RMOST)=ILC(RMOST)
|
|
IF(ITOTPC(RMOST).GT.0)THEN
|
|
WRITE(IDEV,1017)LLINO(RMOST)
|
|
CALL ABEDIN(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ NGELS,NCONTS,SEQ3,LINCON(RMOST),JOINT(RMOST),SEQC2(1,RMOST),
|
|
+ ITOTPC(RMOST),IDOUT(RMOST),IDBSIZ,IDEV,IDEV1,IDEV2,
|
|
+ MAXGEL)
|
|
END IF
|
|
ILC(RMOST)=ILCR
|
|
LTL=LNBR(LINCON(LMOST))
|
|
LTR=LNBR(LINCON(RMOST))
|
|
WRITE(IDEV,1018)LNBR(LINCON(LMOST)),LNBR(LINCON(RMOST))
|
|
1018 FORMAT(' Completing the join between contigs',I6,' and',I6)
|
|
CALL AJOIN2(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+JOINT(LMOST),LTL,LTR,LINCON(LMOST),LINCON(RMOST),IDEV1)
|
|
LLINO(1)=LTL
|
|
IF(ILEFTS(LMOST).GT.ILEFTS(RMOST))THEN
|
|
CALL DELCON(SEQ1,ILEFTS(LMOST),ILC(LMOST),IDIM1)
|
|
CALL DELCON(SEQ1,ILEFTS(RMOST),ILC(RMOST),IDIM1)
|
|
END IF
|
|
IF(ILEFTS(RMOST).GE.ILEFTS(LMOST))THEN
|
|
CALL DELCON(SEQ1,ILEFTS(RMOST),ILC(RMOST),IDIM1)
|
|
CALL DELCON(SEQ1,ILEFTS(LMOST),ILC(LMOST),IDIM1)
|
|
END IF
|
|
LREG=1
|
|
RREG=JOINT(LMOST)
|
|
IGELC=LLINO(1)
|
|
JOB = 1
|
|
CALL ACONSN(RELPG,LNGTHG,LNBR,RNBR,NAMPRO,NGELS,NCONTS,
|
|
+SEQ1,MAXSEQ,SEQ2,IDBSIZ,IDIM1,JOB,IGELC,LREG,RREG,TEMP3,
|
|
+ECHRSZ,MAXGL2,IDEV,IDEV2,IFAIL(1),MAXGEL,IDM,PERCD)
|
|
IF(IFAIL(1).NE.0) THEN
|
|
CALL ERROM(KBOUT,'Error calculating consensus')
|
|
GO TO 900
|
|
END IF
|
|
C CALL FMTDB(SEQ1,IDIM1,1,IDIM1,60,IDEV)
|
|
JNJOIN = JNJOIN + 1
|
|
IF(KFAIL.NE.0) THEN
|
|
C CALL AERROR(IDEV,IDEV8,NAMARC,4)
|
|
C JOINF = JOINF + 1
|
|
END IF
|
|
GO TO 1
|
|
900 CONTINUE
|
|
C CALL FMTDB(SEQ1,IDIM1,1,IDIM1,60,IDEV)
|
|
WRITE(KBOUT,*)'Batch finished'
|
|
WRITE(KBOUT,*)JGEL,' sequences processed'
|
|
WRITE(KBOUT,*)JNGEL,' sequences entered into database'
|
|
WRITE(KBOUT,*)JNJOIN,' joins made'
|
|
WRITE(KBOUT,*)JOINF, ' joins failed'
|
|
END
|
|
SUBROUTINE DBAUTP(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+SEQ2,NAMARC,JOINT,ITYPE,ISENSE,SEQC2,ITOTPC,
|
|
+IDIM2,IDOUT,LLINO,LINCON,IFAIL,IDBSIZ,MAXDB,IDEV,
|
|
+IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,IDEV4,MAXGEL,IMATC,IEMPTY,IOPT)
|
|
INTEGER RELPG(MAXDB)
|
|
INTEGER LNGTHG(MAXDB),LNBR(MAXDB),RNBR(MAXDB)
|
|
CHARACTER SEQ2(MAXGEL),SEQC2(MAXGEL)
|
|
CHARACTER NAMARC*(*)
|
|
C deals with entering all readings into contig 1 (IOPT=2)
|
|
C or all readings into new contigs (IOPT=3)
|
|
IF(IOPT.EQ.2) THEN
|
|
IF(IMATC.EQ.0) THEN
|
|
ITYPE=0
|
|
ISENSE=1
|
|
IDOUT=MAXGEL
|
|
CALL AENTER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ SEQ2,NAMARC,JOINT,ITYPE,ISENSE,SEQC2,ITOTPC,
|
|
+ IDIM2,IDOUT,LLINO,LINCON,IFAIL,IDBSIZ,IDEV,
|
|
+ IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,IDEV4,MAXGEL)
|
|
IF(IFAIL.NE.0) RETURN
|
|
IEMPTY=0
|
|
IMATC = 1
|
|
ELSE
|
|
ITYPE= - 1
|
|
ISENSE=1
|
|
JOINT = 1
|
|
LLINO = 1
|
|
LINCON = IDBSIZ - 1
|
|
ITOTPC = 0
|
|
IDOUT=MAXGEL
|
|
CALL AENTER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ SEQ2,NAMARC,JOINT,ITYPE,ISENSE,SEQC2,ITOTPC,
|
|
+ IDIM2,IDOUT,LLINO,LINCON,IFAIL,IDBSIZ,IDEV,
|
|
+ IDEV1,IDEV2,IDEV3,IDEVT,IDEVT,IDEV4,MAXGEL)
|
|
IF(IFAIL.NE.0) RETURN
|
|
END IF
|
|
ELSE IF(IOPT.EQ.3) THEN
|
|
ITYPE=0
|
|
ISENSE=1
|
|
IDOUT=MAXGEL
|
|
CALL AENTER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ SEQ2,NAMARC,JOINT,ITYPE,ISENSE,SEQC2,ITOTPC,
|
|
+ IDIM2,IDOUT,LLINO,LINCON,IFAIL,IDBSIZ,IDEV,
|
|
+ IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,IDEV4,MAXGEL)
|
|
IF(IFAIL.NE.0) RETURN
|
|
END IF
|
|
END
|
|
SUBROUTINE DBCOPY(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,GEL,
|
|
+NAMPRO,IDEV,IDBSIZ,IERR,KBIN,KBOUT,IDEVR,IDEVW,IDEVN,
|
|
+IDEVT,IDEVC,
|
|
+IHELPS,IHELPE,FILEH,IDEVH,MAXGEL,MAXDB,IDM)
|
|
CHARACTER FILEH*(*)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ)
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
CHARACTER NAMPRO*(*),NAMARC*16,GEL(MAXGEL)
|
|
CHARACTER V2,V1,VIN
|
|
INTEGER IWORD,DELDB,ACTF
|
|
PARAMETER (IWORD=4)
|
|
PARAMETER (MAXPRM = 21)
|
|
CHARACTER PERR(2)*(MAXPRM)
|
|
EXTERNAL DELDB,ACTF
|
|
IERR=1
|
|
LL = INDEX(NAMPRO,'.')
|
|
C
|
|
C save incoming version
|
|
C
|
|
VIN = NAMPRO(LL+1:LL+1)
|
|
1 CONTINUE
|
|
L = 1
|
|
V1='1'
|
|
CALL GTSTR('Make version',V1,V2,L,KBOUT,KBIN,INFLAG)
|
|
CALL CCASE(V2,1)
|
|
IF(INFLAG.EQ.2) RETURN
|
|
IF(INFLAG.EQ.1) THEN
|
|
CALL HELP2(IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
GO TO 1
|
|
END IF
|
|
IF(L.LT.1) V2 = V1
|
|
C
|
|
C check if file is open
|
|
C
|
|
NAMPRO(LL+1:LL+1) = V2
|
|
LM1 = LL - 1
|
|
IOK = ACTF(1,NAMPRO,LM1,V2,KBOUT)
|
|
IF (IOK.NE.0) THEN
|
|
NAMPRO(LL+1:LL+1) = VIN
|
|
RETURN
|
|
END IF
|
|
MN = NGELS + NCONTS + 1
|
|
MX = MAXDB
|
|
NEWSIZ = IDBSIZ
|
|
CALL GETINT(MN,MX,NEWSIZ,
|
|
+'New database size',
|
|
+IVAL,KBIN,KBOUT,IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 40
|
|
CALL BUSY(KBOUT)
|
|
NEWSIZ = IVAL
|
|
C WHERE SHOULD CHARS BE CHANGED ?
|
|
LLL = INDEX(NAMPRO,'.') + 1
|
|
NAMPRO(LLL:)='RL'//V2
|
|
CALL OPENRS(IDEV,NAMPRO,IOK,4,3)
|
|
IF(IOK.NE.0) THEN
|
|
C problem opening file
|
|
IF(IOK.EQ.2) THEN
|
|
CALL ERROM(KBOUT,'File already exists')
|
|
PERR(1) = 'Retype version number'
|
|
PERR(2) = 'Replace database'
|
|
IDO = 1
|
|
CALL RADION('Select action',PERR,2,IDO,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
IF(IDO.LT.1) GO TO 40
|
|
IF(IDO.EQ.1) THEN
|
|
C
|
|
C close busy file (note we dont care what we send for job=2!!!)
|
|
C
|
|
NAMPRO(LL+1:LL+1) = V2
|
|
IOK = ACTF(2,NAMPRO,LLL,V2,KBOUT)
|
|
IF (IOK.NE.0) THEN
|
|
C should never get here
|
|
NAMPRO(LL+1:LL+1) = VIN
|
|
RETURN
|
|
END IF
|
|
GO TO 1
|
|
END IF
|
|
IF(IDO.EQ.2) THEN
|
|
IOK = DELDB(NAMPRO,V2,IDEV,MAXGEL)
|
|
IF(IOK.EQ.0) THEN
|
|
LLL = INDEX(NAMPRO,'.') + 1
|
|
NAMPRO(LLL:)='RL'//V2
|
|
CALL OPENRS(IDEV,NAMPRO,IOK,4,3)
|
|
IF(IOK.EQ.0) GO TO 2
|
|
ELSE
|
|
CALL ERROM(KBOUT,'File delete failed')
|
|
END IF
|
|
END IF
|
|
END IF
|
|
GO TO 100
|
|
END IF
|
|
2 CONTINUE
|
|
CALL WRITER(IDEV,0,MAXDB,NEWSIZ,MAXGEL,IDM)
|
|
CALL WRITER(IDEV,NEWSIZ,NGELS,NCONTS,NGELS,NCONTS)
|
|
DO 10 I=1,NGELS
|
|
CALL WRITER(IDEV,I,RELPG(I),LNGTHG(I),LNBR(I),RNBR(I))
|
|
10 CONTINUE
|
|
M=NEWSIZ-NCONTS
|
|
N=IDBSIZ-NCONTS
|
|
DO 15 I=N,IDBSIZ-1
|
|
CALL WRITER(IDEV,M,RELPG(I),LNGTHG(I),LNBR(I),RNBR(I))
|
|
M=M+1
|
|
15 CONTINUE
|
|
CLOSE(UNIT=IDEV)
|
|
C DO SEQUENCES
|
|
NAMPRO(LLL:)='SQ'//V2
|
|
IREC=MAXGEL/IWORD
|
|
IF(MOD(MAXGEL,IWORD).NE.0)IREC=IREC+1
|
|
CALL OPENRS(IDEV,NAMPRO,IOK,IREC,3)
|
|
IF(IOK.NE.0)GO TO 100
|
|
DO 20 I=1,NGELS
|
|
CALL READW(IDEVW,I,GEL,MAXGEL)
|
|
CALL WRITEW(IDEV,I,GEL,MAXGEL)
|
|
20 CONTINUE
|
|
CLOSE(UNIT=IDEV)
|
|
C DO ARCHIVE NAMES
|
|
NAMPRO(LLL:)='AR'//V2
|
|
CALL OPENRS(IDEV,NAMPRO,IOK,4,3)
|
|
IF(IOK.NE.0)GO TO 100
|
|
DO 30 I=1,NGELS
|
|
CALL READN(IDEVN,I,NAMARC)
|
|
CALL WRITEN(IDEV,I,NAMARC)
|
|
30 CONTINUE
|
|
CLOSE(UNIT=IDEV)
|
|
IF(IDEVT.GT.0.AND.IDEVC.GT.0) THEN
|
|
NAMPRO(LLL:)='TG'//V2
|
|
CALL OPENRS(IDEV,NAMPRO,IOK,5,3)
|
|
IF(IOK.NE.0)GO TO 100
|
|
CALL COPYTG(IDEVT,IDEV,IOK,IDBSIZ,NEWSIZ,NGELS)
|
|
CLOSE(UNIT=IDEV)
|
|
NAMPRO(LLL:)='CC'//V2
|
|
C COMMENT_LENGTH: 11 = (40 + long)/long
|
|
CALL OPENRS(IDEV,NAMPRO,IOK,11,3)
|
|
IF(IOK.NE.0)GO TO 100
|
|
CALL COPYCC(IDEVC,IDEV,IOK)
|
|
CLOSE(UNIT=IDEV)
|
|
ENDIF
|
|
IERR=0
|
|
40 CONTINUE
|
|
C
|
|
C close busy file (note we dont care what we send for job=2!!!)
|
|
C
|
|
LL = INDEX(NAMPRO,'.')
|
|
NAMPRO(LL+1:LL+1) = V2
|
|
IOK = ACTF(2,NAMPRO,LL,V2,KBOUT)
|
|
NAMPRO(LL+1:LL+1) = VIN
|
|
RETURN
|
|
100 CONTINUE
|
|
WRITE(KBOUT,9999)
|
|
9999 FORMAT(' Error opening new database, copy aborted')
|
|
C
|
|
C close busy file
|
|
C
|
|
LL = INDEX(NAMPRO,'.')
|
|
NAMPRO(LL+1:LL+1) = V2
|
|
IOK = ACTF(2,NAMPRO,LL,V2,KBOUT)
|
|
NAMPRO(LL+1:LL+1) = VIN
|
|
END
|
|
INTEGER FUNCTION DELDB(NAMPRO,VERSN,IDEV,MAXGEL)
|
|
CHARACTER NAMPRO*(*),VERSN
|
|
INTEGER DELF
|
|
EXTERNAL DELF
|
|
C
|
|
C delete an xdap database (bap): names are 16 characters!
|
|
C
|
|
C assume relationships are 4 words, names are 4 and seqs are maxgel
|
|
C all recls in BYTES
|
|
C
|
|
DELDB = 1
|
|
LLL = INDEX(NAMPRO,'.') + 1
|
|
NAMPRO(LLL:)='RL'//VERSN
|
|
IF(DELF(NAMPRO,IDEV,16,4).NE.0) RETURN
|
|
NAMPRO(LLL:)='AR'//VERSN
|
|
IF(DELF(NAMPRO,IDEV,16,4).NE.0) RETURN
|
|
NAMPRO(LLL:)='SQ'//VERSN
|
|
IF(DELF(NAMPRO,IDEV,MAXGEL,4).NE.0) RETURN
|
|
NAMPRO(LLL:)='CC'//VERSN
|
|
C COMMENT_LENGTH: 11 = (40 + long)/long
|
|
IF(DELF(NAMPRO,IDEV,44,4).NE.0) RETURN
|
|
NAMPRO(LLL:)='TG'//VERSN
|
|
IF(DELF(NAMPRO,IDEV,20,4).NE.0) RETURN
|
|
DELDB = 0
|
|
END
|
|
SUBROUTINE DBFIX(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+GEL,GEL2,IDBSIZ,KBIN,KBOUT,IDEVR,IDEVW,IDEVN,
|
|
+IHELPS,IHELPE,IHELP1,IHELP2,FILEH,IDEVH,MAXGEL,MAXGLM,
|
|
+IDEVT,IDEVC,TEMP)
|
|
CHARACTER FILEH*(*)
|
|
C AUTHOR: RODGER STADEN
|
|
C 12-12-90 Added function to change raw data parameter file
|
|
C and changed menu routines accordingly
|
|
INTEGER RELPG(IDBSIZ),X,TEMP(IDBSIZ)
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
CHARACTER GEL(MAXGEL)
|
|
CHARACTER NAME*16,NEWNAM*16
|
|
INTEGER L,M,N
|
|
CHARACTER GEL2(MAXGEL)
|
|
PARAMETER (MAXPRM = 32)
|
|
CHARACTER PROMPT(9)*(MAXPRM)
|
|
INTEGER GCLIN,CHAINL
|
|
EXTERNAL GCLIN,CHAINL
|
|
WRITE(KBOUT,1000)
|
|
1000 FORMAT(
|
|
+' Warning:',
|
|
+' make a copy first, and check logical consistency after use')
|
|
10 CONTINUE
|
|
C
|
|
C SELECT OPTION
|
|
C CALL BELL(1,KBOUT)
|
|
C DBMENU now defunct for bap - so we use RADION instead
|
|
C CALL DBMENU(4,NOPT,KOPT,IHELPS,IHELPE,FILEH,IDEVH,
|
|
C +KBIN,KBOUT)
|
|
PROMPT(1) = 'Line change'
|
|
PROMPT(2) = 'Check logical consistency'
|
|
PROMPT(3) = 'Delete contig line'
|
|
PROMPT(4) = 'Shift'
|
|
PROMPT(5) = 'Move gel reading'
|
|
PROMPT(6) = 'Rename gel reading'
|
|
PROMPT(7) = 'Break a contig'
|
|
PROMPT(8) = 'Remove a gel reading'
|
|
PROMPT(9) = 'Alter raw data parameters'
|
|
NOPT = 1
|
|
CALL RADION('Alter relationships', PROMPT, 9, NOPT, IHELPS,
|
|
+ IHELPE, FILEH, IDEVH, KBIN, KBOUT)
|
|
IF(NOPT.LT.1)RETURN
|
|
IF(NOPT.EQ.1)THEN
|
|
C LINE CHANGE
|
|
MN = 0
|
|
MX = IDBSIZ
|
|
LNO = 0
|
|
CALL GETINT(MN,MX,LNO,
|
|
+ 'Number of line to change',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
IF(IVAL.EQ.0) GO TO 18
|
|
LNO = IVAL
|
|
IF(LNO.EQ.IDBSIZ)GO TO 19
|
|
WRITE(KBOUT,*)'Current line'
|
|
WRITE(KBOUT,1001)RELPG(LNO),LNGTHG(LNO),LNBR(LNO),RNBR(LNO)
|
|
1001 FORMAT(' ',4I6)
|
|
MN = 0
|
|
MX = 99999
|
|
X = RELPG(LNO)
|
|
CALL GETINT(MN,MX,X,
|
|
+ 'Relative position',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
X = IVAL
|
|
MN = -MAXGEL
|
|
MX = 999999
|
|
L = LNGTHG(LNO)
|
|
CALL GETINT(MN,MX,L,
|
|
+ 'Length',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
L = IVAL
|
|
MN = 0
|
|
MX = IDBSIZ
|
|
M = LNBR(LNO)
|
|
CALL GETINT(MN,MX,M,
|
|
+ 'Left neighbour',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
M = IVAL
|
|
MN = 0
|
|
MX = IDBSIZ
|
|
N = RNBR(LNO)
|
|
CALL GETINT(MN,MX,N,
|
|
+ 'Right neighbour',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
N = IVAL
|
|
CALL WRITER(IDEVR,LNO,X,L,M,N)
|
|
RELPG(LNO)=X
|
|
LNGTHG(LNO)=L
|
|
LNBR(LNO)=M
|
|
RNBR(LNO)=N
|
|
GO TO 10
|
|
18 CONTINUE
|
|
C
|
|
C deal with record 1/0 which contains: maxdb,actualdbsiz,maxgel,idm
|
|
C
|
|
WRITE(KBOUT,*)
|
|
+ 'Extreme caution: after this record is changed the program'
|
|
WRITE(KBOUT,*)
|
|
+ 'should be restarted, and could malfunction!'
|
|
CALL READR(IDEVR,0,IMSIZ,IASIZ,MXG,IDA)
|
|
WRITE(KBOUT,*)'Current line'
|
|
WRITE(KBOUT,1001)IMSIZ,IASIZ,MXG,IDA
|
|
MN = 0
|
|
MX = 999999
|
|
X = IMSIZ
|
|
CALL GETINT(MN,MX,X,
|
|
+ 'Maximum database size',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
X = IVAL
|
|
MN = 0
|
|
MX = X
|
|
L = IASIZ
|
|
CALL GETINT(MN,MX,L,
|
|
+ 'Actual database size',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
L = IVAL
|
|
MN = 0
|
|
MX = MAXGLM
|
|
M = MXG
|
|
CALL GETINT(MN,MX,M,
|
|
+ 'Maximum reading length',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
M = IVAL
|
|
ID = 0
|
|
IF(IDA.EQ.26) ID = 1
|
|
CALL YESONO(ID,'Database is for DNA',
|
|
+ 'Database is for protein',
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
IF(ID.LT.0) GO TO 10
|
|
N = 5
|
|
IF(ID.EQ.0)N = 5
|
|
IF(ID.EQ.1)N = 26
|
|
CALL WRITER(IDEVR,0,X,L,M,N)
|
|
GO TO 10
|
|
19 CONTINUE
|
|
C NCONTS NGELS LINES
|
|
MN = 0
|
|
MX = IDBSIZ
|
|
LL = NGELS
|
|
CALL GETINT(MN,MX,LL,
|
|
+ 'Number of gel readings',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
LL = IVAL
|
|
MN = 0
|
|
MX = IDBSIZ
|
|
MM = NCONTS
|
|
CALL GETINT(MN,MX,MM,
|
|
+ 'Number of contigs',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
MM = IVAL
|
|
CALL WRITER(IDEVR,IDBSIZ,LL,MM,LL,MM)
|
|
NGELS=LL
|
|
NCONTS=MM
|
|
GO TO 10
|
|
END IF
|
|
C
|
|
C
|
|
IF(NOPT.EQ.4)THEN
|
|
MN = 0
|
|
MX = NGELS
|
|
LNO = 0
|
|
CALL GETINT(MN,MX,LNO,
|
|
+ 'Number of first gel reading to shift',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
IF(IVAL.LT.1) GO TO 10
|
|
LNO = IVAL
|
|
I = CHAINL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,LNO)
|
|
IF(I.EQ.0)THEN
|
|
WRITE(KBOUT,*)
|
|
+ 'Problem with this gel reading. Check logical consistency'
|
|
WRITE(KBOUT,*)'of database. Shift not done'
|
|
GO TO 10
|
|
END IF
|
|
NCONTO = GCLIN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,I)
|
|
IF(NCONTO.EQ.0)THEN
|
|
WRITE(KBOUT,*)
|
|
+ 'No contig line for this contig. Check logical'
|
|
WRITE(KBOUT,*)'consistency of database. Shift not done'
|
|
GO TO 10
|
|
END IF
|
|
MN = 1 - RELPG(LNO)
|
|
MX = RELPG(NCONTO) - RELPG(LNO)
|
|
X = MN
|
|
CALL GETINT(MN,MX,X,
|
|
+ 'Distance to shift',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
X = IVAL
|
|
CALL SHIFTC(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDEVR,
|
|
+ IDBSIZ,LNO,NCONTO,X)
|
|
WRITE(KBOUT,*)'Shift complete'
|
|
GO TO 10
|
|
END IF
|
|
C
|
|
IF (NOPT.EQ.2) THEN
|
|
CALL DBCHEK(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ TEMP,IERR,KBOUT)
|
|
GO TO 10
|
|
END IF
|
|
IF(NOPT.EQ.3)THEN
|
|
WRITE(KBOUT,*)'Remove a contig'
|
|
MN = 0
|
|
MX = NGELS
|
|
LNO = 0
|
|
CALL GETINT(MN,MX,LNO,
|
|
+ 'Number of a reading in the contig to delete',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
IF(IVAL.EQ.0) GO TO 10
|
|
LNO = IVAL
|
|
I = CHAINL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,LNO)
|
|
ICONT = GCLIN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,I)
|
|
IF(ICONT.EQ.0) THEN
|
|
WRITE(KBOUT,*)'No contig line for this reading'
|
|
GO TO 10
|
|
END IF
|
|
CALL REMCON(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ ICONT,GEL,I,IDEVR,IDEVW,IDEVN,MAXGEL,KBOUT,TEMP)
|
|
END IF
|
|
IF(NOPT.EQ.6)THEN
|
|
MN = 0
|
|
MX = NGELS
|
|
LNO = 0
|
|
CALL GETINT(MN,MX,LNO,
|
|
+ 'Number of gel reading to rename',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
IF(IVAL.LT.1) GO TO 10
|
|
LNO = IVAL
|
|
CALL READN(IDEVN,LNO,NAME)
|
|
L = 16
|
|
CALL GTSTR('name for gel reading',
|
|
+ NAME,NEWNAM,L,KBOUT,KBIN,INFLAG)
|
|
IF(L.GT.0)CALL WRITEN(IDEVN,LNO,NEWNAM)
|
|
GO TO 10
|
|
END IF
|
|
IF(NOPT.EQ.5)THEN
|
|
MN = 0
|
|
MX = NGELS
|
|
IFROM = 0
|
|
CALL GETINT(MN,MX,IFROM,
|
|
+ 'Number of gel reading to move',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
IF(IVAL.LT.1) GO TO 10
|
|
IFROM = IVAL
|
|
MN = 0
|
|
MX = NGELS
|
|
ITO = 0
|
|
CALL GETINT(MN,MX,ITO,
|
|
+ 'New number for gel reading',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
IF(IVAL.LT.1) GO TO 10
|
|
ITO = IVAL
|
|
CALL MOVGEL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ GEL,IFROM,ITO,IDEVR,IDEVW,IDEVN,MAXGEL,KBOUT)
|
|
GO TO 10
|
|
END IF
|
|
IF(NOPT.EQ.7)THEN
|
|
CALL BREAKC(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,KBIN,KBOUT,IDEVR,IDEVW,IDEVN,
|
|
+ IHELPS,IHELPE,IHELP1,IHELP2,FILEH,IDEVH,IOK)
|
|
GO TO 10
|
|
END IF
|
|
IF(NOPT.EQ.8) THEN
|
|
CALL REMGD(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ KBIN,KBOUT,GEL,MAXGEL,IDEVR,IDEVW,IDEVN,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
GO TO 10
|
|
END IF
|
|
IF(NOPT.EQ.9) THEN
|
|
CALL FIXRD(IDEVT,IDEVC,IDBSIZ,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH)
|
|
GO TO 10
|
|
END IF
|
|
GO TO 10
|
|
END
|
|
SUBROUTINE REMCON(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ICONT,GEL,LLINO,IDEVR,IDEVW,IDEVN,MAXGEL,KBOUT,TEMP)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),TEMP(IDBSIZ)
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
CHARACTER GEL(MAXGEL)
|
|
C
|
|
C problem is to remove a contig. Strategy is to find all its read numbers
|
|
C and store them in temp. Then go thru and replace them by reads from the
|
|
C end of the list of reads. An earlier version tried to chain thru but
|
|
C ran into complications about replacing reads by reads they pointed to!
|
|
C
|
|
DO 10 I=1,IDBSIZ
|
|
TEMP(I) = 0
|
|
10 CONTINUE
|
|
C
|
|
NDEL = 0
|
|
I = LLINO
|
|
20 CONTINUE
|
|
TEMP(I) = 1
|
|
NDEL = NDEL + 1
|
|
IF (RNBR(I).NE.0) THEN
|
|
I = RNBR(I)
|
|
GO TO 20
|
|
END IF
|
|
C
|
|
C now in temp all the reads in the contig are 1, the rest of temp is zero
|
|
C let i be the read to move and j the last read left in the list of reads
|
|
C so if temp(i) is 1 and temp(j) is zero move read j to i and set j = j - 1
|
|
C then deal with the next i.
|
|
C if temp(i) is 1 and temp(j) is also 1 simply set j = j - 1 and try to move
|
|
C that one
|
|
C we stop when weve gone so far along the list that the next read to
|
|
C delete is equal to the number of reads left in the list
|
|
C
|
|
C what are the difficult cases?
|
|
C 1. only one contig
|
|
C 2. all reads to right of llino are are near the end of the list of reads
|
|
C
|
|
I = 1
|
|
J = NGELS
|
|
30 CONTINUE
|
|
IF (TEMP(I).EQ.1) THEN
|
|
IF (TEMP(J).EQ.0) THEN
|
|
WRITE(*,*)'MOVE ',J,' TO ',I
|
|
CALL MOVGEL(RELPG,LNGTHG,LNBR,RNBR,J,NCONTS,IDBSIZ,
|
|
+ GEL,J,I,IDEVR,IDEVW,IDEVN,MAXGEL,KBOUT)
|
|
J = J - 1
|
|
ELSE
|
|
J = J - 1
|
|
IF (I.LT.J) GO TO 30
|
|
END IF
|
|
END IF
|
|
IF (I.LT.J) THEN
|
|
I = I + 1
|
|
GO TO 30
|
|
END IF
|
|
C
|
|
C fix up number of reads and the contig record
|
|
C
|
|
NGELS = NGELS - NDEL
|
|
CALL REMCNL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ICONT,IDEVR)
|
|
END
|
|
SUBROUTINE REMGD(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+KBIN,KBOUT,GEL,MAXGEL,IDEVR,IDEVW,IDEVN,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
INTEGER RELPG(IDBSIZ),LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
CHARACTER HELPF*(*),GEL(MAXGEL)
|
|
INTEGER REMME,GCLIN,CHAINL
|
|
EXTERNAL GCLIN,CHAINL
|
|
C assumes db is logical consistent
|
|
WRITE(KBOUT,*)'Remove reading from database'
|
|
REMME = NGELS
|
|
CALL GETINT(1,NGELS,REMME,
|
|
+'Number of reading to remove',
|
|
+IVAL,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
REMME = IVAL
|
|
I = CHAINL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,REMME)
|
|
ICONT = GCLIN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,I)
|
|
IF(ICONT.EQ.0) THEN
|
|
WRITE(KBOUT,*)'No contig line for this reading'
|
|
IOK = 1
|
|
RETURN
|
|
END IF
|
|
CALL REMGEL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+REMME,ICONT,KBOUT,GEL,MAXGEL,IDEVR,IDEVW,IDEVN)
|
|
END
|
|
SUBROUTINE REMGEL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+REMME,ICONT,KBOUT,GEL,MAXGEL,IDEVR,IDEVW,IDEVN)
|
|
C Routine to remove a reading from a database
|
|
C Cases: 1 left end
|
|
C 2 right end
|
|
C 3 internal and dispensible
|
|
C 4 internal and indispensible
|
|
C if 1 change contig lnbr, contig length, lnbr of rnbr of remme, relpgs
|
|
C if 2 change contig rnbr, contig length, rnbr of lnbr of remme
|
|
C if 3 change contig length, lnbr of rnbr of remme rnbr of lnbr of remme
|
|
C if 4 need to break contig, then as for 1
|
|
C if 1 and 2 then also remove contig line
|
|
C for all cases move gel ngels to remme (if remme/=ngels)
|
|
C and update line idbsiz
|
|
C
|
|
INTEGER RELPG(IDBSIZ),LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
INTEGER REMME,CLEN
|
|
LOGICAL LEFTE,RIGHTE,DISPEN
|
|
CHARACTER GEL(MAXGEL)
|
|
EXTERNAL CLEN
|
|
LEFTE = .FALSE.
|
|
RIGHTE = .FALSE.
|
|
DISPEN = .FALSE.
|
|
C
|
|
C Left end ?
|
|
C
|
|
IF(LNBR(REMME).EQ.0) LEFTE = .TRUE.
|
|
C
|
|
C Right end ?
|
|
C
|
|
IF(RNBR(REMME).EQ.0) RIGHTE = .TRUE.
|
|
C
|
|
C If both true remove the contig line, then overwrite the gel
|
|
C
|
|
IF(LEFTE.AND.RIGHTE) THEN
|
|
WRITE(KBOUT,*)'Removing reading and contig'
|
|
IFROM = NGELS
|
|
NGELS = NGELS - 1
|
|
CALL REMCNL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ ICONT,IDEVR)
|
|
IF(REMME.NE.IFROM) THEN
|
|
WRITE(KBOUT,*)'Renumbering reading',IFROM,' to',REMME
|
|
CALL MOVGEL(RELPG,LNGTHG,LNBR,RNBR,NGELS,
|
|
+ NCONTS,IDBSIZ,GEL,IFROM,REMME,IDEVR,IDEVW,IDEVN,MAXGEL,KBOUT)
|
|
END IF
|
|
ELSE IF(LEFTE) THEN
|
|
WRITE(KBOUT,*)'Removing reading from left end of contig'
|
|
LNBR(ICONT) = RNBR(REMME)
|
|
I = 1 - RELPG(RNBR(REMME))
|
|
WRITE(KBOUT,*)'Shifting readings in contig by distance=',I
|
|
CALL SHIFTC(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDEVR,
|
|
+ IDBSIZ,RNBR(REMME),ICONT,I)
|
|
I = LNBR(ICONT)
|
|
LNBR(I) = 0
|
|
CALL WRITER(IDEVR,I,RELPG(I),LNGTHG(I),
|
|
+ LNBR(I),RNBR(I))
|
|
IFROM = NGELS
|
|
IF(REMME.NE.IFROM) THEN
|
|
WRITE(KBOUT,*)'Renumbering reading',IFROM,' to',REMME
|
|
CALL MOVGEL(RELPG,LNGTHG,LNBR,RNBR,NGELS,
|
|
+ NCONTS,IDBSIZ,GEL,IFROM,REMME,IDEVR,IDEVW,IDEVN,MAXGEL,KBOUT)
|
|
END IF
|
|
NGELS = NGELS - 1
|
|
CALL WRITER(IDEVR,IDBSIZ,NGELS,NCONTS,NGELS,NCONTS)
|
|
ELSE IF(RIGHTE) THEN
|
|
WRITE(KBOUT,*)'Removing reading from right end of contig'
|
|
RNBR(ICONT) = LNBR(REMME)
|
|
I = RNBR(ICONT)
|
|
RNBR(I) = 0
|
|
CALL WRITER(IDEVR,I,RELPG(I),LNGTHG(I),
|
|
+ LNBR(I),RNBR(I))
|
|
RELPG(ICONT) = CLEN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,LNBR(ICONT))
|
|
CALL WRITER(IDEVR,ICONT,RELPG(ICONT),LNGTHG(ICONT),
|
|
+ LNBR(ICONT),RNBR(ICONT))
|
|
IFROM = NGELS
|
|
IF(REMME.NE.IFROM) THEN
|
|
WRITE(KBOUT,*)'Renumbering reading',IFROM,' to',REMME
|
|
CALL MOVGEL(RELPG,LNGTHG,LNBR,RNBR,NGELS,
|
|
+ NCONTS,IDBSIZ,GEL,IFROM,REMME,IDEVR,IDEVW,IDEVN,MAXGEL,KBOUT)
|
|
END IF
|
|
NGELS = NGELS - 1
|
|
CALL WRITER(IDEVR,IDBSIZ,NGELS,NCONTS,NGELS,NCONTS)
|
|
ELSE
|
|
C Is remme indispensible ?
|
|
NSTART = RELPG(RNBR(REMME))
|
|
I = REMME
|
|
10 CONTINUE
|
|
I = LNBR(I)
|
|
IF(I.NE.0) THEN
|
|
IF((RELPG(I)+ABS(LNGTHG(I))-1).LT.NSTART) GO TO 10
|
|
DISPEN = .TRUE.
|
|
END IF
|
|
IF(DISPEN) THEN
|
|
WRITE(KBOUT,*)
|
|
+ 'Removing dispensible reading from middle of contig'
|
|
I = LNBR(REMME)
|
|
RNBR(I) = RNBR(REMME)
|
|
CALL WRITER(IDEVR,I,RELPG(I),LNGTHG(I),
|
|
+ LNBR(I),RNBR(I))
|
|
I = RNBR(REMME)
|
|
LNBR(I) = LNBR(REMME)
|
|
CALL WRITER(IDEVR,I,RELPG(I),LNGTHG(I),
|
|
+ LNBR(I),RNBR(I))
|
|
IFROM = NGELS
|
|
IF(REMME.NE.IFROM) THEN
|
|
WRITE(KBOUT,*)'Renumbering reading',IFROM,' to',REMME
|
|
CALL MOVGEL(RELPG,LNGTHG,LNBR,RNBR,NGELS,
|
|
+ NCONTS,IDBSIZ,GEL,IFROM,REMME,IDEVR,IDEVW,IDEVN,
|
|
+ MAXGEL,KBOUT)
|
|
END IF
|
|
NGELS = NGELS - 1
|
|
RELPG(ICONT) = CLEN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,LNBR(ICONT))
|
|
CALL WRITER(IDEVR,ICONT,RELPG(ICONT),LNGTHG(ICONT),
|
|
+ LNBR(ICONT),RNBR(ICONT))
|
|
CALL WRITER(IDEVR,IDBSIZ,NGELS,NCONTS,NGELS,NCONTS)
|
|
ELSE
|
|
WRITE(KBOUT,*)
|
|
+ 'Removing indispensible reading from middle of contig'
|
|
WRITE(KBOUT,*)'So breaking contig first'
|
|
IR = REMME
|
|
IL = LNBR(REMME)
|
|
ILO = LNBR(ICONT)
|
|
NCONTO = ICONT
|
|
NCONTR = IDBSIZ - NCONTS - 1
|
|
CALL CBREAK(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ KBOUT,IDEVR,IDEVW,IDEVN,IR,IL,ILO,NCONTO,NCONTR,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
WRITE(KBOUT,*)'Removing reading from left end of contig'
|
|
ICONT = IDBSIZ - NCONTS
|
|
LNBR(ICONT) = RNBR(REMME)
|
|
I = 1 - RELPG(RNBR(REMME))
|
|
WRITE(KBOUT,*)'Shifting readings in contig by distance=',I
|
|
CALL SHIFTC(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDEVR,
|
|
+ IDBSIZ,RNBR(REMME),ICONT,I)
|
|
I = LNBR(ICONT)
|
|
LNBR(I) = 0
|
|
CALL WRITER(IDEVR,I,RELPG(I),LNGTHG(I),
|
|
+ LNBR(I),RNBR(I))
|
|
IFROM = NGELS
|
|
IF(REMME.NE.IFROM) THEN
|
|
WRITE(KBOUT,*)'Renumbering reading',IFROM,' to',REMME
|
|
CALL MOVGEL(RELPG,LNGTHG,LNBR,RNBR,NGELS,
|
|
+ NCONTS,IDBSIZ,GEL,IFROM,REMME,IDEVR,IDEVW,IDEVN,
|
|
+ MAXGEL,KBOUT)
|
|
END IF
|
|
NGELS = NGELS - 1
|
|
CALL WRITER(IDEVR,IDBSIZ,NGELS,NCONTS,NGELS,NCONTS)
|
|
END IF
|
|
END IF
|
|
END
|
|
SUBROUTINE REMCNL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+REMME,IDEVR)
|
|
INTEGER RELPG(IDBSIZ),LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
INTEGER REMME
|
|
C Routine to remove a contig line from a db
|
|
C Loop deals with case of remove top contig
|
|
C Move down all lines from above
|
|
DO 10 I = REMME,IDBSIZ-NCONTS+1,-1
|
|
RELPG(I) = RELPG(I-1)
|
|
LNGTHG(I) = LNGTHG(I-1)
|
|
LNBR(I) = LNBR(I-1)
|
|
RNBR(I) = RNBR(I-1)
|
|
CALL WRITER(IDEVR,I,RELPG(I),LNGTHG(I),LNBR(I),RNBR(I))
|
|
10 CONTINUE
|
|
NCONTS = NCONTS - 1
|
|
CALL WRITER(IDEVR,IDBSIZ,NGELS,NCONTS,NGELS,NCONTS)
|
|
END
|
|
SUBROUTINE MOVGEL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+GEL,FROM,TO,IDEVR,IDEVW,IDEVN,MAXGEL,KBOUT)
|
|
C Subroutine to move a gel from line from to line to
|
|
C Extended 22-5-91
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),FROM,TO
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
CHARACTER NAMGEL*16,GEL(MAXGEL)
|
|
INTEGER GCLIN,CHAINL
|
|
LOGICAL LEFTE,RIGHTE
|
|
EXTERNAL GCLIN,CHAINL
|
|
LEFTE = .FALSE.
|
|
RIGHTE = .FALSE.
|
|
C
|
|
C left end ?
|
|
C
|
|
IF(LNBR(FROM).EQ.0) LEFTE = .TRUE.
|
|
C
|
|
C right end ?
|
|
C
|
|
IF(RNBR(FROM).EQ.0) RIGHTE = .TRUE.
|
|
C
|
|
C if both true remove the contig line, then overwrite the gel
|
|
C
|
|
IF(LEFTE.AND.RIGHTE) THEN
|
|
NCONTO = GCLIN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,FROM)
|
|
IF(NCONTO.EQ.0)THEN
|
|
WRITE(KBOUT,*)
|
|
+ 'This gel has no left neighbour but does not'
|
|
WRITE(KBOUT,*)'appear in a contig line!'
|
|
ELSE
|
|
LNBR(NCONTO) = TO
|
|
RNBR(NCONTO) = TO
|
|
CALL WRITER(IDEVR,NCONTO,RELPG(NCONTO),LNGTHG(NCONTO),
|
|
+ LNBR(NCONTO),RNBR(NCONTO))
|
|
END IF
|
|
ELSE IF(LEFTE) THEN
|
|
NCONTO = GCLIN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,FROM)
|
|
IF(NCONTO.EQ.0)THEN
|
|
WRITE(KBOUT,*)
|
|
+ 'This gel has no left neighbour but does not'
|
|
WRITE(KBOUT,*)'appear in a contig line!'
|
|
ELSE
|
|
LNBR(NCONTO) = TO
|
|
CALL WRITER(IDEVR,NCONTO,RELPG(NCONTO),LNGTHG(NCONTO),
|
|
+ LNBR(NCONTO),RNBR(NCONTO))
|
|
END IF
|
|
ELSE IF(RIGHTE) THEN
|
|
I = CHAINL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,FROM)
|
|
NCONTO = GCLIN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,I)
|
|
IF(NCONTO.EQ.0)THEN
|
|
WRITE(KBOUT,*)
|
|
+ 'This gel has no right neighbour and does not'
|
|
WRITE(KBOUT,*)'appear in a contig!'
|
|
ELSE
|
|
IF(RNBR(NCONTO).NE.FROM)THEN
|
|
WRITE(KBOUT,*)
|
|
+ 'This gel has no right neighbour but does not'
|
|
WRITE(KBOUT,*)'appear in a contig line!'
|
|
ELSE
|
|
RNBR(NCONTO) = TO
|
|
CALL WRITER(IDEVR,NCONTO,RELPG(NCONTO),LNGTHG(NCONTO),
|
|
+ LNBR(NCONTO),RNBR(NCONTO))
|
|
END IF
|
|
END IF
|
|
END IF
|
|
RELPG(TO)=RELPG(FROM)
|
|
LNGTHG(TO)=LNGTHG(FROM)
|
|
LNBR(TO)=LNBR(FROM)
|
|
RNBR(TO)=RNBR(FROM)
|
|
CALL READW(IDEVW,FROM,GEL,MAXGEL)
|
|
CALL WRITEW(IDEVW,TO,GEL,MAXGEL)
|
|
CALL READN(IDEVN,FROM,NAMGEL)
|
|
CALL WRITEN(IDEVN,TO,NAMGEL)
|
|
CALL WRITER(IDEVR,TO,RELPG(TO),LNGTHG(TO),
|
|
+LNBR(TO),RNBR(TO))
|
|
C Do neighbours
|
|
IF(LNBR(FROM).NE.0) THEN
|
|
I=LNBR(FROM)
|
|
RNBR(I)=TO
|
|
CALL WRITER(IDEVR,I,RELPG(I),LNGTHG(I),
|
|
+ LNBR(I),RNBR(I))
|
|
END IF
|
|
IF(RNBR(FROM).NE.0) THEN
|
|
I=RNBR(FROM)
|
|
LNBR(I)=TO
|
|
CALL WRITER(IDEVR,I,RELPG(I),LNGTHG(I),
|
|
+ LNBR(I),RNBR(I))
|
|
END IF
|
|
CALL MOVTAG(FROM,TO)
|
|
END
|
|
SUBROUTINE DBOPEN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,NAMPRO,GEL,
|
|
+IDBSIS,IDBSIZ,IERR,KBIN,KBOUT,
|
|
+IDEVR,IDEVW,IDEVN,IDEVT,IDEVC,
|
|
+MAXGEL,MAXGLM,LLINO,
|
|
+IDM,IHELPS,IHELPE,HELPF,IDEVH)
|
|
C AUTHOR: RODGER STADEN
|
|
CHARACTER GEL(MAXGLM)
|
|
INTEGER RELPG(IDBSIS)
|
|
INTEGER LNGTHG(IDBSIS),LNBR(IDBSIS),RNBR(IDBSIS)
|
|
CHARACTER NAMPRO*(*),COPYNO*4,HELPF*(*)
|
|
INTEGER IWORD,ACTF
|
|
EXTERNAL ACTF
|
|
PARAMETER (IWORD=4)
|
|
C NOTE THIS IS THE MACHINES WORD LENGTH IE HOW MANY CHARS PER WORD
|
|
CALL FILLI(RELPG,IDBSIS,0)
|
|
CALL FILLI(LNGTHG,IDBSIS,0)
|
|
CALL FILLI(LNBR,IDBSIS,0)
|
|
CALL FILLI(RNBR,IDBSIS,0)
|
|
NAMPRO(1:)=' '
|
|
IERR=1
|
|
1 CONTINUE
|
|
L = 0
|
|
CALL GTSTR('Project name',' ',NAMPRO,L,KBOUT,KBIN,INFLAG)
|
|
IF(L.LT.1) RETURN
|
|
LL = L
|
|
CALL CCASE(NAMPRO,1)
|
|
IF(INFLAG.EQ.2) RETURN
|
|
IF(INFLAG.EQ.1) THEN
|
|
CALL HELP2(IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
GO TO 1
|
|
END IF
|
|
L = 1
|
|
CALL GTSTR('Version','0',COPYNO,L,KBOUT,KBIN,INFLAG)
|
|
CALL CCASE(COPYNO,1)
|
|
IF(INFLAG.EQ.2) RETURN
|
|
IF(INFLAG.EQ.1) THEN
|
|
CALL HELP2(IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
GO TO 1
|
|
END IF
|
|
IF(L.LT.1) COPYNO(1:1) = '0'
|
|
C
|
|
C check if file is open
|
|
C
|
|
NAMPRO(LL+1:LL+2) = '.'//COPYNO(1:1)
|
|
IOK = ACTF(1,NAMPRO,LL,COPYNO,KBOUT)
|
|
IF (IOK.NE.0) RETURN
|
|
NAMPRO(LL+1:LL+3)='.AR'
|
|
NAMPRO(LL+4:LL+4)=COPYNO(1:1)
|
|
CALL OPENRS(IDEVN,NAMPRO,IOK,4,4)
|
|
IF(IOK.NE.0)GO TO 100
|
|
NAMPRO(LL+2:LL+3)='RL'
|
|
CALL OPENRS(IDEVR,NAMPRO,IOK,4,4)
|
|
IF(IOK.NE.0)GO TO 100
|
|
CALL READR(IDEVR,0,IDBST,IDBSIZ,MAXGEL,IDM)
|
|
C
|
|
C Do I really need to look at this value? I could simply compare
|
|
C IDBSIS with IDBSIZ which is what really counts. Cannot remember
|
|
C why IDBST is stored!!!!!!!
|
|
C
|
|
IF(IDBST.GT.IDBSIS) THEN
|
|
CALL ERROM(KBOUT,
|
|
+ 'Fatal error: database size too large for program')
|
|
GO TO 100
|
|
END IF
|
|
NAMPRO(LL+2:LL+3)='SQ'
|
|
C DEFINE RECORD LENGTH IN TERMS OF NUMBER OF CHARS PER WORD (4 ON VAX)
|
|
C AND MAXGEL SIZE
|
|
IREC=MAXGEL/IWORD
|
|
IF(MOD(MAXGEL,IWORD).NE.0)IREC=IREC+1
|
|
CALL OPENRS(IDEVW,NAMPRO,IOK,IREC,4)
|
|
IF(IOK.NE.0)GO TO 100
|
|
NAMPRO(LL+2:LL+3) = 'TG'
|
|
CALL OPENRS(IDEVT,NAMPRO,IOK,5,4)
|
|
IF(IOK.NE.0) IDEVT = -1
|
|
NAMPRO(LL+2:LL+3) = 'CC'
|
|
C COMMENT_LENGTH: 11 = (40 + long)/long
|
|
CALL OPENRS(IDEVC,NAMPRO,IOK,11,4)
|
|
IF(IOK.NE.0) IDEVC = -1
|
|
C READ A LINE FOR LUCK
|
|
CALL READW(IDEVW,1,GEL,MAXGEL)
|
|
CALL READR(IDEVR,IDBSIZ,NGELS,NCONTS,IDUM1,IDUM2)
|
|
WRITE(KBOUT,10011)NGELS,NCONTS,IDBSIZ,MAXGEL
|
|
10011 FORMAT(' Number of gel readings ',I6,' Number of contigs ',I6,/,
|
|
+' Database size ',I6,' Maximum gel reading length=',I4)
|
|
LLINO = 0
|
|
IF(NGELS.LT.1)GO TO 5
|
|
DO 3 I=1,NGELS
|
|
CALL READR(IDEVR,I,RELPG(I),LNGTHG(I),LNBR(I),RNBR(I))
|
|
3 CONTINUE
|
|
N=IDBSIZ-NCONTS
|
|
MXT = 0
|
|
DO 4 I=N,IDBSIZ-1
|
|
CALL READR(IDEVR,I,RELPG(I),LNGTHG(I),LNBR(I),RNBR(I))
|
|
IF(RELPG(I).GT.MXT) THEN
|
|
MXT = RELPG(I)
|
|
LLINO = LNBR(I)
|
|
END IF
|
|
4 CONTINUE
|
|
5 CONTINUE
|
|
NAMPRO(LL+2:LL+2) = COPYNO(1:1)
|
|
IERR=0
|
|
RETURN
|
|
100 CONTINUE
|
|
IOK = ACTF(2,NAMPRO,LL,COPYNO,KBOUT)
|
|
WRITE(KBOUT,9999)
|
|
9999 FORMAT(' Error encountered opening database files')
|
|
NAMPRO(LL+2:LL+2) = COPYNO(1:1)
|
|
END
|
|
SUBROUTINE DBSTAR(NAMPRO,GEL,IDBSIS,IDBSIZ,KBIN,KBOUT,
|
|
+IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,
|
|
+IERR,IHELPS,IHELPE,IDEVH,FILEH,
|
|
+MAXGEL,MAXGLM,IDM)
|
|
CHARACTER GEL(MAXGLM),FILEH*(*)
|
|
CHARACTER NAMPRO*(*)
|
|
INTEGER IWORD,ACTF
|
|
PARAMETER (IWORD=4)
|
|
EXTERNAL ACTF
|
|
IERR=1
|
|
3 CONTINUE
|
|
NAMPRO = ' '
|
|
MN = 0
|
|
CALL GTSTR('New project name',' ',NAMPRO,MN,KBOUT,KBIN,INFLAG)
|
|
IF(MN.LT.1) RETURN
|
|
LL = MIN(12,MN)
|
|
CALL CCASE(NAMPRO,1)
|
|
IF(INFLAG.EQ.2) RETURN
|
|
IF(INFLAG.EQ.1) THEN
|
|
CALL HELP2(IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
GO TO 3
|
|
END IF
|
|
NAMPRO(LL+1:LL+2) = '.0'
|
|
IOK = ACTF(1,NAMPRO,LL,'0',KBOUT)
|
|
IF (IOK.NE.0) RETURN
|
|
MN = 10
|
|
MX = IDBSIS
|
|
IDBSIZ = 50
|
|
CALL GETINT(MN,MX,IDBSIZ,
|
|
+'Database size',
|
|
+IVAL,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 40
|
|
IDBSIZ = IVAL
|
|
5 CONTINUE
|
|
MN = 512
|
|
MX = MAXGLM
|
|
MAXIN1 = MIN(512,MAXGEL)
|
|
MAXIN = 1024
|
|
CALL GETINT(MN,MX,MAXIN,
|
|
+'Maximum gel reading length',
|
|
+IVAL,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 40
|
|
MAXGEL = IVAL
|
|
IF(MOD(IVAL,MAXIN1).NE.0) THEN
|
|
MAXGEL = 512 + (IVAL/512)*512
|
|
WRITE(KBOUT,*)'Maximum set to',MAXGEL
|
|
END IF
|
|
CALL YESNO(IDM,'Database is for DNA',
|
|
+IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
IF(IDM.LT.0) GO TO 40
|
|
IF(IDM.EQ.0)IDM = 5
|
|
IF(IDM.EQ.1)IDM = 26
|
|
NAMPRO(LL+1:LL+4)='.RL0'
|
|
CALL OPENRS(IDEV1,NAMPRO,IOK,4,3)
|
|
IF(IOK.NE.0)GO TO 100
|
|
WRITE(KBOUT,1004)
|
|
1004 FORMAT(' Writing new database files')
|
|
J = 0
|
|
K = 0
|
|
N = 0
|
|
M = 0
|
|
C
|
|
C Put maximum possible database size, actual size, max read length and type
|
|
C in record 0 of relationships (actually record 1, as writer adds 1).
|
|
C
|
|
CALL WRITER(IDEV1,0,IDBSIS,IDBSIZ,MAXGEL,IDM)
|
|
C
|
|
C Put zeroes in al other records
|
|
C
|
|
DO 10 I=1,IDBSIZ
|
|
CALL WRITER(IDEV1,I,J,K,M,N)
|
|
10 CONTINUE
|
|
NAMPRO(LL+2:LL+3)='SQ'
|
|
IREC=MAXGEL/IWORD
|
|
IF(MOD(MAXGEL,IWORD).NE.0)IREC=IREC+1
|
|
CALL OPENRS(IDEV2,NAMPRO,IOK,IREC,3)
|
|
IF(IOK.NE.0)GO TO 100
|
|
CALL FILLC(GEL,MAXGEL,' ')
|
|
C CALL WRITEW(IDEV2,IDBSIZ,GEL,MAXGEL)
|
|
C
|
|
C write only the first record into the working versions
|
|
C and assume others will be added when required
|
|
C
|
|
CALL WRITEW(IDEV2,1,GEL,MAXGEL)
|
|
NAMPRO(LL+2:LL+3)='AR'
|
|
CALL OPENRS(IDEV3,NAMPRO,IOK,4,3)
|
|
IF(IOK.NE.0)GO TO 100
|
|
C
|
|
C write only the first record into the names
|
|
C and assume others will be added when required
|
|
C
|
|
CALL WRITEN(IDEV3,1,' ')
|
|
C CREATE TAG FILES (TAGS AND COMMENTS)
|
|
NAMPRO(LL+2:LL+3)='TG'
|
|
CALL OPENRS(IDEVT,NAMPRO,IOK,5,3)
|
|
IF(IOK.NE.0)GO TO 100
|
|
CALL WRITTG(IDEVT,IDBSIZ,IDBSIZ,0,0,0,0)
|
|
NAMPRO(LL+2:LL+3)='CC'
|
|
C COMMENT_LENGTH: 11 = (40 + long)/long
|
|
CALL OPENRS(IDEVC,NAMPRO,IOK,11,3)
|
|
IF(IOK.NE.0)GO TO 100
|
|
CALL WRITCC(IDEVC,1,1,0,' ')
|
|
WRITE(KBOUT,1003)NAMPRO(1:LL),IDBSIZ
|
|
1003 FORMAT(' Database ',A,' version 0, size ',I6,
|
|
+' successfully started')
|
|
IERR=0
|
|
NAMPRO(LL+2:LL+2) = '0'
|
|
RETURN
|
|
40 CONTINUE
|
|
NAMPRO(LL+1:LL+2) = '.0'
|
|
IOK = ACTF(2,NAMPRO,LL,'0',KBOUT)
|
|
RETURN
|
|
100 CONTINUE
|
|
CALL ERROM(KBOUT,'Error writing database files')
|
|
NAMPRO(LL+1:LL+2) = '.0'
|
|
IOK = ACTF(2,NAMPRO,LL,'0',KBOUT)
|
|
END
|
|
SUBROUTINE FIXRD(IDEVT,IDEVC,IDBSIZ,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,FILEH,IDEVH)
|
|
C FILE_NAME_LENGTH
|
|
CHARACTER NAMFIL*18,NEWNAM*18,MTYPE*4,NEWMT*4,FILEH*(*)
|
|
IF(IDEVRD.LT.0) THEN
|
|
WRITE(KBOUT,*)'No raw data file!'
|
|
RETURN
|
|
END IF
|
|
10 CONTINUE
|
|
C Change raw data record
|
|
MN = 0
|
|
MX = IDBSIZ-1
|
|
LNO = 0
|
|
CALL GETINT(MN,MX,LNO,
|
|
+ 'Number of line to change',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
IF(IVAL.EQ.0) RETURN
|
|
LNO = IVAL
|
|
CALL READRD(IDEVT,IDEVC,LNO,LENR,LCUT,LENW,MTYPE,NAMFIL)
|
|
WRITE(KBOUT,*)'Current line'
|
|
WRITE(KBOUT,1001)LENR,LCUT,LENW,MTYPE,NAMFIL
|
|
1001 FORMAT(' ',3I6,' ',A,' ',A)
|
|
MN = 1
|
|
MX = 9999
|
|
LX = LENR
|
|
CALL GETINT(MN,MX,LX,
|
|
+ 'Length raw sequence',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
LX = IVAL
|
|
MN = 1
|
|
MX = LX
|
|
L = LCUT
|
|
CALL GETINT(MN,MX,L,
|
|
+ 'Left cutoff',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 10
|
|
L = IVAL
|
|
MN = 1
|
|
MX = LX
|
|
M = LENW
|
|
CALL GETINT(MN,MX,M,
|
|
+ 'Length of original working sequence',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
M = IVAL
|
|
20 CONTINUE
|
|
LNAM = 4
|
|
CALL GTSTR('Machine type',
|
|
+ MTYPE,NEWMT,LNAM,KBOUT,KBIN,INFLAG)
|
|
IF(INFLAG.EQ.2) RETURN
|
|
IF(INFLAG.EQ.1) THEN
|
|
CALL HELP2(IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
GO TO 20
|
|
END IF
|
|
IF(INFLAG.EQ.3) NEWMT = MTYPE
|
|
30 CONTINUE
|
|
C FILE_NAME_LENGTH
|
|
LNAM = 18
|
|
CALL GTSTR('Name for raw data file',
|
|
+ NAMFIL,NEWNAM,LNAM,KBOUT,KBIN,INFLAG)
|
|
IF(INFLAG.EQ.2) RETURN
|
|
IF(INFLAG.EQ.1) THEN
|
|
CALL HELP2(IHELPS,IHELPE,FILEH,IDEVH,KBIN,KBOUT)
|
|
GO TO 30
|
|
END IF
|
|
IF(INFLAG.EQ.3) NEWNAM = NAMFIL
|
|
WRITE(KBOUT,1001)LX,L,M,NEWMT,NEWNAM
|
|
CALL WRITRD(IDEVT,IDEVC,LNO,LX,L,M,NEWMT,NEWNAM)
|
|
WRITE(KBOUT,*)'New record written to disk'
|
|
GO TO 10
|
|
END
|
|
SUBROUTINE READRD(IDEVT,IDEVC,NGEL,LENR,LCUT,LENW,MTYPE,NAMFIL)
|
|
CHARACTER MTYPE*(*),NAMFIL*(*)
|
|
C COMMENT_LENGTH
|
|
CHARACTER NOTE*40
|
|
IF(IDEVT.GT.0)THEN
|
|
CALL READTG(IDEVT,NGEL,LPOS,LLEN,LCOM,LTYPE,NEXT)
|
|
CALL READCC(IDEVC,LCOM,ICNT,NEXT,NOTE)
|
|
READ(NOTE,1001,ERR=100)LENR,LCUT,LENW,MTYPE,NAMFIL
|
|
1001 FORMAT(3I6,A,A)
|
|
ENDIF
|
|
RETURN
|
|
100 CONTINUE
|
|
LENR = 0
|
|
LCUT = 0
|
|
LENW = 0
|
|
MTYPE = ' '
|
|
NAMFIL = ' '
|
|
END
|
|
SUBROUTINE WRITRD(IDEVT,IDEVC,NGEL,LENR,LCUT,LENW,MTYPE,NAMFIL)
|
|
CHARACTER MTYPE*(*),NAMFIL*(*)
|
|
INTEGER FREECC
|
|
C COMMENT_LENGTH
|
|
CHARACTER NOTE*40
|
|
IF(IDEVT.GT.0)THEN
|
|
CALL READTG(IDEVT,NGEL,LPOS,LLEN,LCOM,LTYPE,NEXT)
|
|
IF(LCOM.EQ.0)THEN
|
|
LCOM = FREECC(IDEVC)
|
|
ENDIF
|
|
WRITE(NOTE,1001,ERR=100)LENR,LCUT,LENW,MTYPE,NAMFIL
|
|
1001 FORMAT(3I6,A,A)
|
|
NEXT = 0
|
|
CALL WRITCC(IDEVC,LCOM,ICNT,NEXT,NOTE)
|
|
ENDIF
|
|
RETURN
|
|
100 CONTINUE
|
|
END
|
|
SUBROUTINE READTG(IDEVT,I,LPOS,LLEN,LCOM,LTYPE,NEXT)
|
|
INTEGER SWAPBO
|
|
EXTERNAL SWAPBO
|
|
IF(IDEVT.GT.0) THEN
|
|
READ(IDEVT,REC=I)LPOS,LLEN,LCOM,LTYPE,NEXT
|
|
LPOS = SWAPBO(LPOS)
|
|
LLEN = SWAPBO(LLEN)
|
|
LCOM = SWAPBO(LCOM)
|
|
NEXT = SWAPBO(NEXT)
|
|
ENDIF
|
|
END
|
|
SUBROUTINE WRITTG(IDEVT,I,LPOS,LLEN,LCOM,LTYPE,NEXT)
|
|
INTEGER SWAPBO
|
|
EXTERNAL SWAPBO
|
|
IF (IDEVT.GT.0) THEN
|
|
WRITE(IDEVT,REC=I)SWAPBO(LPOS),SWAPBO(LLEN),
|
|
+SWAPBO(LCOM),LTYPE,SWAPBO(NEXT)
|
|
ENDIF
|
|
END
|
|
SUBROUTINE READCC(IDEVC,I,ICNT,NEXT,NOTE)
|
|
C COMMENT_LENGTH
|
|
CHARACTER NOTE*40
|
|
C COMMENT_LENGTH - 4
|
|
CHARACTER DUMM*36
|
|
INTEGER SWAPBO
|
|
EXTERNAL SWAPBO
|
|
IF(IDEVC.GT.0)THEN
|
|
READ(IDEVC,REC=1)NEXT,ICNT,DUMM
|
|
NEXT = SWAPBO(NEXT)
|
|
ICNT = SWAPBO(ICNT)
|
|
IF(I.EQ.0.OR.I.GT.ICNT)THEN
|
|
NEXT = 0
|
|
NOTE = ' '
|
|
ELSE
|
|
READ(IDEVC,REC=I)NEXT,NOTE
|
|
NEXT = SWAPBO(NEXT)
|
|
ENDIF
|
|
ENDIF
|
|
END
|
|
SUBROUTINE WRITCC(IDEVC,I,ICNT,NEXT,NOTE)
|
|
C COMMENT_LENGTH
|
|
CHARACTER NOTE*40
|
|
C COMMENT_LENGTH - 4
|
|
CHARACTER DUMM*36
|
|
INTEGER SWAPBO
|
|
EXTERNAL SWAPBO
|
|
IF(IDEVC.GT.0)THEN
|
|
IF(I.EQ.1) THEN
|
|
WRITE(IDEVC,REC=1)SWAPBO(NEXT),SWAPBO(ICNT),DUMM
|
|
ELSE
|
|
READ(IDEVC,REC=1)IDUM,ICNT,DUMM
|
|
ICNT = SWAPBO(ICNT)
|
|
IF(I.LE.ICNT) THEN
|
|
WRITE(IDEVC,REC=I)SWAPBO(NEXT),NOTE
|
|
ENDIF
|
|
ENDIF
|
|
ENDIF
|
|
END
|
|
SUBROUTINE PADCON(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+GEL,LINCON,POSN,NC,IDBSIZ,IDEVR,IDEVW,MAXGEL,KBOUT)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),POSN,X
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
CHARACTER GEL(MAXGEL)
|
|
CHARACTER PAD
|
|
SAVE PAD
|
|
DATA PAD/'*'/
|
|
C NOW FIND FIRST CHAR THAT OVERLAPS REGION
|
|
LLINO=LNBR(LINCON)
|
|
30 CONTINUE
|
|
X=RELPG(LLINO)+ABS(LNGTHG(LLINO))-1
|
|
IF(X.GE.POSN)GO TO 40
|
|
C NOT IN REGION
|
|
LLINO=RNBR(LLINO)
|
|
GO TO 30
|
|
40 CONTINUE
|
|
C NOW GET THIS GEL FROM DISK
|
|
CALL READW(IDEVW,LLINO,GEL,MAXGEL)
|
|
C CALC POSN IN THIS GEL TO EDIT
|
|
X=POSN-RELPG(LLINO)+1
|
|
K=X
|
|
C MOVE THE DATA RIGHT
|
|
M=ABS(LNGTHG(LLINO))
|
|
LNGTHG(LLINO)=LNGTHG(LLINO)+SIGN(NC,LNGTHG(LLINO))
|
|
C CHECK FOR OVER END OF ARRAY
|
|
N=ABS(LNGTHG(LLINO))
|
|
IF(N.GT.MAXGEL)THEN
|
|
WRITE(KBOUT,1000)LLINO
|
|
1000 FORMAT(
|
|
+' Data pushed off end of gel',I4,' during padding')
|
|
M=M-(N-MAXGEL)
|
|
N=MAXGEL
|
|
LNGTHG(LLINO)=SIGN(MAXGEL,LNGTHG(LLINO))
|
|
END IF
|
|
J=M-K+1
|
|
DO 55 I=1,J
|
|
GEL(N)=GEL(M)
|
|
N=N-1
|
|
M=M-1
|
|
55 CONTINUE
|
|
C PERFORM THE INSERTION
|
|
DO 60 I=K,MIN(MAXGEL,K+NC-1)
|
|
GEL(I)=PAD
|
|
60 CONTINUE
|
|
C WRITE BACK TO DISK
|
|
CALL WRITEW(IDEVW,LLINO,GEL,MAXGEL)
|
|
C WRITE NEW LINE
|
|
CALL WRITER(IDEVR,LLINO,RELPG(LLINO),LNGTHG(LLINO),
|
|
+LNBR(LLINO),RNBR(LLINO))
|
|
C NOW UPDATE TAG FILES ACCORDINGLY
|
|
CALL PADTAG(LLINO,K,NC,LNGTHG(LLINO))
|
|
65 CONTINUE
|
|
C NOW GET NEXT GEL
|
|
LLINO=RNBR(LLINO)
|
|
C LAST GEL?
|
|
IF(LLINO.EQ.0)GO TO 70
|
|
C DOES IT HAVE DATA IN REGION?
|
|
C IE DO RELPG AND RELPG+LNGTHG-1 LIE EITHER SIDE OF POSN?
|
|
IF(RELPG(LLINO).GT.POSN)GO TO 70
|
|
X=RELPG(LLINO)+ABS(LNGTHG(LLINO))-1
|
|
IF(X.LT.POSN)GO TO 65
|
|
C WITHIN
|
|
GO TO 40
|
|
70 CONTINUE
|
|
C INSERTS FINISHED SO NEED TO INCREMENT ALL THOSE GELS TO RIGHT
|
|
LLINO=LNBR(LINCON)
|
|
75 CONTINUE
|
|
IF(RELPG(LLINO).GT.POSN)GO TO 80
|
|
76 CONTINUE
|
|
LLINO=RNBR(LLINO)
|
|
IF(LLINO.EQ.0)GO TO 90
|
|
GO TO 75
|
|
80 CONTINUE
|
|
RELPG(LLINO)=RELPG(LLINO)+NC
|
|
C WRITE NEW LINE
|
|
CALL WRITER(IDEVR,LLINO,RELPG(LLINO),LNGTHG(LLINO),
|
|
+LNBR(LLINO),RNBR(LLINO))
|
|
GO TO 76
|
|
90 CONTINUE
|
|
C NEED TO INCREMENT CONTIG LINE
|
|
RELPG(LINCON)=RELPG(LINCON)+NC
|
|
CALL WRITER(IDEVR,LINCON,RELPG(LINCON),LNGTHG(LINCON),
|
|
+LNBR(LINCON),RNBR(LINCON))
|
|
END
|
|
SUBROUTINE AUTOJ(RELPG,LNGTHG,LNBR,RNBR,MAXDB,IDBSIZ,
|
|
+NGELS,NCONTS,MAXGEL,
|
|
+TEMP3,WORDP,WORDN,LPOWRC,POSNS,GELN,
|
|
+SEQ1,MAXSEQ,SEQ2,SEQ3,SEQ4,SEQ5,SEQC2,SEQG2,MATCH,
|
|
+MAXGLM,MAXGL2,CHRSIZ,ECHRSZ,LENGTH,
|
|
+SAV1,SAV2,SAV3,MAXSAV,CENDS,NENDS,MAXCON,CONST,
|
|
+KBIN,KBOUT,IDEV1,IDEV2,IDEV3,IDEV4,IDEV7,IDEV8,IDEV,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,NAMARC,NAMPRO,FILE,
|
|
+PERCD,IOPEN,IDM,SEQG3,SEQC3,IOK,IDEVC,IDEVT)
|
|
INTEGER CHRSIZ,ECHRSZ
|
|
INTEGER RELPG(MAXDB)
|
|
INTEGER LNGTHG(MAXDB),LNBR(MAXDB),RNBR(MAXDB)
|
|
INTEGER JOINT(2),ITOTPC(2),ITOTPG(2),IDIM22(2),IDOUT(2)
|
|
INTEGER LLINO(2),ITYPE(2),IFAIL
|
|
INTEGER ILEFTS(2),ILC(2),IPOSC(2),IPOSG(2),ISENSE(2)
|
|
INTEGER WINDOW
|
|
INTEGER TEMP3(ECHRSZ,MAXGL2),CONST(LENGTH)
|
|
INTEGER POSNS(MAXSEQ),WORDP(LPOWRC),WORDN(LPOWRC),GELN(MAXGLM)
|
|
INTEGER CENDS(MAXCON),NENDS(MAXCON)
|
|
CHARACTER SEQ3(MAXGLM),SEQC2(MAXGLM,2),SEQG2(MAXGLM,2)
|
|
CHARACTER SEQ1(MAXSEQ),SEQ2(MAXGLM),MATCH(MAXGLM),SEQ4(MAXGLM)
|
|
INTEGER SAV1(MAXSAV),SAV2(MAXSAV),SAV3(MAXSAV)
|
|
CHARACTER NAMARC*(*),NAMPRO*(*),FILE*(*)
|
|
CHARACTER SEQ5(MAXGLM),HELPF*(*),SEQG3(MAXGLM),SEQC3(MAXGLM)
|
|
CALL DBCHEK(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ TEMP3,IERR,KBOUT)
|
|
IF(IERR.GT.1) RETURN
|
|
IFAIL = 0
|
|
IF(NGELS.LT.1) RETURN
|
|
MN = LENGTH*2
|
|
MX = MAXGLM + 1
|
|
MINMAT = MAX(15,MN)
|
|
CALL GETINT(MN,MX,MINMAT,
|
|
+'Minimum initial match',
|
|
+IVAL,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
MINMAT = IVAL
|
|
MINSLI = 3
|
|
MN = 0
|
|
MX = 50
|
|
MAXPG = 8
|
|
CALL GETINT(MN,MX,MAXPG,
|
|
+'Maximum pads per sequence',
|
|
+IVAL,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
MAXPG = IVAL
|
|
MAXPC = IVAL
|
|
RMN = 0.
|
|
RMX = 100.
|
|
PERMAX = 75.
|
|
CALL GETRL(RMN,RMX,PERMAX,
|
|
+ 'Maximum percent mismatch after alignment',
|
|
+ VAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
PERMAX = VAL
|
|
MN = MINMAT
|
|
MX = MAXGEL
|
|
WINDOW = 100
|
|
CALL GETINT(MN,MX,WINDOW,
|
|
+'Probe length',
|
|
+IVAL,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
WINDOW = IVAL
|
|
IOK = 1
|
|
I = 0
|
|
CALL YESNO(I,'Use clipped data',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(I.LT.0) RETURN
|
|
IWING = 0
|
|
IF(I.EQ.0) THEN
|
|
MN = 1
|
|
MX = MAXGEL
|
|
IWING = 100
|
|
CALL GETINT(MN,MX,IWING,
|
|
+ 'Window size for good data scan',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
IWING = IVAL
|
|
MN = 1
|
|
MX = MIN(100,IWING)
|
|
NBAD = MIN(IWING,5)
|
|
CALL GETINT(MN,MX,NBAD,
|
|
+ 'Maximum number of dashes in scan window',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
NBAD = IVAL
|
|
END IF
|
|
MAXOVR=MAXGEL-3*MAX(MAXPC,MAXPG)
|
|
IDIM2=MAXGEL
|
|
C
|
|
C Note I am doing something horrible here to save space:
|
|
C temp3 is used for cends, nends and as temp3. I can do this
|
|
C because i call fndcon to reinstate cends ,nends after each
|
|
C use of temp3 for the consensus calc. However iladd and iradd
|
|
C need the same amount of storage as cends and nends so i use
|
|
C the arrays cends and nends for them. iladd and iradd are not
|
|
C updated after use of temp3. The switch in array use is made
|
|
C by the following subroutine call
|
|
C
|
|
C replaced following line 25-9-92
|
|
C +SAV1,SAV2,SAV3,MAXSAV,TEMP3,TEMP3(MAXCON+1,1),MAXCON,
|
|
C
|
|
CALL AUTOJN(SEQ1,MAXSEQ,SEQ2,IDIM2,ILEFTS,ILC,IPOSC,
|
|
+IPOSG,ISENSE,LLINO,IMATC,IFCOMP,MINMAT,POSNS,WORDP,WORDN,
|
|
+CONST,LENGTH,LPOWRC,IDEV,MATCH,MAXGEL,MAXGLM,SEQ5,GELN,
|
|
+SAV1,SAV2,SAV3,MAXSAV,TEMP3,TEMP3(1,1+MAXCON/2),MAXCON,
|
|
+SEQG2,SEQC2,SEQ4,IDOUT,IDIM22,ITOTPG,ITOTPC,JOINT,IFAIL,
|
|
+ITYPE,MAXPC,MAXPG,PERMAX,MINSLI,SEQG3,SEQC3,KFAIL,
|
|
+WINDOW,CENDS,NENDS,RELPG,LNBR,IDBSIZ,NGELS,NCONTS,
|
|
+LNGTHG,RNBR,IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,PERCD,IWING,NBAD,
|
|
+ECHRSZ,SEQ5,TEMP3,MAXGL2,IDM,NAMPRO,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
END
|
|
C
|
|
SUBROUTINE AUTOJN(SEQ1,MAXSEQ,GEL,IDIMGI,ILEFTS,
|
|
+ILC,IPOSC,
|
|
+IPOSG,ISENSE,LLINO,IMATC,IFCOMP,MINMAT,POSNS,WORDP,WORDN,
|
|
+CONST,LENGTH,LPOWRC,IDEV,MATCH,MAXGEL,MAXGLM,GELCOP,GELN,
|
|
+SAVPS,SAVPG,SAVL,MAXSAV,CENDS,NENDS,MAXCON,
|
|
+SEQG2,SEQC2,SEQ4,IDOUT,IDIM22,ITOTPG,ITOTPC,JOINT,IFAIL,
|
|
+ITYPE,MAXPC,MAXPG,PERMAX,MINSLI,SEQG3,SEQC3,KFAIL,
|
|
+WINDOW,ILADD,IRADD,RELPG,LNBR,IDBSIZ,NGELS,NCONTS,
|
|
+LNGTHG,RNBR,IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,PERCD,IWING,NBAD,
|
|
+ECHRSZ,SEQ5,TEMP3,MAXGL2,IDM,NAMPRO,KBIN,KBOUT,
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER ECHRSZ
|
|
INTEGER ILEFTS(2),ILC(2),IPOSC(2),IPOSG(2),ISENSE(2),LLINO(2)
|
|
INTEGER POSNS(MAXSEQ),GELN(MAXGLM),WORDP(LPOWRC),SAVPS(MAXSAV)
|
|
INTEGER SAVPG(MAXSAV),SAVL(MAXSAV)
|
|
INTEGER WORDN(LPOWRC)
|
|
CHARACTER GELCOP(MAXGLM),MATCH(MAXGLM),NAMPRO*(*),HELPF*(*)
|
|
INTEGER CENDS(MAXCON),NENDS(MAXCON),ILADD(MAXCON),IRADD(MAXCON)
|
|
INTEGER CONST(LENGTH)
|
|
CHARACTER SEQ1(MAXSEQ),GEL(MAXGLM),SEQ5(MAXGLM)
|
|
CHARACTER SEQG2(MAXGLM,2),SEQC2(MAXGLM,2),SEQ4(MAXGLM)
|
|
INTEGER IDOUT(2),IDIM22(2),ITOTPG(2),ITOTPC(2),JOINT(2)
|
|
INTEGER IFAIL,ITYPE(2)
|
|
PARAMETER (MAXC = 100)
|
|
CHARACTER SEQG3(MAXGLM),SEQC3(MAXGLM)
|
|
INTEGER JLEFTS(MAXC),JLC(MAXC),JPOSC(MAXC),JPOSG(MAXC),MCON(MAXC)
|
|
INTEGER JSENSE(MAXC),JLLINO(MAXC),WINDOW,CSTART
|
|
INTEGER RELPG(IDBSIZ),LNBR(IDBSIZ),LNGTHG(IDBSIZ),RNBR(IDBSIZ)
|
|
INTEGER TEMP3(ECHRSZ,MAXGL2)
|
|
INTEGER CLINNO,CHAINL,XVERSN
|
|
EXTERNAL CLINNO,CHAINL,XVERSN
|
|
IEXT = 0
|
|
IF (IWING.NE.0) IEXT = 1
|
|
IFAIL = 1
|
|
KFAIL = 0
|
|
C
|
|
C calc the consensus and add the unused data if required
|
|
C save the lengths of the unused data in iladd, iradd
|
|
C
|
|
CALL FILLI(ILADD,NCONTS,0)
|
|
CALL FILLI(IRADD,NCONTS,0)
|
|
IDIM = 0
|
|
JOB = 0
|
|
CALL JCONS(RELPG,LNGTHG,LNBR,RNBR,NAMPRO,NGELS,NCONTS,
|
|
+SEQ1,MAXSEQ,GEL,IDBSIZ,IDIM,JOB,KDUMM,KDUMM,KDUMM,KDUMM,
|
|
+TEMP3,
|
|
+ECHRSZ,MAXGL2,KBOUT,IDEV2,IFAIL,MAXGEL,IDM,PERCD,SEQ5,
|
|
+ILADD,IRADD,MAXCON,KDUMM,IWING,NBAD)
|
|
IDCEND=MAXCON
|
|
CALL FNDCON(SEQ1,IDIM,CENDS,NENDS,IDCEND,MAXCON,KBOUT)
|
|
C WRITE(*,*)'IDIM',IDIM
|
|
C WRITE(*,*)(CENDS(K),NENDS(K),K=1,IDCEND)
|
|
C
|
|
C find possible missed joins
|
|
C we have consensus in seq1 in order first contig,second contig etc
|
|
C compare the ends in reverse order, simultaneously shortening the consensus
|
|
C the total consensus length is denoted by idimin
|
|
C but we compare the last contig to the rest so denote the length
|
|
C of the rest by idim
|
|
C
|
|
NFOUND = 0
|
|
NMADE = 0
|
|
C
|
|
C init hashing constants
|
|
C
|
|
CALL INITE(CONST,CSTART,LENGTH)
|
|
C
|
|
C idim is the length up to the end of the last but 1 consensus
|
|
C (we get from the start of the last contig)
|
|
C jcon is the last contig (ie the one we are currently comparing)
|
|
C
|
|
JCON = IDCEND + 1
|
|
C
|
|
C come here for the next contig
|
|
C
|
|
10 CONTINUE
|
|
C
|
|
C
|
|
IDIMIN = IDIM
|
|
C WRITE(*,*)'IDIM,IDIMIN',IDIM,IDIMIN
|
|
JCON = JCON - 1
|
|
C
|
|
C come here if weve just made a join
|
|
C
|
|
20 CONTINUE
|
|
C
|
|
C
|
|
C WRITE(*,*)'*********COMPARING CONTIG',JCON
|
|
IF (JCON.GT.1) THEN
|
|
C WRITE(*,*)'IDCEND,JCON'
|
|
C WRITE(*,*)IDCEND,JCON
|
|
C WRITE(*,1222)
|
|
C + (K,NENDS(K),CENDS(K),ILADD(K),IRADD(K),K=1,JCON)
|
|
C 1222 FORMAT(5I6)
|
|
C WRITE(*,*)(ILADD(K),IRADD(K),K=1,JCON)
|
|
IDIM = CENDS(JCON) - 1
|
|
C WRITE(*,*)'IDIM,IDIMIN',IDIM,IDIMIN
|
|
CALL ENCOF(SEQ1,IDIM,CONST,CSTART,LENGTH,POSNS)
|
|
CALL ENCONN(POSNS,IDIM,WORDP,WORDN,LPOWRC,LENGTH,1)
|
|
C
|
|
C point to probes at the ends of the current contig
|
|
C
|
|
JS = CENDS(JCON) + 20
|
|
JE = CENDS(JCON+1) - 1
|
|
IEND = 1
|
|
IDIMG = MIN(WINDOW,JE-JS+1)
|
|
C
|
|
C check for case where contig is shorter than probe (window)
|
|
C in which case only compare the left hand end
|
|
C
|
|
IF(JE-JS+1.LE.WINDOW) IEND = 2
|
|
C
|
|
C come here for left end of contig
|
|
C
|
|
30 CONTINUE
|
|
C
|
|
C
|
|
C WRITE(*,*)'JCON,JS,JE,IDIMG'
|
|
C WRITE(*,*)JCON,JS,JE,IDIMG
|
|
IMATC = 0
|
|
IF(IEND.EQ.1) THEN
|
|
CALL SQCOPY(SEQ1(JE-IDIMG+1),GEL,IDIMG)
|
|
ELSE
|
|
CALL SQCOPY(SEQ1(JS),GEL,IDIMG)
|
|
END IF
|
|
CALL SQCOPY(GEL,GELCOP,IDIMG)
|
|
ISTRAN=1
|
|
C
|
|
C come here for strand 2
|
|
C
|
|
40 CONTINUE
|
|
C
|
|
C
|
|
CALL ENCOF(GEL,IDIMG,CONST,CSTART,LENGTH,GELN)
|
|
IDSAV=MAXSAV
|
|
CALL CFGEL(GELN,IDIMG,POSNS,IDIM,WORDP,WORDN,LENGTH,LPOWRC,
|
|
+ SAVPG,SAVPS,SAVL,
|
|
+ IDSAV,SEQ1,GEL,MINMAT,IFCOMP,IDEV)
|
|
IF(IFCOMP.NE.0) THEN
|
|
CALL ERROM(KBOUT,'Error in CFGEL')
|
|
RETURN
|
|
END IF
|
|
IF(IDSAV.NE.0)THEN
|
|
CALL ADISM5(IDIM,IDIMG,SAVPS,SAVPG,IDSAV,CENDS,NENDS,
|
|
+ IDCEND,MAXCON,JLEFTS,JLC,JPOSC,JPOSG,JSENSE,JLLINO,
|
|
+ IMATC,ISTRAN,MAXC,KBOUT,MCON)
|
|
END IF
|
|
ISTRAN=ISTRAN+1
|
|
C
|
|
C have we done both strands for this probe ?
|
|
C
|
|
IF(ISTRAN.EQ.2) THEN
|
|
CALL SQCOPY(GELCOP,GEL,IDIMG)
|
|
CALL SQREV(GEL,IDIMG)
|
|
CALL SQCOM(GEL,IDIMG)
|
|
GO TO 40
|
|
END IF
|
|
CALL SQCOPY(GELCOP,GEL,IDIMG)
|
|
KSENSE = 0
|
|
C WRITE(KBOUT,*)'Total matches found',IMATC
|
|
C CALL FMTDB(SEQ1,IDIMIN,1,IDIMIN,60,KBOUT)
|
|
IF(IMATC.NE.0) THEN
|
|
JMATC = 0
|
|
C
|
|
C loop for all matches on both strands
|
|
C
|
|
DO 100 I = 1,IMATC
|
|
IF(JSENSE(I).EQ.-1) THEN
|
|
IF(KSENSE.EQ.0) THEN
|
|
CALL SQREV(GEL,IDIMG)
|
|
CALL SQCOM(GEL,IDIMG)
|
|
KSENSE = 1
|
|
END IF
|
|
END IF
|
|
JDIM22 = IDIMG
|
|
JDOUT = MAXGEL
|
|
IDSAV = MAXSAV
|
|
CALL ALINEJ(SEQ1(JLEFTS(I)),GEL,SEQG3,SEQC3,
|
|
+ SAVPS,SAVPG,SAVL,IDSAV,JLC(I),JDIM22,JDOUT,
|
|
+ JPOSC(I),JPOSG(I),MINSLI,JJOINT,JTOTPC,JTOTPG,
|
|
+ JFAIL,JTYPE,MAXPC,MAXPG,PERMAX,KBOUT,IDEV,
|
|
+ SEQ4,MAXGEL,PERMS,
|
|
+ NENDS(JCON),IEND,JLLINO(I),JSENSE(I),ILADD,IRADD,MAXCON,
|
|
+ MCON(I),
|
|
+ JCON,RELPG,LNBR,IDBSIZ,NCONTS,NFOUND,IENDCT,IENDGT)
|
|
IF (JFAIL.NE.0) GO TO 100
|
|
C
|
|
C if this is the x version and the results are being shown on the screen
|
|
C we can offer the join editor and make the join
|
|
C
|
|
IF ((XVERSN().NE.0).AND.(IDEV.EQ.KBOUT)) THEN
|
|
C
|
|
C another horrible thing here: sending posns for use by dbchek after
|
|
C a join is made (only used then)
|
|
C
|
|
CALL USJED(
|
|
+ SEQ1,MAXSEQ,IDIMIN,IDIM,GEL,KBIN,KBOUT,MAXGEL,MAXGLM,
|
|
+ CENDS,NENDS,MAXCON,IFAIL,MCON(I),JCON,NAMPRO,
|
|
+ JSENSE(I), JLLINO(I),
|
|
+ ILADD,IRADD,RELPG,LNBR,IDBSIZ,NGELS,
|
|
+ NCONTS,
|
|
+ LNGTHG,RNBR,IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,
|
|
+ PERCD,IWING,NBAD,
|
|
+ IENDGT,IENDCT,SEQ4,IDCEND,IEXT,
|
|
+ ECHRSZ,SEQ5,TEMP3,MAXGL2,IDM,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,POSNS,NMADE,IOK)
|
|
IF (IFAIL.NE.0) GO TO 200
|
|
IF (NCONTS.EQ.1) GO TO 200
|
|
IF (IOK.EQ.2) GO TO 200
|
|
C
|
|
C if join has been made, assume everything is tidy and do this contig again
|
|
C
|
|
IF (IOK.EQ.0) GO TO 20
|
|
END IF
|
|
100 CONTINUE
|
|
END IF
|
|
IEND = IEND + 1
|
|
C
|
|
C have we done both ends of this contig ?
|
|
C
|
|
IF(IEND.EQ.2) GO TO 30
|
|
C
|
|
C do the next contig
|
|
C
|
|
GO TO 10
|
|
END IF
|
|
C
|
|
C jump here if only one contig left ! or user wants to escape
|
|
C
|
|
200 CONTINUE
|
|
WRITE(KBOUT,*)'Number of potential joins found',NFOUND
|
|
WRITE(KBOUT,*)'Number of joins made',NMADE
|
|
END
|
|
SUBROUTINE USJED(
|
|
+SEQ1,MAXSEQ,IDIMIN,IDIM,GEL,KBIN,KBOUT,MAXGEL,MAXGLM,
|
|
+CENDS,NENDS,MAXCON,IFAIL,MCON,JCON,NAMPRO,JSENSE,JLLINO,
|
|
+ILADD,IRADD,RELPG,LNBR,IDBSIZ,NGELS,NCONTS,
|
|
+LNGTHG,RNBR,IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,PERCD,IWING,NBAD,
|
|
+IENDGT,IENDCT,SEQ2,IDCEND,IEXT,
|
|
+ECHRSZ,SEQ5,TEMP3,MAXGL2,IDM,IHELPS,IHELPE,HELPF,IDEVH,
|
|
+TEMP1,NMADE,IOK)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER ECHRSZ
|
|
INTEGER CENDS(MAXCON),NENDS(MAXCON),ILADD(MAXCON),IRADD(MAXCON)
|
|
CHARACTER SEQ1(MAXSEQ),GEL(MAXGLM),SEQ5(MAXGLM),SEQ2(MAXGEL)
|
|
INTEGER JSENSE,JLLINO,MCON,TEMP1(IDBSIZ)
|
|
INTEGER RELPG(IDBSIZ),LNBR(IDBSIZ),LNGTHG(IDBSIZ),RNBR(IDBSIZ)
|
|
INTEGER TEMP3(ECHRSZ,MAXGL2)
|
|
INTEGER CLINNO,CHAINL
|
|
CHARACTER HELPF*(*),NAMPRO*(*)
|
|
EXTERNAL CLINNO,CHAINL
|
|
C WRITE(*,*)'CURRENT MAX CONSENSUS LENGTH',MAXSEQ
|
|
C WRITE(*,*)'CURRENT CONSENSUS LENGTH',IDIM
|
|
C WRITE(*,*)'MAXGEL,MAXGLM,MAXCON',MAXGEL,MAXGLM,MAXCON
|
|
C WRITE(*,*)'NGELS,NCONTS',NGELS,NCONTS
|
|
C
|
|
C routine to offer use of join editor and if join is made tidy up
|
|
C
|
|
IFAIL = 0
|
|
IOK = 2
|
|
IJOIN = 0
|
|
CALL YESNO(IJOIN,'Use the editor',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(IJOIN.LT.0) RETURN
|
|
IOK = 1
|
|
IF (IJOIN.EQ.0) THEN
|
|
C
|
|
C weve called the editor
|
|
C
|
|
C
|
|
C the contig names are nends(jcon) and jllino(i)
|
|
C if jsense is -1 the contig to complement is nends(jcon)
|
|
C the position of the overlap is iendct in jllino(i) and iendgt in nends(jcon)
|
|
C we are going to need to know the endpoints of the consensuses to delete
|
|
C the element number in nends of one contig is mcon(i) and the other is jcon
|
|
C cends(i) stores the position relative to the consensus of the left end
|
|
C of contig i
|
|
C
|
|
C WRITE(*,*)'NENDS(JCON)',NENDS(JCON)
|
|
C WRITE(*,*)'IDBSIZ,NCONTS',IDBSIZ,NCONTS
|
|
IF (JSENSE.EQ.-1) THEN
|
|
LINCON = CLINNO(LNBR,IDBSIZ,NCONTS,NENDS(JCON))
|
|
IF (LINCON.EQ.0) THEN
|
|
C
|
|
C major cockup here: weve lost the contig number
|
|
C
|
|
CALL ERROM(KBOUT,
|
|
+ 'Error: contig line not found for complementing')
|
|
IFAIL = 1
|
|
RETURN
|
|
END IF
|
|
CALL CMPLMT(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ LINCON,NENDS(JCON),GEL,IDBSIZ,KBOUT,IDEV1,IDEV2,MAXGEL)
|
|
NENDS(JCON) = LNBR(LINCON)
|
|
END IF
|
|
C
|
|
C note the number of contigs
|
|
C
|
|
NCIN = NCONTS
|
|
C WRITE(*,*)'JCON,NENDS(JCON)',JCON,NENDS(JCON)
|
|
LNCONL = CLINNO(LNBR,IDBSIZ,NCONTS,NENDS(JCON))
|
|
IF (LNCONL.EQ.0) THEN
|
|
C
|
|
C major cockup here: weve lost the contig number
|
|
C
|
|
CALL ERROM(KBOUT,
|
|
+ 'Error: contig line not found for left contig')
|
|
IFAIL = 1
|
|
RETURN
|
|
END IF
|
|
C WRITE(*,*)'I,MCON,NENDS(MCON)',I,MCON,NENDS(MCON)
|
|
LNCONR = CLINNO(LNBR,IDBSIZ,NCONTS,NENDS(MCON))
|
|
IF (LNCONR.EQ.0) THEN
|
|
C
|
|
C major cockup here: weve lost the contig number
|
|
C
|
|
CALL ERROM(KBOUT,
|
|
+ 'Error: contig line not found for right contig')
|
|
IFAIL = 1
|
|
RETURN
|
|
END IF
|
|
C WRITE(*,*)IDBSIZ,LNCONL,NENDS(JCON),LNCONR,NENDS(MCON)
|
|
C WRITE(*,*)IENDGT,IENDCT,PERCD,NGELS,NCONTS
|
|
C
|
|
C call the editor
|
|
C
|
|
CALL JXEDIT(IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,
|
|
+ RELPG,LNGTHG,LNBR,RNBR,MAXGEL,
|
|
+ IDBSIZ,LNCONR,NENDS(MCON),LNCONL,NENDS(JCON),
|
|
+ 0,IENDCT,0,IENDGT,PERCD,NGELS,NCONTS,IDM,IEXT,JOK)
|
|
C
|
|
C jok 0 = nothing done
|
|
C 1 = left contig saved
|
|
C 2 = right contig saved
|
|
C 3 = both contigs saved
|
|
C 7 = join made
|
|
C we could use these values but we dont: we only count contigs.
|
|
C if a contig is saved its consensus may be different. For now dont worry!!!!
|
|
C
|
|
C
|
|
C note we may not have joined but may have either
|
|
C complemented a contig or edited one! We shall ignore the edited
|
|
C contig, but must uncomplement the contig. See below.
|
|
C Remove the 2 consensus sequences
|
|
C Add the new one to the right end
|
|
C Update nends, cends
|
|
C Update iladd iradd, and add the new data to the end of the list
|
|
C to update iladd, iradd move all data left one element
|
|
C
|
|
IF (NCIN.GT.NCONTS) THEN
|
|
C
|
|
C weve made a join
|
|
C
|
|
NMADE = NMADE + 1
|
|
C
|
|
C check for consistency
|
|
C
|
|
CALL DBCHEK(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ TEMP1,IFAIL,KBOUT)
|
|
IF (IFAIL.NE.0) RETURN
|
|
C
|
|
C sometimes there may only be one contig now!
|
|
C
|
|
IF (NCONTS.EQ.1) RETURN
|
|
C
|
|
C get contig left gel number
|
|
C
|
|
LLINOC = CHAINL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ IDBSIZ,NENDS(JCON))
|
|
IF (LLINOC.EQ.0) THEN
|
|
CALL ERROM(KBOUT,'Gel number not in a contig')
|
|
IFAIL = 1
|
|
RETURN
|
|
END IF
|
|
C
|
|
C get contig line number
|
|
C
|
|
LINCON = CLINNO(LNBR,IDBSIZ,NCONTS,LLINOC)
|
|
IF (LINCON.EQ.0) THEN
|
|
CALL ERROM(KBOUT,'Error: contig line not found')
|
|
IFAIL = 1
|
|
RETURN
|
|
END IF
|
|
C WRITE(*,*)'DELETE CONSENSUS FROM',CENDS(MCON)
|
|
C WRITE(*,*)'TO',CENDS(MCON+1)
|
|
C WRITE(*,*)'SHORTEN THE CONSENSUS BY'
|
|
C WRITE(*,*)CENDS(MCON+1)-CENDS(MCON)
|
|
C WRITE(*,*)CENDS(JCON+1)-CENDS(JCON)
|
|
C
|
|
C shift the consensus to delete the internal consensus then
|
|
C reduce idim,idimin accordingly (for both contigs). Finally
|
|
C move the extension sizes into the appropriate elements.
|
|
C
|
|
CALL SHFTLA(SEQ1,MAXSEQ,CENDS(MCON+1),CENDS(MCON),
|
|
+ IDIMIN)
|
|
IDIM = IDIM - (CENDS(JCON+1)-CENDS(JCON))
|
|
IDIMIN = IDIMIN - (CENDS(JCON+1)-CENDS(JCON))
|
|
IDIM = IDIM - (CENDS(MCON+1)-CENDS(MCON))
|
|
IDIMIN = IDIMIN - (CENDS(MCON+1)-CENDS(MCON))
|
|
DO 50 J=MCON,NCONTS-1
|
|
ILADD(J) = ILADD(J+1)
|
|
IRADD(J) = IRADD(J+1)
|
|
50 CONTINUE
|
|
C
|
|
C make a consensus for the joined contig, putting it at the end
|
|
C
|
|
JOB = 1
|
|
JCON = JCON - 1
|
|
CALL JCONS(RELPG,LNGTHG,LNBR,RNBR,NAMPRO,NGELS,NCONTS,
|
|
+ SEQ1,MAXSEQ,SEQ2,IDBSIZ,IDIMIN,JOB,LLINOC,LINCON,1,
|
|
+ RELPG(LINCON),TEMP3,
|
|
+ ECHRSZ,MAXGL2,KBOUT,IDEV2,IFAIL,MAXGEL,IDM,PERCD,SEQ5,
|
|
+ ILADD,IRADD,MAXCON,JCON,IWING,NBAD)
|
|
IF (IFAIL.NE.0) RETURN
|
|
C
|
|
C now restore the contig end positions because jcons has overwritten them
|
|
C (they have changed anyway)
|
|
C
|
|
CALL FNDCON(SEQ1,IDIMIN,CENDS,NENDS,IDCEND,MAXCON,KBOUT)
|
|
IOK = 0
|
|
C
|
|
C end of join made
|
|
C
|
|
ELSE
|
|
C
|
|
C if we get here we have called the editor but not made a join
|
|
C so we might need to tidy up
|
|
C
|
|
IF (JSENSE.EQ.-1) THEN
|
|
LINCON = CLINNO(LNBR,IDBSIZ,NCONTS,NENDS(JCON))
|
|
IF (LINCON.EQ.0) THEN
|
|
C
|
|
C major cockup here: weve lost the contig number
|
|
C
|
|
CALL ERROM(KBOUT,
|
|
+ 'Error: contig line not found for complementing')
|
|
IFAIL = 1
|
|
RETURN
|
|
END IF
|
|
CALL CMPLMT(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ LINCON,NENDS(JCON),GEL,IDBSIZ,KBOUT,IDEV1,IDEV2,
|
|
+ MAXGEL)
|
|
NENDS(JCON) = LNBR(LINCON)
|
|
END IF
|
|
END IF
|
|
END IF
|
|
END
|
|
SUBROUTINE ALINEJ(SEQ1,SEQ2,SEQG2,SEQC2,ISAV1,ISAV2,ISAV3,
|
|
+IDSAV,IDC,IDIM2,IDOUT,IC1,IG1,MINSLI,JOINT,
|
|
+ITOTPC,ITOTPG,IFAIL,ITYPE,MAXPC,MAXPG,PERMAX,KBOUT,IDEV,
|
|
+SEQ3,MAXGEL,
|
|
+PERCM,JCONN,IEND,NCON,JSENSE,ILADD,IRADD,MAXCON,MCON,JCON,
|
|
+RELPG,LNBR,IDBSIZ,NCONTS,NFOUND,IENDCT,IENDGT)
|
|
C AUTHOR: RODGER STADEN
|
|
CHARACTER SEQ1(IDC),SEQ2(IDIM2),SEQG2(IDOUT),SEQC2(IDOUT)
|
|
CHARACTER SEQ3(MAXGEL)
|
|
INTEGER ISAV1(IDSAV),ISAV2(IDSAV),ISAV3(IDSAV)
|
|
INTEGER ILADD(MAXCON),IRADD(MAXCON)
|
|
INTEGER RELPG(IDBSIZ),LNBR(IDBSIZ),WINDOW
|
|
MINSLT=MINSLI
|
|
C SAVE SEQ2
|
|
CALL SQCOPY(SEQ2,SEQ3,IDIM2)
|
|
IFAIL=1
|
|
C FIND MATCHES
|
|
IPP=IDSAV
|
|
CALL SLIDER(SEQ1,IDC,SEQ3,IDIM2,IC1,IG1,MAXPG,MAXPC,MINSLT,
|
|
+ISAV1,ISAV2,ISAV3,IPP)
|
|
IF(IPP.GT.IDSAV)RETURN
|
|
IF(IPP.LT.1)RETURN
|
|
CALL REMOVL(ISAV2,ISAV3,ISAV1,IPP)
|
|
CALL BUB3AS(ISAV2,ISAV3,ISAV1,IPP)
|
|
C DO TOPOLOGICAL CHECK
|
|
CALL TPCHEK(ISAV2,ISAV3,ISAV1,IPP)
|
|
IFAIL=0
|
|
C
|
|
C save idim2 as window
|
|
C
|
|
WINDOW = IDIM2
|
|
CALL LINEUP(SEQ2,SEQ1,SEQG2,SEQC2,IDC,IDIM2,IDOUT,ISAV3,ISAV2,
|
|
+ISAV1,IPP,ITOTPC,ITOTPG,JOINT,ITYPE,KBOUT,MAXGEL,IFAIL)
|
|
IF(IFAIL.NE.0)RETURN
|
|
C IDIM2 IS NOW LENGTH OF ALIGNED GEL
|
|
CALL JALIGN(SEQC2,SEQG2,SEQ3,MAXGEL,IDOUT,IDIM2,JOINT,
|
|
+ITYPE,PERCM,KBOUT,IDEV,IFAIL,PERMAX,JCONN,IEND,NCON,JSENSE,
|
|
+ILADD,IRADD,MAXCON,MCON,JCON,RELPG,LNBR,IDBSIZ,NCONTS,WINDOW,
|
|
+NFOUND,IENDCT,IENDGT)
|
|
IF (ITOTPC.GT.MAXPC) IFAIL = 1
|
|
IF (ITOTPG.GT.MAXPG) IFAIL = 1
|
|
IF (PERCM.GT.PERMAX) IFAIL = 1
|
|
END
|
|
C
|
|
C COUNTS MISMATCHES AND DISPLAYS OVERLAP.
|
|
SUBROUTINE JALIGN(SEQC2,SEQG2,SEQ3,MAXGEL,IDOUT,IDIM2,
|
|
+JOINT,ITYPE,X,KBOUT,IDEV,IFAIL,PERMAX,JCONN,IEND,NCON,JSENSE,
|
|
+ILADD,IRADD,MAXCON,MCON,JCON,RELPG,LNBR,IDBSIZ,NCONTS,WINDOW,
|
|
+NFOUND,IENDCT,IENDGT)
|
|
C AUTHOR: RODGER STADEN
|
|
CHARACTER SEQC2(MAXGEL),SEQG2(MAXGEL),SEQ3(MAXGEL)
|
|
CHARACTER PAD,DASH,STRAND,NAME1*6,NAME2*6
|
|
INTEGER ILADD(MAXCON),IRADD(MAXCON),RELPG(IDBSIZ),LNBR(IDBSIZ)
|
|
INTEGER CLINNO,WINDOW
|
|
EXTERNAL CLINNO
|
|
SAVE PAD,DASH
|
|
DATA PAD,DASH/',','-'/
|
|
C
|
|
C where are the overlaps?
|
|
C
|
|
C if ITYPE is 1 the overlap starts within the reading at JOINT
|
|
C else it starts at the left end of the reading at JOINT in the contig
|
|
C
|
|
C
|
|
IENDG=1
|
|
IENDC=JOINT
|
|
C ONLY LOOK AT OVERLAP WHICH IS FROM JOINT FOR LEFT TYPE JOIN
|
|
IF(ITYPE.EQ.1)THEN
|
|
IENDG=JOINT
|
|
IENDC=1
|
|
END IF
|
|
C LENGTH OF OVERLAP?
|
|
LG=IDIM2-IENDG+1
|
|
LO=MIN(IDOUT,LG)
|
|
C SAVE RAW DATA
|
|
CALL SQCOPY(SEQG2,SEQ3,IDIM2)
|
|
X=REAL(LO)
|
|
Y=X
|
|
K=IENDG+LO-1
|
|
C POINT TO CONSENSUS
|
|
J=0
|
|
C CHECK FOR OVERFLOW
|
|
IF(K.GT.MAXGEL)THEN
|
|
CALL ERROM(KBOUT,'Matching region too long in JALIGN')
|
|
IFAIL=1
|
|
RETURN
|
|
END IF
|
|
DO 200 I=IENDG,K
|
|
J=J+1
|
|
IF(SEQC2(J).EQ.SEQ3(I))GO TO 200
|
|
C IF(SEQ3(I).EQ.DASH)GO TO 200
|
|
C IF(SEQC2(J).EQ.DASH)GO TO 200
|
|
C IF(SEQC2(J).EQ.PAD)GO TO 200
|
|
X=X-1.
|
|
200 CONTINUE
|
|
X=(Y-X)*100./Y
|
|
IF(X.LT.PERMAX) THEN
|
|
IC = CLINNO(LNBR,IDBSIZ,NCONTS,JCONN)
|
|
IF(IC.EQ.0) THEN
|
|
CALL ERROM(KBOUT,'Error: contig line not found')
|
|
IFAIL = 99
|
|
RETURN
|
|
END IF
|
|
IF(JSENSE.EQ.-1) THEN
|
|
STRAND = '-'
|
|
IF(IEND.EQ.1) THEN
|
|
C
|
|
C probe is complement of right hand end of contig. Give posns assuming
|
|
C this contig is complemented.
|
|
C
|
|
IENDGT = IENDG - IRADD(JCON)
|
|
C
|
|
C Next line gives posns relative to original orientation
|
|
C IENDGT = RELPG(IC) + IRADD(JCON) - IENDG + 1
|
|
ELSE
|
|
C
|
|
C probe is complement of left end of contig. Give posns assuming this
|
|
C contig is going to be complemented.
|
|
C
|
|
IENDGT = RELPG(IC) - WINDOW + ILADD(JCON) + IENDG
|
|
END IF
|
|
ELSE
|
|
STRAND = '+'
|
|
IF(IEND.EQ.2) THEN
|
|
C
|
|
C probe is left hand end of contig in original sense
|
|
C
|
|
IENDGT = IENDG - ILADD(JCON)
|
|
ELSE
|
|
C
|
|
C probe is right hand end of contig in original sense
|
|
C
|
|
IENDGT = RELPG(IC) + IRADD(JCON) - WINDOW + IENDG
|
|
END IF
|
|
END IF
|
|
IENDCT = IENDC - ILADD(MCON)
|
|
WRITE(IDEV,*)
|
|
+ ' Possible join between contig ',JCONN,' in the ',
|
|
+ STRAND,' sense and contig ',NCON
|
|
WRITE(IDEV,1000)X
|
|
1000 FORMAT(' Percentage mismatch after alignment = ',F4.1)
|
|
WRITE(NAME1,1002)JCONN
|
|
WRITE(NAME2,1002)NCON
|
|
1002 FORMAT(I6)
|
|
CALL FMT4LP(SEQC2(1),SEQG2(IENDG),LO,IENDCT,IENDGT,IDEV,
|
|
+ NAME2,NAME1)
|
|
NFOUND = NFOUND + 1
|
|
END IF
|
|
IFAIL=0
|
|
END
|
|
SUBROUTINE ADISM5(IDIM,IDIMG,SAVPS,SAVPG,IDSAV,
|
|
+CENDS,NENDS,IDCEND,MAXCON,ILEFTS,ILC,IPOSC,IPOSG,ISENSE,
|
|
+LLINO,IMATC,ISTRAN,MAXC,KBOUT,MCON)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER ILEFTS(MAXC),ILC(MAXC),IPOSC(MAXC),IPOSG(MAXC)
|
|
INTEGER ISENSE(MAXC),LLINO(MAXC),MCON(MAXC)
|
|
INTEGER CENDS(MAXCON)
|
|
INTEGER NENDS(MAXCON)
|
|
INTEGER SAVPS(IDSAV),SAVPG(IDSAV)
|
|
NEXTC=IDIM+1
|
|
CALL BUB2AS(SAVPS,SAVPG,IDSAV)
|
|
IMATC=IMATC+1
|
|
CALL ADISM6(SAVPS(1),SAVPG(1),CENDS,NENDS,IDCEND,MAXCON,
|
|
+ ILEFTS,ILC,IPOSC,IPOSG,ISENSE,LLINO,IMATC,ISTRAN,NEXTC,MAXC,
|
|
+ KBOUT,MCON)
|
|
LEND=IDIMG-SAVPG(1)+SAVPS(1)
|
|
DO 10 I=2,IDSAV
|
|
IF((SAVPS(I).LT.LEND).AND.(SAVPS(I).LT.NEXTC))GO TO 10
|
|
IMATC=IMATC+1
|
|
CALL ADISM6(SAVPS(I),SAVPG(I),CENDS,NENDS,IDCEND,MAXCON,
|
|
+ ILEFTS,ILC,IPOSC,IPOSG,ISENSE,LLINO,IMATC,ISTRAN,NEXTC,MAXC,
|
|
+ KBOUT,MCON)
|
|
LEND=IDIMG-SAVPG(I)+SAVPS(I)
|
|
10 CONTINUE
|
|
IMATC = MIN(IMATC,MAXC)
|
|
END
|
|
SUBROUTINE ADISM6(ISAVPS,SAVPG,CENDS,NENDS,
|
|
+IDCEND,MAXCON,ILEFTS,ILC,IPOSC,IPOSG,ISENSE,LLINO,IMATC,ISTRAN,
|
|
+NEXTC,MAXC,KBOUT,MCON)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER ILEFTS(MAXC),ILC(MAXC),IPOSC(MAXC),IPOSG(MAXC)
|
|
INTEGER ISENSE(MAXC),LLINO(MAXC),MCON(MAXC)
|
|
INTEGER SAVPS,SAVPG,CENDS(MAXCON)
|
|
INTEGER NENDS(MAXCON)
|
|
SAVPS=ISAVPS-19
|
|
JJ=1
|
|
DO 5 J=2,IDCEND
|
|
IF(SAVPS.GT.CENDS(J))GO TO 5
|
|
JJ=J-1
|
|
GO TO 6
|
|
5 CONTINUE
|
|
JJ=IDCEND
|
|
6 CONTINUE
|
|
C new bit need to save contig number for alinej
|
|
SAVPS=SAVPS-1
|
|
LCL=SAVPS-CENDS(JJ)
|
|
LCR=CENDS(JJ+1)-ISAVPS-1
|
|
NEXTC=CENDS(JJ+1)+20
|
|
IF(IMATC.LE.MAXC) THEN
|
|
MCON(IMATC) = JJ
|
|
ILEFTS(IMATC)=CENDS(JJ)+20
|
|
ILC(IMATC)=LCL+LCR+1
|
|
IPOSC(IMATC)=LCL+1
|
|
IPOSG(IMATC)=SAVPG
|
|
LLINO(IMATC)=NENDS(JJ)
|
|
ISENSE(IMATC)=1
|
|
IF(ISTRAN.EQ.2)ISENSE(IMATC)=-1
|
|
ELSE
|
|
CALL ERROM(KBOUT,'Warning: too many overlaps')
|
|
END IF
|
|
END
|
|
C JCONS
|
|
SUBROUTINE JCONS(RELPG,LNGTHG,LNBR,RNBR,NAMPRO,NGELS,NCONTS,
|
|
+SEQ1,IDIM1,GEL,IDBSIZ,ISTART,JOB,LLINO,LINCON,LREG,RREG,TEMP,
|
|
+CHRSIZ,MAXGL2,KBOUT,
|
|
+IDEVW,IFAIL,MAXGEL,IDM,PERCD,TGEL,ILADD,IRADD,MAXCON,JCON,
|
|
+IWIN,NBAD)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),CHRSIZ,ILADD(MAXCON),IRADD(MAXCON)
|
|
INTEGER LREG,RREG,X,Y,TEMP(CHRSIZ,MAXGL2)
|
|
CHARACTER SEQ1(IDIM1)
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
CHARACTER GEL(MAXGEL),TGEL(MAXGEL)
|
|
CHARACTER NAMPRO*(*)
|
|
CALL BUSY(KBOUT)
|
|
C WRITE(*,*)'NGELS,NCONTS,IDIM1,ISTART,LREG,RREG'
|
|
C WRITE(*,*)NGELS,NCONTS,IDIM1,ISTART,LREG,RREG
|
|
C JOB = 0 do all contigs, = 1 do one contig llino,lincon
|
|
IFAIL=0
|
|
IF (JOB.EQ.0) THEN
|
|
N=IDBSIZ-NCONTS
|
|
NCONS = 0
|
|
DO 110 I=N,IDBSIZ-1
|
|
J=LNBR(I)
|
|
X=1
|
|
Y=RELPG(I)
|
|
ISTART=ISTART+1
|
|
IF((ISTART+19+Y+2*MAXGEL).GT.IDIM1)THEN
|
|
WRITE(KBOUT,1009)IDIM1
|
|
1009 FORMAT(
|
|
+ ' Maximum consensus length(',I6,') exceeded',/,
|
|
+ ' calculation aborted')
|
|
IFAIL=1
|
|
RETURN
|
|
END IF
|
|
CALL ADDTIT(SEQ1(ISTART),NAMPRO,J,ISTART)
|
|
NCONS = NCONS + 1
|
|
IDIN = 0
|
|
IF(IWIN.GT.0) THEN
|
|
IDIN = MAXGEL
|
|
CALL GETEX(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,I,1,
|
|
+ GEL,TGEL,IDIN,IWIN,NBAD)
|
|
IF(IDIN.GT.0)CALL SQCOPY(GEL,SEQ1(ISTART),IDIN)
|
|
END IF
|
|
ILADD(NCONS) = IDIN
|
|
ISTART = ISTART + IDIN
|
|
CALL SUMMER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ SEQ1(ISTART),Y,GEL,X,Y,J,IDBSIZ,TEMP,CHRSIZ,MAXGL2,IDEVW,
|
|
+ MAXGEL,IDM,PERCD)
|
|
ISTART=ISTART+Y-1
|
|
IDIN = 0
|
|
IF(IWIN.GT.0) THEN
|
|
IDIN = MAXGEL
|
|
CALL GETEX(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,I,2,
|
|
+ GEL,TGEL,IDIN,IWIN,NBAD)
|
|
IF(IDIN.GT.0)CALL SQCOPY(GEL,SEQ1(ISTART+1),IDIN)
|
|
END IF
|
|
IRADD(NCONS) = IDIN
|
|
ISTART = ISTART + IDIN
|
|
110 CONTINUE
|
|
ELSE
|
|
C
|
|
C do it for one contig
|
|
C
|
|
J = LLINO
|
|
I = LINCON
|
|
X = LREG
|
|
Y = RREG
|
|
ISTART = ISTART + 1
|
|
C WRITE(*,*)
|
|
C +'LLINO,LINCON,NCONTS,ISTART',LLINO,LINCON,NCONTS,ISTART
|
|
IF((ISTART+19+Y+2*MAXGEL).GT.IDIM1)THEN
|
|
WRITE(KBOUT,1009)IDIM1
|
|
IFAIL=1
|
|
RETURN
|
|
END IF
|
|
CALL ADDTIT(SEQ1(ISTART),NAMPRO,J,ISTART)
|
|
NCONS = NCONTS
|
|
IDIN = 0
|
|
IF(IWIN.GT.0) THEN
|
|
IDIN = MAXGEL
|
|
CALL GETEX(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,I,1,
|
|
+ GEL,TGEL,IDIN,IWIN,NBAD)
|
|
IF(IDIN.GT.0)CALL SQCOPY(GEL,SEQ1(ISTART),IDIN)
|
|
END IF
|
|
C WRITE(*,*)'ILADD(JCON)',IDIN
|
|
ILADD(JCON) = IDIN
|
|
ISTART = ISTART + IDIN
|
|
CALL SUMMER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ SEQ1(ISTART),Y,GEL,X,Y,J,IDBSIZ,TEMP,CHRSIZ,MAXGL2,
|
|
+ IDEVW,MAXGEL,IDM,PERCD)
|
|
C repleced by next line 26-1-93 ISTART = ISTART + Y - 1
|
|
ISTART = ISTART + RREG - LREG
|
|
IDIN = 0
|
|
IF(IWIN.GT.0) THEN
|
|
IDIN = MAXGEL
|
|
CALL GETEX(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,I,2,
|
|
+ GEL,TGEL,IDIN,IWIN,NBAD)
|
|
IF(IDIN.GT.0)CALL SQCOPY(GEL,SEQ1(ISTART+1),IDIN)
|
|
END IF
|
|
C WRITE(*,*)'IRADD(JCON)',IDIN
|
|
IRADD(JCON) = IDIN
|
|
ISTART = ISTART + IDIN
|
|
END IF
|
|
END
|
|
SUBROUTINE GETEX(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,ICONT,IEND,
|
|
+GEL,GELT,ID,IWIN,NBAD)
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ),RELPG(IDBSIZ)
|
|
CHARACTER GEL(ID),GELT(ID)
|
|
EXTERNAL NOK2
|
|
MAXGEL = ID
|
|
IDT = 0
|
|
C routine to find a possible extension to a contig by looking in a tag file
|
|
C contig ICONT end IEND = 1 =left 2=right
|
|
C return data in GEL, of length ID
|
|
C The worst aspect of this is that if we dont choose the very end reading
|
|
C we dont know where it lies relative to the consensus. For now just assume
|
|
C there are no length difference between the extension and the consensus
|
|
C and just add it on the end
|
|
IF(IEND.EQ.1) THEN
|
|
LMOST = 1
|
|
IGEL = LNBR(ICONT)
|
|
10 CONTINUE
|
|
IF(IGEL.EQ.0) GO TO 100
|
|
IF(RELPG(IGEL).GT.MAXGEL) GO TO 100
|
|
IF(LNGTHG(IGEL).LT.0) THEN
|
|
ID = MAXGEL
|
|
CALL GETEXT(IGEL,GELT,ID,IOK)
|
|
IF(IOK.EQ.0) THEN
|
|
C look for number of n's here and change id accordingly
|
|
K = NOK2(GELT,ID,IWIN,NBAD)
|
|
LT = MIN(LMOST,RELPG(IGEL)-K)
|
|
IF(LT.LT.LMOST) THEN
|
|
IS = RELPG(IGEL)
|
|
N = K - IS + 1
|
|
CALL SQCOPY(GELT(IS),GEL,N)
|
|
CALL SQREV(GEL,N)
|
|
CALL SQCOM(GEL,N)
|
|
IDT = N
|
|
LMOST = LT
|
|
END IF
|
|
ELSE
|
|
C WRITE(*,*)'COCKUP IN GETEXT, gel',IGEL
|
|
END IF
|
|
END IF
|
|
IGEL = RNBR(IGEL)
|
|
GO TO 10
|
|
ELSE
|
|
IGEL = RNBR(ICONT)
|
|
LMOST = RELPG(ICONT)
|
|
IDC = RELPG(ICONT)
|
|
20 CONTINUE
|
|
IF(IGEL.EQ.0) GO TO 100
|
|
IF(LMOST-RELPG(IGEL).GT.MAXGEL) GO TO 100
|
|
IF(LNGTHG(IGEL).GT.0) THEN
|
|
ID = MAXGEL
|
|
CALL GETEXT(IGEL,GELT,ID,IOK)
|
|
IF(IOK.EQ.0) THEN
|
|
K = NOK2(GELT,ID,IWIN,NBAD)
|
|
LT = MAX(LMOST,RELPG(IGEL)+LNGTHG(IGEL)+K-1)
|
|
IF(LT.GT.LMOST) THEN
|
|
IS = RELPG(ICONT) - (RELPG(IGEL) + LNGTHG(IGEL)) + 2
|
|
N = K - IS + 1
|
|
CALL SQCOPY(GELT(IS),GEL,N)
|
|
IDT = N
|
|
LMOST = LT
|
|
END IF
|
|
ELSE
|
|
C WRITE(*,*)'COCKUP IN GETEXT, GEL',IGEL
|
|
END IF
|
|
END IF
|
|
IGEL = LNBR(IGEL)
|
|
GO TO 20
|
|
END IF
|
|
100 CONTINUE
|
|
ID = IDT
|
|
END
|
|
INTEGER FUNCTION NOK2(GEL,ID,IWIN,NBADIN)
|
|
CHARACTER GEL(ID)
|
|
PARAMETER (MAXPOS = 101)
|
|
INTEGER POSNS(MAXPOS),R
|
|
EXTERNAL KWRAP
|
|
C count N's over a window of iwin, return position
|
|
C when over NBAD
|
|
C INIT
|
|
NBAD = NBADIN + 1
|
|
I = 0
|
|
N = 0
|
|
NOK2 = ID
|
|
L = 1
|
|
R = 0
|
|
IF(NBAD.GT.MAXPOS)THEN
|
|
WRITE(*,*)'Scream: nok2 not happy'
|
|
RETURN
|
|
END IF
|
|
10 CONTINUE
|
|
I = I + 1
|
|
IF(I.GT.ID) RETURN
|
|
IF(GEL(I).EQ.'-') THEN
|
|
N = N + 1
|
|
R = KWRAP(R,NBAD)
|
|
POSNS(R) = I
|
|
IF(N.GE.NBAD) THEN
|
|
IF(POSNS(R)-POSNS(L)+1.LT.IWIN) THEN
|
|
NOK2 = POSNS(L)
|
|
RETURN
|
|
END IF
|
|
L = KWRAP(L,NBAD)
|
|
END IF
|
|
END IF
|
|
GO TO 10
|
|
END
|
|
INTEGER FUNCTION KWRAP(I,J)
|
|
IT = I + 1
|
|
IF(IT.GT.J) IT = 1
|
|
KWRAP = IT
|
|
END
|
|
SUBROUTINE PADSHF(RELPG,LNGTHG,LNBR,RNBR,
|
|
+READS,IDIM1,LREG,RREG,LLINO,IDBSIZ,FREQS,CHRSIZ,MAXGL2,
|
|
+IDEVW,MAXGEL,READPS,READNS,PSTART,NPADS,MXG,PADPOS,MAXP,
|
|
+KBOUT)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),CHRSIZ
|
|
INTEGER LREG,RREG
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
CHARACTER READS(IDIM1),PAD,DASH
|
|
C PARAMETER (MXG = 1000)
|
|
C PARAMETER (MAXP = MXG * 20)
|
|
INTEGER FREQS(CHRSIZ,MAXGL2),READPS(MXG),READNS(MXG)
|
|
INTEGER CR,SHIFTF,READST,PCF,WINS,WINE,PG,PC,PCR,CHNRP1
|
|
INTEGER PADPOS(MAXP),PSTART(MXG),NPADS(MXG),PADST,PGE
|
|
LOGICAL FIRSTP,PASS1
|
|
EXTERNAL INDEXS,CHNRP1,NSAMPL
|
|
SAVE PAD,DASH
|
|
DATA PAD/'*'/,DASH/'-'/
|
|
C
|
|
C region to process is lreg to rreg
|
|
C current active region is wins to wine
|
|
C current reading is cr
|
|
C we process all the readings covering wins to wine,
|
|
C write them to disk, then move wins and wine, and get the first
|
|
C reading that covers them, etc
|
|
C we save readings in reads and their numbers in readns
|
|
C and their start positions in reads in readps
|
|
C we accumulate frequencies in freqs
|
|
C
|
|
NEDITS = 0
|
|
NSWAPT = 0
|
|
INDPAD = INDEXS(PAD,JUNK)
|
|
WINS = LREG
|
|
WINE = MIN(WINS + MAXGL2 - 1,RREG)
|
|
1 CONTINUE
|
|
CR = CHNRP1(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,
|
|
+ LLINO,WINS)
|
|
IF (CR.EQ.0) RETURN
|
|
DO 10 I = 1,MAXGL2
|
|
DO 10 J = 1,CHRSIZ
|
|
FREQS(J,I) = 0
|
|
10 CONTINUE
|
|
NREADS = 0
|
|
READST = 1
|
|
PADST = 1
|
|
PASS1 = .TRUE.
|
|
20 CONTINUE
|
|
CALL READW(IDEVW,CR,READS(READST),MAXGEL)
|
|
IF (NREADS.EQ.MXG) THEN
|
|
CALL ERROM(KBOUT,'Too many reads for buffer')
|
|
RETURN
|
|
END IF
|
|
NREADS = NREADS + 1
|
|
READNS(NREADS) = CR
|
|
READPS(NREADS) = READST
|
|
READST = READST + ABS(LNGTHG(CR))
|
|
SHIFTF = RELPG(CR) - WINS
|
|
DO 30 I = MAX(1,WINS-RELPG(CR) + 1),
|
|
+ MIN(ABS(LNGTHG(CR)),WINE-RELPG(CR)+1)
|
|
JJ = INDEXS(READS(READPS(NREADS)+I-1),JSCORE)
|
|
JJJ = SHIFTF + I
|
|
FREQS(JJ,JJJ) = FREQS(JJ,JJJ) + 1
|
|
30 CONTINUE
|
|
IF (RNBR(CR).NE.0) THEN
|
|
CR = RNBR(CR)
|
|
IF (RELPG(CR).LE.WINE) GO TO 20
|
|
END IF
|
|
DO 35 I=1,NREADS
|
|
NPADS(I) = 0
|
|
35 CONTINUE
|
|
40 CONTINUE
|
|
NSWAP = 0
|
|
C pg is position in reading ir
|
|
C PCR is position in reads
|
|
C pc is position in contig
|
|
C pcf is position in freqs
|
|
C
|
|
C we swap if chars to left and right are not already pads and
|
|
C theres at least one of the char to the left in the current column
|
|
C or the char to the left is a dash (not always beneficial) and
|
|
C the pads are not already aligned with only 1 other symbol
|
|
C
|
|
C we put the element number of the first pad for each reading at padst
|
|
C we count the number of pads in npads
|
|
C we store the pad positions in padpos
|
|
C
|
|
DO 60 IR = 1,NREADS
|
|
SHIFTF = WINS - 1
|
|
FIRSTP = .TRUE.
|
|
DO 50 PG = MAX(2,WINS-RELPG(READNS(IR))+2),
|
|
+ MIN(ABS(LNGTHG(READNS(IR))),WINE-RELPG(READNS(IR))+1)
|
|
PCR = PG + READPS(IR) - 1
|
|
C
|
|
C horrible exception for case where base 1 in read is pad
|
|
C we wont see it this pass but will when we write out tags
|
|
C so make sure we do see it now or we cockup later
|
|
C
|
|
IF (PASS1) THEN
|
|
IF (PG.EQ.2) THEN
|
|
IF (READS(PCR-1).EQ.PAD) THEN
|
|
PSTART(IR) = PADST
|
|
NPADS(IR) = NPADS(IR) + 1
|
|
PADPOS(PADST) = PG - 1
|
|
IF (PADST.EQ.MAXP) THEN
|
|
CALL ERROM(KBOUT,'Too many pads for buffer')
|
|
RETURN
|
|
END IF
|
|
PADST = PADST + 1
|
|
FIRSTP = .FALSE.
|
|
END IF
|
|
END IF
|
|
END IF
|
|
IF (READS(PCR).EQ.PAD) THEN
|
|
IF (PASS1) THEN
|
|
IF (FIRSTP) PSTART(IR) = PADST
|
|
NPADS(IR) = NPADS(IR) + 1
|
|
PADPOS(PADST) = PG
|
|
IF (PADST.EQ.MAXP) THEN
|
|
CALL ERROM(KBOUT,'Too many pads for buffer')
|
|
RETURN
|
|
END IF
|
|
PADST = PADST + 1
|
|
FIRSTP = .FALSE.
|
|
END IF
|
|
PC = PG + RELPG(READNS(IR)) - 1
|
|
INPCRM = INDEXS(READS(PCR-1),JSCORM)
|
|
IF (INPCRM.NE.INDPAD) THEN
|
|
PCF = PC - SHIFTF
|
|
IF (((FREQS(INPCRM,PCF).NE.0).OR.
|
|
+ (READS(PCR-1).EQ.DASH)).AND.
|
|
+ (FREQS(INDPAD,PCF).LT.(NSAMPL(FREQS(1,PCF),CHRSIZ)-1)))
|
|
+ THEN
|
|
C
|
|
C swap
|
|
C
|
|
READS(PCR) = READS(PCR-1)
|
|
READS(PCR-1) = PAD
|
|
INPCR = INDEXS(READS(PCR),JSCORE)
|
|
FREQS(INPCRM,PCF-1) = FREQS(INPCRM,PCF-1) - 1
|
|
FREQS(INPCR,PCF) = FREQS(INPCR,PCF) - 1
|
|
FREQS(INPCR,PCF-1) = FREQS(INPCR,PCF-1) + 1
|
|
FREQS(INPCRM,PCF) = FREQS(INPCRM,PCF) + 1
|
|
NSWAP = NSWAP + 1
|
|
END IF
|
|
END IF
|
|
END IF
|
|
50 CONTINUE
|
|
60 CONTINUE
|
|
C WRITE(*,*)'NUMBER OF SWAPS',NSWAP
|
|
NSWAPT = NSWAPT + NSWAP
|
|
PASS1 = .FALSE.
|
|
IF (NSWAP.GT.0) GO TO 40
|
|
C
|
|
C no more movement so write back wot weve changed so far
|
|
C
|
|
DO 70 IR = 1,NREADS
|
|
K = 0
|
|
NPAD = 0
|
|
DO 65 PG = MAX(1,WINS-RELPG(READNS(IR))+1),
|
|
+ MIN(ABS(LNGTHG(READNS(IR))),WINE-RELPG(READNS(IR))+1)
|
|
PCR = PG + READPS(IR) - 1
|
|
IF (READS(PCR).EQ.PAD) THEN
|
|
L = PADPOS(PSTART(IR)+NPAD)
|
|
IF (PG.LT.L) THEN
|
|
PC = PG + RELPG(READNS(IR)) - 1
|
|
LE = L
|
|
PGE = PG
|
|
IF (LNGTHG(READNS(IR)).LT.0) THEN
|
|
LE = 1 - LNGTHG(READNS(IR)) - LE
|
|
PGE = 1 - LNGTHG(READNS(IR)) - PGE
|
|
C WRITE(*,*)READNS(IR),' insert at',PGE,PC,' del at',LE
|
|
CALL DELEDT(READNS(IR),LE,'*')
|
|
CALL INSEDT(READNS(IR),PGE,'*')
|
|
ELSE
|
|
C WRITE(*,*)READNS(IR),' del at',LE,' insert at',PGE,PC
|
|
CALL DELEDT(READNS(IR),LE,'*')
|
|
CALL INSEDT(READNS(IR),PGE,'*')
|
|
END IF
|
|
NEDITS = NEDITS + 1
|
|
END IF
|
|
NPAD = NPAD + 1
|
|
END IF
|
|
65 CONTINUE
|
|
CALL WRITEW(IDEVW,READNS(IR),READS(READPS(IR)),MAXGEL)
|
|
70 CONTINUE
|
|
C
|
|
C we allow the next chunk to overlap the previous one by 20
|
|
C in case some pads at the right edge of it could still move
|
|
C further left. Alternatively we could look for 2 adjacent
|
|
C positions which have 100% consensus and are different
|
|
C
|
|
IF (WINE.NE.RREG) THEN
|
|
WINS = MAX(LREG,WINE-19)
|
|
WINE = MIN(WINS + MAXGL2 - 1,RREG)
|
|
GO TO 1
|
|
END IF
|
|
C WRITE(*,*)'Total swaps',NSWAPT
|
|
WRITE(KBOUT,*)'Number of pads moved',NEDITS
|
|
END
|
|
INTEGER FUNCTION NSAMPL(COUNTS,IDM)
|
|
C
|
|
C count the number of reads covering this position
|
|
C
|
|
INTEGER IDM
|
|
INTEGER COUNTS(IDM)
|
|
ISUM = 0
|
|
DO 5 I=1,IDM
|
|
ISUM = ISUM + COUNTS(I)
|
|
5 CONTINUE
|
|
NSAMPL = ISUM
|
|
END
|
|
INTEGER FUNCTION CONOK(
|
|
+RELPG,LNGTHG,LNBR,RNBR,NAMPRO,NGELS,NCONTS,
|
|
+SEQ,MAXSEQ,SEQ2,IDBSIZ,TEMP3,
|
|
+ECHRSZ,MAXGL2,KBOUT,IDEVW,IDEV,
|
|
+MAXGEL,IDM,PERCD,CENDS,NENDS,MAXCON)
|
|
INTEGER RELPG(IDBSIZ),LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
INTEGER ECHRSZ,TEMP3(ECHRSZ,MAXGL2),CENDS(MAXCON),NENDS(MAXCON)
|
|
CHARACTER SEQ(MAXSEQ),SEQ2(MAXGEL),NAMPRO*(*)
|
|
EXTERNAL NOK2
|
|
C
|
|
C routine to calc a consensus and look for ibad dashes in iwin bases
|
|
C
|
|
C returns number of cocked-up contigs or -1 for error in consensus calc
|
|
C for each contig it stops where it find the first problem.
|
|
C
|
|
CONOK = 0
|
|
JOB = 0
|
|
IDIM1 = 0
|
|
CALL ACONSN(RELPG,LNGTHG,LNBR,RNBR,NAMPRO,NGELS,NCONTS,
|
|
+SEQ,MAXSEQ,SEQ2,IDBSIZ,IDIM1,JOB,KDUMM,KDUMM,KDUMM,TEMP3,
|
|
+ECHRSZ,MAXGL2,KBOUT,IDEVW,IFAIL,MAXGEL,IDM,PERCD)
|
|
IF(IFAIL.NE.0) THEN
|
|
CALL ERROM(KBOUT,'Error calculating consensus')
|
|
CONOK = -1
|
|
RETURN
|
|
END IF
|
|
CALL FNDCON(SEQ,IDIM1,CENDS,NENDS,IDCEND,MAXCON,KBOUT)
|
|
IWIN = 20
|
|
NBAD = 15
|
|
DO 10 I=1,IDCEND
|
|
J = CENDS(I) + 20
|
|
IDIM = CENDS(I+1) - J
|
|
K = NOK2(SEQ(J),IDIM,IWIN,NBAD)
|
|
IF (K.LT.IDIM) THEN
|
|
WRITE(IDEV,*)'Problem at position',K,' In contig',NENDS(I)
|
|
CONOK = CONOK + 1
|
|
END IF
|
|
10 CONTINUE
|
|
IF (CONOK.EQ.0) THEN
|
|
WRITE(KBOUT,*)
|
|
+' Consensus has no segments with',NBAD,' dashes in',IWIN
|
|
ELSE
|
|
WRITE(KBOUT,*)
|
|
+ CONOK,' contigs have segments with',NBAD,' dashes in',IWIN
|
|
END IF
|
|
END
|
|
SUBROUTINE AUTOM(RELPG,LNGTHG,LNBR,RNBR,MAXDB,IDBSIZ,
|
|
+NGELS,NCONTS,LLINO,LINCON,MAXGEL,
|
|
+TEMP3,SEQ1,MAXSEQ,SEQ2,SEQ5,SEQC2,SEQG2,
|
|
+MAXGLM,MAXGL2,ECHRSZ,
|
|
+SAV1,SAV2,SAV3,MAXSAV,
|
|
+KBIN,KBOUT,IDEV1,IDEV2,IDEV3,IDEV8,IDEV,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,NAMARC,NAMPRO,FILE,
|
|
+PERCD,IOPEN,IDM,SEQG3,SEQC3,IOK,IDEVC,IDEVT)
|
|
INTEGER ECHRSZ,SELCON
|
|
INTEGER RELPG(MAXDB)
|
|
INTEGER LNGTHG(MAXDB),LNBR(MAXDB),RNBR(MAXDB)
|
|
INTEGER TEMP3(ECHRSZ,MAXGL2)
|
|
CHARACTER SEQC2(MAXGLM,2),SEQG2(MAXGLM,2)
|
|
CHARACTER SEQ1(MAXSEQ),SEQ2(MAXGLM)
|
|
INTEGER SAV1(MAXSAV),SAV2(MAXSAV),SAV3(MAXSAV)
|
|
CHARACTER NAMARC*(*),NAMPRO*(*),FILE*(*)
|
|
CHARACTER SEQ5(MAXGLM),HELPF*(*),SEQG3(MAXGLM),SEQC3(MAXGLM)
|
|
C
|
|
C take the poor data from reads to see if it aligns with the consensus
|
|
C
|
|
CALL DBCHEK(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,
|
|
+ TEMP3,IERR,KBOUT)
|
|
IF(IERR.GT.1) RETURN
|
|
IFAIL = 0
|
|
IF(NGELS.LT.1) RETURN
|
|
C MN = 0
|
|
C MX = 50
|
|
C MAXPG = 20
|
|
C CALL GETINT(MN,MX,MAXPG,
|
|
C +'Maximum pads per sequence',
|
|
C +IVAL,KBIN,KBOUT,
|
|
C +IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
C IF(IOK.NE.0) RETURN
|
|
MAXPG = IVAL
|
|
MAXPC = IVAL
|
|
RMN = 0.
|
|
RMX = 100.
|
|
PERMAX = 13.
|
|
CALL GETRL(RMN,RMX,PERMAX,
|
|
+ 'Maximum percent mismatch after alignment',
|
|
+ VAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
PERMAX = VAL
|
|
MN = 1
|
|
MX = MAXGEL
|
|
IWING = 100
|
|
CALL GETINT(MN,MX,IWING,
|
|
+ 'Window size for good data scan',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
IWING = IVAL
|
|
MN = 1
|
|
MX = MIN(100,IWING)
|
|
NBAD = MIN(IWING,5)
|
|
CALL GETINT(MN,MX,NBAD,
|
|
+ 'Maximum number of dashes in scan window',
|
|
+ IVAL,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) RETURN
|
|
NBAD = IVAL
|
|
CALL YESNO(IREPN,'Save failed names in a file',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF (IREPN.LT.0) RETURN
|
|
IF (IREPN.EQ.0) THEN
|
|
CALL OPENF1(IDEV8,FILE,1,IOK,KBIN,KBOUT,
|
|
+ 'File for names of failed reads',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH)
|
|
IF (IOK.NE.0) RETURN
|
|
END IF
|
|
CALL YESNO(SELCON,'Select contigs',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(SELCON.LT.0) GO TO 100
|
|
IF(SELCON.EQ.0) THEN
|
|
CALL GETLN3(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,LINCON,
|
|
+ LLINO,NULGEL,IOK,IDBSIZ,KBIN,KBOUT,IDEV3,
|
|
+ 'Contig identifier',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH)
|
|
IF(IOK.NE.0) GO TO 100
|
|
CALL GETREG(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ 1,RELPG(LINCON),LREG,RREG,LINCON,LLINO,IDBSIZ,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,IOK)
|
|
IF(IOK.NE.0) GO TO 100
|
|
END IF
|
|
IDIM2=MAXGEL
|
|
CALL AUTOMN(SEQ1,MAXSEQ,SEQ2,IDIM2,
|
|
+LLINO,LINCON,IFCOMP,LREG,RREG,
|
|
+IDEV,MAXGEL,MAXGLM,
|
|
+SAV1,SAV2,SAV3,MAXSAV,
|
|
+SEQG2,SEQC2,
|
|
+PERMAX,SEQG3,SEQC3,
|
|
+RELPG,LNBR,IDBSIZ,NGELS,NCONTS,
|
|
+LNGTHG,RNBR,IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,PERCD,IWING,NBAD,
|
|
+ECHRSZ,TEMP3,MAXGL2,IDM,NAMPRO,KBIN,KBOUT,SELCON,IREPN,IDEV8,
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
100 CONTINUE
|
|
IF(IREPN.EQ.0) CLOSE(UNIT=IDEV8)
|
|
END
|
|
SUBROUTINE AUTOMN(SEQ1,MAXSEQ,GEL,IDIMGI,
|
|
+LLINO,LINCON,IFCOMP,LREG,RREG,
|
|
+IDEV,MAXGEL,MAXGLM,
|
|
+SAVPS,SAVPG,SAVL,MAXSAV,
|
|
+SEQG2,SEQC2,
|
|
+PERMAX,SEQG3,SEQC3,
|
|
+RELPG,LNBR,IDBSIZ,NGELS,NCONTS,
|
|
+LNGTHG,RNBR,IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,PERCD,IWING,NBAD,
|
|
+ECHRSZ,TEMP3,MAXGL2,IDM,NAMPRO,KBIN,KBOUT,SELCON,IREPN,IDEV8,
|
|
+IHELPS,IHELPE,HELPF,IDEVH)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER ECHRSZ
|
|
INTEGER SAVPS(MAXSAV),SAVPG(MAXSAV),SAVL(MAXSAV)
|
|
CHARACTER NAMPRO*(*),HELPF*(*)
|
|
CHARACTER SEQ1(MAXSEQ),GEL(MAXGLM)
|
|
CHARACTER SEQG2(MAXGLM,2),SEQC2(MAXGLM,2)
|
|
CHARACTER SEQG3(MAXGLM),SEQC3(MAXGLM)
|
|
CHARACTER NAMT*16
|
|
PARAMETER (MAXCON = 1,MING = 10,MAXG = 500)
|
|
INTEGER ILADD(MAXCON),IRADD(MAXCON)
|
|
INTEGER RREG,FIRSTR,SELCON,FIRSTC
|
|
INTEGER RELPG(IDBSIZ),LNBR(IDBSIZ),LNGTHG(IDBSIZ),RNBR(IDBSIZ)
|
|
INTEGER TEMP3(ECHRSZ,MAXGL2)
|
|
INTEGER CHNRP1,XVERSN
|
|
EXTERNAL CHNRP1,XVERSN,NOK2
|
|
JOB = 1
|
|
JCON = 1
|
|
NPROBS = 0
|
|
IF (SELCON.EQ.0) THEN
|
|
FIRSTC = LINCON
|
|
LASTC = LINCON
|
|
ELSE IF (SELCON.EQ.1) THEN
|
|
FIRSTC = IDBSIZ-NCONTS
|
|
LASTC = IDBSIZ-1
|
|
ELSE
|
|
CALL ERROM(KBOUT,'Error in AUTOM: unexpected option')
|
|
RETURN
|
|
END IF
|
|
DO 50 ICONT = FIRSTC,LASTC
|
|
IF (SELCON.EQ.1) THEN
|
|
LINCON = ICONT
|
|
LLINO = LNBR(LINCON)
|
|
LREG = 1
|
|
RREG = RELPG(LINCON)
|
|
END IF
|
|
IDIM = 0
|
|
CALL JCONS(RELPG,LNGTHG,LNBR,RNBR,NAMPRO,NGELS,NCONTS,
|
|
+ SEQ1,MAXSEQ,GEL,IDBSIZ,IDIM,JOB,LLINO,LINCON,LREG,RREG,
|
|
+ TEMP3,
|
|
+ ECHRSZ,MAXGL2,KBOUT,IDEV2,IOK,MAXGEL,IDM,PERCD,SEQG3,
|
|
+ ILADD,IRADD,MAXCON,JCON,0,NBAD)
|
|
IF (IOK.NE.0) RETURN
|
|
C CALL FMTDB(SEQ1,IDIM,1,IDIM,60,KBOUT)
|
|
LENCON = RREG - LREG + 1
|
|
ISTRAN = 1
|
|
5 CONTINUE
|
|
C
|
|
C main loop
|
|
C
|
|
FIRSTR = CHNRP1(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,
|
|
+ LNBR(LINCON),LREG)
|
|
IF(FIRSTR.EQ.0) GO TO 50
|
|
IRNO = FIRSTR
|
|
10 CONTINUE
|
|
IF (LNGTHG(IRNO).GT.0) THEN
|
|
IDIN = MIN(MAXG-LNGTHG(IRNO)+1,MAXGLM)
|
|
IF(IDIN.GT.MING) THEN
|
|
CALL GETEXT(IRNO,GEL,IDIN,IOK)
|
|
IF(IOK.EQ.0) THEN
|
|
K = NOK2(GEL,IDIN,IWING,NBAD)
|
|
IDIN = K
|
|
ELSE
|
|
IDIN = 0
|
|
END IF
|
|
ELSE
|
|
IDIN = 0
|
|
END IF
|
|
C
|
|
C fudge idin (it seems to return 1 less char than its value!)
|
|
C
|
|
IDIN = IDIN - 2
|
|
C WRITE(*,*)'IDIN',IDIN
|
|
IF(IDIN.GT.MING) THEN
|
|
C
|
|
C trim gel if it goes off the end of the consensus!
|
|
C
|
|
IDIN = MIN((LENCON-(RELPG(IRNO)+LNGTHG(IRNO)-1)),IDIN)
|
|
C CALL FMTDB(GEL,IDIN,1,IDIN,60,KBOUT)
|
|
JDIM22 = IDIN
|
|
JDOUT = MAXGEL
|
|
IDSAV = MAXSAV
|
|
JPOSC = RELPG(IRNO)+LNGTHG(IRNO) - LREG + 1
|
|
LC = LENCON - JPOSC + 1
|
|
IPC = RELPG(IRNO) + LNGTHG(IRNO) - 1
|
|
C
|
|
C fudge the consensus position as we dont seem to get alignments for
|
|
C cases where a deletion is required at the start of the consensus
|
|
C
|
|
C CALL ALINEM(SEQ1(20+JPOSC),GEL,SEQG3,SEQC3,
|
|
CALL ALINEM(SEQ1(18+JPOSC),GEL,SEQG3,SEQC3,
|
|
+ SAVPS,SAVPG,SAVL,IDSAV,LC,JDIM22,JDOUT,
|
|
+ 1,1,
|
|
+ JFAIL,PERMAX,KBOUT,IDEV,
|
|
+ MAXGEL,PERMS,IPC,IRNO,LNBR(LINCON),IFAIL)
|
|
IF (IFAIL.EQ.0) THEN
|
|
IF (JFAIL.NE.0) NPROBS = NPROBS + 1
|
|
IF ((JFAIL.NE.0).AND.(IREPN.EQ.0)) THEN
|
|
CALL READN(IDEV3,IRNO,NAMT)
|
|
WRITE(IDEV8,1000)NAMT,IRNO,LLINO,IPC
|
|
1000 FORMAT(' ',A,' ',3I6)
|
|
END IF
|
|
IF((JFAIL.NE.0).AND.(XVERSN().NE.0).AND.(IDEV.EQ.KBOUT))
|
|
+ THEN
|
|
CALL BPAUSE(KBIN,KBOUT,IOK)
|
|
IF (IOK.NE.0) GO TO 100
|
|
CALL YESNO(IJOIN,'Use the editor',
|
|
+ IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
IF(IJOIN.LT.0) GO TO 100
|
|
IF (IJOIN.EQ.0) THEN
|
|
CALL CXEDIT(IDEV1,IDEV2,IDEV3,IDEVT,IDEVC,
|
|
+ RELPG,LNGTHG,LNBR,RNBR,MAXGEL,
|
|
+ IDBSIZ,LINCON,LNBR(LINCON),IRNO,LNGTHG(IRNO),
|
|
+ PERCD,IDM,1,IOK)
|
|
END IF
|
|
CALL BUSY(KBOUT)
|
|
END IF
|
|
END IF
|
|
END IF
|
|
END IF
|
|
IRNO = RNBR(IRNO)
|
|
IF(IRNO.NE.0) THEN
|
|
IF((RELPG(IRNO)+LNGTHG(IRNO)).LT.RREG) GO TO 10
|
|
END IF
|
|
IF (ISTRAN.EQ.1) THEN
|
|
CALL CMPLMT(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,LINCON,
|
|
+ LNBR(LINCON),GEL,IDBSIZ,KBOUT,IDEV1,IDEV2,MAXGEL)
|
|
LREG = RELPG(LINCON) - RREG + 1
|
|
RREG = LREG + LENCON - 1
|
|
CALL SQREV(SEQ1(21),LENCON)
|
|
CALL SQCOM(SEQ1(21),LENCON)
|
|
ISTRAN = 2
|
|
GO TO 5
|
|
END IF
|
|
CALL CMPLMT(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,LINCON,
|
|
+ LNBR(LINCON),GEL,IDBSIZ,KBOUT,IDEV1,IDEV2,MAXGEL)
|
|
ISTRAN = 3
|
|
50 CONTINUE
|
|
100 CONTINUE
|
|
WRITE(KBOUT,*)'Number of possible problems ',NPROBS
|
|
C
|
|
C if required return the contig to original sense
|
|
C
|
|
IF (ISTRAN.EQ.2) THEN
|
|
CALL CMPLMT(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,LINCON,
|
|
+ LNBR(LINCON),GEL,IDBSIZ,KBOUT,IDEV1,IDEV2,MAXGEL)
|
|
END IF
|
|
END
|
|
SUBROUTINE ALINEM(SEQ1,SEQ2,SEQG2,SEQC2,ISAV1,ISAV2,ISAV3,
|
|
+IDSAV,IDC,IDIM2,IDOUT,IC1,IG1,
|
|
+IFAIL,PERMAX,KBOUT,IDEV,MAXGEL,PERCM,IPC,IRNO,LLINO,JFAIL)
|
|
C AUTHOR: RODGER STADEN
|
|
CHARACTER SEQ1(IDC),SEQ2(IDIM2),SEQG2(IDOUT),SEQC2(IDOUT)
|
|
INTEGER ISAV1(IDSAV),ISAV2(IDSAV),ISAV3(IDSAV)
|
|
C INTEGER PATH(5000)
|
|
CHARACTER NAME1*6,NAME2*6
|
|
PARAMETER (MINOVR = 10, MINSLT = 3)
|
|
EXTERNAL PMIS
|
|
C
|
|
C the crap must be at least minovr bases in length or we dont bother
|
|
C
|
|
JFAIL = 1
|
|
IF((IDC.LT.MINOVR).OR.(IDIM2.LT.MINOVR)) RETURN
|
|
WRITE(NAME1,1002)LLINO
|
|
WRITE(NAME2,1002)IRNO
|
|
1002 FORMAT(I6)
|
|
NAME1(1:1) = 'C'
|
|
NAME2(1:1) = 'R'
|
|
IDC1 = MIN(IDC,IDIM2)
|
|
X = MAX(1.0,REAL(IDC1)/100.)
|
|
C
|
|
C allowing 10% and 7% of overlap to be pads: consensus and read respectively
|
|
C
|
|
MAXPC = NINT(10.0*X)
|
|
MAXPG = NINT(7.0*X)
|
|
IPP=IDSAV
|
|
CALL SLIDER(SEQ1,IDC,SEQ2,IDIM2,IC1,IG1,MAXPG,MAXPC,MINSLT,
|
|
+ISAV1,ISAV2,ISAV3,IPP)
|
|
IF((IPP.GT.IDSAV).OR.(IPP.LT.1)) THEN
|
|
CALL ERROM(KBOUT,'Warning: slider failed')
|
|
CALL ERROM(KBOUT,NAME2)
|
|
RETURN
|
|
END IF
|
|
CALL REMOVL(ISAV2,ISAV3,ISAV1,IPP)
|
|
CALL BUB3AS(ISAV2,ISAV3,ISAV1,IPP)
|
|
CALL TPCHEK(ISAV2,ISAV3,ISAV1,IPP)
|
|
CALL UPCHEK(ISAV2,ISAV3,ISAV1,IPP)
|
|
CALL LINEUP(SEQ2,SEQ1,SEQG2,SEQC2,IDC,IDIM2,IDOUT,ISAV3,ISAV2,
|
|
+ISAV1,IPP,ITOTPC,ITOTPG,JOINT,ITYPE,KBOUT,MAXGEL,JFAIL)
|
|
IF(JFAIL.NE.0) THEN
|
|
CALL ERROM(KBOUT,'Warning: lineup failed')
|
|
CALL ERROM(KBOUT,NAME2)
|
|
RETURN
|
|
END IF
|
|
C IDC1 = MIN(IDC,IDIM2+MAX(MAXPG,MAXPC))
|
|
C IDC1 = MIN(IDC,IDIM2)
|
|
C WRITE(*,*)'IDC,IDIM2',IDC,IDIM2
|
|
C MMAXPC = -1*(MAXPC+MAXPC/2)
|
|
C X = MAX(1.0,REAL(IDC1)/100.)
|
|
C MMAXPC = NINT(-10.*X)
|
|
C MMAXPG = NINT(7.*X)
|
|
C MAXPC = -MMAXPC
|
|
C MAXPG = MMAXPG
|
|
C IDIM3 = IDC1
|
|
C WRITE(*,*)'PADS ALLOWED',MAXPC,MAXPG
|
|
C CALL GALIGN(SEQ1,IDC1,SEQ2,IDIM3,MMAXPC,MMAXPG,PATH)
|
|
C CALL DISP(SEQ1,IDC1,SEQ2,IDIM3,PATH)
|
|
C CALL EXPAND(SEQ1,IDC1,SEQ2,IDIM3,SEQC2,IDC2,SEQG2,IDG2,PATH)
|
|
C LO = MIN(IDC2,IDG2)
|
|
C WRITE(*,*)'IDC2,IDG2',IDC2,IDG2
|
|
C LO = MIN(NOTIRL(SEQC2,IDC2,','),NOTIRL(SEQG2,IDG2,','))
|
|
C WRITE(*,*)'LO',LO
|
|
C END IF
|
|
LO = MIN(IDOUT,IDIM2)
|
|
PERCM = PMIS(SEQC2,SEQG2,LO)
|
|
ITOTPC = NSYM(SEQC2,LO,',')
|
|
ITOTPG = NSYM(SEQG2,LO,',')
|
|
IFAIL = 0
|
|
IF (ITOTPC.GT.MAXPC) IFAIL = 1
|
|
IF (ITOTPG.GT.MAXPG) IFAIL = 1
|
|
IF (PERCM.GT.PERMAX) IFAIL = 1
|
|
IF(IFAIL.NE.0) THEN
|
|
WRITE(IDEV,1000)PERCM,ITOTPC,ITOTPG
|
|
1000 FORMAT(' Percentage mismatch ',F4.1,', Pads',I4,I4)
|
|
CALL FMT4LP(SEQC2(1),SEQG2,LO,IPC,1,IDEV,
|
|
+ NAME1,NAME2)
|
|
CALL BUSY(KBOUT)
|
|
END IF
|
|
END
|
|
INTEGER FUNCTION NSYM(SEQ,ID,SYM)
|
|
CHARACTER SEQ(ID),SYM
|
|
N = 0
|
|
DO 10 I=1,ID
|
|
IF(SEQ(I).EQ.SYM) N = N + 1
|
|
10 CONTINUE
|
|
NSYM = N
|
|
END
|
|
REAL FUNCTION PMIS(SEQ1,SEQ2,ID)
|
|
CHARACTER SEQ1(ID),SEQ2(ID),PAD
|
|
PARAMETER (PAD = ',')
|
|
C
|
|
C count mismatch. padding in consensus ignored
|
|
C
|
|
J = 0
|
|
DO 10 I=1,ID
|
|
IF(SEQ1(I).NE.SEQ2(I)) THEN
|
|
IF(SEQ1(I).NE.PAD) J = J + 1
|
|
END IF
|
|
C IF(SEQ1(I).NE.SEQ2(I)) J = J + 1
|
|
10 CONTINUE
|
|
PMIS = 100. * REAL(J)/REAL(ID)
|
|
END
|
|
SUBROUTINE UPCHEK(PC,PG,L,N)
|
|
INTEGER PC(N),PG(N),L(N),DC,DG
|
|
C AUTHOR RODGER STADEN
|
|
C
|
|
C only allow gaps that are shorter than the next block of identity
|
|
C
|
|
K1 = 2
|
|
1 CONTINUE
|
|
DO 10 I = K1,N
|
|
J1 = I
|
|
DC = PC(I) - PC(I-1) - L(I-1)
|
|
DG = PG(I) - PG(I-1) - L(I-1)
|
|
IF(ABS(DC-DG).GE.L(I)) GO TO 20
|
|
10 CONTINUE
|
|
RETURN
|
|
20 CONTINUE
|
|
C WRITE(*,*)'REMOVING!!'
|
|
CALL ML(PC,PG,L,N,J1)
|
|
C IF(L(J1-1).GT.L(J1)) THEN
|
|
C CALL ML(PC,PG,L,N,J1)
|
|
C ELSE
|
|
C CALL ML(PC,PG,L,N,J1-1)
|
|
C END IF
|
|
K1 = MAX(2,J1-1)
|
|
N = N - 1
|
|
GO TO 1
|
|
END
|
|
SUBROUTINE ARRFIN(IDEV,SEQNCE,J,KBOUT,ICRAP)
|
|
CHARACTER TEMP(80),SEQNCE(J)
|
|
CHARACTER SPACE,TITCHR,GT
|
|
EXTERNAL NOK2
|
|
SAVE SPACE,TITCHR,GT
|
|
PARAMETER (IWIN = 100, NBAD = 5)
|
|
DATA SPACE/' '/
|
|
DATA TITCHR/';'/
|
|
DATA GT/'>'/
|
|
C
|
|
C routine to read reads, including their poor data if required
|
|
C allowing 5 n's in 100 bases
|
|
C
|
|
IDMX=J
|
|
J=0
|
|
1 CONTINUE
|
|
READ(IDEV,1001,END=30,ERR=80)TEMP
|
|
1001 FORMAT(80A1)
|
|
IF(TEMP(1).EQ.TITCHR) GO TO 1
|
|
DO 20 I=1,80
|
|
IF(TEMP(I).NE.SPACE)THEN
|
|
IF(J.EQ.IDMX)THEN
|
|
CALL ERROM(KBOUT,'Error in arrfin: too much data')
|
|
RETURN
|
|
END IF
|
|
J=J+1
|
|
SEQNCE(J)=TEMP(I)
|
|
END IF
|
|
20 CONTINUE
|
|
GO TO 1
|
|
30 CONTINUE
|
|
C WRITE(*,*)' GOOD',J
|
|
IF(ICRAP.NE.0) RETURN
|
|
C
|
|
C read in the crap
|
|
C
|
|
REWIND(IDEV)
|
|
JG = J
|
|
50 CONTINUE
|
|
READ(IDEV,1001,END=70,ERR=80)TEMP
|
|
IF(TEMP(1).NE.TITCHR) GO TO 50
|
|
IF(TEMP(2).EQ.GT) THEN
|
|
DO 60 I=3,80
|
|
IF(TEMP(I).NE.SPACE)THEN
|
|
IF(J.EQ.IDMX)THEN
|
|
CALL ERROM(KBOUT,'Error in arrfin: too much data')
|
|
GO TO 70
|
|
END IF
|
|
J=J+1
|
|
SEQNCE(J)=TEMP(I)
|
|
END IF
|
|
60 CONTINUE
|
|
END IF
|
|
GO TO 50
|
|
70 CONTINUE
|
|
C WRITE(*,*)' GOOD+BAD',J
|
|
K = NOK2(SEQNCE(JG+1),J-JG+1,IWIN,NBAD)
|
|
J = JG + K
|
|
C WRITE(*,*)' GOOD+BAD-CRAP',J
|
|
C WRITE(*,*)(SEQNCE(K),K=1,J)
|
|
RETURN
|
|
80 CONTINUE
|
|
CALL ERROM(KBOUT,'Error reading file in arrfin')
|
|
J = 0
|
|
END
|
|
SUBROUTINE DBLSTR(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,SEQ1,
|
|
+IDIM1,GEL,IDBSIZ,TEMP3,ID1,CHRSIZ,MAXGL2,KBIN,KBOUT,IDEVW,
|
|
+IDEVR,IDEV,LINLEN,PERCD,
|
|
+IHELPS,IHELPE,FILEH,IDEVH,MAXGEL,IDEVN,
|
|
+ LLINO,LINCON,LREG,RREG,MXGOOD,SEQ2,IDM,JOB,IERR,LSTRT,LEND)
|
|
CHARACTER FILEH*(*)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RREG, RELPG(IDBSIZ),CHRSIZ
|
|
INTEGER LREG,TEMP3(ID1,CHRSIZ,MAXGL2)
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
CHARACTER GEL(MAXGEL)
|
|
CHARACTER SEQ1(IDIM1),SEQ2(IDIM1)
|
|
C DBLSTR
|
|
C JOB 0 = Ask initial questions
|
|
C JOB 1 = Generate quality information
|
|
C JOB 2 = Complement contig and region info.
|
|
C
|
|
IF (JOB.EQ.0) THEN
|
|
CALL GETLN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+LINCON,LLINO,IERR,IDBSIZ,KBIN,KBOUT,IDEVN,
|
|
+IHELPS,IHELPE,FILEH,IDEVH)
|
|
IF(IERR.NE.0) RETURN
|
|
CALL GETREG(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+1,RELPG(LINCON),LREG,RREG,LINCON,LLINO,IDBSIZ,KBIN,KBOUT,
|
|
+ IHELPS,IHELPE,FILEH,IDEVH,IERR)
|
|
IF(IERR.NE.0) RETURN
|
|
ELSE IF (JOB.EQ.1) THEN
|
|
CALL SUMMAN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+SEQ1,IDIM1,GEL,LSTRT,LEND,LLINO,PERCD,IDBSIZ,
|
|
+TEMP3,ID1,CHRSIZ,MAXGL2,IDEVW,
|
|
+MAXGEL,MXGOOD)
|
|
IDIM1=LEND-LSTRT+1
|
|
CALL SUMMER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+SEQ2(LSTRT),IDIM1,GEL,LSTRT,LEND,LLINO,IDBSIZ,TEMP3,
|
|
+CHRSIZ,MAXGL2,IDEVW,MAXGEL,IDM,PERCD)
|
|
ELSE IF (JOB.EQ.2) THEN
|
|
C
|
|
C for complementary strand need to juggle the active region
|
|
C
|
|
LSTRT = LREG-LSTRT
|
|
LEND = LREG-LEND
|
|
IREG = RREG
|
|
RREG = RELPG(LINCON) - LREG + 1
|
|
LREG = RELPG(LINCON) - IREG + 1
|
|
LSTRT = LREG-LSTRT
|
|
LEND = LREG-LEND
|
|
IF (LSTRT.LT.1) LSTRT=1
|
|
IF (LEND.LT.1) LEND=1
|
|
C
|
|
CALL CMPLMT(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ LINCON,LNBR(LINCON),GEL,IDBSIZ,KBOUT,IDEVR,IDEVW,
|
|
+ MAXGEL)
|
|
LLINO = LNBR(LINCON)
|
|
CALL SUMMAN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ SEQ1,IDIM1,GEL,LSTRT,LEND,LLINO,PERCD,IDBSIZ,
|
|
+ TEMP3,ID1,CHRSIZ,MAXGL2,IDEVW,
|
|
+ MAXGEL,MXGOOD)
|
|
CALL SUMMER(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+ SEQ2(LSTRT),IDIM1,GEL,LSTRT,LEND,LLINO,IDBSIZ,TEMP3,
|
|
+ CHRSIZ,MAXGL2,IDEVW,MAXGEL,IDM,PERCD)
|
|
END IF
|
|
END
|
|
SUBROUTINE SUMMAN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,
|
|
+SEQ1,IDIM1,GEL,LREG,RREG,IGELC,PERCD,IDBSIZ,CHARS,
|
|
+ID1,CHRSIZ,MAXGL2,IDEVW,MAXGEL,MXGOOD)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),CHRSIZ
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
CHARACTER SEQ1(IDIM1)
|
|
CHARACTER GEL(MAXGEL)
|
|
INTEGER LREG,RREG,LSEQNO,POSN,Y,POSN1
|
|
INTEGER GELC
|
|
INTEGER CHARS(CHRSIZ,ID1,MAXGL2)
|
|
CHARACTER SYMS(7,7)
|
|
EXTERNAL INDEXS,LWRAPS,IGTCON
|
|
SAVE SYMS
|
|
DATA SYMS/
|
|
+'0','4','4','4','5','5','1',
|
|
+'4','0','4','4','5','5','1',
|
|
+'4','4','0','4','5','5','1',
|
|
+'4','4','4','0','5','5','1',
|
|
+'6','6','6','6','3','5','7',
|
|
+'6','6','6','6','6','3','7',
|
|
+'2','2','2','2','8','8','9'/
|
|
C
|
|
C codes are: using g = good consensus, b = bad consensus, n = no data
|
|
C
|
|
C + -
|
|
C 0 g g =
|
|
C 1 g n
|
|
C 2 n g
|
|
C 3 b b
|
|
C 4 g g !=
|
|
C 5 g b
|
|
C 6 b g
|
|
C 7 b n
|
|
C 8 n b
|
|
C 9 n n (should never occur!!!)
|
|
C
|
|
C hard to understand this very old code! rewrite it.
|
|
C we have a summing array of twice the length of the longest sequence
|
|
C posn is posn in contig of next to write, lseqno is current posn in contig
|
|
C we write when lseqno-posn ge the length of the longest seq
|
|
POSN=LREG
|
|
GELC=IGELC
|
|
LINLEN=MAXGEL
|
|
LSEQNO=RELPG(GELC)
|
|
IEND=0
|
|
DO 40 I=1,MAXGL2
|
|
DO 40 J=1,ID1
|
|
DO 40 K=1,CHRSIZ
|
|
CHARS(K,J,I)=0
|
|
40 CONTINUE
|
|
50 CONTINUE
|
|
ISS=1
|
|
IF(LNGTHG(GELC).LT.0)ISS=2
|
|
CALL READW(IDEVW,GELC,GEL,MAXGEL)
|
|
C LOOP FOR RELEVANT ELEMENTS THIS GEL
|
|
C only use mxgood characters from start of read
|
|
C
|
|
IF(ISS.EQ.1) THEN
|
|
N = MIN(MXGOOD,ABS(LNGTHG(GELC)))
|
|
IF(LSEQNO.LT.LREG)LSEQNO=LREG
|
|
IS = LSEQNO-RELPG(GELC)+1
|
|
ELSE
|
|
C LOOP FOR RELEVANT ELEMENTS THIS GEL
|
|
C only use mxgood characters from start of read (right end for these)
|
|
C
|
|
IREND = RELPG(GELC) - LNGTHG(GELC) + 1
|
|
IF (MXGOOD.LT.MAXGEL) THEN
|
|
LSEQNO = IREND - MXGOOD + 1
|
|
ELSE
|
|
LSEQNO = RELPG(GELC)
|
|
END IF
|
|
LSEQNO = MAX(LSEQNO,LREG)
|
|
IS = LSEQNO - RELPG(GELC) + 1
|
|
N = ABS(LNGTHG(GELC))
|
|
END IF
|
|
DO 70 I=IS,N
|
|
JJ = INDEXS(GEL(I),JSCORE)
|
|
JJJ = LWRAPS(LSEQNO,MAXGL2)
|
|
CHARS(JJ,ISS,JJJ) = CHARS(JJ,ISS,JJJ) + JSCORE
|
|
LSEQNO = LSEQNO + 1
|
|
70 CONTINUE
|
|
IF(RNBR(GELC).EQ.0)GO TO 200
|
|
GELC=RNBR(GELC)
|
|
LSEQNO=RELPG(GELC)
|
|
IF(LSEQNO.GT.RREG)GO TO 200
|
|
C ENOUGH TO OUTPUT?
|
|
Y=LSEQNO-POSN
|
|
IF(Y.GE.MAXGEL)GO TO 210
|
|
GO TO 50
|
|
200 CONTINUE
|
|
C SET FLAG TO SHOW END REACHED
|
|
IEND=1
|
|
LINLEN=MAXGEL
|
|
Y=RREG-POSN
|
|
IF(Y.LT.MAXGEL)LINLEN=Y+1
|
|
210 CONTINUE
|
|
C SET POINTER TO SEQ1
|
|
POSN1=POSN-1
|
|
C PREPARE NEXT SECTION OF CHARS FOR OUTPUT
|
|
DO 230 I=1,LINLEN
|
|
JJJ = LWRAPS(POSN,MAXGL2)
|
|
L1 = IGTCON(CHARS(1,1,JJJ),CHRSIZ,PERCD)
|
|
L2 = IGTCON(CHARS(1,2,JJJ),CHRSIZ,PERCD)
|
|
POSN1 = POSN1 + 1
|
|
SEQ1(POSN1) = SYMS(L2,L1)
|
|
DO 250 J=1,CHRSIZ
|
|
CHARS(J,1,JJJ)=0
|
|
CHARS(J,2,JJJ)=0
|
|
250 CONTINUE
|
|
POSN=POSN+1
|
|
230 CONTINUE
|
|
IF(POSN.GT.RREG)RETURN
|
|
IF((IEND.EQ.1).AND.(POSN.LE.RREG))GO TO 200
|
|
C ANY MORE MAXGEL CHAR LENGTHS TO OUTPUT
|
|
Y=LSEQNO-POSN
|
|
IF(Y.LT.MAXGEL)GO TO 50
|
|
C FINISHED COMPLETELY?
|
|
GO TO 210
|
|
END
|
|
INTEGER FUNCTION IGTCON(COUNTS,IDM,CUT)
|
|
INTEGER IDM
|
|
INTEGER COUNTS(IDM)
|
|
CHARACTER CHARSU
|
|
EXTERNAL CHARSU
|
|
C
|
|
C returns values 1 to 7: 7 means no data, 1-5 means a,c,g,t,*
|
|
C 6 means no consensus. if we change gtconc we should change this too
|
|
C
|
|
IGTCON = 7
|
|
ISUM = 0
|
|
DO 5 I=1,IDM
|
|
ISUM = ISUM + COUNTS(I)
|
|
5 CONTINUE
|
|
IF(ISUM.EQ.0) RETURN
|
|
Y = ISUM
|
|
DO 10 I = 1,IDM - 1
|
|
X = REAL(COUNTS(I))/Y
|
|
IF(X.GE.CUT) THEN
|
|
IGTCON = I
|
|
RETURN
|
|
END IF
|
|
10 CONTINUE
|
|
IGTCON = 6
|
|
END
|
|
SUBROUTINE CHKREV(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,NGELS,NCONTS,
|
|
+IDEVN,IDEV,KBIN,KBOUT,DONEIT,MATES,IGELS,DIST,RNAMES)
|
|
INTEGER RELPG(IDBSIZ),LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
INTEGER DONEIT(IDBSIZ),MATES(IDBSIZ),IGELS(IDBSIZ),DIST(IDBSIZ)
|
|
INTEGER GCLIN,DTOEND,CHAINL
|
|
CHARACTER*(*) RNAMES(IDBSIZ)
|
|
LOGICAL SAMEC
|
|
EXTERNAL GCLIN,DTOEND,NMATE,SAMEC,CHAINL
|
|
PARAMETER (MAXD = 99999)
|
|
C
|
|
C need to know which ones have already been mated
|
|
C
|
|
CALL FILLI(DONEIT,NGELS,0)
|
|
CALL FILLI(IGELS,NGELS,0)
|
|
CALL FILLI(MATES,NGELS,0)
|
|
CALL GNAMES(IDEVN,RNAMES,NGELS)
|
|
IGEL = 0
|
|
IDONE = 0
|
|
10 CONTINUE
|
|
IGEL = IGEL + 1
|
|
IF (IGEL.EQ.NGELS) THEN
|
|
IF (IDONE.GT.0) THEN
|
|
CALL SHOWP(DIST,IGELS,MATES,IDONE,RNAMES,RELPG,IDBSIZ,IDEV)
|
|
ELSE
|
|
WRITE(KBOUT,*)'None found'
|
|
END IF
|
|
RETURN
|
|
END IF
|
|
IF (DONEIT(IGEL).EQ.1) GO TO 10
|
|
MATE = NMATE(IGEL,NGELS,RNAMES,2)
|
|
IF (MATE.EQ.0) GO TO 10
|
|
DONEIT(MATE) = 1
|
|
IDONE = IDONE + 1
|
|
C
|
|
C we have a pair igel and mate. deal separately with cases of 1 or 2 contigs
|
|
C
|
|
IF (SAMEC(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,NGELS,NCONTS,IGEL,MATE))
|
|
+ THEN
|
|
C
|
|
C same contig
|
|
C
|
|
IF (((LNGTHG(IGEL).GT.0).AND.(LNGTHG(MATE).GT.0)).OR.
|
|
+ ((LNGTHG(IGEL).LT.0).AND.(LNGTHG(MATE).LT.0))) THEN
|
|
C
|
|
C >>>>>>> >>>>>>>>> or <<<<<<<< <<<<<<<<<<
|
|
C
|
|
IGELS(IDONE) = IGEL
|
|
MATES(IDONE) = MATE
|
|
DIST(IDONE) = MAXD
|
|
C WRITE(KBOUT,*)IGEL,MATE,
|
|
C + ' are in the same contig face the same way'
|
|
GO TO 10
|
|
END IF
|
|
IF (LNGTHG(IGEL).GT.0) THEN
|
|
IF (RELPG(IGEL).GT.RELPG(MATE)+ABS(LNGTHG(MATE)) -1 ) THEN
|
|
C
|
|
C <<<<<<<222 111>>>>>>>>
|
|
C
|
|
IGELS(IDONE) = IGEL
|
|
MATES(IDONE) = MATE
|
|
DIST(IDONE) = MAXD
|
|
C WRITE(KBOUT,*)IGEL,MATE,
|
|
C + ' are in the same contig and face away'
|
|
ELSE
|
|
C
|
|
C 111>>>>> <<<<<<222
|
|
C
|
|
IGELS(IDONE) = IGEL
|
|
MATES(IDONE) = MATE
|
|
DIST(IDONE) = RELPG(MATE)+ABS(LNGTHG(MATE))-RELPG(IGEL)
|
|
C WRITE(KBOUT,*)IGEL,RELPG(IGEL),MATE,RELPG(MATE)
|
|
END IF
|
|
ELSE
|
|
C
|
|
C <<<<<<111
|
|
C
|
|
IF (RELPG(MATE).GT.RELPG(IGEL)+ABS(LNGTHG(IGEL)) -1 ) THEN
|
|
C
|
|
C <<<<<111 222>>>>>
|
|
C
|
|
IGELS(IDONE) = IGEL
|
|
MATES(IDONE) = MATE
|
|
DIST(IDONE) = MAXD
|
|
C WRITE(KBOUT,*)IGEL,MATE,
|
|
C + ' are in the same contig and face away'
|
|
ELSE
|
|
C
|
|
C >>>>>222 <<<<<111
|
|
C
|
|
IGELS(IDONE) = IGEL
|
|
MATES(IDONE) = MATE
|
|
DIST(IDONE) = RELPG(IGEL)+ABS(LNGTHG(IGEL))-RELPG(MATE)
|
|
C WRITE(KBOUT,*)'same ',IGEL,RELPG(IGEL),MATE,RELPG(MATE)
|
|
END IF
|
|
END IF
|
|
ELSE
|
|
C
|
|
C different contigs
|
|
C
|
|
ICONTL = CHAINL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,IGEL)
|
|
ICONT = GCLIN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,ICONTL)
|
|
MCONTL = CHAINL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,MATE)
|
|
MCONT = GCLIN(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,MCONTL)
|
|
LENC1 = RELPG(ICONT)
|
|
LENC2 = RELPG(MCONT)
|
|
LENI = DTOEND(RELPG(IGEL),LNGTHG(IGEL),LENC1)
|
|
LENM = DTOEND(RELPG(MATE),LNGTHG(MATE),LENC2)
|
|
LENIM = LENI + LENM
|
|
IGELS(IDONE) = -IGEL
|
|
MATES(IDONE) = MATE
|
|
DIST(IDONE) = LENIM
|
|
C WRITE(KBOUT,*)'different contigs',IGEL,MATE,' separation',LENIM
|
|
END IF
|
|
GO TO 10
|
|
END
|
|
SUBROUTINE SHOWP(DIST,IGELS,MATES,IDONE,RNAMES,RELPG,IDBSIZ,IDEV)
|
|
INTEGER DIST(IDONE),IGELS(IDONE),MATES(IDONE),RELPG(IDBSIZ)
|
|
CHARACTER*(*) RNAMES(IDBSIZ)
|
|
CALL BUBBL3(DIST,IGELS,MATES,IDONE)
|
|
DO 10 I=1,IDONE
|
|
IF (IGELS(I).GT.0) THEN
|
|
J = ABS(IGELS(I))
|
|
WRITE(IDEV,1000)RNAMES(J),J,RELPG(J),DIST(I)
|
|
J = MATES(I)
|
|
WRITE(IDEV,1000)RNAMES(J),J,RELPG(J),DIST(I)
|
|
END IF
|
|
10 CONTINUE
|
|
DO 20 I=1,IDONE
|
|
IF (IGELS(I).LT.1) THEN
|
|
J = ABS(IGELS(I))
|
|
WRITE(IDEV,1000)RNAMES(J),J,RELPG(J),DIST(I),'*'
|
|
J = MATES(I)
|
|
WRITE(IDEV,1000)RNAMES(J),J,RELPG(J),DIST(I),'*'
|
|
END IF
|
|
20 CONTINUE
|
|
1000 FORMAT(' ',A,3I6,A)
|
|
END
|
|
INTEGER FUNCTION DTOEND(RELPG,LNGTHG,LENCON)
|
|
INTEGER RELPG
|
|
IF (LNGTHG.GT.0) THEN
|
|
DTOEND = LENCON - RELPG + 1
|
|
ELSE
|
|
DTOEND = RELPG + ABS(LNGTHG) - 1
|
|
END IF
|
|
END
|
|
LOGICAL FUNCTION SAMEC(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,
|
|
+NGELS,NCONTS,IGEL,JGEL)
|
|
INTEGER RELPG(IDBSIZ),LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
INTEGER CHAINL
|
|
EXTERNAL CHAINL
|
|
I = CHAINL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,IGEL)
|
|
J = CHAINL(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,IDBSIZ,JGEL)
|
|
SAMEC = .FALSE.
|
|
IF (I.EQ.J) SAMEC = .TRUE.
|
|
END
|
|
SUBROUTINE GNAMES(IDEVN,RNAMES,NGELS)
|
|
CHARACTER*(*) RNAMES(NGELS)
|
|
DO 10 I=1,NGELS
|
|
CALL READN(IDEVN,I,RNAMES(I))
|
|
10 CONTINUE
|
|
END
|
|
INTEGER FUNCTION NMATE(IGEL,NGELS,RNAMES,JOB)
|
|
PARAMETER (NAMLEN = 16)
|
|
CHARACTER*(NAMLEN) NAME,NAMEM,RNAMES(NGELS)
|
|
C
|
|
C assume s1 and r1 are mates
|
|
C
|
|
C if job = 1 we look at all the names
|
|
C if job = 2 we look only at the ones with higher numbers
|
|
NMATE = 0
|
|
IF(JOB.EQ.1) THEN
|
|
JGEL = 1
|
|
ELSE IF(JOB.EQ.2) THEN
|
|
JGEL = IGEL +1
|
|
ELSE
|
|
WRITE(*,*)'SCREAM: NMATE'
|
|
RETURN
|
|
END IF
|
|
NAME = RNAMES(IGEL)
|
|
I = INDEX(NAME,'.')
|
|
IF (NAME(I+1:I+1).EQ.'s') THEN
|
|
NAME(I+1:I+1) = 'r'
|
|
ELSE
|
|
NAME(I+1:I+1) = 's'
|
|
END IF
|
|
DO 10 I=JGEL,NGELS
|
|
IF (I.NE.IGEL) THEN
|
|
NAMEM = RNAMES(I)
|
|
IF (NAME.EQ.NAMEM) THEN
|
|
NMATE = I
|
|
RETURN
|
|
END IF
|
|
END IF
|
|
10 CONTINUE
|
|
END
|
|
SUBROUTINE LMATES(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,NGELS,
|
|
+IDEV,IGELS,MATES,DIST,RNAMES,NMATES)
|
|
INTEGER RELPG(IDBSIZ),LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
INTEGER MATES(IDBSIZ),IGELS(IDBSIZ),DIST(IDBSIZ)
|
|
CHARACTER*(*) RNAMES(IDBSIZ)
|
|
C
|
|
C list mates
|
|
C
|
|
C WRITE(*,*)NMATES
|
|
DO 10 I=1,NMATES
|
|
J = IGELS(I)
|
|
K = MATES(I)
|
|
IF(DIST(I).GT.0) THEN
|
|
WRITE(IDEV,1000)J,RNAMES(J),RELPG(J),LNGTHG(J),DIST(I)
|
|
WRITE(IDEV,1000)K,RNAMES(K),RELPG(K),LNGTHG(K)
|
|
1000 FORMAT(' ',I6,' ',A,' ',3I6)
|
|
END IF
|
|
10 CONTINUE
|
|
DO 20 I=1,NMATES
|
|
J = IGELS(I)
|
|
K = MATES(I)
|
|
IF(DIST(I).LT.0) THEN
|
|
WRITE(IDEV,1000)J,RNAMES(J),RELPG(J),LNGTHG(J),DIST(I)
|
|
WRITE(IDEV,1000)K,RNAMES(K),RELPG(K),LNGTHG(K)
|
|
END IF
|
|
20 CONTINUE
|
|
END
|
|
SUBROUTINE CHKREW(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,NGELS,
|
|
+IDEVN,IGELS,MATES,DIST,RNAMES,LLINO,NMATES)
|
|
INTEGER RELPG(IDBSIZ),LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
INTEGER MATES(IDBSIZ),IGELS(IDBSIZ),DIST(IDBSIZ)
|
|
CHARACTER*(*) RNAMES(IDBSIZ)
|
|
LOGICAL SAMEC
|
|
EXTERNAL NMATE,SAMEC
|
|
C
|
|
C need to know which ones have already been mated
|
|
C
|
|
CALL FILLI(IGELS,NGELS,0)
|
|
CALL FILLI(MATES,NGELS,0)
|
|
CALL GNAMES(IDEVN,RNAMES,NGELS)
|
|
IGEL = LLINO
|
|
NMATES = 0
|
|
10 CONTINUE
|
|
MATE = NMATE(IGEL,NGELS,RNAMES,1)
|
|
IF (MATE.NE.0) THEN
|
|
NMATES = NMATES + 1
|
|
IGELS(NMATES) = IGEL
|
|
MATES(NMATES) = MATE
|
|
C WRITE(*,*)IGELS(NMATES),MATES(NMATES)
|
|
C
|
|
C we have a pair igel and mate. deal separately with cases of 1 or 2 contigs
|
|
C
|
|
IF (SAMEC(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,NGELS,NCONTS,IGEL,MATE))
|
|
+ THEN
|
|
C
|
|
C same contig
|
|
C
|
|
DIST(NMATES) = RELPG(MATE) + ABS(LNGTHG(MATE)) - RELPG(IGEL)
|
|
IF(RELPG(MATE).LT.RELPG(IGEL)) THEN
|
|
C
|
|
C seen it before so skip it
|
|
C
|
|
NMATES = NMATES - 1
|
|
C IGELS(NMATES) = MATE
|
|
C MATES(NMATES) = IGEL
|
|
END IF
|
|
C WRITE(*,*)IGELS(NMATES),MATES(NMATES)
|
|
C WRITE(*,*)DIST(NMATES)
|
|
ELSE
|
|
C
|
|
C different contigs
|
|
C
|
|
DIST(NMATES) = -ABS(LNGTHG(IGELS(NMATES)))
|
|
C WRITE(*,*)DIST(NMATES)
|
|
END IF
|
|
END IF
|
|
IF (RNBR(IGEL).NE.0) THEN
|
|
IGEL = RNBR(IGEL)
|
|
GO TO 10
|
|
END IF
|
|
END
|
|
SUBROUTINE GDEPTH(RELPG,LNGTHG,
|
|
+IDBSIZ,LREG,RREG,IGELS,MATES,DIST,NMATES,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),IGELS(NMATES),MATES(NMATES),DIST(NMATES)
|
|
INTEGER LNGTHG(IDBSIZ)
|
|
INTEGER RREG,DEPTH
|
|
C WRITE(*,*)(IGELS(K),MATES(K),DIST(K),K=1,NMATES)
|
|
CALL GDPTH(RELPG,
|
|
+IDBSIZ,LREG,RREG,IGELS,MATES,DIST,NMATES,DEPTH)
|
|
IF(NMATES.LT.1) RETURN
|
|
CALL PLTCOG(RELPG,LNGTHG,IDBSIZ,IGELS,MATES,DIST,NMATES,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX,LREG,RREG,DEPTH)
|
|
END
|
|
SUBROUTINE GDPTH(RELPG,
|
|
+IDBSIZ,LREG,RREG,IGELS,MATES,DIST,NMATES,DEPTH)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),MATES(NMATES),IGELS(NMATES),DIST(NMATES)
|
|
INTEGER RREG,DEPTH
|
|
EXTERNAL NCDEPM
|
|
C LREG = left contig position
|
|
C RREG = right '' ''
|
|
C
|
|
C have nmates reading numbers in igels and mates, their spacings are dist
|
|
C dist <0 means they are in different contigs
|
|
C the igels numbers are in left to right order
|
|
C
|
|
I = 1
|
|
5 CONTINUE
|
|
IF(I.LE.NMATES) THEN
|
|
IF((RELPG(IGELS(I))+ABS(DIST(I))-1).LE.LREG) THEN
|
|
I = I + 1
|
|
GO TO 5
|
|
END IF
|
|
END IF
|
|
DEPTH = 0
|
|
10 CONTINUE
|
|
IF(I.LT.NMATES) THEN
|
|
IF(RELPG(IGELS(I)).LE.RREG) THEN
|
|
K = RELPG(IGELS(I)) + ABS(DIST(I))
|
|
DEPTH = MAX(NCDEPM(RELPG,IDBSIZ,I,IGELS,NMATES,K),DEPTH)
|
|
I = I + 1
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
C WRITE(*,*)'DEPTH',DEPTH
|
|
END
|
|
SUBROUTINE PLTCOG(RELPG,LNGTHG,IDBSIZ,IGELS,MATES,DIST,NMATES,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX,LREG,RREG,DEPTH)
|
|
INTEGER DEPTH
|
|
INTEGER RELPG(IDBSIZ),IGELS(NMATES),MATES(NMATES),DIST(NMATES)
|
|
INTEGER LNGTHG(IDBSIZ)
|
|
INTEGER RREG
|
|
C have window size margt starting at margb
|
|
C depths depthp, depthm
|
|
CALL VECTOM
|
|
CALL FRAME(MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
CALL PLTCG(RELPG,LNGTHG,
|
|
+IDBSIZ,LREG,RREG,IGELS,MATES,DIST,NMATES,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX,DEPTH)
|
|
CALL VT100M
|
|
END
|
|
SUBROUTINE PLTCG(RELPG,LNGTHG,
|
|
+IDBSIZ,LREG,RREG,IGELS,MATES,DIST,NMATES,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX,DEPTH)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),IGELS(NMATES),MATES(NMATES),DIST(NMATES)
|
|
INTEGER LNGTHG(IDBSIZ)
|
|
INTEGER RREG,DEPTH
|
|
CHARACTER SYM
|
|
LOGICAL MATEY
|
|
EXTERNAL MATEY
|
|
YINC = REAL(ISYMAX) / (DEPTH + 1)
|
|
XMIN = LREG
|
|
XMAX = RREG
|
|
YMAX = ISYMAX
|
|
YMIN = 0.
|
|
XSLOPE = REAL(MARGR - MARGL)/100.
|
|
YSLOPE = YINC/10.
|
|
C WRITE(*,*)XMIN,XMAX
|
|
I = 1
|
|
5 CONTINUE
|
|
IF(I.LE.NMATES) THEN
|
|
IF((RELPG(IGELS(I))+ABS(DIST(I))-1).LT.LREG) THEN
|
|
I = I + 1
|
|
GO TO 5
|
|
END IF
|
|
END IF
|
|
N = 0
|
|
10 CONTINUE
|
|
IF(I.LE.NMATES) THEN
|
|
IF(RELPG(IGELS(I)).LE.RREG) THEN
|
|
XF = MAX(RELPG(IGELS(I)),LREG)
|
|
XT = MIN(ABS(DIST(I))+RELPG(IGELS(I))-1,RREG)
|
|
N = N + 1
|
|
IF(N.GT.DEPTH) N = 1
|
|
YF = N * YINC
|
|
C WRITE(*,*)I,IGELS(I),MATES(I),DIST(I)
|
|
C WRITE(*,*)RELPG(IGELS(I)),LNGTHG(IGELS(I))
|
|
C WRITE(*,*)RELPG(MATES(I)),LNGTHG(MATES(I))
|
|
C WRITE(*,*)I,XF,XT
|
|
CALL LINE(XT,XF,YF,YF,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
IF (.NOT.(MATEY(IGELS(I),MATES(I),
|
|
+ RELPG,LNGTHG,IDBSIZ)).AND.(DIST(I).GT.0)) THEN
|
|
XF1 = XF + (XT-XF)/2.
|
|
YT1 = YF - YSLOPE
|
|
YT2 = YF + YSLOPE
|
|
CALL LINE(XF1,XF1,YT1,YT2,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
END IF
|
|
IF (LNGTHG(IGELS(I)).GT.0) THEN
|
|
XT = RELPG(IGELS(I))+LNGTHG(IGELS(I))
|
|
CALL RARROW(XT,YF,XSLOPE,YSLOPE,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
ELSE
|
|
XT = MAX(RELPG(IGELS(I)),LREG)
|
|
CALL LARROW(XT,YF,XSLOPE,YSLOPE,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
END IF
|
|
IF (DIST(I).GT.0) THEN
|
|
IF (LNGTHG(MATES(I)).GT.0) THEN
|
|
XT = RELPG(MATES(I))+LNGTHG(MATES(I))
|
|
CALL RARROW(XT,YF,XSLOPE,YSLOPE,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
ELSE
|
|
XT = MIN(RELPG(MATES(I)),RREG)
|
|
CALL LARROW(XT,YF,XSLOPE,YSLOPE,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
END IF
|
|
ELSE
|
|
SYM = '*'
|
|
CALL TEXT(XT,YF,SYM,1,ISIZE,XMAX,XMIN,YMAX,YMIN,
|
|
+ MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
END IF
|
|
I = I + 1
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
END
|
|
SUBROUTINE LARROW(XT,YF,XSLOPE,YSLOPE,XMAX,XMIN,YMAX,YMIN,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
XF1 = XT + XSLOPE
|
|
YT1 = YF + YSLOPE
|
|
XF2 = XT + XSLOPE
|
|
YT2 = YF - YSLOPE
|
|
CALL LINE(XT,XF1,YF,YT1,XMAX,XMIN,YMAX,YMIN,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
CALL LINE(XT,XF2,YF,YT2,XMAX,XMIN,YMAX,YMIN,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
END
|
|
SUBROUTINE RARROW(XT,YF,XSLOPE,YSLOPE,XMAX,XMIN,YMAX,YMIN,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
XF1 = XT - XSLOPE
|
|
YT1 = YF + YSLOPE
|
|
XF2 = XT - XSLOPE
|
|
YT2 = YF - YSLOPE
|
|
CALL LINE(XT,XF1,YF,YT1,XMAX,XMIN,YMAX,YMIN,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
CALL LINE(XT,XF2,YF,YT2,XMAX,XMIN,YMAX,YMIN,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX)
|
|
END
|
|
LOGICAL FUNCTION MATEY(IGEL,MATE,RELPG,LNGTHG,NGELS)
|
|
INTEGER RELPG(NGELS),LNGTHG(NGELS)
|
|
MATEY = .FALSE.
|
|
IF (RELPG(IGEL).LT.RELPG(MATE)) THEN
|
|
IF ((LNGTHG(IGEL).GT.0).AND.(LNGTHG(MATE).LT.0)) MATEY = .TRUE.
|
|
ELSE
|
|
IF ((LNGTHG(MATE).GT.0).AND.(LNGTHG(IGEL).LT.0)) MATEY = .TRUE.
|
|
END IF
|
|
END
|
|
INTEGER FUNCTION NCDEPM(RELPG,IDBSIZ,IGEL,
|
|
+IGELS,NMATES,RREG)
|
|
INTEGER RELPG(IDBSIZ)
|
|
INTEGER IGELS(NMATES)
|
|
INTEGER RREG
|
|
NCDEPM = 0
|
|
N = 0
|
|
I = IGEL
|
|
10 CONTINUE
|
|
IF(I.LE.NMATES) THEN
|
|
IF(RELPG(IGELS(I)).LE.RREG) THEN
|
|
N = N + 1
|
|
I = I + 1
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
NCDEPM = N
|
|
END
|
|
SUBROUTINE XHGAP(RELPG,LNGTHG,LNBR,RNBR,
|
|
+IDBSIZ,NCONTS,LLINOI,LINCNI,LREG,RREG,
|
|
+WINDOW,GWIND,LENCON,DEPTHP,DEPTHM,
|
|
+MARGL,MARGR,MARGB,MARGT,MAXOPT,ISXMAX,ISYMAX,KBIN,IDEV,
|
|
+KBOUT,GEL,GEL2,IDEV2,IDEV3,LINLEN,PERCD,MAXGEL,IDM,
|
|
+SEQ1,IDIM1,NGELS,TEMP3,CHRSIZ,MAXGL2,LINOU1,LINOU2,
|
|
+NOPT1,NOPT2,NOPT3,NOPT4,
|
|
+IHELPS,IHELPE,HELPF,IDEVH,MXGOOD,RNAMES,IGELS,MATES,DIST)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),WINDOW,CHRSIZ,GWIND
|
|
INTEGER LNGTHG(IDBSIZ),LNBR(IDBSIZ),RNBR(IDBSIZ)
|
|
INTEGER MARGB(MAXOPT),MARGT(MAXOPT)
|
|
INTEGER RREG,DEPTHP,DEPTHM,STRAND,CHNRP1,HQN
|
|
INTEGER TEMP3(2,CHRSIZ,MAXGL2)
|
|
CHARACTER GEL(MAXGEL),GEL2(MAXGEL)
|
|
CHARACTER TERM,TUPPER,NAMARC*16,HELPF*(*)
|
|
CHARACTER SEQ1(IDIM1),LINOU1(MAXGEL),LINOU2(MAXGEL)
|
|
CHARACTER*(*) RNAMES(IDBSIZ)
|
|
INTEGER MATES(IDBSIZ),IGELS(IDBSIZ),DIST(IDBSIZ),DEPTHR
|
|
EXTERNAL NOPWIN,CWORLD,TUPPER,CHNRP1,HQN
|
|
C nopt1 = single contig
|
|
C nopt2 = all contigs
|
|
C nopt3 = scan
|
|
C nopt4 = plot templates with pairs of reads
|
|
CALL CHKREW(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,NGELS,
|
|
+IDEV3,IGELS,MATES,DIST,RNAMES,LLINOI,NMATES)
|
|
10 CONTINUE
|
|
LLINO = LLINOI
|
|
LINCON = LINCNI
|
|
LOCLR = 0
|
|
LOCRR = 0
|
|
CALL BPAUSE(KBIN,KBOUT,IOK)
|
|
CALL CLEARV
|
|
CALL XHAIRR(ISXMAX,ISYMAX,IX,IY,TERM,DBTDUX,DBTDUY)
|
|
CALL VT100M
|
|
INFLAG = HQN(TERM)
|
|
IF(INFLAG.EQ.1) THEN
|
|
CALL HELP2(IHELPS,IHELPE,HELPF,IDEVH,KBIN,KBOUT)
|
|
GO TO 10
|
|
END IF
|
|
IF(INFLAG.EQ.2) RETURN
|
|
IF(INFLAG.EQ.3) RETURN
|
|
NOPT = NOPWIN(IY,MARGB,MARGT,MAXOPT)
|
|
C WRITE(*,*)NOPT
|
|
TERM = TUPPER(TERM)
|
|
IF(NOPT.EQ.0) RETURN
|
|
IF(NOPT.EQ.NOPT3) THEN
|
|
IF(TERM.EQ.'S') THEN
|
|
XMIN = LREG
|
|
XMAX = RREG
|
|
X = CWORLD(IX,MARGL,MARGR,XMIN,XMAX)
|
|
LOCLR = MAX(LREG,NINT(X)-WINDOW)
|
|
LOCRR = MIN(RREG,NINT(X)+WINDOW-1)
|
|
IF(LOCLR.NE.0) THEN
|
|
CALL DSPLAY(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ GEL,LLINO,LINCON,LOCLR,LOCRR,GEL2,I1,I2,0,I,
|
|
+ IDBSIZ,IDEV,KBOUT,
|
|
+ IDEV2,IDEV3,LINLEN,PERCD,MAXGEL,IDM)
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
IF((TERM.EQ.'N').OR.(TERM.EQ.'Z').OR.(TERM.EQ.'I')) GO TO 10
|
|
END IF
|
|
IF(NOPT.EQ.NOPT1) THEN
|
|
STRAND = 1
|
|
CALL FDPTH(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ IDBSIZ,LLINO,LREG,RREG,LENCON,STRAND,DEPTHP)
|
|
IF(DEPTHP.LT.0) RETURN
|
|
STRAND = -1
|
|
CALL FDPTH(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ IDBSIZ,LLINO,LREG,RREG,LENCON,STRAND,DEPTHM)
|
|
IF(DEPTHM.LT.0) RETURN
|
|
YMAX = ISYMAX
|
|
YMIN = 0.
|
|
XMIN = LREG
|
|
XMAX = RREG
|
|
RINC = ISYMAX / (DEPTHP + DEPTHM + 2)
|
|
RMID =(DEPTHM+1) * RINC
|
|
X = CWORLD(IX,MARGL,MARGR,XMIN,XMAX)
|
|
Y = CWORLD(IY,MARGB(NOPT),MARGT(NOPT),YMIN,YMAX)
|
|
C WRITE(*,*)X,Y,IX,IY
|
|
IF(TERM.EQ.'I') THEN
|
|
STRAND = 1
|
|
YINC = RINC * STRAND
|
|
CALL IPLTC(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ IDBSIZ,LLINO,LREG,RREG,STRAND,
|
|
+ MARGL,MARGR,MARGB(NOPT),MARGT(NOPT),
|
|
+ ISXMAX,ISYMAX,RMID,YINC,DEPTHP,X,Y,
|
|
+ KBOUT,IGEL,ICLOSE)
|
|
IF(ICLOSE.EQ.1) THEN
|
|
STRAND = -1
|
|
YINC = RINC * STRAND
|
|
CALL IPLTC(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ IDBSIZ,LLINO,LREG,RREG,STRAND,
|
|
+ MARGL,MARGR,MARGB(NOPT),MARGT(NOPT),
|
|
+ ISXMAX,ISYMAX,RMID,YINC,DEPTHM,X,Y,
|
|
+ KBOUT,IGEL,ICLOSE)
|
|
END IF
|
|
IF(ICLOSE.EQ.1) GO TO 10
|
|
CALL READN(IDEV3,IGEL,NAMARC)
|
|
WRITE(IDEV,1006)NAMARC,IGEL,RELPG(IGEL),LNGTHG(IGEL)
|
|
1006 FORMAT
|
|
+ ( ' Name ',A,' Number ',I6,' Rel. Posn. ',I7,' Length ',I5)
|
|
GO TO 10
|
|
END IF
|
|
IF(TERM.EQ.'Z') THEN
|
|
STRAND = 1
|
|
YINC = RINC * STRAND
|
|
CALL IPLTC(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ IDBSIZ,LLINO,LREG,RREG,STRAND,
|
|
+ MARGL,MARGR,MARGB(NOPT),MARGT(NOPT),
|
|
+ ISXMAX,ISYMAX,RMID,YINC,DEPTHP,X,Y,
|
|
+ KBOUT,IGEL,ICLOSE)
|
|
IF(ICLOSE.EQ.1) THEN
|
|
STRAND = -1
|
|
YINC = RINC * STRAND
|
|
CALL IPLTC(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ IDBSIZ,LLINO,LREG,RREG,STRAND,
|
|
+ MARGL,MARGR,MARGB(NOPT),MARGT(NOPT),
|
|
+ ISXMAX,ISYMAX,RMID,YINC,DEPTHM,X,Y,
|
|
+ KBOUT,IGEL,ICLOSE)
|
|
END IF
|
|
IF(IGEL.EQ.0) GO TO 10
|
|
CALL CLEARG
|
|
CALL PLC(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,LINCON,IGEL,
|
|
+ NCONTS,MARGL,MARGR,MARGB(NOPT2),MARGT(NOPT2),ISXMAX,ISYMAX)
|
|
LREG = MAX(1,RELPG(IGEL)-GWIND)
|
|
RREG = MIN(RELPG(LINCON),RELPG(IGEL)+GWIND)
|
|
LLINO = LNBR(LINCON)
|
|
LLINOI = LLINO
|
|
LINCNI = LINCON
|
|
LENCON = RREG - LREG + 1
|
|
CALL FDEPTH(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ IDBSIZ,LLINO,LREG,RREG,LENCON,
|
|
+ MARGL,MARGR,MARGB(NOPT1),MARGT(NOPT1),ISXMAX,ISYMAX)
|
|
GO TO 10
|
|
END IF
|
|
IF(TERM.EQ.'S') THEN
|
|
LOCLR = MAX(LREG,NINT(X)-WINDOW)
|
|
LOCRR = MIN(RREG,NINT(X)+WINDOW-1)
|
|
IF(LOCLR.NE.0) THEN
|
|
CALL DSPLAY(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ GEL,LLINO,LINCON,LOCLR,LOCRR,GEL2,I1,I2,0,I,
|
|
+ IDBSIZ,IDEV,KBOUT,
|
|
+ IDEV2,IDEV3,LINLEN,PERCD,MAXGEL,IDM)
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
IF(TERM.EQ.'N') THEN
|
|
LOCLR = MAX(LREG,NINT(X)-WINDOW)
|
|
LOCRR = MIN(RREG,NINT(X)+WINDOW-1)
|
|
IGEL = CHNRP1(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,
|
|
+ LLINO,LREG)
|
|
IF(LOCLR.NE.0) THEN
|
|
CALL LSTCON(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,IGEL,
|
|
+ LOCRR,IDEV,IDEV3,NAMARC)
|
|
END IF
|
|
GO TO 10
|
|
END IF
|
|
IF(TERM.EQ.'Q') THEN
|
|
CALL DBSCNP(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,SEQ1,
|
|
+ IDIM1,GEL,IDBSIZ,TEMP3,2,CHRSIZ,MAXGL2,IDEV2,LLINO,
|
|
+ PERCD,MAXGEL,LINOU1,LINOU2,LREG,RREG,
|
|
+ MARGL,MARGR,MARGB(NOPT3),MARGT(NOPT3),ISXMAX,ISYMAX,
|
|
+ MXGOOD)
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
IF(NOPT.EQ.NOPT2) THEN
|
|
CALL IDPLC(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,
|
|
+ NCONTS,IX,IY,MARGL,MARGR,MARGB(NOPT),MARGT(NOPT),
|
|
+ ISXMAX,ISYMAX,DBTDUX,DBTDUY,
|
|
+ LINCON,IGEL,IS)
|
|
IF(IGEL.EQ.0) RETURN
|
|
IF(TERM.EQ.'Z') THEN
|
|
CALL CLEARG
|
|
CALL PLC(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,LINCON,IGEL,
|
|
+ NCONTS,MARGL,MARGR,MARGB(NOPT),MARGT(NOPT),ISXMAX,ISYMAX)
|
|
LREG = 1
|
|
RREG = RELPG(LINCON)
|
|
LLINO = LNBR(LINCON)
|
|
LLINOI = LLINO
|
|
LINCNI = LINCON
|
|
LENCON = RREG - LREG + 1
|
|
CALL FDEPTH(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ IDBSIZ,LLINO,LREG,RREG,LENCON,
|
|
+ MARGL,MARGR,MARGB(NOPT1),MARGT(NOPT1),ISXMAX,ISYMAX)
|
|
GO TO 10
|
|
END IF
|
|
IF(TERM.EQ.'Q') THEN
|
|
CALL CLEARG
|
|
CALL PLC(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,LINCON,IGEL,
|
|
+ NCONTS,MARGL,MARGR,MARGB(NOPT),MARGT(NOPT),ISXMAX,ISYMAX)
|
|
LREG = 1
|
|
RREG = RELPG(LINCON)
|
|
LLINO = LNBR(LINCON)
|
|
LLINOI = LLINO
|
|
LINCNI = LINCON
|
|
LENCON = RREG - LREG + 1
|
|
CALL FDEPTH(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ IDBSIZ,LLINO,LREG,RREG,LENCON,
|
|
+ MARGL,MARGR,MARGB(NOPT1),MARGT(NOPT1),ISXMAX,ISYMAX)
|
|
CALL DBSCNP(RELPG,LNGTHG,LNBR,RNBR,NGELS,NCONTS,SEQ1,
|
|
+ IDIM1,GEL,IDBSIZ,TEMP3,2,CHRSIZ,MAXGL2,IDEV2,LLINO,
|
|
+ PERCD,MAXGEL,LINOU1,LINOU2,LREG,RREG,
|
|
+ MARGL,MARGR,MARGB(NOPT3),MARGT(NOPT3),ISXMAX,ISYMAX,
|
|
+ MXGOOD)
|
|
GO TO 10
|
|
END IF
|
|
IF(TERM.EQ.'I') THEN
|
|
CALL READN(IDEV3,IGEL,NAMARC)
|
|
WRITE(IDEV,1006)NAMARC,IGEL,RELPG(IGEL),LNGTHG(IGEL)
|
|
GO TO 10
|
|
END IF
|
|
IF(TERM.EQ.'S') THEN
|
|
LOCLR = MAX(1,IS-WINDOW)
|
|
LOCRR = MIN(RELPG(LINCON),IS+WINDOW-1)
|
|
LLINO = LNBR(LINCON)
|
|
IF(LOCLR.NE.0) THEN
|
|
CALL DSPLAY(RELPG,LNGTHG,LNBR,RNBR,
|
|
+ GEL,LLINO,LINCON,LOCLR,LOCRR,GEL2,I1,I2,0,I,
|
|
+ IDBSIZ,IDEV,KBOUT,
|
|
+ IDEV2,IDEV3,LINLEN,PERCD,MAXGEL,IDM)
|
|
END IF
|
|
GO TO 10
|
|
END IF
|
|
IF(TERM.EQ.'N') THEN
|
|
LOCLR = MAX(1,IS-WINDOW)
|
|
LOCRR = MIN(RELPG(LINCON),IS+WINDOW-1)
|
|
LLINO = LNBR(LINCON)
|
|
IF(LOCLR.NE.0) THEN
|
|
CALL LSTCON(RELPG,LNGTHG,LNBR,RNBR,IDBSIZ,IGEL,
|
|
+ LOCRR,IDEV,IDEV3,NAMARC)
|
|
END IF
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
IF(NOPT.EQ.NOPT4) THEN
|
|
CALL GDPTH(RELPG,
|
|
+ IDBSIZ,LREG,RREG,IGELS,MATES,DIST,NMATES,DEPTHR)
|
|
C + MARGL,MARGR,MARGB(NOPT),MARGT(NOPT),ISXMAX,ISYMAX)
|
|
IF(DEPTHR.LT.0) RETURN
|
|
YMAX = ISYMAX
|
|
YMIN = 0.
|
|
XMIN = LREG
|
|
XMAX = RREG
|
|
X = CWORLD(IX,MARGL,MARGR,XMIN,XMAX)
|
|
Y = CWORLD(IY,MARGB(NOPT),MARGT(NOPT),YMIN,YMAX)
|
|
C WRITE(*,*)X,Y,IX,IY
|
|
IF(TERM.EQ.'I') THEN
|
|
CALL IPLTCG(RELPG,LNGTHG,
|
|
+ IDBSIZ,LREG,RREG,IGELS,MATES,DIST,NMATES,
|
|
+ MARGL,MARGR,MARGB(NOPT),MARGT(NOPT),ISXMAX,ISYMAX,DEPTHR,X,Y,
|
|
+ IGEL,ICLOSE)
|
|
IF(ICLOSE.EQ.1) GO TO 10
|
|
C CALL READN(IDEV3,IGEL,NAMARC)
|
|
J = IGELS(IGEL)
|
|
WRITE(IDEV,1007)RNAMES(J),J,RELPG(J),LNGTHG(J)
|
|
J = MATES(IGEL)
|
|
WRITE(IDEV,1007)RNAMES(J),J,RELPG(J),LNGTHG(J)
|
|
1007 FORMAT(' 'A,3I6)
|
|
C WRITE(IDEV,1006)NAMARC,IGEL,RELPG(IGEL),LNGTHG(IGEL)
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
END
|
|
SUBROUTINE IPLTCG(RELPG,LNGTHG,
|
|
+IDBSIZ,LREG,RREG,IGELS,MATES,DIST,NMATES,
|
|
+MARGL,MARGR,MARGB,MARGT,ISXMAX,ISYMAX,DEPTH,X,Y,IGEL,IOK)
|
|
C AUTHOR: RODGER STADEN
|
|
INTEGER RELPG(IDBSIZ),IGELS(NMATES),MATES(NMATES),DIST(NMATES)
|
|
INTEGER LNGTHG(IDBSIZ)
|
|
INTEGER RREG,DEPTH
|
|
LOGICAL MATEY
|
|
EXTERNAL MATEY
|
|
YINC = REAL(ISYMAX) / (DEPTH + 1)
|
|
YINCO2 = YINC/2.
|
|
XMIN = LREG
|
|
XMAX = RREG
|
|
YMAX = ISYMAX
|
|
YMIN = 0.
|
|
C WRITE(*,*)X,Y
|
|
C WRITE(*,*)XMIN,XMAX,YMIN,YMAX,YINC
|
|
I = 1
|
|
5 CONTINUE
|
|
IF(I.LE.NMATES) THEN
|
|
IF((RELPG(IGELS(I))+ABS(DIST(I))-1).LT.LREG) THEN
|
|
I = I + 1
|
|
GO TO 5
|
|
END IF
|
|
END IF
|
|
N = 0
|
|
10 CONTINUE
|
|
IF(I.LE.NMATES) THEN
|
|
IF(RELPG(IGELS(I)).LE.RREG) THEN
|
|
XF = MAX(RELPG(IGELS(I)),LREG)
|
|
XT = MIN(ABS(DIST(I))+RELPG(IGELS(I))-1,RREG)
|
|
N = N + 1
|
|
IF(N.GT.DEPTH) N = 1
|
|
YF = N * YINC
|
|
C WRITE(*,*)XF,XT,YF
|
|
C WRITE(*,*)I,IGELS(I),MATES(I),DIST(I)
|
|
C WRITE(*,*)RELPG(IGELS(I)),LNGTHG(IGELS(I))
|
|
C WRITE(*,*)RELPG(MATES(I)),LNGTHG(MATES(I))
|
|
C WRITE(*,*)I,XF,XT
|
|
IF((X.GE.XF).AND.(X.LE.XT)) THEN
|
|
IF((Y.GE.YF-YINCO2).AND.(Y.LE.YF+YINCO2)) THEN
|
|
IOK = 0
|
|
IGEL = I
|
|
RETURN
|
|
END IF
|
|
END IF
|
|
I = I + 1
|
|
GO TO 10
|
|
END IF
|
|
END IF
|
|
IOK = 1
|
|
END
|