Nothing
'#
Authors
Torsten Pook, torsten.pook@uni-goettingen.de
Copyright (C) 2017 -- 2020 Torsten Pook
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 3
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
'#
#' Simulation of a breeding program based on a JSON-file from MoBPSweb
#'
#' Function to simulate a breeding program based on a JSON-file from MoBPSweb
#' @param file Path to a json-file generated by the user-interface
#' @param log Provide Path where to write a log-file of your simulation (or false to not write a log-file)
#' @param total Json-file imported via jsonlite::read_json
#' @param fast.mode Set to TRUE work on a small genome with few markers
#' @param progress.bars Set to TRUE to display progress bars
#' @param size.scaling Scale the size of nodes by this factor (especially for testing smaller examples)
#' @param rep.max Maximum number of repeats to use in fast.mode
#' @param verbose Set to FALSE to not display any prints
#' @param miraculix.cores Number of cores used in miraculix applications (default: 1)
#' @param skip.population Set to TRUE to not execute breeding actions (only cost/time estimation will be performed)
#' @param miraculix.chol Set to FALSE to manually deactive the use of miraculix for any cholesky decompostion even though miraculix is actived
#' @param time.check Set to TRUE to automatically check simulation run-time before executing breeding actions
#' @param time.max Maximum length of the simulation in seconds when time.check is active
#' @param export.population Path were to export the population to (at state selected in export.gen/timepoint)
#' @param export.gen Last generation to simulate before exporting population to file
#' @param export.timepoint Last timepoint to simulate before exporting population to file
#' @param fixed.generation.order Vector containing the order of cohorts to generate (Advanced // Testing Parameter!)
#' @examples
#' data(ex_json)
#' \donttest{population <- json.simulation(total=ex_json)}
#' @return Population-list
#' @export
#'
json.simulation <- function(file=NULL, log=NULL, total=NULL, fast.mode=FALSE,
progress.bars=FALSE, size.scaling=NULL, rep.max=1,
verbose=TRUE, miraculix.cores=NULL,
miraculix.chol = NULL,
skip.population = FALSE,
time.check = FALSE,
time.max = 7200,
export.population=FALSE,
export.gen=NULL,
export.timepoint=NULL,
fixed.generation.order=NULL){
if(length(log)==0 && length(file)>0){
log <- paste0(file, ".log")
} else if(length(log)==0 && length(file)==0){
log <- FALSE
}
if(log != FALSE){
zz <- file(log, open="wt")
sink(zz, append = TRUE, type = c("output"), split=TRUE)
cat("#############################################################
############ Modular Breeding Program Simulator #############
#############################################################")
cat("\n")
cat(paste0("Simulation started: ", Sys.time(), "\n"))
cat(paste0("MoBPS version used: ", utils::sessionInfo()$otherPkgs$MoBPS$Version, "\n"))
cat("Copyright (C) 2017-2020 Torsten Pook\n\n")
}
if(requireNamespace("miraculix", quietly = TRUE)){
miraculix <- TRUE
miraculix.dataset <- TRUE
if(length(miraculix.chol)==0){
miraculix.chol <- TRUE
}
} else{
miraculix <- FALSE
miraculix.dataset <- FALSE
if(length(miraculix.chol)==0){
miraculix.chol <- FALSE
}
}
if(length(file)>0){
if(requireNamespace("jsonlite", quietly = TRUE)){
total <- jsonlite::read_json(path=file)
} else{
stop("Use of jsonlite without being installed!")
}
} else if(length(total)==0){
stop("No dataset provided in file or total \n")
}
{
{
if(requireNamespace("RandomFieldsUtils", quietly = TRUE)){
RandomFieldsUtils::RFoptions(helpinfo=FALSE)
if(length(miraculix.cores)>0){
RandomFieldsUtils::RFoptions(cores=miraculix.cores)
} else if( length(total$`Genomic Info`$`number-simulations-parallel`)>0 &&
total$`Genomic Info`$`number-simulations-parallel`!=""){
RandomFieldsUtils::RFoptions(cores=as.numeric(total$`Genomic Info`$`number-simulations-parallel`))
}
}
nodes <- total$Nodes
edges <- total$Edges
geninfo <- total$`Genomic Info`
traitinfo <- total$`Trait Info`
cullinginfo <- total$Culling
subpopulations <- total$Subpopulation$subpopulation_list
litter_size <- total$litter_size
first_pop <- subpopulations[[1]]$Name
if(length(edges)==0){
stop("No edges detected. No valid simulation!")
}
if(length(geninfo$mutation_rate)>0){
mutation.rate <- as.numeric(geninfo$mutation_rate)
} else{
mutation.rate <- 1e-08
}
major_table <- major <- list()
n_traits <- length(traitinfo)
map <- NULL
bp <- NULL
if(length(geninfo$'advanced_miraculix')>0 && geninfo$'advanced_miraculix'==FALSE){
if(verbose) cat("Miraculix has been manually disabled!\n")
miraculix.dataset <- FALSE
miraculix <- FALSE
miraculix.chol <- FALSE
}
if(length(geninfo$`speed-mode`)>0 && geninfo$`speed-mode`=="Yes"){
fast.mode <- TRUE
}
if(length(geninfo$`individual-scaling`)>0 && nchar(geninfo$`individual-scaling`) > 0 && length(size.scaling)==0){
size.scaling <- as.numeric(geninfo$`individual-scaling`)
} else if(length(size.scaling)==0){
size.scaling <- 1
}
fixed_costs <- as.numeric(total$Economy$`Fixed Cost`)
fixed_annual_costs <- 0
geninfo$`Time Unit`
interest_rate <- 1+ as.numeric(total$Economy$`Interest Rate`)/100
if(geninfo$`Time Unit`=="Weeks"){
interest_rate <- interest_rate ^(1/52)
} else if(geninfo$`Time Unit`=="Month"){
interest_rate <- interest_rate ^(1/12)
} else if(geninfo$`Time Unit`=="Days"){
interest_rate <- interest_rate ^(1/365)
}
genotyping_costs <- as.numeric(total$Economy$`Genotyping Cost`)
selection_index_raw <- total$`Selection Index`
selection_index_scaling_raw <- total$`Selection Index Scaling`
if(n_traits>0){
selection_index_name <- numeric(length(selection_index_raw))
selection_index <- matrix(0, nrow=length(selection_index_raw), ncol=n_traits)
selection_index_miesenberger <- rep(FALSE,length(selection_index_raw))
selection_index_miesenberger_wscaling <- rep(FALSE,length(selection_index_raw))
for(index in 1:length(selection_index_raw)){
selection_index[index,] <- as.numeric(unlist(selection_index_raw[[index]])[1:n_traits+1])
selection_index_name[index] <- selection_index_raw[[index]]$'Name'
if(length(selection_index_scaling_raw[[index]]$w_scaling)>0){
if(length(selection_index_scaling_raw[[index]]$miesenberger)>0){
selection_index_miesenberger[index] <- selection_index_scaling_raw[[index]]$miesenberger
}
if(selection_index_scaling_raw[[index]]$w_scaling=="Per Unit"){
selection_index_miesenberger_wscaling[index] <- "unit"
} else if(selection_index_scaling_raw[[index]]$w_scaling=="Per Phenotypic SD"){
selection_index_miesenberger_wscaling[index] <- "pheno_sd"
} else if(selection_index_scaling_raw[[index]]$w_scaling=="Per Breeding Value SD"){
selection_index_miesenberger_wscaling[index] <- "bve_sd"
} else if(selection_index_scaling_raw[[index]]$w_scaling=="Per Genomic Value SD"){
selection_index_miesenberger_wscaling[index] <- "bv_sd"
}
}
}
if(sum(is.na(selection_index))>0){
selection_index[is.na(selection_index)] <- 0
}
} else{
selection_index <- NULL
selection_index_miesenberger_wscaling <- NULL
selection_index_miesenberger <- NULL
selection_index_name <- NULL
}
if(length(litter_size)>0 && length(litter_size$litter_info)>0){
repeat.mating <- matrix(NA, nrow=length(litter_size$litter_info), 2)
for(index in 1:length(litter_size$litter_info)){
repeat.mating[index,1] <- as.numeric(litter_size$litter_info[[index]]$name)
repeat.mating[index,2] <- as.numeric(litter_size$litter_info[[index]]$prob)
}
repeat.mating[is.na(repeat.mating[,1]),1] <- 1
repeat.mating[is.na(repeat.mating[,2]),2] <- 0
} else{
repeat.mating <- cbind(1,1)
}
repeat.mating.copy <- cbind(1,1)
if(length(geninfo$'history_baseline')>0){
base.cycle <- as.numeric(geninfo$'history_baseline')
} else{
base.cycle <- Inf
}
if(length(geninfo$'history_delete')>0){
delete.previous.gen <- geninfo$'history_delete'
} else{
delete.previous.gen <- FALSE
}
pheno_index_raw <- total$`Phenotyping Info`
if(n_traits>0){
pheno_index_name <- pheno_index_costs <- numeric(length(pheno_index_raw))
pheno_index <- matrix(0, nrow=length(pheno_index_raw), ncol=n_traits)
for(index in 1:length(pheno_index_raw)){
pheno_index[index,] <- as.numeric(unlist(pheno_index_raw[[index]])[1:n_traits+2])
pheno_index_costs[index] <- as.numeric(pheno_index_raw[[index]][[2]])
pheno_index_name[index] <- pheno_index_raw[[index]]$'Name'
}
} else{
pheno_index_name <- pheno_index_costs <- numeric(0)
pheno_index <- NULL
}
culling_reason <- matrix(NA, nrow=length(cullinginfo$culling_reasons), ncol=8)
if(length(culling_reason)>0){
for(index in 1:nrow(culling_reason)){
culling_reason[index,] <- unlist(cullinginfo$culling_reasons[[index]])
}
rearrange <- sort(as.numeric(culling_reason[,2]), index.return=TRUE)
culling_reason <- culling_reason[rearrange$ix, ,drop=FALSE]
}
housing_index_raw <- total$Economy$`Animal Housing Costs`
housing_index_name <- housing_index <- numeric(length(housing_index_raw))
if(length(housing_index_raw)>0){
for(index in 1:length(housing_index_raw)){
housing_index[index] <- as.numeric(housing_index_raw[[index]][2])
housing_index_name[index] <- housing_index_raw[[index]]$'Name'
}
housing_cost <- cbind(housing_index_name, housing_index)
}
if(n_traits>0){
trafos <- list()
combi_weights <- list()
trait_matrix <- matrix(0, nrow=n_traits, ncol=18)
for(index in 1:n_traits){
trait_matrix[index,1] <- traitinfo[[index]]$`Trait Name`
trait_matrix[index,2] <- traitinfo[[index]]$`Trait Unit`
trait_matrix[index,3] <- traitinfo[[index]]$`Trait Mean`
trait_matrix[index,4] <- traitinfo[[index]]$`Trait Std Deviation`
trait_matrix[index,5] <- traitinfo[[index]]$`Trait Heritability`
trait_matrix[index,6] <- traitinfo[[index]]$`Trait Number of Polygenic Loci`
trait_matrix[index,7] <- traitinfo[[index]]$`Trait Major QTL`
trait_matrix[index,8] <- traitinfo[[index]]$`Trait Value per Unit`
if(length(traitinfo[[index]]$`dominant_qtl`)==1){
trait_matrix[index,9] <- traitinfo[[index]]$`dominant_qtl`
}
if(length(traitinfo[[index]]$`qualitative_qtl`)==1){
trait_matrix[index,10] <- traitinfo[[index]]$`qualitative_qtl`
}
if(length(traitinfo[[index]]$`quantitative_qtl`)==1){
trait_matrix[index,11] <- traitinfo[[index]]$`quantitative_qtl`
}
if(length(traitinfo[[index]]$'Trait Repeatability')==1 && suppressWarnings(!is.na(as.numeric(traitinfo[[index]]$'Trait Repeatability')))){
trait_matrix[index,12] <- traitinfo[[index]]$'Trait Repeatability'
} else{
trait_matrix[index,12] <- trait_matrix[index,5]
}
if(length(traitinfo[[index]]$is_maternal)==1){
trait_matrix[index,13] <- traitinfo[[index]]$is_maternal
} else{
trait_matrix[index,13] <- FALSE
}
if(length(traitinfo[[index]]$is_paternal)==1){
trait_matrix[index,14] <- traitinfo[[index]]$is_paternal
} else{
trait_matrix[index,14] <- FALSE
}
if(length(traitinfo[[index]]$is_combi)==1){
trait_matrix[index,15] <- traitinfo[[index]]$is_combi
if(traitinfo[[index]]$is_combi){
combi_weights[[index]] <- unlist(as.numeric(traitinfo[[index]]$combi_weights))
}
} else{
trait_matrix[index,15] <- FALSE
}
if(length(traitinfo[[index]]$`additive_equal`)==1){
trait_matrix[index,16] <- traitinfo[[index]]$`additive_equal`
}
if(length(traitinfo[[index]]$`dominant_equal`)==1){
trait_matrix[index,17] <- traitinfo[[index]]$`dominant_equal`
}
if(length(traitinfo[[index]]$`dominant_positive`)==1){
trait_matrix[index,18] <- traitinfo[[index]]$`dominant_positive`
} else{
trait_matrix[index,18] <- FALSE
}
if(length(traitinfo[[index]]$'trafo')>0 && traitinfo[[index]]$'trafo' != "function(x){return(x)}" && nchar(traitinfo[[index]]$'trafo')>11){
f_temp <- NULL
eval(parse(text=paste0("f_temp <- ", traitinfo[[index]]$'trafo')))
trafos[[index]] <- f_temp
}
}
# derive genetic variance based on total variance
pheno_var <- as.numeric(trait_matrix[,4])
trait_matrix[,4] <- as.numeric(trait_matrix[,4]) * sqrt(as.numeric(trait_matrix[,5]))
#traitmean <- as.numeric(c(geninfo[[1]]$MilkYield_Mean, geninfo[[1]]$NonReturnRate_Mean, geninfo[[1]]$SomaticCellScore_Mean))
traitmean <- as.numeric(trait_matrix[,3])
trait_weigths <- as.numeric(trait_matrix[,8])
#heritability <- as.numeric(c(geninfo[[1]]$MilkYield_Hertit, geninfo[[1]]$NonReturnRate_Hertit, geninfo[[1]]$SomaticCellScore_Hertit))
heritability <- as.numeric(trait_matrix[,5])
if(sum(is.na(heritability))>0){
if(verbose) cat("No heritability entered - set heritability to 0.5 for all those traits\n")
heritability[is.na(heritability)]==0
}
repeatability <- as.numeric(trait_matrix[,12])
} else{
heritability <- NULL
}
if(n_traits>0){
cor_gen <- cor_pheno <- matrix(0, nrow=n_traits, ncol= n_traits)
nr <- 1
vector_pheno <- as.numeric(unlist(total$`Phenotypic Correlation`))
vector_gen <- as.numeric(unlist(total$`Genetic Correlation`))
for(index1 in 1:n_traits){
for(index2 in 1:index1){
cor_pheno[index2,index1] <- cor_pheno[index1,index2] <- vector_pheno[nr]
cor_gen[index2,index1] <- cor_gen[index1,index2] <- vector_gen[nr]
nr <- nr + 1
}
}
eigen_gen <- eigen(cor_gen)
if(sum(eigen_gen$values<0)>0){
if(verbose) cat("Genetic covariance matrix is not positive definit.\n")
if(verbose) cat("Generate projection on the set of positive definit matrices:")
test <- eigen_gen
test$values[test$values<0] <- 0
M <- diag(test$values)
S <- test$vectors
newA <- S %*% M %*% solve(S)
diag(newA) <- diag(newA) + 0.0000001 # Avoid numerical issues with inversion
newA <- newA * matrix(1/sqrt(diag(newA)), nrow=nrow(newA), ncol=nrow(newA), byrow=TRUE) * matrix(1/sqrt(diag(newA)), nrow=nrow(newA), ncol=nrow(newA), byrow=FALSE)
if(verbose) cat("new suggested genetic correlation matrix:\n")
cor_gen <- newA
if(verbose) print(round(cor_gen, digits=3))
}
if(length(total$`PhenotypicResidual`)>0 && total$`PhenotypicResidual`){
gen_var <- (diag(trait_matrix[,4])) %*% cor_gen %*% (diag(trait_matrix[,4]))
pheno_var<- (diag(pheno_var)) %*% cor_pheno %*% (diag(pheno_var))
residual_var <- pheno_var - gen_var
cor_pheno <- sqrt(diag(diag(1/residual_var))) %*% residual_var %*% sqrt(diag(diag(1/residual_var)))
cor_pheno[is.na(cor_pheno)] = 0
if(verbose) cat("Used genetic covariance:\n")
if(verbose) print(cor_pheno)
}
eigen_pheno <- eigen(cor_pheno)
if(sum(eigen_pheno$values<0)>0){
if(verbose) cat("Residual correlation matrix is not positive definit.\n")
if(verbose) cat("Generate projection on the set of positive definit matrices:\n")
test <- eigen_pheno
test$values[test$values<0] <- 0
M <- diag(test$values)
S <- test$vectors
newA <- S %*% M %*% solve(S)
diag(newA) <- diag(newA) + 0.0000001 # Avoid numerical issues with inversion
newA <- newA * matrix(1/sqrt(diag(newA)), nrow=nrow(newA), ncol=nrow(newA), byrow=TRUE) * matrix(1/sqrt(diag(newA)), nrow=nrow(newA), ncol=nrow(newA), byrow=FALSE)
if(verbose) cat("new suggested residual correlation matrix:\n")
cor_pheno <- newA
if(verbose) print(round(cor_pheno, digits=3))
}
}
# Creating-types:
# 0 - Founder
# 1 - Selection
# 2 - Reproduction
# 3 - Recombination
# 4 - Selfing
# 5 - DH-Production
# 6 - Cloning
# 7 - Combine
# 8 - Aging
# 9 - Split
#### Character - Numeric trafo ######
for(index in 1:length(nodes)){
nodes[[index]]$'Number of Individuals' <- as.numeric(nodes[[index]]$'Number of Individuals')
nodes[[index]]$'Proportion of Male' <- as.numeric(nodes[[index]]$'Proportion of Male')
}
}
#### MANUEL MODIFICATION:
{
modify_path1 <- sum(dir()=="UserMaps")==0
modify_path2 <- sum(dir()=="UserGenos")==0
if(modify_path1){
if(length(geninfo$`Own Map Path` )>0 && nchar(geninfo$`Own Map Path` )>0){
geninfo$`Own Map Path` <- substr(geninfo$`Own Map Path`, start = gregexpr(pattern = "_" , text = geninfo$`Own Map Path`)[[1]][1] + 1, stop = 1000)
}
}
for(index in 1:length(nodes)){
if(length(nodes[[index]]$delete_info)==0){
nodes[[index]]$delete_info <- FALSE
}
if(modify_path2 && length(nodes[[index]]$Path)>0 && nchar(nodes[[index]]$Path)>0){
nodes[[index]]$Path <- substr(nodes[[index]]$Path, start = gregexpr(pattern = "_" , text = nodes[[index]]$Path)[[1]][1] + 1, stop = 1000)
}
}
for(index in 1:length(edges)){
if(length(edges[[index]]$phenotype_used)==0){
if(length(edges[[index]]$`Use Offspring for BVE`)>0 && edges[[index]]$`Use Offspring for BVE`=="Yes"){
edges[[index]]$phenotype_used <- "Avg. offspring phenotype"
} else{
edges[[index]]$phenotype_used <- "Own phenotype"
}
}
if(edges[[index]]$phenotype_used == "Own phenotype"){
edges[[index]]$`Use Offspring for BVE` <- "No"
} else{
edges[[index]]$`Use Offspring for BVE` <- "Yes"
}
if(length(edges[[index]]$skip)==0){
edges[[index]]$skip <- "No"
}
}
for(index in 1:length(edges)){
if(length(edges[[index]]$'half_sib')==0){
edges[[index]]$'half_sib' <- FALSE
}
if(length(edges[[index]]$'full_sib')==0){
edges[[index]]$'full_sib' <- FALSE
}
if(edges[[index]]$'Breeding Type'=="Semen-collection"){
edges[[index]]$'Breeding Type' <- 'DH-Production'
}
}
if(fast.mode){
if(verbose) cat("Reduce length of genome! Maximum number of Repeat set to ", rep.max, ". I like fast simulations!\n")
geninfo$`Use Ensembl Map` <- "No"
geninfo$`Use Own Map` <- "No"
geninfo$`Chromosomes of Equal Length` <- "Yes"
geninfo$'Chromosomes Info' <- list()
geninfo$'Chromosomes Info'[[1]] <- list()
geninfo$'Chromosomes Info'[[1]]$MD <- 10
geninfo$'Chromosomes Info'[[1]]$Length <- 20
geninfo$'Chromosomes Info'[[1]]$Recombination <- 1
geninfo$'Number of Chromosomes' <- 5
for(index in 1:length(nodes)){
nodes[[index]]$'Path' <- NULL
}
for(index in 1:length(edges)){
if(edges[[index]]$`Breeding Type`=="Repeat"){
edges[[index]]$'Number of Repeat' <- min(as.numeric(edges[[index]]$'Number of Repeat'),rep.max)
}
}
}
for(index in 1:length(edges)){
if(length(edges[[index]]$last_available)==0){
edges[[index]]$last_available <- FALSE
}
}
if(size.scaling!=1){
ids <- numeric(length(nodes))
for(index in 1:length(nodes)){
if(!(length(nodes[[index]]$active_scaling)> 0 && nodes[[index]]$active_scaling)){
nodes[[index]]$`Number of Individuals` <- round(as.numeric(nodes[[index]]$`Number of Individuals`) * size.scaling, digits=0)
}
if( nodes[[index]]$`Number of Individuals`==0){
nodes[[index]]$`Number of Individuals` <- 1
}
ids[index] <- nodes[[index]]$id
}
to_node <- NULL
to_number <- NULL
for(index in 1:length(edges)){
if(edges[[index]]$`Breeding Type`=="Combine"){
to_node <- c(to_node, edges[[index]]$to)
to_number <- c(to_number, nodes[[which(ids==edges[[index]]$from)]]$`Number of Individuals`)
}
}
for(index in unique(to_node)){
take <- which(to_node==index)
if(nodes[[which(ids==index)]]$`Number of Individuals` != sum(to_number[take])){
if(verbose) cat(paste0("Changed number of individuals in ", index, " from ",nodes[[which(ids==index)]]$`Number of Individuals`, " to ",
sum(to_number[take]), " because of Combine rounding in size scaling.\n"))
nodes[[which(ids==index)]]$`Number of Individuals` <- sum(to_number[take])
}
}
to_node <- NULL
to_number <- NULL
for(index in 1:length(edges)){
if(edges[[index]]$`Breeding Type`=="Split"){
to_node <- c(to_node, edges[[index]]$from)
to_number <- c(to_number, nodes[[which(ids==edges[[index]]$to)]]$`Number of Individuals`)
}
}
for(index in unique(to_node)){
take <- which(to_node==index)
if(nodes[[which(ids==index)]]$`Number of Individuals` != sum(to_number[take])){
if(verbose) cat(paste0("Changed number of individuals in ", index, " from ",nodes[[which(ids==index)]]$`Number of Individuals`, " to ",
sum(to_number[take]), " because of Split rounding in size scaling.\n"))
nodes[[which(ids==index)]]$`Number of Individuals` <- sum(to_number[take])
}
}
}
remove.effect.position <- FALSE
for(index in 1:length(nodes)){
if(length(nodes[[index]]$Path)>0 && nchar(nodes[[index]]$Path)>0 && nodes[[index]]$`Genotype generation`!="Upload Genotypes"){
nodes[[index]]$Path <- ""
}
}
}
######################## REALITY - CHECKS ############################
{
for(index in 1:length(edges)){
if(length(edges[[index]]$`Time Needed`)>0){
edges[[index]]$`Time Needed` <- as.numeric(gsub(",", ".", edges[[index]]$`Time Needed`))
}
if((edges[[index]]$`Breeding Type`=="Reproduction" || edges[[index]]$`Breeding Type`=="Cloning" || edges[[index]]$`Breeding Type`=="Selfing" || edges[[index]]$`Breeding Type`=="DH-Production") && (is.na(as.numeric(edges[[index]]$`Time Needed`)) || as.numeric(edges[[index]]$`Time Needed`)==0 )){
edges[[index]]$`Time Needed` <- 0.0000001
if(verbose) cat("Edges to generate new individuals are not supposed to take 0 time. Automatically set to 0.001 time unit.\n")
}
}
for(index in 1:length(nodes)){
if(nodes[[index]]$Sex=="Both" && (nodes[[index]]$`Proportion of Male`==0 ||nodes[[index]]$`Proportion of Male`==1)){
nodes[[index]]$`Proportion of Male` <- 0.5
if(verbose) cat("Both node used with only one sex - set Proportion of male to 0.5. \n")
}
}
for(index in 1:length(nodes)){
if(length(nodes[[index]]$`Genotype generation`)==0 && nodes[[index]]$Founder=="Yes"){
nodes[[index]]$`Genotype generation` <- "Random-sampling"
if(length(nodes[[index]]$Path)>0 && nchar(nodes[[index]]$Path)>0){
nodes[[index]]$`Genotype generation` <- "Upload Genotypes"
}
if(length(nodes[[index]]$`Genotype generation subpopulation`)==0){
nodes[[index]]$`Genotype generation subpopulation` <- first_pop
}
} else if(length(nodes[[index]]$`Genotype generation subpopulation`)==0){
nodes[[index]]$`Genotype generation subpopulation` <- first_pop
}
}
for(index in 1:length(edges)){
if(length(edges[[index]]$'Use Offspring for BVE')==0){
if(verbose) cat("Manually entered Use offspring for BVE \n")
edges[[index]]$'Use Offspring for BVE' <- "No"
}
}
{
for(index in 1:length(nodes)){
if(length(nodes[[index]]$earliest_time)==0){
nodes[[index]]$earliest_time <- 0
}
}
}
{
for(index in 1:length(edges)){
if(length(edges[[index]]$'Time Needed')==0){
if(verbose) cat("Manually entered time_needed\n")
if(edges[[index]]$'Breeding Type'=="Selection" ||edges[[index]]$'Breeding Type'=="Combine" || edges[[index]]$'Breeding Type'=="Split"){
edges[[index]]$'Time Needed' <- 0
} else{
edges[[index]]$'Time Needed' <- 1
}
}
}
}
{
for(index in 1:length(edges)){
if(edges[[index]]$`Breeding Type`=="Selection" && (length(edges[[index]]$`Selection Index`)==0 || sum(selection_index_name==edges[[index]]$`Selection Index`)==0)){
if(verbose) cat(paste0("Invalid Selection Index selected - use ", selection_index_name[1], "\n"))
edges[[index]]$`Selection Index` <- selection_index_name[1]
}
}
}
# Diagonal of pheno/cor-matrix
if(n_traits>0 && sum(diag(cor_gen)==1)!= nrow(cor_gen)){
diag(cor_gen) <- 1
if(verbose) cat("Diagonal of cor-matrix must be 1\n")
}
if(n_traits>0 && sum(diag(cor_pheno)==1)!= nrow(cor_pheno)){
diag(cor_pheno) <- 1
if(verbose) cat("Diagonal of cor-matrix must be 1\n")
}
# Correct nodes are Founders
ids <- possible_founder <- earliest_time <- numeric(length(nodes))
for(index in 1:length(nodes)){
ids[index] <- nodes[[index]]$id
earliest_time[index] <- nodes[[index]]$earliest_time
}
for(index in 1:length(edges)){
possible_founder[which(edges[[index]]$to==ids)] <- 1
}
for(index in which(possible_founder==0)){
if(nodes[[index]]$Founder=="No"){
nodes[[index]]$Founder <- "Yes"
if(verbose) cat(paste0("Changed Note ", nodes[[index]]$id, " to Founder-Note! No incoming edge\n"))
}
}
# All edges
for(index in length(edges):1){
if(sum(edges[[index]]$to==ids)==0 || sum(edges[[index]]$from==ids)==0){
if(verbose) cat("Remove illegal edge. Connected Node not present\n")
edges[[index]] <- NULL
}
}
for(index in 1:length(edges)){
if(edges[[index]]$'Breeding Type'==""){
if(verbose) cat(paste0("Edge ", edges[[index]]$Name, " without breeding type. Invalid breeding program!\n"))
}
}
# Selected cohort is not bigger than founder cohort
for(index in 1:length(edges)){
from <- which(edges[[index]]$from==ids)
to <- which(edges[[index]]$to==ids)
if((nodes[[from]]$'Number of Individuals' < nodes[[to]]$'Number of Individuals') && (edges[[index]]$'Breeding Type'=="Selection")){
if(verbose) cat(paste0("more individuals selected than in founding note in ", edges[[index]]$Name, ". Automatically changed to selected all individuals\n"))
nodes[[to]]$'Number of Individuals' <- nodes[[from]]$'Number of Individuals'
}
}
}
################### Preparation ############################
{
# check for Ind
for(index in 1:length(nodes)){
if(nodes[[index]]$'Sex'=="Indefinit"){
for(index2 in 1:length(nodes)){
nodes[[index2]]$'Sex' <- "Male"
}
}
}
splits <- NULL
n_male <- NULL
for(index in 1:length(nodes)){
if(nodes[[index]]$'Sex'=="Both"){
splits <- c(splits, nodes[[index]]$id)
node_id <- nodes[[index]]$id
nodes[[length(nodes)+1]] <- nodes[[index]]
nodes[[index]]$'Sex' <- "Male"
nodes[[length(nodes)]]$'Sex' <- "Female"
nodes[[index]]$id <- paste0(nodes[[index]]$id, "_M")
nodes[[length(nodes)]]$id <- paste0(nodes[[length(nodes)]]$id, "_F")
nodes[[index]]$label <- paste0(nodes[[index]]$label, "_M")
nodes[[length(nodes)]]$label <- paste0(nodes[[length(nodes)]]$label, "_F")
nodes[[index]]$'Number of Individuals' <- round(as.numeric(nodes[[index]]$'Number of Individuals') * as.numeric(nodes[[index]]$'Proportion of Male'))
n_male <- c(n_male,nodes[[index]]$'Number of Individuals')
nodes[[length(nodes)]]$'Number of Individuals' <- as.numeric(nodes[[length(nodes)]]$'Number of Individuals') - nodes[[index]]$'Number of Individuals'
nodes[[index]]$'Proportion of Male' <- 1
nodes[[length(nodes)]]$'Proportion of Male' <- 0
for(index2 in 1:length(edges)){
if(length(edges[[index2]]$`Manuel selected cohorts`)>0){
for(index3 in 1:length(edges[[index2]]$`Manuel selected cohorts`)){
split_up <- strsplit(edges[[index2]]$`Manuel selected cohorts`[[index3]],":")[[1]]
if(split_up[1]==node_id){
edges[[index2]]$`Manuel selected cohorts`[[index3]] <- paste0(nodes[[index]]$id,":", split_up[2])
edges[[index2]]$`Manuel selected cohorts`[[length(edges[[index2]]$`Manuel selected cohorts`)+1]] <- paste0(nodes[[length(nodes)]]$id,":", split_up[2])
}
}
}
}
}
}
for(index in 1:length(edges)){
wfrom <- which(edges[[index]]$from==splits)
wto <- which(edges[[index]]$to==splits)
if(length(wfrom)>0 && length(wto)>0){
edges[[length(edges)+1]] <- edges[[index]]
edges[[index]]$to <- paste0(edges[[index]]$to, "_M")
edges[[length(edges)]]$to <- paste0(edges[[length(edges)]]$to, "_F")
edges[[index]]$from <- paste0(edges[[index]]$from, "_M")
edges[[length(edges)]]$from <- paste0(edges[[length(edges)]]$from, "_F")
} else if(length(wto)>0){
edges[[length(edges)+1]] <- edges[[index]]
edges[[index]]$to <- paste0(edges[[index]]$to, "_M")
edges[[length(edges)]]$to <- paste0(edges[[length(edges)]]$to, "_F")
} else if(length(wfrom)>0){
edges[[length(edges)+1]] <- edges[[index]]
edges[[index]]$from <- paste0(edges[[index]]$from, "_M")
edges[[length(edges)]]$from <- paste0(edges[[length(edges)]]$from, "_F")
}
}
ids <- earliest_time <- numeric(length(nodes))
sex <- numeric(length(nodes))
for(index in 1:length(nodes)){
earliest_time[index] <- nodes[[index]]$earliest_time
ids[index] <- nodes[[index]]$id
sex[index] <- as.numeric(nodes[[index]]$'Sex'=="Female")
}
# Remove sex-switch
for(index in length(edges):1){
sex1 <- sex[which(edges[[index]]$from==ids)]
sex2 <- sex[which(edges[[index]]$to==ids)]
if(sex1!=sex2 && !(edges[[index]]$'Breeding Type'=="Reproduction" || edges[[index]]$'Breeding Type'=="Selfing" ||
edges[[index]]$'Breeding Type'=="Cloning" || edges[[index]]$'Breeding Type'=="DH-Production")){
edges[[index]] <- NULL
}
}
# Identify founder cohorts
# Add "Founder" as a breeding type - If selected allow for a dataset to be imported!
founder <- NULL
ids <- earliest_time <- numeric(length(nodes))
repeat_node <- numeric(length(nodes))
incoming_repeat_node <- numeric(length(nodes))
for(index in 1:length(nodes)){
if(nodes[[index]]$Founder=="Yes"){
founder <- c(founder,index)
}
ids[index] <- nodes[[index]]$id
earliest_time[index] <- nodes[[index]]$earliest_time
}
for(index in 1:length(edges)){
if(edges[[index]]$'Breeding Type'=="Repeat"){
repeat_node[which(edges[[index]]$from==ids)] <- 1
incoming_repeat_node[which(edges[[index]]$to==ids)] <- 1
}
}
}
## REPEATS
{
repeat_links <- NULL
for(index in 1:length(edges)){
if(edges[[index]]$`Breeding Type`=="Repeat"){
repeat_links <- rbind(repeat_links, c(edges[[index]]$from, edges[[index]]$to))
}
}
while(sum(repeat_node)>0){
to_early_to_repeat <- 0
while(length(to_early_to_repeat)>0){
required_nodes <- list()
required_nodes[[length(ids)+1]] <- "placeholder"
for(index in 1:length(edges)){
if(edges[[index]]$`Breeding Type`!="Repeat"){
new_node <- which(edges[[index]]$to==ids)
required_nodes[[new_node]] <- c(required_nodes[[new_node]], which(edges[[index]]$from==ids))
}
}
n_nodes <- length(nodes)
n_edges <- length(edges)
link <- NULL
link2 <- NULL
if(length(to_early_to_repeat)==1 || to_early_to_repeat==0){
start <- founder
} else{
start <- intersect(founder, (1:length(ids))[-to_early_to_repeat])
}
step <- start
step_vali <- intersect(start, which(incoming_repeat_node==1))
nodes_to_repeat <- NULL
edges_to_repeat <- NULL
prev <- NULL
while(length(step)!=length(prev) || length(prev)==0){
prev <- step
for(index2 in 1:length(edges)){
desti <- which(edges[[index2]]$to==ids)
ori <- which(edges[[index2]]$from==ids)
if(length(intersect(ori, step))>0 & length(required_nodes[[desti]])== length(intersect(required_nodes[[desti]], step))){
if(repeat_node[ori]==1){
if(edges[[index2]]$`Breeding Type`=="Repeat"){
start <- unique(c(start, desti))
step <- unique(c(step, desti))
if(incoming_repeat_node[desti] || length(intersect(step_vali,ori))>0){
step_vali <- unique(c(step_vali, desti))
}
link <- rbind(link, c(ori, desti))
link2 <- unique(c(link2, index2))
}
} else{
if(sum(ori==step)>0){
step <- unique(c(step, desti))
if(incoming_repeat_node[desti]==1 || length(intersect(step_vali,ori))>0){
step_vali <- unique(c(step_vali, desti))
}
}
}
}
}
}
nodes_to_repeat <- step_vali[!duplicated(c(start, step_vali))[-(1:length(start))]]
incoming_repeat_node2 <- numeric(length(incoming_repeat_node))
for(index in 1:length(edges)){
if(edges[[index]]$'Breeding Type'=="Repeat"){
temp2 <- which(edges[[index]]$from==ids)
if(length(intersect(temp2, prev))>0){
temp1 <- which(edges[[index]]$to==ids)
incoming_repeat_node2[temp1] <- 1
}
}
}
# check if all founder nodes for repeat nodes are available
to_early_to_repeat_temp <- intersect(intersect(ids[prev],ids[incoming_repeat_node2==1]), ids[-link[,2]])
if(length(to_early_to_repeat_temp)>0){
for(index in 1:length(to_early_to_repeat_temp)){
to_early_to_repeat_temp[index] <- which(ids==to_early_to_repeat_temp[index])
}
to_early_to_repeat <- as.numeric(c(to_early_to_repeat_temp, to_early_to_repeat))
} else{
to_early_to_repeat <- NULL
}
}
nodes_to_repeat <- sort(nodes_to_repeat)
n_rep <- 0
for(index in 1:length(edges)){
if(sum(edges[[index]]$to==ids[nodes_to_repeat])>0 && edges[[index]]$`Breeding Type`!="Repeat"){
edges_to_repeat <- c(edges_to_repeat, index)
}
if(sum(edges[[index]]$from==ids[nodes_to_repeat]) && edges[[index]]$`Breeding Type`=="Repeat"){
n_rep <- as.numeric(edges[[index]]$`Number of Repeat`)
}
}
# edges_to_repeat <- sort(edges_to_repeat)
link <- unique(link)
# for(index4 in start){
# edges[[index4]]$'Breeding Type' <- "PerformedRepeat"
# }
for(index4 in 1:n_rep){
for(nodrep in nodes_to_repeat){
nodes[[length(nodes)+1]] <- nodes[[nodrep]]
nodes[[length(nodes)]]$Founder <- "No"
nodes[[length(nodes)]]$id <- paste0(nodes[[length(nodes)]]$id, "_",index4)
nodes[[length(nodes)]]$label <- paste0(nodes[[length(nodes)]]$label, "_",index4)
}
for(edgesrep in edges_to_repeat){
edges[[length(edges)+1]] <- edges[[edgesrep]]
test <- ids[start]==edges[[length(edges)]]$from
if(sum(edges[[edgesrep]]$from==ids[step_vali])==0){
if(verbose) cat(paste0("Repeat from edge not included in the repeat: ", edges[[edgesrep]]$from, "\n"))
} else if(sum(test)){
take <- link[which(link[,2]==start[which(test)]),1]
if(index4==1){
edges[[length(edges)]]$from <- paste0(ids[take])
} else{
edges[[length(edges)]]$from <- paste0(ids[take],"_", index4-1)
}
} else{
edges[[length(edges)]]$from <- paste0(edges[[length(edges)]]$from,"_", index4)
}
edges[[length(edges)]]$to <- paste0(edges[[length(edges)]]$to,"_", index4)
edges[[length(edges)]]$id <- paste0(edges[[length(edges)]]$id,"_", index4)
}
}
for(sonstcheck in (1:length(ids))[-step_vali]){
for(edgeindex in 1:length(edges)){
if(edges[[edgeindex]]$to==ids[sonstcheck] && sum(edges[[edgeindex]]$from==ids[nodes_to_repeat])>0){
edges[[edgeindex]]$from <- paste0(edges[[edgeindex]]$from, "_", n_rep)
}
}
}
for(index in 1:length(edges)){
if(edges[[index]]$`Breeding Type`=="Repeat" && sum(ids[step_vali]==edges[[index]]$to)>0){
edges[[index]]$to <- paste0(edges[[index]]$to, "_", n_rep)
}
}
for(index in 1:nrow(link)){
if(verbose) cat(paste0("Successful repeat for ", ids[link[index,1]] ," to ",ids[link[index,2]],"\n" ))
}
if(verbose) cat(paste0(length(nodes)-n_nodes, " Nodes & ", length(edges) - n_edges, " Edges were added to generation script.\n"))
ids <- earliest_time <- numeric(length(nodes))
founder <- NULL
for(index5 in 1:length(nodes)){
if(nodes[[index5]]$Founder=="Yes"){
founder <- c(founder,index5)
}
ids[index5] <- nodes[[index5]]$id
earliest_time[index5] <- nodes[[index5]]$earliest_time
}
for(changeto in sort(link2, decreasing = TRUE)){
edges[[changeto]] <- NULL
}
repeat_node <- incoming_repeat_node <- numeric(length(nodes))
for(index in 1:length(edges)){
if(edges[[index]]$'Breeding Type'=="Repeat"){
repeat_node[which(edges[[index]]$from==ids)] <- 1
incoming_repeat_node[which(edges[[index]]$to==ids)] <- 1
}
}
}
## Change manuel selected according to repeats
for(index in 1:length(edges)){
if(length(edges[[index]]$`Manuel selected cohorts`)>0){
for(index2 in 1:length(edges[[index]]$`Manuel selected cohorts`)){
split_up <- strsplit(edges[[index]]$`Manuel selected cohorts`[[index2]], ":")[[1]]
repeat_change <- which(split_up[1]==repeat_links[,2])
if(length(repeat_change)==1 && split_up[2] != "X"){
edges[[index]]$`Manuel selected cohorts`[[index2]] <- paste0(repeat_links[repeat_change,1], ":", as.numeric(split_up[2])-1)
}
}
}
}
## Determine genetic architecture
if(geninfo$'Use Ensembl Map'=="Yes"){
if (requireNamespace("MoBPSmaps", quietly = TRUE)) {
if(geninfo$`Ensembl Dataset`=="Axiom Genotyping Array"){
map <- MoBPSmaps::map_pig1
} else if(geninfo$`Ensembl Dataset`=="GGP Porcine HD"){
map <- MoBPSmaps::map_pig2
} else if(geninfo$`Ensembl Dataset`=="GGP Porcine LD"){
map <- MoBPSmaps::map_pig3
} else if(geninfo$`Ensembl Dataset`=="Illumina_PorcineSNP60"){
map <- MoBPSmaps::map_pig4
} else if(geninfo$`Ensembl Dataset`=="Affymetrix Chicken600K Array"){
map <- MoBPSmaps::map_chicken1
} else if(geninfo$`Ensembl Dataset`=="Affymetrix Chicken600K Array_subset_diversity"){
map <- MoBPSmaps::map_chicken2
} else if(geninfo$`Ensembl Dataset`=="Affymetrix Chicken600K Array_subset_50k"){
map <- MoBPSmaps::map_chicken3
} else if(geninfo$`Ensembl Dataset`=="Illumina BovineSNP50 BeadChip"){
map <- MoBPSmaps::map_cattle1
} else if(geninfo$`Ensembl Dataset`=="Illumina BovineHD BeadChip"){
map <- MoBPSmaps::map_cattle2
} else if(geninfo$`Ensembl Dataset`=="Illumina BovineLD BeadChip"){
map <- MoBPSmaps::map_cattle3
}else if(geninfo$`Ensembl Dataset`=="Genotyping chip variations"){
map <- MoBPSmaps::map_cattle4
}else if(geninfo$`Ensembl Dataset`=="Illumina EquineSNP50 BeadChip"){
map <- MoBPSmaps::map_horse1
}else if(geninfo$`Ensembl Dataset`=="IlluminaOvineHDSNP"){
map <- MoBPSmaps::map_sheep1
}else if(geninfo$`Ensembl Dataset`=="IlluminaOvineSNP50"){
map <- MoBPSmaps::map_sheep2
} else if(geninfo$`Ensembl Dataset`=="Genotyping chip variants"){
map <- MoBPSmaps::map_sheep3
} else if(geninfo$`Ensembl Dataset`=="Illumina_GoatSNP50"){
map <- MoBPSmaps::map_goat1
} else if(geninfo$`Ensembl Dataset`=="Affy GeneChip500K"){
map <- MoBPSmaps::map_human1
} else if(geninfo$`Ensembl Dataset`=="Illumina_1M-duo"){
map <- MoBPSmaps::map_human2
} else if(geninfo$`Ensembl Dataset`== "Illumina_HumanHap550"){
map <- MoBPSmaps::map_human3
} else if(geninfo$`Ensembl Dataset`== "Affymetrix Axiom Maize Genotyping Array"){
map <- MoBPSmaps::map_maize1
} else if(geninfo$`Ensembl Dataset`== "3k_Sorghum_Bekele"){
map <- MoBPSmaps::map_sorghum1
} else if(geninfo$`Ensembl Dataset`== "12k_Wheat_Bekele"){
map <- MoBPSmaps::map_wheat1
} else if(geninfo$`Ensembl Dataset`== "29k_Wheat_Wen"){
map <- MoBPSmaps::map_wheat2
} else if(geninfo$`Ensembl Dataset`== "96k_AtlanticSalmon_male_Tsai"){
map <- MoBPSmaps::map_salmon1
} else if(geninfo$`Ensembl Dataset`== "96k_AtlanticSalmon_female_Tsai"){
map <- MoBPSmaps::map_salmon2
} else if(geninfo$`Ensembl Dataset`== "68k_NileTilapia_Penaloza"){
map <- MoBPSmaps::map_salmon2
} else{
map <- MoBPSmaps::ensembl.map(dataset = geninfo$'Ensembl Dataset',
filter = geninfo$'Ensembl Filter',
filter.values = geninfo$'Ensembl Filter Values')
}
} else{
stop("Use of Ensembl-Maps without MoBPSmaps R-package.
## To Install:
## devtools::install_github('tpook92/MoBPS', subdir='pkg-maps')
## Or download from https://github.com/tpook92/MoBPS")
}
if(length(geninfo$'Max Number of SNPs')>0 && as.numeric(geninfo$'Max Number of SNPs')<nrow(map)){
map <- map[sort(sample(1:nrow(map), as.numeric(geninfo$'Max Number of SNPs'))),]
}
chr.opt <- unique(map[,1])
nchromo <- length(chr.opt)
nsnp <- chromo.length <- numeric(nchromo)
for(index in 1:nchromo){
nsnp[index] <- sum(map[,1]==chr.opt[index])
chromo.length[index] <- max(as.numeric(map[map[,1]==chr.opt[index],4])) / 100000000
}
if(verbose) cat("Assume 100.000.000 bp/M in Ensembl Map \n")
} else if(geninfo$'Use Own Map'=="Yes" || length(map)>0){
map_path <- geninfo$`Own Map Path`
if(length(map_path)>0){
map_type <- substr(map_path, start= nchar(map_path)-2, stop= nchar(map_path))
if(map_type=="vcf"){
if(verbose) cat("Map identified as vcf-file - extract map information")
if(requireNamespace("vcfR", quietly = TRUE)){
vcf_file <- vcfR::read.vcfR(map_path)
map <- cbind(vcf_file@fix[,c(1,3)],NA,vcf_file@fix[,2], NA)
} else{
stop("Use of vcfR without being installed!")
}
} else if(map_type=="map"){
if(verbose) cat("Map identified as Ped-map-file - extract map information")
map_file <- utils::read.table(map_path)
map <- cbind(map_file[,c(1,2,3,4)],NA)
if(sum(is.na(map[,3]))>0 || sum(map[,3]==0)>0){
map[,3] <- NA
}
} else if(map_type=="ata"){
map_store <- load(map_path)
if(length(map_store)>1){
if(verbose) cat("More than one object contain in .Rdata - unless the map is named 'map' import will fail! \n")
} else{
map <- eval(parse(text=map_store[[1]]))
}
}
map <- as.matrix(map)
}
chr.opt <- unique(map[,1])
nchromo <- length(chr.opt)
nsnp <- chromo.length <- numeric(nchromo)
for(index in 1:nchromo){
nsnp[index] <- sum(map[,1]==chr.opt[index])
chromo.length[index] <- max(as.numeric(map[map[,1]==chr.opt[index],4])) / 100000000
}
if(verbose) cat("Assume 100.000.000 bp/M in Imported Map \n")
} else if(geninfo$'Chromosomes of Equal Length'=="Yes"){
nchromo <- as.numeric(geninfo$'Number of Chromosomes')
nsnp <- ceiling(rep(as.numeric(geninfo$'Chromosomes Info'[[1]]$Length) * as.numeric(geninfo$'Chromosomes Info'[[1]]$MD), nchromo))
chromo.length <- rep(as.numeric(geninfo$'Chromosomes Info'[[1]]$Length) * as.numeric(geninfo$'Chromosomes Info'[[1]]$Recombination, nchromo) /100, nchromo)
bp <- ceiling((0.5:nsnp[1]) / nsnp[1] * as.numeric(geninfo$'Chromosomes Info'[[1]]$Length) * 1000000)
} else{
nchromo <- as.numeric(geninfo$'Number of Chromosomes')
chromo.length <- nsnp <- numeric(nchromo)
bp <- numeric(sum(nsnp))
prev <- 0
for(index in 1:nchromo){
nsnp[index] <- ceiling(as.numeric(geninfo$'Chromosomes Info'[[index]]$Length) * as.numeric(geninfo$'Chromosomes Info'[[index]]$MD))
chromo.length[index] <- as.numeric(geninfo$'Chromosomes Info'[[index]]$Length) * as.numeric(geninfo$'Chromosomes Info'[[index]]$Recombination, nchromo) /100
bp[(1:nsnp[index])+prev] <- ceiling((0.5:nsnp[index]) / nsnp[index] * as.numeric(geninfo$'Chromosomes Info'[[index]]$Length) * 1000000)
prev <- prev + nsnp[index]
}
}
if(length(map)==0){
for(index in 1:nchromo){
map <- rbind(map, cbind(index, paste0("Chr", index,"SNP", 1:nsnp[index]), NA, NA, NA))
}
if(length(bp)>0){
map[,4] <- bp
} else{
if(verbose) cat("Why is there no base-pair in auto generated map?")
}
}
groups <- length(nodes)
# Species not needed, traitvariance not needed,
position <- matrix(0, nrow=length(nodes), ncol=4)
rownames(position) <- ids
founding_a <- c(0,0)
mig_m <- numeric(0)
mig_f <- numeric(0)
sex.s <- NULL
genotyped.s <- NULL
# CHECK FOR ME THAN 2 FOUNDERS (MORE THAN 1 of a sex- Migration level...)
new_mig <- c(0,0)
founder_data <- FALSE
for(index in 1:length(founder)){
sex <- as.numeric(nodes[[founder[index]]]$'Sex'=="Female")+1
size <- as.numeric(nodes[[founder[index]]]$'Number of Individuals')
founding_a[sex] <- founding_a[sex] + size
sex.s <- c(sex.s, rep(sex, size))
genotyped.s <- c(genotyped.s, stats::rbinom(as.numeric(nodes[[founder[index]]]$'Number of Individuals'), 1,
as.numeric(nodes[[founder[index]]]$`Proportion of genotyped individuals`)))
if(sex==1){
mig_m <- c(mig_m, rep(new_mig[1], size))
mig <- new_mig[1]
new_mig[1] <- new_mig[1] + 1
}
if(sex==2){
mig_f <- c(mig_f, rep(new_mig[2], size))
mig <- new_mig[2]
new_mig[2] <- new_mig[2] + 1
}
position[founder[index],] <- c(1, sex, mig, size)
if(nodes[[index]]$`Genotype generation` == "Upload Genotypes" && length(nodes[[founder[index]]]$Path)>0 && nchar(nodes[[founder[index]]]$Path)>0){
founder_data <- TRUE
}
}
}
################### Import genetic data + major QTL import ####################
{
founder_pedigree <- NULL
if(!skip.population){
path_list <- NULL
phasing_list <- NULL
no_data <- numeric(length(founder))
founder_pop <- numeric(length(founder))
for(index in 1:length(founder)){
path_list <- c(path_list, nodes[[founder[index]]]$Path )
if(length(nodes[[founder[index]]]$'phasing_required') >0 && nodes[[founder[index]]]$'phasing_required'=="Yes"){
phasing_list <- c(phasing_list, nodes[[founder[index]]]$Path )
}
founder_pop[index] <- nodes[[founder[index]]]$`Genotype generation subpopulation`
}
p_i <- list()
major_sub <- list()
majorp <- NULL
major_table_sub <- list()
subpopulation_info <- matrix(0, nrow=length(subpopulations), ncol=10)
for(index in 1:length(subpopulations)){
subpopulation_info[index,] <- unlist(subpopulations[[index]])[1:10]
}
if( (length(subpopulations)>1 || sum(as.numeric(subpopulation_info[,9])>0) > 0) & miraculix.dataset){
miraculix.dataset <- FALSE
if(verbose) cat("miraculix dataset generation only for single subpopulation and no LD build up.\n")
}
if(length(unique(c("",path_list)))>1 && miraculix.dataset){
miraculix.dataset <- FALSE
if(verbose) cat("Detected data-path for founders. Deactive miraculix.dataset.\n")
}
if(!miraculix.dataset){
dataset <- matrix(0L, nrow= sum(nsnp), ncol = sum(position[,4])*2)
}
starts <- cumsum(c(1, position[founder,4]))*2-1
for(path in unique(c("",path_list))[-1]){
data_path <- path
if(length(data_path)>0){
data_type <- substr(data_path, start= nchar(data_path)-2, stop= nchar(data_path))
if(data_type=="vcf"|| data_type==".gz"){
if(verbose) cat("Data input identified as vcf-file - extract genomic information \n")
if(requireNamespace("vcfR", quietly = TRUE)){
if(sum(path==phasing_list)>0){
pedmap.to.phasedbeaglevcf(vcf_path = data_path)
data_path <- "/home/nha/Plink/DB/temp_vcf_phased.vcf.gz"
}
vcf_file <- vcfR::read.vcfR(data_path)
haplo1 <- substr(vcf_file@gt[,-1], start=1, stop=1)
haplo2 <- substr(vcf_file@gt[,-1], start=3, stop=3)
storage.mode(haplo1) <- storage.mode(haplo2) <- "integer"
haplo_temp <- cbind(haplo1, haplo2)
haplo_temp <- haplo_temp[,c(0,ncol(haplo1)) + rep(1:ncol(haplo1), each=2)]
order <- numeric(nrow(vcf_file@fix))
prev <- 0
max_print1 <- TRUE
max_print2 <- TRUE
it <- nrow(map)
skips <- numeric(nrow(vcf_file@fix))
for(index in 1:nrow(vcf_file@fix)){
if(prev < it && map[prev+1,2]==vcf_file@fix[index,3]){
prev <- prev +1
} else{
prev <- which(map[,2]==vcf_file@fix[index,3])
}
if(length(prev)==1){
order[index] <- prev
} else if(length(prev)>1){
skips[index] <- index
if(max_print1){
warning("Marker names in vcf are not unique!")
max_print1 <- FALSE
}
} else if(length(prev)==0 ){
skips[index] <- index
if(max_print2){
warning("Marker in VCF that do not appear in the map!\n")
max_print2 <- FALSE
}
}
if(length(prev)!=1){
prev <- 0
}
}
haplo <- matrix(0,nrow=nrow(map), ncol=ncol(haplo_temp))
if(sum(skips)>0){
haplo[order,] <- haplo_temp[-skips,]
} else{
haplo[order,] <- haplo_temp
}
if(sum(is.na(haplo))>0){
warning("Your dataset contains missing elements! \nEverything missing is set to 0. Consider Phasing!\n")
if(verbose) cat(paste0(sum(is.na(haplo)), " missing entries in haplotype data!\n"))
haplo[is.na(haplo)] <- 0
}
} else{
stop("Use of vcfR without being installed!")
}
} else if(data_type=="ped"){
if(verbose) cat("Data input identified as Ped-map-file - extract genomic information \n")
if(sum(path==phasing_list)>0){
pedmap.to.phasedbeaglevcf(ped_path=path, map_path=geninfo$`Own Map Path`)
if(requireNamespace("vcfR", quietly = TRUE)){
vcf_file <- vcfR::read.vcfR("/home/nha/Plink/DB/temp_vcf_phased.vcf.gz")
haplo1 <- substr(vcf_file@gt[,-1], start=1, stop=1)
haplo2 <- substr(vcf_file@gt[,-1], start=3, stop=3)
storage.mode(haplo1) <- storage.mode(haplo2) <- "integer"
haplo_temp <- cbind(haplo1, haplo2)
haplo_temp <- haplo_temp[,c(0,ncol(haplo1)) + rep(1:ncol(haplo1), each=2)]
order <- numeric(nrow(vcf_file@fix))
prev <- 0
max_print1 <- TRUE
max_print2 <- TRUE
it <- nrow(map)
skips <- numeric(nrow(vcf_file@fix))
for(index in 1:nrow(vcf_file@fix)){
if(prev < it && map[prev+1,2]==vcf_file@fix[index,3]){
prev <- prev +1
} else{
prev <- which(map[,2]==vcf_file@fix[index,3])
}
if(length(prev)==1){
order[index] <- prev
} else if(length(prev)>1){
skips[index] <- index
if(max_print1){
warning("Marker names in vcf are not unique!")
max_print1 <- FALSE
}
} else if(length(prev)==0 ){
skips[index] <- index
if(max_print2){
warning("Marker in VCF that do not appear in the map!\n")
max_print2 <- FALSE
}
}
if(length(prev)!=1){
prev <- 0
}
}
haplo <- matrix(0,nrow=nrow(map), ncol=ncol(haplo_temp))
if(sum(skips)>0){
haplo[order,] <- haplo_temp[-skips,]
} else{
haplo[order,] <- haplo_temp
}
if(sum(is.na(haplo))>0){
warning("Your dataset contains missing elements.\n Everything missing is set to 0. Consider Phasing!")
if(verbose) cat(paste0(sum(is.na(haplo)), " missing entries in haplotype data!\n"))
haplo[is.na(haplo)] <- 0
}
} else{
stop("Use of vcfR without being installed!")
}
} else{
if(verbose) cat("Haplotype phase is assumed to be by colum - No internal phasing performed! \n")
ped_file <- utils::read.table(data_path)
haplo12 <- t(ped_file[,-(1:6)])
haplo <- matrix(0, ncol = ncol(haplo12)*2, nrow=nrow(haplo12)/2)
for(index1 in 1:ncol(haplo12)){
haplo[,index1*2+c(-1,0)] <- matrix(haplo12[,index1], ncol=2, byrow=TRUE)
}
}
} else if(data_type=="ata"){
Map <- map ## Lisas subpos contain a element named map"
data_store <- load(data_path)
map <- Map
if(length(data_store)>1){
if(verbose) cat("More than one object contain in .Rdata - unless the data object is named 'haplo' import will fail! \n")
} else{
haplo <- eval(parse(text=map_store[[1]]))
}
}
}
takes <- which(path_list==path)
no_data[takes] <- 1
take <- NULL
for(index in takes){
take <- c(take, starts[index]:(starts[index+1]-1))
}
if(nrow(haplo)!=sum(nsnp)){
if(data_type=="vcf" || data_type==".gz"){
chr.opt <- unique(vcf_file@fix[,1])
nsnp <- numeric(length(chr.opt))
for(index in 1:length(chr.opt)){
nsnp[index] <- sum(vcf_file@fix[,1]==chr.opt[index])
}
} else{
nsnp <- floor(nsnp*nrow(haplo)/sum(nsnp))
nsnp[1] <- nsnp[1] + nrow(haplo)-sum(nsnp)
}
dataset <- matrix(0L, nrow= sum(nsnp), ncol = sum(position[,4])*2)
map <- NULL
for(index in 1:nchromo){
map <- rbind(map, cbind(index, paste0("SNP", 1:nsnp[index]), NA, NA, NA))
}
}
if(length(haplo) != length(dataset[,take])){
if(verbose) cat("Size of Founder-dataset not in concorance with node size!\n")
haplo <- matrix(rep(haplo, length.out=length(dataset[,take])), ncol=length(take))
}
if(storage.mode(haplo)!="integer"){
storage.mode(haplo) <- "integer"
}
dataset[,take] <- haplo
}
if((n_traits)>0){
for(index in 1:n_traits){
major[[index]] <- traitinfo[[index]]$'Trait QTL Info'
if(length(major[[index]])>0 && length(major[[index]][[1]])>1){
ncols <- length(major[[index]][[1]])
mqtl <- matrix(unlist(major[[index]]), ncol=ncols, byrow=TRUE)
to_enter <- mqtl[,c(1,4:8), drop=FALSE]
to_enter_name <- mqtl[,c(2,3), drop=FALSE]
suppressWarnings(storage.mode(to_enter) <- "numeric")
if(sum(is.na(to_enter[,2]))>0 || sum(is.na(to_enter[,1]))>0){
check_qtl <- unique(c(which(is.na(to_enter[,1])), which(is.na(to_enter[,2]))))
for(sample_index in check_qtl){
take_qtl <- which(to_enter_name[sample_index,1] == map[,2])
snp_name_identifier <- to_enter_name[sample_index,1]
if(length(take_qtl)>0){
to_enter[sample_index, 2] <- as.numeric(map[take_qtl[1], 1])
to_enter[sample_index, 1] <- sum(map[1:take_qtl[1], 1] == map[take_qtl[1], 1])
if(verbose) cat("SNP-position assigned via SNP-name\n")
} else{
if(is.na(to_enter[sample_index,2])){
to_enter[sample_index,2] <- 1
if(verbose) cat("Illegal Chromosom-name. Simulate major SNP-effect on chromosome 1.\n")
}
if(!is.na(as.numeric(to_enter_name[sample_index,2]))){
diff_to <- - abs(as.numeric(map[,4]) - as.numeric(to_enter_name[sample_index,2]))
diff_to[map[,1]!=to_enter[sample_index,2]] <- -Inf
take_qtl <- which.max(diff_to)
to_enter[sample_index, 2] <- as.numeric(map[take_qtl[1], 1])
to_enter[sample_index, 1] <- sum(map[1:take_qtl[1], 1] == map[take_qtl[1], 1])
if(verbose) cat("SNP-position assigned via bp position\n")
} else{
to_enter[sample_index,1] <- sample(1:nsnp[to_enter[sample_index,2]],1)
if(verbose) cat(paste0("Illegal SNP-name/bp/position. SNP-effect generated on SNP ", to_enter[sample_index,1],"\n"))
}
if(nchar(snp_name_identifier)>0){
map[to_enter[sample_index,1],2] <- snp_name_identifier
if(verbose) cat(paste0("Major QTL SNP has been renamed to ", snp_name_identifier))
}
}
}
}
to_enter[to_enter[,6]=="",6] <- NA
major_table[[index]] <- to_enter
}
}
}
# Change Allele frequencies according to selection in major QTLs
if(length(major_table)>0){
for(index in 1:length(major_table)){
if(length(major_table[[index]])>0){
for(index2 in 1:nrow(major_table[[index]])){
take <- which(cumsum((map[,1]==major_table[[index]][index2,2]))==major_table[[index]][index2,1])[1]
map[take,5] <- major_table[[index]][index2,6]
majorp <- c(majorp, take)
}
}
}
}
majorp <- unique(majorp)
founder_pedigree <- diag((max(starts)-1) / 2)
prior <- 0
for(subpop in 1:length(subpopulations)){
if(sum(no_data==0 & founder_pop==subpopulation_info[subpop,1])>0){
takes <- which(no_data==0 & founder_pop==subpopulation_info[subpop,1])
take <- NULL
for(index in takes){
take <- c(take, starts[index]:(starts[index+1]-1))
}
if(sum(no_data==1 & founder_pop==subpopulation_info[subpop,1])>0){
takes2 <- which(no_data==1 & founder_pop==subpopulation_info[subpop,1])
take2 <- NULL
for(index in takes2){
take2 <- c(take2, starts[index]:(starts[index+1]-1))
}
if(miraculix.dataset){
p_i[[subpop]] <- miraculix::allele_freq(miraculix::genomicmatrix(dataset[,take2]))*2
} else{
p_i[[subpop]] <- rowMeans(dataset[,take2])
}
} else{
p_i[[subpop]] <- stats::rbeta(sum(nsnp), shape1=as.numeric(subpopulation_info[subpop,2]),
shape2=as.numeric(subpopulation_info[subpop,3]))
set_zero <- sample(1:sum(nsnp), sum(nsnp) * as.numeric(subpopulation_info[subpop,4]))
if(length(set_zero)>0){
set_one <- sample((1:sum(nsnp))[-set_zero], sum(nsnp) * as.numeric(subpopulation_info[subpop,5]))
} else{
set_one <- sample((1:sum(nsnp)), sum(nsnp) * as.numeric(subpopulation_info[subpop,5]))
}
if(length(set_zero)>0){
p_i[[subpop]][set_zero] <- 0
}
if(length(set_one)>0){
p_i[[subpop]][set_one] <- 1
}
}
if(sum(!is.na(map[,5]))){
p_i[[subpop]][!is.na(map[,5])] <- as.numeric(map[!is.na(map[,5]),5])
}
## Manually selected replaces!
major_sub[[index]] <- subpopulations[[subpop]]$'QTL Info'
if(length(major_sub[[index]])>0 && length(major_sub[[index]][[1]])>1){
ncols <- length(major_sub[[index]][[1]])
mqtl_sub <- matrix(unlist(major_sub[[index]]), ncol=ncols, byrow=TRUE)
to_enter_sub <- mqtl_sub[,c(1,4:6), drop=FALSE]
to_enter_name_sub <- mqtl[,c(2,3), drop=FALSE]
suppressWarnings(storage.mode(to_enter_sub) <- "numeric")
if(sum(is.na(to_enter_sub[,2]))>0 || sum(is.na(to_enter_sub[,1]))>0){
check_qtl <- unique(c(which(is.na(to_enter_sub[,1])), which(is.na(to_enter_sub[,2]))))
for(sample_index in check_qtl){
take_qtl <- which(to_enter_name_sub[sample_index,1] == map[,2])
if(length(take_qtl)>0){
to_enter_sub[sample_index, 2] <- as.numeric(map[take_qtl[1], 1])
to_enter_sub[sample_index, 1] <- sum(map[1:take_qtl[1], 1] == map[take_qtl[1], 1])
if(verbose) cat("SNP-position assigned via SNP-name\n")
} else{
if(is.na(to_enter_sub[sample_index,2])){
to_enter_sub[sample_index,2] <- 1
if(verbose) cat("Illegal Chromosom-name. Simulate major SNP-effect on chromosome 1")
}
if(!is.na(as.numeric(to_enter_name_sub[sample_index,2]))){
diff_to <- - abs(as.numeric(map[,4]) - as.numeric(to_enter_name_sub[sample_index,2]))
diff_to[map[,1]!=to_enter_sub[sample_index,2]] <- -Inf
take_qtl <- which.max(diff_to)
to_enter_sub[sample_index, 2] <- as.numeric(map[take_qtl[1], 1])
to_enter_sub[sample_index, 1] <- sum(map[1:take_qtl[1], 1] == map[take_qtl[1], 1])
if(verbose) cat("SNP-position assigned via bp position\n")
} else{
to_enter_sub[sample_index,1] <- sample(1:nsnp[to_enter_sub[sample_index,2]],1)
if(verbose) cat(paste0("Illegal SNP-name/bp/position. SNP-effect generated on SNP ", to_enter_sub[sample_index,1],"\n"))
}
}
}
}
to_enter_sub[to_enter_sub[,4]=="",4] <- NA
major_table_sub[[index]] <- to_enter_sub
position_change <- cumsum(c(0, nsnp))[to_enter_sub[,2]] + to_enter_sub[,1]
p_i[[subpop]][position_change] <- to_enter_sub[,3]
majorp <- c(majorp, position_change)
}
} else{
p_i[[subpop]] <- rep(0, nrow(map))
takes <- which(no_data==0 & founder_pop==subpopulation_info[subpop,1])
take <- NULL
}
majorp <- unique(majorp)
if(miraculix.dataset){
if(length(takes)>0){
nindi <- (starts[max(takes)+1] - starts[min(takes)])/2
chr.nr <- map[,1]
chr.opt <- unique(chr.nr)
dataset <- list()
for(chr_index in 1:length(chr.opt)){
dataset[[chr_index]] <- miraculix::rhaplo(p_i[[subpop]][which(chr.nr==chr.opt[chr_index])], indiv = nindi, loci = nsnp[chr_index])
}
}
} else{
nfinal <- 0
for(index2 in takes){
nfinal <- nfinal + (starts[index2+1]) - starts[index2]
}
if(verbose){cat(paste0("Start LD build-up for subpopulation ", subpopulation_info[subpop,1]))}
ntemp <- as.numeric(subpopulation_info[subpop,9])
if(ntemp>0){
dataset_temp1 <- founder.simulation(nindi = as.numeric(subpopulation_info[subpop,7]),
nsnp = 0,
sex.quota = as.numeric(subpopulation_info[subpop,8]),
n.gen = ntemp,
nfinal = nfinal, verbose=FALSE,
map = map[,1:4],
freq = p_i[[subpop]],
big.output = TRUE,plot=FALSE)
dataset_temp <- dataset_temp1[[1]]
if(subpopulation_info[subpop,10]=="TRUE"){
for(change_p in majorp){
dataset_temp[change_p,] <- stats::rbinom(ncol(dataset_temp),1, prob=p_i[[subpop]][change_p])
}
}
}
for(index in takes){
# Integer * Integer > 2^32 --> NA! conversion to numerics
take <- starts[index]:(starts[index+1]-1)
input_type <- nodes[[founder[index]]]$`Genotype generation`
if(ntemp > 0){
founder_pedigree[take[(1:(length(take)/2)) *2]/2,take[(1:(length(take)/2)) *2]/2] <- dataset_temp1[[4]][1:(length(take)/2)+prior, 1:(length(take)/2)+prior]
}
if(input_type=="Random-sampling"){
if(ntemp>0){
dataset[,take] <- dataset_temp[,prior + 1:length(take)]
} else{
dataset[,take] <- stats::rbinom(as.numeric(nrow(dataset))*as.numeric(length(take)),1, prob=p_i[[subpop]])
}
} else if(input_type=="Fully-homozygous"){
if(ntemp>0){
dataset[,take[c(1:(length(take)/2)*2-1,1:(length(take)/2)*2)]] <- dataset_temp[,prior + 1:(length(take)/2)]
} else{
dataset[,take[c(1:(length(take)/2)*2-1,1:(length(take)/2)*2)]] <- stats::rbinom(as.numeric(nrow(dataset))*as.numeric(length(take))/2,1, prob=p_i[[subpop]])
}
# generate half as many alleles as spots to fill - rest is autocompleted
} else if(input_type=="Fully-heterozygous"){
dataset[,c(1:(length(take)/2)*2)-1] <- 0L # just for safety but not necessary
dataset[,c(1:(length(take)/2)*2)] <- 1L
} else if(input_type=="All-B-Allele"){
dataset[,take] <- 1L
} else if(input_type=="All-A-Allele "){
dataset[,take] <- 0L
} else if(input_type!="Upload Genotypes"){
stop("Invalid input type for founder node!")
}
prior <- prior + length(take)
}
}
}
}
}
############### Generate Base-Population
{
if(!skip.population){
population <- NULL
population <- creating.diploid(dataset = dataset, nindi=length(sex.s),
sex.s = sex.s, genotyped.s = genotyped.s,
chromosome.length = chromo.length,
#snps.equidistant = if(is.na(map[1,4])) {TRUE} else {FALSE},
snps.equidistant = TRUE,
verbose=FALSE,
miraculix = miraculix,
miraculix.dataset = miraculix.dataset,
chr.nr = map[,1], bp=map[,4], snp.name = map[,2],
freq = map[,5], snp.position = if(is.na(map[1,3])) {NULL} else {map[,3]})
if(length(mutation.rate)>0){
population <- set.default(population, parameter.name ="mutation.rate", parameter.value = mutation.rate)
population <- set.default(population, parameter.name ="remutation.rate", parameter.value = mutation.rate)
}
if(length(founder_pedigree)>0 && sum(founder_pedigree) > sum(diag(founder_pedigree))){
population <- add.founder.kinship(population, founder.kinship = founder_pedigree)
}
# Cohort names
gender_founder <- numeric(length(founder))
nr <- 1
for(index in 1:length(founder)){
gender_founder[index] <- sex.s[nr]
nr <- nr + position[founder[index],4]
}
founder_temp2 <- sort(gender_founder,index.return=TRUE)$ix
founder_temp <- founder[sort(gender_founder,index.return=TRUE)$ix]
cohort_info <- cbind(ids[founder_temp], 1, position[founder_temp,4], position[founder_temp,4], position[founder_temp,3], 1+cumsum(c(0,position[founder_temp,4]))[-(length(founder_temp)+1)], 1+cumsum(c(0,position[founder_temp,4]))[-(length(founder_temp)+1)] - sum(sex.s==1),
earliest_time[founder_temp], 0)
cohort_info[(cohort_info[,7]>=1),c(3,6)] <- 0
cohort_info[(cohort_info[,7]<1),c(4,7)] <- 0
founder_pop <- founder_pop[founder_temp2]
population$info$cohorts <- cohort_info
colnames(population$info$cohorts) <- c("name","generation", "male individuals", "female individuals", "class", "position first male", "position first female",
"time point", "creating.type")
if(n_traits>0){
population <- creating.trait(population, n.additive = as.numeric(trait_matrix[,6]),
n.dominant = as.numeric(trait_matrix[,9]),
n.qualitative = as.numeric(trait_matrix[,10]),
n.quantitative = as.numeric(trait_matrix[,11]),
n.equal.additive = as.numeric(trait_matrix[,16]),
n.equal.dominant = as.numeric(trait_matrix[,17]),
dominant.only.positive = as.logical(trait_matrix[,18]),
shuffle.cor = cor_gen, new.phenotype.correlation = cor_pheno,
shuffle.traits=1:n_traits,
trait.name = trait_matrix[,1], verbose=verbose,
is.maternal = as.logical(trait_matrix[,13]),
is.paternal = as.logical(trait_matrix[,14]))
for(index in which(as.logical(trait_matrix[,15]))){
population <- add.combi(population, trait = index, combi.weights = combi_weights[[index]], trait.name = trait_matrix[index,1])
}
# Correct Scaling
snp.before <- cumsum(c(0,population$info$snp))
## REASON?!
if(TRUE){
for(index in 1:n_traits){
if(length(population$info$real.bv.add[[index]])>0){
t <- population$info$real.bv.add[[index]]
take <- sort(t[,1]+ snp.before[t[,2]], index.return=TRUE)
t <- t[take$ix,,drop=FALSE]
take <- sort(t[,1]+ t[,2] * 10^10)
keep <- c(0,which(diff(take)!=0), length(take))
if(length(keep) <= nrow(t)){
for(index2 in 2:(length(keep))){
t[keep[index2],3:5] <- colSums(t[(keep[index2-1]+1):keep[index2],3:5, drop=FALSE])
}
population$info$real.bv.add[[index]] <- t[keep,]
}
}
}
}
## Variance Standardization
for(index in 1:n_traits){
new_var <- as.numeric(trait_matrix[index,4])^2
if(population$info$bv.calculated==FALSE){
population <- breeding.diploid(population, verbose=verbose)
}
active_sub <- subpopulation_info[1,1]
standard_cohort <- population$info$cohorts[which(founder_pop==active_sub),1]
var_test <- stats::var(get.bv(population, cohorts= standard_cohort)[index,])
test1 <- TRUE
if(length(population$info$real.bv.add[[index]])>0){
population$info$real.bv.add[[index]][,3:5] <- population$info$real.bv.add[[index]][,3:5] * sqrt( new_var / var_test)
test1 <- FALSE
}
if(length(population$info$real.bv.mult[[index]])>0){
population$info$real.bv.mult[[index]][,5:13] <- population$info$real.bv.mult[[index]][,5:13] * sqrt( new_var / var_test)
test1 <- FALSE
}
if(test1 && verbose) cat("You entered a trait without quantitative loci. Is this intentional?\n")
}
population$info$bv.calculated <- FALSE
## Mean Standardization
for(index in 1:n_traits){
if(population$info$bv.calculated==FALSE){
population <- breeding.diploid(population, verbose=verbose)
}
active_sub <- subpopulation_info[1,1]
standard_cohort <- population$info$cohorts[which(founder_pop==active_sub),1]
mean_test <- mean(get.bv(population, cohorts= standard_cohort)[index,])
population$info$base.bv[index] <- traitmean[index] + population$info$base.bv[index] - mean_test
}
population$info$bv.calculated <- FALSE
## Add Major QTL
for(index in 1:n_traits){
if(length(major_table)>=index && length(major_table[[index]])>0){
to_enter <- major_table[[index]][,1:5,drop=FALSE]
} else{
to_enter <- rbind(NULL)
}
if(length(to_enter)>0 || length(population$info$real.bv.add[[index]])>0){
population$info$real.bv.add[[index]] <- rbind(to_enter, population$info$real.bv.add[[index]])
if(population$info$real.bv.length[1]<index && length(to_enter)>0){
population$info$real.bv.length[1] <- index
}
}
}
## Subpopulation differences
if(length(subpopulations)>1){
for(subpop in 1:length(subpopulations)){
if(population$info$bv.calculated==FALSE){
population <- breeding.diploid(population, verbose=verbose)
}
cutoffsub <- max(which(names(subpopulations[[subpop]])=="majorfreq"), which(names(subpopulations[[subpop]])=="QTL Info"))
if(length(subpopulations[[subpop]])>cutoffsub){
for(indexmod in (cutoffsub+1):length(subpopulations[[subpop]])){
if(subpopulations[[subpop]][[indexmod]]!=""){
active_trait <- as.numeric(substr(names(subpopulations[[subpop]])[[indexmod]], start=2, stop=nchar(names(subpopulations[[subpop]])[[indexmod]])))
valid_markers0 <- valid_markers1 <- which(p_i[[subpop]]>0 & p_i[[subpop]]<1)
for(reduction in (1:length(subpopulations))[-subpop]){
valid_markers0 <- intersect(valid_markers0, which(p_i[[reduction]]==0))
}
for(reduction in (1:length(subpopulations))[-subpop]){
valid_markers1 <- intersect(valid_markers1, which(p_i[[reduction]]==1))
}
valid_markers <- c(valid_markers0, valid_markers1)
subpop_name <- subpopulation_info[subpop,1]
sub_cohort <- population$info$cohorts[which(founder_pop==subpop_name),1]
if(length(sub_cohort)>0){
mean_ref <- mean(get.bv(population, cohorts= standard_cohort )[active_trait,])
mean_sub <- mean(get.bv(population, cohorts= sub_cohort )[active_trait,])
current_diff <- mean_sub - mean_ref
target_diff <- as.numeric(subpopulations[[subpop]][[indexmod]])
if(subpop==1){
change <- target_diff
} else{
change <- target_diff - current_diff
}
diff_freq <- c(p_i[[subpop]][valid_markers0], - p_i[[subpop]][valid_markers1])
effect_size <- change / sum(abs(diff_freq)) / 2
direction <- diff_freq > 0
snp_index <- chromo_index <- numeric(length(valid_markers))
for(index2 in 1:length(valid_markers)){
chromo_index[index2] <- max(which(c(0,population$info$cumsnp)<=valid_markers[index2]))
snp_index[index2] <- valid_markers[index2] - c(0,population$info$cumsnp)[chromo_index[index2]]
}
add.effects <- cbind(snp_index, chromo_index, 2* effect_size, effect_size, 0)
add.effects[direction==TRUE,3:5] <- cbind(rep(0, sum(direction)), rep(effect_size, sum(direction)), rep(2*effect_size, sum(direction)))
population$info$real.bv.add[[active_trait]] <- rbind(population$info$real.bv.add[[active_trait]], add.effects)
population$info$bv.calculated <- FALSE
if(subpop==1){
population <- breeding.diploid(population, verbose=verbose)
population$info$bv.calculated <- FALSE
}
}
}
}
}
}
}
if(length(trafos)>0){
for(index in 1:length(trafos)){
if(length(trafos[[index]])>0){
population <- creating.phenotypic.transform(population, phenotypic.transform.function = trafos[[index]], trait = index)
}
}
}
}
population$breeding[[1]][[5]] <- mig_m
population$breeding[[1]][[6]] <- mig_f
}
}
################## Last preps ########################
{
# Add edges info to nodes
for(index in 1:length(edges)){
to_node <- which(edges[[index]]$to==ids)
nodes[[to_node]]$'Breeding Type' <- edges[[index]]$'Breeding Type'
if(nodes[[to_node]]$'Breeding Type'=="Selection" || nodes[[to_node]]$'Breeding Type'=="Split" || nodes[[to_node]]$'Breeding Type'=="Aging" ){
nodes[[to_node]]$'Selection Type' <- edges[[index]]$'Selection Type'
if(length(nodes[[to_node]]$'Selection Type')>0 && nodes[[to_node]]$'Selection Type'=="Pseudo-BVE"){
pseudo <- NULL
for(t in as.character(1:n_traits-1)){
pseudo <- c(pseudo, edges[[index]][[t]])
}
nodes[[to_node]]$'PseudoAcc' <- as.numeric(pseudo)
}
# nodes[[to_node]]$proportion <- edges[[index]]$proportion # not needed?
nodes[[to_node]]$origin <- edges[[index]]$from
nodes[[to_node]]$last_available <- edges[[index]]$last_available
if(edges[[index]]$last_available){
temp1 <- unlist(strsplit(edges[[index]]$from, split="_"))
temp1 <- temp1[1:max(1,length(temp1)-1)]
temp1 <- temp2 <- paste0(temp1)
if(sum(repeat_links[,1]==temp1)>0){
temp1 <- repeat_links[which(repeat_links[,1]==temp1),2]
}
if(sum(founder == which(ids==temp1))==0){
nodes[[to_node]]$require <- c(nodes[[to_node]]$require, temp2)
}
}
nodes[[to_node]]$'Relationship Matrix' <- edges[[index]]$'Relationship Matrix'
nodes[[to_node]]$skip <- edges[[index]]$skip
nodes[[to_node]]$'BVE Method' <- edges[[index]]$'BVE Method'
nodes[[to_node]]$'MAS_marker' <- edges[[index]]$'MAS_marker'
nodes[[to_node]]$'Use Offspring for BVE' <- edges[[index]]$'Use Offspring for BVE'
nodes[[to_node]]$phenotype_used <- edges[[index]]$phenotype_used
nodes[[to_node]]$edge.nr <- c(nodes[[to_node]]$edge.nr,index)
nodes[[to_node]]$'Time Needed' <- c(nodes[[to_node]]$'Time Needed',edges[[index]]$'Time Needed')
nodes[[to_node]]$'Selection Index' <- edges[[index]]$'Selection Index'
nodes[[to_node]]$'Cohorts used in BVE' <- edges[[index]]$'Cohorts used in BVE'
nodes[[to_node]]$'Selection Type Mean' <- edges[[index]]$'Selection Type Mean'
if(nodes[[to_node]]$'Cohorts used in BVE'=="Manual select"){
edges[[index]]$`Manuel selected cohorts`
repeat_iter_test <- strsplit(nodes[[to_node]]$id, "_")[[1]]
repeat_iter <- suppressWarnings(as.numeric(repeat_iter_test[length(repeat_iter_test)]))
if(is.na(repeat_iter)){
repeat_iter <- 0
}
id_bve <- NULL
for(index3 in 1:length(edges[[index]]$`Manuel selected cohorts`)){
cohort_split <- strsplit(edges[[index]]$`Manuel selected cohorts`[[index3]], ":")[[1]]
if(length(cohort_split)==1){
cohort_iter <- -(0:20) + repeat_iter
cohort_iter <- cohort_iter[cohort_iter>=0]
} else{
if(!is.na(cohort_split[length(cohort_split)]) && cohort_split[length(cohort_split)]=="X"){
cohort_iter <- (- repeat_iter)
} else{
cohort_iter <- as.numeric(cohort_split[length(cohort_split)])
if(is.na(cohort_iter)){
cohort_iter <- 0
}
}
cohort_iter <- repeat_iter + cohort_iter
}
for(index4 in cohort_iter){
if(index4< (-1)){
} else if(index4 == (-1)){
take_repeat <- which(repeat_links[,1]==cohort_split[1])
if(length(take_repeat)==1){
id_bve <- c(id_bve, repeat_links[take_repeat,2])
}
} else if(index4==0){
id_bve <- c(id_bve, cohort_split[1])
} else{
id_bve <- c(id_bve, paste0(cohort_split[1],"_", cohort_iter))
}
}
}
clean_up <- which(duplicated(c(ids, id_bve))[-(1:length(ids))])
id_bve <- id_bve[clean_up]
id_bve <- intersect(ids, id_bve)
nodes[[to_node]]$'Manuel selected cohorts' <- id_bve
}
nodes[[to_node]]$'threshold' <- edges[[index]]$'threshold'
nodes[[to_node]]$'bve_solve' <- edges[[index]]$'BVE Method_solver'
nodes[[to_node]]$'threshold_sign' <- edges[[index]]$'threshold_sign'
nodes[[to_node]]$'Depth of Pedigree' <- edges[[index]]$'Depth of Pedigree'
}
if(nodes[[to_node]]$'Breeding Type'=="Reproduction" || nodes[[to_node]]$'Breeding Type'=="Selfing" || nodes[[to_node]]$'Breeding Type'=="DH-Production" || nodes[[to_node]]$'Breeding Type'=="Cloning"){
nodes[[to_node]]$origin <- c(nodes[[to_node]]$origin,edges[[index]]$from)
nodes[[to_node]]$last_available <- c(nodes[[to_node]]$last_available, edges[[index]]$last_available)
if(edges[[index]]$last_available){
temp1 <- unlist(strsplit(edges[[index]]$from, split="_"))
temp1 <- temp1[1:max(1,length(temp1)-1)]
temp1 <- temp2 <- paste0(temp1)
if(sum(repeat_links[,1]==temp1)>0){
temp1 <- repeat_links[which(repeat_links[,1]==temp1),2]
}
if(sum(founder == which(ids==temp1))==0){
nodes[[to_node]]$require <- c(nodes[[to_node]]$require, temp2)
}
}
if(nodes[[to_node]]$'Breeding Type' == "Reproduction"){
if(length(nodes[[to_node]]$full_sib)==0 || nodes[[to_node]]$full_sib == FALSE){
nodes[[to_node]]$full_sib <- edges[[index]]$full_sib
}
if(length(nodes[[to_node]]$half_sib)==0 || nodes[[to_node]]$half_sib == FALSE){
nodes[[to_node]]$half_sib <- edges[[index]]$half_sib
}
}
nodes[[to_node]]$edge.nr <- c(nodes[[to_node]]$edge.nr,index)
nodes[[to_node]]$'Time Needed' <- c(nodes[[to_node]]$'Time Needed',edges[[index]]$'Time Needed')
if(length(edges[[index]]$OGC)>0 && edges[[index]]$OGC=="Yes"){
nodes[[to_node]]$OGC <- TRUE
if(length(edges[[index]]$'ogc_target')>0 || length(nodes[[to_node]]$ogc_cAc)==0){
nodes[[to_node]]$ogc_target <- (edges[[index]]$'ogc_target')
nodes[[to_node]]$ogc_relation <- (edges[[index]]$'ogc_relation')
if(length(nodes[[to_node]]$ogc_relation )==0){
nodes[[to_node]]$ogc_relation <- "vanRaden"
}
if(length(nodes[[to_node]]$ogc_relation) && nodes[[to_node]]$ogc_relation == "Pedigree"){
nodes[[to_node]]$ogc_relation <- "pedigree"
}
if(length( nodes[[to_node]]$constrain) == 0 ){
nodes[[to_node]]$constrain <- list()
}
constrains <- list(edges[[index]]$'ogc_constrain1', edges[[index]]$'ogc_constrain2', edges[[index]]$'ogc_constrain3')
constrains_value <- list(edges[[index]]$'ogc_constrain1_value', edges[[index]]$'ogc_constrain2_value', edges[[index]]$'ogc_constrain3_value')
for(index in 1:3){
if(constrains[[index]]=="ub.BV"){
nodes[[to_node]]$constrain[[1]] <- as.numeric(constrains_value[[index]])
}
if(constrains[[index]]=="eq.BV"){
nodes[[to_node]]$constrain[[2]] <- as.numeric(constrains_value[[index]])
}
if(constrains[[index]]=="lb.BV"){
nodes[[to_node]]$constrain[[3]] <- as.numeric(constrains_value[[index]])
}
if(constrains[[index]]=="ub.sKin"){
nodes[[to_node]]$constrain[[4]] <- as.numeric(constrains_value[[index]])
}
if(constrains[[index]]=="uniform"){
nodes[[to_node]]$constrain[[5]] <- constrains_value[[index]]
}
if(constrains[[index]]=="lb.BV.increase"){
nodes[[to_node]]$constrain[[6]] <- as.numeric(constrains_value[[index]])
}
if(constrains[[index]]=="ub.sKin.increase"){
nodes[[to_node]]$constrain[[7]] <- as.numeric(constrains_value[[index]])
}
if(length( nodes[[to_node]]$constrain)>0){
nodes[[to_node]]$constrain[[8]] <- "placeholder"
}
}
}
} else if(length(nodes[[to_node]]$OGC)==0){
nodes[[to_node]]$OGC <- FALSE
}
if(length(nodes[[to_node]]$selection_ratio)==0){
nodes[[to_node]]$selection_ratio <- c(1,1)
}
if(length(nodes[[to_node]]$selection_ratio_type)==0){
nodes[[to_node]]$selection_ratio_type <- c("bv","bv")
}
if(length(nodes[[to_node]]$selection_ratio_index)==0){
if(n_traits>0){
nodes[[to_node]]$selection_ratio_index <- c(selection_index_name[1], selection_index_name[1])
}
}
if(length(edges[[index]]$'selection_ratio')>0){
sex <- as.numeric(nodes[[which(ids==edges[[index]]$from)]]$Sex=="Female") + 1
if(length(edges[[index]]$'selection_ratio')==1){
nodes[[to_node]]$selection_ratio[sex] <- as.numeric(edges[[index]]$'selection_ratio')
if(length(edges[[index]]$'selection_ratio_type')>0){
if(edges[[index]]$'selection_ratio_type'=="Genomic Value"){
nodes[[to_node]]$selection_ratio_type[sex] <- "bv"
} else if(edges[[index]]$'selection_ratio_type'=="Breeding Value"){
nodes[[to_node]]$selection_ratio_type[sex] <- "bve"
} else if(edges[[index]]$'selection_ratio_type'=="Phenotype"){
nodes[[to_node]]$selection_ratio_type[sex] <- "pheno"
}
}
if(is.na(nodes[[to_node]]$selection_ratio[sex])){
nodes[[to_node]]$selection_ratio[sex] <- 1
}
if(is.na(nodes[[to_node]]$selection_ratio_type[sex])){
nodes[[to_node]]$selection_ratio[sex] <- "bv"
}
}
}
if(length(edges[[index]]$'selection_ratio_index')>0){
sex <- as.numeric(nodes[[which(ids==edges[[index]]$from)]]$Sex=="Female") + 1
nodes[[to_node]]$selection_ratio_index[sex] <- (edges[[index]]$'selection_ratio_index'$Name)
}
}
if(nodes[[to_node]]$'Breeding Type'=="Combine"){
nodes[[to_node]]$origin <- c(nodes[[to_node]]$origin,edges[[index]]$from)
nodes[[to_node]]$last_available <- c(nodes[[to_node]]$last_available, edges[[index]]$last_available)
if(edges[[index]]$last_available){
temp1 <- unlist(strsplit(edges[[index]]$from, split="_"))
temp1 <- temp1[1:max(1,length(temp1)-1)]
temp1 <- temp2 <- paste0(temp1)
if(sum(repeat_links[,1]==temp1)>0){
temp1 <- repeat_links[which(repeat_links[,1]==temp1),2]
}
if(sum(founder == which(ids==temp1))==0){
nodes[[to_node]]$require <- c(nodes[[to_node]]$require, temp2)
}
}
nodes[[to_node]]$edge.nr <- c(nodes[[to_node]]$edge.nr,index)
nodes[[to_node]]$'Time Needed' <- c(nodes[[to_node]]$'Time Needed',edges[[index]]$'Time Needed')
}
if(length(nodes[[to_node]]$max_offspring)!=2){
nodes[[to_node]]$max_offspring <- c(Inf, Inf)
}
if(length(nodes[[to_node]]$max_offspring_pair)!=2){
nodes[[to_node]]$max_offspring_pair <- c(Inf, Inf)
}
if(length(edges[[index]]$'Max Offspring')>0 && nchar(edges[[index]]$'Max Offspring')>0){
from_node <- which(ids==edges[[index]]$from)
switch <- as.numeric(nodes[[from_node]]$Sex=="Female") +1
nodes[[to_node]]$max_offspring[switch] <- as.numeric(edges[[index]]$'Max Offspring')
}
if(length(edges[[index]]$'Max Offspring Pair')>0 && nchar(edges[[index]]$'Max Offspring Pair')>0){
from_node <- which(ids==edges[[index]]$from)
switch <- as.numeric(nodes[[from_node]]$Sex=="Female") +1
nodes[[to_node]]$max_offspring_pair[switch] <- as.numeric(edges[[index]]$'Max Offspring Pair')
}
if(length(edges[[index]]$repeat_mating)>0){
nodes[[to_node]]$repeat_mating <- as.numeric(edges[[index]]$repeat_mating)
} else{
if(length(nodes[[to_node]]$repeat_mating)==0){
nodes[[to_node]]$repeat_mating <- 1
}
}
}
phenotype_groups <- numeric(length(nodes))
for(index in 1:length(edges)){
phenotype_groups[which(ids==edges[[index]]$from)] <- 1
}
priority_breeding <- ids[(1-phenotype_groups)*1:length(ids)]
n_tester <- n_tester_generated <- numeric(length(nodes))
for(index in 1:length(edges)){
if(length(intersect(priority_breeding, edges[[index]]$to))){
n_tester[which(ids==edges[[index]]$from)] <- n_tester[which(ids==edges[[index]]$from)] +1
}
}
for(index in 1:length(nodes)){
exp_genotyped <- as.numeric(nodes[[index]]$`Proportion of genotyped individuals`)
real_genotyped <- 0
n_animals <- 0
for(index2 in nodes[[index]]$origin){
real_genotyped <- real_genotyped + as.numeric(nodes[[which(index2==ids)]]$`Proportion of genotyped individuals`) * as.numeric(nodes[[which(index2==ids)]]$`Number of Individuals`)
n_animals <- n_animals + nodes[[which(index2==ids)]]$`Number of Individuals`
}
if(n_animals>0){
nodes[[index]]$'Proportion of added genotypes' <- max(0, exp_genotyped - real_genotyped/n_animals)
} else{
nodes[[index]]$'Proportion of added genotypes' <- exp_genotyped
}
}
}
# Check for Split nodes
{
to_split <- NULL
split_info <- list()
split_part <- list()
for(index in 1:length(edges)){
if(edges[[index]]$'Breeding Type'=="Split"){
to_split <- unique(c(to_split, edges[[index]]$from))
split_nr <- which(to_split==edges[[index]]$from)
nodes_nr <- which(edges[[index]]$from==ids)
split_info[[split_nr]] <- 1:nodes[[nodes_nr]]$'Number of Individuals'
if(length(split_part)>= split_nr){
split_part[[split_nr]] <- c(split_part[[split_nr]], edges[[index]]$to)
} else{
split_part[[split_nr]] <- edges[[index]]$to
}
}
}
}
############## Founder Phenotypes ########################
{
if(n_traits>0 & !skip.population){
for(index in 1:nrow(population$info$cohorts)){
population <- breeding.diploid(population, heritability = heritability,
repeatability = repeatability,
new.bv.observation.cohorts = population$info$cohorts[index,1],
sigma.e.database = cbind(1,(1:2)[population$info$size[1,]>0]),
n.observation = pheno_index[which(pheno_index_name==nodes[[which(ids==population$info$cohorts[index,1])]]$'Phenotyping Class'),],
verbose=verbose)
}
}
}
#### Add litter size
population$info$repeat.mating <- repeat.mating
population$info$repeat.mating.copy <- repeat.mating.copy
### derive order of generation
{
last_avail <- rep(FALSE, length(edges))
for(index in 1:length(edges)){
if(length(edges[[index]]$last_available)>0){
last_avail[index] <- edges[[index]]$last_available
}
}
ids_type <- ids
ids_rep <- numeric(length(ids))
temp1 <- strsplit(ids_type, split="_")
for(index in 1:length(ids_type)){
if(length(temp1[[index]])==1){
ids_rep[index] <- 0
ids_type[index] <- (temp1[[index]][1])
} else if(suppressWarnings(is.na(as.numeric(temp1[[index]][length(temp1[[index]])])))){
ids_rep[index] <- 0
text1 <- NULL
for(merge in 1:length(temp1[[index]])){
if(merge==length(temp1[[index]])){
text1 <- paste0(text1, temp1[[index]][merge])
} else{
text1 <- paste0(text1 ,temp1[[index]][merge], "_")
}
}
ids_type[index] <- text1
} else {
ids_rep[index] <- as.numeric(temp1[[index]][length(temp1[[index]])])
text1 <- NULL
for(merge in 1:(length(temp1[[index]])-1)){
if(merge==(length(temp1[[index]])-1)){
text1 <- paste0(text1, temp1[[index]][merge])
} else{
text1 <- paste0(text1 ,temp1[[index]][merge], "_")
}
}
ids_type[index] <- text1
}
}
ids_rep <- as.numeric(ids_rep)
{
simulated <- founder
left <- (1:groups)[-simulated]
generation <- 1
time.point.list <- numeric(length(ids))
for(index in founder){
time.point.list[index] <- nodes[[index]]$earliest_time
}
possible <- NULL
while(length(left)>0){
generation <- generation + 1
if(TRUE || length(possible)==0){
possible <- ids[left]
left1 <- left
} else{
suppressWarnings(test <- as.numeric(unlist(strsplit(possible, split="_"))))
min_rep <- min(test, na.rm=TRUE)
max_rep <- max(test, na.rm=TRUE)
if(min_rep==Inf){
min_rep <- 0
}
if(max_rep==(-Inf)){
max_rep <- Inf
}
possible <- ids[left]
rep_nr <- ids_rep[left]
keep <- rep_nr < 10 | (rep_nr > (min_rep-15) & rep_nr < (max_rep+15))
left1 <- left[keep]
possible <- possible[left1]
}
stock <- ids[-left]
# require-check
require_remove <- rep(TRUE, length(left1))
for(index in 1:length(left1)){
if(length(nodes[[left1[index]]]$require)>0){
for(check_coh in nodes[[left1[index]]]$require){
if(sum(stock==check_coh)==0){
require_remove[index] <- FALSE
}
}
}
}
possible <- possible[require_remove]
for(index in (1:length(edges))[!last_avail]){
there <- which(edges[[index]]$to==possible)
if(length(there)>0){
if(sum(edges[[index]]$from==stock)==0){
possible <- possible[-there]
next
}
to_node <- which(ids==edges[[index]]$to)
if(length(nodes[[to_node]]$'Manuel selected cohorts')>0){
manu <- setdiff(nodes[[to_node]]$'Manuel selected cohorts', edges[[index]]$from)
if(length(intersect(manu, stock))<length(manu)){
possible <- possible[-there]
}
}
}
}
if(length(intersect(possible, priority_breeding))>0){
possible <- intersect(possible, priority_breeding)
} else{
stock <- ids[-unique(c(left, (n_tester>n_tester_generated)*(1:length(n_tester))))]
for(index in 1:length(edges)){
there <- which(edges[[index]]$to==possible)
if(length(there)>0){
if(sum(edges[[index]]$from==stock)==0){
possible <- possible
}
}
}
}
# Remove group for which not all testers are generated
if(length(possible)==0){
stop("invalite breeding program")
}
for(group in possible){
groupnr <- which(ids==group)
simulated <- c(simulated, groupnr)
time.point <- 0
origins <- nodes[[which(ids==group)]]$origin
time_needed <- as.numeric(nodes[[groupnr]]$'Time Needed')
for(temp1 in (1:length(origins))[!nodes[[groupnr]]$last_available]){
time.point <- max(time.point.list[simulated][which(origins[temp1]==ids[simulated])] + time_needed[temp1],time.point)
}
if(length(nodes[[groupnr]]$require)>0){
for(coh in nodes[[groupnr]]$require){
time.point <- max(time.point, time.point.list[simulated][which(coh==ids[simulated])]+0.0000001)
}
}
for(temp1 in nodes[[groupnr]]$'Manuel selected cohorts'){
groupnr2 <- which(ids==temp1)
time.point <- max(time.point.list[groupnr2],time.point)
}
time.point <- max(nodes[[which(ids==group)]]$earliest_time, time.point)
time.point.list[groupnr] <- time.point
}
left <- (1:groups)[-simulated]
}
generation_times <- sort(unique(time.point.list))
generation_group <- list()
generation_bv_size <- list()
for(index in 1:length(generation_times)){
nrs <- setdiff(which(time.point.list==generation_times[[index]]), founder)
btype <- numeric(length(nrs))
if(length(nrs)>0){
for(index2 in 1:length(nrs)){
btype[index2] <- nodes[[nrs[index2]]]$'Breeding Type'
}
}
prio <- which(btype=="Reproduction"| btype=="Selfing" | btype=="DH-Production" | btype=="Cloning" | btype=="Combine")
if(length(prio)>0){
nrs <- c(nrs[prio], nrs[-prio])
btype <- c(btype[prio], btype[-prio])
}
# Both sorting ###
if(length(splits)>0){
reorder <- numeric(length(nrs))
next1 <- 1
for(index2 in nrs){
split1 <- strsplit(ids[index2], split = "_")[[1]][1]
split3 <- strsplit(ids[index2], split = "_")[[1]][3]
if(sum( split1 == splits)>0){
if(strsplit(ids[index2], split = "_")[[1]][2]=="F"){
next
} else if(strsplit(ids[index2], split = "_")[[1]][2]=="M"){
reorder[next1] <- which(nrs==index2)
for(index3 in nrs){
if(length(strsplit(ids[index3], split = "_")[[1]])>1 && strsplit(ids[index3], split = "_")[[1]][2]=="F" && strsplit(ids[index3], split = "_")[[1]][1]==split1 &&
((is.na(split3) && length(strsplit(ids[index3], split = "_")[[1]])<3) ||
((!is.na(split3) && length(strsplit(ids[index3], split = "_")[[1]])==3) && strsplit(ids[index3], split = "_")[[1]][3] == strsplit(ids[index2], split = "_")[[1]][3])
)){
partner <- index3
}
}
reorder[next1+1] <- which(nrs==partner)
next1 <- next1 + 2
} else{
stop("Something wrong in Both-nodes")
}
} else{
reorder[next1] <- which(nrs==index2)
next1 <- next1 + 1
}
}
nrs <- nrs[reorder]
btype <- btype[reorder]
}
if(sum(btype=="Split")>1){
activ_type <- which(btype=="Split")
activ_node <- nrs[activ_type]
ns <- length(activ_node)
activ_prior <- activ_time <- numeric(ns)
for(index2 in 1:ns){
activ_prior[index2] <- nodes[[activ_node[index2]]]$origin
activ_time[index2] <- nodes[[activ_node[index2]]]$`Time Needed`
}
avail_prior <- unique(activ_prior)
for(index2 in 1:length(avail_prior)){
activ <- which(activ_prior==avail_prior[index2])
new_order <- activ[sort(activ_time[activ], index.return=TRUE)$ix]
nrs[activ_type][activ] <- nrs[activ_type][new_order]
btype[activ_type][activ] <- btype[activ_type][new_order]
}
}
generation_group[[index]] <- ids[nrs]
changes <- 1
while(changes>0){
changes <- 0
if(length(nrs)>0){
for(index2 in 1:length(nrs)){
if(length(intersect(c(nodes[[nrs[index2]]]$origin, nodes[[nrs[index2]]]$'Manuel selected cohorts'), generation_group[[index]][index2:length(nrs)]))>0){
back <- 0
for(i_origin in c(nodes[[nrs[index2]]]$origin, nodes[[nrs[index2]]]$'Manuel selected cohorts')){
back1 <- which(generation_group[[index]]==i_origin)
back <- max(back, back1[length(back1)])
}
store_temp <- nrs[index2]
nrs[index2:(back-1)] <- nrs[(index2+1):back]
nrs[back] <- store_temp
generation_group[[index]] <- ids[nrs]
changes <- 1
}
}
}
}
btype <- numeric(length(nrs))
if(length(nrs)>0){
for(index2 in 1:length(nrs)){
btype[index2] <- nodes[[nrs[index2]]]$'Breeding Type'
}
}
database_add <- c(0,0)
while(length(btype)>0 && (btype[1]=="Reproduction"| btype[1]=="Selfing" | btype[1]=="DH-Production" | btype[1]=="Cloning" | btype[1]=="Combine")){
sex <- as.numeric(nodes[[nrs[1]]]$'Sex'=="Female") + 1
database_add <- database_add + as.numeric(nodes[[nrs[1]]]$`Number of Individuals`) * c(sex==1, sex==2)
btype <- btype[-1]
nrs <- nrs[-1]
}
generation_bv_size[[index]] <- database_add
}
}
{
for(index in length(generation_group):1){
if(length(generation_group[[index]])==0){
generation_group[[index]] <- NULL
generation_bv_size[[index]] <- NULL
}
}
}
generation_both <- list()
for(index in 1:length(generation_group)){
generation_both[[index]] <- rep(0, length(generation_group[[index]]))
for(index2 in 1:length(generation_group[[index]])){
if(sum(strsplit(generation_group[[index]][index2], split = "_")[[1]][1] == splits)>0){
test1 <- nodes[[which(ids==generation_group[[index]][index2])]]$'Breeding Type'
if(test1=="Reproduction" || test1 =="Selfing" || test1 == "DH-Production" || test1 =="Cloning"){
if(strsplit(generation_group[[index]][index2], split = "_")[[1]][2]=="M"){
generation_both[[index]][index2] <- 1
} else if(strsplit(generation_group[[index]][index2], split = "_")[[1]][2]=="F"){
generation_both[[index]][index2] <- 2
}
}
}
}
if(sum(generation_both[[index]]==1) != sum(generation_both[[index]]==2)){
stop("Something wrong in Both generation!")
}
}
for(index1 in (1:length(nodes))[-founder]){
for(index2 in which(nodes[[index1]]$last_available)){
cohort_type <- ids_type[which(nodes[[index1]]$origin[index2]==ids)]
class_time <- time.point.list[which(ids_type==cohort_type)] - time.point.list[index1]
class_time[class_time>=0] <- -Inf
take_repeat <- which.max(class_time)
if(max(class_time) == (-Inf) && length(class_time)>1){
if(length(which(repeat_links[,1]==cohort_type))!=1){
warning("Problem in UseLastAvailable!\n")
stop("No automatic fix possible")
}
nodes[[index1]]$origin[index2] <- repeat_links[which(repeat_links[,1]==cohort_type),2]
} else{
nodes[[index1]]$origin[index2] <- ids[which(cohort_type==ids_type)][take_repeat]
}
}
}
#### Check if same breeding value estimation is done multiple times
{
bv_cohort <- NULL
bv_cohort_ori <- NULL
bv_included <- list()
for(index in 1:length(nodes)){
bve.breeding.type <- nodes[[index]]$`Breeding Type`=="Selection" || nodes[[index]]$`Breeding Type`=="Aging" || nodes[[index]]$`Breeding Type`=="Split"
if(bve.breeding.type && nodes[[index]]$`Selection Type`=="BVE" && length(nodes[[index]]$'Cohorts used in BVE')>0 && nodes[[index]]$'Cohorts used in BVE'=="Manual select"){
bv_cohort <- c(bv_cohort, nodes[[index]]$id)
bv_cohort_ori <- c(bv_cohort_ori, nodes[[index]]$origin)
bv_included[[length(bv_included)+1]] <- sort(unique( c(nodes[[index]]$origin, nodes[[index]]$'Manuel selected cohorts')))
}
}
if(length(bv_included)>0){
bv_class <- unique(bv_included)
bv_cohort_class <- numeric(length(bv_cohort))
bv_class_member <- list()
for(index in 1:length(bv_cohort)){
for(index2 in 1:length(bv_class)){
if(length(bv_class[[index2]])==length(bv_included[[index]]) && prod(bv_class[[index2]] == bv_included[[index]])){
bv_cohort_class[index] <- index2
if(length(bv_class_member)<index2){
bv_class_member[[index2]] <- bv_cohort_ori[index]
} else{
bv_class_member[[index2]] <- c(bv_class_member[[index2]], bv_cohort_ori[index])
}
}
}
}
class_executed <- rep(FALSE,length(bv_class))
}
}
}
############## Attach json-infos ########################
housing <- list(housing_index, housing_index_name)
phenotyping <- list(pheno_index_costs, pheno_index_name, pheno_index)
###### Write List of Cohorts + costs #####################
if(TRUE){
costdata <- NULL
# Name of Cohort , Time-point, total money spent, cost of genotyping, cost of phenotyping, cost of housing
for(index in founder){
group_name <- ids[index]
groupnr <- which(ids==group_name)
group_time <- time.point.list[groupnr]
group_size <- nodes[[groupnr]]$`Number of Individuals`
share_geno <- as.numeric(nodes[[groupnr]]$`Proportion of genotyped individuals`)
if(length(share_geno)==0) share_geno <- 1
cost_geno <- nodes[[groupnr]]$`Number of Individuals` * genotyping_costs * as.numeric(share_geno)
pheno_class <- which(nodes[[groupnr]]$`Phenotyping Class`==phenotyping[[2]])
if(length(pheno_class)==0){
pheno_class <- 1
}
cost_pheno <- nodes[[groupnr]]$`Number of Individuals` * phenotyping[[1]][pheno_class]
housing_class <- which(nodes[[groupnr]]$`Housing Cost Class`==housing[[2]])
if(length(housing_class)==0){
housing_class <- 1
}
cost_housing <- nodes[[groupnr]]$`Number of Individuals` * housing[[1]][housing_class]
if(n_traits==0){
costdata <- rbind(costdata, c(group_name, group_size, group_time, cost_geno+cost_housing, cost_geno, cost_housing))
} else{
costdata <- rbind(costdata, c(group_name, group_size, group_time, cost_geno+cost_housing+cost_pheno, cost_geno, cost_pheno, cost_housing))
}
}
for(index in 1:length(generation_group)){
if(length(generation_group[[index]])){
for(index2 in 1:length(generation_group[[index]])){
group_name <- generation_group[[index]][[index2]]
groupnr <- which(ids==group_name)
group_time <- time.point.list[groupnr]
group_size <- nodes[[groupnr]]$`Number of Individuals`
share_geno <- as.numeric(nodes[[groupnr]]$`Proportion of genotyped individuals`)
if(length(share_geno)==0) share_geno <- 1
cost_geno <- nodes[[groupnr]]$`Number of Individuals` * genotyping_costs * share_geno
pheno_class <- which(nodes[[groupnr]]$`Phenotyping Class`==phenotyping[[2]])
if(length(pheno_class)==0){
pheno_class <- 1
}
cost_pheno <- nodes[[groupnr]]$`Number of Individuals` * phenotyping[[1]][pheno_class]
housing_class <- which(nodes[[groupnr]]$`Housing Cost Class`==housing[[2]])
if(length(housing_class)==0){
housing_class <- 1
}
cost_housing <- nodes[[groupnr]]$`Number of Individuals` * housing[[1]][housing_class]
### Subtract costs from earlier
if(nodes[[groupnr]]$`Breeding Type`=="Selection" || nodes[[groupnr]]$`Breeding Type`=="Aging" || nodes[[groupnr]]$`Breeding Type`=="Split"){
prior_node <- which(ids==nodes[[groupnr]]$`origin`)
share_geno <- as.numeric(nodes[[prior_node]]$`Proportion of genotyped individuals`)
cost_geno <- cost_geno - nodes[[groupnr]]$`Number of Individuals` * genotyping_costs * share_geno
if(length(share_geno)==0) share_geno <- 1
if(cost_geno<0) cost_geno <- 0
pheno_class_old <- which(nodes[[prior_node]]$`Phenotyping Class`==phenotyping[[2]])
if(length(pheno_class_old)==0){
pheno_class_old <- 1
}
if(pheno_class_old==pheno_class){
cost_pheno <- 0
}
}
if(n_traits==0){
costdata <- rbind(costdata, c(group_name, group_size, group_time, cost_geno+cost_housing, cost_geno, cost_housing))
} else{
costdata <- rbind(costdata, c(group_name, group_size, group_time, cost_geno+cost_housing+cost_pheno, cost_geno, cost_pheno, cost_housing))
}
}
}
}
costdata[,-(1:3)] <- round(as.numeric(costdata[,-(1:3)]) * interest_rate^(as.numeric(costdata[,3])), digits=2)
if(n_traits==0){
colnames(costdata) <- c("Cohort name", "Nr. of individuals", "Time-point", "Total costs", "Cost genotyping", "Cost housing")
} else{
colnames(costdata) <- c("Cohort name", "Nr. of individuals", "Time-point", "Total costs", "Cost genotyping", "Cost phenotyping", "Cost housing")
}
if(FALSE){
as.data.frame(costdata)
time_point_plot <- unique(sort(time.point.list))
cost_plot <- numeric(length(time_point_plot))
for(index in 1:length(time_point_plot)){
cost_plot[index] <- sum(as.numeric(costdata[costdata[,3]==time_point_plot[index],4]))
}
plot(time_point_plot, cost_plot, main="Cost - overview", ylab="cost in Euro", xlab="time point", ylim=c(0, max(cost_plot)))
barplot(cost_plot, names=time_point_plot, ylab="cost in Euro", xlab="time point")
cost_plot_sex <- matrix(0, ncol=length(time_point_plot), nrow=2)
for(index in 1:length(time.point.list)){
index2 <- which(costdata[index,1]==ids)
sex <- as.numeric(nodes[[index2]]$"Sex"=="Female") + 1
index3 <- which(time_point_plot==costdata[index,3])
cost_plot_sex[sex,index3] <- cost_plot_sex[sex,index3] + as.numeric(costdata[index,4])
}
barplot(cost_plot_sex, names=time_point_plot, ylab="cost in Euro", xlab="time point", col=c("blue", "red"))
ntype <- length(unique(ids_type))
type_list <- unique(ids_type)
cost_plot_ntype <- matrix(0, ncol = length(time_point_plot), nrow=ntype)
for(index in 1:length(time.point.list)){
index2 <- which(costdata[index,1]==ids)
ntype1 <- which(type_list == ids_type[index2])
index3 <- which(time_point_plot==costdata[index,3])
cost_plot_ntype[ntype1,index3] <- cost_plot_ntype[ntype1,index3] + as.numeric(costdata[index,4])
}
keep1 <- which(rowSums(cost_plot_ntype)>0)
cost_plot_ntype <- cost_plot_ntype[keep1,]
type_list <- type_list[keep1]
qual_col_pals = RColorBrewer::brewer.pal.info[RColorBrewer::brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(RColorBrewer::brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
used_color <- sample(col_vector, length(type_list))
cost_plot_type <- matrix(0, ncol=length(time_point_plot), nrow=3)
for(index in 1:length(time_point_plot)){
cost_plot_type[1,index] <- sum(as.numeric(costdata[costdata[,3]==time_point_plot[index],5]))
cost_plot_type[2,index] <- sum(as.numeric(costdata[costdata[,3]==time_point_plot[index],6]))
cost_plot_type[3,index] <- sum(as.numeric(costdata[costdata[,3]==time_point_plot[index],7]))
}
barplot(cost_plot_ntype, names=time_point_plot,
ylab="cost in Euro", xlab="time point", col=used_color,
ylim=c(0, max(cost_plot_type)*1.5))
X11()
w <- 20
a <- barplot(cost_plot_ntype, names=time_point_plot, width=w, space=(c(0,diff(time_point_plot))-w)/w,
ylab="cost in Euro", xlab="time point", col=used_color,
ylim=c(0, max(cost_plot_ntype)*2),
axes=FALSE, axisnames = FALSE)
axis(1, at=seq(0,600, by=100)-20, label=seq(0,600, by=100) )
axis(2)
legend("topleft", type_list, col=used_color, lty=0, pch=15, ncol = 5)
w <- 20
order <- c(17,2,1,18,4,3,11,5,6,10,7,8,9, 12,13,15,14,16)
palette = c(
col=RColorBrewer::brewer.pal(5,"Greens")[-(c(1,5))],
col=RColorBrewer::brewer.pal(7,"Blues")[-(c(1,7))],
col=RColorBrewer::brewer.pal(5,"Purples")[-(c(1,5))],
col=RColorBrewer::brewer.pal(5,"Greys")[-(c(1,5))],
col=RColorBrewer::brewer.pal(7,"Reds")[-(c(1,7))]
)
used_color <- palette[c(3,2,1,
6,5,4,
19,10,9,
8,13,12,
7, 18,17,16,15,14)]
options("scipen"=100)
a <- barplot(cost_plot_ntype[order,], names=time_point_plot, width=w, space=(c(0,diff(time_point_plot))-w)/w,
ylab="cost in Euro", xlab="time point", col=used_color,
ylim=c(0, max(cost_plot_ntype)*2),
axes=FALSE, axisnames = FALSE)
axis(1, at=seq(0,600, by=100)-20, label=seq(0,600, by=100) )
axis(2)
# 16 = 11
type_list
legend("topleft", type_list[order][-18], col=used_color, lty=0, pch=15, ncol = 6, cex=1.5)
w <- 20
barplot(cost_plot_type, names=time_point_plot, ylab="cost in Euro", xlab="time point",
width=w, space=(c(0,diff(time_point_plot))-w)/w,
col=c("red", "blue", "green"), ylim=c(0, max(cost_plot_type)*1.35),
axes=FALSE, axisnames = FALSE)
axis(1, at=seq(0,600, by=100)-20, label=seq(0,600, by=100) )
axis(2)
legend("topleft", c("Genotyping", "Phenotyping", "Housing"), col=c("red", "blue", "green"), lty=0, pch=15)
barplot(cost_plot_type, names=time_point_plot, ylab="cost in Euro", xlab="time point",
col=c("red", "blue", "green"), ylim=c(0, max(cost_plot_type)*1.35))
legend("topleft", c("Genotyping", "Phenotyping", "Housing"), col=c("red", "blue", "green"), lty=0, pch=15)
plot(time_point_plot, cumsum(cost_plot_type[1,]), type="l", col="red",
ylab="cost in Euro", xlab="time point", lwd=2, ylim=c(0, sum(cost_plot_type)))
lines(time_point_plot, cumsum(cost_plot_type[2,]), col="blue", lwd=2)
lines(time_point_plot, cumsum(cost_plot_type[3,]), col = "green", lwd=2)
legend("topleft", c("Genotyping", "Phenotyping", "Housing"), col=c("red", "blue", "green"), lty=1, lwd=2)
write.csv(file="C:/Users/pook/Desktop/Cost_overview.csv", costdata, row.names = FALSE, quote=FALSE)
}
}
# Clean up everything we do not need anymore
if(TRUE & !skip.population){
rm(dataset)
rm(map)
}
############## Actual simulations ########################
{
if(export.population != FALSE && length(export.timepoint)>0 && export.timepoint!=""){
export.gen <- sum(generation_times<export.timepoint) +1
}
if(verbose) {
if(length(export.gen)>0){
if(export.gen==1){
cat("Exported population contains all cohorts till:", ids[founder])
} else{
cat("Exported population contains all cohorts till:", generation_group[[export.gen-1]])
}
}
}
expected_time <- NULL
if(skip.population || time.check){
## Estimate computing time:
if(n_traits>0){
pop_check <- creating.diploid(nindi = 2, nsnp = 100, chromosome.length = chromo.length,
n.additive = as.numeric(trait_matrix[,6]),
n.dominant = as.numeric(trait_matrix[,9]),
n.qualitative = as.numeric(trait_matrix[,10]),
n.quantitative = as.numeric(trait_matrix[,11]),
verbose = FALSE)
pop_check <- breeding.diploid(pop_check, breeding.size = 1000, verbose =FALSE)
pop_check <- breeding.diploid(pop_check, bve=TRUE, new.bv.observation = "all", verbose = FALSE)
pop_check <- breeding.diploid(pop_check, breeding.size = 1000, copy.individual = TRUE, verbose = FALSE)
time_bve1000 <- pop_check$info$comp.times.bve[2,10]
time_gen1000 <- pop_check$info$comp.times.generation[1,4]
time_copy1000 <- pop_check$info$comp.times.generation[3,4]
} else{
pop_check <- creating.diploid(nindi = 2, nsnp = 100, chromosome.length = chromo.length,
verbose = FALSE)
pop_check <- breeding.diploid(pop_check, breeding.size = 1000, verbose =FALSE)
pop_check <- breeding.diploid(pop_check)
pop_check <- breeding.diploid(pop_check, breeding.size = 1000, copy.individual = TRUE, verbose = FALSE)
time_bve1000 <- 0.25 #pop_check$info$comp.times.bve[2,10]
time_gen1000 <- pop_check$info$comp.times.generation[1,4]
time_copy1000 <- pop_check$info$comp.times.generation[3,4]
}
expected_time <- matrix(0, nrow=length(ids), ncol=4)
index <- 1 + length(founder)
expected_time[1:length(founder),1] <- ids[founder]
for(generation in 1:(length(generation_group)) +1){
for(group in generation_group[[generation-1]]){
expected_time[index,1] <- group
groupnr <- which(ids==group)
if(length(nodes[[groupnr]]$'Manuel selected cohorts')>0){
takes <- c(nodes[[groupnr]]$'Manuel selected cohorts')
n_bve <- sum(as.numeric(costdata[duplicated(c(takes, costdata[,1]))[-(1:length(takes))],2]))
expected_time[index,2] <- max(time_bve1000 * (n_bve/1000)^3, time_bve1000*0.1)
}
if(nodes[[groupnr]]$`Breeding Type` %in% c("Aging", "Selection", "Split", "Cloning", "Combine")){
expected_time[index,3] <- max(time_copy1000 * as.numeric(costdata[which(costdata[,1]==group),2]) / 1000, time_copy1000*0.1)
} else{
expected_time[index,3] <- max(time_gen1000 * as.numeric(costdata[which(costdata[,1]==group),2]) / 1000, time_gen1000*0.1)
}
index <- index + 1
}
}
expected_time[,4] <- as.numeric(expected_time[,2]) + as.numeric(expected_time[,3])
t1 <- round(sum(as.numeric(expected_time[,2])) / 60, digits=1)
t2 <- round(sum(as.numeric(expected_time[,3])) / 60, digits=1)
t <- round(sum(as.numeric(expected_time[,4])) / 60, digits=1)
if(verbose){
if(t > 200){
cat(paste0("Your simulation is expected to take ", round(t/60, digits=2), " hours on this operating system!!!\n" ))
cat(paste0(round(t1/60, digits=2), " hours are due to breeding value estimation.\n" ))
cat(paste0(round(t2/60, digits=2), " hours are due to generation of new individuals.\n" ))
} else{
cat(paste0("Your simulation is expected to take ", t, "minutes on this operating system.\n" ))
cat(paste0(t1, " minutes are due to breeding value estimation.\n" ))
cat(paste0(t2, " minutes are due to generation of new individuals.\n" ))
}
}
rbind(c("TOTAL", sum(as.numeric(expected_time[,2])), sum(as.numeric(expected_time[,3])),
sum(as.numeric(expected_time[,4]))), expected_time)
colnames(expected_time) <- c("Cohort name", "BVEtime", "Gentime", "Totaltime")
}
if(!skip.population){
if(time.check){
if(t*60 > time.max){
stop("Expected run-time exceeds your limits!\n Simulation stopped!!")
}
}
# Derive Founders that are alive
alive_cohorts <- population$info$cohorts[,1]
alive_numbers <- as.numeric(population$info$cohorts[,3])
alive_numbers[alive_numbers==0] <- as.numeric(population$info$cohorts[alive_numbers==0,4])
death_to <- matrix(0, ncol=nrow(culling_reason), nrow=length(alive_cohorts))
if(length(fixed.generation.order)>0){
generation_times <- 1:length(fixed.generation.order)
generation_group <- list()
for(index in 1:length(fixed.generation.order)){
generation_group[[index]] <- fixed.generation.order[index]
}
}
generation_times <- round(generation_times, digits = 4)
for(generation in 1:(length(generation_group)+1) +1){
if(verbose) cat(paste0("Start simulation of generation: ", generation," (time point: ", c(generation_times, Inf)[generation-1], ")\n"))
t1 <- as.numeric(Sys.time())
if(generation != (length(generation_group)+2)){
if(length(geninfo$'trait_rescaling') > 0 && sum(generation_times[generation-1]==as.numeric(unlist(strsplit(geninfo$'trait_rescaling', split=";"))))>0){
population <- bv.standardization(population, mean.target = as.numeric(trait_matrix[,3]),
var.target = as.numeric(trait_matrix[,4])^2,
gen = nrow(population$info$size),
adapt.bve = TRUE,
adapt.pheno = TRUE)
population <- breeding.diploid(population, verbose=verbose)
population <- breeding.diploid(population, sigma.e.gen = nrow(population$info$size),
heritability = as.numeric(trait_matrix[,5]), verbose=verbose)
}
for(group in generation_group[[generation-1]]){
{
if(sum((as.numeric(population$info$cohorts[,3]) + as.numeric(population$info$cohorts[,4]))==0)>0){stop("SF")}
groupnr <- which(ids==group)
sex <- as.numeric(nodes[[groupnr]]$'Sex'=="Female") + 1
breeding.size <- as.numeric(nodes[[groupnr]]$'Number of Individuals') * c(sex==1, sex==2)
name_cohort = nodes[[groupnr]]$id
both_mode <- generation_both[[generation-1]][which(generation_group[[generation-1]]==group)]
if(both_mode==1){
second <- gsub("_M_", "_F_",group)
if(substr(group, start=nchar(group)-1, stop= nchar(group))=="_M"){
second <- paste0(substr(group, start=1, stop = nchar(group)-2), "_F")
}
groupnr2 <- which(ids==second)
breeding.size[2] <- as.numeric(nodes[[groupnr2]]$'Number of Individuals')
name_cohort <- gsub("_M_", "_",group)
if(substr(name_cohort, start=nchar(name_cohort)-1, stop= nchar(name_cohort))=="_M"){
name_cohort <- substr(name_cohort, start=1, stop = nchar(name_cohort)-2)
}
}
if(both_mode==2){
next
}
involved_cohorts <- nodes[[groupnr]]$origin
cohort_data <- population$info$cohorts[involved_cohorts,,drop=FALSE]
sex_cohorts <- (as.numeric(cohort_data[,3])==0) +1
selection.size <- selection.size.max <- c(sum(as.numeric(cohort_data[,3])), sum(as.numeric(cohort_data[,4])))
if(selection.size[1]>0){
selection.size.max[1] <- sum(get.class(population, cohorts=cohort_data[as.numeric(cohort_data[,3])>0,1])!=(-1))
}
if(selection.size[2]>0){
selection.size.max[2] <-sum(get.class(population, cohorts=cohort_data[as.numeric(cohort_data[,4])>0,1])!=(-1))
}
if(nodes[[groupnr]]$'Breeding Type'=="Split"){
split_nr <- which(nodes[[groupnr]]$origin==to_split)
if(selection.size[1]>0){
selection.size.max[1] <- sum(get.class(population, cohorts=cohort_data[as.numeric(cohort_data[,3])>0,1])[split_info[[split_nr]]]!=(-1), na.rm = TRUE)
}
if(selection.size[2]>0){
selection.size.max[2] <-sum(get.class(population, cohorts=cohort_data[as.numeric(cohort_data[,4])>0,1])[split_info[[split_nr]]]!=(-1), na.rm=TRUE)
}
}
if(selection.size[1]>selection.size.max[1] || selection.size[2]>selection.size.max[2] ){
if(verbose) cat(paste0("Less individuals available than designed for cohort: ", group,".\n"))
if(verbose) cat(paste0(selection.size.max[1], " male individuals & ", selection.size.max[2], " female individuals.\n"))
selection.size[selection.size>selection.size.max] <- selection.size.max[selection.size>selection.size.max]
}
if(nodes[[groupnr]]$'Breeding Type'=="Selection" ||
nodes[[groupnr]]$'Breeding Type'=="Aging" ||
nodes[[groupnr]]$'Breeding Type'=="Split" ||
nodes[[groupnr]]$'Breeding Type'=="Combine"){
for(temp1 in 1:2){
if(breeding.size[temp1]>selection.size[temp1]){
breeding.size[temp1] <- selection.size[temp1]
if(verbose) cat(paste0("Reduce size of cohort ", group," to ", breeding.size[temp1],".\n"))
if(sum(group==to_split)>0){
split_info[[which(group==to_split)]] <- 1:breeding.size[temp1]
}
}
}
}
#
if(nodes[[groupnr]]$'Breeding Type'=="Reproduction" ||
nodes[[groupnr]]$'Breeding Type'=="Selfing" ||
nodes[[groupnr]]$'Breeding Type'=="DH-Production" ||
nodes[[groupnr]]$'Breeding Type'=="Cloning"
){
if(length(nodes[[groupnr]]$max_offspring)>0 && nodes[[groupnr]]$max_offspring[1] < Inf && sum(breeding.size)>(selection.size[1] *nodes[[groupnr]]$max_offspring[1])){
breeding.size[breeding.size>0] <- selection.size[1] * nodes[[groupnr]]$max_offspring[1]
if(verbose) cat(paste0("Reduce size of cohort ", group," to ", sum(breeding.size),".\n"))
}
if(length(nodes[[groupnr]]$max_offspring)>0 && nodes[[groupnr]]$max_offspring[2] < Inf && sum(breeding.size)>(selection.size[2] *nodes[[groupnr]]$max_offspring[2])){
breeding.size[breeding.size>0] <- selection.size[2] * nodes[[groupnr]]$max_offspring[2]
if(verbose) cat(paste0("Reduce size of cohort ", group," to ", sum(breeding.size),".\n"))
}
}
share.genotyped <- as.numeric(nodes[[groupnr]]$`Proportion of genotyped individuals`)
cohorts.m <- involved_cohorts[sex_cohorts==1]
cohorts.f <- involved_cohorts[sex_cohorts==2]
involved_groups <- cbind(as.numeric(cohort_data[,2]), sex_cohorts)
# Derive time.point
time.point <- 0
origins <- nodes[[which(ids==group)]]$origin
time.point <- generation_times[generation-1]
bve.database <- NULL
bve_exe <- TRUE
bve_insert_cohorts = c(cohorts.m, cohorts.f)
bve.breeding.type <- nodes[[groupnr]]$`Breeding Type`=="Selection" || nodes[[groupnr]]$`Breeding Type`=="Aging" || nodes[[groupnr]]$`Breeding Type`=="Split"
if(length(nodes[[groupnr]]$'Cohorts used in BVE') || bve.breeding.type){
if(length(nodes[[groupnr]]$'Cohorts used in BVE')==0 || nodes[[groupnr]]$'Cohorts used in BVE'=="Only this cohort" || nodes[[groupnr]]$`Selection Type`=="Pseudo-BVE"){
bve.database <- involved_groups[,1:2, drop=FALSE]
} else if(nodes[[groupnr]]$'Cohorts used in BVE'=="Last 2 Generations"){
bve.database <- get.database(population, gen=max(1,generation-2):(generation-1))
if(generation_bv_size[[generation-1]][1]>0){
bve.database <- rbind(bve.database, c(generation,1,1,generation_bv_size[[generation-1]][1]))
}
if(generation_bv_size[[generation-1]][2]>0){
bve.database <- rbind(bve.database, c(generation,2,1,generation_bv_size[[generation-1]][2]))
}
if(verbose){
cat("Last 2 generations for ", group, "includes:\n")
cat(population$info$cohorts[as.numeric(population$info$cohorts[,2]) >= (generation-2),1])
cat("\n")
}
} else if(nodes[[groupnr]]$'Cohorts used in BVE'=="Last 3 Generations"){
bve.database <- get.database(population, gen=max(1,generation-3):(generation-1))
if(generation_bv_size[[generation-1]][1]>0){
bve.database <- rbind(bve.database, c(generation,1,1,generation_bv_size[[generation-1]][1]))
}
if(generation_bv_size[[generation-1]][2]>0){
bve.database <- rbind(bve.database, c(generation,2,1,generation_bv_size[[generation-1]][2]))
}
if(verbose){
cat("Last 3 generations for ", group, "includes:\n")
cat(population$info$cohorts[as.numeric(population$info$cohorts[,2]) >= (generation-3),1])
cat("\n")
}
} else if(nodes[[groupnr]]$'Cohorts used in BVE'=="Last Generation"){
bve.database <- get.database(population, gen=max(1,generation-1):(generation-1))
if(generation_bv_size[[generation-1]][1]>0){
bve.database <- rbind(bve.database, c(generation,1,1,generation_bv_size[[generation-1]][1]))
}
if(generation_bv_size[[generation-1]][2]>0){
bve.database <- rbind(bve.database, c(generation,2,1,generation_bv_size[[generation-1]][2]))
}
if(verbose){
cat("Last generation for ", group, "includes:\n")
cat(population$info$cohorts[as.numeric(population$info$cohorts[,2]) >= (generation-1),1])
cat("\n")
}
} else if(nodes[[groupnr]]$'Cohorts used in BVE'=="All"){
bve.database <- get.database(population, gen=1:(generation-1))
} else if(nodes[[groupnr]]$'Cohorts used in BVE'=="Manual select"){
bve.database <- get.database(population, cohorts=unique(nodes[[groupnr]]$'Manuel selected cohorts'))
bve.database <- bve.database[!is.na(bve.database[,1]),, drop=FALSE]
bv_cohort_index <- which(bv_cohort==group)
bv_cohort_class_index <- bv_cohort_class[bv_cohort_index]
bve_insert_cohorts <- bv_class_member[[bv_cohort_class_index]]
if(class_executed[bv_cohort_class_index]){
bve_exe <- FALSE
} else{
class_executed[bv_cohort_class_index] <- TRUE
}
}
}
if(bve.breeding.type){
activemmreml <-activesommer <- activemultisommer <- activbglr <- activerrblup <- FALSE
bve_solve <- "exact"
singlestep.active <- FALSE
depth <- 0
parent_average <- FALSE
grandparent_average <- FALSE
mean_between <- NULL
phenotype.bv <- FALSE
pseudo_bve <- FALSE
computeA <- "vanRaden"
input_phenotype <- "own"
pseudo_acc <- NULL
bglrmodel <- "RKHS"
bvemas <- FALSE
masmarker <- 0
threshold <- NULL
threshold_sign <- ">"
if(length(nodes[[groupnr]]$skip)>0 && nodes[[groupnr]]$skip=="Yes"){
skip <- "zero"
} else{
skip <- NULL
}
if(length(nodes[[groupnr]]$'threshold')>0){
threshold <- as.numeric(nodes[[groupnr]]$'threshold')
}
if(length(nodes[[groupnr]]$'bve_solve')>0){
if(nodes[[groupnr]]$'bve_solve'=="PCG"){
bve_solve <- "pcg"
}
}
if(length(nodes[[groupnr]]$'threshold_sign')>0){
threshold_sign <- nodes[[groupnr]]$'threshold_sign'
}
if(length(nodes[[groupnr]]$'Selection Type')==0){
if(verbose) cat("No selection type selected in some edges. Assume selection type 'Random'")
bve <- FALSE
selection <- "random"
}else if(nodes[[groupnr]]$'Selection Type'=="Parent_Mean"){
bve <- FALSE
parent_average <- TRUE
selection <- "function"
if(nodes[[groupnr]]$'Selection Type Mean'=="BVE"){
mean_between <- "bve"
} else if(nodes[[groupnr]]$'Selection Type Mean'=="Phenotype"){
mean_between <- "pheno"
} else if(nodes[[groupnr]]$'Selection Type Mean'=="BV"){
mean_between <- "bv"
} else {
mean_between <- "bvepheno"
}
} else if(nodes[[groupnr]]$'Selection Type'=="Grandparent_Mean"){
bve <- FALSE
grandparent_average <- TRUE
selection <- "function"
if(nodes[[groupnr]]$'Selection Type Mean'=="BVE"){
mean_between <- "bve"
} else if(nodes[[groupnr]]$'Selection Type Mean'=="Phenotype"){
mean_between <- "pheno"
} else if(nodes[[groupnr]]$'Selection Type Mean'=="BV"){
mean_between <- "bv"
} else {
mean_between <- "bvepheno"
}
} else if(nodes[[groupnr]]$'Selection Type'=="Phenotypic"){
bve <- FALSE
selection <- "function"
phenotype.bv <- TRUE
} else if(nodes[[groupnr]]$'Selection Type' == "Pseudo-BVE"){
bve <- pseudo_bve <- TRUE
pseudo_acc <- nodes[[groupnr]]$'PseudoAcc'
selection <- "function"
} else if(nodes[[groupnr]]$'Selection Type'=="BVE"){
bve <- TRUE
selection <- "function"
if(nodes[[groupnr]]$'Relationship Matrix'=="Pedigree"){
computeA <- "kinship"
depth <- as.numeric(nodes[[groupnr]]$'Depth of Pedigree')
} else if(nodes[[groupnr]]$'Relationship Matrix'=="Single Step"){
depth <- as.numeric(nodes[[groupnr]]$'Depth of Pedigree')
singlestep.active <- TRUE
}
if(nodes[[groupnr]]$'BVE Method'=="REML-GBLUP" || nodes[[groupnr]]$'BVE Method'=="REML-GBLUP (EMMREML)"){
activemmreml <- TRUE
} else if(nodes[[groupnr]]$'BVE Method' =="Last BVE"){
bve <- FALSE
} else if(nodes[[groupnr]]$'BVE Method'=="REML-GBLUP (sommer)") {
activesommer <- TRUE
} else if(nodes[[groupnr]]$'BVE Method'=="Multi-trait REML-GBLUP (sommer)") {
activemultisommer <- TRUE
} else if(nodes[[groupnr]]$'BVE Method'=="REML-GBLUP (rrBLUP)"){
activerrblup <- TRUE
} else if(nodes[[groupnr]]$'BVE Method'=="BayesA (BGLR)"){
activbglr <- TRUE
bglrmodel <- "BayesA"
} else if(nodes[[groupnr]]$'BVE Method'=="BayesB (BGLR)"){
activbglr <- TRUE
bglrmodel <- "BayesB"
} else if(nodes[[groupnr]]$'BVE Method'=="BayesC (BGLR)"){
activbglr <- TRUE
bglrmodel <- "BayesC"
} else if(nodes[[groupnr]]$'BVE Method'=="RKHS (BGLR)"){
activbglr <- TRUE
bglrmodel <- "RKHS"
} else if(nodes[[groupnr]]$'BVE Method'=="BL (BGLR)"){
activbglr <- TRUE
bglrmodel <- "BL"
} else if(nodes[[groupnr]]$'BVE Method'=="BRR (BGLR)"){
activbglr <- TRUE
bglrmodel <- "BRR"
} else if(nodes[[groupnr]]$'BVE Method'=="Marker assisted selection (lm)"){
bvemas <-TRUE
if(length(nodes[[groupnr]]$'MAS_marker')==1){
masmarker <- as.numeric(nodes[[groupnr]]$'MAS_marker')
} else{
masmarker <- 10
}
}
} else if(nodes[[groupnr]]$'Selection Type'=="Random"){
bve <- FALSE
selection <- "random"
} else{
if(verbose) cat("No selection type selected in some edges. Assume selection type 'Random'")
bve <- FALSE
selection <- "random"
}
if(length(involved_cohorts)>1){
stop("Only one cohort to select from allowed in selection - check for error")
} else{
add.observation <- pheno_index[which(pheno_index_name==nodes[[groupnr]]$'Phenotyping Class'),]
add.observation[add.observation<0] <- 0
}
reduced.selection.panel.m <- NULL
reduced.selection.panel.f <- NULL
if(nodes[[groupnr]]$'Breeding Type'=="Selection"){
creating.type <- 1
} else if(nodes[[groupnr]]$'Breeding Type'=="Aging"){
creating.type <- 8
} else if(nodes[[groupnr]]$'Breeding Type'=="Split"){
creating.type <- 9
split_nr <- which(nodes[[groupnr]]$origin==to_split)
reduced.selection.panel.m <- split_info[[split_nr]]
reduced.selection.panel.f <- split_info[[split_nr]]
}
if(nodes[[groupnr]]$'Use Offspring for BVE'=="Yes"){
offspring.bve.parents.database <- get.database(population, cohorts=c(cohorts.m, cohorts.f))
if(nodes[[groupnr]]$phenotype_used=="Avg. offspring phenotye"){
input_phenotype <- "off"
} else if(nodes[[groupnr]]$phenotype_used=="Mean own/offspring phenotype"){
input_phenotype <- "mean"
} else if(nodes[[groupnr]]$phenotype_used=="Weighted own/offspring phenotype"){
input_phenotype <- "weighted"
}
} else{
offspring.bve.parents.database <- NULL
}
population <- breeding.diploid(population, breeding.size=breeding.size,
bve=(bve&bve_exe),
bve.solve = bve_solve,
computation.A = computeA,
bve.pseudo = pseudo_bve,
bve.pseudo.accuracy = pseudo_acc,
offspring.bve.parents.database=offspring.bve.parents.database,
BGLR.bve = activbglr,
BGLR.model = bglrmodel,
emmreml.bve = activemmreml,
rrblup.bve = activerrblup,
sommer.bve = activesommer,
sommer.multi.bve = activemultisommer,
mas.bve = bvemas,
mas.number = masmarker,
selection.size= breeding.size,
copy.individual = TRUE,
added.genotyped = nodes[[groupnr]]$`Proportion of added genotypes`,
max.offspring = c(1,1),
heritability = heritability,
sigma.e.database = cbind(1,(1:2)[population$info$size[1,]>0]),
new.bv.child="addobs",
selection.m = selection,
selection.f = selection,
phenotype.bv = phenotype.bv,
add.gen = generation,
bve.database = bve.database,
selfing.mating=TRUE,
selfing.sex=(sex-1),
selection.m.cohorts = cohorts.m,
selection.f.cohorts = cohorts.f,
new.class = new_mig[sex],
selection.m.miesenberger = if(n_traits>0){selection_index_miesenberger[which(selection_index_name==nodes[[groupnr]]$'Selection Index')]} else{FALSE} ,
multiple.bve.scale.m = selection_index_miesenberger_wscaling[which(selection_index_name==nodes[[groupnr]]$'Selection Index')],
n.observation = add.observation,
remove.effect.position = remove.effect.position,
name.cohort = name_cohort,
time.point = time.point,
creating.type = creating.type,
depth.pedigree = depth,
store.breeding.totals = TRUE,
reduced.selection.panel.m = reduced.selection.panel.m,
reduced.selection.panel.f = reduced.selection.panel.f,
bve.cohorts = bve_insert_cohorts,
bve.insert.cohorts = bve_insert_cohorts,
display.progress=progress.bars,
singlestep.active=singlestep.active,
share.genotyped = share.genotyped,
multiple.bve.weights.m = selection_index[which(selection_index_name==nodes[[groupnr]]$'Selection Index'),],
verbose=verbose,
bve.parent.mean=parent_average,
bve.grandparent.mean=grandparent_average,
bve.mean.between=mean_between,
threshold.selection = threshold,
threshold.sign = threshold_sign,
input.phenotype = input_phenotype,
bve.ignore.traits = skip,
miraculix.chol = miraculix.chol
)
if(nodes[[groupnr]]$'Breeding Type'=="Split"){
split_info[[split_nr]] <- sort(setdiff(split_info[[split_nr]], split_info[[split_nr]][population$info$breeding.totals[[length(population$info$breeding.totals)]][[7]][[sex]]]))
}
} else if(nodes[[groupnr]]$'Breeding Type'=="Reproduction"){
generation.sex <- as.numeric(selection.size[2]>0)- 0.5 * as.numeric((selection.size[1]>0)*(selection.size[2]>0))
if(generation.sex==0 || generation.sex==1){
same.sex.activ <- TRUE
same.sex.sex <- generation.sex
} else{
same.sex.activ <- FALSE
same.sex.sex <- generation.sex
}
half_sib <- nodes[[groupnr]]$half_sib
full_sib <- nodes[[groupnr]]$full_sib
population <- breeding.diploid(population, breeding.size=breeding.size,
selection.size= selection.size,
heritability = heritability,
sigma.e.database = cbind(1,(1:2)[population$info$size[1,]>0]),
new.bv.child="obs",
selection.m = if(sum(nodes[[groupnr]]$selection_ratio==1)==2){"random"} else{"function"},
add.gen = generation,
bve.database = bve.database,
selection.m.cohorts = cohorts.m,
selection.f.cohorts = cohorts.f,
avoid.mating.halfsib = half_sib,
avoid.mating.fullsib = full_sib,
n.observation = pheno_index[which(pheno_index_name==nodes[[groupnr]]$'Phenotyping Class'),],
new.class = new_mig[sex],
same.sex.activ = same.sex.activ,
same.sex.sex = same.sex.sex,
same.sex.selfing = FALSE,
name.cohort = name_cohort,
time.point = time.point,
creating.type = 2,
store.breeding.totals = TRUE,
display.progress=progress.bars,
share.genotyped = share.genotyped,
added.genotyped = nodes[[groupnr]]$`Proportion of added genotypes`,
ogc = nodes[[groupnr]]$OGC,
ogc.target = nodes[[groupnr]]$ogc_target,
relationship.matrix.ogc = nodes[[groupnr]]$ogc_relation,
ogc.ub.BV = nodes[[groupnr]]$constrain[[1]],
ogc.eq.BV = nodes[[groupnr]]$constrain[[2]],
ogc.lb.BV = nodes[[groupnr]]$constrain[[3]],
ogc.ub.sKin = nodes[[groupnr]]$constrain[[4]],
ogc.uniform = nodes[[groupnr]]$constrain[[5]],
ogc.lb.BV.increase = nodes[[groupnr]]$constrain[[6]],
ogc.ub.sKin.increase = nodes[[groupnr]]$constrain[[7]],
repeat.mating = repeat.mating * nodes[[groupnr]]$repeat_mating,
repeat.mating.overwrite = FALSE,
max.offspring = nodes[[groupnr]]$max_offspring,
max.mating.pair = min(nodes[[groupnr]]$max_offspring_pair),
best.selection.ratio.m = nodes[[groupnr]]$selection_ratio[1],
best.selection.ratio.f = nodes[[groupnr]]$selection_ratio[2],
best.selection.criteria.m = nodes[[groupnr]]$selection_ratio_type[1],
best.selection.criteria.f = nodes[[groupnr]]$selection_ratio_type[2],
multiple.bve.weights.m = selection_index[max(1,which(selection_index_name==nodes[[groupnr]]$selection_ratio_index[1])),],
multiple.bve.weights.f = selection_index[max(1,which(selection_index_name==nodes[[groupnr]]$selection_ratio_index[2])),],
multiple.bve.scale.m = selection_index_miesenberger_wscaling[max(1,which(selection_index_name==nodes[[groupnr]]$selection_ratio_index[1]))],
multiple.bve.scale.f = selection_index_miesenberger_wscaling[max(1,which(selection_index_name==nodes[[groupnr]]$selection_ratio_index[2]))],
verbose=verbose,
miraculix.chol = miraculix.chol
)
} else if(nodes[[groupnr]]$'Breeding Type'=="Selfing"){
selfing.sex <- as.numeric(selection.size[2]>0)- 0.5 * as.numeric((selection.size[1]>0)*(selection.size[2]>0))
population <- breeding.diploid(population, breeding.size=breeding.size,
selection.size= selection.size,
selfing.mating = TRUE,
selfing.sex = selfing.sex,
heritability = heritability,
sigma.e.database = cbind(1,(1:2)[population$info$size[1,]>0]),
new.bv.child="obs",
selection.m = "random",
add.gen = generation,
bve.database = bve.database,
selection.m.cohorts = cohorts.m,
selection.f.cohorts = cohorts.f,
n.observation = pheno_index[which(pheno_index_name==nodes[[groupnr]]$'Phenotyping Class'),],
new.class = new_mig[sex],
name.cohort = name_cohort,
time.point = time.point,
creating.type = 4,
store.breeding.totals = TRUE,
display.progress=progress.bars,
share.genotyped = share.genotyped,
added.genotyped = nodes[[groupnr]]$`Proportion of added genotypes`,
repeat.mating = repeat.mating* nodes[[groupnr]]$repeat_mating,
repeat.mating.overwrite = FALSE,
max.offspring = nodes[[groupnr]]$max_offspring,
max.mating.pair = min(nodes[[groupnr]]$max_offspring_pair),
verbose=verbose,
miraculix.chol = miraculix.chol)
} else if(nodes[[groupnr]]$'Breeding Type'=="DH-Production"){
dh.sex <- as.numeric(selection.size[2]>0)- 0.5 * as.numeric((selection.size[1]>0)*(selection.size[2]>0))
population <- breeding.diploid(population, breeding.size=breeding.size,
selection.size= selection.size,
dh.mating = TRUE,
dh.sex = dh.sex,
selfing.mating = TRUE,
selfing.sex = dh.sex,
heritability = heritability,
sigma.e.database = cbind(1,(1:2)[population$info$size[1,]>0]),
new.bv.child="obs",
selection.m = "random",
add.gen = generation,
bve.database = bve.database,
selection.m.cohorts = cohorts.m,
selection.f.cohorts = cohorts.f,
n.observation = pheno_index[which(pheno_index_name==nodes[[groupnr]]$'Phenotyping Class'),],
new.class = new_mig[sex],
name.cohort = name_cohort,
time.point = time.point,
creating.type = 5,
store.breeding.totals = TRUE,
display.progress=progress.bars,
share.genotyped = share.genotyped,
added.genotyped = nodes[[groupnr]]$`Proportion of added genotypes`,
repeat.mating = repeat.mating * nodes[[groupnr]]$repeat_mating,
repeat.mating.overwrite = FALSE,
max.offspring = nodes[[groupnr]]$max_offspring,
max.mating.pair = min(nodes[[groupnr]]$max_offspring_pair),
verbose=verbose,
miraculix.chol = miraculix.chol)
} else if(nodes[[groupnr]]$'Breeding Type'=="Recombination"){
population <- breeding.diploid(population, breeding.size=breeding.size,
mutation.rate = nodes[[groupnr]]$mutation,
remutation.rate = nodes[[groupnr]]$remutation,
recombination.rate = nodes[[groupnr]]$recom,
selection.size= selection.size,
heritability = heritability,
sigma.e.database = cbind(1,(1:2)[population$info$size[1,]>0]),
new.bv.child="obs",
selection.m = "random",
add.gen = generation,
bve.database = bve.database,
selection.m.cohorts = cohorts.m,
selection.f.cohorts = cohorts.f,
n.observation = pheno_index[which(pheno_index_name==nodes[[groupnr]]$'Phenotyping Class'),],
new.class = new_mig[sex],
name.cohort = name_cohort,
time.point = time.point,
creating.type = 3,
store.breeding.totals = TRUE,
display.progress=progress.bars,
share.genotyped = share.genotyped,
added.genotyped = nodes[[groupnr]]$`Proportion of added genotypes`,
repeat.mating = repeat.mating * nodes[[groupnr]]$repeat_mating,
repeat.mating.overwrite = FALSE,
max.offspring = nodes[[groupnr]]$max_offspring,
max.mating.pair = min(nodes[[groupnr]]$max_offspring_pair),
verbose=verbose,
miraculix.chol = miraculix.chol)
} else if(nodes[[groupnr]]$'Breeding Type'=="Cloning"){
selfing.sex <- as.numeric(selection.size[2]>0)- 0.5 * as.numeric((selection.size[1]>0)*(selection.size[2]>0))
if(length(involved_cohorts)>1){
stop("Only one cohort to select from allowed in selection - check for error")
} else{
add.observation <- pheno_index[which(pheno_index_name==nodes[[groupnr]]$'Phenotyping Class'),]
add.observation[add.observation<0] <- 0
}
population <- breeding.diploid(population, breeding.size=breeding.size,
selection.size= selection.size,
copy.individual = TRUE,
added.genotyped = nodes[[groupnr]]$`Proportion of added genotypes`,
selfing.mating = TRUE,
selfing.sex = selfing.sex,
heritability = heritability,
sigma.e.database = cbind(1,(1:2)[population$info$size[1,]>0]),
new.bv.child="addobs",
selection.m = "random",
add.gen = generation,
bve.database = bve.database,
selection.m.cohorts = cohorts.m,
selection.f.cohorts = cohorts.f,
n.observation = pheno_index[which(pheno_index_name==nodes[[groupnr]]$'Phenotyping Class'),],
new.class = new_mig[sex],
name.cohort = name_cohort,
time.point = time.point,
creating.type = 6,
store.breeding.totals = TRUE,
display.progress=progress.bars,
share.genotyped = share.genotyped,
repeat.mating.copy = repeat.mating.copy * nodes[[groupnr]]$repeat_mating,
repeat.mating.overwrite = FALSE,
max.offspring = nodes[[groupnr]]$max_offspring,
max.mating.pair = min(nodes[[groupnr]]$max_offspring_pair),
verbose=verbose,
miraculix.chol = miraculix.chol)
} else if(nodes[[groupnr]]$'Breeding Type'=="Combine"){
selfing.sex <- (as.numeric(selection.size[2])>0)- 0.5 * as.numeric((selection.size[1]>0)*(selection.size[2]>0))
if(FALSE){
stop("Only one cohort to select from allowed in selection - check for error")
} else{
add.observation <- pheno_index[which(pheno_index_name==nodes[[groupnr]]$'Phenotyping Class'),]
add.observation[add.observation<0] <- 0
}
population <- breeding.diploid(population, breeding.size=breeding.size,
selection.size= breeding.size,
copy.individual = TRUE,
added.genotyped = nodes[[groupnr]]$`Proportion of added genotypes`,
max.offspring = c(1,1),
new.bv.child="addobs",
selection.m = "random",
selfing.mating=TRUE,
selfing.sex=selfing.sex,
add.gen = generation,
selection.m.cohorts = cohorts.m,
selection.f.cohorts = cohorts.f,
n.observation = add.observation,
new.class = new_mig[sex],
name.cohort = name_cohort,
time.point = time.point,
creating.type = 7,
store.breeding.totals = TRUE,
display.progress=progress.bars,
share.genotyped = share.genotyped,
verbose=verbose,
miraculix.chol = miraculix.chol)
}
if(both_mode==1 & breeding.size[1]>0 & breeding.size[2]>0){
population$info$cohorts[c(-1,0)+nrow(population$info$cohorts),1] <- c(group, second)
rownames(population$info$cohorts)[c(-1,0)+nrow(population$info$cohorts)]<- c(group, second)
}
position[groupnr,] <- c(generation, sex, new_mig[sex], sum(breeding.size))
new_mig[sex] <- new_mig[sex] + 1
if(phenotype_groups[groupnr]==0){
tested <- which(duplicated(c(nodes[[groupnr]]$origin, ids))[-(1:length(nodes[[groupnr]]$origin))])
n_tester_generated[tested] <- n_tester_generated[tested] + 1
}
}
}
}
if(length(geninfo$advanced_history)>0 && geninfo$advanced_history && base.cycle < Inf && generation%%base.cycle == 0){
population <- new.base.generation(population, base.gen = generation,
delete.previous.gen = delete.previous.gen)
}
if(generation != (length(generation_group)+2)){
alive_cohorts <- c(alive_cohorts, generation_group[[generation-1]])
new_cohorts <- population$info$cohorts[(nrow(population$info$cohorts)-length(generation_group[[generation-1]])+1):nrow(population$info$cohorts),,drop=FALSE]
new_numbers <- as.numeric(new_cohorts[,3])
new_numbers[new_numbers==0] <- as.numeric(new_cohorts[new_numbers==0,4])
# Creating-types:
# 0 - Founder
# 1 - Selection
# 2 - Reproduction
# 3 - Recombination
# 4 - Selfing
# 5 - DH-Production
# 6 - Cloning
# 7 - Combine
# 8 - Aging
# 9 - Split
non_copy <- new_cohorts[,9] == 0 | new_cohorts[,9] == 2 | new_cohorts[,9] == 3 | new_cohorts[,9] == 4 | new_cohorts[,9] == 5 | new_cohorts[,9] == 6
new_numbers[!non_copy] <- 0
alive_numbers <- c(alive_numbers, new_numbers)
}
# death_to <- rbind(death_to, matrix(0, nrow=length(generation_group[[generation-1]]), ncol=nrow(culling_reason)))
t2 <- as.numeric(Sys.time())
from <- if(generation==2){0} else{generation_times[generation-2]}
if(generation != (length(generation_group)+2)){
to <- generation_times[generation-1]
} else{
to <- Inf
}
if(verbose & nrow(culling_reason)>0) cat(paste0("Start culling modul for time span ", from ," to ", to,".\n"))
for(cohort_index in alive_cohorts[alive_numbers>0]){
individual_time <- get.age.point(population, cohorts = cohort_index)
age_start <- min(from - individual_time) ## min()
age_end <- max(to - individual_time) ## max()
sex <- if(population$info$cohorts[cohort_index,4]=="0"){"Male"} else{"Female"}
sex1 <- as.numeric(sex=="Female") + 1
# Time - period, sex
active_culling <- as.numeric(culling_reason[,2]) > age_start & as.numeric(culling_reason[,2]) <= age_end & (culling_reason[,3] == sex | culling_reason[,3]=="Both")
active_cohort <- which(population$info$cohorts[,1]==cohort_index)
for(culling_index in which(active_culling)){
age_start_single <- from - individual_time
age_end_single <- to - individual_time
active_single <- (as.numeric(culling_reason[culling_index,2]) > age_start_single) & (as.numeric(culling_reason[culling_index,2]) <= age_end_single) & rep(culling_reason[culling_index,3] == sex | culling_reason[culling_index,3]=="Both", length(age_end_single))
death_counter <- get.class(population, cohorts=cohort_index)
active_single <- active_single & (death_counter != (-1))
if(sum(active_single)>0){
population <- breeding.diploid(population,
culling.cohort = cohort_index,
culling.time = as.numeric(culling_reason[culling_index,2]),
culling.name = culling_reason[culling_index,1],
culling.bv1 = as.numeric(culling_reason[culling_index,6]),
culling.bv2 = as.numeric(culling_reason[culling_index,8]),
culling.share1 = as.numeric(culling_reason[culling_index,5]),
culling.share2 = as.numeric(culling_reason[culling_index,7]),
culling.index = selection_index[which(selection_index_name==culling_reason[culling_index,4]),],
culling.single = active_single,
culling.all.copy=TRUE,
verbose=verbose
)
if(length(population$info$culling.stats)>=active_cohort){
new_death <- sum(death_counter != get.class(population, cohorts=cohort_index))
} else{
new_death <- 0
}
if(length(new_death)>0){
alive_numbers[active_cohort] <- alive_numbers[active_cohort] - new_death
}
}
}
if(alive_numbers[active_cohort]==0){
if(verbose) cat(paste0("All individuals in cohort ", cohort_index, " are death.\n"))
if(verbose) cat("Reason for death:\n")
for(index in 1:nrow(population$info$culling.stats[[active_cohort]])){
if(verbose) cat(population$info$culling.stats[[active_cohort]][index,])
if(verbose) cat("\n")
}
}
}
t3 <- as.numeric(Sys.time())
if(generation != (length(generation_group)+2)){
if(verbose){
cat("Generated groups: \n")
cat(generation_group[[generation-1]])
cat("\n")
cat("Time spend: ")
cat(round(t3-t1, digits=2))
cat(" seconds.\n")
cat("For generation: ")
cat(round(t2-t1, digits=2))
cat(" seconds.\n")
cat("For culling: ")
cat(round(t3-t2, digits=2))
cat(" seconds.\n")
}
}
if(TRUE){
temp1 <- length(population$breeding)
population <- clean.up(population, gen = temp1)
for(index2 in 3:30){
test <- population$breeding[[temp1]][[index2]]==as.integer(population$breeding[[temp1]][[index2]])
test[is.na(test)] <- 0
test[is.na(population$breeding[[temp1]][[index2]])] <- 1
if(prod(test)==1){
storage.mode(population$breeding[[temp1]][[index2]]) <- "integer"
}
}
}
if(sum(alive_numbers<0)>0){ stop("Some multi-death?!")}
if(export.population!=FALSE && export.gen == generation){
if(verbose) cat("Export population list to ", export.population)
population_demiraculix <- demiraculix(population)
save(file=export.population, list=c("population", "population_demiraculix"))
}
}
}
}
if(!skip.population){
delete_stuff <- rep(FALSE, length(nodes))
for(index in 1:length(nodes)){
delete_stuff[index] <- nodes[[index]]$delete_info
}
if(sum(delete_stuff)>0){
del_database <- get.database(population, cohorts = ids[delete_stuff])
for(index in 1:nrow(del_database)){
for(index2 in del_database[index,3]:del_database[index,4]){
population$breeding[[del_database[index,1]]][[del_database[index,2]]][[index2]] <- "removed"
}
}
}
}
if(log != FALSE){
sink(zz, append = TRUE, type = c("message"))
warnings()
sink(NULL)
sink(NULL, type=c("message"))
}
if(skip.population){
return(list(costdata, expected_time))
} else{
population$info$json <- list(nodes, edges, geninfo, traitinfo, major, housing, phenotyping, ids)
population$info$cost.data <- as.data.frame(costdata)
population$info$expected.time <- as.data.frame(expected_time)
return(population)
}
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.