Provides functions to simulate single or multiple, independent or correlated phenotype(s) with additive, dominance effects and their interactions. Also includes functions to generate phenotype(s) with specific heritability. Flexible and user-friendly options for simulation.
Bug report git repository for SimPhe
SimPhe is an R package dedicated to simulate multiple phenotyes based on genotyping data. The main feature of this package is the possibility to take genetic epistatic effects, not only additive-additive interaction but also additive-dominance and dominance-dominance, heritability and multiple phenotypes into account. Moreover, we provide a variety of convenient functions, including suggestion for set the variation of random variable according to user specified genetic effects and flexible support for input formats, as well as output formats. It also supports the input of plink formats.
install.packages(SimPhe)# load packagelibrary(SimPhe)First, we show how easy to get phenotype(s) by using sim.phe. Before running sim.phe, we need specify the parameter file and genotype file for simulation. After installing SimPhe, these two toy files exist in package folder. To get the path, type
# get file path of simulation parameters# (two shared SNP pairs and one independent SNP pair for each phenotype)fpar.path <- system.file("extdata", "simupars.txt", package="SimPhe") # get file path of genotype file: rows are individuals and columns are SNPsfgeno.path <- system.file("extdata", "10SNP.txt", package="SimPhe")Then simulate the phenotypes as designed in the parameter file after loading package:
# simulate phenotype(s)phe <- sim.phe(sim.pars = fpar.path, fgeno = fgeno.path, ftype = "snp.head", seed = 123, fwrite = FALSE)In the parameter file, we conduct two phenotypes contrbituted by two common SNP pairs with epistatic effects and one independent SNP pair with epistatic effects, the example of simulation parameters, same as the example file, is also included in SimPhe. Users could touch it by typing genepars.
geneparsPhenotype 1 has been set with certain heritability but phenotype 2 has not. We will show whether the heritability of the simulated phenotype 1 is same as the set. SimPhe includes the coefficients and the allele frequencies for simulating phenotype 1: gene.coef and allele.freq, which is extracted from the simulation parameters.
# regression coefficients for simulating phenotype 1gene.coefficients# allele frequencies of SNPs for simulating phenotype 1allele.freq # calculate genetic variancegenevar <- calc.gene.var(gene.coefficients, allele.freq) # the variance of simulated phenotype 1phe1var <- var(phe[, "p1"]) # heritability of simulated phenotype 1simuht <- genevar / phe1varsimuhtThe result is not the exactly 0.6 due to the (pseudo)random numbers generated in R. To get phenotype 2 with a specific heritability, for example, 0.45, we could proceed as:
# get regression coefficients settings for phenotype 2genecoef <- get.gene.coef(main.pars = specify.pars(genetic.pars = genepars, effect.type = "main", phe.index = 2), epi.pars = specify.pars(genetic.pars = genepars, effect.type = "epistasis", phe.index = 2))# get genotype datagenotype <- read.geno(fname = fgeno.path, ftype = "snp.head") # get allele frequencies of the SNPs set for phenotype 2freq2 <- get.freq(geno = genotype, epi.pars = specify.pars(genetic.pars = genepars, effect.type = "epistasis", phe.index = 2)) # get noise variationexp.noise.var <- get.noise.var(gene.coef = genecoef, freq = freq2, heritability = 0.45)Then when simulating a phenotype, just give this value as argument noise.var to function sim.phe. It will generate a phenotype which has a heritability close to 0.45.
As we mentioned earlier in this article, buiding the correlation by setting the shared interactive SNP pairs cannot be controlled. We can take a look at the correlation between the simulated phenotype 1 and phenotype 2:
# correlation testcor.test(phe[, "p1"], phe[, "p2"])According to the result of correlation test, the two simulated phenotypes are significantly correlated but the correlation coefficient is small and the amount cannot easily be predicted. To get a certain amount or higher correlation, we can conduct correlation by applying the correlation matrix to two independent variables. For two phenotypes, if we set different SNP pairs for each, we assume these two phenotypes are independent. Here we give another parameter file to but same genotype file to the arguments in sim.phe, fgenetic.pars and fgeno:
# get file path of simulation parameters# (different SNP pairs for each phenotype)fpar.path <- system.file("extdata", "sep_simupars.txt", package="SimPhe") # simulate phenotype(s)indphe <- sim.phe(sim.pars = fpar.path, fgeno = fgeno.path, ftype = "snp.head", seed = 123, fwrite = FALSE)We can test the correlation between initial phenotypes with seperated SNP pairs settings:
# test original correlationcor.test(indphe[, "p1"], indphe[, "p2"])Apparently, these two phenotypes are not related. We could continue our work on coverting them to be correlated. First, conduct a correlation matrix:
# correlation matrixcorm <- matrix(c(1, 0.6, 0.6, 1), ncol = 2)cormBefore applying correlation matrix to simulated phenotypes, we would like to know what the data looks like:
apply(indphe, 2, mean)apply(indphe, 2, sd)Then we can build correlation between the two initial phenotypes:
# build correlationcorphe <- build.cor.phe(indphe, corMtr = corm)Now let's test the correlation between the two new phenotypes and see if there is anything difference:
# check mean and standard deviation of new data setapply(corphe, 2, mean)apply(corphe, 2, sd)# test correlationcor.test(corphe[, "p1"], corphe[, "p2"])Obviously, there is no significant difference on means and standard deviations between the initial phenotypes and the new phenotypes.