Nothing
#####################################################################
#
# xchr.R
#
# copyright (c) 2004-2018, Karl W Broman
# last modified Mar, 2018
# first written Apr, 2004
#
# 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: Utilities for dealing with the X chromosome.
# getsex, getgenonames, reviseXdata, scanoneXnull
# revisecovar, dropXcol
# [See also fixXgeno.bc & fixXgeno.f2 in read.cross.R]
#
######################################################################
# get sex and pgm columns from phenotype data
getsex <-
function(cross)
{
type <- crosstype(cross)
if(type != "bc" && type != "f2" && type != "4way") return(list(sex=NULL, pgm=NULL))
phe.names <- names(cross$pheno)
sex.column <- grep("^[Ss][Ee][Xx]$", phe.names)
pgm.column <- grep("^[Pp][Gg][Mm]$", phe.names)
if(length(sex.column)==0) { # no sex included
sex <- NULL
}
else {
if(length(sex.column)>1)
warning("'sex' included multiple times. Using the first one.")
temp <- cross$pheno[,sex.column[1]]
if(is.numeric(temp)) {
if(any(!is.na(temp) & temp != 0 & temp != 1)) {
warning("Sex column should be coded as 0=female 1=male; sex ignored.")
sex <- NULL
}
else sex <- temp
}
else {
if(!is.factor(temp)) temp <- as.factor(temp)
if(length(levels(temp)) == 1) {
if(levels(temp) == "F" || levels(temp)=="f" ||
toupper(levels(temp)) == "FEMALE") sex <- rep(0,nind(cross))
else if(levels(temp) == "M" || levels(temp)=="m" ||
toupper(levels(temp)) == "MALE") sex <- rep(1,nind(cross))
else
warning("Sex column should be coded as 0=female 1=male; sex ignored.")
}
else if(length(levels(temp)) > 2) {
warning("Sex column should be coded as a two-level factor; sex ignored.")
sex <- NULL
}
else { # is a factor with two levels
lev <- levels(temp)
if(length(grep("^[Ff]",lev))>0 &&
length(males <- grep("^[Mm]",lev))>0) {
temp <- as.character(temp)
sex <- rep(0,length(temp))
sex[is.na(temp)] <- NA
sex[!is.na(temp) & temp==lev[males]] <- 1
}
else
warning("Don't understand levels in sex column; sex ignored.")
}
}
}
if(length(pgm.column)==0 || type=="4way") { # no pgm included
pgm <- NULL
}
else {
if(length(pgm.column)>1)
warning("'pgm' included multiple times. Using the first one.")
temp <- cross$pheno[,pgm.column[1]]
if(!is.numeric(temp))
temp <- as.numeric(temp)-1
if(any(!is.na(temp) & temp != 0 & temp != 1)) {
warning("pgm column should be coded as 0/1; pgm ignored.")
pgm <- NULL
}
else pgm <- temp
}
if(!is.null(sex) && any(is.na(sex))) {
if(all(sex[!is.na(sex)]==1)) {
warning(sum(is.na(sex)), " individuals with missing sex; assuming they're male like the others")
sex[is.na(sex)] <- 1
}
else if(all(sex[!is.na(sex)]==0)) {
warning(sum(is.na(sex)), " individuals with missing sex; assuming they're female like the others")
sex[is.na(sex)] <- 0
}
else {
warning(sum(is.na(sex)), " individuals with missing sex; assuming they're female")
sex[is.na(sex)] <- 0
}
}
if(!is.null(pgm) && any(is.na(pgm))) {
if(all(pgm[!is.na(pgm)]==1)) {
warning(sum(is.na(pgm)), " individuals with missing pgm; assuming pgm==1 like the others")
pgm[is.na(pgm)] <- 1
}
else if(all(pgm[!is.na(pgm)]==0)) {
warning(sum(is.na(pgm)), " individuals with missing pgm; assuming pgm==0 like the others")
pgm[is.na(pgm)] <- 0
}
else {
warning(sum(is.na(pgm)), " individuals with missing pgm; assuming pgm==0")
pgm[is.na(pgm)] <- 0
}
}
list(sex=sex,pgm=pgm)
}
# get names of genotypes
# used in discan, effectplot, plotPXG, scanone, scantwo, vbscan, reviseXdata
# cross.attr gives the cross attributes
getgenonames <-
function(type=c("f2","bc","riself","risib","4way","dh","haploid","special","bcsft"),
chrtype=c("A","X"), expandX=c("simple","standard","full"),
sexpgm, cross.attr)
{
type <- match.arg(type)
chrtype <- match.arg(chrtype)
expandX <- match.arg(expandX)
## Treat bcsft as bc if no intercross generations; otherwise as f2.
if(type == "bcsft") {
if(cross.attr$scheme[2] == 0)
type <- "bc"
else
type <- "f2"
}
if(chrtype=="X") {
sex <- sexpgm$sex
pgm <- sexpgm$pgm
}
if(type=="special") return(cross.attr$genotypes)
if(missing(cross.attr) || !("alleles" %in% names(cross.attr))) {
if(type == "4way") alleles <- LETTERS[1:4]
else alleles <- LETTERS[1:2]
}
else {
alleles <- cross.attr$alleles
}
if(type=="4way") { # ensure that we have enough allele codes
if(length(alleles) < 4) alleles <- LETTERS[1:4]
} else {
if(length(alleles) < 2) alleles <- LETTERS[1:2]
}
tempgn <- c(paste(rep(alleles[1],2),collapse=""),
paste(alleles,collapse=""),
paste(rep(alleles[2],2),collapse=""),
paste(alleles[1],"Y",sep=""),
paste(alleles[2],"Y",sep=""))
# get rid of missing sex and pgm values, if there are any
if(chrtype=="X") {
if(length(sex)>0) sex <- sex[!is.na(sex)]
if(length(pgm)>0) pgm <- pgm[!is.na(pgm)]
}
if(type=="riself" || type=="risib" || type=="dh")
gen.names <- tempgn[c(1,3)]
else if(type=="haploid")
gen.names <- alleles
else if(type == "4way") {
if(chrtype=="A")
gen.names <- c(paste(alleles[1],alleles[3],sep=""),
paste(alleles[2],alleles[3],sep=""),
paste(alleles[1],alleles[4],sep=""),
paste(alleles[2],alleles[4],sep=""))
else
gen.names <- c(paste(alleles[1],alleles[3],sep=""),
paste(alleles[2],alleles[3],sep=""),
paste(alleles[1],"Y",sep=""),
paste(alleles[2],"Y",sep=""))
}
else if(type == "bc") {
if(chrtype=="A") # autosome
gen.names <- tempgn[1:2] # AA AB
else { # X chromosome
# simple standard full
# -both sexes A-/AB/BY AA/AB/AY/BY same as std
# -all females AA/AB same same
# -all males AY/BY same same
if(length(sex)==0 || all(sex==0)) # all females
gen.names <- tempgn[1:2] # AA AB
else if(all(sex==1)) # all males
gen.names <- tempgn[4:5] # AY BY
else { # some of each
if(expandX == "simple")
gen.names <- c(paste(alleles[1], "-", sep=""),
tempgn[c(2,5)]) # A-, AB, BY
else gen.names <- tempgn[c(1,2,4,5)] # AA,AB,AY,BY
}
}
}
else { # intercross
if(chrtype == "A") # autosomal
gen.names <- tempgn[1:3]
else { # X chromsome
# both crosses simple standard full
# -both sexes A-/AB/B- AA/AB/BB/AY/BY AA/AB1/AB2/BB/AY/BY
# -all females AA/AB/BB same as simple AA/AB1/AB2/BB
# -all males AY/BY same same
# forw cross
# -both sexes A-/AB/BY AA/AB/AY/BY same as std
# -all females AA/AB same same
# -all males AY/BY same same
# backw cross
# -both sexes B-/AB/AY BB/AB/AY/BY same as std
# -all females BB/AB same same
# -all males AY/BY same same
if(length(sex)==0 || all(sex==0)) { # all females
if(length(pgm)==0 || all(pgm==0)) # all forw dir
gen.names <- tempgn[1:2] # AA AB
else if(all(pgm==1)) # all backw dir
gen.names <- tempgn[3:2] # BB AB
else { # some of each direction
if(expandX=="full")
gen.names <- c(tempgn[1],
paste(tempgn[2],c("f","r"), sep=""),
tempgn[3])
else gen.names <- tempgn[1:3]
}
}
else if(all(sex==1)) # all males
gen.names <- tempgn[4:5]
else { # some of each sex
if(length(pgm)==0 || all(pgm==0)) { # all forw
if(expandX=="simple")
gen.names <- c(paste(alleles[1],"-", sep=""),
tempgn[c(2,5)])
else gen.names <- tempgn[c(1,2,4,5)]
}
else if (all(pgm==1)) { # all backw
if(expandX=="simple")
gen.names <- c(paste(alleles[2], "-",sep=""),
tempgn[c(2,4)])
else gen.names <- tempgn[c(3,2,4,5)]
}
else { # some of each dir
if(expandX=="simple")
gen.names <- c(paste(alleles[1],"-",sep=""),
tempgn[2],
paste(alleles[2],"-",sep=""))
else if(expandX=="standard")
gen.names <- tempgn
else
gen.names <- c(tempgn[1],
paste(tempgn[2],c("f","r"),sep=""),
tempgn[3:5])
}
}
}
}
gen.names
}
# revise genotype data, probabilities or imputations for the X chromosome
reviseXdata <-
function(type=c("f2","bc","bcsft"), expandX=c("simple","standard","full"),
sexpgm, geno, prob, draws, pairprob, cross.attr, force=FALSE)
{
type <- match.arg(type)
expandX <- match.arg(expandX)
## Treat bcsft as bc if no intercross generations; otherwise as f2.
if(type == "bcsft") {
if(cross.attr$scheme[2] == 0)
type <- "bc"
else
type <- "f2"
}
sex <- sexpgm$sex
pgm <- sexpgm$pgm
notmissing <- (!missing(geno)) + (!missing(prob)) + (!missing(draws)) +
(!missing(pairprob))
if(notmissing == 0)
stop("Provide one of geno, prob, draws, pairprob.")
if(notmissing > 1)
stop("Provide just one of geno, prob, draws, pairprob.")
# get genonames
genonames <- getgenonames(type, "X", expandX, sexpgm, cross.attr)
if(type == "bc") { # backcross
if(length(sex)==0 || ((all(sex==0) || all(sex==1)) && !force)) { # all one sex
# no changes necessary
if(!missing(geno)) return(geno)
else if(!missing(prob)) {
dimnames(prob)[[3]] <- genonames
return(prob)
}
else if(!missing(draws))
return(draws)
else # pairprob
return(pairprob)
}
else { # both sexes
if(!missing(geno)) {
gmale <- geno[sex==1,]
if(expandX=="simple")
gmale[!is.na(gmale) & gmale==2] <- 3
else {
gmale[!is.na(gmale) & gmale==1] <- 3
gmale[!is.na(gmale) & gmale==2] <- 4
}
geno[sex==1,] <- gmale
return(geno)
}
else if(!missing(draws)) {
gmale <- draws[sex==1,,]
if(expandX=="simple")
gmale[gmale==2] <- 3
else {
gmale[gmale==1] <- 3
gmale[gmale==2] <- 4
}
draws[sex==1,,] <- gmale
return(draws)
}
else if(!missing(prob)) {
dimprob <- dim(prob)
dimprob[3] <- length(genonames)
newprob <- array(0,dim=dimprob)
dimnames(newprob) <- c(dimnames(prob)[1:2],list(genonames))
newprob[sex==0,,1:2] <- prob[sex==0,,1:2]
if(expandX=="simple") {
newprob[sex==1,,1] <- prob[sex==1,,1]
newprob[sex==1,,3] <- prob[sex==1,,2]
}
else {
newprob[sex==1,,3] <- prob[sex==1,,1]
newprob[sex==1,,4] <- prob[sex==1,,2]
}
return(newprob)
}
else { # pairprob
dimpairprob <- dim(pairprob)
dimpairprob[3] <- dimpairprob[4] <- length(genonames)
newpairprob <- array(0,dim=dimpairprob)
newpairprob[sex==0,,1:2,1:2] <- pairprob[sex==0,,,]
if(expandX=="simple") {
newpairprob[sex==1,,1,1] <- pairprob[sex==1,,1,1]
newpairprob[sex==1,,1,3] <- pairprob[sex==1,,1,2]
newpairprob[sex==1,,3,1] <- pairprob[sex==1,,2,1]
newpairprob[sex==1,,3,3] <- pairprob[sex==1,,2,2]
}
else {
newpairprob[sex==1,,3,3] <- pairprob[sex==1,,1,1]
newpairprob[sex==1,,3,4] <- pairprob[sex==1,,1,2]
newpairprob[sex==1,,4,3] <- pairprob[sex==1,,2,1]
newpairprob[sex==1,,4,4] <- pairprob[sex==1,,2,2]
}
return(newpairprob)
}
} # end of "both sexes" / backcross
} # end of backcross
else { # intercross
if(length(sex)==0 || all(sex==0)) { # all females
if(length(pgm)==0 || ((all(pgm==0) || all(pgm==1)) && !force)) { # one dir, females
if(!missing(geno)) return(geno)
else if(!missing(draws)) return(draws)
else if(!missing(pairprob)) return(pairprob)
else {
dimnames(prob)[[3]] <- genonames
return(prob)
}
}
else { # both dir, females
if(!missing(geno)) {
gback <- geno[pgm==1,]
if(expandX!="full") {
gback[!is.na(gback) & gback==1] <- 3
geno[pgm==1,] <- gback
}
else {
gback[!is.na(gback) & gback==1] <- 4
gback[!is.na(gback) & gback==2] <- 3
geno[pgm==1,] <- gback
}
return(geno)
}
else if(!missing(draws)) {
gback <- draws[pgm==1,,]
if(expandX!="full") {
gback[!is.na(gback) & gback==1] <- 3
}
else {
gback[!is.na(gback) & gback==1] <- 4
gback[!is.na(gback) & gback==2] <- 3
}
draws[pgm==1,,] <- gback
return(draws)
}
else if(!missing(prob)) {
dimprob <- dim(prob)
dimprob[3] <- length(genonames)
newprob <- array(0,dim=dimprob)
dimnames(newprob) <- c(dimnames(prob)[1:2],list(genonames))
newprob[pgm==0,,1:2] <- prob[pgm==0,,1:2]
if(expandX!="full") { # simple/standard
newprob[pgm==1,,3] <- prob[pgm==1,,1]
newprob[pgm==1,,2] <- prob[pgm==1,,2]
}
else {
newprob[pgm==1,,4] <- prob[pgm==1,,1]
newprob[pgm==1,,3] <- prob[pgm==1,,2]
}
return(newprob)
}
else { # pairprob
dimpairprob <- dim(pairprob)
dimpairprob[3] <- dimpairprob[4] <- length(genonames)
newpairprob <- array(0,dim=dimpairprob)
newpairprob[pgm==0,,1:2,1:2] <- pairprob[pgm==0,,,]
if(expandX!="full") { # simple/standard
newpairprob[pgm==1,,3,3] <- pairprob[pgm==1,,1,1]
newpairprob[pgm==1,,3,2] <- pairprob[pgm==1,,1,2]
newpairprob[pgm==1,,2,3] <- pairprob[pgm==1,,2,1]
newpairprob[pgm==1,,2,2] <- pairprob[pgm==1,,2,2]
}
else {
newpairprob[pgm==1,,4,4] <- pairprob[pgm==1,,1,1]
newpairprob[pgm==1,,4,3] <- pairprob[pgm==1,,1,2]
newpairprob[pgm==1,,3,4] <- pairprob[pgm==1,,2,1]
newpairprob[pgm==1,,3,3] <- pairprob[pgm==1,,2,2]
}
return(newpairprob)
}
}
}
else if(all(sex==1) && !force) { # all males
if(!missing(geno)) return(geno)
else if(!missing(draws)) return(draws)
else if(!missing(pairprob)) return(pairprob)
else {
dimnames(prob)[[3]] <- genonames
return(prob)
}
}
else { # both sexes
if(length(pgm)==0 || all(pgm==0)) { # both sexes, forw dir
if(!missing(geno)) {
gmale <- geno[sex==1,]
if(expandX=="simple")
gmale[!is.na(gmale) & gmale==2] <- 3
else {
gmale[!is.na(gmale) & gmale==1] <- 3
gmale[!is.na(gmale) & gmale==2] <- 4
}
geno[sex==1,] <- gmale
return(geno)
}
else if(!missing(draws)) {
gmale <- draws[sex==1,,]
if(expandX=="simple")
gmale[gmale==2] <- 3
else {
gmale[gmale==1] <- 3
gmale[gmale==2] <- 4
}
draws[sex==1,,] <- gmale
return(draws)
}
else if(!missing(prob)) {
dimprob <- dim(prob)
dimprob[3] <- length(genonames)
newprob <- array(0,dim=dimprob)
dimnames(newprob) <- c(dimnames(prob)[1:2],list(genonames))
newprob[sex==0,,1:2] <- prob[sex==0,,1:2]
if(expandX=="simple") {
newprob[sex==1,,1] <- prob[sex==1,,1]
newprob[sex==1,,3] <- prob[sex==1,,2]
}
else {
newprob[sex==1,,3] <- prob[sex==1,,1]
newprob[sex==1,,4] <- prob[sex==1,,2]
}
return(newprob)
}
else { # pairprob
dimpairprob <- dim(pairprob)
dimpairprob[3] <- dimpairprob[4] <- length(genonames)
newpairprob <- array(0,dim=dimpairprob)
newpairprob[sex==0,,1:2,1:2] <- pairprob[sex==0,,,]
if(expandX=="simple") {
newpairprob[sex==1,,1,1] <- pairprob[sex==1,,1,1]
newpairprob[sex==1,,1,3] <- pairprob[sex==1,,1,2]
newpairprob[sex==1,,3,1] <- pairprob[sex==1,,2,1]
newpairprob[sex==1,,3,3] <- pairprob[sex==1,,2,2]
}
else {
newpairprob[sex==1,,3,3] <- pairprob[sex==1,,1,1]
newpairprob[sex==1,,3,4] <- pairprob[sex==1,,1,2]
newpairprob[sex==1,,4,3] <- pairprob[sex==1,,2,1]
newpairprob[sex==1,,4,4] <- pairprob[sex==1,,2,2]
}
return(newpairprob)
}
} # both sexes, forw dir
if(all(pgm==1) && !force) { # both sexes, backw dir
if(!missing(geno)) {
gmale <- geno[sex==1,]
if(expandX=="simple") {
gmale[!is.na(gmale) & gmale==1] <- 3
gmale[!is.na(gmale) & gmale==2] <- 1
}
else {
gmale[!is.na(gmale) & gmale==1] <- 3
gmale[!is.na(gmale) & gmale==2] <- 4
}
geno[sex==1,] <- gmale
return(geno)
}
else if(!missing(draws)) {
gmale <- draws[sex==1,,]
if(expandX=="simple") {
gmale[gmale==1] <- 3
gmale[gmale==2] <- 1
}
else {
gmale[gmale==1] <- 3
gmale[gmale==2] <- 4
}
draws[sex==1,,] <- gmale
return(draws)
}
else if(!missing(prob)) {
dimprob <- dim(prob)
dimprob[3] <- length(genonames)
newprob <- array(0,dim=dimprob)
dimnames(newprob) <- c(dimnames(prob)[1:2],list(genonames))
newprob[sex==0,,1:2] <- prob[sex==0,,1:2]
if(expandX=="simple") {
newprob[sex==1,,3] <- prob[sex==1,,1]
newprob[sex==1,,1] <- prob[sex==1,,2]
}
else {
newprob[sex==1,,3] <- prob[sex==1,,1]
newprob[sex==1,,4] <- prob[sex==1,,2]
}
return(newprob)
}
else { # pairprob
dimpairprob <- dim(pairprob)
dimpairprob[3] <- dimpairprob[4] <- length(genonames)
newpairprob <- array(0,dim=dimpairprob)
newpairprob[sex==0,,1:2,1:2] <- pairprob[sex==0,,,]
if(expandX=="simple") {
newpairprob[sex==1,,3,3] <- pairprob[sex==1,,1,1]
newpairprob[sex==1,,1,3] <- pairprob[sex==1,,2,1]
newpairprob[sex==1,,3,1] <- pairprob[sex==1,,1,2]
newpairprob[sex==1,,1,1] <- pairprob[sex==1,,2,2]
}
else {
newpairprob[sex==1,,3,3] <- pairprob[sex==1,,1,1]
newpairprob[sex==1,,3,4] <- pairprob[sex==1,,1,2]
newpairprob[sex==1,,4,3] <- pairprob[sex==1,,2,1]
newpairprob[sex==1,,4,4] <- pairprob[sex==1,,2,2]
}
return(newpairprob)
}
} # both sexes, backw dir
else { # both dir, both sexes
if(!missing(geno)) {
gmale <- geno[sex==1,]
gfemaler <- geno[sex==0 & pgm==1,]
if(expandX=="simple") {
gmale[!is.na(gmale) & gmale==2] <- 3
gfemaler[!is.na(gfemaler) & gfemaler==1] <- 3
}
else if(expandX=="standard") {
gmale[!is.na(gmale) & gmale==1] <- 4
gmale[!is.na(gmale) & gmale==2] <- 5
gfemaler[!is.na(gfemaler) & gfemaler==1] <- 3
}
else {
gmale[!is.na(gmale) & gmale==1] <- 5
gmale[!is.na(gmale) & gmale==2] <- 6
gfemaler[!is.na(gfemaler) & gfemaler==1] <- 4
gfemaler[!is.na(gfemaler) & gfemaler==2] <- 3
}
geno[sex==1,] <- gmale
geno[sex==0 & pgm==1,] <- gfemaler
return(geno)
}
else if(!missing(draws)) {
gmale <- draws[sex==1,,]
gfemaler <- draws[sex==0 & pgm==1,,]
if(expandX=="simple") {
gmale[gmale==2] <- 3
gfemaler[gfemaler==1] <- 3
}
else if(expandX=="standard") {
gmale[gmale==1] <- 4
gmale[gmale==2] <- 5
gfemaler[gfemaler==1] <- 3
}
else {
gmale[gmale==1] <- 5
gmale[gmale==2] <- 6
gfemaler[gfemaler==1] <- 4
gfemaler[gfemaler==2] <- 3
}
draws[sex==1,,] <- gmale
draws[sex==0 & pgm==1,,] <- gfemaler
return(draws)
}
else if(!missing(prob)) {
dimprob <- dim(prob)
dimprob[3] <- length(genonames)
newprob <- array(0,dim=dimprob)
dimnames(newprob) <- c(dimnames(prob)[1:2],list(genonames))
newprob[sex==0 & pgm==0,,1:2] <- prob[sex==0 & pgm==0,,1:2]
if(expandX=="simple") {
newprob[sex==1,,1] <- prob[sex==1,,1]
newprob[sex==1,,3] <- prob[sex==1,,2]
newprob[sex==0 & pgm==1,,3] <- prob[sex==0 & pgm==1,,1]
newprob[sex==0 & pgm==1,,2] <- prob[sex==0 & pgm==1,,2]
}
else if(expandX=="standard") {
newprob[sex==1,,4] <- prob[sex==1,,1]
newprob[sex==1,,5] <- prob[sex==1,,2]
newprob[sex==0 & pgm==1,,3] <- prob[sex==0 & pgm==1,,1]
newprob[sex==0 & pgm==1,,2] <- prob[sex==0 & pgm==1,,2]
}
else {
newprob[sex==1,,5] <- prob[sex==1,,1]
newprob[sex==1,,6] <- prob[sex==1,,2]
newprob[sex==0 & pgm==1,,4] <- prob[sex==0 & pgm==1,,1]
newprob[sex==0 & pgm==1,,3] <- prob[sex==0 & pgm==1,,2]
}
return(newprob)
}
else { # pairprob
dimpairprob <- dim(pairprob)
dimpairprob[3] <- dimpairprob[4] <- length(genonames)
newpairprob <- array(0,dim=dimpairprob)
newpairprob[sex==0 & pgm==0,,1:2,1:2] <- pairprob[sex==0 & pgm==0,,,]
male <- (sex==1)
femaler <- (sex==0) & (pgm==1)
if(expandX=="simple") {
newpairprob[male,,1,1] <- pairprob[male,,1,1]
newpairprob[male,,1,3] <- pairprob[male,,1,2]
newpairprob[male,,3,1] <- pairprob[male,,2,1]
newpairprob[male,,3,3] <- pairprob[male,,2,2]
newpairprob[femaler,,3,3] <- pairprob[femaler,,1,1]
newpairprob[femaler,,3,2] <- pairprob[femaler,,1,2]
newpairprob[femaler,,2,3] <- pairprob[femaler,,2,1]
newpairprob[femaler,,2,2] <- pairprob[femaler,,2,2]
}
else if(expandX=="standard") {
newpairprob[male,,4,4] <- pairprob[male,,1,1]
newpairprob[male,,4,5] <- pairprob[male,,1,2]
newpairprob[male,,5,4] <- pairprob[male,,2,1]
newpairprob[male,,5,5] <- pairprob[male,,2,2]
newpairprob[femaler,,3,3] <- pairprob[femaler,,1,1]
newpairprob[femaler,,3,2] <- pairprob[femaler,,1,2]
newpairprob[femaler,,2,3] <- pairprob[femaler,,2,1]
newpairprob[femaler,,2,2] <- pairprob[femaler,,2,2]
}
else {
newpairprob[male,,5,5] <- pairprob[male,,1,1]
newpairprob[male,,5,6] <- pairprob[male,,1,2]
newpairprob[male,,6,5] <- pairprob[male,,2,1]
newpairprob[male,,6,6] <- pairprob[male,,2,2]
newpairprob[femaler,,4,4] <- pairprob[femaler,,1,1]
newpairprob[femaler,,4,3] <- pairprob[femaler,,1,2]
newpairprob[femaler,,3,4] <- pairprob[femaler,,2,1]
newpairprob[femaler,,3,3] <- pairprob[femaler,,2,2]
}
return(newpairprob)
}
}
}
} # end of intercross
}
######################################################################
# scanoneXnull
#
# figure out null hypothesis business for scanone on X chromosome
######################################################################
scanoneXnull <-
function(type, sexpgm, cross.attr)
{
sex <- sexpgm$sex
pgm <- sexpgm$pgm
if(type == "risib" || type=="riself" || type=="dh" || type=="haploid") type <- "bc"
if(type == "bcsft") {
if(cross.attr$scheme[2] == 0)
type <- "bc"
else
type <- "f2"
}
### first figure out sex/pgm pattern
# sex
if(length(sex)==0 || all(sex==0)) { # all female
onesex <- allfemale <- TRUE
}
else if(all(sex==1)) { # all male
onesex <- TRUE
allfemale <- FALSE
}
else { # both sexes
onesex <- allfemale <- FALSE
}
# pgm
if(length(pgm)==0 || all(pgm==0) || all(pgm==1)) # one direction
onedir <- TRUE
else onedir <- FALSE
allmale <- onesex && !allfemale
bothsex <- !onesex
bothdir <- !onedir
### now figure out the null hypothesis and pull out appropriate
### covariates for the null
# backcross, one sex
# OR intercross, one dir and one sex
# OR intercross, both dir and all male
if((type=="bc" && onesex) ||
(type=="f2" && ((onedir && onesex) || (bothdir && allmale)))) {
adjustX <- FALSE
parX0 <- 1
sexpgmcovar <- sexpgmcovar.alt <- NULL
}
# backcross, both sexes
# OR intercross, one direction and both sexes
else if((type=="bc" && bothsex) ||
(type=="f2" && onedir && bothsex)) {
adjustX <- TRUE
parX0 <- 2
sexpgmcovar <- cbind(sex)
sexpgmcovar.alt <- sex+1
}
# intercross, both dir and all female
else if(type=="f2" && bothdir && allfemale) {
adjustX <- TRUE
parX0 <- 2
sexpgmcovar <- cbind(pgm)
sexpgmcovar.alt <- pgm+1
}
# intercross, both dir and both sexes
else {
adjustX <- TRUE
parX0 <- 3
sexpgmcovar <- cbind(sex,as.numeric(sex==0 & pgm==1))
sexpgmcovar.alt <- rep(3,length(sex))
sexpgmcovar.alt[sex==0 & pgm==0] <- 1
sexpgmcovar.alt[sex==0 & pgm==1] <- 2
}
list(adjustX=adjustX, parX0=parX0, sexpgmcovar=sexpgmcovar,
sexpgmcovar.alt=sexpgmcovar.alt)
}
######################################################################
# revisecovar
#
# Drop sex and pgm and their interxn as covariates for the X chr.
######################################################################
revisecovar <-
function(sexpgm, covar)
{
if(is.null(covar) || (is.null(sexpgm$sex) && is.null(sexpgm$pgm))) {
if(!is.null(covar)) attr(covar, "n.dropped") <- 0
return(covar)
}
covar <- as.matrix(covar)
sex <- sexpgm$sex
pgm <- sexpgm$pgm
if(!is.null(pgm) && length(unique(pgm))==1) pgm <- NULL
allfemale <- FALSE
if(is.null(sex)) allfemale <- TRUE
else {
if(all(sex==0)) {
allfemale <- TRUE
sex <- NULL
}
else if(all(sex==1)) {
allfemale <- FALSE
sex <- NULL
}
}
if(!is.null(pgm)) { # some of each direction
if(!is.null(sex)) { # some of each sex
femf <- as.numeric(pgm==0 & sex==0)
femr <- as.numeric(pgm==1 & sex==0)
mal <- sex
X <- cbind(femf, femr, mal)
}
else { # all of one sex
if(allfemale)
X <- cbind(1-pgm, pgm)
else
X <- cbind(rep(1, nrow(covar)))
}
}
else { # all of one direction
if(!is.null(sex)) # some of each sex
X <- cbind(sex, 1-sex)
else X <- cbind(rep(1, nrow(covar)))
}
nc <- ncol(X)
keep <- rep(TRUE,ncol(covar))
for(i in 1:ncol(covar)) {
if(qr(cbind(X,covar[,i]))$rank <= nc)
keep[i] <- FALSE
}
if(!any(keep))
covar <- numeric(0)
else
covar <- covar[,keep,drop=FALSE]
attr(covar, "n.dropped") <- sum(!keep)
covar
}
######################################################################
# dropXcol: for use with scantwo() for the X chromosome:
# figure out what columns to drop...both for the full model
# and for the additive model.
######################################################################
dropXcol <-
function(type=c("f2","bc", "riself", "risib", "4way", "dh", "haploid", "special","bcsft"),
sexpgm, cross.attr)
{
type <- match.arg(type)
## Treat bcsft as bc if no intercross generations; otherwise as f2.
if(type == "bcsft") {
if(cross.attr$scheme[2] == 0)
type <- "bc"
else
type <- "f2"
}
gn <- getgenonames(type, "X", "full", sexpgm, cross.attr)
if(length(gn)==2) return(rep(0,4))
if(length(gn)==4) return( c(0,0,0,0,0,1,0, 0,1,1,1,1,1,1,1,0) )
if(length(gn)==6) {
todrop <- c(rep(0,11), rep(1,25))
todrop[c(8,10)] <- 1
todrop[11+c(1,13,25)] <- 0
return(todrop)
}
return(rep(0,length(gn)^2))
}
# end of xchr.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.