Replies: 1 comment
-
Change ev to ebv. Only noticed this question now. |
Beta Was this translation helpful? Give feedback.
0 replies
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
-
Please Professor, I would like to ask abouk this error message in R with the package. I have used a genomic model to estimate the breeding value and apply it to the course's scenario eg. But it is like the package does not store the ebv or there is a problem related to my parameters. I wondering what is the reason for the error message: Error in pop@ev[, trait, drop = FALSE] : subscript out of bounds.
Here is the script :
Simulate founder genomes
founderPop = runMacs(nInd = 2000,
nChr = 30,
segSites = 4000,
species = "CATTLE")
Global simulation parameters
SP = SimParam$new(founderPop)
SP$setSexes("yes_sys")
SP$addTraitA(nQtlPerChr = 1000, mean = 158)
SP$addSnpChip(3000)
heritability = 0.3
Base population
pop = newPop(founderPop, simParam = SP)
Phenotype the base population
pop = setPheno(pop, varE=0.3)
#Run GS model and set EBV
ans = RRBLUP(pop, simParam = SP)
pop = setEBV(pop, ans, simParam = SP)
pop@ebv= matrix(rnorm(pop@nInd), nrow = pop@nInd, ncol = 1)
#Evaluate accuracy of the predicted ebv
cor(gv(pop), ebv(pop))
Assign year of birth for the founders
year = 0
founders = setMisc(x = pop,
node = "yearOfBirth",
value = year)
head(getMisc(x = pop, node = "yearOfBirth"))
Initial parent populations
males = selectInd(pop = founders, nInd = 50, use = "ebv", sex = "M")
baseSires2 = males[ 1:10]
baseSires2 = setMisc(x = baseSires2,
node = "yearOfBirth",
value = -1)
baseSires1 = males[11:50]
baseSires1 = setMisc(x = baseSires1,
node = "yearOfBirth",
value = 0)
baseSires = c(baseSires2, baseSires1)
nInd(baseSires)
Inspect year of birth fr the sires
table(unlist(getMisc(x = baseSires, node = "yearOfBirth")))
#Female founder selection
cat("Founder females\n")
(nFemales = sum(founders@sex == "F"))
females = selectInd(pop = founders, nInd = nFemales, use = "ebv", sex = "F")
Here we define how many dams are kept per age group
nDams1 = 500
nDams2 = 250
nDams3 = 150
nDams4 = 75
nDams5 = 25
sum(nDams1, nDams2, nDams3, nDams4, nDams5)
#oldest group of dams
cat("Dams5\n")
(start = 1)
(end = nDams5)
baseDams5 = females[start:end]
baseDams5 = setMisc(x = baseDams5,
node = "yearOfBirth",
value = -4)
nInd(baseDams5)
cat("Dams4\n")
(start = end + 1)
(end = start - 1 + nDams4)
baseDams4 = females[start:end]
baseDams4 = setMisc(x = baseDams4,
node = "yearOfBirth",
value = -3)
nInd(baseDams4)
cat("Dams3\n")
(start = end + 1)
(end = start - 1 + nDams3)
baseDams3 = females[start:end]
baseDams3 = setMisc(x = baseDams3,
node = "yearOfBirth",
value = -2)
nInd(baseDams3)
cat("Dams2\n")
(start = end + 1)
(end = start - 1 + nDams2)
baseDams2 = females[start:end]
baseDams2 = setMisc(x = baseDams2,
node = "yearOfBirth",
value = -1)
nInd(baseDams2)
cat("Dams1\n")
(start = end + 1)
(end = start - 1 + nDams1)
baseDams1 = females[start:end]
baseDams1 = setMisc(x = baseDams1,
node = "yearOfBirth",
value = 0)
nInd(baseDams1)
baseDams = c(baseDams5, baseDams4, baseDams3, baseDams2, baseDams1)
nInd(baseDams)
Year of birth for the dams inspection
table(unlist(getMisc(x = baseDams, node = "yearOfBirth")))
#Data recording
recordData <- function(data = NULL, pop, yearOfUse = NA, scenario = NA) {
popData = data.frame(id = pop@id,
father = pop@father,
mother = pop@mother,
sex = pop@sex,
gv = pop@gv[, "Trait1"],
pheno = pop@pheno[, "Trait1"],
ebv = pop@ebv[, "Trait1"],
yearOfBirth = unlist(getMisc(x = pop, node ="yearOfBirth")),
yearOfUse = yearOfUse,
scenario = scenario)
Manage first instance of calling this function, when data is NULL
if (is.null(data)) {
ret = popData
} else {
ret = rbind(data, popData)
}
return(ret)
}
#dataframes creation
data4AllAnimals = recordData(pop = founders,
scenario = "Founders")
head(data4AllAnimals)
data4AllParents = recordData(pop = c(baseSires, baseDams),
yearOfUse = year,
scenario = "Founders") # year is 0 at this stage in the script
head(data4AllParents)
data4NewParents = recordData(pop = c(baseSires1, baseDams1),
scenario = "Founders")
head(data4NewParents)
Multiple years of the two scenarios
Initiate the sires and dams for the Scenario 1
sires = baseSires
sires1 = baseSires1
dams = baseDams
dams1 = baseDams1
dams2 = baseDams2
dams3 = baseDams3
dams4 = baseDams4
Set the current scenario name
currentScenario = "Scenario_1"
for (year in 1:20) {
cat("Working on the year:", year, "\n")
Generate progeny from current dams and sires
candidates = randCross2(males = sires, females = dams, nCrosses = nInd(dams))
candidates = setMisc(x = candidates, node = "yearOfBirth", value = year)
candidates = setPheno(candidates, h2 = heritability)
Record data for all newborn animals
data4AllAnimals = recordData(data = data4AllAnimals,
pop = candidates,
scenario = currentScenario)
write.csv(data4AllAnimals, file = "data4AllAnimals-ebv1-S1.csv")
Record data for the used sires and dams (young and old)
data4AllParents = recordData(data = data4AllParents,
pop = c(sires, dams),
yearOfUse = year,
scenario = currentScenario)
write.csv(data4AllParents, file = "data4AllParents-ebv1-S1.csv")
Update and select sires
sires2 = selectInd(pop = sires1, nInd = 1, use = "ebv")
sires1 = selectInd(pop = candidates, nInd = 2, use = "ebv", sex = "M")
sires = c(sires2, sires1)
Update and select dams
dams5 = selectInd(pop = dams4, nInd = nDams5, use = "ebv")
dams4 = selectInd(pop = dams3, nInd = nDams4, use = "ebv")
dams3 = selectInd(pop = dams2, nInd = nDams3, use = "ebv")
dams2 = selectInd(pop = dams1, nInd = nDams2, use = "ebv")
dams1 = selectInd(pop = candidates, nInd = nDams1, use = "ebv", sex = "F")
dams = c(dams5, dams4, dams3, dams2, dams1)
Record data for the newly selected sires and dams (just the new ones)
data4NewParents = recordData(data = data4NewParents,
pop = c(sires1, dams1),
scenario = currentScenario)
write.csv(data4NewParents, file = "data4NewParents-ebv1-s1.csv")
}
Error in pop@ev[, trait, drop = FALSE] : subscript out of bounds
Beta Was this translation helpful? Give feedback.
All reactions