WOMBAT – A program for Mixed Model Analyses by Restricted Maximum Likelihood

Differences

This shows you the differences between two versions of the page.

@@ Line 1: / Line 1: @@
+Click tab to download
+<file Fortran Gin4Imprint.f90>
+!------------------------------------------------------------------------------
+!                        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
+</file>

Trace: • referrers • sexlinked • sidebar • wombatauxf90 • wombatcites

Differences

Search

Navigation

Print/export

Tools

QR Code