Design and Maintenance: Dong Yin, Xuanning
Zhang, Lilin Yin ,Haohao Zhang, and Xiaolei
Liu.
Contributors: Zhenshuang Tang, Jingya Xu,
Xiaohui Yuan, Xinyun Li, and Shuhong Zhao.
If you have any bug reports or questions, please feed back :point_right:here:point_left:.
WE STRONGLY RECOMMEND TO INSTALL SIMER ON Microsoft R Open (https://mran.microsoft.com/download/).
::install_github("xiaolei-lab/SIMER") devtools
After installed successfully, SIMER
can
be loaded by typing
> library(simer)
Typing ?simer
could get the details of all
parameters.
Genotype data should be
Numeric format (m
rows and n columns,
m is the number of SNPs,
n is the number of individuals). Other
genotype data such as PLINK
Binary format (details see
http://zzz.bwh.harvard.edu/plink/data.shtml#bed),
VCF or Hapmap can be
converted to Numeric format using
MVP.Data
function in the rMVP
(https://github.com/xiaolei-lab/rMVP).
genotype.txt
2 | 1 | 0 | 1 | 0 | … | 0 |
1 | 2 | 0 | 1 | 0 | … | 0 |
1 | 1 | 2 | 1 | 0 | … | 0 |
1 | 1 | 0 | 2 | 1 | … | 0 |
0 | 0 | 0 | 0 | 2 | … | 0 |
back to top
Genotypic Map is necessary in
SIMER
. The first column is SNP
name, the second column is Chromosome
ID, the third column is physical
position, the fourth column is
REF, and the fifth column is
ALT. It will be used to generate
annotation data, genotype
data, and phenotype data.
map.txt
SNP | Chrom | BP | REF | ALT |
---|---|---|---|---|
1_10673082 | 1 | 10673082 | T | C |
1_10723065 | 1 | 10723065 | A | G |
1_11407894 | 1 | 11407894 | A | G |
1_11426075 | 1 | 11426075 | T | C |
1_13996200 | 1 | 13996200 | T | C |
1_14638936 | 1 | 14638936 | T | C |
back to top
SIMER
supports user designed
pedigree to control mating process. User
designed pedigree is useful only in userped
reproduction. The first column is sample id,
the second column is paternal id, and the
third column is maternal id. Please make sure
that paternal id and maternal
id can match to genotype
data.
userped.txt
Index | Sire | Dam |
---|---|---|
41 | 1 | 11 |
42 | 1 | 11 |
43 | 1 | 11 |
44 | 1 | 11 |
45 | 2 | 12 |
46 | 2 | 12 |
back to top
At least users should prepare two datasets: genotypic
map and genotype data.
genotype data,
Numeric format (m
rows and n columns,
m is the number of SNPs,
n is the number of individuals)
genotypic map, SNP map information, the first
column is SNP name, the second column is
Chromosome ID, the third column is
physical position, the fourth column is
REF, and the fifth column is
ALT.
<- read.table("genotype.txt")
pop.geno <- read.table("map.txt" , head = TRUE) pop.map
back to top
Mating process can be designed by user designed
pedigree.
pedigree, pedigree information, the first column is sample id, the second column is paternal id, and the third column is maternal id. Note that the individuals in the pedigree do not need to be sorted by the date of birth, and the missing value can be replaced by NA or 0.
<- read.table("userped.txt", header = TRUE) userped
All simulation processes can be devided into 2 steps: 1) generate simulation parameters; 2) run simulation process.
A quick start for Complete Simulation is shown below:
# Generate all simulation parameters
<- param.simer(out = "simer")
SP
# Run Simer
<- simer(SP) SP
A quick start for Annotation Simulation is shown below:
# Real genotypic map
# pop.map <- read.table("Real_Genotypic_map.txt", header = TRUE)
# Simulated genotypic map
<- generate.map(pop.marker = 1e4)
pop.map
# Generate annotation simulation parameters
<- param.annot(pop.map = pop.map, qtn.num = 10)
SP
# Run annotation simulation
<- annotation(SP) SP
A quick start for Genotype Simulation is shown below:
# Generate genotype simulation parameters
<- param.geno(pop.marker = 1e4, pop.ind = 1e2)
SP
# Run genotype simulation
<- genotype(SP) SP
A quick start for Phenotype Simulation is shown below:
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
Genotype data in
SIMER
will be generated randomly or from
outside genotype matrix. Chromosome crossovers and base mutations depend
on block information and recombination information of
Annotation data.
genotype
, main function of Genotype
Simulation:
Paramater | Default | Options | Description |
pop.geno | NULL | big.matrix or matrix | the genotype data. |
incols | 1 | 1 or 2 | ‘1’: one-column genotype represents an individual; ‘2’: two-column genotype represents an individual. |
pop.marker | 1e4 | num | the number of markers. |
pop.ind | 1e2 | num | the number of individuals in the base population. |
prob | NULL | num vector | the genotype code probability. |
rate.mut | 1e-8 | num | the mutation rate of the genotype data. |
annotation
, main function of Annotation
Simulation:
Paramater | Default | Options | Description |
recom.spot | FALSE | TRUE or FALSE | whether to generate recombination events. |
range.hot | 4:6 | num vector | the recombination times range in the hot spot. |
range.cold | 1:5 | num vector | the recombination times range in the cold spot. |
There are two different ways to generate genotype matrix of base population.
### 01 Use Genotype Data from Outside ###
# Create a genotype matrix
<- matrix(0, nrow = 1e4, ncol = 1e2)
pop.geno
# Generate genotype simulation parameters
<- param.geno(pop.geno = pop.geno)
SP
# Run genotype simulation
<- genotype(SP)
SP
### 02 Create Genotype Data Randomly ###
# Generate genotype simulation parameters
<- param.geno(pop.marker = 1e4, pop.ind = 1e2)
SP
# Run genotype simulation
<- genotype(SP) SP
With annotation data, chromosome crossovers and mutations can be added to genotype matrix.
# Generate annotation simulation parameters
# If recom.spot = TRUE, chromsome crossovers will be added to genotype matrix
<- param.annot(recom.spot = TRUE)
SP # Generate genotype simulation parameters
# Base mutation rate is 1e8
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2, rate.mut = 1e-8)
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP) SP
Note that recombination only exists in meiosis. Therefore, some
reproduction methods such as clone
do not have
recombination process. Users can set recom.spot = FALSE
to
add only mutations to genotype matrix.
# Generate annotation simulation parameters
# If recom.spot = FALSE, chromsome crossovers will not be added to genotype matrix
<- param.annot(recom.spot = FALSE)
SP # Generate genotype simulation parameters
# Base mutation rate is 1e8
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2, rate.mut = 1e-8)
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP) SP
Phenotype data in
SIMER
will be generated according to
different models, including:
(1) Single-Trait Model
(2) Multiple-Trait Model
(3) Repeated Record Model
(4) Genetic Effect Model
(Additive effect,
Dominant effect, and
Genetic-Genetic
interaction effect)
(5) Effect Distribution Model (QTN effect distribution:
****Normal distribution, Geometric
distribution, Gamma*** distribution, and
Beta distribution)
(6) Linear Mixed Model (Fixed
effect, Environmental
Random effect,
Genetic Random
effect, and
Genetic-Environmental
effect)
phenotype
, main function of Phenotype
Simulation:
Paramater | Default | Options | Description |
pop | NULL | data.frame | the population information containing environmental factors and other effects. |
pop.ind | 100 | num | the number of individuals in the base population. |
pop.rep | 1 | num | the repeated times of repeated records. |
pop.rep.bal | TRUE | TRUE or FALSE | whether repeated records are balanced. |
pop.env | NULL | list | a list of environmental factors setting. |
phe.model | list(tr1 = “T1 = A + E”) | list | a list of genetic model of phenotype such as “T1 = A + E”. |
phe.h2A | list(tr1 = 0.3) | list | a list of additive heritability. |
phe.h2D | list(tr1 = 0.1) | list | a list of dominant heritability. |
phe.h2GxG | NULL | list | a list of GxG interaction heritability. |
phe.h2GxE | NULL | list | a list of GxE interaction heritability. |
phe.h2PE | NULL | list | a list of permanent environmental heritability. |
phe.var | NULL | list | a list of phenotype variance. |
phe.corA | NULL | matrix | the additive genetic correlation matrix. |
phe.corD | NULL | matrix | the dominant genetic correlation matrix. |
phe.corGxG | NULL | list | a list of the GxG genetic correlation matrix. |
phe.corPE | NULL | matrix | the permanent environmental correlation matrix. |
phe.corE | NULL | matrix | the residual correlation matrix. |
annotation
, main function of Annotation
Simulation:
Paramater | Default | Options | Description |
pop.map | NULL | data.frame | the map data with annotation information. |
qtn.num | 10 | num, vector or matrix | integer: the QTN number of single trait; vector: the multiple group QTN number of single trait; matrix: the QTN number of multiple traits. |
qtn.model | ‘A’ | character | the genetic model of QTN such as ‘A + D’. |
qtn.dist | list(tr1 = ‘norm’) | list | the QTN distribution containing ‘norm’, ‘geom’, ‘gamma’ or ‘beta’. |
qtn.sd | list(tr1 = 1) | list | the standard deviations for normal distribution. |
qtn.prob | NULL | list | the probability of success for geometric distribution. |
qtn.shape | NULL | list | the shape parameter for gamma distribution. |
qtn.scale | NULL | list | the scale parameter for gamma distribution. |
qtn.shape1 | NULL | list | the shape1 parameter for beta distribution. |
qtn.shape2 | NULL | list | the shape2 parameter for beta distribution. |
qtn.ncp | NULL | list | the ncp parameter for beta distribution. |
qtn.spot | NULL | list | the QTN distribution probability in each block. |
len.block | 5e7 | num | the block length. |
maf | NULL | num | the maf threshold, markers less than this threshold will be exclude. |
In A model,
SIMER
only considers
Additive effect as genetic effect. Users
should prepare Additive
QTN effect in the Annotation
data for generating Additive
Individual effect.
Additive single-trait simulation is displayed
as follows:
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10, qtn.model = "A") # Additive effect
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(tr1 = "T1 = A + E"), # "T1" (Trait 1) consists of Additive effect and Residual effect
# phe.var = list(tr1 = 100),
phe.h2A = list(tr1 = 0.3)
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
In multiple-trait simulation, SIMER
can
build accurate Additive genetic correlation
between multiple traits. Additive
multiple-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(qtn.num = matrix(c(6, 4, 4, 6), 2, 2), qtn.model = "A") # Additive effect
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(
tr1 = "T1 = A + E", # "T1" (Trait 1) consists of Additive effect and Residual effect
tr2 = "T2 = A + E" # "T2" (Trait 2) consists of Additive effect and Residual effect
),# phe.var = list(tr1 = 100, tr2 = 100),
phe.h2A = list(tr1 = 0.3, tr2 = 0.3),
phe.corA = matrix(c(1, 0.5, 0.5, 1), 2, 2) # Additive genetic correlation
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
In single-trait simulation, the trait can be controlled by multiple-group QTNs. An example of single-trait controlled by two-group QTNs is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(qtn.num = c(2, 8), qtn.model = "A") # Group1: 2 QTNs; Group 2: 8 QTNs
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(tr1 = "T1 = A + E"), # "T1" (Trait 1) consists of Additive effect and Residual effect
# phe.var = list(tr1 = 100),
phe.h2A = list(tr1 = 0.3)
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
In AD model,
SIMER
considers
Additive effect and
Dominant effect as genetic effect. Users
should prepare Additive
QTN effect and
Dominant QTN effect
in the Annotation data for generating
Additive Individual
effect and Dominant
Individual effect.
Additive and
Dominant single-trait simulation is displayed
as follows:
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10, qtn.model = "A + D") # Additive effect and Dominant effect
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(tr1 = "T1 = A + D + E"), # "T1" (Trait 1) consists of Additive effect, Dominant effect, and Residual effect
# phe.var = list(tr1 = 100),
phe.h2A = list(tr1 = 0.3),
phe.h2D = list(tr1 = 0.1)
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
In multiple-trait simulation, SIMER
can
build accurate Additive genetic correlation
and accurate Dominant genetic correlation
between multiple traits. Additive and
Dominant multiple-trait simulation is
displayed as follows:
# Generate annotation simulation parameters
<- param.annot(qtn.num = matrix(c(6, 4, 4, 6), 2, 2), qtn.model = "A + D") # Additive effect and Dominant effect
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(
tr1 = "T1 = A + D + E", # "T1" (Trait 1) consists of Additive effect, Dominant effect, and Residual effect
tr2 = "T2 = A + D + E" # "T2" (Trait 2) consists of Additive effect, Dominant effect, and Residual effect
),# phe.var = list(tr1 = 100, tr2 = 100),
phe.h2A = list(tr1 = 0.3, tr2 = 0.3),
phe.h2D = list(tr1 = 0.1, tr2 = 0.1),
phe.corA = matrix(c(1, 0.5, 0.5, 1), 2, 2), # Additive genetic correlation
phe.corD = matrix(c(1, 0.5, 0.5, 1), 2, 2) # Dominant genetic correlation
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
In GxG model,
SIMER
considers
Genetic-Genetic
effect as genetic effect. Users should prepare
Genetic-Genetic
QTN effect in the Annotation
data for generating
Genetic-Genetic
Individual effect. An example of
Additive-Dominant
interaction in single-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10, qtn.model = "A + D + A:D") # Additive effect, Dominant effect, and Additive-Dominant interaction effect
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(tr1 = "T1 = A + D + A:D + E"), # "T1" (Trait 1) consists of Additive effect, Dominant effect, Additive-Dominant interaction effect, and Residual effect
# phe.var = list(tr1 = 100),
phe.h2A = list(tr1 = 0.3),
phe.h2D = list(tr1 = 0.1),
phe.h2GxG = list(tr1 = list("A:D" = 0.1))
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
In multiple-trait simulation, SIMER
can
build accurate Genetic-Genetic interaction
correlation between multiple traits. An example of
Additive-Dominant
interaction in multiple-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(qtn.num = matrix(c(6, 4, 4, 6), 2, 2), qtn.model = "A + D + A:D") # Additive effect, Dominant effect, and Additive-Dominant interaction effect
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(
tr1 = "T1 = A + D + A:D + E", # "T1" (Trait 1) consists of Additive effect, Dominant effect, Additive-Dominant interaction effect, and Residual effect
tr2 = "T2 = A + D + A:D + E" # "T2" (Trait 2) consists of Additive effect, Dominant effect, Additive-Dominant interaction effect, and Residual effect
),# phe.var = list(tr1 = 100, tr2 = 100),
phe.h2A = list(tr1 = 0.3, tr2 = 0.3),
phe.h2D = list(tr1 = 0.1, tr2 = 0.1),
phe.h2GxG = list(tr1 = list("A:D" = 0.1), tr2 = list("A:D" = 0.1)),
phe.corA = matrix(c(1, 0.5, 0.5, 1), 2, 2), # Additive genetic correlation
phe.corD = matrix(c(1, 0.5, 0.5, 1), 2, 2), # Dominant genetic correlation
phe.corGxG = list("A:D" = matrix(c(1, 0.5, 0.5, 1), 2, 2)) # Additive-Dominant interaction genetic correlation
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
In Repeated Record model,
SIMER
adds PE
(Permanent
Environmental) effect to the phenotype. The
number of repeated records can be set by pop.rep
. In the
meantime, pop.rep.bal
can be used to determine whether
repeated records are balanced. Repeated Record
in single-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10, qtn.model = "A") # Additive effect
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
pop.rep = 2, # The number of repeated records is 2
pop.rep.bal = TRUE, # Repeated records are balanced
phe.model = list(tr1 = "T1 = A + E"), # "T1" (Trait 1) consists of Additive effect and Residual effect
# phe.var = list(tr1 = 100),
phe.h2A = list(tr1 = 0.3)
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
In multiple-trait simulation, SIMER
can
build accurate Permanent Environmental
correlation between multiple traits. Repeated
Record in multiple-trait simulation is displayed as
follows:
# Generate annotation simulation parameters
<- param.annot(qtn.num = matrix(c(6, 4, 4, 6), 2, 2), qtn.model = "A") # Additive effect
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
pop.rep = 2, # The number of repeated records is 2
pop.rep.bal = TRUE, # Repeated records are balanced
phe.model = list(
tr1 = "T1 = A + E", # "T1" (Trait 1) consists of Additive effect and Residual effect
tr2 = "T2 = A + E" # "T2" (Trait 2) consists of Additive effect and Residual effect
),# phe.var = list(tr1 = 100, tr2 = 100),
phe.h2A = list(tr1 = 0.3, tr2 = 0.3),
phe.corA = matrix(c(1, 0.5, 0.5, 1), 2, 2), # Additive genetic correlation
phe.corPE = matrix(c(1, 0.5, 0.5, 1), 2, 2) # Permanent Environmental correlation
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
Normal distribution is the most common QTN effect distribution. Phenotype controlled by QTNs subject to Normal distribution in single-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(
SP qtn.num = 10,
qtn.model = "A",
qtn.dist = list(tr1 = "norm"),
qtn.sd = list(tr1 = 1)
)# Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(tr1 = "T1 = A + E"), # "T1" (Trait 1) consists of Additive effect and Residual effect
# phe.var = list(tr1 = 100),
phe.h2A = list(tr1 = 0.3)
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
Phenotype controlled by QTNs subject to Normal distribution in multiple-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(
SP qtn.num = matrix(c(6, 4, 4, 6), 2, 2),
qtn.model = "A",
qtn.dist = list(tr1 = "norm", tr2 = "norm"),
qtn.sd = list(tr1 = 1, tr2 = 1)
)# Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(
tr1 = "T1 = A + E", # "T1" (Trait 1) consists of Additive effect and Residual effect
tr2 = "T2 = A + E" # "T2" (Trait 2) consists of Additive effect and Residual effect
),# phe.var = list(tr1 = 100, tr2 = 100),
phe.h2A = list(tr1 = 0.3, tr2 = 0.3),
phe.corA = matrix(c(1, 0.5, 0.5, 1), 2, 2) # Additive genetic correlation
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
Geometric distribution is the probability of success for the first time obtained only after K trials among the N Bernoulli trials. Geometric distribution can be used as a QTN effect distribution. Phenotype controlled by QTNs subject to Geometric distribution in single-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(
SP qtn.num = 10,
qtn.model = "A",
qtn.dist = list(tr1 = "geom"),
qtn.prob = list(tr1 = 0.5)
)# Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(tr1 = "T1 = A + E"), # "T1" (Trait 1) consists of Additive effect and Residual effect
# phe.var = list(tr1 = 100),
phe.h2A = list(tr1 = 0.3)
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
Phenotype controlled by QTNs subject to Geometric distribution in multiple-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(
SP qtn.num = matrix(c(6, 4, 4, 6), 2, 2),
qtn.model = "A",
qtn.dist = list(tr1 = "geom", tr2 = "geom"),
qtn.prob = list(tr1 = 0.5, tr2 = 0.5)
)# Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(
tr1 = "T1 = A + E", # "T1" (Trait 1) consists of Additive effect and Residual effect
tr2 = "T2 = A + E" # "T2" (Trait 2) consists of Additive effect and Residual effect
),# phe.var = list(tr1 = 100, tr2 = 100),
phe.h2A = list(tr1 = 0.3, tr2 = 0.3),
phe.corA = matrix(c(1, 0.5, 0.5, 1), 2, 2) # Additive genetic correlation
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
Gamma distribution is the sum of N
independent exponential random variables. Note that
Exponential distribution is a special form of
Gamma distribution when
qtn.shape = 1
and qtn.scale = 1
. Phenotype
controlled by QTNs subject to Gamma
distribution in single-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(
SP qtn.num = 10,
qtn.model = "A",
qtn.dist = list(tr1 = "gamma"),
qtn.shape = list(tr1 = 1),
qtn.scale = list(tr1 = 1)
)# Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(tr1 = "T1 = A + E"), # "T1" (Trait 1) consists of Additive effect and Residual effect
# phe.var = list(tr1 = 100),
phe.h2A = list(tr1 = 0.3)
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
Phenotype controlled by QTNs subject to Gamma distribution in multiple-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(
SP qtn.num = matrix(c(6, 4, 4, 6), 2, 2),
qtn.model = "A",
qtn.dist = list(tr1 = "gamma", tr2 = "gamma"),
qtn.shape = list(tr1 = 1, tr2 = 1),
qtn.scale = list(tr1 = 1, tr2 = 1)
)# Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(
tr1 = "T1 = A + E", # "T1" (Trait 1) consists of Additive effect and Residual effect
tr2 = "T2 = A + E" # "T2" (Trait 2) consists of Additive effect and Residual effect
),# phe.var = list(tr1 = 100, tr2 = 100),
phe.h2A = list(tr1 = 0.3, tr2 = 0.3),
phe.corA = matrix(c(1, 0.5, 0.5, 1), 2, 2) # Additive genetic correlation
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
Beta distribution is a density function of conjugate prior distribution as Bernoulli distribution and Binomial distribution. Phenotype controlled by QTNs subject to Beta distribution in single-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(
SP qtn.num = 10,
qtn.model = "A",
qtn.dist = list(tr1 = "beta"),
qtn.shape1 = list(tr1 = 1),
qtn.shape2 = list(tr1 = 1),
qtn.ncp = list(tr1 = 0)
)# Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(tr1 = "T1 = A + E"), # "T1" (Trait 1) consists of Additive effect and Residual effect
# phe.var = list(tr1 = 100),
phe.h2A = list(tr1 = 0.3)
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
Phenotype controlled by QTNs subject to Beta distribution in multiple-trait simulation is displayed as follows:
# Generate annotation simulation parameters
<- param.annot(
SP qtn.num = matrix(c(6, 4, 4, 6), 2, 2),
qtn.model = "A",
qtn.dist = list(tr1 = "beta", tr2 = "beta"),
qtn.shape1 = list(tr1 = 1, tr2 = 1),
qtn.shape2 = list(tr1 = 1, tr2 = 1),
qtn.ncp = list(tr1 = 0, tr2 = 0)
)# Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
phe.model = list(
tr1 = "T1 = A + E", # "T1" (Trait 1) consists of Additive effect and Residual effect
tr2 = "T2 = A + E" # "T2" (Trait 2) consists of Additive effect and Residual effect
),# phe.var = list(tr1 = 100, tr2 = 100),
phe.h2A = list(tr1 = 0.3, tr2 = 0.3),
phe.corA = matrix(c(1, 0.5, 0.5, 1), 2, 2) # Additive genetic correlation
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
SIMER
supports add
Fixed effects and
Environmental Random
effects to phenotype. Users should prepare a list of environmental
factors setting. The effect value of Fixed
effect is determined by eff
and the variance ratio to
phenotype variance (similar to heritability) of
Environmental Random
is set by ratio
. Phenotype with
Fixed effect and
Environmental Random
effect in single-trait simulation is displayed as follows:
# Prepare environmental factor list
<- list(
pop.env F1 = list( # fixed effect 1
level = c("1", "2"),
eff = list(tr1 = c(50, 30))
), F2 = list( # fixed effect 2
level = c("d1", "d2", "d3"),
eff = list(tr1 = c(10, 20, 30))
),R1 = list( # random effect 1
level = c("l1", "l2", "l3"),
ratio = list(tr1 = 0.1)
)
)
# Generate genotype simulation parameters
<- param.annot(qtn.num = 10, qtn.model = "A")
SP # Generate annotation simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
pop.env = pop.env,
phe.model = list(tr1 = "T1 = A + F1 + F2 + R1 + E"), # "T1" (Trait 1) consists of Additive effect, F1, F2, R1, and Residual effect
# phe.var = list(tr1 = 100),
phe.h2A = list(tr1 = 0.3)
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
Phenotype with Fixed effect and Environmental Random effect in multiple-trait simulation is displayed as follows:
# Prepare environmental factor list
<- list(
pop.env F1 = list( # fixed effect 1
level = c("1", "2"),
eff = list(tr1 = c(50, 30), tr2 = c(50, 30))
), F2 = list( # fixed effect 2
level = c("d1", "d2", "d3"),
eff = list(tr1 = c(10, 20, 30), tr2 = c(10, 20, 30))
),R1 = list( # random effect 1
level = c("l1", "l2", "l3"),
ratio = list(tr1 = 0.1, tr2 = 0.1)
)
)
# Generate genotype simulation parameters
<- param.annot(qtn.num = matrix(c(6, 4, 4, 6), 2, 2), qtn.model = "A")
SP # Generate annotation simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
pop.env = pop.env,
phe.model = list(
tr1 = "T1 = A + F1 + F2 + R1 + E", # "T1" (Trait 1) consists of Additive effect, F1, F2, R1, and Residual effect
tr2 = "T2 = A + F1 + F2 + R1 + E" # "T2" (Trait 1) consists of Additive effect, F1, F2, R1, and Residual effect
),# phe.var = list(tr1 = 100, tr2 = 100),
phe.h2A = list(tr1 = 0.3, tr2 = 0.3),
phe.corA = matrix(c(1, 0.5, 0.5, 1), 2, 2) # Additive genetic correlation
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
In GxE model,
SIMER
adds
Genetic-Environmental
interaction effect to phenotype. Users should prepare
Genetic QTN effect
in the Annotation data and environmental
factor by pop.env
for generating
Genetic-Environmental
Individual effect. An example of
Genetic-Environmental
interaction in single-trait simulation is displayed as follows:
# Prepare environmental factor list
<- list(
pop.env F1 = list( # fixed effect 1
level = c("1", "2"),
eff = list(tr1 = c(50, 30))
), F2 = list( # fixed effect 2
level = c("d1", "d2", "d3"),
eff = list(tr1 = c(10, 20, 30))
),R1 = list( # random effect 1
level = c("l1", "l2", "l3"),
ratio = list(tr1 = 0.1)
)
)
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10, qtn.model = "A") # Additive effect
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
pop.env = pop.env,
phe.model = list(
tr1 = "T1 = A + F1 + F2 + R1 + A:F1 + E" # "T1" (Trait 1) consists of Additive effect, F1, F2, R1, Additive-F1 interaction effect, and Residual effect
),# phe.var = list(tr1 = 100),
phe.h2A = list(tr1 = 0.3),
phe.h2GxE = list(tr1 = list("A:F1" = 0.1))
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
An example of Genetic-Environmental interaction in multiple-trait simulation is displayed as follows:
# Prepare environmental factor list
<- list(
pop.env F1 = list( # fixed effect 1
level = c("1", "2"),
eff = list(tr1 = c(50, 30), tr2 = c(50, 30))
), F2 = list( # fixed effect 2
level = c("d1", "d2", "d3"),
eff = list(tr1 = c(10, 20, 30), tr2 = c(10, 20, 30))
),R1 = list( # random effect 1
level = c("l1", "l2", "l3"),
ratio = list(tr1 = 0.1, tr2 = 0.1)
)
)
# Generate annotation simulation parameters
<- param.annot(qtn.num = matrix(c(6, 4, 4, 6), 2, 2), qtn.model = "A") # Additive effect
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
pop.env = pop.env,
phe.model = list(
tr1 = "T1 = A + F1 + F2 + R1 + A:F1 + E", # "T1" (Trait 1) consists of Additive effect, F1, F2, R1, Additive-F1 interaction effect, and Residual effect
tr2 = "T2 = A + F1 + F2 + R1 + A:F1 + E" # "T2" (Trait 2) consists of Additive effect, F1, F2, R1, Additive-F1 interaction effect, and Residual effect
),# phe.var = list(tr1 = 100, tr2 = 100),
phe.h2A = list(tr1 = 0.3, tr2 = 0.3),
phe.corA = matrix(c(1, 0.5, 0.5, 1), 2, 2) # Additive genetic correlation
)
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP) SP
SIMER
imitates the reproductive process
of organisms to generate Multiple-Generation
population. The genotype data and
phenotype data of the population are screened
by single-trait selection or
multiple-trait selection, and then amplified
by species-specific reproduction.
selects
, main function of
Selection:
Paramater | Default | Options | Description |
pop.sel | NULL | list | the selected males and females. |
ps | c(0.8, 0.8) | num vector | if ps <= 1, fraction selected in selection of males and females; if ps > 1, ps is number of selected males and females. |
decr | TRUE | TRUE or FALSE | whether the sort order is decreasing. |
sel.crit | ‘pheno’ | character | the selection criteria, it can be ‘TBV’, ‘TGV’, and ‘pheno’. |
sel.single | ‘comb’ | character | the single-trait selection method, it can be ‘ind’, ‘fam’, ‘infam’, and ‘comb’. |
sel.multi | ‘index’ | character | the multiple-trait selection method, it can be ‘index’, ‘indcul’, and ‘tmd’. |
index.wt | c(0.5, 0.5) | num vector | the weight of each trait for multiple-trait selection. |
index.tdm | 1 | num | the index of tandem selection for multiple-trait selection. |
goal.perc | 0.1 | num | the percentage of goal more than the mean of scores of individuals. |
pass.perc | 0.9 | num | the percentage of expected excellent individuals. |
reproduces
, main function of
Reproduction:
Paramater | Default | Options | Description |
pop.gen | 2 | num | the generations of simulated population. |
reprod.way | ‘randmate’ | character | reproduction method, it consists of ‘clone’, ‘dh’, ‘selfpol’, ‘randmate’, ‘randexself’, ‘2waycro’, ‘3waycro’, ‘4waycro’, ‘backcro’, and ‘userped’. |
sex.rate | 0.5 | num | the male rate in the population. |
prog | 2 | num | the progeny number of an individual. |
Individual selection is a selecting method according to the phenotype of individual traits, also known as mixed selection or collective selection. This selection method is simple and easy to used for traits with high heritability.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "ind")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP) SP
Family selection is a selecting method by family based on the average of the family. This selection method is used for traits with low heritability.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "fam")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP) SP
Within-family selection is a selecting method according to the deviation of individual phenotype and family mean value in each family. This selection method is used for traits with low heritability and small family.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "infam")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP) SP
Combined selection is a selecting method according to weighed combination of the deviation of individual phenotype and family mean value.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "comb")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP) SP
Tandem selection is a method for
sequentially selecting a plurality of target traits one by
one. The index of selected trait is index.tdm
and this parameter should not controlled by
Users.
# Generate genotype simulation parameters
<- param.annot(qtn.num = matrix(c(6, 4, 4, 6), 2, 2), qtn.model = "A")
SP # Generate annotation simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
# phe.var = list(tr1 = 100, tr2 = 100),
phe.model = list(
tr1 = "T1 = A + E",
tr2 = "T2 = A + E"
)
)# Generate selection parameters
<- param.sel(SP = SP, sel.multi = "tdm")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP) SP
After setting a minimum selection criterion for each target trait. Independent culling selection will eliminate this individual when the candidate’s performance on any trait is lower than the corresponding criteria.
# Generate genotype simulation parameters
<- param.annot(qtn.num = matrix(c(6, 4, 4, 6), 2, 2), qtn.model = "A")
SP # Generate annotation simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
# phe.var = list(tr1 = 100, tr2 = 100),
phe.model = list(
tr1 = "T1 = A + E",
tr2 = "T2 = A + E"
)
)# Generate selection parameters
<- param.sel(SP = SP, sel.multi = "indcul")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP) SP
Index selection is a comprehensive
selection that will consider several traits based on their respective
heritabilities, phenotypic
variances, economic weights,
corresponding genetic correlations, and
phenotypes. Then calculate the
index value of each trait, and eliminate or
select it according to its level. Users can set the weight of each trait
by index.wt
.
# Generate genotype simulation parameters
<- param.annot(qtn.num = matrix(c(6, 4, 4, 6), 2, 2), qtn.model = "A")
SP # Generate annotation simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(
SP SP = SP,
pop.ind = 100,
# phe.var = list(tr1 = 100, tr2 = 100),
phe.model = list(
tr1 = "T1 = A + E",
tr2 = "T2 = A + E"
)
)# Generate selection parameters
<- param.sel(SP = SP, sel.multi = "index")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP) SP
Clone is an sexual reproduction method
which does not involve germ cells, and does not require a process of
fertilization, directly forming a new individual’s reproductive mode
from a part of the mother. Sex of offspring
will be 0 in clone
.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "comb")
SP # Generate reproduction parameters
<- param.reprod(SP = SP, reprod.way = "clone")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP)
SP # Run reproduction
<- reproduces(SP) SP
Double haploid is a reproduction method for
breeding workers to obtain haploid plants. It induced to double the
number of chromosomes and restore the number of chromosomes in normal
plants. Sex of offspring will be
0 in dh
.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "comb")
SP # Generate reproduction parameters
<- param.reprod(SP = SP, reprod.way = "dh")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP)
SP # Run reproduction
<- reproduces(SP) SP
Self pollination refers to the combination
of male and female gametes from the same individual or between
individuals from the same clonal breeding line.
Sex of offspring will be
0 in selfpol
.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "comb")
SP # Generate reproduction parameters
<- param.reprod(SP = SP, reprod.way = "selfpol")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP)
SP # Run reproduction
<- reproduces(SP) SP
In random mating, any female or male
individual have the same probability to mate with any opposite sex in a
sexually reproducing organism. Sex of
offspring in random mating is controlled by sex.ratio
in
randmate
.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "comb")
SP # Generate reproduction parameters
<- param.reprod(SP = SP, reprod.way = "randmate")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP)
SP # Run reproduction
<- reproduces(SP) SP
In random mating excluding self
pollination, an individual cannot mate to itself.
Sex of offspring in random mating is
controlled by sex.ratio
in randexself
.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "comb")
SP # Generate reproduction parameters
<- param.reprod(SP = SP, reprod.way = "randexself")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run selection
<- selects(SP)
SP # Run reproduction
<- reproduces(SP) SP
Two-way cross method needs to use sex to distinguish two different breeds, in which the first breed is sire and the second breed is dam.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "comb")
SP # Generate reproduction parameters
<- param.reprod(SP = SP, reprod.way = "2waycro")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Two different breeds are cut by sex
$pheno$pop$gen1$sex <- rep(c(1, 2), c(50, 50))
SP# Run selection
<- selects(SP)
SP # Run reproduction
<- reproduces(SP) SP
Three-way cross method needs to use sex to distinguish three different breeds, in which the first breed is sire and the second breed is dam in the first two-way cross, the third breed is termimal sire.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "comb")
SP # Generate reproduction parameters
<- param.reprod(SP = SP, reprod.way = "3waycro")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Three different breeds are cut by sex
$pheno$pop$gen1$sex <- rep(c(1, 2, 1), c(30, 30, 40))
SP# Run selection
<- selects(SP)
SP # Run reproduction
<- reproduces(SP) SP
Four-way cross method needs to use sex to distinguish four different breeds, in which the first breed is sire and the second breed is dam in the first two-way cross, the third breed is sire and the fourth breed is dam in the second two-way cross.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "comb")
SP # Generate reproduction parameters
<- param.reprod(SP = SP, reprod.way = "4waycro")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Three different breeds are cut by sex
$pheno$pop$gen1$sex <- rep(c(1, 2, 1, 2), c(25, 25, 25, 25))
SP# Run selection
<- selects(SP)
SP # Run reproduction
<- reproduces(SP) SP
Back cross method needs to use sex to distinguish two different breeds, in which the first breed is always sire in each generation and the second breed is dam in the first two-way cross.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate selection parameters
<- param.sel(SP = SP, sel.single = "comb")
SP # Generate reproduction parameters
<- param.reprod(SP = SP, reprod.way = "backcro")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Two different breeds are cut by sex
$pheno$pop$gen1$sex <- rep(c(1, 2), c(50, 50))
SP# Run selection
<- selects(SP)
SP # Run reproduction
<- reproduces(SP) SP
User designed pedigree mating needs a specific user designed pedigree to control mating process. The first column is sample id, the second column is paternal id, and the third column is maternal id. Please make sure that paternal id and maternal id can match to genotype data.
# Generate annotation simulation parameters
<- param.annot(qtn.num = 10)
SP # Generate genotype simulation parameters
<- param.geno(SP = SP, pop.marker = 1e4, pop.ind = 1e2)
SP # Generate phenotype simulation parameters
<- param.pheno(SP = SP, pop.ind = 100)
SP # Generate reproduction parameters
<- param.reprod(SP = SP, reprod.way = "userped")
SP
# Run annotation simulation
<- annotation(SP)
SP # Run genotype simulation
<- genotype(SP)
SP # Run phenotype simulation
<- phenotype(SP)
SP # Run reproduction
<- reproduces(SP) SP
The above methods are to generate population step by
step, which are easy to
understand. Actually, SIMER
can directly generate a population by a MORE
CONVENIENT way.
# Generate all simulation parameters
<- param.simer(qtn.num = 10, pop.marker = 1e4, pop.ind = 1e2, sel.single = "comb", reprod.way = "randmate")
SP
# Run Simer
<- simer(SP) SP
After generating a population, further work can be done. Breeders
wish to evaluate their Breeding Program
Design. To save a lot of money and time,
SIMER
can assist breeders to evaluate
their Breeding Program Design.
read.selgeno()
, function to make comparison on breeding
plans:
pop, total population information
selPath, the path of breeding plans
outpath, path of output files
Breeding plans should be stored on different files respectively. Filenames must begin with breeding_plan, such like breeding_plan01.txt. For now, SIMER supports different breeding plans on generation, family_index, within_family_index, and sex.
> # get path of breeding plan
> selPath <- system.file("extdata", "01breeding_plan", package = "simer")
> # show filename format
> dir(selPath)
1] "breeding_plan01.txt" "breeding_plan02.txt" "breeding_plan03.txt" "ReadMe.txt"
[
> # show file contents format
> fn1 <- file.path(selPath, "breeding_plan01.txt")
> bp1 <- read.delim(fn1, header = TRUE)
> # ###NOTE###
> # 1. any generation can be choosed within num.gen.
> # 2. any sex(1 represents sir 2 represents dam) can be choosed.
> # 3. any family index can be choosed in all family. 1:5 means
> # that the first 5 family will be chosen.
> # 4. any infamily index can be choosed in every family.
> # 5. numbers should be separated by comma(",").
> # 6. ":" can choose serial number, "1:3,7" represents "1 2 3 7".
> # 7. "all" will choose all of options in a category.
> bp1
generation family_index within_family_index sex1 1 1:3,7 all 1
2 2 1:3,7 all 1
3 3 1:3,7 all 1
After a total reproduction process, breeding plans comparison can be done.
# random-mating
<-
simer.list simer(num.gen = 10,
verbose = verbose,
outpath = outpath,
input.map = input.map,
rawgeno1 = rawgeno, # use your own genotype matrix
# num.ind = 100,
mtd.reprod = "randmate",
num.prog = 2,
ratio = 0.5)
<- simer.list$pop
pop <- read.selgeno(pop = pop, selPath = selPath, outpath = outpath)
out.pop
# make comparison on breeding plans
<- out.pop$plan1
plan1 <- out.pop$plan2
plan2 <- out.pop$plan3
plan3 summary(plan1)
summary(plan2)
summary(plan3)
In this part, calculation of population size and different ourput methods will be introduced.
simer()
, main function:
num.gen, number of generations in simulation
replication, replication index of simulation
verbose, whether to print detail
out, prefix of output file name
outpath, path of output files
out.format, format of output, “numeric” or
“plink”
seed.sim, random seed of a simulation process
seed.map, random seed of map file
out.geno.gen, indice of generations of output
genotype
out.pheno.gen, indice of generations of output
phenotype
The following is the method of obtaining population size of every generation. Every elements in cound.ind are population size in this generation respectively.
# parameters that controls population size of every generations
<- 40
nind <- 10
num.gen <- 0.5
ratio <- TRUE # whether selection exsits
sel.on <- 0.8
ps <- ifelse(sel.on, ps, 1)
ps <- 2
num.prog <- "randmate"
mtd.reprod
# populations of the first generation
<- getpop(nind = nind, from = 1, ratio = ratio)
basepop <- getpop(nind = nind, from = 41, ratio = ratio)
pop2 <- getpop(nind = nind, from = 81, ratio = ratio)
pop3 <- getpop(nind = nind, from = 121, ratio = ratio)
pop4
# calculate number of individuals in every generation
<- rep(nind, num.gen)
count.ind if (mtd.reprod == "clone" || mtd.reprod == "dh" || mtd.reprod == "selfpol") {
if (num.gen > 1) {
for(i in 2:num.gen) {
<- round(count.ind[i-1] * (1-ratio) * ps) * num.prog
count.ind[i]
}
}
else if (mtd.reprod == "randmate" || mtd.reprod == "randexself") {
} if (num.gen > 1) {
for(i in 2:num.gen) {
<- round(count.ind[i-1] * (1-ratio) * ps) * num.prog
count.ind[i]
}
}
else if (mtd.reprod == "singcro") {
} <- round(nrow(pop2) * ps) * num.prog
sing.ind <- c(nrow(basepop), nrow(pop2), sing.ind)
count.ind
else if (mtd.reprod == "tricro") {
} <- round(nrow(pop2) * ps) * prog.tri
dam21.ind <- round(dam21.ind * (1-ratio) * ps) * num.prog
tri.ind <- c(nrow(basepop), nrow(pop2), nrow(pop3), dam21.ind, tri.ind)
count.ind
else if (mtd.reprod == "doubcro") {
} <- round(nrow(pop2) * ps) * prog.doub
sir11.ind <- round(nrow(pop4) * ps) * prog.doub
dam22.ind <- round(dam22.ind * (1-ratio) * ps) * num.prog
doub.ind <- c(nrow(basepop), nrow(pop2), nrow(pop3), nrow(pop4), sir11.ind, dam22.ind, doub.ind)
count.ind
else if (mtd.reprod == "backcro") {
} 1] <- nrow(basepop) + nrow(pop2)
count.ind[if (num.gen > 1) {
2] <- round(nrow(pop2) * ps) * num.prog
count.ind[for(i in 3:num.gen) {
<- round(count.ind[i-1] * (1-ratio) * ps) * num.prog
count.ind[i]
}
} }
Simulation of multiple populations can be realized by “for” in R and replication. Random seed of simulation is random in every replication and random seed of map is fixed. In every replication, you can set your own random seed of simulation and random seed of map by seed.geno| and seed.map respectively.
# random-mating
<- 2
rep for (i in 1:rep) {
<-
simer.list simer(num.gen = 3,
replication = i, # set index of replication
verbose = verbose,
outpath = outpath,
input.map = input.map,
# rawgeno1 = rawgeno, # use your own genotype matrix
num.ind = 100,
mtd.reprod = "randmate",
num.prog = 2,
ratio = 0.5)
}
SIMER won’t output files by default. A series of
files with prefix out = simer
output when specify a exsited
out path by outpath. File output format can be
“numeric” or “plink” by out.format.
# set a output path
= getwd()
outpath
# "numeric" format
<-
simer.list simer(num.gen = 3,
verbose = verbose,
outpath = outpath,
out.format = "numeric",
input.map = input.map,
# rawgeno1 = rawgeno, # use your own genotype matrix
num.ind = 100,
mtd.reprod = "randmate",
num.prog = 2,
ratio = 0.5)
# "plink" format
<-
simer.list simer(num.gen = 3,
verbose = verbose,
outpath = outpath,
out.format = "plink",
input.map = input.map,
# rawgeno1 = rawgeno, # use your own genotype matrix
num.ind = 100,
mtd.reprod = "randmate",
num.prog = 2,
ratio = 0.5)
Output of genotype and phenotype can be generation-selective by out.geno.gen and out.pheno.gen. For example, out.geno.gen = 3:5 and out.pheno.gen = 1:5 represent outputting genotype from generation3 to generation5 and outputting phenotype from generation1 to generation5.
# set output generations of genotype and phenotype
<- 3:5
out.geno.gen <- 1:5
out.pheno.gen <-
simer.list simer(num.gen = 5,
verbose = verbose,
outpath = outpath,
out.geno.gen = out.geno.gen,
out.pheno.gen = out.pheno.gen,
input.map = input.map,
# rawgeno1 = rawgeno, # use your own genotype matrix
num.ind = 100,
mtd.reprod = "randmate",
num.prog = 2,
ratio = 0.5)
SIMER outputs data including population information,
marker effects, trait information, genotype, genotypic id, genotypic
map, and selection intensity. Note that several result files with prefix
out = "simer"
will be generated. Firstly, a result path
will be generated. The number at the beginning represents the total
individuals number and the ending character is the format of output
(“num_Simer_Data_format”). Secondly, you will see a path named
“replication1”. It is the first replication of simulation and you can
get numerous different replications by setting parameter
replication. In “replication1”, results are
following:
simer.geno.id
contains indice of genotyped
individuals
simer.geno.bin
and simer.geno.desc
contain
genotype matrix of all individuals
simer.map
contains input map with block information and
recombination information
simer.ped
contains pedigree of individuals
simer.phe
contains phenotype of individuals
Population information contains generation, individual indice, family indice, within-family indice, sire indice, dam indice, sex, and phenotpye.
> pop <- simer.list$pop
> head(pop)
gen index fam infam sir dam sex pheno1 1 1 1 1 0 0 1 -1.043716
2 1 2 2 2 0 0 1 0.551076
3 1 3 3 3 0 0 1 -5.727581
4 1 4 4 4 0 0 1 5.825591
5 1 5 5 5 0 0 1 4.430751
6 1 6 6 6 0 0 1 15.059641
Marker effects is a list with selected markers and effects of markers.
> effs <- simer.list$effs
> str(effs)
2
List of $ mrk1: int [1:18] 4006 4950 10416 18070 21899 23817 26038 28886 29437 30012 ...
$ eff1:List of 1
$ eff.a: num [1:18] -0.0144 -5.7868 -1.2436 1.2436 -3.4893 ... ..
Trait information is displayed by generation in single-breed reproduction or by breed in multiple-breeds reproduction. In every generation (or breed), it contains trait information (variance components and heritability), effect information (phenotype decomposition), phenotype information (TBV, TGV, pEBVs, gEBVs, ssEBVs, phenotype), and others.
> trait <- simer.list$trait
> str(trait)
5
List of $ gen1:List of 4
$ info.tr :List of 3
..$ Vg: num 23.3
.. ..$ Ve: num 63
.. ..$ h2: num 0.27
.. ..$ info.eff :'data.frame': 40 obs. of 2 variables:
..$ ind.a : num [1:40] -6.46 6.46 6.46 6.46 0 ...
.. ..$ ind.env: num [1:40] 5.419 -5.912 -12.191 -0.638 4.431 ...
.. ..$ info.pheno:'data.frame': 40 obs. of 3 variables:
..$ TBV : num [1:40] -6.46 6.46 6.46 6.46 0 ...
.. ..$ TGV : num [1:40] -6.46 6.46 6.46 6.46 0 ...
.. ..$ pheno: num [1:40] -1.044 0.551 -5.728 5.826 4.431 ...
.. ..$ qtn.a : num [1:18, 1:40] 1 1 1 1 1 1 1 1 1 1 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. .. ..$ : chr [1:40] "V1" "V3" "V5" "V7" ...
.. .. ..$ gen2:List of 4
$ info.tr :List of 3
..$ Vg: num 41.7
.. ..$ Ve: num 113
.. ..$ h2: num 0.27
.. ..$ info.eff :'data.frame': 32 obs. of 2 variables:
..$ ind.a : num [1:32] 5.17 1.29 6.46 6.46 6.46 ...
.. ..$ ind.env: num [1:32] 10.8 -8.61 -2.41 1.8 -5.04 ...
.. ..$ info.pheno:'data.frame': 32 obs. of 3 variables:
..$ TBV : num [1:32] 5.17 1.29 6.46 6.46 6.46 ...
.. ..$ TGV : num [1:32] 5.17 1.29 6.46 6.46 6.46 ...
.. ..$ pheno: num [1:32] 15.97 -7.31 4.06 8.26 1.43 ...
.. ..$ qtn.a : num [1:18, 1:32] 0 -1 0 0 0 0 0 0 0 0 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. .. ..$ : chr [1:32] "V17" "V18" "t1" "t2" ...
.. .. ..$ gen3:List of 4
$ info.tr :List of 3
..$ Vg: num 17.5
.. ..$ Ve: num 36.5
.. ..$ h2: num 0.324
.. ..$ info.eff :'data.frame': 26 obs. of 2 variables:
..$ ind.a : num [1:26] -1.525 -1.525 0 0.609 6.463 ...
.. ..$ ind.env: num [1:26] 3.29 2.31 -6.94 5.97 15.28 ...
.. ..$ info.pheno:'data.frame': 26 obs. of 3 variables:
..$ TBV : num [1:26] -1.525 -1.525 0 0.609 6.463 ...
.. ..$ TGV : num [1:26] -1.525 -1.525 0 0.609 6.463 ...
.. ..$ pheno: num [1:26] 1.77 0.781 -6.936 6.576 21.743 ...
.. ..$ qtn.a : num [1:18, 1:26] 0 0 0 0 0 0 0 0 0 0 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. .. ..$ : chr [1:26] "t2" "t2" "t2" "t1" ...
.. .. ..$ gen4:List of 4
$ info.tr :List of 3
..$ Vg: num 15.8
.. ..$ Ve: num 32
.. ..$ h2: num 0.33
.. ..$ info.eff :'data.frame': 20 obs. of 2 variables:
..$ ind.a : num [1:20] -2.808 -0.242 9.085 -2.013 4.052 ...
.. ..$ ind.env: num [1:20] 3.96 2.64 -2.9 -2.65 2.2 ...
.. ..$ info.pheno:'data.frame': 20 obs. of 3 variables:
..$ TBV : num [1:20] -2.808 -0.242 9.085 -2.013 4.052 ...
.. ..$ TGV : num [1:20] -2.808 -0.242 9.085 -2.013 4.052 ...
.. ..$ pheno: num [1:20] 1.15 2.4 6.19 -4.66 6.26 ...
.. ..$ qtn.a : num [1:18, 1:20] 1 1 1 0 0 1 1 1 0 1 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. .. ..$ : chr [1:20] "V18" "t2" "t1" "t1" ...
.. .. ..$ gen5:List of 4
$ info.tr :List of 3
..$ Vg: num 12.1
.. ..$ Ve: num 28.4
.. ..$ h2: num 0.299
.. ..$ info.eff :'data.frame': 16 obs. of 2 variables:
..$ ind.a : num [1:16] 3.79 4.377 0.551 0.551 7.68 ...
.. ..$ ind.env: num [1:16] -1.22 -5.5 -7.99 8.79 3.1 ...
.. ..$ info.pheno:'data.frame': 16 obs. of 3 variables:
..$ TBV : num [1:16] 3.79 4.377 0.551 0.551 7.68 ...
.. ..$ TGV : num [1:16] 3.79 4.377 0.551 0.551 7.68 ...
.. ..$ pheno: num [1:16] 2.57 -1.12 -7.44 9.34 10.78 ...
.. ..$ qtn.a : num [1:18, 1:16] -1 -1 -1 0 0 0 -1 -1 -1 -1 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. .. ..$ : chr [1:16] "t2" "t1" "t2" "t2" ... .. .. ..
In genotype matrix, each row represents a marker and each column represents a individual.
> geno <- simer.list$geno
> geno[1:6, 1:6]
1] [,2] [,3] [,4] [,5] [,6]
[,1,] 1 1 1 1 0 0
[2,] 1 1 1 1 0 0
[3,] 1 1 1 1 0 0
[4,] 1 1 1 1 0 0
[5,] 1 1 1 1 0 0
[6,] 1 1 1 1 0 0 [
Genotypic id is the indice of genotyped individuals.
> genoid <- simer.list$genoid
> head(genoid)
1] 73 74 75 76 77 78 [
In SNP map information, the first column is SNP name, the second column is Chromosome ID, the third column is phsical position, the fourth column is REF, the fifth column is ALT, the sixth column is block ID, and the seventh is recombination information (“1” represents making recombination and “0” represents not).
> map <- simer.list$map
> head(map)
SNP Chrom BP REF ALT block recom1,] "1_10673082" "1" " 10673082" "T" "C" "1" "1"
[2,] "1_10723065" "1" " 10723065" "A" "G" "1" "1"
[3,] "1_11407894" "1" " 11407894" "A" "G" "1" "1"
[4,] "1_11426075" "1" " 11426075" "T" "C" "1" "1"
[5,] "1_13996200" "1" " 13996200" "T" "C" "1" "1"
[6,] "1_14638936" "1" " 14638936" "T" "C" "1" "1" [
Selection intensity is calculated according to ratio of selected individuals.
> si <- simer.list$si
> si
1] 0.3499524 [
For SIMER:
Hope it will be coming soon!
For ADI model:
Kao, Chenhung, et al. "Modeling Epistasis of Quantitative Trait Loci Using Cockerham's Model." Genetics 160.3 (2002): 1243-1261.
For build.cov:
B. D. Ripley "Stochastic Simulation." Wiley-Interscience (1987): Page 98.
:sos: Question1: Failing to install “devtools”:
ERROR: configuration failed for package ‘git2r’
removing ‘/Users/acer/R/3.4/library/git2r’
ERROR: dependency ‘git2r’ is not available for package ‘devtools’
removing ‘/Users/acer/R/3.4/library/devtools’
:yum: Answer: Please try following codes in terminal:
apt-get install libssl-dev/unstable
:sos: Question2: When installing packages from Github with “devtools”, an error occurred:
Error in curl::curl_fetch_disk(url, x$path, handle = handle): Problem with the SSL CA cert (path? access rights?)
:yum: Answer: Please try following codes and then try agian.
library(httr)
set_config(config(ssl_verifypeer = 0L))
Questions, suggestions, and bug reports are welcome and appreciated. :arrow_right: