# Pedigrees and Inbreeding 

# Fun1.d - small example data.

library(MASS)
# Retrieve data for the small example
zz=file.choose()
funped <- read.table(file=zz,header=FALSE,
col.names=c("anim","sire","dam"))
funped

# Find out how many unique animals there are

sl = levels(as.factor(funped$sire))
dl = levels(as.factor(funped$dam))
al = levels(as.factor(funped$anim))
cl = c(sl,dl,al)
aid = levels(as.factor(cl))
aid
nanim = length(aid)
nanim

# create new data frame with all animals
k = match(aid,funped$anim)
sid = funped$sire[k]
did = funped$dam[k]
ped = data.frame(aid,sid,did)
ped

# Order animals: Parents before progeny
gen = rep(1,nanim)
for(L in 1:55){
k = 0
 for(i in 1:nanim) {
      kgen = gen[i] + 1;
	ks = sid[i];
	if(!is.na(ks)){
         ksl = match(ks,aid);
	   if(kgen > gen[ksl]){
            gen[ksl] = kgen;
            k = k + 1;
           }
         }
	kd = did[i];
	if(!is.na(kd)){
	   kdl = match(kd,aid);
	   if(kgen > gen[kdl]){
            gen[kdl] = kgen;
            k = k + 1;
        }
     }
    }
if(k == 0 )break;
  }
  
gen
kord = order(-gen)
#  order pedigrees and number animals consecutively
rm(funped)
anim=aid[kord]
sire=sid[kord]
dam=did[kord]
funped=data.frame(anim,sire,dam)
funped   # Parents always come before progeny 

# Number animals consecutively
aid = match(anim,anim)
sid = match(sire,anim)
did = match(dam,anim)
cbind(aid,sid,did)

# Get Bill Szkotnicki's C-functions, in wrappers
source("Bills.R")

# THE SMALL EXAMPLE DATA
xpdinit(nanim) # room for nanim 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

for(i in 1:nanim) { 
   ks = sid[i];
	kd = did[i];
	if(is.na(ks))ks = 0;
      if(is.na(kd))kd = 0;
  inbc[i]=xpdadd(ks,kd);
  inbb[i] = xpdd(i); } 
inbc
inbb
xpdfree()


# Larger data file - animals already numbered consecutively, chronologically
# Except parents without records are not in list
# Read in 'Fun2.d' data file - you must first save it on to the Desktop

zz=file.choose()
dprog <- read.table(file=zz,header=FALSE,
col.names=c("anim","sire","dam","dobirth"))
nrow(dprog)

# Get the birth year of each animal,
#  make it a factor.
bdat = as.integer(dprog$dobirth/10000)
bdatfac = factor(bdat)
dprog[1:10,]
dprog$bdfac = bdatfac
dprog[100:110,]


# Pull out the animal identifications
aid = dprog$anim
sid = dprog$sire
did = dprog$dam
mmm = aid[1]-1
# mmm equals the first animal ID in the data minus 1
#  Animals 1 to mmm were ancestors to those in your data
#  The routines must include these animals, and we
#  must assume that their parents are unknown
aaid = c( 1:mmm, aid)
ssid = c( rep(0,mmm), sid)
ddid = c( rep(0,mmm), did)
ssid[1]
ddid[mmm]
mtotal = length(aaid)
# mtotal = total number of animals in your pedigree
mtotal
mt3=mtotal*3
# initialize space for those animals, and
#  initialize arrays for the inbreeding coefficient and b values
xpdinit(mt3) # room for mtotal
inbc=rep(0,mt3)
inbb=rep(0,mt3)

for(i in 1:mtotal) { inbc[i]=xpdadd(ssid[i],ddid[i]);
  inbb[i] = xpdd(i); } 
inbc[mtotal]
inbb[mtotal]

# Now just grab the inbreeding coefficients for 
# the animals in the data set (ignore animals 1 to mmm)
subf = inbc[aid[1]:mtotal]
subb = length(inbc[inbc > 0])
subb  # Number of animals that are inbred
length(subf)  # should be number of records in 'Fun2.d'

# 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(subf,bdatfac,mean)
aves
plot(aves)
plot(aves, type="b", xaxis = bdatfac)
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,.05),
   xlab="Year of Birth", ylab="Inbreeding Level")
title(main="Change in Inbreeding By Year of Birth")

subb = subf
subb[subb > 0] = 1
nbyr = tapply(subb,bdatfac,sum)
nbyr
plot(nbyr, type="b", lwd=3, col="green",xlab="Year of Birth",
  ylab="Number Inbred")
title(main="Number Inbred By Year of Birth")

# Generate another generation of animals
# Pick 100 sires from aid[1] to aid[1]+1000
# Pick 2000 dams from aid[1]+1001 to mtotal
ks = (sample(1000,100)) + aid[1]
ks
xl = mtotal - aid[1] - 1000
kd = (sample(xl,2000))+aid[1]+1000
kd
# generate one progeny per dam
man = mtotal
sires = sample(ks,2000,replace=TRUE)
for(i in 1:2000){
  jd = kd[i]
  js = sires[i]
  man = man + 1
  inbc[man]=xpdadd(js,jd);
  inbb[man] = xpdd(man); 
}
mtotalp = mtotal + 1
bb=c(aaid,mtotalp:man)
cc=c(ssid,sires)
length(cc)
dd=c(ddid,kd)
length(dd)
aaid = bb
ssid = cc
ddid = dd 

mtotal = man
mtotal

xpdfree()

dyn.unload("rclib.dll")
inbb = inbb[1:mtotal]
length(inbb)

nanim=mtotal
nanimp = nanim+1
bii = 1.0/inbb

# Create the diagonals of A-inverse
adiag = rep(0,nanimp)       # diagonals of A-inverse
bii = c(bii,1)               # bii values for animals
ssid[ssid==0] = nanimp
ddid[ddid==0] = nanimp

#  Set up diagonals of A-inverse
for(i in 1:nanim){
  kis = ssid[i]
  kid = ddid[i]
  d = bii[i]
  x = 0.25*d
  adiag[i] = adiag[i] + d
  adiag[kis] = adiag[kis] + x
  adiag[kid] = adiag[kid] + x
}
adiag[nanimp]=0
adiag

# Set up a Coded pedigree file
# A coded pedigree file is needed for Gibbs Sampling, (not for solving MME)
# But MME will converge faster with coded pedigree file

nped = 0
aaa = rep(0,3*nanim)
sss = aaa
ddd = aaa
cods = aaa
for(i in 1:nanim){
   ka = aaid[i]
   ks = ssid[i]
   kd = ddid[i]
	nped = nped + 1
      cods[nped]=0
      aaa[nped] = ka
      sss[nped] = ks
      ddd[nped] = kd
    if(ks < nanimp){
       nped = nped + 1
       cods[nped]=1
       aaa[nped] = ks
       sss[nped] = ka
       ddd[nped] = kd }
    if(kd < nanimp){
       nped = nped + 1
       cods[nped]=2
       aaa[nped] = kd
       sss[nped] = ka
       ddd[nped] = ks }
  }

aaa = aaa[1:nped]
cods = cods[1:nped]
sss = sss[1:nped]
ddd = ddd[1:nped]
kord = order(aaa)   # Sorting the coded pedigree data - in kord
aaa[kord]