#  Genetic Relationships, Inbreeding

library(MASS)

# new lrscrips.R
zlrs = file.choose()
source(file=zlrs)

# Order pedigrees function





# Enter the data for the small example
#  aid = animal ID (1 to 8)
#  sid = sire IDs - use 0 for unknown sires
#  did = dam IDs  - use 0 for unknown dams

aid = c(1:8)
aid
sid = c(0, 0, 1, 1, 4, 1, 5, 7)
did = c(0, 0, 2, 2, 3, 3, 6, 6)

#  Compute A matrix
numer8 = function(sid,did) {
  N = length(sid)
   A = id(N)
    for(i in 1:N){
      for(j in i:N) {
        ks = sid[j]
        kd = did[j]
        if( ks > 0 && kd > 0 ){
           if(i == j){
               A[i,i] = 1 + 0.5*A[ks,kd] } else
           { A[i,j] = 0.5*(A[i,ks] + A[i,kd])
              A[j,i] = A[i,j] }
         }
        if( ks == 0 && kd > 0 ){
           if(i == j){ A[i,i] = 1 }
            else { A[i,j] = 0.5*A[i,kd] 
                   A[j,i] = A[i,j] }
         }
        if( kd == 0 && ks > 0 ){
           if( i == j){ A[i,i] = 1 } else
              { A[i,j] = 0.5*A[i,ks] 
                 A[j,i] = A[i,j] }
         }
       }
      }
   return(A)
   }
    
A = numer8(sid,did)

Fi = diag(A) - 1
Fi
         
Lp = chol(A)
Lp
L = t(Lp)
L
Bh = diag(Lp)
Bh
B = diag(Bh)
BI = ginv(B)
BI
T = L %*% BI
Tp = BI %*% Lp
T
BB = B %*% B
BB
bi = diag(BB)
bi

peds = data.frame(aid,sid,did,Fi,bi)
peds

TI = ginv(T)
TI = round(TI*10)
TI = TI/10
TI

#  Routines written by Bill Szkotnicki
#  for large pedigrees following the Meuwissen and Luo (1992) algorithm
#  The rclib.dll should be retrieved from course website and stored on DESKTOP.
#  Also, Bills.R and dped.d should be stored on DESKTOP.

zrlib=file.choose()  # open the rclib.dll file
dyn.load(zrlib) # open source file and demo compiled routines
getLoadedDLLs()  # This shows what DLLs are being used

# The "wrapper" functions need to be defined, as given below
#  The routines were written in the C language, and a "wrapper"
#  is necessary for R to understand and use it.
#  You can also make wrappers for FORTRAN programs.

zbill = file.choose()  # open the Bills.R file
source(file=zbill)

# To use these routines ANIMALS MUST BE NUMBERED CONSECUTIVELY
# and they MUST BE IN PARENTS BEFORE PROGENY

# THE SMALL EXAMPLE DATA
nanim = length(aid)
xpdinit(nanim)    # initialize for the number of animals
inbc=rep(0,nanim) # create an array to store the inbreeding coefficients
inbb=rep(0,nanim) # create an array to store the b-values

# Loops are possible in R

for(i in 1:nanim) { 
  inbc[i]=xpdadd(sid[i],did[i])
  inbb[i] = xpdd(i) } 

AI = AINV(sid,did,inbb)
AI
A = solve(AI)
A

xpdfree()

# 1(c) Genetic covariance between full-sibs 
axy=0.5
dxy=0.25
v10 = 1600
v01 = 1000
v11 = 400
v20 = 600
v02 = 200

gcov = axy*v10 + dxy*v01 + axy*dxy*v11 + axy*axy*v20 + dxy*dxy*v02
gcov

#  The assignment asked only for full-sibs, but I compared it to
#    half-sibs in the lab, so some students may show it.
#  covariance between half sibs
axy = 0.25
dxy = 0.0
gcov = axy*v10 + dxy*v01 + axy*dxy*v11 + axy*axy*v20 + dxy*dxy*v02
gcov


#  2) Using R routines on larger data file

# Read in 'dped.d' data file 

