Skip to content

Commit

Permalink
Merge pull request #64 from gaynorr/devel
Browse files Browse the repository at this point in the history
Devel
  • Loading branch information
gaynorr committed Jul 6, 2022
2 parents f0227c2 + fe6426c commit d91c02b
Show file tree
Hide file tree
Showing 12 changed files with 333 additions and 317 deletions.
6 changes: 3 additions & 3 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -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 = "[email protected]",
role = c("aut", "cre"), comment = c(ORCID = "0000-0003-0558-6656")),
person("Gregor", "Gorjanc", role = "ctb",
Expand Down Expand Up @@ -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
6 changes: 6 additions & 0 deletions NEWS
Original file line number Diff line number Diff line change
@@ -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
Expand Down
11 changes: 6 additions & 5 deletions R/Class-SimParam.R
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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))
Expand All @@ -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){
Expand Down
4 changes: 0 additions & 4 deletions R/RcppExports.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
}
Expand Down
104 changes: 0 additions & 104 deletions R/misc.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions R/pullGeno.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
138 changes: 138 additions & 0 deletions R/writePlink.R
Original file line number Diff line number Diff line change
@@ -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())
}
Loading

0 comments on commit d91c02b

Please sign in to comment.