C ***************************************************************** C * * C * FORTRAN-77 PROGRAM * C * * C * EVALUATING STATISTICAL SIGNIFICANCE * C * IN TWO-STAGE GENOMEWIDE ASSOCIATION STUDIES * C * * C * JANUARY 5, 2006 * C * * C ***************************************************************** C C C C C ***************************************************************** C * * C * MAIN ROUTINE * C * * C ***************************************************************** C C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*70 TITLE,DAT,OUT DIMENSION A(530000000),ID(1000000) OPEN (5,FILE='gas2.dat') C READ (5,12) TITLE READ (5,12) DAT READ (5,12) OUT READ (5,*) N READ (5,*) M READ (5,*) N1 READ (5,*) C1 READ (5,*) NR C 12 FORMAT (A70) C ALPHA=0.5 IX=940 IY=439 IZ=1299 C OPEN (6,FILE=DAT) OPEN (7,FILE=OUT) C WRITE (7,100) WRITE (7,110) WRITE (7,120) WRITE (7,110) WRITE (7,130) WRITE (7,110) WRITE (7,100) WRITE (7,990) WRITE (7,990) WRITE (7,210) TITLE WRITE (7,990) WRITE (7,220) DAT WRITE (7,230) OUT WRITE (7,240) N WRITE (7,250) M WRITE (7,260) N1 WRITE (7,270) C1 WRITE (7,280) NR WRITE (7,990) WRITE (7,990) WRITE (7,300) WRITE (7,310) WRITE (7,320) C 100 FORMAT('* * * * * * * * * * * * * * * * * * * * * * * * * * * *') 110 FORMAT('* *') 120 FORMAT('* EVALUATING STATISTICAL SIGNIFICANCE *') 130 FORMAT('* IN TWO-STAGE GENOMEWIDE ASSOCIATION STUDIES *') 210 FORMAT(A50) 220 FORMAT('DATA FILE: ',A30) 230 FORMAT('OUTPUT FILE: ',A30) 240 FORMAT('TOTAL SAMPLE SIZE =',I8) 250 FORMAT('TOTAL NUMBER OF MARKERS =',I8) 260 FORMAT('SAMPLE SIZE IN STAGE 1 =',I8) 270 FORMAT('CRITICAL VALUE FOR STAGE 1 =',F10.5) 280 FORMAT('NUMBER OF MONTE CARLO SAMPLES =',I6) 300 FORMAT(T28,'UNADJUSTED',T48,'ADJUSTED P-VALUES') 310 FORMAT(T5,'MARKER',T13,'TEST STAT',T30,'P-VALUE', & T47,'HOLM',T62,'LIN') 320 FORMAT(T5,'------',T13,'---------',T27,'------------', & T44,'----------',T59,'----------') 990 FORMAT(' ') C LY=1+N*M LX0=LY+N LXB=LX0+M LU=LXB+M LUU=LU+N1*M LV=LUU+N1 LT=LV+M LP=LT+M LG1=LP+M LG2=LG1+N1 LTH1=LG2+N1 LTH2=LTH1+M*NR C CALL INPUT(N,M,A,A(LY)) CALL TEST(N1,N,M,NR,IX,IY,IZ,C1,ALPHA,A,A(LY),A(LXB), & A(LU),A(LUU),A(LV),A(LT), & A(LP),A(LG1),A(LG2),A(LTH1),A(LTH2),ID) C STOP END C C C C ***************************************************************** C * * C * INPUT SUBROUTINE * C * * C ***************************************************************** C C SUBROUTINE INPUT(N,M,X,Y) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION X(N,M),Y(N) C DO 100 I=1,N READ (6,*) Y(I),(X(I,J),J=1,M) 100 CONTINUE C RETURN END C C C C ***************************************************************** C * * C * TEST SUBROUTINE * C * * C ***************************************************************** C C SUBROUTINE TEST(N1,N,M,NR,IX,IY,IZ,C1,ALPHA,X,Y,XB, & U,UU,V,T,P,G1,G2,TH1,TH2,ID) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION X(N,M),Y(N),XB(M),U(N1,M),UU(N1),V(M),T(M),P(M), & G1(N1),G2(N1),TH1(M,NR),TH2(M,NR),ID(M) C N2=N-N1 C YB=0.D0 C DO 10 I=1,N1 YB=YB+Y(I) 10 CONTINUE C YB=YB/DFLOAT(N1) C DO 12 J=1,M X0=0.D0 XB(J)=0.D0 DO 11 I=1,N1 IF (X(I,J).GT.8.D0) GOTO 11 X0=X0+1.D0 XB(J)=XB(J)+X(I,J) 11 CONTINUE XB(J)=XB(J)/X0 12 CONTINUE C DO 20 J=1,M ID(J)=J T(J)=0.D0 V(J)=0.D0 DO 15 I=1,N1 IF (X(I,J).GT.8.D0) THEN U(I,J)=0.D0 ELSE U(I,J)=(Y(I)-YB)*(X(I,J)-XB(J)) T(J)=T(J)+U(I,J) V(J)=V(J)+U(I,J)*U(I,J) ENDIF 15 CONTINUE IF (V(J).LT.1.D-10) V(J)=1.D-10 T(J)=T(J)*T(J)/V(J) 20 CONTINUE C YB=0.D0 DO 30 I=1,N YB=YB+Y(I) 30 CONTINUE YB=YB/DFLOAT(N) C DO 50 J=1,M IF (T(J).LT.C1) THEN T(J)=0.D0 ELSE X0=0.D0 XB(J)=0.D0 DO 40 I=1,N IF (X(I,J).LT.8.D0) THEN X0=X0+1.D0 XB(J)=XB(J)+X(I,J) ENDIF 40 CONTINUE XB(J)=XB(J)/X0 T(J)=0.D0 V(J)=0.D0 DO 45 I=1,N IF (X(I,J).GT.8.D0) GOTO 45 UJ=(Y(I)-YB)*(X(I,J)-XB(J)) T(J)=T(J)+UJ V(J)=V(J)+UJ*UJ 45 CONTINUE IF (V(J).LT.1.D-10) V(J)=1.D-10 T(J)=T(J)*T(J)/V(J) ENDIF 50 CONTINUE C CALL SORT(M,N1,T,U,ID,UU) C DO 60 J=1,M P(J)=GAMMQ(0.5D0,T(J)/2.D0) V(J)=0.D0 DO 55 I=1,N1 V(J)=V(J)+U(I,J)*U(I,J) 55 CONTINUE IF (V(J).LT.1.D-10) V(J)=1.D-10 60 CONTINUE C DO 250 L=1,NR C DO 160 I=1,N1 G1(I)=GASDEV(IX,IY,IZ,GSET) G2(I)=DSQRT(DFLOAT(N2)/DFLOAT(N1))*GASDEV(IX,IY,IZ,GSET) 160 CONTINUE C DO 200 J=1,M TH1(J,L)=0.D0 TH2(J,L)=0.D0 DO 180 I=1,N1 TH1(J,L)=TH1(J,L)+U(I,J)*G1(I) TH2(J,L)=TH2(J,L)+U(I,J)*G2(I) 180 CONTINUE TH2(J,L)=TH1(J,L)+TH2(J,L) TH1(J,L)=TH1(J,L)*TH1(J,L)/V(J) TH2(J,L)=TH2(J,L)*TH2(J,L)/(DFLOAT(N)*V(J)/DFLOAT(N1)) 200 CONTINUE C 250 CONTINUE C PH=0.D0 PM=0.D0 C DO 380 NJ=M,1,-1 R=0.D0 DO 360 L=1,NR DO 340 J=1,NJ IF ((TH1(J,L).GT.C1).AND.(TH2(J,L).GT.T(NJ))) THEN R=R+1.D0 GOTO 360 ENDIF 340 CONTINUE 360 CONTINUE PT=P(NJ)*DFLOAT(NJ) IF (PT.GT.PH) PH=PT IF (PH.GT.1.D0) PH=1.D0 PT=R/DFLOAT(NR) IF (PT.GT.PM) PM=PT IF (PM.GT.ALPHA) GOTO 990 WRITE (7,900) ID(NJ),T(NJ),P(NJ),PH,PM C 380 CONTINUE C 900 FORMAT(I10,F10.3,F18.10,2F15.8) C 990 RETURN END C C C C ***************************************************************** C * * C * SORT SUBROUTINE * C * * C ***************************************************************** C C C SUBROUTINE SORT(N,NP,X,Z,ID,ZZ) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION X(N),ID(N),Z(NP,N),ZZ(NP) L=N/2+1 IR=N 10 CONTINUE IF (L.GT.1) THEN L=L-1 XX=X(L) IDD=ID(L) DO 100 K=1,NP ZZ(K)=Z(K,L) 100 CONTINUE ELSE XX=X(IR) IDD=ID(IR) DO 110 K=1,NP ZZ(K)=Z(K,IR) 110 CONTINUE X(IR)=X(1) ID(IR)=ID(1) DO 120 K=1,NP Z(K,IR)=Z(K,1) 120 CONTINUE IR=IR-1 IF (IR.EQ.1) THEN X(1)=XX ID(1)=IDD DO 130 K=1,NP Z(K,1)=ZZ(K) 130 CONTINUE RETURN ENDIF ENDIF I=L J=L+L 20 IF (J.LE.IR) THEN IF (J.LT.IR) THEN IF (X(J).LT.X(J+1)) J=J+1 ENDIF IF (XX.LT.X(J)) THEN X(I)=X(J) ID(I)=ID(J) DO 150 K=1,NP Z(K,I)=Z(K,J) 150 CONTINUE I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF X(I)=XX ID(I)=IDD DO 160 K=1,NP Z(K,I)=ZZ(K) 160 CONTINUE GO TO 10 END C C C C ***************************************************************** C * * C * UNIFORM DERIVATIVE SUBROUTINE * C * * C ***************************************************************** C C C DOUBLE PRECISION FUNCTION UNIDEV(IX,IY,IZ) IX = 171*MOD(IX,177) - 2*(IX/177) IY = 172*MOD(IY,176) - 35*(IY/176) IZ = 170*MOD(IZ,178) - 63*(IZ/178) IF(IX.LT.0) IX = IX+30269 IF(IY.LT.0) IY = IY+30307 IF(IZ.LT.0) IZ = IZ+30323 UNIDEV = DMOD(DFLOAT(IX)/30269.D0 + DFLOAT(IY)/30307.D0 + & DFLOAT(IZ)/30323.D0, 1.D0) RETURN END C C C C **************************************************************** C * * C * GAUSSIAN DERIVATIVE SUBROUTINE * C * * C **************************************************************** C C DOUBLE PRECISION FUNCTION GASDEV(IX,IY,IZ,GSET) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DATA ISET/0/ IF (ISET.EQ.0) THEN 1 V1=2.D0*UNIDEV(IX,IY,IZ)-1.D0 V2=2.D0*UNIDEV(IX,IY,IZ)-1.D0 R=V1**2+V2**2 IF (R.GE.1.D0) GO TO 1 FAC=DSQRT(-2.D0*DLOG(R)/R) GSET=V1*FAC GASDEV=V2*FAC ISET=1 ELSE GASDEV=GSET ISET=0 ENDIF RETURN END C C C C ***************************************************************** C * * C * GAMMQ FUNCTION SUBROUTINE * C * * C ***************************************************************** C C C DOUBLE PRECISION FUNCTION GAMMQ(A,X) IMPLICIT DOUBLE PRECISION(A-H,O-Z) IF ((X.LT.0.D0).OR.(A.LE.0.D0)) PAUSE IF (X.LT.(A+1.D0)) THEN CALL GSER(GAMSER,A,X,GLN) GAMMQ=1.D0-GAMSER ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMQ=GAMMCF ENDIF RETURN END C C C C ***************************************************************** C * * C * GSER FUNCTION SUBROUTINE * C * * C ***************************************************************** C C C SUBROUTINE GSER(GAMSER,A,X,GLN) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER (ITMAX=100,EPS=3.D-7) GLN=GAMMLN(A) IF (X.LE.0.D0) THEN IF (X.LT.0.D0) PAUSE GAMSER=0.D0 RETURN ENDIF AP=A SUM=1.D0/A DEL=SUM DO 11 N=1,ITMAX AP=AP+1.D0 DEL=DEL*X/AP SUM=SUM+DEL IF (DABS(DEL).LT.(DABS(SUM)*EPS)) GO TO 1 11 CONTINUE PAUSE 'A too large, ITMAX too small' 1 GAMSER=SUM*DEXP(-X+A*DLOG(X)-GLN) RETURN END C C C C ***************************************************************** C * * C * GCF FUNCTION SUBROUTINE * C * * C ***************************************************************** C C C SUBROUTINE GCF(GAMMCF,A,X,GLN) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER (ITMAX=100,EPS=3.D-7) GLN=GAMMLN(A) GOLD=0.D0 A0=1.D0 A1=X B0=0.D0 B1=1.D0 FAC=1.D0 DO 11 N=1,ITMAX AN=DFLOAT(N) ANA=AN-A A0=(A1+A0*ANA)*FAC B0=(B1+B0*ANA)*FAC ANF=AN*FAC A1=X*A0+ANF*A1 B1=X*B0+ANF*B1 IF (A1.NE.0.D0) THEN FAC=1.D0/A1 G=B1*FAC IF (DABS((G-GOLD)/G).LT.EPS) GO TO 1 GOLD=G ENDIF 11 CONTINUE PAUSE 'A too large, ITMAX too small' 1 GAMMCF=DEXP(-X+A*DLOG(X)-GLN)*G RETURN END C C C C ***************************************************************** C * * C * GAMMLN FUNCTION SUBROUTINE * C * * C ***************************************************************** C C C DOUBLE PRECISION FUNCTION GAMMLN(XX) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION COF(6) DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, * -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/ DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/ X=XX-ONE TMP=X+FPF TMP=(X+HALF)*DLOG(TMP)-TMP SER=ONE DO 11 J=1,6 X=X+ONE SER=SER+COF(J)/X 11 CONTINUE GAMMLN=TMP+DLOG(STP*SER) RETURN END