!------------------------------------------------------------------------------ ! Program GIN4IMPRINT ! ! Auxiliary program for WOMBAT to set up inverse of gametic relationship matrix ! Files written out allow imprinting effects to be fitted with a GIN covariance ! structure: ! 1. The .gin file -> default name: imprint.gin ! 2. The corresponding .codes file -> default name: imprint.codes ! Gametes are coded as: animal ID * 10 + 1 (paternal) and ! animal ID * 10 + 2 (maternal) ! => need to add corresponding column to data file! ! Input required: ! Pedigree file (after pruning) with animal, sire and dam ID recoded in ! running order and animals' inbreeding coefficients; ! this can easily be generated by WOMBAT by running a set-up step for a ! corresponding model without the imprinting effect(s) fitted ! ! Code for gametic relationship matrix written by Bruce Tier !-------------------------------------------------------------km----02/2011----- ! =================== program Gin4Imprint ! =================== implicit none integer, parameter :: ibig = 2147483647 ! Output files for WOMBAT, edit names as desired character(len=13) :: ginfile = 'imprint.gin ', & & codesfile = 'imprint.codes' ! Pruned pedigree file with inbreeding coefficients written out by WOMBAT character(len=18) :: ppedfile = 'PrunedPedFile1.dat' integer :: nanis, ii, ian, idad, imum, jan, jdad, jmum real(8) :: xx write(*,*)'*********** Program "Gin4Imprinting" **************' ! read pedigree file & count no. of animals open(31, file = ppedfile, status = 'old', form ='formatted', action ='read') do read(31, *, iostat=ii ) ian, idad, imum, jan, jdad, jmum, xx if( ii < 0 ) then exit ! end of file else if( ii > 0 ) then write(*,*) 'Some sort of error: reading file "',ppedfile,'"' write(*,*) 'Last running no. read correctly was', ian stop 'Check file!' end if end do write(*,*)'Highest animal ID found =', ian, jan ! set up inverse of gametic relationship matrix nanis = ian; rewind(31) call fbt write(*,*)'"Gametic inverse" complete' write(*,*)'Files written out: "',trim(ginfile),'" and "',trim(codesfile),'"' write(*,*)'******************************************************' stop contains ! ============== subroutine fbt ! ============== real(8), dimension(:), allocatable :: eff, diag, pxp, pxg integer, dimension(:,:), allocatable :: ged, ped real(8) :: rlogl, ras, das, rad, ff, slogl integer :: npxp, i, ia, is, id, ka1, ka2, & & ks1, ks2, kd1, kd2, ngams ! memory requirement = nanis*14 words ngams = nanis*2 allocate( pxg(0:ngams), ped(2,0:nanis), ged(0:2,0:ngams), & eff(0:nanis), diag(0:ngams), pxp(0:ngams), stat=ii ) ped = 0; ged = 0; eff = 0.d0; diag = 0.d0; pxp = 0.d0; pxg = 0.d0 ! read in pedigree & inbreeding and store elements eff(0)= - 1.d0; rlogl = 0.d0 ! read in pedigree and store elements do i = 1, nanis read(31,*) ia, is, id, jan, jmum, jdad, ff xx = jan * 10 + 2 if( xx > dble(ibig) ) then write(*,*)'Original ID too big for coding gametes as "ID*10+i"' write(*,*) jan stop end if ped(1,i) = is ped(2,i) = id eff(i) = ff ! codes for gametes -> original ID*10+i; i=1: paternal, i=2: maternal ka2 = ia * 2; ka1 = ka2 - 1 ged(0,ka1) = 10 * jan + 1 ged(0,ka2) = 10 * jan + 2 ! gametes contributions ras = 2.d0 / (1.d0 - eff(is) ) rad = 2.d0 / (1.d0 - eff(id) ) rlogl = rlogl + log(ras) + log(rad) diag(ka1) = ras diag(ka2) = rad if( is > 0 ) then ! sire gametes ks2 = is * 2 ks1 = ks2 - 1 ged(1,ka1) = ks1 ged(2,ka1) = ks2 ras = -0.5d0*ras pxg(ka1) = ras ras = -0.5d0*ras pxp(ks2) = pxp(ks2) + ras diag(ks1) = diag(ks1) + ras diag(ks2) = diag(ks2) + ras end if if( id > 0 ) then ! dam gametes kd2 = id * 2 kd1 = kd2 - 1 ged(1,ka2) = kd1 ged(2,ka2) = kd2 rad = -0.5d0*rad pxg(ka2) = rad rad = -0.5d0*rad pxp(kd2) = pxp(kd2)+rad diag(kd1) = diag(kd1)+rad diag(kd2) = diag(kd2)+rad end if if( mod(ia,50000) == 0 ) write(*,'(2hAt,3(i12,2i8))') ia, ped(:,ia), ka1, & & ged(1:,ka1), ka2, ged(1:,ka2) end do close(31) ! write out files for WOMBAT open(408, file = codesfile, status = 'unknown', form = 'formatted', & & action = 'write' ) open(409, file = ginfile, status = 'unknown', form = 'formatted', & & action = 'write' ) write(409,*) -rlogl, -rlogl/log(10.0) ! log determinant npxp = 0 do i = 1, ngams write(408,'(i8,i15)') i, ged(0,i) if( pxg(i) /= 0.d0 ) then write(409,'(2i8,f16.8)') ged(1,i), i, pxg(i) write(409,'(2i8,f16.8)') ged(2,i), i, pxg(i) npxp = npxp + 2 end if if( mod(i,2) == 0 .and. pxp(i) /= 0.d0) then write(409,'(2i8,f16.8)') i-1, i, pxp(i) npxp = npxp + 1 end if write(409,'(2i8,f16.8)') i, i, diag(i) end do close(408); close(409) write(*,*)'No. of gametes =', ngams write(*,*)'No. of elements in inverse G =', ngams + npxp return end subroutine fbt end program Gin4Imprint