######################################################################
#
# argmax.geno.R
#
# copyright (c) 2001-2019, Karl W Broman
# last modified Dec, 2019
# first written Nov, 2001
#
# 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: argmax.geno
#
######################################################################
######################################################################
#
# argmax.geno: Use Viterbi algorithm to find most likely sequence of
# underlying genotypes, given observed marker data
#
######################################################################
argmax.geno <-
function(cross, step=0, off.end=0, error.prob=0.0001,
map.function=c("haldane","kosambi","c-f","morgan"),
stepwidth=c("fixed", "variable", "max"))
{
if(!inherits(cross, "cross"))
stop("cross should have class \"cross\".")
# map function
map.function <- match.arg(map.function)
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
stepwidth <- match.arg(stepwidth)
# 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)
n.chr <- nchr(cross)
n.mar <- nmar(cross)
type <- crosstype(cross)
# loop over chromosomes
for(i in 1:n.chr) {
if(n.mar[i]==1) temp.offend <- max(c(off.end,5))
else temp.offend <- off.end
chr_type <- chrtype(cross$geno[[i]])
if(chr_type=="X") xchr <- TRUE
else xchr <- FALSE
# which type of cross is this?
if(type=="f2") {
one.map <- TRUE
if(!xchr) # autosomal
cfunc <- "argmax_geno_f2"
else # X chromsome
cfunc <- "argmax_geno_bc"
}
else if(type=="bc" || type=="dh" || type=="riself" || type=="risib" || type=="haploid") {
cfunc <- "argmax_geno_bc"
one.map <- TRUE
}
else if(type == "4way") {
cfunc <- "argmax_geno_4way"
one.map <- FALSE
}
else if(type == "ri8sib" || type=="ri4sib" || type=="ri8self" || type=="ri4self" || type=="bgmagic16") {
cfunc <- paste("argmax_geno_", type, sep="")
one.map <- TRUE
if(xchr)
warning("argmax.geno not working properly for the X chromosome for 4- or 8-way RIL.")
}
else if(type == "bcsft") {
one.map <- TRUE
cfunc <- "argmax_geno_bcsft"
cross.scheme <- attr(cross, "scheme") ## c(s,t) for BC(s)F(t)
if(xchr) { ## X chr
cross.scheme[1] <- cross.scheme[1] + cross.scheme[2] - (cross.scheme[1] == 0)
cross.scheme[2] <- 0
}
}
else
stop("argmax.geno not available for cross type ", type, ".")
# genotype data
gen <- cross$geno[[i]]$data
gen[is.na(gen)] <- 0
# recombination fractions
if(one.map) {
# recombination fractions
map <- create.map(cross$geno[[i]]$map,step,temp.offend,stepwidth)
rf <- mf(diff(map))
if(type=="risib" || type=="riself")
rf <- adjust.rf.ri(rf, sub("^ri", "", type), chr_type)
rf[rf < 1e-14] <- 1e-14
# new genotype matrix with pseudomarkers filled in
newgen <- matrix(ncol=length(map),nrow=nrow(gen))
dimnames(newgen) <- list(NULL,names(map))
newgen[,colnames(gen)] <- gen
newgen[is.na(newgen)] <- 0
n.pos <- ncol(newgen)
}
else {
map <- create.map(cross$geno[[i]]$map,step,temp.offend,stepwidth)
rf <- mf(diff(map[1,]))
rf[rf < 1e-14] <- 1e-14
rf2 <- mf(diff(map[2,]))
rf2[rf2 < 1e-14] <- 1e-14
# new genotype matrix with pseudomarkers filled in
newgen <- matrix(ncol=ncol(map),nrow=nrow(gen))
dimnames(newgen) <- list(NULL,dimnames(map)[[2]])
newgen[,colnames(gen)] <- gen
newgen[is.na(newgen)] <- 0
n.pos <- ncol(newgen)
}
if(any(is.na(rf))) { # this occurs when there is only one marker
rf <- rf2 <- 0
warn <- paste("Only one marker on chr ", names(cross$geno)[i],
": argmax results tenuous.", sep="")
warning(warn)
}
# call the C function
if(one.map) {
## Hide cross scheme in genoprob to pass to routine. BY
temp <- newgen
if(type == "bcsft")
temp[1:2] <- cross.scheme
z <- .C(cfunc,
as.integer(n.ind), # number of individuals
as.integer(n.pos), # number of markers
as.integer(newgen), # genotype data
as.double(rf), # recombination fractions
as.double(error.prob),
argmax=as.integer(temp), # the output
PACKAGE="qtl")
cross$geno[[i]]$argmax <- matrix(z$argmax,ncol=n.pos)
dimnames(cross$geno[[i]]$argmax) <- list(NULL, names(map))
}
else {
z <- .C(cfunc,
as.integer(n.ind), # number of individuals
as.integer(n.pos), # number of markers
as.integer(newgen), # genotype data
as.double(rf), # recombination fractions
as.double(rf2), # recombination fractions
as.double(error.prob),
argmax=as.integer(newgen), # the output
PACKAGE="qtl")
cross$geno[[i]]$argmax <- matrix(z$argmax,ncol=n.pos)
dimnames(cross$geno[[i]]$argmax) <- list(NULL, colnames(map))
}
# attribute set to the error.prob value used, for later
# reference
attr(cross$geno[[i]]$argmax, "map") <- map
attr(cross$geno[[i]]$argmax,"error.prob") <- error.prob
attr(cross$geno[[i]]$argmax,"step") <- step
attr(cross$geno[[i]]$argmax,"off.end") <- temp.offend
attr(cross$geno[[i]]$argmax,"map.function") <- map.function
attr(cross$geno[[i]]$argmax,"stepwidth") <- stepwidth
}
# store argmax values as integers
for(i in 1:nchr(cross))
storage.mode(cross$geno[[i]]$argmax) <- "integer"
# 4- and 8-way RIL: reorganize the results
if(type=="ri4self" || type=="ri4sib" || type=="ri8self" || type=="ri8sib" || type=="bgmagic16")
cross <- reorgRIargmax(cross)
cross
}
# end of argmax.geno.R
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.