# Author: Kosuke Hamazaki hamazaki@ut-biomet.org
# 2020 The University of Tokyo
#
# Description:
# Definition of traitInfo class
#' R6 Class representing information of traits (especially simulated ones)
#'
#' @description traitInfo object store specific information of traits (especially simulated ones).
#'
#' @export
#' @import R6
traitInfo <- R6::R6Class(
"traitInfo",
public = list(
#' @field lociInfo [lociInfo class] lociInfo of the specie
lociInfo = NULL,
#' @field nMarkers [numeric] Number of markers for each chromosome
nMarkers = NULL,
#' @field nTraits [numeric] Number of traits assumed in this breeding scheme
nTraits = NULL,
#' @field traitNames [character] Names of traits
traitNames = NULL,
#' @field nQTLs [matrix] Number of QTLs of each chromosome for each trait
nQTLs = NULL,
#' @field qtlVarList [list] A list of variances of QTLs for additive/dominance, pleiotropy, and epistasis
#'
#' \itemize{
#' \item{\code{$qtlPle}:} {A variance of QTLs for pleiotropy}
#' \item{\code{$qtlEach}:} {A variance of QTLs for each trait (except epistasis)}
#' \item{\code{$qtlEpi}:} {A variance of QTLs for each trait (epistasis)}
#' }
qtlVarList = NULL,
#' @field qtlOverlap [logical] If TRUE, pleiotropy across traits will be assumed
qtlOverlap = NULL,
#' @field nOverlap [numeric] Number of QTLs common across traits for each chromosome
nOverlap = NULL,
#' @field effCor [numeric] Correlation of pleiotropic QTL effects between traits
effCor = NULL,
#' @field propDomi [numeric] Proportion of dominant effects for each trait
propDomi = NULL,
#' @field interactionMean [numeric] The expected number that each QTL interacts with others for each trait
interactionMean = NULL,
#' @field actionTypeEpiSimple [logical] If TRUE, `actionTypeEpi` will be same as
#' `actionType` for each QTL.
actionTypeEpiSimple = NULL,
#' @field seedSelectMarker [numeric] Random seed for selecting marker information
seedSelectMarker = NA,
#' @field seedSelectQTLEach [numeric] Random seed for selecting each QTL and its effect (it should have the same length as nTraits)
seedSelectQTLEach = NA,
#' @field seedSelectQTLPle [numeric] Random seed for selecting pleiotropic QTLs and their effects
seedSelectQTLPle = NA,
#' @field seedSelectQTLDom [numeric] Random seed for selecting QTL positions with dominant effects (it should have the same length as nTraits)
seedSelectQTLDom = NA,
#' @field seedSelectQTLEpi [numeric] Random seed for selecting epistasis and its effect (it should have the same length as nTraits)
seedSelectQTLEpi = NA,
#' @field mrkPos [numeric] Marker positions
mrkPos = NULL,
#' @field mrkPosEachChr [list] Marker positions for each chromosome
mrkPosEachChr = NULL,
#' @field qtlPos [list] QTL positions for each trait
qtlPos = NULL,
#' @field qtlEachPos [list] Positions of QTLs specific to one trait for each trait
qtlEachPos = NULL,
#' @field qtlPlePos [list] Pleiotropic QTL positions for each trait
qtlPlePos = NULL,
#' @field qtlEpiPos [list] Epistatic QTL positions for each trait
qtlEpiPos = NULL,
#' @field qtlPosEachChrList [list] QTL positions for each trait for each chromosome
qtlPosEachChrList = NULL,
#' @field qtlEachPosEachChrList [list] Positions of QTLs specific to one trait for each trait for each chromosome
qtlEachPosEachChrList = NULL,
#' @field qtlPlePosEachChrList [list] Pleiotropic QTL positions for each trait for each chromosome
qtlPlePosEachChrList = NULL,
#' @field qtlNamesList [list] QTL names for each trait
qtlNamesList = NULL,
#' @field qtlEachNamesList [list] Names of QTLs common to one trait for each trait
qtlEachNamesList = NULL,
#' @field qtlPleNamesList [list] Pleiotropic QTL names for each trait
qtlPleNamesList = NULL,
#' @field qtlEpiNamesList [list] Epistatic QTL names for each trait
qtlEpiNamesList = NULL,
#' @field qtlAllNamesList [list] All QTL names for each trait (including epistatic effects)
qtlAllNamesList = NULL,
#' @field actionType [list] A list of vectors representing additive by 0 and dominance by 1 for each trait
actionType = NULL,
#' @field actionTypeEpi [list] A list of matrix representing additive by 0 and dominance by 1 for epitasis for each trait
actionTypeEpi = NULL,
#' @field qtlEff [list] QTL effects for each trait
qtlEff = NULL,
#' @field qtlPleEff [list] Pleiotropic QTL effects for each trait
qtlPleEff = NULL,
#' @field qtlEachEff [list] Effects of QTLs specific to one trait for each trait
qtlEachEff = NULL,
#' @field qtlEpiEff [list] Epistatic QTL effects for each trait
qtlEpiEff = NULL,
#' @field qtlAllEff [list] All QTL effects for each trait (including epistasis)
qtlAllEff = NULL,
#' @field nQTLsEach [numeric] Number of QTLs specific to each trait
nQTLsEach = NULL,
#' @field nEpi [numeric] Number of epistatic effects
nEpi = NULL,
#' @description Create a new traitInfo object.
#' @param lociInfo [lociInfo class] lociInfo of the specie
#' @param nMarkers [numeric] Number of markers for each chromosome
#' @param nTraits [numeric] Number of traits assumed in this breeding scheme
#' @param traitNames [character] Names of traits
#' @param nQTLs [matrix] Number of QTLs of each chromosome for each trait
#' @param qtlVarList [list] A list of variances of QTLs for additive/dominance, pleiotropy, and epistasis
#'
#' \itemize{
#' \item{\code{$qtlPle}:} {A variance of QTLs for pleiotropy}
#' \item{\code{$qtlEach}:} {A variance of QTLs for each trait (except epistasis)}
#' \item{\code{$qtlEpi}:} {A variance of QTLs for each trait (epistasis)}
#' }
#' @param qtlOverlap [logical] If TRUE, pleiotropy across traits will be assumed
#' @param nOverlap [numeric] Number of QTLs common across traits for each chromosome
#' @param effCor [numeric] Correlation of pleiotropic QTL effects between traits
#' @param propDomi [numeric] Proportion of dominant effects for each trait
#' @param interactionMean [numeric] The expected number that each QTL interacts with others for each trait
#' @param actionTypeEpiSimple [logical] If TRUE, `actionTypeEpi` will be same as
#' `actionType` for each QTL.
#' @param qtlPleNamesList [list] Pleiotropic QTL names for each trait
#' @param qtlEachNamesList [list] Names of QTLs common to one trait for each trait
#' @param seedSelectMarker [numeric] Random seed for selecting marker information
#' @param seedSelectQTLEach [numeric] Random seed for selecting each QTL and its effect (it should have the same length as nTraits)
#' @param seedSelectQTLPle [numeric] Random seed for selecting pleiotropic QTLs and their effects
#' @param seedSelectQTLDom [numeric] Random seed for selecting QTL positions with dominant effects (it should have the same length as nTraits)
#' @param seedSelectQTLEpi [numeric] Random seed for selecting epistasis and its effect (it should have the same length as nTraits)
#' @return A new `traitInfo` object.
#' @examples
#' mySimInfo <- simInfo$new(simName = "Simulation Example",
#' simGeno = TRUE,
#' simPheno = TRUE,
#' nSimGeno = 1,
#' nSimPheno = 3,
#' nCoreMax = 4,
#' nCorePerGeno = 1,
#' nCorePerPheno = 3,
#' saveDataFileName = NULL)
#'
#' ### create specie information
#' mySpec <- specie$new(nChr = 3,
#' lChr = c(100, 150, 200),
#' specName = "Example 1",
#' ploidy = 2,
#' mutRate = 10^-8,
#' recombRate = 10^-7,
#' recombRateOneVal = FALSE,
#' chrNames = c("C1", "C2", "C3"),
#' nLoci = 100,
#' effPopSize = 100,
#' simInfo = mySimInfo,
#' verbose = TRUE)
#'
#' ### create SNP information
#' myLoci <- lociInfo$new(genoMap = NULL,
#' specie = mySpec,
#' seedSimGeno = NA,
#' lociNames = NULL,
#' founderNames = NULL)
#'
#'
#' ### create traitInfo object
#' myTrait <- traitInfo$new(lociInfo = myLoci,
#' nMarkers = 80,
#' nTraits = 3,
#' nQTLs = c(4, 8, 3),
#' actionTypeEpiSimple = TRUE,
#' qtlOverlap = TRUE,
#' nOverlap = c(2, 3, 0),
#' effCor = 0.1,
#' propDomi = 0.2,
#' interactionMean = c(4, 1, 2))
#' print(myTrait)
#'
initialize = function(lociInfo,
nMarkers,
nTraits = NULL,
traitNames = NULL,
nQTLs = NULL,
qtlVarList = NULL,
qtlOverlap = TRUE,
nOverlap = NULL,
effCor = NULL,
propDomi = NULL,
interactionMean = NULL,
actionTypeEpiSimple = TRUE,
qtlPleNamesList = NULL,
qtlEachNamesList = NULL,
seedSelectMarker = NA,
seedSelectQTLEach = NA,
seedSelectQTLPle = NA,
seedSelectQTLDom = NA,
seedSelectQTLEpi = NA
) {
# CHECKS:
# parameters classes
if (class(lociInfo)[1] != "lociInfo") {
stop('"class(lociInfo)" must be "lociInfo"')
}
nChr <- lociInfo$specie$nChr
chrNames <- lociInfo$specie$chrNames
nLoci <- lociInfo$specie$nLoci
qtlPosCand <- sapply(1:nChr, function(chrNo) {
qtlPosCandEachChr <- 1:nLoci[chrNo]
names(qtlPosCandEachChr) <- lociInfo$ids[[chrNo]]
return(qtlPosCandEachChr)
},
simplify = FALSE)
names(qtlPosCand) <- chrNames
if (!is.numeric(nMarkers)) stop("nMarkers must be numeric.")
if (any(floor(nMarkers) - nMarkers != 0)) stop("nMarkers must be integers.")
if (length(nMarkers) != 1 && length(nMarkers) != nChr) {
stop(paste("length(nMarkers) must be equal to 1 (all chr have the same",
"size) or equal to nChr."))
} else if (length(nMarkers) == 1) {
nMarkers <- rep(nMarkers, nChr)
}
stopifnot(!any(nMarkers < 0))
names(nMarkers) <- chrNames
if (lociInfo$specie$simInfo$simPheno) {
if (is.null(nTraits)) {
message('"nTraits" was not specified. The nTraits had been set to "1"')
nTraits <- 1
}
if (!is.numeric(nTraits)) stop("nTraits must be numeric.")
if (any(floor(nTraits) - nTraits != 0)) stop("nTraits must be integers.")
stopifnot(nTraits >= 1)
if (is.null(traitNames)) {
traitNames <- .charSeq(paste0("Trait_"), seq(nTraits))
}
stopifnot(length(traitNames) == nTraits)
if (is.null(propDomi)) {
message('"propDomi" was not specified. The propDomi had been set to "0".')
propDomi <- rep(0, nTraits)
}
if (!is.numeric(propDomi)) stop("propDomi must be numeric.")
if (length(propDomi) != 1 && length(propDomi) != nTraits) {
stop(paste("length(propDomi) must be equal to 1 (all traits have the same",
"value) or equal to nTraits."))
} else if (length(propDomi) == 1) {
propDomi <- rep(propDomi, nTraits)
}
if (any((propDomi > 1) | (propDomi < 0))) {
stop("All elements in `propDomi` should be within the interval [0, 1].")
}
if (is.null(interactionMean)) {
message('"interactionMean" was not specified. The interactionMean had been set to "0".')
interactionMean <- rep(0, nTraits)
}
if (!is.numeric(interactionMean)) stop("interactionMean must be numeric.")
if (any(floor(interactionMean) - interactionMean != 0)) stop("interactionMean must be integers.")
if (length(interactionMean) != 1 && length(interactionMean) != nTraits) {
stop(paste("length(interactionMean) must be equal to 1 (all traits have the same",
"value) or equal to nTraits."))
} else if (length(interactionMean) == 1) {
interactionMean <- rep(interactionMean, nTraits)
}
stopifnot(!any(propDomi < 0))
stopifnot(!any(interactionMean < 0))
names(propDomi) <- names(interactionMean) <- traitNames
if (!is.null(seedSelectMarker)) {
if (!is.na(seedSelectMarker)) {
stopifnot(is.numeric(seedSelectMarker))
seedSelectMarker <- floor(seedSelectMarker)
} else {
seedSelectMarker <- sample(x = 1e9, size = 1)
}
}
if (!is.null(seedSelectQTLEach)) {
if (!is.na(seedSelectQTLEach)) {
stopifnot(is.numeric(seedSelectQTLEach))
stopifnot(length(seedSelectQTLEach) == nTraits)
seedSelectQTLEach <- floor(seedSelectQTLEach)
} else {
seedSelectQTLEach <- sample(x = 1e9, size = nTraits)
}
}
if (!is.null(seedSelectQTLPle)) {
if (!is.na(seedSelectQTLPle)) {
stopifnot(is.numeric(seedSelectQTLPle))
seedSelectQTLPle <- floor(seedSelectQTLPle)
} else {
seedSelectQTLPle <- sample(x = 1e9, size = 1)
}
}
if (!is.null(seedSelectQTLDom)) {
if (!is.na(seedSelectQTLDom)) {
stopifnot(is.numeric(seedSelectQTLDom))
stopifnot(length(seedSelectQTLDom) == nTraits)
seedSelectQTLDom <- floor(seedSelectQTLDom)
} else {
seedSelectQTLDom <- sample(x = 1e9, size = nTraits)
}
}
if (!is.null(seedSelectQTLEpi)) {
if (!is.na(seedSelectQTLEpi)) {
stopifnot(is.numeric(seedSelectQTLEpi))
stopifnot(length(seedSelectQTLEpi) == nTraits)
seedSelectQTLEpi <- floor(seedSelectQTLEpi)
} else {
seedSelectQTLEpi <- sample(x = 1e9, size = nTraits)
}
}
if (nTraits == 1) {
if (is.null(nQTLs)) {
message('"nQTLs" was not specified. The nQTLs had been set to "10" for all chromosome & traits.')
nQTLs <- rep(10, nChr)
}
if (!is.numeric(nQTLs)) stop("nQTLs must be numeric.")
if (any(floor(nQTLs) - nQTLs != 0)) stop("nQTLs must be integers.")
if (length(nQTLs) != 1 && length(nQTLs) != nChr) {
stop(paste("length(nQTLs) must be equal to 1 (all chr have the same",
"size) or equal to nChr."))
} else if (length(nQTLs) == 1) {
nQTLs <- rep(nQTLs, nChr)
}
names(nQTLs) <- chrNames
if (is.null(qtlVarList)) {
qtlVarList <- list(qtlEach = 1,
qtlEpi = 1)
}
if (!all(c("qtlEach", "qtlEpi") %in% names(qtlVarList))) {
stop("`qtlVarList` should contain objects names `qtlPle`, `qtlEach`, `qtlEpi`.")
} else {
effVarEach <- qtlVarList$qtlEach
effVarEpi <- qtlVarList$qtlEpi
if (!is.numeric(effVarEach)) stop("qtlEach must be numeric.")
if (!is.numeric(effVarEpi)) stop("qtlEpi must be numeric.")
}
names(effVarEach) <- names(effVarEpi) <- traitNames
qtlVarList <- list(qtlEach = effVarEach,
qtlEpi = effVarEpi)
qtlOverlap <- FALSE
nOverlap <- NULL
effCor <- NULL
nLociUsed <- nMarkers + nQTLs
nQTLsEach <- nQTLs
qtlPleNamesList <- list(NULL)
if (is.null(qtlEachNamesList)) {
qtlEachNamesList <- list(paste0("QTL_", 1:sum(nQTLsEach)))
}
stopifnot(length(qtlEachNamesList) == 1)
stopifnot(length(qtlEachNamesList[[1]]) == sum(nQTLsEach))
names(qtlPleNamesList) <- names(qtlEachNamesList) <- traitNames
qtlPleEff <- list(NULL)
qtlPlePos <- list(NULL)
qtlEachNames <- qtlEachNamesList[[1]]
set.seed(seed = seedSelectQTLEach)
qtlEachPosEachChr <- sapply(1:nChr, function(chrNo) {
sampleVec(x = qtlPosCand[[chrNo]], size = nQTLsEach[chrNo])
}, simplify = FALSE)
names(qtlEachPosEachChr) <- chrNames
set.seed(seed = seedSelectQTLEach)
qtlEachPos <- list(sampleVec(unlist(sapply(1:nChr, function(chrNo) {
qtlEachPosEachChr[[chrNo]] + cumsum(c(0, nLoci))[chrNo]
}, simplify = FALSE))))
names(qtlEachPos[[1]]) <- qtlEachNames
qtlPos <- qtlEachPos
qtlPosCandNow <- sapply(1:nChr, function(chrNo) {
qtlPosCandNowEachChr <- (qtlPosCand[[chrNo]])[!(names(qtlPosCand[[chrNo]]) %in%
names(qtlEachPosEachChr[[chrNo]]))]
return(qtlPosCandNowEachChr)
}, simplify = FALSE)
qtlEachPosEachChrList <- list(qtlEachPosEachChr)
qtlPlePosEachChrList <- list(NULL)
names(qtlEachPosEachChrList) <-
names(qtlPlePosEachChrList) <- traitNames
qtlPosEachChrList <- qtlEachPosEachChrList
set.seed(seed = seedSelectQTLEach)
qtlEachEff0 <- rnorm(n = sum(nQTLsEach), mean = 0,
sd = sqrt(effVarEach / sum(nQTLsEach)))
qtlEachEff <- list(qtlEachEff0[order(abs(qtlEachEff0), decreasing = TRUE)])
names(qtlEachEff[[1]]) <- qtlEachNames
set.seed(seed = seedSelectQTLEpi)
nEpi <- round(sum(rpois(n = sum(nQTLs), lambda = interactionMean)) / 2, 0)
names(nEpi) <- traitNames
if (nEpi > 0) {
set.seed(seed = seedSelectQTLEpi)
qtlEpiPos0 <- t(sapply(1:nEpi, function(x) sort(sampleVec(qtlPos[[1]], 2))))
qtlEpiPos0 <- qtlEpiPos0[!duplicated(qtlEpiPos0), , drop = FALSE]
nEpi <- nrow(qtlEpiPos0)
qtlEpiNames <- paste0(qtlEachNames[match((qtlEpiPos0)[, 1], qtlPos[[1]])], " x ",
qtlEachNames[match((qtlEpiPos0)[, 2], qtlPos[[1]])])
rownames(qtlEpiPos0) <- qtlEpiNames
colnames(qtlEpiPos0) <- paste0("Epi_", 1:2)
qtlEpiPos <- list(qtlEpiPos0)
qtlAllNames <- c(qtlEachNames, qtlEpiNames)
set.seed(seed = seedSelectQTLDom)
actionType <- list(rbinom(sum(nQTLs), 1, propDomi))
names(actionType[[1]]) <- qtlEachNames
if (!actionTypeEpiSimple) {
set.seed(seed = seedSelectQTLDom)
actionTypeEpi <- list(matrix(rbinom(2 * nEpi, 1, propDomi), nrow = nEpi))
} else {
actionTypeEpi <- list(matrix((actionType[[1]])[apply(qtlEpiPos0, 2, function(x) {
match(x, qtlPos[[1]])
})], nrow = nEpi))
}
rownames(actionTypeEpi[[1]]) <- qtlEpiNames
colnames(actionTypeEpi[[1]]) <- paste0("Epi_", 1:2)
set.seed(seed = seedSelectQTLEpi)
qtlEpiEff0 <- rnorm(n = nEpi, mean = 0, sd = sqrt(effVarEpi / nEpi))
qtlEpiEff <- list(qtlEpiEff0[order(abs(qtlEpiEff0), decreasing = TRUE)])
names(qtlEpiEff[[1]]) <- qtlEpiNames
qtlAllEff <- list(c(qtlEachEff[[1]], qtlEpiEff[[1]]))
} else {
qtlEpiPos <- list(NULL)
qtlEpiNames <- NULL
qtlAllNames <- qtlEachNames
set.seed(seed = seedSelectQTLDom)
actionType <- list(rbinom(sum(nQTLs), 1, propDomi))
names(actionType[[1]]) <- qtlEachNames
actionTypeEpi <- list(NULL)
qtlVarList$qtlEpi <- 0
qtlEpiEff <- list(NULL)
qtlAllEff <- qtlEachEff
}
qtlNamesList <- qtlEachNamesList
qtlEpiNamesList <- list(qtlEpiNames)
qtlAllNamesList <- list(qtlAllNames)
qtlEff <- qtlEachEff
names(qtlPos) <- names(qtlEachPos) <- names(qtlPlePos) <- names(qtlEpiPos) <-
names(qtlEpiNamesList) <- names(qtlAllNamesList) <-
names(actionType) <- names(actionTypeEpi) <-
names(qtlPleEff) <- names(qtlEachEff) <- names(qtlEff) <-
names(qtlEpiEff) <- names(qtlAllEff) <- traitNames
} else {
if (is.null(nQTLs)) {
message('"nQTLs" was not specified. The nQTLs had been set to "10" for all chromosome & traits.')
nQTLs <- matrix(10, nrow = nTraits, ncol = nChr)
}
if (is.vector(nQTLs)) {
if ((length(nQTLs) != 1) && (length(nQTLs) != nChr) &&
(length(nQTLs) != nTraits)) {
stop(paste("length(nQTLs) must be equal to 1 (all chr have the same",
"size) or equal to nChr."))
} else if (length(nQTLs) == 1) {
nQTLs <- matrix(nQTLs, nrow = nTraits, ncol = nChr)
} else if (length(nQTLs) == nTraits) {
nQTLs <- matrix(rep(nQTLs, nChr), nrow = nTraits, ncol = nChr)
} else if (length(nQTLs) == nChr) {
nQTLs <- matrix(rep(nQTLs, nTraits), nrow = nTraits,
ncol = nChr, byrow = TRUE)
}
} else if (is.matrix(nQTLs)) {
if (length(dim(nQTLs)) == 2) {
if (!all(dim(nQTLs) == c(nTraits, nChr))) {
stop(paste0("The dimension of `nQTLs` must be equal to the number of traits ",
"and the number of chromosomes respectively."))
}
} else {
stop("`nQTLs` must be a matrix (2-dimensional array) class object.")
}
} else {
stop("`nQTLs` must be a matrix (2-dimensional array) class object.")
}
if (!is.numeric(nQTLs)) stop("nQTLs must be numeric.")
if (any(floor(nQTLs) - nQTLs != 0)) stop("nQTLs must be integers.")
stopifnot(!any(nQTLs < 0))
rownames(nQTLs) <- traitNames
colnames(nQTLs) <- chrNames
if (is.null(qtlVarList)) {
qtlVarList <- list(qtlPle = rep(1, nTraits),
qtlEach = rep(1, nTraits),
qtlEpi = rep(1, nTraits))
}
if (!all(c("qtlPle", "qtlEach", "qtlEpi") %in% names(qtlVarList))) {
stop("`qtlVarList` should contain objects names `qtlPle`, `qtlEach`, `qtlEpi`.")
} else {
effVarPle <- qtlVarList$qtlPle
effVarEach <- qtlVarList$qtlEach
effVarEpi <- qtlVarList$qtlEpi
if (!is.numeric(effVarPle)) stop("qtlPle must be numeric.")
if (length(effVarPle) != 1 && length(effVarPle) != nTraits) {
stop(paste("length(qtlPle) must be equal to 1 (all trait have the same",
"size) or equal to nTraits."))
} else if (length(effVarPle) == 1) {
effVarPle <- rep(effVarPle, nTraits)
}
if (!is.numeric(effVarEach)) stop("qtlEach must be numeric.")
if (length(effVarEach) != 1 && length(effVarEach) != nTraits) {
stop(paste("length(qtlEach) must be equal to 1 (all trait have the same",
"size) or equal to nTraits."))
} else if (length(effVarEach) == 1) {
effVarEach <- rep(effVarEach, nTraits)
}
if (!is.numeric(effVarEpi)) stop("qtlEpi must be numeric.")
if (length(effVarEpi) != 1 && length(effVarEpi) != nTraits) {
stop(paste("length(qtlEpi) must be equal to 1 (all trait have the same",
"size) or equal to nTraits."))
} else if (length(effVarEpi) == 1) {
effVarEpi <- rep(effVarEpi, nTraits)
}
}
names(effVarPle) <- names(effVarEach) <- names(effVarEpi) <- traitNames
qtlVarList <- list(qtlPle = effVarPle,
qtlEach = effVarEach,
qtlEpi = effVarEpi)
if (qtlOverlap) {
if (is.null(nOverlap)) {
message('"nOverlap" was not specified. The nOverlap had been set to "2" for all chromosome & traits.')
nOverlap <- rep(2, nChr)
}
if (!is.numeric(nOverlap)) stop("nOverlap must be numeric.")
if (any(floor(nOverlap) - nOverlap != 0)) stop("nOverlap must be integers.")
if (length(nOverlap) != 1 && length(nOverlap) != nChr) {
stop(paste("length(nOverlap) must be equal to 1 (all chr have the same",
"size) or equal to nChr."))
} else if (length(nOverlap) == 1) {
nOverlap <- rep(nOverlap, nChr)
}
stopifnot(!any(nOverlap < 0))
if (sum(nOverlap) == 0) {
stop("If nOverlap = 0 for all chromosomes, you should set 'qtlOverlap = FALSE'!!!")
}
if (any(apply(nQTLs, 1, function(x) any(x < nOverlap)))) {
stop("The number of overlapping QTLs across traits (`nOverlap`) should be equal to or less than nQTLs for any chromosome or any trait.")
}
if (is.null(effCor)) {
message('"effCor" was not specified. The effCor had been set to "0.5" for all chromosome & traits.')
effCor <- 0.5
}
if (!is.numeric(effCor)) stop("effCor must be numeric.")
if (is.vector(effCor)) {
if (length(effCor) != 1) {
stop("`effCor` should be scalar or 2-dimensional matrix with `nTraits x nTraits` dimension.")
}
effCorMat <- diag(nTraits)
effCorMat[upper.tri(x = effCorMat)] <- effCor
effCorMat[lower.tri(x = effCorMat)] <- effCor
effCor <- effCorMat
} else if (is.matrix(effCor)) {
if (length(dim(effCor)) != 2) {
stop("`effCor` should be scalar or 2-dimensional matrix with `nTraits x nTraits` dimension.")
} else {
if (!all(dim(effCor) == c(nTraits, nTraits))) {
stop("`effCor` should be scalar or 2-dimensional matrix with `nTraits x nTraits` dimension.")
}
stopifnot(all(diag(effCor) == 1))
}
} else {
stop("`effCor` should be scalar or 2-dimensional matrix with `nTraits x nTraits` dimension.")
}
if (any((effCor > 1) | (effCor < -1))) {
stop("All elements in `effCor` should be within the interval [-1, 1].")
}
if (is.null(qtlPleNamesList)) {
qtlPleNamesList <- sapply(X = rep(sum(nOverlap), nTraits), function(i) {
paste0("QTLPle_", 1:i)
}, simplify = FALSE)
}
stopifnot(length(qtlPleNamesList) == nTraits)
stopifnot(all(unlist(lapply(qtlPleNamesList, length)) == rep(sum(nOverlap), nTraits)))
} else {
nOverlap <- rep(0, nChr)
effCor <- diag(nTraits)
qtlPleNamesList <- rep(list(NULL), nTraits)
}
names(nOverlap) <- chrNames
rownames(effCor) <- colnames(effCor) <- traitNames
nLociUsed <- nMarkers + (apply(nQTLs, 2, sum) - nOverlap)
nQTLsEach <- nQTLs - matrix(rep(nOverlap, nTraits), nrow = nTraits, byrow = TRUE)
if (is.null(qtlEachNamesList)) {
qtlEachNamesList <- sapply(X = apply(nQTLsEach, 1, sum), function(i) {
if (i >= 1) {
return(paste0("QTL_", 1:i))
} else {
return(NULL)
}
}, simplify = FALSE)
}
stopifnot(length(qtlEachNamesList) == nTraits)
stopifnot(all(unlist(lapply(qtlEachNamesList, length)) == apply(nQTLsEach, 1, sum)))
names(qtlPleNamesList) <- names(qtlEachNamesList) <- traitNames
nEpi <- rep(NA, nTraits)
names(nEpi) <- traitNames
qtlPos <- qtlEachPos <- qtlEpiPos <-
qtlPosEachChrList <- qtlEachPosEachChrList <-
qtlNamesList <- qtlEpiNamesList <- qtlAllNamesList <-
actionType <- actionTypeEpi <- qtlEff <-
qtlPleEff <- qtlEachEff <- qtlEpiEff <-
qtlAllEff <- rep(list(NULL), nTraits)
if (qtlOverlap) {
set.seed(seed = seedSelectQTLPle)
qtlPlePosEachChr <- sapply(1:nChr, function(chrNo) {
sampleVec(x = qtlPosCand[[chrNo]], size = nOverlap[chrNo])
}, simplify = FALSE)
names(qtlPlePosEachChr) <- chrNames
qtlPlePosEachChrList <- rep(list(qtlPlePosEachChr), nTraits)
set.seed(seed = seedSelectQTLPle)
qtlPlePos <- rep(list(sampleVec(unlist(sapply(1:nChr, function(chrNo) {
qtlPlePosEachChr[[chrNo]] + cumsum(c(0, nLoci))[chrNo]
}, simplify = FALSE)))), nTraits)
qtlPlePos <- lapply(1:nTraits, function(traitNo) {
names(qtlPlePos[[traitNo]]) <- qtlPleNamesList[[traitNo]]
return(qtlPlePos[[traitNo]])
})
qtlPosCandExcPle <- sapply(1:nChr, function(chrNo) {
qtlPosCandExcPleEachChr <- (qtlPosCand[[chrNo]])[!(names(qtlPosCand[[chrNo]]) %in%
names(qtlPlePosEachChr[[chrNo]]))]
return(qtlPosCandExcPleEachChr)
}, simplify = FALSE)
effVarPleMat <- diag(sqrt(effVarPle)) %*% effCor %*% diag(sqrt(effVarPle))
rownames(effVarPleMat) <- colnames(effVarPleMat) <- traitNames
set.seed(seed = seedSelectQTLPle)
qtlPleEffs0 <- matrix(MASS::mvrnorm(n = sum(nOverlap), mu = rep(0, nTraits),
Sigma = effVarPleMat / sum(nOverlap)),
nrow = sum(nOverlap), ncol = nTraits)
# qtlPleEffs0 <- matrix(mvrnorm(n = 1, mu = rep(0, nTrait * nOverlap),
# Sigma = kronecker(diag(nOverlap), effVarAll)),
# nrow = nOverlap, byrow = TRUE)
qtlPleEffSize <- apply(qtlPleEffs0, 1, function(x) sum(x ^ 2))
qtlPleEffs <- qtlPleEffs0[order(qtlPleEffSize, decreasing = TRUE), , drop = FALSE]
rownames(qtlPleEffs) <- qtlPleNamesList[[1]]
} else {
qtlPlePosEachChr <- rep(list(NULL), nChr)
names(qtlPlePosEachChr) <- chrNames
qtlPlePosEachChrList <- rep(list(qtlPlePosEachChr), nTraits)
qtlPlePos <- rep(list(NULL), nTraits)
qtlPosCandExcPle <- qtlPosCand
qtlPleEffs <- NULL
}
for (traitNo in 1:nTraits) {
qtlPleNames <- qtlPleNamesList[[traitNo]]
qtlEachNames <- qtlEachNamesList[[traitNo]]
qtlNames <- c(qtlPleNames, qtlEachNames)
set.seed(seed = seedSelectQTLEpi[traitNo])
nEpi[traitNo] <- round(sum(rpois(n = sum(nQTLs), lambda = interactionMean[traitNo])) / 2, 0)
if (traitNo == 1) {
qtlPosCandNow <- qtlPosCandExcPle
} else {
qtlPosCandNow <- sapply(1:nChr, function(chrNo) {
qtlPosCandNowEachChr <- (qtlPosCandNow[[chrNo]])[!(names(qtlPosCandNow[[chrNo]]) %in%
names(qtlEachPosEachChr[[chrNo]]))]
return(qtlPosCandNowEachChr)
}, simplify = FALSE)
}
set.seed(seed = seedSelectQTLEach[traitNo])
qtlEachPosEachChr <- sapply(1:nChr, function(chrNo) {
sampleVec(x = qtlPosCandNow[[chrNo]], size = nQTLsEach[traitNo, chrNo])
}, simplify = FALSE)
names(qtlEachPosEachChr) <- chrNames
qtlEachPosEachChrList[[traitNo]] <- qtlEachPosEachChr
set.seed(seed = seedSelectQTLEach[traitNo])
qtlEachPos[[traitNo]] <- sampleVec(unlist(sapply(1:nChr, function(chrNo) {
qtlEachPosEachChr[[chrNo]] + cumsum(c(0, nLoci))[chrNo]
}, simplify = FALSE)))
names(qtlEachPos[[traitNo]]) <- qtlEachNames
qtlPos[[traitNo]] <- c(qtlPlePos[[traitNo]], qtlEachPos[[traitNo]])
names(qtlPos[[traitNo]]) <- c(qtlPleNames, qtlEachNames)
qtlPosEachChr <- sapply(1:nChr, function(chrNo) {
c(qtlPlePosEachChr[[chrNo]], qtlEachPosEachChr[[chrNo]])
}, simplify = FALSE)
qtlPosEachChrList[[traitNo]] <- qtlPosEachChr
if (qtlOverlap) {
qtlPleEffNow <- qtlPleEffs[, traitNo]
names(qtlPleEffNow) <- qtlPleNames
} else {
qtlPleEffNow <- NULL
}
set.seed(seed = seedSelectQTLEach[traitNo])
qtlEachEffNow0 <- rnorm(n = sum(nQTLsEach[traitNo, ]), mean = 0,
sd = sqrt(effVarEach[traitNo] / sum(nQTLsEach[traitNo, ])))
qtlEachEffNow <-
qtlEachEffNow0[order(abs(qtlEachEffNow0), decreasing = TRUE)]
names(qtlEachEffNow) <- qtlEachNames
qtlEffNow <- c(qtlPleEffNow, qtlEachEffNow)
if (nEpi[traitNo] > 0) {
set.seed(seed = seedSelectQTLEpi[traitNo])
qtlEpiPos0 <- t(sapply(1:nEpi[traitNo], function(x) sort(sampleVec(qtlPos[[traitNo]], 2))))
qtlEpiPos0 <- qtlEpiPos0[!duplicated(qtlEpiPos0), , drop = FALSE]
nEpi[traitNo] <- nrow(qtlEpiPos0)
qtlEpiNames <- paste0(qtlNames[match((qtlEpiPos0)[, 1], qtlPos[[traitNo]])], " x ",
qtlNames[match((qtlEpiPos0)[, 2], qtlPos[[traitNo]])])
rownames(qtlEpiPos0) <- qtlEpiNames
colnames(qtlEpiPos0) <- paste0("Epi_", 1:2)
qtlEpiNamesList[[traitNo]] <- qtlEpiNames
qtlEpiPos[[traitNo]] <- qtlEpiPos0
qtlAllNames <- c(qtlNames, qtlEpiNames)
set.seed(seed = seedSelectQTLDom[traitNo])
actionType[[traitNo]] <- rbinom(sum(nQTLs[traitNo, ]), 1, propDomi[traitNo])
names(actionType[[traitNo]]) <- qtlNames
set.seed(seed = seedSelectQTLDom[traitNo])
if (!actionTypeEpiSimple) {
actionTypeEpi[[traitNo]] <- matrix(rbinom(2 * nEpi[traitNo], 1, propDomi[traitNo]),
nrow = nEpi[traitNo])
} else {
actionTypeEpi[[traitNo]] <- matrix((actionType[[traitNo]])[apply(qtlEpiPos0, 2, function(x) {
match(x, qtlPos[[traitNo]])
})], nrow = nEpi[traitNo])
}
rownames(actionTypeEpi[[traitNo]]) <- qtlEpiNames
colnames(actionTypeEpi[[traitNo]]) <- paste0("Epi_", 1:2)
set.seed(seed = seedSelectQTLEpi[traitNo])
qtlEpiEff0 <- rnorm(n = nEpi[traitNo], mean = 0,
sd = sqrt(effVarEpi[traitNo] / nEpi[traitNo]))
qtlEpiEffNow <- qtlEpiEff0[order(abs(qtlEpiEff0), decreasing = TRUE)]
names(qtlEpiEffNow) <- qtlEpiNames
} else {
qtlAllNames <- c(qtlPleNames, qtlEachNames)
set.seed(seed = seedSelectQTLDom)
actionType[[traitNo]] <- rbinom(sum(nQTLs[traitNo, ]), 1, propDomi[traitNo])
names(actionType[[traitNo]]) <- qtlNames
qtlVarList$qtlEpi[traitNo] <- 0
qtlEpiEffNow <- NULL
}
qtlNamesList[[traitNo]] <- qtlNames
qtlAllNamesList[[traitNo]] <- qtlAllNames
qtlAllEffNow <- c(qtlPleEffNow, qtlEachEffNow, qtlEpiEffNow)
if (!is.null(qtlPleEffNow)) {
qtlPleEff[[traitNo]] <- qtlPleEffNow
}
qtlEachEff[[traitNo]] <- qtlEachEffNow
qtlEff[[traitNo]] <- qtlEffNow
if (!is.null(qtlEpiEffNow)) {
qtlEpiEff[[traitNo]] <- qtlEpiEffNow
}
qtlAllEff[[traitNo]] <- qtlAllEffNow
} ### traitNo
names(qtlPos) <- names(qtlEachPos) <- names(qtlPlePos) <-
names(qtlEpiPos) <- names(qtlPosEachChrList) <-
names(qtlEachPosEachChrList) <- names(qtlPlePosEachChrList) <-
names(qtlNamesList) <- names(qtlEpiNamesList) <- names(qtlAllNamesList) <-
names(actionType) <- names(actionTypeEpi) <-
names(qtlPleEff) <- names(qtlEachEff) <- names(qtlEff) <-
names(qtlEpiEff) <- names(qtlAllEff) <- traitNames
}
if (any(lociInfo$specie$nLoci < nLociUsed)) {
stop("The number of loci used (sum of the number of markers & QTLs) must be lower than nLoci for any chromosome.")
}
mrkPosCand <- sapply(1:nChr, function(chrNo) {
qtlPosCandNowEachChr <- (qtlPosCandNow[[chrNo]])[!(names(qtlPosCandNow[[chrNo]]) %in%
names(qtlEachPosEachChr[[chrNo]]))]
}, simplify = FALSE)
} else {
nQTLs <- qtlVarList <- qtlOverlap <- nOverlap <- effCor <- propDomi <-
interactionMean <- qtlPos <- qtlEachPos <- qtlPlePos <- qtlEpiPos <-
qtlPosEachChrList <- qtlEachPosEachChrList <- qtlPlePosEachChrList <-
qtlNamesList <- qtlEachNamesList <- qtlPleNamesList <- qtlEpiNamesList <-
qtlAllNamesList <- actionType <- actionTypeEpi <- qtlPleEff <-
qtlEachEff <- qtlEff <- qtlEpiEff <- qtlAllEff <- seedSelectQTLEach <-
seedSelectQTLPle <- seedSelectQTLDom <- seedSelectQTLEpi <-
nQTLsEach <- nEpi <- NULL
mrkPosCand <- qtlPosCand
}
if (!is.null(seedSelectMarker)) {
if (!is.na(seedSelectMarker)) {
stopifnot(is.numeric(seedSelectMarker))
seedSelectMarker <- floor(seedSelectMarker)
} else {
seedSelectMarker <- sample(x = 1e9, size = 1)
}
}
set.seed(seed = seedSelectMarker)
mrkPosEachChr <- sapply(1:nChr, function(chrNo) {
sort(sampleVec(x = mrkPosCand[[chrNo]], size = nMarkers[chrNo], replace = FALSE))
}, simplify = FALSE)
names(mrkPosEachChr) <- chrNames
mrkPos <- unlist(sapply(1:nChr, function(chrNo) {
mrkPosEachChr[[chrNo]] + cumsum(c(0, nLoci))[chrNo]
}, simplify = FALSE))
self$lociInfo <- lociInfo
self$nMarkers <- nMarkers
self$nTraits <- nTraits
self$traitNames <- traitNames
self$nQTLs <- nQTLs
self$qtlVarList <- qtlVarList
self$qtlOverlap <- qtlOverlap
self$nOverlap <- nOverlap
self$effCor <- effCor
self$propDomi <- propDomi
self$interactionMean <- interactionMean
self$actionTypeEpiSimple <- actionTypeEpiSimple
self$seedSelectMarker <- seedSelectMarker
self$seedSelectQTLEach <- seedSelectQTLEach
self$seedSelectQTLPle <- seedSelectQTLPle
self$seedSelectQTLDom <- seedSelectQTLDom
self$seedSelectQTLEpi <- seedSelectQTLEpi
self$mrkPos <- mrkPos
self$mrkPosEachChr <- mrkPosEachChr
self$qtlPos <- qtlPos
self$qtlEachPos <- qtlEachPos
self$qtlPlePos <- qtlPlePos
self$qtlEpiPos <- qtlEpiPos
self$qtlPosEachChrList <- qtlPosEachChrList
self$qtlEachPosEachChrList <- qtlEachPosEachChrList
self$qtlPlePosEachChrList <- qtlPlePosEachChrList
self$qtlNamesList <- qtlNamesList
self$qtlEachNamesList <- qtlEachNamesList
self$qtlPleNamesList <- qtlPleNamesList
self$qtlEpiNamesList <- qtlEpiNamesList
self$qtlAllNamesList <- qtlAllNamesList
self$actionType <- actionType
self$actionTypeEpi <- actionTypeEpi
self$qtlPleEff <- qtlPleEff
self$qtlEachEff <- qtlEachEff
self$qtlEff <- qtlEff
self$qtlEpiEff <- qtlEpiEff
self$qtlAllEff <- qtlAllEff
self$nQTLsEach <- nQTLsEach
self$nEpi <- nEpi
},
#' @description Detail information of map for marker genome
#' @param lociNames [character] loci ids
#' @examples
#' myLoci$genoMapDetail()
#' myLoci$genoMapDetail("Loci_6")
genoMapDetail = function(lociNames = NULL) {
genoMap <- self$lociInfo$genoMap
if (self$lociInfo$specie$simInfo$simPheno) {
for (traitNo in 1:self$nTraits) {
traitNameNow <- self$traitNames[traitNo]
qtlNamesNow <- self$qtlNamesList[[traitNo]]
qtlPosNow <- self$qtlPos[[traitNo]]
qtlNameForMapNow <- rep(NA, sum(self$lociInfo$specie$nLoci))
qtlNameForMapNow[qtlPosNow] <- qtlNamesNow
genoMap <- data.frame(genoMap, qtlNameForMapNow)
colnames(genoMap)[ncol(genoMap)] <- traitNameNow
}
} else {
genoMapNA <- matrix(NA, nrow = nrow(genoMap), ncol = self$nTraits)
rownames(genoMapNA) <- rownames(genoMap)
colnames(genoMapNA) <- self$traitNames
genoMap <- cbind(genoMap, genoMapNA)
}
if (is.null(lociNames)) {
return(genoMap)
} else {
genoMap[genoMap$lociNames %in% lociNames]
}
},
#' @description Display summary information about the object: traitInfo
print = function() {
cat(paste0("Number of traits: ", self$nTraits, "\n",
"Number of markers: \n"))
print(self$nMarkers)
cat(paste0("Number of QTLs (except epistasis) in total: \n"))
print(self$nQTLs)
cat(paste0("----------------------------------------- \n"))
cat(paste0("Number of QTLs specific to each trait: \n"))
print(self$nQTLsEach)
cat(paste0("Number of pleiotropic QTLs across traits: \n"))
print(self$nOverlap)
cat(paste0("Number of epistatic effects for each trait: \n"))
print(self$nEpi)
},
#' @description
#' plot QTL map using the \pkg{plotly} package
#' @param alphaMarker [numeric] transparency for markers see \link[plotly]{plot_ly}
#' @param alphaQTL [numeric] transparency for QTLs see \link[plotly]{plot_ly}
#' @param sizeMarker [numeric] size for markers
#' @param sizeQTLmin [numeric] size for QTL with minimum effect
#' @param sizeQTLmax [numeric] size for QTL with maximum effect
#'
#' @import plotly
#'
#' @examples
#' myLoci$plot(alphaMarker = 0.1, alphaQTL = 0.9,
#' sizeMarker = 5, sizeQTLmin = 6, sizeQTLmax = 12)
plot = function(alphaMarker = 0.05,
alphaQTL = 0.8,
sizeMarker = 5,
sizeQTLmin = 6,
sizeQTLmax = 12) {
ends <- self$lociInfo$specie$lChr
chrNames <- self$lociInfo$specie$chrNames
if (self$lociInfo$specie$simInfo$simPheno) {
genoDetail <- self$genoMapDetail()
genoDetail$rec <- round(genoDetail$rec, 3)
genoDetailQTLs <- genoDetail[apply(genoDetail[, -c(1:4), drop = FALSE], 1, function(x) any(!is.na(x))), ]
genoDetailMarkers <- genoDetail[self$mrkPos, ]
effQTLs <- sapply(X = 1:self$nTraits, function(traitNo) {
effQTLNow <- rep(NA, nrow(genoDetailQTLs))
effQTLNow[match(names(self$qtlEff[[traitNo]]), genoDetailQTLs[, traitNo + 4])] <-
self$qtlEff[[traitNo]]
return(effQTLNow)
})
effQTLsWithNA <- effQTLs
effQTLs[is.na(effQTLs)] <- 0
rownames(effQTLs) <- genoDetailQTLs$lociNames
colnames(effQTLs) <- self$traitNames
effQTLSize <- apply(effQTLs, 1, function(x) sqrt(sum(x ^ 2)))
effQTLSizeScaled <- sizeQTLmin + (effQTLSize - min(effQTLSize)) *
(sizeQTLmax - sizeQTLmin) / diff(range(effQTLSize))
for (traitNo in 1:self$nTraits) {
nonNANow <- !is.na(genoDetailQTLs[, traitNo + 4])
actionTypeNow <- self$actionType[[traitNo]]
actionTypeNow <- factor(ifelse(actionTypeNow == 1, "Dom", "Add"), levels = c("Add", "Dom"))
genoDetailQTLs[nonNANow, traitNo + 4] <-
paste0(genoDetailQTLs[nonNANow, traitNo + 4], ", ",
round(effQTLsWithNA[nonNANow, traitNo], 3), ", ",
actionTypeNow[genoDetailQTLs[nonNANow, traitNo + 4]])
}
} else {
genoDetailMarkers <- self$genoMapDetail()[self$mrkPos, ]
}
genoDetailMarkers$chr <- factor(genoDetailMarkers$chr, levels = chrNames)
plt <- plotly::plot_ly(data = genoDetailMarkers,
x = ~ chr,
y = ~ pos,
type = "scatter",
mode = "markers",
alpha = alphaMarker,
marker = list(size = sizeMarker),
name = "Markers",
hoverinfo = 'text',
text = apply(genoDetailMarkers, 1, function(l) {
paste(names(l), ":", l, collapse = "\n")
})) %>%
plotly::add_markers(x = rep(names(ends), 2),
y = c(ends, rep(0, length(ends))),
alpha = 1,
name = "Chromosome's edges",
hoverinfo = 'text',
text = paste(rep(names(ends), 2),
": length =",
c(ends, rep(0, length(ends)))))
if (self$lociInfo$specie$simInfo$simPheno) {
genoDetailQTLs$chr <- factor(genoDetailQTLs$chr, levels = chrNames)
plt <- plt %>%
plotly::add_markers(data = genoDetailQTLs,
x = ~ chr,
y = ~ pos,
type = "scatter",
mode = "markers",
alpha = alphaQTL,
marker = list(size = effQTLSizeScaled),
name = "QTLs",
hoverinfo = 'text',
text = apply(genoDetailQTLs, 1, function(l) {
paste(names(l), ":", l, collapse = "\n")
}))
}
print(plt)
}
)
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.