zz=file.choose()
peds <- read.table(file=zz,header=FALSE,
col.names=c("anim","sire","dam","birth","sex"))

# Get the birth year of each animal,
#  make it a factor.
bdat = as.integer(peds$birth/10000)
bdatfac = factor(bdat)
nlevels(bdatfac)  # number of birth years in the data


# Pull out the animal identifications
aid = peds$anim
sid = peds$sire
did = peds$dam
nanim = length(aid)
# nanim = total number of animals in your pedigree
nanim

# initialize space for those animals, and
#  initialize arrays for the inbreeding coefficient and b values
xpdinit(nanim) # room for mtotal
inbc=rep(0,nanim)
inbb=rep(0,nanim)

for(i in 1:nanim) { inbc[i]=xpdadd(sid[i],did[i])
  inbb[i] = xpdd(i) } 
inbc[nanim]
inbb[nanim]

# Compute the average inbreeding coefficients by year of birth
#  and plot the results - add titles and labels
#  print the plot and attach to your lab
aves = tapply(inbc,bdatfac,mean)
aves
plot(aves)
plot(aves, type="b", xaxis = bdatfac)  # this statement did not work???
plot(aves, type="b", lwd=3)
plot(aves, type="b", lwd=3, col="blue")
plot(aves, type="b", lwd=3, col="blue", ylim=c(0,.03),
   xlab="Year of Birth", ylab="Inbreeding Level")
title(main="Change in Inbreeding By Year of Birth")

# The students were told to only show the final graph, not all of the
#  intermediate steps

xpdfree()

# Percentage of animals that are inbred
sum(inbc > 0)/nanim

# range of inbreeding coefficients
summary(inbc)
hist(inbc,breaks=100)

# Day 4 Problem

aid = c(1:14)
sid = c(0,0,0,0,1,1,3,3,5,6,6,5,10,11)
did = c(0,0,0,0,2,4,2,4,7,8,7,8,9,8)
anwr=c(5:14)
year=c(1,1,1,1,1,2,2,2,2,2)
cg=c(1,1,2,2,2,3,3,3,4,4)
y=c(16,34,27,23,38,44,25,39,14,26)
nrec=length(y)
nrec
A=numer8(sid,did)
AI=solve(A)
AI = AINV(sid,did,bi) # if bi is known
AI
X=desgn(year,0)
Zc=desgn(cg,0)
Zc
Za=desgn(anwr,0)
Za
Z=cbind(Zc,Za)
Gc=diag(c(6,6,6,6))
Gc
Ga=AI*2
GI=block(Gc,Ga)
RI=id(10)
qqq=MME(X,Z,GI,RI,y)
bhat=qqq$SOLNS
kyr=c(1,2)
kcg=c(3:6)
kan=c(7:20)
length(bhat)
sst=sum(y*y)
ssr=qqq$SSR
vres=as.numeric((sst-ssr)/(nrec-2))
vres
vpe=qqq$VPE
vpe
sep=sqrt(vpe[kan]*vres)
dA=diag(A)
rti=(dA - vpe[kan]*2)/dA
ebvs=cbind(bhat[kan, ],sep,vpe[kan],dA,rti)
ebvs

ebv=bhat[kan]
ebvwr=ebv[anwr]
t(X)%*%ebvwr
yr=factor(year)
tapply(ebvwr,yr,mean)

# Problem 2
bi=c(1,1,1,1,.5,.5,.75,.5,.5)
aid=c(1:9)
sid=c(0,0,0,0,1,3,5,6,3)
did=c(0,0,0,0,2,4,0,7,7)
y=c(15,32,27,8,33,19)
anwr=c(3:8)
X=jd(6,1)
Za=desgn(anwr,9)
A=numer8(sid,did)
A
AI=solve(A)

GI = AI*1.5
RI = id(6)
qqq=MME(X,Za,GI,RI,y)
bhat=qqq$SOLNS
sst=sum(y*y)
ssr=qqq$SSR
vres=as.numeric((sst-ssr)/5)
vpe=qqq$VPE
sep=sqrt(vpe*vres)
cbind(bhat,sep)