Nothing
convert2bcsft <- function(cross, BC.gen = 0, F.gen = 0, estimate.map = TRUE,
error.prob=0.0001, map.function=c("haldane","kosambi","c-f","morgan"),
verbose=FALSE)
{
cross.class <- crosstype(cross)
if(cross.class=="bcsft") cross.class <- "f2"
if((cross.class %in% c("bc","f2"))) {
class(cross) <- c("bcsft", "cross")
## If BC.gen = 0 and F.gen = 0, then set to BC1F0 (bc) or BC0F2 (f2).
if(cross.class == "bc" & F.gen > 0) {
stop("input cross has only 2 genotypes--cannot have F.gen > 0")
if(BC.gen == 0)
BC.gen <- 1
}
if(cross.class == "f2") {
if(F.gen == 0) {
if(BC.gen == 0)
F.gen <- 2
else
stop("input cross has 3 genotypes--cannot have F.gen = 0")
}
}
if(BC.gen < 0 | F.gen < 0)
stop("BC.gen and F.gen cannot be negative")
attr(cross, "scheme") <- c(BC.gen, F.gen)
cross
}
else stop("cross object has to be of class bc or f2 to be converted to bcsft")
# re-estimate map?
if(estimate.map) {
cat(" --Estimating genetic map\n")
newmap <- est.map(cross, error.prob=error.prob, map.function=map.function, verbose=verbose)
cross <- replace.map(cross, newmap)
}
cross
}
read.cross.bcsft <- function(..., BC.gen = 0, F.gen = 0, cross = NULL, force.bcsft = FALSE,
estimate.map=TRUE)
{
## Must specify s = BC.gen and t = F.gen.
## Later: Could import in clever way from qtlcart? See qtlcart_io.R and their software.
## Make sure we only estimate map once!
if(is.null(cross)) # Estimate map at end of this routine (called read.cross.bcsft directly).
cross <- read.cross(..., estimate.map = FALSE)
else # Estimate map in parent read.cross() call (read.cross.bcsft is pass-through from read.cross).
estimate.map <- FALSE
force.bcsft <- force.bcsft | (BC.gen > 0 | F.gen > 0)
if((crosstype(cross) %in% c("bc","f2","bcsft")) && force.bcsft) {
# deal with ... args
dots <- list(...)
if("verbose" %in% names(dots)) verbose <- dots$verbose
else verbose <- TRUE
if("error.prob" %in% names(dots)) error.prob <- dots$error.prob
else error.prob <- 0.0001
if("map.function" %in% names(dots)) map.function <- dots$map.function
else map.function <- "haldane"
cross <- convert2bcsft(cross, BC.gen, F.gen, estimate.map = estimate.map,
error.prob=error.prob, map.function=map.function,
verbose=verbose)
}
cross
}
######################################################################
sim.cross.bcsft <- function(map,model,n.ind,error.prob,missing.prob,
partial.missing.prob,keep.errorind,
m,p,map.function,
cross.scheme)
{
if(map.function=="kosambi") mf <- mf.k
else if(map.function=="c-f") mf <- mf.cf
else if(map.function=="morgan") mf <- mf.m
else mf <- mf.h
if(any(sapply(map,is.matrix)))
stop("Map must not be sex-specific.")
## cross.scheme = c(s,t) for bcsft.
if(missing(cross.scheme))
stop("must specify cross.scheme for bcsft")
if(length(cross.scheme) != 2)
stop("cross.scheme for bcsft must have 2 values")
cross.scheme <- round(cross.scheme)
if(min(cross.scheme) < 0)
stop("cross.scheme for bcsft must have 2 non-negative integers")
n.eff <- 3 + (cross.scheme[2] > 0)
## chromosome types
chr.type <- sapply(map,chrtype)
n.chr <- length(map)
if(is.null(model)) n.qtl <- 0
else {
if(!((!is.matrix(model) && length(model) == n.eff) ||
(is.matrix(model) && ncol(model) == n.eff))) {
stop(paste("Model must be a matrix with ", n.eff, " columns (chr, pos and effect",
ifelse(n.eff == 4, "s", ""), ").", sep = ""))
}
if(!is.matrix(model)) model <- rbind(model)
n.qtl <- nrow(model)
if(any(model[,1] < 0 | model[,1] > n.chr))
stop("Chromosome indicators in model matrix out of range.")
model[,2] <- model[,2]+1e-14 # so QTL not on top of marker
}
# if any QTLs, place qtls on map
if(n.qtl > 0) {
for(i in 1:n.qtl) {
temp <- map[[model[i,1]]]
if(model[i,2] < min(temp)) {
temp <- c(model[i,2],temp)
names(temp)[1] <- paste("QTL",i,sep="")
}
else if(model[i,2] > max(temp)) {
temp <- c(temp,model[i,2])
names(temp)[length(temp)] <- paste("QTL",i,sep="")
}
else {
j <- max((seq(along=temp))[temp < model[i,2]])
temp <- c(temp[1:j],model[i,2],temp[(j+1):length(temp)])
names(temp)[j+1] <- paste("QTL",i,sep="")
}
map[[model[i,1]]] <- temp
}
}
geno <- vector("list", n.chr)
names(geno) <- names(map)
n.mar <- sapply(map,length)
mar.names <- lapply(map,names)
BC.gen <- cross.scheme[1]
F.gen <- cross.scheme[2] - (BC.gen == 0)
for(i in 1:n.chr) {
# simulate genotype data
bcallele1 <- sim.bcg(n.ind, map[[i]], m, p, map.function) - 1
## BCs: multiply independent instances of meiosis together.
if(BC.gen > 0) {
if(BC.gen > 1) for(j in seq(2, BC.gen))
bcallele1 <- bcallele1 * (sim.bcg(n.ind, map[[i]], m, p, map.function) - 1)
}
if(F.gen == 0) ## BCs only.
thedata <- bcallele1 + 1
else {
if(chr.type[i] == "X") {
if(F.gen > 1) for(j in seq(F.gen))
bcallele1 <- bcallele1 * (sim.bcg(n.ind, map[[i]], m, p, map.function) - 1)
thedata <- bcallele1 + 1
}
else { ## chr.type[i] != "X"
if(BC.gen > 0) { ## Two unique alleles from BC(s).
bcallele2 <- bcallele1 * (sim.bcg(n.ind, map[[i]], m, p, map.function) - 1)
bcallele1 <- bcallele1 * (sim.bcg(n.ind, map[[i]], m, p, map.function) - 1)
}
else ## Starting from F(1) with two unique alleles.
bcallele2 <- sim.bcg(n.ind, map[[i]], m, p, map.function) - 1
if(F.gen > 1) for(j in seq(2, F.gen)) {
## need two instances.
allelemask1 <- sim.bcg(n.ind, map[[i]], m, p, map.function) - 1
allelemask2 <- sim.bcg(n.ind, map[[i]], m, p, map.function) - 1
bcallele1 <- bcallele1 * allelemask1 + bcallele2 * (1 - allelemask1)
bcallele2 <- bcallele2 * allelemask2 + bcallele1 * (1 - allelemask2)
}
thedata <- bcallele1 + bcallele2 + 1
}
}
dimnames(thedata) <- list(NULL,mar.names[[i]])
geno[[i]] <- list(data = thedata, map = map[[i]])
class(geno[[i]]) <- chr.type[i]
class(geno[[i]]$map) <- NULL
} # end loop over chromosomes
# simulate phenotypes
pheno <- rnorm(n.ind,0,1)
if(n.qtl > 0) {
# find QTL positions in genotype data
QTL.chr <- QTL.loc <- NULL
for(i in 1:n.chr) {
o <- grep("^QTL[0-9]+",mar.names[[i]])
if(length(o)>0) {
QTL.chr <- c(QTL.chr,rep(i,length(o)))
QTL.loc <- c(QTL.loc,o)
}
}
# incorporate QTL effects
for(i in 1:n.qtl) {
QTL.geno <- geno[[QTL.chr[i]]]$data[,QTL.loc[i]]
pheno[QTL.geno==2] <- pheno[QTL.geno==2] + model[i,n.eff]
if(n.eff == 4) {
pheno[QTL.geno==1] <- pheno[QTL.geno==1] - model[i,3]
pheno[QTL.geno==3] <- pheno[QTL.geno==3] + model[i,3]
}
}
} # end simulate phenotype
n.mar <- sapply(geno, function(a) length(a$map))
# add errors
if(error.prob > 0) {
for(i in 1:n.chr) {
if(chr.type[i]=="X") {
a <- sample(0:1,n.mar[i]*n.ind,replace=TRUE,
prob=c(1-error.prob,error.prob))
geno[[i]]$data[a == 1] <- 3 - geno[[i]]$data[a == 1]
}
else {
a <- sample(0:2,n.mar[i]*n.ind,replace=TRUE,
prob=c(1-error.prob,error.prob/2,error.prob/2))
if(any(a>0 & geno[[i]]$data==1))
geno[[i]]$data[a>0 & geno[[i]]$data==1] <-
(geno[[i]]$data+a)[a>0 & geno[[i]]$data==1]
if(any(a>0 & geno[[i]]$data==2)) {
geno[[i]]$data[a>0 & geno[[i]]$data==2] <-
(geno[[i]]$data+a)[a>0 & geno[[i]]$data==2]
geno[[i]]$data[geno[[i]]$data>3] <- 1
}
if(any(a>0 & geno[[i]]$data==3))
geno[[i]]$data[a>0 & geno[[i]]$data==3] <-
(geno[[i]]$data-a)[a>0 & geno[[i]]$data==3]
}
if(keep.errorind) {
errors <- matrix(0,n.ind,n.mar[i])
errors[a>0] <- 1
colnames(errors) <- colnames(geno[[i]]$data)
geno[[i]]$errors <- errors
}
} # end loop over chromosomes
} # end simulate genotyping errors
## add partial missing
if(partial.missing.prob > 0) {
for(i in 1:n.chr) {
if(chr.type[i] != "X") {
o <- sample(c(TRUE,FALSE),n.mar[i],replace=TRUE,
prob=c(partial.missing.prob,1-partial.missing.prob))
if(any(o)) {
o2 <- grep("^QTL[0-9]+",mar.names[[i]])
if(length(o2)>0)
x <- geno[[i]]$data[,o2]
m <- (1:n.mar[i])[o]
for(j in m) {
if(runif(1) < 0.5)
geno[[i]]$data[geno[[i]]$data[,j]==1 | geno[[i]]$data[,j]==2,j] <- 4
else
geno[[i]]$data[geno[[i]]$data[,j]==3 | geno[[i]]$data[,j]==2,j] <- 5
}
if(length(o2)>0)
geno[[i]]$data[,o2] <- x
}
}
} # end loop over chromosomes
} # end simulate partially missing data
# add missing
if(missing.prob > 0) {
for(i in 1:n.chr) {
o <- grep("^QTL[0-9]+",mar.names[[i]])
if(length(o)>0)
x <- geno[[i]]$data[,o]
geno[[i]]$data[sample(c(TRUE,FALSE),n.mar[i]*n.ind,replace=TRUE,
prob=c(missing.prob,1-missing.prob))] <- NA
if(length(o)>0)
geno[[i]]$data[,o] <- x
}
}
pheno <- data.frame(phenotype=pheno)
cross <- list(geno=geno,pheno=pheno)
class(cross) <- c("bcsft","cross")
attr(cross, "scheme") <- cross.scheme
cross
}
## End bcsft.R
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.