pValue <- function(dist = seq(25,40, by = 5), pop.size = 100:500, map.function = "kosambi", LOD = FALSE){ colour_hue <- function(n) { hues = seq(15, 375, length = n + 1) hcl(h = hues, l = 65, c = 100)[1:n] } if(max(dist) > 100) stop("Genetic distance exceeds maximum allowed.") mf <- switch(map.function, kosambi = mf.k, haldane = mf.h, morgan = mf.m, cf = mf.cf) rf <- mf(dist) val <- list() for(i in 1:length(rf)){ rec <- rf[i]*pop.size if(LOD){ val[[i]] <- (pop.size - rec)*log10(2*(1 - rf[i])) + rec*log10(2*rf[i]) ylab <- list("LOD score of linkage", cex = 1.5) } else { val[[i]] <- -log10(exp(-2*((pop.size/2 - rec)^2)/pop.size)) ylab <- list(expression(paste("-log10 ",epsilon, sep = "")), cex = 1.5) } } dat <- cbind.data.frame(val = unlist(val)) dat$pop.size <- rep(pop.size, length(dist)) dat$dist <- rep(dist, each = length(pop.size)) cols <- colour_hue(length(dist)) labs <- paste(dist, " cM", sepm = "") print(xyplot(val ~ pop.size, type = "l", data = dat, groups = dat$dist, lwd = 2, col = cols, xlab = "Number of genotypes in population", ylab = ylab, key = list(x = 0.05, y = 0.9, text = list(labs, cex = 2), lines = list(col = cols, lwd = 3)))) } createRQTL.map <- function(pop.size, no.chr = 1, no.mark = 200, chr.length = 200, random = TRUE, mf = "kosambi", QTL = NULL, error.rate = NULL, miss.rate = NULL, pop = "DH"){ if(no.chr > 1){ if(length(no.mark) != no.chr) if(length(no.mark) == 1) no.mark <- rep(no.mark, no.chr) else stop("Need to know how many markers in each linkage group.") if(length(chr.length) != no.chr) if(length(chr.length) == 1) chr.length <- rep(chr.length, no.chr) else stop("Need to know how long you want each linkage group.") } gn <- paste("g", 1:pop.size, sep = "") mapl <- list() mapl$geno <- list() for(j in 1:no.chr){ if(random) mp <- sort(sample(seq(0.01, chr.length[j] - 0.01, by = 0.01), no.mark[j] - 2, replace = FALSE)) else mp <- seq(0, chr.length[j], length.out = no.mark[j])[2:(no.mark[j] - 1)] if(!is.null(QTL[[j]])){ mp <- unique(sort(c(mp, QTL[[j]]))) no.mark[j] <- length(mp) + 2 } dist <- c(0, diff(c(0, mp, chr.length[j]))) mm <- matrix(0, pop.size, no.mark[j]) # r <- mm if(mf == "kosambi") rfrac <- tanh(dist/50)/2 else rfrac <- (1 - exp(-2 * dist/100))/2 if(pop %in% c("DH","BC")){ mm[, 1] <- sample(0:1, pop.size, T) for(i in 2:no.mark[j]) { r <- sample(0:1, pop.size, T, prob = c(1 - rfrac[i], rfrac[i])) rind <- (1:pop.size)[r == 1] mt <- mm[,i - 1] mt[rind] <- as.numeric(!as.logical(mt[rind])) mm[, i] <- mt } cl <- "bc" } else if (length(grep("F", pop))) { gen <- as.numeric(gsub("F","", pop)) het <- 1/2^(gen - 1) hom <- (1 - het)/2 mm[,1] <- sample(c(1,2,3), pop.size, T, prob = c(hom,het,hom)) for(i in 2:no.mark[j]) { r <- sample(0:1, pop.size, T, prob = c(1 - 2*rfrac[i], 2*rfrac[i])) r <- (1:pop.size)[r == 1] mm[,i] <- mm[,i - 1] if(length(r)){ mc <- letters[mm[, i][r]] for(k in 1:length(mc)) mc[k] <- switch(mc[k], a = 2, b = sample(c(1,3), 1), c = 2) mm[, i][r] <- as.numeric(mc) } } cl <- "f2" } else stop("Population not allowed.") if(!is.null(error.rate)){ err <- sample(0:1, pop.size*no.mark[j], TRUE, prob = c(1 - error.rate, error.rate)) mmv <- c(mm) mmv[as.logical(err)] <- as.numeric(!as.logical(mmv[as.logical(err)])) mm <- matrix(mmv, ncol = no.mark[j]) } if(!is.null(miss.rate)){ miss <- sample(0:1, pop.size*no.mark[j], TRUE, prob = c(1 - miss.rate, miss.rate)) mmv <- c(mm) mmv[as.logical(miss)] <- NA mm <- matrix(mmv, ncol = no.mark[j]) } mm[mm == 0] <- 2 chn <- paste("c", as.character(j), sep = "") mapl$geno[[chn]] <- list() mapl$geno[[chn]]$map <- cumsum(dist) mapl$geno[[chn]]$data <- mm chrn <- rep(paste("c", as.character(j), sep = ""), no.mark[j]) markn <- paste(chrn, paste("m", 1:no.mark[j], sep = ""), sep = ".") if(!is.null(QTL[[j]])){ wh <- (1:no.mark[j])[mapl$geno[[chn]]$map %in% QTL[[j]]] no.mark[j] <- no.mark[j] - length(QTL[[j]]) markn[wh] <- paste(chrn[wh], paste("q",1:length(QTL[[j]]), sep = ""),sep = ".") markn[-wh] <- paste(chrn[-wh], paste("m",1:no.mark[j], sep = ""), sep = ".") } dimnames(mapl$geno[[chn]]$data)[[2]] <- names(mapl$geno[[chn]]$map) <- markn class(mapl$geno[[chn]]) <- "A" } mapl$pheno <- data.frame(id = factor(gn, levels = gn)) class(mapl) <- c(cl, "cross") mapl }