diff --git a/DESCRIPTION b/DESCRIPTION index 0c3f13be..090b3ac7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: AlphaSimR Type: Package Title: Breeding Program Simulations -Version: 1.2 -Date: 2022-5-15 +Version: 1.2.1 +Date: 2022-7-5 Authors@R: c(person("Chris", "Gaynor", email = "gaynor.robert@hotmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-0558-6656")), person("Gregor", "Gorjanc", role = "ctb", @@ -32,7 +32,7 @@ Encoding: UTF-8 Depends: R (>= 4.0.0), methods, R6 Imports: Rcpp (>= 0.12.7) LinkingTo: Rcpp, RcppArmadillo (>= 0.7.500.0.0), BH -RoxygenNote: 7.1.2 +RoxygenNote: 7.2.0 Suggests: knitr, rmarkdown, testthat VignetteBuilder: knitr NeedsCompilation: true diff --git a/NEWS b/NEWS index d3d0a001..c2755b95 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,9 @@ +Changes in version 1.2.1 + + Bug fixes + -fixed bugs relating to importData functions + -fixed writePlink errors and no longer requires equal length chromosomes + Changes in version 1.2.0 New features diff --git a/R/Class-SimParam.R b/R/Class-SimParam.R index e9c6914a..d9a50943 100644 --- a/R/Class-SimParam.R +++ b/R/Class-SimParam.R @@ -51,7 +51,7 @@ SimParam = R6Class( #' #Set simulation parameters #' SP = SimParam$new(founderPop) initialize = function(founderPop){ - stopifnot(class(founderPop)=="MapPop") + stopifnot(is(founderPop, "MapPop")) # Public items self$nThreads = getNumThreads() @@ -1212,16 +1212,17 @@ SimParam = R6Class( if(!is.null(varE)){ varE = as.numeric(varE) stopifnot(length(varE)==nTraits) + }else{ + varE = rep(NA_real_, nTraits) } # Prepare trait names - if(is.null(names)){ + if(is.null(name)){ name = paste0("Trait",1:nTraits+self$nTraits) }else{ - stopifnot(length(names)==nTraits) + stopifnot(length(name)==nTraits) } - # Extract genetic map and check if marker names are on the map genMapMarkerNames = unlist(lapply(private$.femaleMap, names)) stopifnot(all(markerNames%in%genMapMarkerNames)) @@ -1242,7 +1243,7 @@ SimParam = R6Class( lociLoc[[i]] = integer() # Find matches if they exist - take = match(names(genMap[[i]]), markerNames) + take = match(names(private$.femaleMap[[i]]), markerNames) lociPerChr[i] = length(na.omit(take)) if(lociPerChr[i]>0L){ diff --git a/R/RcppExports.R b/R/RcppExports.R index 2c36e63d..2f5606c4 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -349,10 +349,6 @@ packHaplo <- function(haplo, ploidy, inbred) { .Call(`_AlphaSimR_packHaplo`, haplo, ploidy, inbred) } -writePlinkPed <- function(fam, haplo, nInd, ploidy, nLoc, file) { - invisible(.Call(`_AlphaSimR_writePlinkPed`, fam, haplo, nInd, ploidy, nLoc, file)) -} - MaCS <- function(args, maxSites, inbred, ploidy, nThreads, seed) { .Call(`_AlphaSimR_MaCS`, args, maxSites, inbred, ploidy, nThreads, seed) } diff --git a/R/misc.R b/R/misc.R index a0e504d5..be913c44 100644 --- a/R/misc.R +++ b/R/misc.R @@ -506,110 +506,6 @@ usefulness = function(pop,trait=1,use="gv",p=0.1, return(mean(response)) } -#' @title Writes a Pop-class as PLINK files -#' -#' @description -#' Writes a Pop-class as PLINK PED and MAP files -#' -#' @param pop an object of \code{\link{Pop-class}} -#' @param baseName a character. Basename of PED and MAP files. -#' @param trait an integer. Which phenotype trait should be used. -#' @param snpChip an integer. Which SNP array should be used. -#' @param simParam an object of \code{\link{SimParam}} -#' @param chromLength an integer. The size of chromosomes in base -#' pairs; assuming all chromosomes are of the same size. -#' -#' @examples -#' \dontrun{ -#' #Create founder haplotypes -#' founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) -#' -#' #Set simulation parameters -#' SP = SimParam$new(founderPop) -#' SP$setSexes(sex = "yes_rand") -#' SP$addTraitA(nQtlPerChr = 10) -#' SP$addSnpChip(nSnpPerChr = 5) -#' -#' #Create population -#' pop = newPop(rawPop = founderPop) -#' pop = setPheno(pop, varE = SP$varA) -#' writePlink(pop, baseName="test") -#' -#' #Test -#' test = read.table(file = "test.ped") -#' #...sex -#' if (!identical(x = c("M", "F")[test[[5]]], y = pop@sex)) { stop() } -#' #...pheno (issues with rounding) -#' # if (!identical(x = test[[6]], y = pop@pheno[, 1])) { stop() } -#' #...genotypes -#' x = test[, -(1:6)] - 1 -#' x[, 1] = x[, 1] + x[, 2] -#' x[, 2] = x[, 3] + x[, 4] -#' x[, 3] = x[, 5] + x[, 6] -#' x[, 4] = x[, 7] + x[, 8] -#' x[, 5] = x[, 9] + x[, 10] -#' y = pullSnpGeno(pop) -#' if (sum(x[, 1:5] - y) != 0) { stop() } -#' } -#' @export -writePlink = function(pop, baseName, trait = 1L, snpChip = 1L, simParam = NULL, - chromLength = 10L^8) { - if (is.null(simParam)) { - simParam = get(x = "SP", envir = .GlobalEnv) - } - if (pop@ploidy != 2L) { - stop(paste0("writePlink() will write ", pop@ploidy, " alleles for each locus!")) - } - - # ---- Map ---- - - # This assumes equal number of markers per chromosome! - map = data.frame(chr = rep(x = 1L:simParam$nChr, - each = simParam$snpChips[[snpChip]]@lociPerChr[1L]), - loc = paste0("SNP_", 1L:simParam$snpChips[[snpChip]]@nLoci), - posGenetic = 0L, - pos = simParam$snpChips[[snpChip]]@lociLoc) - for (chr in 1L:simParam$nChr) { - # chr = 1 - sel = map$chr == chr - map$posGenetic[sel] = simParam$genMap[[chr]][map$pos[sel]] - } - map$pos = round(map$posGenetic * chromLength) - write.table(x = map, file = paste0(baseName, ".map"), - col.names = FALSE, row.names = FALSE, quote = FALSE) - - # ---- Ped ---- - - # First the FAM format, which covers the first 6 columns of the PED format - fam = data.frame(family = rep(x = 1L, times = pop@nInd), - id = as.integer(pop@id), - father = as.integer(pop@father), - mother = as.integer(pop@mother), - sex = 0L, - pheno = 0) - if (!any(pop@sex == "H")) { - fam$sex = (pop@sex == "F") + 1L - } - if (!anyNA(pop@pheno[, trait])) { - fam$pheno = pop@pheno[, trait] - } - - # Select loci on the SNP array - tmp = selectLoci(chr = 1L:simParam$nChr, - inLociPerChr = simParam$snpChips[[snpChip]]@lociPerChr, - inLociLoc = simParam$snpChips[[snpChip]]@lociLoc) - # Add loci alleles to fam and write to file directly from C++ - writePlinkPed(fam = fam, - haplo = getHaplo(geno = pop@geno, - lociPerChr = tmp$lociPerChr, - lociLoc = tmp$lociLoc, - nThreads = simParam$nThreads), - nInd = nrow(fam), - ploidy = pop@ploidy, - nLoc = sum(tmp$lociPerChr), - file = paste0(baseName, ".ped")) -} - #' @title Linear transformation matrix #' #' @description diff --git a/R/pullGeno.R b/R/pullGeno.R index 5df012cc..d56fb7ad 100644 --- a/R/pullGeno.R +++ b/R/pullGeno.R @@ -129,7 +129,7 @@ getSnpMap = function(snpChip=1, sex="A", simParam=NULL){ #Create a data.frame with SNP postions on genetic map output = data.frame(id=getLociNames(snp@lociPerChr, snp@lociLoc, genMap), - chr=rep(1:simParam$nChr,snp@lociPerChr), + chr=rep(names(genMap),snp@lociPerChr), site=snp@lociLoc, pos=do.call("c",snpMap)) return(output) @@ -215,7 +215,7 @@ getQtlMap = function(trait=1, sex="A", simParam=NULL){ #Create a data.frame with QTL positions on genetic map output = data.frame(id=getLociNames(qtl@lociPerChr, qtl@lociLoc, genMap), - chr=rep(1:simParam$nChr,qtl@lociPerChr), + chr=rep(names(genMap),qtl@lociPerChr), site=qtl@lociLoc, pos=do.call("c",qtlMap)) return(output) diff --git a/R/writePlink.R b/R/writePlink.R new file mode 100644 index 00000000..6898c027 --- /dev/null +++ b/R/writePlink.R @@ -0,0 +1,138 @@ +#' @title Writes a Pop-class as PLINK files +#' +#' @description +#' Writes a Pop-class to PLINK PED and MAP files. The arguments +#' for this function were chosen for consistency with +#' \code{\link{RRBLUP2}}. The base pair coordinate will the locus +#' position as stored in AlphaSimR and not an actual base pair +#' position. This is because AlphaSimR doesn't track base pair +#' positions, only relative positions for the loci used in the +#' simulation. +#' +#' @param pop an object of \code{\link{Pop-class}} +#' @param baseName basename for PED and MAP files. +#' @param traits an integer indicating the trait to write, a trait name, or a +#' function of the traits returning a single value. +#' @param use what to use for PLINK's phenotype field. Either phenotypes "pheno", +#' genetic values "gv", estimated breeding values "ebv", breeding values "bv", +#' or random values "rand". +#' @param snpChip an integer indicating which SNP chip genotype +#' to use +#' @param useQtl should QTL genotypes be used instead of a SNP chip. +#' If TRUE, snpChip specifies which trait's QTL to use, and thus these +#' QTL may not match the QTL underlying the phenotype supplied in traits. +#' @param simParam an object of \code{\link{SimParam}} +#' @param ... additional arguments if using a function for +#' traits +#' +#' @examples +#' \dontrun{ +#' #Create founder haplotypes +#' founderPop = quickHaplo(nInd=10, nChr=1, segSites=15) +#' +#' #Set simulation parameters +#' SP = SimParam$new(founderPop) +#' SP$setSexes(sex="yes_rand") +#' SP$addTraitA(nQtlPerChr=10) +#' SP$addSnpChip(nSnpPerChr=5) +#' SP$setVarE(h2=0.5) +#' +#' #Create population +#' pop = newPop(rawPop = founderPop) +#' +#' # Write out PLINK files +#' writePlink(pop, baseName="test") +#' } +#' @export +writePlink = function(pop, baseName, traits=1, use="pheno", + snpChip=1, useQtl=FALSE, simParam=NULL, + ...){ + if(pop@ploidy!=2L){ + stop("writePlink() only supports ploidy=2") + } + + if(is.null(simParam)){ + simParam = get(x="SP", envir=.GlobalEnv) + } + + # Pull "phenotype" data indicated by traits + y = getResponse(pop=pop, trait=traits, use=use, + simParam=simParam, ...) + + # Pull QTL/SNP data indicated by snpChip and useQtl + if(useQtl){ + H1 = pullQtlHaplo(pop=pop, trait=snpChip, + haplo=1, asRaw=TRUE, + simParam=simParam) + + H2 = pullQtlHaplo(pop=pop, trait=snpChip, + haplo=2, asRaw=TRUE, + simParam=simParam) + + map = getQtlMap(trait=snpChip, simParam=simParam) + }else{ + H1 = pullSnpHaplo(pop=pop, snpChip=snpChip, + haplo=1, asRaw=TRUE, + simParam=simParam) + + H2 = pullSnpHaplo(pop=pop, snpChip=snpChip, + haplo=2, asRaw=TRUE, + simParam=simParam) + + map = getSnpMap(snpChip=snpChip, simParam=simParam) + } + + ## Make .ped file + + # Format pop data for a .fam (first columns of .ped) + # Format sex for PLINK + sex = pop@sex + sex[which(sex=="H")] = "0" + sex[which(sex=="M")] = "1" + sex[which(sex=="F")] = "2" + + # Determine within-family ID of father, "0" if not present + father = pop@id[match(pop@father, pop@id)] + father[is.na(father)] = "0" + + # Determine within-family ID of mother, "0" if not present + mother = pop@id[match(pop@mother, pop@id)] + mother[is.na(mother)] = "0" + + fam = rbind(rep("1", pop@nInd), # Family ID + pop@id, # Within-family ID + father, # Within-family ID of father + mother, # Within-family ID of mother + sex, # Sex + as.character(c(y))) # Phenotype + + # Weave together haplotype data for writing to a file with + # the write function (requires a transposed matrix) + H = unname(rbind(t(H1), t(H2))) + H = H[c(matrix(1:nrow(H),nrow=2,byrow=T)),] + + # Free up some memory + rm(H1, H2) + + # Convert to characters for writing to file + H = ifelse(H, "2", "1") + + # Append .fam + H = rbind(fam,H) + + # Write .ped file + write(H, file=paste0(baseName,".ped"), ncolumns=nrow(H)) + + ## Make .map file + map = rbind(map$chr, # Chromosome + map$id, # Variant id + as.character(map$pos*10), # Genetic map position (cM) + as.character(map$site) # Physical map position + ) + + # Write .map file + write(map, file=paste0(baseName,".map"), ncolumns=nrow(map)) + + # Don't return anything + return(invisible()) +} diff --git a/man/SimParam.Rd b/man/SimParam.Rd index 3a909f1d..b5057640 100644 --- a/man/SimParam.Rd +++ b/man/SimParam.Rd @@ -338,43 +338,43 @@ genetic map} \section{Methods}{ \subsection{Public methods}{ \itemize{ -\item \href{#method-new}{\code{SimParam$new()}} -\item \href{#method-setTrackPed}{\code{SimParam$setTrackPed()}} -\item \href{#method-setTrackRec}{\code{SimParam$setTrackRec()}} -\item \href{#method-resetPed}{\code{SimParam$resetPed()}} -\item \href{#method-restrSegSites}{\code{SimParam$restrSegSites()}} -\item \href{#method-setSexes}{\code{SimParam$setSexes()}} -\item \href{#method-addSnpChip}{\code{SimParam$addSnpChip()}} -\item \href{#method-addStructuredSnpChip}{\code{SimParam$addStructuredSnpChip()}} -\item \href{#method-addTraitA}{\code{SimParam$addTraitA()}} -\item \href{#method-addTraitAD}{\code{SimParam$addTraitAD()}} -\item \href{#method-addTraitAG}{\code{SimParam$addTraitAG()}} -\item \href{#method-addTraitADG}{\code{SimParam$addTraitADG()}} -\item \href{#method-addTraitAE}{\code{SimParam$addTraitAE()}} -\item \href{#method-addTraitADE}{\code{SimParam$addTraitADE()}} -\item \href{#method-addTraitAEG}{\code{SimParam$addTraitAEG()}} -\item \href{#method-addTraitADEG}{\code{SimParam$addTraitADEG()}} -\item \href{#method-manAddTrait}{\code{SimParam$manAddTrait()}} -\item \href{#method-importTrait}{\code{SimParam$importTrait()}} -\item \href{#method-switchTrait}{\code{SimParam$switchTrait()}} -\item \href{#method-removeTrait}{\code{SimParam$removeTrait()}} -\item \href{#method-setVarE}{\code{SimParam$setVarE()}} -\item \href{#method-setCorE}{\code{SimParam$setCorE()}} -\item \href{#method-rescaleTraits}{\code{SimParam$rescaleTraits()}} -\item \href{#method-setRecombRatio}{\code{SimParam$setRecombRatio()}} -\item \href{#method-switchGenMap}{\code{SimParam$switchGenMap()}} -\item \href{#method-switchFemaleMap}{\code{SimParam$switchFemaleMap()}} -\item \href{#method-switchMaleMap}{\code{SimParam$switchMaleMap()}} -\item \href{#method-addToRec}{\code{SimParam$addToRec()}} -\item \href{#method-ibdHaplo}{\code{SimParam$ibdHaplo()}} -\item \href{#method-updateLastId}{\code{SimParam$updateLastId()}} -\item \href{#method-addToPed}{\code{SimParam$addToPed()}} -\item \href{#method-clone}{\code{SimParam$clone()}} +\item \href{#method-SimParam-new}{\code{SimParam$new()}} +\item \href{#method-SimParam-setTrackPed}{\code{SimParam$setTrackPed()}} +\item \href{#method-SimParam-setTrackRec}{\code{SimParam$setTrackRec()}} +\item \href{#method-SimParam-resetPed}{\code{SimParam$resetPed()}} +\item \href{#method-SimParam-restrSegSites}{\code{SimParam$restrSegSites()}} +\item \href{#method-SimParam-setSexes}{\code{SimParam$setSexes()}} +\item \href{#method-SimParam-addSnpChip}{\code{SimParam$addSnpChip()}} +\item \href{#method-SimParam-addStructuredSnpChip}{\code{SimParam$addStructuredSnpChip()}} +\item \href{#method-SimParam-addTraitA}{\code{SimParam$addTraitA()}} +\item \href{#method-SimParam-addTraitAD}{\code{SimParam$addTraitAD()}} +\item \href{#method-SimParam-addTraitAG}{\code{SimParam$addTraitAG()}} +\item \href{#method-SimParam-addTraitADG}{\code{SimParam$addTraitADG()}} +\item \href{#method-SimParam-addTraitAE}{\code{SimParam$addTraitAE()}} +\item \href{#method-SimParam-addTraitADE}{\code{SimParam$addTraitADE()}} +\item \href{#method-SimParam-addTraitAEG}{\code{SimParam$addTraitAEG()}} +\item \href{#method-SimParam-addTraitADEG}{\code{SimParam$addTraitADEG()}} +\item \href{#method-SimParam-manAddTrait}{\code{SimParam$manAddTrait()}} +\item \href{#method-SimParam-importTrait}{\code{SimParam$importTrait()}} +\item \href{#method-SimParam-switchTrait}{\code{SimParam$switchTrait()}} +\item \href{#method-SimParam-removeTrait}{\code{SimParam$removeTrait()}} +\item \href{#method-SimParam-setVarE}{\code{SimParam$setVarE()}} +\item \href{#method-SimParam-setCorE}{\code{SimParam$setCorE()}} +\item \href{#method-SimParam-rescaleTraits}{\code{SimParam$rescaleTraits()}} +\item \href{#method-SimParam-setRecombRatio}{\code{SimParam$setRecombRatio()}} +\item \href{#method-SimParam-switchGenMap}{\code{SimParam$switchGenMap()}} +\item \href{#method-SimParam-switchFemaleMap}{\code{SimParam$switchFemaleMap()}} +\item \href{#method-SimParam-switchMaleMap}{\code{SimParam$switchMaleMap()}} +\item \href{#method-SimParam-addToRec}{\code{SimParam$addToRec()}} +\item \href{#method-SimParam-ibdHaplo}{\code{SimParam$ibdHaplo()}} +\item \href{#method-SimParam-updateLastId}{\code{SimParam$updateLastId()}} +\item \href{#method-SimParam-addToPed}{\code{SimParam$addToPed()}} +\item \href{#method-SimParam-clone}{\code{SimParam$clone()}} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-new}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-new}{}}} \subsection{Method \code{new()}}{ Starts the process of building a new simulation by creating a new SimParam object and assigning a founder @@ -409,8 +409,8 @@ SP = SimParam$new(founderPop) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-setTrackPed}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-setTrackPed}{}}} \subsection{Method \code{setTrackPed()}}{ Sets pedigree tracking for the simulation. By default pedigree tracking is turned off. When turned on, @@ -447,8 +447,8 @@ SP$setTrackPed(TRUE) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-setTrackRec}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-setTrackRec}{}}} \subsection{Method \code{setTrackRec()}}{ Sets recombination tracking for the simulation. By default recombination tracking is turned off. When turned @@ -485,8 +485,8 @@ SP$setTrackRec(TRUE) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-resetPed}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-resetPed}{}}} \subsection{Method \code{resetPed()}}{ Resets the internal lastId, the pedigree and recombination tracking (if in use) to the @@ -527,8 +527,8 @@ pop2@id # 1:10 } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-restrSegSites}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-restrSegSites}{}}} \subsection{Method \code{restrSegSites()}}{ Sets restrictions on which segregating sites can serve as SNP and/or QTL. @@ -574,8 +574,8 @@ SP$restrSegSites(minQtlPerChr=5, minSnpPerChr=5) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-setSexes}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-setSexes}{}}} \subsection{Method \code{setSexes()}}{ Changes how sexes are determined in the simulation. The default sexes is "no", indicating all individuals are hermaphrodites. @@ -615,8 +615,8 @@ SP$setSexes("yes_sys") } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addSnpChip}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addSnpChip}{}}} \subsection{Method \code{addSnpChip()}}{ Randomly assigns eligible SNPs to a SNP chip \subsection{Usage}{ @@ -654,8 +654,8 @@ SP$addSnpChip(10) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addStructuredSnpChip}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addStructuredSnpChip}{}}} \subsection{Method \code{addStructuredSnpChip()}}{ Randomly selects the number of snps in structure and then assigns them to chips based on structure @@ -679,8 +679,8 @@ ignored. Only set to TRUE if you know what you are doing.} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addTraitA}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addTraitA}{}}} \subsection{Method \code{addTraitA()}}{ Randomly assigns eligible QTLs for one or more additive traits. If simulating more than one trait, all traits will be pleiotrophic @@ -735,8 +735,8 @@ SP$addTraitA(10) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addTraitAD}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addTraitAD}{}}} \subsection{Method \code{addTraitAD()}}{ Randomly assigns eligible QTLs for one or more traits with dominance. If simulating more than one trait, all traits will be pleiotrophic @@ -804,8 +804,8 @@ SP$addTraitAD(10, meanDD=0.5) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addTraitAG}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addTraitAG}{}}} \subsection{Method \code{addTraitAG()}}{ Randomly assigns eligible QTLs for one ore more additive GxE traits. If simulating more than one trait, all traits will be pleiotrophic @@ -869,8 +869,8 @@ SP$addTraitAG(10, varGxE=2) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addTraitADG}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addTraitADG}{}}} \subsection{Method \code{addTraitADG()}}{ Randomly assigns eligible QTLs for a trait with dominance and GxE. \subsection{Usage}{ @@ -945,8 +945,8 @@ SP$addTraitADG(10, meanDD=0.5, varGxE=2) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addTraitAE}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addTraitAE}{}}} \subsection{Method \code{addTraitAE()}}{ Randomly assigns eligible QTLs for one or more additive and epistasis traits. If simulating more than one trait, all traits will be pleiotrophic @@ -1008,8 +1008,8 @@ founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addTraitADE}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addTraitADE}{}}} \subsection{Method \code{addTraitADE()}}{ Randomly assigns eligible QTLs for one or more traits with dominance and epistasis. If simulating more than one trait, all traits will be pleiotrophic @@ -1084,8 +1084,8 @@ SP$addTraitADE(10) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addTraitAEG}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addTraitAEG}{}}} \subsection{Method \code{addTraitAEG()}}{ Randomly assigns eligible QTLs for one or more additive and epistasis GxE traits. If simulating more than one trait, all traits will be pleiotrophic @@ -1160,8 +1160,8 @@ SP$addTraitAEG(10, varGxE=2) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addTraitADEG}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addTraitADEG}{}}} \subsection{Method \code{addTraitADEG()}}{ Randomly assigns eligible QTLs for a trait with dominance, epistasis and GxE. @@ -1244,8 +1244,8 @@ SP$addTraitADEG(10, meanDD=0.5, varGxE=2) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-manAddTrait}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-manAddTrait}{}}} \subsection{Method \code{manAddTrait()}}{ Manually add a new trait to the simulation. Trait must be formatted as a \code{\link{LociMap-class}}. If the @@ -1269,8 +1269,8 @@ ignored. Only set to TRUE if you know what you are doing} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-importTrait}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-importTrait}{}}} \subsection{Method \code{importTrait()}}{ Manually add a new trait(s) to the simulation. Unlike the manAddTrait function, this function does not require @@ -1314,8 +1314,8 @@ ignored. Only set to TRUE if you know what you are doing} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-switchTrait}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-switchTrait}{}}} \subsection{Method \code{switchTrait()}}{ Switch a trait in the simulation. \subsection{Usage}{ @@ -1339,8 +1339,8 @@ ignored. Only set to TRUE if you know what you are doing} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-removeTrait}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-removeTrait}{}}} \subsection{Method \code{removeTrait()}}{ Remove a trait from the simulation \subsection{Usage}{ @@ -1359,8 +1359,8 @@ ignored. Only set to TRUE if you know what you are doing} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-setVarE}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-setVarE}{}}} \subsection{Method \code{setVarE()}}{ Defines a default value for error variances in the simulation. @@ -1395,8 +1395,8 @@ SP$setVarE(h2=0.5) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-setCorE}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-setCorE}{}}} \subsection{Method \code{setCorE()}}{ Defines a correlation structure for default error variances. You must call \code{setVarE} first to define @@ -1430,8 +1430,8 @@ SP$setCorE(E) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-rescaleTraits}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-rescaleTraits}{}}} \subsection{Method \code{rescaleTraits()}}{ Linearly scales all traits to achieve desired values of means and variances in the founder population. @@ -1485,8 +1485,8 @@ meanG(pop) } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-setRecombRatio}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-setRecombRatio}{}}} \subsection{Method \code{setRecombRatio()}}{ Set the relative recombination rates between males and females. This allows for sex-specific recombination rates, @@ -1519,8 +1519,8 @@ SP$setRecombRatio(2) #Twice as much recombination in females } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-switchGenMap}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-switchGenMap}{}}} \subsection{Method \code{switchGenMap()}}{ Replaces existing genetic map. \subsection{Usage}{ @@ -1542,8 +1542,8 @@ be metacentric.} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-switchFemaleMap}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-switchFemaleMap}{}}} \subsection{Method \code{switchFemaleMap()}}{ Replaces existing female genetic map. \subsection{Usage}{ @@ -1565,8 +1565,8 @@ be metacentric.} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-switchMaleMap}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-switchMaleMap}{}}} \subsection{Method \code{switchMaleMap()}}{ Replaces existing male genetic map. \subsection{Usage}{ @@ -1588,8 +1588,8 @@ be metacentric.} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addToRec}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addToRec}{}}} \subsection{Method \code{addToRec()}}{ For internal use only. \subsection{Usage}{ @@ -1617,8 +1617,8 @@ For internal use only. } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-ibdHaplo}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-ibdHaplo}{}}} \subsection{Method \code{ibdHaplo()}}{ For internal use only. \subsection{Usage}{ @@ -1634,8 +1634,8 @@ For internal use only. } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-updateLastId}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-updateLastId}{}}} \subsection{Method \code{updateLastId()}}{ For internal use only. \subsection{Usage}{ @@ -1651,8 +1651,8 @@ For internal use only. } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-addToPed}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-addToPed}{}}} \subsection{Method \code{addToPed()}}{ For internal use only. \subsection{Usage}{ @@ -1676,8 +1676,8 @@ For internal use only. } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-clone}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-SimParam-clone}{}}} \subsection{Method \code{clone()}}{ The objects of this class are cloneable with this method. \subsection{Usage}{ diff --git a/man/writePlink.Rd b/man/writePlink.Rd index ce4d540e..9e1efeb7 100644 --- a/man/writePlink.Rd +++ b/man/writePlink.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/misc.R +% Please edit documentation in R/writePlink.R \name{writePlink} \alias{writePlink} \title{Writes a Pop-class as PLINK files} @@ -7,59 +7,63 @@ writePlink( pop, baseName, - trait = 1L, - snpChip = 1L, + traits = 1, + use = "pheno", + snpChip = 1, + useQtl = FALSE, simParam = NULL, - chromLength = 10L^8 + ... ) } \arguments{ \item{pop}{an object of \code{\link{Pop-class}}} -\item{baseName}{a character. Basename of PED and MAP files.} +\item{baseName}{basename for PED and MAP files.} -\item{trait}{an integer. Which phenotype trait should be used.} +\item{traits}{an integer indicating the trait to write, a trait name, or a +function of the traits returning a single value.} -\item{snpChip}{an integer. Which SNP array should be used.} +\item{use}{what to use for PLINK's phenotype field. Either phenotypes "pheno", +genetic values "gv", estimated breeding values "ebv", breeding values "bv", +or random values "rand".} + +\item{snpChip}{an integer indicating which SNP chip genotype +to use} + +\item{useQtl}{should QTL genotypes be used instead of a SNP chip. +If TRUE, snpChip specifies which trait's QTL to use, and thus these +QTL may not match the QTL underlying the phenotype supplied in traits.} \item{simParam}{an object of \code{\link{SimParam}}} -\item{chromLength}{an integer. The size of chromosomes in base -pairs; assuming all chromosomes are of the same size.} +\item{...}{additional arguments if using a function for +traits} } \description{ -Writes a Pop-class as PLINK PED and MAP files +Writes a Pop-class to PLINK PED and MAP files. The arguments +for this function were chosen for consistency with +\code{\link{RRBLUP2}}. The base pair coordinate will the locus +position as stored in AlphaSimR and not an actual base pair +position. This is because AlphaSimR doesn't track base pair +positions, only relative positions for the loci used in the +simulation. } \examples{ \dontrun{ #Create founder haplotypes -founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) +founderPop = quickHaplo(nInd=10, nChr=1, segSites=15) #Set simulation parameters SP = SimParam$new(founderPop) -SP$setSexes(sex = "yes_rand") -SP$addTraitA(nQtlPerChr = 10) -SP$addSnpChip(nSnpPerChr = 5) +SP$setSexes(sex="yes_rand") +SP$addTraitA(nQtlPerChr=10) +SP$addSnpChip(nSnpPerChr=5) +SP$setVarE(h2=0.5) #Create population pop = newPop(rawPop = founderPop) -pop = setPheno(pop, varE = SP$varA) -writePlink(pop, baseName="test") -#Test -test = read.table(file = "test.ped") -#...sex -if (!identical(x = c("M", "F")[test[[5]]], y = pop@sex)) { stop() } -#...pheno (issues with rounding) -# if (!identical(x = test[[6]], y = pop@pheno[, 1])) { stop() } -#...genotypes -x = test[, -(1:6)] - 1 -x[, 1] = x[, 1] + x[, 2] -x[, 2] = x[, 3] + x[, 4] -x[, 3] = x[, 5] + x[, 6] -x[, 4] = x[, 7] + x[, 8] -x[, 5] = x[, 9] + x[, 10] -y = pullSnpGeno(pop) -if (sum(x[, 1:5] - y) != 0) { stop() } +# Write out PLINK files +writePlink(pop, baseName="test") } } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index e4da679a..c3089404 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -790,21 +790,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// writePlinkPed -void writePlinkPed(const Rcpp::DataFrame& fam, const arma::Mat& haplo, const arma::uword& nInd, const arma::uword& ploidy, const arma::uword& nLoc, const Rcpp::String& file); -RcppExport SEXP _AlphaSimR_writePlinkPed(SEXP famSEXP, SEXP haploSEXP, SEXP nIndSEXP, SEXP ploidySEXP, SEXP nLocSEXP, SEXP fileSEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const Rcpp::DataFrame& >::type fam(famSEXP); - Rcpp::traits::input_parameter< const arma::Mat& >::type haplo(haploSEXP); - Rcpp::traits::input_parameter< const arma::uword& >::type nInd(nIndSEXP); - Rcpp::traits::input_parameter< const arma::uword& >::type ploidy(ploidySEXP); - Rcpp::traits::input_parameter< const arma::uword& >::type nLoc(nLocSEXP); - Rcpp::traits::input_parameter< const Rcpp::String& >::type file(fileSEXP); - writePlinkPed(fam, haplo, nInd, ploidy, nLoc, file); - return R_NilValue; -END_RCPP -} // MaCS Rcpp::List MaCS(Rcpp::String args, arma::uvec maxSites, bool inbred, arma::uword ploidy, int nThreads, Rcpp::StringVector seed); RcppExport SEXP _AlphaSimR_MaCS(SEXP argsSEXP, SEXP maxSitesSEXP, SEXP inbredSEXP, SEXP ploidySEXP, SEXP nThreadsSEXP, SEXP seedSEXP) { @@ -872,7 +857,6 @@ static const R_CallMethodDef CallEntries[] = { {"_AlphaSimR_calcCoef", (DL_FUNC) &_AlphaSimR_calcCoef, 2}, {"_AlphaSimR_getNumThreads", (DL_FUNC) &_AlphaSimR_getNumThreads, 0}, {"_AlphaSimR_packHaplo", (DL_FUNC) &_AlphaSimR_packHaplo, 3}, - {"_AlphaSimR_writePlinkPed", (DL_FUNC) &_AlphaSimR_writePlinkPed, 6}, {"_AlphaSimR_MaCS", (DL_FUNC) &_AlphaSimR_MaCS, 6}, {NULL, NULL, 0} }; diff --git a/src/plink.cpp b/src/plink.cpp deleted file mode 100644 index e3531be9..00000000 --- a/src/plink.cpp +++ /dev/null @@ -1,56 +0,0 @@ -#include "alphasimr.h" - -// Write PLINK PED file -// [[Rcpp::export]] -void writePlinkPed(const Rcpp::DataFrame & fam, - const arma::Mat & haplo, - const arma::uword & nInd, - const arma::uword & ploidy, - const arma::uword & nLoc, - const Rcpp::String & file) { - - // FAM file stuff - Rcpp::IntegerVector family = fam["family"]; - Rcpp::IntegerVector id = fam["id"]; - Rcpp::IntegerVector father = fam["father"]; - Rcpp::IntegerVector mother = fam["mother"]; - Rcpp::IntegerVector sex = fam["sex"]; - Rcpp::NumericVector pheno = fam["pheno"]; - - std::ofstream plinkPed; - plinkPed.open(file, std::ofstream::trunc); - arma::Col alleles(nLoc * ploidy); - for (arma::uword ind = 0; ind < nInd; ind++) { - // Family & Individual data - plinkPed << family[ind] << " " << id[ind] << " " << father[ind] << " " - << mother[ind] << " " << sex[ind] << " " << pheno[ind]; - // Haplotype data - // ... first gather locus alleles from all chromosomes - arma::uword indLastChrom = (ind + 1) * ploidy; - for (arma::uword p = ploidy; p > 0; p--) { - arma::uword indChrom = indLastChrom - p; - for (arma::uword loc = 0; loc < nLoc; loc++) { - arma::uword locLastChrom = (loc + 1) * ploidy; - /* - std::cout << " Ind " << ind + 1 - << " Chrom " << ploidy - p + 1 - << " Locus " << loc + 1 - << " Haplo row " << indChrom - << " Allele " << haplo(indChrom, loc) + 49 - << " Allele " << char(haplo(indChrom, loc) + 49) - << "\n"; - */ - // haplo stores 0 and 1 as "unsigned char" [0 - 255] - // to get ASCII characters 1 and 2 we need values 49 and 50 http://www.asciitable.com - // and cast these "unsigned char" values as "char" - alleles(locLastChrom - p) = char(haplo(indChrom, loc) + 49); - } - } - // ... now print locus alleles from all chromosomes - for (arma::uword loc = 0; loc < (nLoc * ploidy); loc++) { - plinkPed << " " << alleles(loc); - } - plinkPed << "\n"; - } - plinkPed.close(); -} diff --git a/tests/testthat/test-importData.R b/tests/testthat/test-importData.R new file mode 100644 index 00000000..49a9ed39 --- /dev/null +++ b/tests/testthat/test-importData.R @@ -0,0 +1,47 @@ +context("importData") + +test_that("importTrait",{ + # Create haplotype data + haplo = rbind(c(1,1,0,1,0), + c(1,1,0,1,0), + c(0,1,1,0,0), + c(0,1,1,0,0)) + colnames(haplo) = letters[1:5] + + # Create genetic map + genMap = data.frame(markerName=letters[1:5], + chromosome=c(1,1,1,2,2), + position=c(0,0.5,1,0.15,0.4)) + + # Create pedigree + ped = data.frame(id=c("a","b"), + mother=c(0,0), + father=c(0,0)) + + # Generate an external trait + myTrait = data.frame(marker = c("a","c","d"), + a = c(1,-1,1)) + + founderPop = importHaplo(haplo = haplo, + genMap = genMap, + ploidy = 2L, + ped = ped) + + SP = SimParam$new(founderPop=founderPop) + + # Import trait + SP$importTrait(markerNames = myTrait$marker, + addEff = myTrait$a, + name = "myTrait") + + expect_equal(SP$traits[[1]]@addEff, myTrait$a, tolerance=1e-6) + + expect_equal(SP$traits[[1]]@intercept, 0, tolerance=1e-6) + + pop = newPop(founderPop,simParam=SP) + + expect_equal(unname(pop@gv[1,1]), 3, tolerance=1e-6) + + expect_equal(unname(pop@gv[2,1]), -3, tolerance=1e-6) +}) +