#####################################################################
#
# mqmaugment.R
#
# Copyright (c) 2009-2019, Danny Arends
#
# Modified by Pjotr Prins; slight modification by Karl Broman
#
#
# first written Februari 2009
# last modified Dec 2019
#
# 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: mqmaugment
#
#
#####################################################################
######################################################################
#
# mqmaugment_on_cofactors: Data Augmentation routine for MQM only using the cofactors, the other markers are filled by fill.geno()
#
######################################################################
mqmaugment_on_cofactors <- function(cross, cofactors, maxaugind=82, minprob=0.1, strategy=c("default","impute","drop"), verbose=FALSE, ...){
markernames <- as.character(unlist(lapply(pull.map(cross),names)))
todrop <- markernames[which(cofactors==0)]
toaugment <- drop.markers(cross,todrop)
cross <- fill.geno(cross, ...)
augmented <- mqmaugment(toaugment,maxaugind,minprob,strategy,verbose)
newgenomatrix <- NULL
for(marker in markernames){
if(marker %in% todrop){
newgenomatrix <- cbind(newgenomatrix,pull.geno(cross)[augmented$mqm[[3]],marker])
}else{
newgenomatrix <- cbind(newgenomatrix,pull.geno(augmented)[,marker])
}
}
colnames(newgenomatrix) <- markernames
augmented$geno <- vector("list",nchr(cross))
for(chr in 1:length(cross$geno)){
cat(chr,"\n")
augmented$geno[[chr]]$map <- cross$geno[[chr]]$map
augmented$geno[[chr]]$data <- newgenomatrix[,names(cross$geno[[chr]]$map)]
names(augmented$geno) <- names(cross$geno)
class(augmented$geno[[chr]]) <- chrtype(cross$geno[[chr]])
}
augmented
}
######################################################################
#
# mqmaugment: dataaugmentation routine for MQM
#
######################################################################
mqmaugment <- function(cross,maxaugind=82, minprob=0.1, strategy=c("default","impute","drop"), verbose=FALSE) {
# ---- check input supplied by user, start a timer
if(minprob <= 0 || minprob > 1){ stop("Error minprob should be a value between 0 and 1.") }
starttime <- proc.time()
maxaug <- nind(cross) * maxaugind # maxaug is the maximum of individuals to augment to
supported <- c("default","impute","drop")
strategy <- pmatch(strategy, supported)
# ---- check for supported crosses and set ctype
isF2 = 1
isBC = 2
isRIL = 3
isAA = 1
isAB = 2
isH = 2
isBB = 3
isNOTBB = 4
isNOTAA = 5
isMISSING = 9
crosstype <- crosstype(cross)
if (crosstype == "f2") {
ctype = isF2
}
else if (crosstype == "bc" || crosstype == "dh" || crosstype=="haploid") {
ctype = isBC
}
else if (crosstype == "riself") {
ctype = isRIL
}
else {
stop("Currently only F2, BC, and selfed RIL crosses can be analyzed by MQM.")
}
if (verbose) cat("INFO: Received a valid cross file type:", crosstype,".\n")
# ---- Check sex chromosome
# check whether the X chromosome should be dropped
# (backcross with one sex should be fine)
chr_type <- sapply(cross$geno, chrtype)
# Drop the X chromosome in F2 and related crosses
if (any(chr_type == "X") && (ctype == isF2 || length(getgenonames(crosstype, "X", "full", getsex(cross), attributes(cross))) != 2)) {
warning("MQM not yet available for the X chromosome; omitting chr ",
paste(names(cross$geno)[chr_type == "X"], collapse=" "))
cross <- subset(cross, chr=(chr_type != "X"))
}
# ---- Count
n.ind <- nind(cross)
n.chr <- nchr(cross)
n.aug <- maxaug
if (verbose) {
cat("INFO: Number of individuals:",n.ind,".\n")
cat("INFO: Number of chr:",n.chr,".\n")
}
# ---- Genotype
geno <- pull.geno(cross)
chr <- rep(1:nchr(cross), nmar(cross))
dist <- unlist(pull.map(cross))
# Create a fake phenotype for augmentation
pheno <- rep(1:n.ind)
n.mark <- ncol(geno)
if (verbose) cat("INFO: Number of markers:",n.mark,".\n")
# Check for NA genotypes and replace them with a 9
geno[is.na(geno)] <- isMISSING
if (ctype==isRIL) {
nH = sum(geno==isH)
if (nH>0) {
#warning("RIL dataset contains ", nH," heterozygous genotypes")
if (any(geno==isBB)) { # have 3/BB's, so replace 2/H's with missing values
geno[geno==isH] <- isMISSING
#warning("Removed heterozygous genotypes from RIL set")
} else {
#warning("Converting heterozygous genotypes to BB from RIL set")
geno[geno==isH] <- isBB
}
}
} # end if(RIL)
# ---- Call data augmentation
result <- .C("R_mqmaugment",
as.integer(geno),
as.double(dist),
as.double(pheno),
augGeno=as.integer(rep(0,n.mark*maxaug)),
augPheno=as.double(rep(0,maxaug)),
augIND=as.integer(rep(0,maxaugind*n.ind)),
nind=as.integer(n.ind),
naug=as.integer(n.aug),
as.integer(n.mark),
as.integer(1), # 1 phenotype
as.integer(maxaug),
as.integer(maxaugind),
as.double(minprob),
as.integer(chr),
as.integer(ctype),
as.integer(strategy),
as.integer(verbose),
PACKAGE="qtl")
n.indold = n.ind
n.ind = result$nind
n.aug = result$naug
markONchr <- 0
markdone <- 0
pheno <- NULL
pnames <- phenames(cross)
oldpheno <- pull.pheno(cross)
result$augIND <- result$augIND+1
for(x in result$augPheno[1:n.aug]){
if(nphe(cross)>1){
pheno <- rbind(pheno,oldpheno[x,])
}else{
pheno <- c(pheno,oldpheno[x])
}
}
for(c in 1:n.chr){
#print(paste("Cromosome",c,"\n",sep=""))
matri <- NULL
markONchr <- dim(cross$geno[[c]]$data)[2]
#print(paste("# markers",markONchr,"\n",sep=""))
for(j in markdone:(markdone+markONchr-1)){
#print(paste("Start",markdone,":End",(markdone+markONchr-1),"\n",sep=""))
ind2 <- NULL
ind2 <- result$augGeno[(1+(j*maxaug)):(n.aug+(j*maxaug))]
matri <- rbind(matri,ind2)
}
matri <- t(matri)
if(markdone==0){
colnames(matri) <- colnames(geno)[markdone:(markdone+markONchr)]
} else {
#print(paste("Markdone",markdone,"End",(markdone+markONchr-1)))
colnames(matri) <- colnames(geno)[(markdone+1):(markdone+markONchr)]
}
cross$geno[[c]]$data <- matri
markdone <- (markdone+markONchr)
}
if(nphe(cross) > 1){ colnames(pheno) <- colnames(cross$pheno) } # Add the phenotype names if we have multiple phenotypes
cross$pheno <- as.data.frame(pheno, stringsAsFactors=TRUE)
# Store extra information (needed by the MQM algorithm) which individual was which original etc..
cross$mqm$Nind <- n.ind
cross$mqm$Naug <- n.aug
result$augIND <- result$augIND-1
cross$mqm$augIND <- result$augIND[1:n.aug]
# ---- RESULTS
endtime <- proc.time()
if(n.ind != n.indold){
if(verbose) warning("SERIOUS WARNING: Dropped ",abs(n.ind - n.indold)," original individuals.\n Information lost, please increase minprob.")
}
if(verbose) cat("INFO: DATA-Augmentation took: ",round((endtime-starttime)[3], digits=3)," seconds\n")
colnames(cross$pheno) <- pnames # fix up phenotype names problem
cross # return cross type
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.