R/calc.genoprob.R
In kbroman/qtl: Tools for Analyzing QTL Experiments
Defines functions
calc.genoprob.special
calc.genoprob
Documented in
calc.genoprob
######################################################################
#
# calc.genoprob.R
#
# copyright (c) 2001-2019, Karl W Broman
# last modified Dec, 2019
# first written Feb, 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: calc.genoprob
#
######################################################################
######################################################################
#
# calc.genoprob: calculate genotype probabilities conditional on
# observed marker genotypes
#
######################################################################
calc.genoprob <-
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("Input 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)
# calculate genotype probabilities one chromosome at a time
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 <- "calc_genoprob_f2"
n.gen <- 3
gen.names <- getgenonames("f2", "A", cross.attr=attributes(cross))
}
else { # X chromsome
cfunc <- "calc_genoprob_bc"
n.gen <- 2
gen.names <- c("g1","g2")
}
}
else if(type == "bc") {
cfunc <- "calc_genoprob_bc"
n.gen <- 2
if(!xchr)
gen.names <- getgenonames("bc", "A", cross.attr=attributes(cross))
else gen.names <- c("g1","g2")
one.map <- TRUE
}
else if(type == "riself" || type=="risib" || type=="dh" || type=="haploid") {
cfunc <- "calc_genoprob_bc"
n.gen <- 2
gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
one.map <- TRUE
}
else if(type == "4way") {
cfunc <- "calc_genoprob_4way"
n.gen <- 4
one.map <- FALSE
gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
}
else if(type=="ri8sib" || type=="ri4sib" || type=="ri8self" || type=="ri4self" || type=="bgmagic16") {
cfunc <- paste("calc_genoprob_", type, sep="")
if(type=="bgmagic16") n.gen <- 16
else n.gen <- as.numeric(substr(type, 3, 3))
one.map <- TRUE
gen.names <- LETTERS[1:n.gen]
if(xchr)
warning("calc.genoprob not working properly for the X chromosome for 4- or 8-way RIL.")
}
else if(type == "bcsft") {
one.map <- TRUE
cfunc <- "calc_genoprob_bcsft"
cross.scheme <- attr(cross, "scheme") ## c(s,t) for BC(s)F(t)
if(!xchr) { # autosomal
gen.names <- getgenonames("bcsft", "A", cross.attr=attributes(cross))
n.gen <- 2 + (cross.scheme[2] > 0)
}
else { ## X chr
cross.scheme[1] <- cross.scheme[1] + cross.scheme[2] - (cross.scheme[1] == 0)
cross.scheme[2] <- 0
gen.names <- c("g1","g2")
n.gen <- 2
}
}
else if(type == "ri8selfIRIP1") {
one.map <- TRUE
cfunc <- "calc_genoprob_ri8selfIRIP1"
n.gen <- 8
gen.names <- LETTERS[1:n.gen]
}
else
stop("calc.genoprob 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)
marnames <- names(map)
}
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)
marnames <- colnames(map)
}
## Cross scheme being added for Ft and BCs.
## cross_scheme = c(BC = s, F = t).
## BC has cross_scheme = c(1,0)
## F2 has cross_scheme = c(0,2)
## BCsFt has cross_scheme = c(s,t)
## Other designs such as Ft with test cross don't quite fit in (yet).
## *** Need to change all the C routines!!
# call the C function
if(one.map) {
## Hide cross scheme in genoprob to pass to routine. BY
temp <- as.double(rep(0,n.gen*n.ind*n.pos))
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), #
genoprob=as.double(temp),
PACKAGE="qtl")
}
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), #
genoprob=as.double(rep(0,n.gen*n.ind*n.pos)),
PACKAGE="qtl")
}
# re-arrange marginal probabilites
cross$geno[[i]]$prob <- array(z$genoprob,dim=c(n.ind,n.pos,n.gen))
dimnames(cross$geno[[i]]$prob) <- list(NULL, marnames, gen.names)
# attribute set to the error.prob value used, for later
# reference, especially by calc.errorlod()
attr(cross$geno[[i]]$prob, "map") <- map
attr(cross$geno[[i]]$prob,"error.prob") <- error.prob
attr(cross$geno[[i]]$prob,"step") <- step
attr(cross$geno[[i]]$prob,"off.end") <- temp.offend
attr(cross$geno[[i]]$prob,"map.function") <- map.function
attr(cross$geno[[i]]$prob,"stepwidth") <- stepwidth
} # end loop over chromosomes
# 4- and 8-way RIL: reorganize the results
if(type=="ri4self" || type=="ri4sib" || type=="ri8self" || type=="ri8sib" || type=="bgmagic16")
cross <- reorgRIgenoprob(cross)
cross
}
######################################################################
#
# calc.genoprob.special:
# special version used by calc.errorlod
# for each individual and marker, calculate probabilities allowing
# that genotype to be in error but assuming that all other genotypes
# are correct
#
######################################################################
calc.genoprob.special <-
function(cross, error.prob=0.0001,
map.function=c("haldane","kosambi","c-f","morgan"))
{
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
step <- 0
off.end <- 0
stepwidth <- "fixed"
# 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
# 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)
# calculate genotype probabilities one chromosome at a time
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 <- "calc_genoprob_special_f2"
n.gen <- 3
gen.names <- getgenonames("f2", "A", cross.attr=attributes(cross))
}
else { # X chromsome
cfunc <- "calc_genoprob_special_bc"
n.gen <- 2
gen.names <- c("g1","g2")
}
}
else if(type == "bc") {
cfunc <- "calc_genoprob_special_bc"
n.gen <- 2
if(!xchr)
gen.names <- getgenonames("bc", "A", cross.attr=attributes(cross))
else gen.names <- c("g1","g2")
one.map <- TRUE
}
else if(type == "riself" || type=="risib" || type=="dh" || type=="haploid") {
cfunc <- "calc_genoprob_special_bc"
n.gen <- 2
gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
one.map <- TRUE
}
else if(type == "4way") {
cfunc <- "calc_genoprob_special_4way"
n.gen <- 4
one.map <- FALSE
gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
}
else if(type=="ri8sib" || type=="ri4sib" || type=="ri8self" || type=="ri4self" || type=="bgmagic16") {
cfunc <- paste("calc_genoprob_special_", type, sep="")
if(type=="bgmagic16") n.gen <- 16
else n.gen <- as.numeric(substr(type, 3, 3))
one.map <- TRUE
gen.names <- LETTERS[1:n.gen]
if(xchr)
warning("calc.genoprob.special not working properly for the X chromosome for 4- or 8-way RIL.")
}
else if(type == "bcsft") {
one.map <- TRUE
cfunc <- "calc_genoprob_bcsft"
cross.scheme <- attr(cross, "scheme") ## c(s,t) for BC(s)F(t)
if(!xchr) { # autosomal
if(cross.scheme[2] == 0) {
gen.names <- getgenonames("bc", "A", cross.attr=attributes(cross))
n.gen <- 2
}
else {
gen.names <- getgenonames("f2", "A", cross.attr=attributes(cross))
n.gen <- 3
}
}
else { ## X chr
cross.scheme[1] <- cross.scheme[1] + cross.scheme[2] - (cross.scheme[1] == 0)
cross.scheme[2] <- 0
gen.names <- c("g1","g2")
n.gen <- 2
}
}
else
stop("calc.genoprob.special 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)
marnames <- names(map)
}
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)
marnames <- colnames(map)
}
# call the C function
if(one.map) {
## Hide cross scheme in genoprob to pass to routine. BY
temp <- as.double(rep(0,n.gen*n.ind*n.pos))
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), #
genoprob=as.double(temp),
PACKAGE="qtl")
}
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), #
genoprob=as.double(rep(0,n.gen*n.ind*n.pos)),
PACKAGE="qtl")
}
# re-arrange marginal probabilites
cross$geno[[i]]$prob <- array(z$genoprob,dim=c(n.ind,n.pos,n.gen))
dimnames(cross$geno[[i]]$prob) <- list(NULL, marnames, gen.names)
# attribute set to the error.prob value used, for later
# reference, especially by calc.errorlod()
attr(cross$geno[[i]]$prob, "map") <- map
attr(cross$geno[[i]]$prob,"error.prob") <- error.prob
attr(cross$geno[[i]]$prob,"step") <- step
attr(cross$geno[[i]]$prob,"off.end") <- temp.offend
attr(cross$geno[[i]]$prob,"map.function") <- map.function
attr(cross$geno[[i]]$prob,"stepwidth") <- stepwidth
} # end loop over chromosomes
cross
}
# end of calc.genoprob.R
kbroman/qtl documentation built on Jan. 13, 2024, 10:14 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions calc.genoprob.special calc.genoprob
Documented in calc.genoprob
######################################################################
#
# calc.genoprob.R
#
# copyright (c) 2001-2019, Karl W Broman
# last modified Dec, 2019
# first written Feb, 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: calc.genoprob
#
######################################################################
######################################################################
#
# calc.genoprob: calculate genotype probabilities conditional on
# observed marker genotypes
#
######################################################################
calc.genoprob <-
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("Input 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)
# calculate genotype probabilities one chromosome at a time
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 <- "calc_genoprob_f2"
n.gen <- 3
gen.names <- getgenonames("f2", "A", cross.attr=attributes(cross))
}
else { # X chromsome
cfunc <- "calc_genoprob_bc"
n.gen <- 2
gen.names <- c("g1","g2")
}
}
else if(type == "bc") {
cfunc <- "calc_genoprob_bc"
n.gen <- 2
if(!xchr)
gen.names <- getgenonames("bc", "A", cross.attr=attributes(cross))
else gen.names <- c("g1","g2")
one.map <- TRUE
}
else if(type == "riself" || type=="risib" || type=="dh" || type=="haploid") {
cfunc <- "calc_genoprob_bc"
n.gen <- 2
gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
one.map <- TRUE
}
else if(type == "4way") {
cfunc <- "calc_genoprob_4way"
n.gen <- 4
one.map <- FALSE
gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
}
else if(type=="ri8sib" || type=="ri4sib" || type=="ri8self" || type=="ri4self" || type=="bgmagic16") {
cfunc <- paste("calc_genoprob_", type, sep="")
if(type=="bgmagic16") n.gen <- 16
else n.gen <- as.numeric(substr(type, 3, 3))
one.map <- TRUE
gen.names <- LETTERS[1:n.gen]
if(xchr)
warning("calc.genoprob not working properly for the X chromosome for 4- or 8-way RIL.")
}
else if(type == "bcsft") {
one.map <- TRUE
cfunc <- "calc_genoprob_bcsft"
cross.scheme <- attr(cross, "scheme") ## c(s,t) for BC(s)F(t)
if(!xchr) { # autosomal
gen.names <- getgenonames("bcsft", "A", cross.attr=attributes(cross))
n.gen <- 2 + (cross.scheme[2] > 0)
}
else { ## X chr
cross.scheme[1] <- cross.scheme[1] + cross.scheme[2] - (cross.scheme[1] == 0)
cross.scheme[2] <- 0
gen.names <- c("g1","g2")
n.gen <- 2
}
}
else if(type == "ri8selfIRIP1") {
one.map <- TRUE
cfunc <- "calc_genoprob_ri8selfIRIP1"
n.gen <- 8
gen.names <- LETTERS[1:n.gen]
}
else
stop("calc.genoprob 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)
marnames <- names(map)
}
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)
marnames <- colnames(map)
}
## Cross scheme being added for Ft and BCs.
## cross_scheme = c(BC = s, F = t).
## BC has cross_scheme = c(1,0)
## F2 has cross_scheme = c(0,2)
## BCsFt has cross_scheme = c(s,t)
## Other designs such as Ft with test cross don't quite fit in (yet).
## *** Need to change all the C routines!!
# call the C function
if(one.map) {
## Hide cross scheme in genoprob to pass to routine. BY
temp <- as.double(rep(0,n.gen*n.ind*n.pos))
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), #
genoprob=as.double(temp),
PACKAGE="qtl")
}
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), #
genoprob=as.double(rep(0,n.gen*n.ind*n.pos)),
PACKAGE="qtl")
}
# re-arrange marginal probabilites
cross$geno[[i]]$prob <- array(z$genoprob,dim=c(n.ind,n.pos,n.gen))
dimnames(cross$geno[[i]]$prob) <- list(NULL, marnames, gen.names)
# attribute set to the error.prob value used, for later
# reference, especially by calc.errorlod()
attr(cross$geno[[i]]$prob, "map") <- map
attr(cross$geno[[i]]$prob,"error.prob") <- error.prob
attr(cross$geno[[i]]$prob,"step") <- step
attr(cross$geno[[i]]$prob,"off.end") <- temp.offend
attr(cross$geno[[i]]$prob,"map.function") <- map.function
attr(cross$geno[[i]]$prob,"stepwidth") <- stepwidth
} # end loop over chromosomes
# 4- and 8-way RIL: reorganize the results
if(type=="ri4self" || type=="ri4sib" || type=="ri8self" || type=="ri8sib" || type=="bgmagic16")
cross <- reorgRIgenoprob(cross)
cross
}
######################################################################
#
# calc.genoprob.special:
# special version used by calc.errorlod
# for each individual and marker, calculate probabilities allowing
# that genotype to be in error but assuming that all other genotypes
# are correct
#
######################################################################
calc.genoprob.special <-
function(cross, error.prob=0.0001,
map.function=c("haldane","kosambi","c-f","morgan"))
{
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
step <- 0
off.end <- 0
stepwidth <- "fixed"
# 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
# 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)
# calculate genotype probabilities one chromosome at a time
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 <- "calc_genoprob_special_f2"
n.gen <- 3
gen.names <- getgenonames("f2", "A", cross.attr=attributes(cross))
}
else { # X chromsome
cfunc <- "calc_genoprob_special_bc"
n.gen <- 2
gen.names <- c("g1","g2")
}
}
else if(type == "bc") {
cfunc <- "calc_genoprob_special_bc"
n.gen <- 2
if(!xchr)
gen.names <- getgenonames("bc", "A", cross.attr=attributes(cross))
else gen.names <- c("g1","g2")
one.map <- TRUE
}
else if(type == "riself" || type=="risib" || type=="dh" || type=="haploid") {
cfunc <- "calc_genoprob_special_bc"
n.gen <- 2
gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
one.map <- TRUE
}
else if(type == "4way") {
cfunc <- "calc_genoprob_special_4way"
n.gen <- 4
one.map <- FALSE
gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
}
else if(type=="ri8sib" || type=="ri4sib" || type=="ri8self" || type=="ri4self" || type=="bgmagic16") {
cfunc <- paste("calc_genoprob_special_", type, sep="")
if(type=="bgmagic16") n.gen <- 16
else n.gen <- as.numeric(substr(type, 3, 3))
one.map <- TRUE
gen.names <- LETTERS[1:n.gen]
if(xchr)
warning("calc.genoprob.special not working properly for the X chromosome for 4- or 8-way RIL.")
}
else if(type == "bcsft") {
one.map <- TRUE
cfunc <- "calc_genoprob_bcsft"
cross.scheme <- attr(cross, "scheme") ## c(s,t) for BC(s)F(t)
if(!xchr) { # autosomal
if(cross.scheme[2] == 0) {
gen.names <- getgenonames("bc", "A", cross.attr=attributes(cross))
n.gen <- 2
}
else {
gen.names <- getgenonames("f2", "A", cross.attr=attributes(cross))
n.gen <- 3
}
}
else { ## X chr
cross.scheme[1] <- cross.scheme[1] + cross.scheme[2] - (cross.scheme[1] == 0)
cross.scheme[2] <- 0
gen.names <- c("g1","g2")
n.gen <- 2
}
}
else
stop("calc.genoprob.special 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)
marnames <- names(map)
}
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)
marnames <- colnames(map)
}
# call the C function
if(one.map) {
## Hide cross scheme in genoprob to pass to routine. BY
temp <- as.double(rep(0,n.gen*n.ind*n.pos))
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), #
genoprob=as.double(temp),
PACKAGE="qtl")
}
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), #
genoprob=as.double(rep(0,n.gen*n.ind*n.pos)),
PACKAGE="qtl")
}
# re-arrange marginal probabilites
cross$geno[[i]]$prob <- array(z$genoprob,dim=c(n.ind,n.pos,n.gen))
dimnames(cross$geno[[i]]$prob) <- list(NULL, marnames, gen.names)
# attribute set to the error.prob value used, for later
# reference, especially by calc.errorlod()
attr(cross$geno[[i]]$prob, "map") <- map
attr(cross$geno[[i]]$prob,"error.prob") <- error.prob
attr(cross$geno[[i]]$prob,"step") <- step
attr(cross$geno[[i]]$prob,"off.end") <- temp.offend
attr(cross$geno[[i]]$prob,"map.function") <- map.function
attr(cross$geno[[i]]$prob,"stepwidth") <- stepwidth
} # end loop over chromosomes
cross
}
# end of calc.genoprob.R
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.