PROGRAM impt_anct IMPLICIT NONE INTEGER, PARAMETER :: nt=3682500,ng=25,nimp=10000 INTEGER, PARAMETER :: itg=SELECTED_INT_KIND (2) !!$ INTEGER, PARAMETER :: cod=SELECTED_INT_KIND (4) !!$ INTEGER, PARAMETER :: id=SELECTED_INT_KIND (6) INTEGER, PARAMETER :: doub=SELECTED_REAL_KIND (14) INTEGER :: dob1,dob2,ir4,ir5,ir6,nr,pgrp INTEGER :: a,i,j,k,l,s, nl,nbit,ip,im,ii, & npop,ng1, nfond, io INTEGER, ALLOCATABLE :: act(:) REAL(KIND=doub) :: x,y,r,z,sum,xr1,xr2,xmin,xmaj !!$ real :: t1,t2,gtimer INTEGER :: indi(nimp),pp(nimp),mm(nimp) INTEGER (KIND=itg) :: nxgen(0:nt),ncgen(0:nt) INTEGER :: dob,ann(nt),ndes(nt),eff(0:ng),prc(nt),mrc(nt),ord(nt) REAL(KIND=doub) :: prob(nt), pro(nt),pr1(nimp),pr2(nimp),cum,aj(nt)!,dgr(nt) CHARACTER (LEN=5) :: refbd,sub !, breed(nt) LOGICAL :: t, pop(nt), donact(0:nt) CHARACTER :: pedfile*50,nfile*50, date*8,time*10 1101 FORMAT (/'No. of individuals:',i8) 1102 FORMAT ('No. of non parents:',i8) 1103 FORMAT ('No. of whole founders in entire pop :',i8) 1104 FORMAT ('Problems in pedigree file ',i8) 1105 FORMAT ('No. of animals in ref pop :', i8/ & 'Max. no. of gen. traced : ', i2) 1106 FORMAT ('No. of founders of ref pop :',i8) 1107 FORMAT ('Eff. no. of founders : ',f10.6) OPEN(5, FILE='anctpgrp.in', STATUS='OLD') READ(5,*) pgrp,pedfile !,popfile,nfile PRINT '(a/a,a)','------------------------------------------------', & 'Input files: pedigree= ',TRIM(pedfile) CALL DATE_AND_TIME (date,time) OPEN(9, FILE='anctSARD.'//date, STATUS='new') PRINT *, 'Eff. no. of Ancestors from anctSARD.f90 on: ',date, ' ', time WRITE(9,*) 'Eff. no. of Ancestors from anctSARD.f90 on: ',date, ' ', time WRITE(9,'(a/a,a)') '------------------------------------------------', & 'Input pedigree file: ',TRIM(pedfile) READ(5,*) dob1,dob2,refbd IF (dob1>dob2) THEN PRINT *, 'Error in parameter file:', dob1, dob2 PRINT *, 'First birthday later than last birthday' STOP ENDIF PRINT 109,nimp,pgrp,refbd,dob1,dob2 WRITE(9,109) nimp,pgrp,refbd,dob1,dob2 109 FORMAT (/'Analysis of important ancestors' & /'No. of ancestors to be searched :',i6 & /'No. of parent groups (nullified):',i6 & /'Breed of the reference pop : ',a5 & /'Birthday range of reference pop :',2i9) READ(5,*) a ALLOCATE (act(1:a)) DO i=1,a READ(5,*) act(i) ENDDO PRINT *, 'Number of carrier ancestors declared:',a ir6=0 indi=0 pp=0 mm=0 nxgen=0 ncgen=0 ndes=0 donact=.TRUE. eff=0 ord=0 pro=0.d0 pr1=0.d0 pr2=0.d0 cum=0.d0 aj=0.d0 !!$Reading pedigree and ref pop OPEN (1,file=pedfile,STATUS='old',FORM='FORMATTED',IOSTAT=io) IF (io > 0) THEN PRINT *,'Cannot find input pedigree file, error code:',io PRINT *,'Pls check location/name of file: ',TRIM(pedfile) STOP ENDIF !1001 FORMAT (4(i8)) nl=0 npop=0 ng1=0 pop=.FALSE. prc=0 mrc=0 ann=0 prob=0.d0 DO READ(1,*,IOSTAT=io) i,ip,im,dob,sub IF (io==-1) THEN WRITE (9,*)'End of record in ',pedfile WRITE (9,*) i,'ip=',ip,'im=',im,'dob=',dob,'sub=',sub EXIT ENDIF ! i=i-pgrp ! ip=ip-pgrp ! im=im-pgrp ! IF (i<=0) CYCLE IF (ip<=pgrp) ip=0 IF (im<=pgrp) im=0 nl=nl+1 IF (nl>nt) THEN WRITE(5,*) 'Memory limited to nt:',nt EXIT ! STOP 'Increase nt' ENDIF IF (i/=nl) THEN PRINT *, i,'ip=',ip,'im=',im,'dob=',dob PRINT *,'File not sorted according to individuals','i=',i,'nl=',nl STOP ENDIF prc(nl)=ip mrc(nl)=im ann(nl)=dob !!$ breed(nl)=sub !!$ Define ref population and initialise prob IF (ann(i)>=dob1 .AND. ann(i)<=dob2 & .AND. prc(i)/=0 .AND. mrc(i)/=0 .AND. TRIM(sub)==TRIM(refbd)) THEN prob(i)=1.d0 pop(i)=.TRUE. npop=npop+1 ng1=MAX(ng1,nxgen(i)) ENDIF ENDDO CLOSE(1) PRINT *, 'Last line from ',pedfile PRINT *, sub,nl,prc(nl),mrc(nl),ann(nl) PRINT 1101, nl WRITE(9,*) 'Last line from ',pedfile WRITE(9,*) nl,prc(nl),mrc(nl),ann(nl) WRITE(9,1101) nl PRINT 1105, npop,ng1 WRITE(9,1105) npop,ng1 PRINT *, 'Identities and pedigree of carrier ancestors' WRITE(9,*) 'Identities and pedigree of carrier ancestors' DO i=1,a donact(act(i))=.FALSE. PRINT *, act(i),prc(act(i)),mrc(act(i)) WRITE(9,'(3i7)') act(i),prc(act(i)),mrc(act(i)) ENDDO !!$ compute number of progeny k=0 ndes=0 DO i=1,nl IF (prc(i)/=0) ndes(prc(i))=ndes(prc(i)) + 1 IF (mrc(i)/=0) ndes(mrc(i))=ndes(mrc(i)) + 1 ENDDO DO i=1,nl IF (ndes(i)==0) k=k+1 ENDDO PRINT 1102, k WRITE(9,1102) k !!$ Assignment of values for nxgen and ord !!$ nxgen(i)=0 for 'whole' founders, -1 for otherwise !!$ ord(k) is for reordering of pedigree (ord(k) not necessary for us) k=0 DO i=1,nl nxgen(i)=-1 IF (prc(i)==0 .AND. mrc(i)==0) THEN nxgen(i)=0 k=k+1 ord(k)=i ENDIF ENDDO PRINT 1103, k WRITE(9,1103) k !!$ Reordering of pedigree (Max no. of loops needed is max. no. of gen., ng) PRINT *, 'Reordering pedigree' nbit=0 DO WHILE(k0) STOP !!$ Count the no. of individuals, eff(nxgen(i)), tracing nxgen(i) no. of gen. !!$ to its earliest founder. ! OPEN (10,file=nfile,FORM='FORMATTED',IOSTAT=io) ! IF (io > 0) THEN ! PRINT *,'Cannot find output pedigree file, error code:',io ! PRINT *,'Pls check file access: ',TRIM(nfile) ! STOP ! ENDIF eff=0 DO i=1,nl ii=nxgen(i) eff(ii)=eff(ii)+1 ! WRITE(10,2010) i,ann(i),ncgen(i),nxgen(i),ndes(i) ENDDO !2010 FORMAT(i7,i5,2i3,i6) PRINT *, 'Max no. of generations traced by individual' WRITE(9,*) 'Max no. of generations traced by individual' DO i=1, MAXVAL(nxgen) PRINT *, i,eff(i) WRITE(9,*) i,eff(i) ENDDO !!$ Calculate the gene contribution of parents of individuals in ref pop !!$ First DO loop points to the individual id of the last animal ordered DO k=nl,1,-1 i=ord(k) IF (pop(i)) THEN ip=prc(i) im=mrc(i) IF (ip/=0) prob(ip)=prob(ip) + prob(i)*.5d0 IF (im/=0) prob(im)=prob(im) + prob(i)*.5d0 ENDIF ENDDO !!$ Tracing the pedigree further back DO k=nl,1,-1 i=ord(k) IF (.NOT.pop(i)) THEN ip=prc(i) im=mrc(i) IF (ip/=0) prob(ip)=prob(ip) + prob(i)*.5d0 IF (im/=0) prob(im)=prob(im) + prob(i)*.5d0 ENDIF ENDDO !!$ Calculating and summing the eff gene contrib. of each founder to the ref pop nfond=0 x=0.d0 DO i=1,nl IF (prc(i)==0 .AND. mrc(i)==0 .AND. .NOT.pop(i) .AND. prob(i)>0.d0) THEN nfond=nfond+1 y=prob(i)/REAL((npop),KIND=doub) x=x + y*y ELSEIF ((prc(i)==0 .OR. mrc(i)==0) .AND. .NOT.pop(i) .AND. prob(i)>0.d0) THEN nfond=nfond+1 y=prob(i)*.5d0/REAL((npop),KIND=doub) x=x + y*y ENDIF ENDDO x=1.d0/x PRINT 1106,nfond PRINT 1107,x WRITE(9,1106) nfond WRITE(9,1107) x !!$ *********Analysis of Important Ancestors********** 1108 FORMAT (/' # Id Contrib sire dam Nxgen Progeny Lower Upper' & /' total marginal cumulated bound') PRINT 1108 WRITE (9,1108) !!$ Initialisation cum=0.d0 sum=0.d0 pro=prob l=0 x=0.d0 t=.TRUE. !!$ $$$$$$$$$$$$ analyse du lme fondateur $$$$$$$$$$$$$$$$$$$$ DO WHILE(l<=nimp .AND. cum<1.d0) IF (t .AND. a>0) THEN DO i=1,a IF (.NOT. donact(act(i))) THEN t=.FALSE. PRINT *,l,donact(act(i)),act(i),i EXIT ENDIF ENDDO ENDIF IF (t) THEN PRINT *, 'All carrier ancestors identified' donact=.TRUE. DO i=1,a donact(act(i))=.FALSE. ENDDO i=a DO j=1,nl IF ((.NOT. donact(prc(j)) .OR. .NOT. donact(mrc(j))) & .AND. .NOT. pop(j)) THEN i=i+1 donact(j)=.FALSE. ! IF (prc(j)==act(1)) PRINT *,j,donact(j) ENDIF ENDDO PRINT *, '# potential carrier ancestors',i WRITE(9,*) '# potential carrier ancestors',i DO i=1,a donact(act(i))=.TRUE. ENDDO t=.FALSE. a=0 ELSEIF (a>0) THEN t=.TRUE. ENDIF l=l+1 !!$ IF (l==11) then !!$ gtimer is a unix specific function !!$ t2=gtimer() !!$ t1=0. !!$ print *,'CPU per founder : ', (t2-t1)*.1 !!$ ENDIF !!$ a:recherche de pro maxi DO i=1,l-1 pro(indi(i))=0.d0 ENDDO x=0.d0 DO i=1,nl IF (.NOT. donact(i) .AND. x1) THEN r=x/REAL((npop),KIND=doub) nr=INT((1.d0-cum)/r) z=1.d0 - cum - REAL((nr),KIND=doub)*r xmin=1.d0 / (sum + REAL((nr),KIND=doub)*r*r + z*z) !!$ print *,x,npop,r,cum,nr,z,xmin nr=nfond - l + 1 r=(1.d0-cum)/ REAL((nr),KIND=doub) xmaj=1.d0/(sum + REAL((nr),KIND=doub)*r*r) ! DO i=1,a ! IF(indi(l-1)==act(i)) THEN ! donact(act(i))=.TRUE. ! ENDIF ! ENDDO PRINT *, l-1,indi(l-1),xr1,xr2,cum,pp(l-1),mm(l-1),ndes(indi(l-1)) WRITE (9,1109) l-1,indi(l-1),xr1,xr2,cum,pp(l-1),mm(l-1), & nxgen(indi(l-1)),ndes(indi(l-1)),xmin,xmaj IF (indi(l-1)==indi(l-2)) EXIT ENDIF 1109 FORMAT(i4,i8,3f9.6,2i8,i3,i6,2f7.2) !!$ b: sauvegarde, annulation de g‚p‚alogie, reinitialisation de pro x=prob(j)/ REAL((npop),KIND=doub) y=pro(j)/ REAL((npop),KIND=doub) cum=cum + y sum=sum + y*y !!$ mise en reserve pour l inpression a l+1 pp(l)=prc(j) mm(l)=mrc(j) ir6=ndes(j) xr1=x xr2=y indi(l)=j pr1(l)=x pr2(l)=y prc(j)=0 mrc(j)=0 !!$ initialisation DO i=1,nl IF (pop(i)) THEN pro(i)=1.d0 ELSE pro(i)=0.d0 ENDIF ENDDO !!$ c: calcul des pro !!$ on traite d'abord la population de depart DO k=nl,1,-1 i=ord(k) IF (pop(i)) THEN ip=prc(i) im=mrc(i) IF (ip/=0) pro(ip)=pro(ip) + pro(i)*.5d0 IF (im/=0) pro(im)=pro(im) + pro(i)*.5d0 ENDIF ENDDO !!$ on traite ensuite les autres DO k=nl,1,-1 i=ord(k) IF (.NOT.pop(i)) THEN ip=prc(i) im=mrc(i) IF (ip/=0) pro(ip)=pro(ip) + pro(i)*.5d0 IF (im/=0) pro(im)=pro(im) + pro(i)*.5d0 ENDIF ENDDO !!$ d: ajustements pour les ancetres deja choisis DO i=1,nl aj(i)=0.d0 ENDDO DO i=1,l ii=indi(i) aj(ii)=1.d0 ENDDO DO k=1,nl i=ord(k) ip=prc(i) IF (ip/=0 .AND. aj(ip)/=0.d0) aj(i)=aj(i) + .5d0*aj(ip) im=mrc(i) IF (im/=0 .AND. aj(im)/=0.d0) aj(i)=aj(i) + .5d0*aj(im) ENDDO DO i=1,nl pro(i)=pro(i)*(1.d0 - aj(i)) ENDDO ENDDO !!$ a:recherche de pro maxi pour l+1 (pour le grand minorant) l=l+1 DO i=1,l-1 pro(indi(i))=0.d0 ENDDO r=0.d0 DO i=1,nl IF (.NOT. donact(i) .AND. r