######################################################################
#
# tryallpositions.R
#
# copyright (c) 2007-2019, Karl W Broman
# last modified Dec, 2019
# first written Oct, 2007
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License,
# version 3, as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but without any warranty; without even the implied warranty of
# merchantability or fitness for a particular purpose. See the GNU
# General Public License, version 3, for more details.
#
# A copy of the GNU General Public License, version 3, is available
# at http://www.r-project.org/Licenses/GPL-3
#
# Part of the R/qtl package
# Contains: tryallpositions, markerloglik, allchrsplits
#
######################################################################
######################################################################
# tryallpositions
#
# Place a given marker in all possible positions on a selected set of
# chromosomes, keeping the positions of all other markers fixed, and
# evaluate the likelihood and estimate the chromosome length
######################################################################
tryallpositions <-
function(cross, marker, chr, error.prob=0.0001,
map.function=c("haldane","kosambi","c-f","morgan"),
m=0, p=0, maxit=4000, tol=1e-6, sex.sp=TRUE, verbose=TRUE)
{
map.function <- match.arg(map.function)
if(missing(chr)) chr <- names(cross$geno)
else chr <- matchchr(chr, names(cross$geno))
thechr <- find.markerpos(cross, marker)[1,1]
if(is.na(thechr))
stop("Marker ", marker, " not found.")
markerll <- markerloglik(cross, marker, error.prob)
allchr <- names(subset(cross, chr=chr)$geno)
results <- NULL
for(i in allchr) {
if(i == thechr) { # marker already on this chromosome
pos <- cross$geno[[i]]$map
if(is.matrix(pos)) {
pos <- pos[1,]
matrixmap <- TRUE
}
else matrixmap <- FALSE
n.mar <- ncol(cross$geno[[i]]$data)
if(n.mar == 1) { # this is the only marker
pos <- 0
length <- length.male <- 0
llik <- 0
if(verbose) cat(i, pos, llik/log(10), "\n")
interval <- "---"
} # just this marker
else if(n.mar == 2) { # just two markers
pos <- 0
themar <- colnames(cross$geno[[i]]$data)
initialloglik <- sum(sapply(themar, function(a, b, c) markerloglik(b, a, c), cross, error.prob))
nm <- est.map(cross, chr=i, error.prob=error.prob,
map.function=map.function, m=m, p=p, maxit=maxit,
tol=tol, sex.sp=sex.sp, verbose=FALSE,
omit.noninformative=FALSE)[[1]]
llik <- attr(nm, "loglik") - initialloglik
if(verbose) cat(i, pos, llik/log(10), "\n")
if(is.matrix(nm)) {
length <- diff(range(nm[1,]))
length.male <- diff(range(nm[2,]))
}
else length <- diff(range(nm))
interval <- "---"
} # just two markers
else { # >2 markers
temp <- drop.markers(cross, marker)
pos <- temp$geno[[i]]$map
if(is.matrix(pos)) {
pos <- pos[1,]
matrixmap <- TRUE
}
else matrixmap <- FALSE
pos <- c(min(pos)-10, pos, max(pos)+10)
pos <- (pos[-1] + pos[-length(pos)])/2
themarkers <- colnames(temp$geno[[i]]$data)
int2 <- c("*", match(themarkers, colnames(cross$geno[[i]]$data)), "*")
themarkers <- c("pter", themarkers, "qter")
interval <- paste(themarkers[-length(themarkers)], themarkers[-1],
sep="-")
int2 <- paste("(", int2[-length(int2)], "-", int2[-1], ")", sep="")
interval <- paste(interval, int2)
initialmap <- est.map(temp, chr=i,
error.prob=error.prob,
map.function=map.function, m=m, p=p, maxit=maxit,
tol=tol, sex.sp=sex.sp, verbose=FALSE,
omit.noninformative=FALSE)[[1]]
initialloglik <- attr(initialmap, "loglik") + markerll
llik <- length <- length.male <- rep(NA, length(pos))
for(j in seq(along=pos)) {
temp <- movemarker(cross, marker, i, pos[j])
nm <- est.map(temp, chr=i, error.prob=error.prob,
map.function=map.function, m=m, p=p, maxit=maxit,
tol=tol, sex.sp=sex.sp, verbose=FALSE,
omit.noninformative=FALSE)[[1]]
llik[j] <- attr(nm, "loglik") - initialloglik
if(verbose) cat(i, pos[j], llik[j]/log(10), "\n")
if(is.matrix(nm)) {
length[j] <- diff(range(nm[1,]))
length.male[j] <- diff(range(nm[2,]))
}
else length[j] <- diff(range(nm))
}
} # >2 markers
}
else { # marker not on this chromosome
pos <- cross$geno[[i]]$map
if(is.matrix(pos)) {
pos <- pos[1,]
matrixmap <- TRUE
}
else matrixmap <- FALSE
n.mar <- ncol(cross$geno[[i]]$data)
if(n.mar > 1) {
pos <- c(min(pos)-10, pos, max(pos)+10)
pos <- (pos[-1] + pos[-length(pos)])/2
themarkers <- colnames(cross$geno[[i]]$data)
int2 <- c("*", 1:length(themarkers), "*")
themarkers <- c("pter", themarkers, "qter")
interval <- paste(themarkers[-length(themarkers)], themarkers[-1],
sep="-")
int2 <- paste("(", int2[-length(int2)], "-", int2[-1], ")", sep="")
interval <- paste(interval, int2)
initialmap <- est.map(cross, chr=i, error.prob=error.prob,
map.function=map.function, m=m, p=p, maxit=maxit,
tol=tol, sex.sp=sex.sp, verbose=FALSE,
omit.noninformative=FALSE)[[1]]
initialloglik <- attr(initialmap, "loglik") + markerll
llik <- length <- length.male <- rep(NA, length(pos))
for(j in seq(along=pos)) {
temp <- movemarker(cross, marker, i, pos[j])
nm <- est.map(temp, chr=i, error.prob=error.prob,
map.function=map.function, m=m, p=p, maxit=maxit,
tol=tol, sex.sp=sex.sp, verbose=FALSE,
omit.noninformative=FALSE)[[1]]
llik[j] <- attr(nm, "loglik") - initialloglik
if(verbose) cat(i, pos[j], llik[j]/log(10), "\n")
if(is.matrix(nm)) {
length[j] <- diff(range(nm[1,]))
length.male[j] <- diff(range(nm[2,]))
}
else length[j] <- diff(range(nm))
}
} # >1 marker on chromosome
else {
initialloglik <- markerloglik(cross, colnames(cross$geno[[i]]$data), error.prob) + markerll
pos <- pos+10
interval <- "---"
temp <- movemarker(cross, marker, i, pos)
nm <- est.map(temp, chr=i, error.prob=error.prob,
map.function=map.function, m=m, p=p, maxit=maxit,
tol=tol, sex.sp=sex.sp, verbose=FALSE,
omit.noninformative=FALSE)[[1]]
llik <- attr(nm, "loglik") - initialloglik
if(verbose) cat(i, pos, llik/log(10), "\n")
if(is.matrix(nm)) {
length <- diff(range(nm[1,]))
length.male <- diff(range(nm[2,]))
}
else length <- diff(range(nm))
} # one marker on chromosome
} # marker not on this chr
if(matrixmap && sex.sp)
tempres <- data.frame(chr=rep(i, length(pos)),
pos=pos,
lod=llik/log(10),
length.female=length,
length.male=length.male,
interval=interval,
stringsAsFactors=FALSE)
else
tempres <- data.frame(chr=rep(i, length(pos)),
pos=pos,
lod=llik/log(10),
length=length,
interval=interval,
stringsAsFactors=FALSE)
results <- rbind(results, tempres)
} # loop over chromosomes
rownames(results) <- results$interval
results <- results[,-ncol(results)]
class(results) <- c("scanone", "data.frame")
results
}
######################################################################
#
# markerloglik: Calculate log likelihood for a given marker
#
######################################################################
markerloglik <-
function(cross, marker, error.prob=0.0001)
{
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
type <- crosstype(cross)
if(length(marker) > 1) {
ll <- sapply(marker, function(a,b,d) markerloglik(b, a, d), cross, error.prob)
names(ll) <- marker
return(ll)
}
# don't let error.prob be exactly zero (or >1)
if(error.prob < 1e-50) error.prob <- 1e-50
if(error.prob > 1) {
error.prob <- 1-1e-50
warning("error.prob shouldn't be > 1!")
}
n.ind <- nind(cross)
thechr <- find.markerpos(cross, marker)[1,1]
if(is.na(thechr))
stop("Marker ", marker, " not found.")
chr_type <- chrtype(cross$geno[[thechr]])
g <- pull.geno(cross, chr=thechr)
m <- match(marker, colnames(g))
if(is.na(m))
stop("Marker ", marker, " not found.")
g <- g[,m]
g[is.na(g)] <- 0
# which type of cross is this?
if(type == "f2") {
if(chr_type == "A") # autosomal
cfunc <- "marker_loglik_f2"
else # X chromsome
cfunc <- "marker_loglik_bc"
}
else if(type == "bc" || type=="riself" || type=="risib" || type=="dh" || type=="haploid") {
cfunc <- "marker_loglik_bc"
}
else if(type == "4way") {
cfunc <- "marker_loglik_4way"
}
else if(type=="ri4sib" || type=="ri4self" || type=="ri8sib" || type=="ri8self" || type=="bgmagic16") {
cfunc <- paste("marker_loglik_", type, sep="")
if(chr_type=="X")
warning("markerloglik not working properly for the X chromosome for 4- or 8-way RIL.")
}
else if(type == "bcsft") {
cfunc <- "marker_loglik_bcsft"
cross.scheme <- attr(cross, "scheme") ## c(s,t) for BC(s)F(t)
if(chr_type != "A") { ## X chromosome
cross.scheme[1] <- cross.scheme[1] + cross.scheme[2] - (cross.scheme[1] == 0)
cross.scheme[2] <- 0
}
}
else
stop("markerloglik not available for cross type ", type, ".")
## Hide cross scheme in genoprob to pass to routine. BY
temp <- 0
if(type == "bcsft")
temp[1] <- cross.scheme[1] * 1000 + cross.scheme[2]
# call the C function
z <- .C(cfunc,
as.integer(n.ind), # number of individuals
as.integer(g), # genotype data
as.double(error.prob),
loglik=as.double(temp), # log likelihood
PACKAGE="qtl")
z$loglik
}
######################################################################
# allchrsplits
#
# get LOD scores for each possible split of each chromosome into
# two pieces
#
######################################################################
allchrsplits <-
function(cross, chr, error.prob=0.0001,
map.function=c("haldane","kosambi","c-f","morgan"),
m=0, p=0, maxit=4000, tol=1e-6, sex.sp=TRUE, verbose=TRUE)
{
map.function <- match.arg(map.function)
if(!missing(chr)) cross <- subset(cross, chr=chr)
biggap <- imf.h(0.5 - 1e-14)
n.mar <- nmar(cross)
chrnam <- names(cross$geno)
result <- NULL
for(i in seq(along=cross$geno)) {
if(n.mar[i] == 1) {
# temp <- data.frame(chr=chrnam[i], pos=pos, lod=NA, gap=0)
# rownames(temp) <- themarkers
# result <- cbind(result, temp)
next
}
thischr <- subset(cross, chr=chrnam[i])
if(verbose) cat("Chr ", chrnam[i], " (", n.mar[i], " markers)\n", sep="")
pos <- cross$geno[[i]]$map
themarkers <- colnames(cross$geno[[i]]$data)
if(is.matrix(pos)) {
pos <- pos[1,]
matrixmap <- TRUE
}
else matrixmap <- FALSE
gap <- (pos[-1] - pos[-length(pos)])
pos <- (pos[-1] + pos[-length(pos)])/2
int2 <- match(themarkers, colnames(cross$geno[[i]]$data))
interval <- paste(themarkers[-length(themarkers)], themarkers[-1],
sep="-")
int2 <- paste("(", int2[-length(int2)], "-", int2[-1], ")", sep="")
interval <- paste(interval, int2)
initialmap <- est.map(thischr, error.prob=error.prob,
map.function=map.function, m=m, p=p, maxit=maxit, tol=tol, sex.sp=sex.sp)
thischr$geno[[1]]$map <- initialmap[[1]]
initialloglik <- attr(initialmap[[1]], "loglik")
if(n.mar[i] == 2) { # 2 markers
mmll <- markerloglik(thischr, markernames(thischr), error.prob=error.prob)
temp <- data.frame(chr=chrnam[i], pos=pos,
lod=(initialloglik - sum(mmll))/log(10), gap=gap, stringsAsFactors=TRUE)
rownames(temp) <- interval
}
else { # >2 markers
mn <- markernames(thischr)
mmll <- markerloglik(thischr, mn[c(1,length(mn))], error.prob=error.prob)
lod <- rep(NA, length(mn)-1)
if(verbose) cat(" interval 1\n")
# first interval
lod[1] <- initialloglik - mmll[1] -
attr(est.map(drop.markers(thischr, mn[1]), error.prob=error.prob,
map.function=map.function, m=m, p=p, maxit=maxit, tol=tol, sex.sp=sex.sp)[[1]], "loglik")
if(n.mar[i] > 3) {
for(j in 2:(n.mar[i]-2)) {
if(verbose) cat(" interval", j, "\n")
temp1 <- est.map(pull.markers(thischr, mn[1:j]), error.prob=error.prob,
map.function=map.function, m=m, p=p, maxit=maxit, tol=tol, sex.sp=sex.sp)[[1]]
temp2 <- est.map(drop.markers(thischr, mn[1:j]), error.prob=error.prob,
map.function=map.function, m=m, p=p, maxit=maxit, tol=tol, sex.sp=sex.sp)[[1]]
# lod[j] <- initialloglik - attr(temp1, "loglik") - attr(temp2, "loglik")
if(any(is.na(temp1)) || any(is.na(temp2)))
stop("Missing values in estimated map on chr ", chrnam[i], " with split at interval ", j, "\n")
# the likelihoods aren't adding properly, so I'll use the following kluge:
temp3 <- thischr
if(is.matrix(temp1))
temp3$geno[[1]]$map <- cbind(temp1, biggap+temp2)
else
temp3$geno[[1]]$map <- c(temp1, biggap+temp2)
temp3 <- est.map(temp3, error.prob=error.prob, map.function=map.function, m=m, p=p, maxit=0, tol=tol,
sex.sp=sex.sp)[[1]]
lod[j] <- initialloglik - attr(temp3, "loglik")
}
}
if(verbose) cat(" interval", n.mar[i]-1, "\n")
# last interval
lod[length(mn)-1] <- initialloglik - mmll[2] -
attr(est.map(drop.markers(thischr, mn[length(mn)]), error.prob=error.prob,
map.function=map.function, m=m, p=p, maxit=maxit, tol=tol, sex.sp=sex.sp)[[1]], "loglik")
temp <- data.frame(chr=rep(chrnam[i], length(interval)), pos=pos, lod=lod/log(10), gap=gap, stringsAsFactors=TRUE)
rownames(temp) <- interval
}
result <- rbind(result, temp)
}
class(result) <- c("scanone", "data.frame")
result
}
# end of tryallpositions.R
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.