# NOT RUN {
data(map)
data(Cattle)
dir <- system.file("extdata", package = "optiSel")
files<- file.path(dir, paste("Chr", 1:2, ".phased", sep=""))
### Compute genomic kinship and genomic kinship at native segments
G <- segIBD(files, map, minSNP=20, minL=3.0)
GN <- segIBDatN(files, Cattle, map, thisBreed="Angler", refBreeds="others",
ubFreq=0.02, minSNP=20, minL=3.0, lowMem=TRUE)
Kin <- kinlist(G=G, GN=GN)
### Compute migrant contributions of selection candidates
Haplo<- haplofreq(files, Cattle, map, thisBreed="Angler", refBreeds="others",
minSNP=20, minL=3.0, ubFreq=0.02, what="match")
Comp <- segBreedComp(Haplo$match, map)
Cattle$MC <- NA
Cattle[rownames(Comp), "MC"] <- 1-Comp$native
apply(Comp[,-1],2,mean)
# native F H R
#0.551844104 0.009739393 0.202216271 0.236200232
########################################
# Find optimum breed contributions #
########################################
lb <- c(Angler=0.10, Holstein=0.20, Fleckvieh=0.20)
bc <- opticomp(G, Breed=Cattle$Breed, obj.fun="NGD", lb=lb)$bc
round(bc,3)
# Angler Fleckvieh Holstein Rotbunt
# 0.355 0.445 0.200 0.000
########################################
# Check available objective functions #
# and constraints #
########################################
help.opticont4mb(Kin, Cattle)
##################################################################
# Compute the minimum segment based kinship achievable #
# across breeds while constraining it within the breed #
##################################################################
con <- list(ub.G=0.05, ub=c(M=NA, F=-1))
minG <- opticont4mb("min.G.acrossBreeds", Kin, Cattle, bc, thisBreed="Angler", con=con, trace=FALSE)
minG.s <- summary(minG)
minG.s[,c("G.acrossBreeds")]
#[1] 0.02289039
##################################################################
# Compute the genetic progress achievable while constraining #
# segment based kinship within and across breeds #
# and migrant contributions #
##################################################################
con <- list(ub.G=0.05, ub.G.acrossBreeds=0.026, ub.MC=0.32, ub=c(M=NA, F=-1))
maxBV.G.MC <- opticont4mb("max.BV", Kin, Cattle, bc, thisBreed="Angler", con=con, trace=FALSE)
maxBV.G.MC.s <- summary(maxBV.G.MC)
maxBV.G.MC.s$meanBV
# [1] 0.2851194
##################################################################
# Compute the minimum achievable kinship at native alleles #
# while constraining kinship within and across breeds #
# and migrant contributions #
##################################################################
con <- list(ub.G=0.05, ub.G.acrossBreeds=0.026, ub.MC=0.32, ub=c(M=NA, F=-1))
minGN <- opticont4mb("min.GN", Kin, Cattle, bc, thisBreed="Angler", con=con, solver="slsqp")
minGN.s <- summary(minGN)
minGN.s$GN
#[1] 0.04114953
##################################################################
# Summary statistics from different optimizations #
# can be combined in a data frame. The most important parameters #
# are printed for comparison: #
##################################################################
Res <- rbind(minG.s, maxBV.G.MC.s, minGN.s)
format(Res[,c("valid","meanBV", "meanMC", "G.acrossBreeds", "G", "GN")],digits=4)
# valid meanBV meanMC G.acrossBreeds G GN
#minG TRUE -0.4329 0.3379 0.02289 0.03451 0.03986
#maxBV.G.MC TRUE 0.2851 0.3200 0.02475 0.05000 0.06830
#minGN TRUE -0.5321 0.3200 0.02312 0.03600 0.04115
cor(cbind(minG$parent$oc, maxBV.G.MC$parent$oc, minGN$parent$oc))
# [,1] [,2] [,3]
#[1,] 1.0000000 0.2741736 0.8540842
#[2,] 0.2741736 1.0000000 0.3608900
#[3,] 0.8540842 0.3608900 1.0000000
# }
Run the code above in your browser using DataLab