Nothing
R/new-restrict.R
In OncoSimulR: Forward Genetic Simulation of Cancer Progression with Epistasis
Defines functions
list_g_matches_fixed
print.sampledGenotypes
sampledGenotypes
detectionProbCheckParse
matchGeneIDs
geneCounts
get.gene.counts
allNamedGenes
countGenesFe
nr_oncoSimul.internal
plot.fitnessEffects
fitnessEffectsToIgraph
evalAllGenotypesFitAndMut
generateAllGenotypes
evalAllGenotypes
evalAllGenotypesORMut
evalGenotypeFitAndMut
evalGenotype
evalGenotypeORMut
allFitnessEffects
allFitnessORMutatorEffects
getGeneIDNum
getNamesID
checkRT
addIntID.epist.order
to.long.epist.order
to.pairs.modules
epist.order.to.pairs.modules
oe.to.df
epist.order.element
to.long.rt
list.of.deps
geneModuleNull
gm.to.geneModuleL
nice.vector.eo
gtm2
check.gm
Documented in
allFitnessEffects
evalAllGenotypes
evalAllGenotypesFitAndMut
evalGenotype
evalGenotypeFitAndMut
plot.fitnessEffects
print.sampledGenotypes
sampledGenotypes
## Copyright 2013, 2014, 2015, 2016 Ramon Diaz-Uriarte
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
## 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 for more details.
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.
## Say which are drivers: populate the drv vector.
## // In R, the user says which are the drivers. If does not say anuthing,
## // the default (NULL) then drivers are all in poset, epist, restrict. The
## // user can pass a vector with the names of the genes (not modules). Allow
## // also for empty, so this is faster if not needed. And check that if we
## // use a stopping rule based on drivers that drv vectors is not empty.
## - Say that a user can use a "0" as a gene name, but that is BAD idea.
## - Modules and order effects can be kind of funny?
## Gene names can contain no spaces, commas, or ">"
check.gm <- function(gm) {
## Yes, Root: we want no ambiguities.
## Actually, that sucks. So we do not require it, but check for
## consistency.
if(any(gm == "Root") && (gm[1] != "Root") )
stop("If Root is in the module table, it must be the first")
if(any(names(gm) == "Root") && (names(gm)[1] != "Root") )
stop("If the name Root is in the module table, it must be the first")
if( (names(gm)[1] == "Root") && (gm[1] != "Root") )
stop("The name Root is in the module table, but is not of Root")
if( (gm[1] == "Root") && (names(gm)[1] != "Root") )
stop("Root is in the module table, but with a different name")
## if(gm[1] != "Root")
## stop("First value of a module table must be Root")
## if(names(gm)[1] != "Root")
## stop("First name of a module table must be Root")
if(length(unique(names(gm))) != length(gm))
stop("Number of unique module names different from length of vector")
if(gm[1] != "Root")
gm <- c("Root" = "Root", gm)
return(gm)
}
gtm2 <- function(x) {
data.frame(cbind(nice.vector.eo(x, ","), x), stringsAsFactors = TRUE)
}
## nice.vector.eo <- function(z, sep) {
## ## with epistasis, maybe we want sorted?
## setdiff(unlist(lapply(strsplit(z, " "),
## function(u) strsplit(u, sep))), "")
## }
nice.vector.eo <- function(z, sep, rm.sign = FALSE) {
## with epistasis, maybe we want sorted?
if(! rm.sign)
return(setdiff(unlist(lapply(strsplit(z, " "),
function(u) strsplit(u, sep))), ""))
else ## remove the " ", the -, and the sep
return(setdiff(unlist(lapply(strsplit(z, " "), function(u)
strsplit(unlist(strsplit(u, "-")), sep))), ""))
}
gm.to.geneModuleL <- function(gm, one.to.one) {
## the table will be sorted by gene name
gm <- check.gm(gm)
## the named vector with the mapping into the long geneModule df
geneMod <- as.data.frame(rbindlist(lapply(gm, gtm2)),
stringsAsFactors = TRUE) ## this stringsAsFactors is key
geneMod$Module <- names(gm)[geneMod[, 2]] ## reverse lookup table
colnames(geneMod)[1] <- c("Gene")
geneMod <- geneMod[, -2]
geneMod$Gene <- as.character(geneMod$Gene)
## geneMod$Module <- as.character(geneMod$Module) ## already a char
geneMod <- geneMod[c(1, order(geneMod$Gene[-1]) + 1), ]
geneMod$GeneNumID <- 0:(nrow(geneMod) - 1)
## this assumes sorted! and they need not be
## rlem <- rle(geneMod$Module)
## geneMod$ModuleNumID <- rep( 0:(length(rlem$values) - 1), rlem$lengths)
if(one.to.one) {
geneMod$ModuleNumID <- geneMod$GeneNumID
if(!(identical(geneMod$Gene, geneMod$Module)))
stop("Impossible: we are supposed to be in one-to-one for Module-Gene.")
} else {
## Works, but does not give the most ordered module names. But
## keeps implicit order given by user.
idm <- seq.int(length(gm) - 1)
idm <- c("Root" = 0L, idm)
names(idm) <- names(gm)
## This sorts by the names but is not optimal either
## idm1 <- seq.int(length(gm) - 1)
## idm <- c(0L, idm1)
## names(idm) <- c("Root", sort(setdiff(names(gm), "Root")))
geneMod$ModuleNumID <- idm[geneMod[, "Module"]]
}
if(length(unique(geneMod$Gene)) != nrow(geneMod)) {
stop("Are there identical gene names in different modules?")
}
## I think this is unreacheable now. Caught earlier.
if(length(unique(geneMod$GeneNumID)) != nrow(geneMod)) {
stop("Are there identical gene names in different modules?")
}
rownames(geneMod) <- 1:nrow(geneMod)
geneMod
}
geneModuleNull <- function(namesM) {
v <- c("Root", setdiff(namesM, "Root"))
names(v) <- v
return(v)
}
list.of.deps <- function(x) {
## lookupTypeDep <- c("MN" = 1, "monotone" = 1,
## "SM" = 2, "semimonotone" = 2)
lookupTypeDep <- c("MN" = "monotone",
"AND" = "monotone",
"monotone" = "monotone",
"CMPN" = "monotone",
"OR" = "semimonotone",
"SM" = "semimonotone",
"semimonotone" = "semimonotone",
"DMPN" = "semimonotone",
"XOR" = "xmpn",
"xmpn" = "xmpn",
"XMPN" = "xmpn",
"--" = "--",
"-" = "--")
## FIXME: check values of typeDep
if(length(x) > 1) {
if(length(unique(x$s))!= 1)
stop("Not all s identical within a child")
if(length(unique(x$sh))!= 1)
stop("Not all sh identical within a child")
if(length(unique(x$typeDep))!= 1)
stop("Not all typeDep identical within a child")
if(length(unique(x$child))!= 1)
stop("child not unique")
}
typeDep <- lookupTypeDep[unique(x$typeDep)]
if(any(is_null_na(typeDep)))
stop("typeDep value incorrect. Check spelling.")
return(list(
child = unique(x$child),
s = unique(x$s),
sh = unique(x$sh),
typeDep = typeDep,
parents = unlist(x$parent)))
}
to.long.rt <- function(rt, idm) {
## We now do this inconditionally, so that we do not need to use the
## "stringsAsFactors = FALSE". This is now done before
## if(is.numeric(rt$parent))
## rt$parent <- as.character(rt$parent)
## if(is.numeric(rt$child))
## rt$child <- as.character(rt$child)
if(!("Root" %in% rt$parent))
stop("Root must be one parent node")
## rt$parent <- unlist(lapply(rt$parent, nice.string))
## rt$child <- unlist(lapply(rt$child, nice.string))
srt <- rt[order(rt$child), ]
## Not relevant if we allow non-numeric names
## all.child.genes <- as.integer(
## unlist(lapply(rt[, 2],
## function(x) strsplit(x, ","))))
## ## check all childs
## if(!identical(sort(unique(all.child.genes)),
## seq.int(max(all.child.genes))))
## stop("Not all children present")
long.rt <- lapply(split(srt, srt$child), list.of.deps)
## geneModule <- gene.to.module(srt)
## idm <- seq.int(length(names(long.rt)))
## names(idm) <- names(long.rt)
## idm <- c("0" = 0L, idm)
## geneModule$ModuleNumID <- idm[geneModule[, "Module"]]
## idm is just a look up table for the id of the module
## idm <- unique(geneModule$ModuleNumID)
## names(idm) <- unique(geneModule$Module)
## add integer IDs
addIntID <- function(z, idm) {
z$childNumID <- idm[z$child]
z$parentsNumID <- idm[z$parents]
if( any(is.na(z$parentsNumID)) ||
any(is.na(z$childNumID)) ) {
## I think this can no longer be reached ever. Caught before.
stop(paste("An ID is NA:",
"Is a gene part of two different modules?",
"(That includes being by itself and part",
"of a module.)"))
}
## These checks could be somewhere else instead of here.
if(length(unique(z$parentsNumID)) != length(z$parentsNumID))
stop("No node can have the same parent multiple times")
if(length(unique(z$parents)) != length(z$parents))
stop("No node can have the same parent multiple times")
if(length(z$child) > 1)
stop("Child nodes can have one, and only one, member")
## sort the parents, always.
o <- order(z$parentsNumID)
z$parentsNumID <- z$parentsNumID[o]
z$parents <- as.character(z$parents[o])
return(z)
}
long.rt <- lapply(long.rt, function(x) addIntID(x, idm = idm))
## if(verbosity >= 4) {
## message(paste("Number of drivers: ",
## length(unique(geneModule[, "Gene"]))))
## message(paste("Number of modules: ",
## length(unique(geneModule[, "Module"]))))
## }
return(long.rt)
## return(list(long.rt = long.rt, geneModule = geneModule))
}
epist.order.element <- function(x, y, sep, rm.sign = FALSE) {
list(ids = nice.vector.eo(x, sep = sep, rm.sign = rm.sign), s = y)
}
oe.to.df <- function(x) {
ma <- matrix(ncol = 2, nrow = 1 + length(x) - 2)
if(length(x) == 2) {
ma[1, 1] <- x[1]
ma[1, 2] <- x[2]
} else {
ma[, 1] <- x[-length(x)]
ma[, 2] <- x[-1]
}
return(data.frame(ma, stringsAsFactors = FALSE))
}
epist.order.to.pairs.modules <- function(x, sep, rm.sign = TRUE) {
## We discard, do not even accept, the coefficient
tmp <- epist.order.element(x, -99, sep = sep, rm.sign = rm.sign)$ids
if(length(tmp) > 1) {
## if a single gene, as when we specify all genotypes, we do not
## want this
if(sep == ":")
return(data.frame(combinations(n = length(tmp), r = 2, v = tmp),
stringsAsFactors = FALSE))
else if(sep == ">") {
return(oe.to.df(tmp))
}
}
}
to.pairs.modules <- function(x, sep, rm.sign = TRUE) {
df <- data.frame(rbindlist(
lapply(names(x),
function(z) epist.order.to.pairs.modules(z, sep))),
stringsAsFactors = TRUE)
if(nrow(df) == 0L) { ## if only single genes in epist, we get nothing here.
return(df)
} else {
colnames(df) <- c("parent", "child")
if(sep == ":")
df$typeDep <- "epistasis"
else if(sep == ">")
df$typeDep <- "orderEffect"
return(unique(df))
}
}
to.long.epist.order <- function(epor, sep, rm.sign = FALSE) {
## just vectors for now
long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
names(epor), epor)
## if(is.vector(epor))
## long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
## names(epor), epor)
## else if(is.data.frame(epor))
## long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
## as.character(epor$ids),
## epor$s)
names(long) <- NULL
return(long)
}
## addIntID.epist.order <- function(z, idm, sort) {
## z$NumID <- idm[z$ids]
## if(sort) {
## ## essential for epistasis, but never do for order effects
## o <- order(z$NumID)
## z$NumID <- z$NumID[o]
## z$ids <- z$ids[o]
## }
## return(z)
## }
addIntID.epist.order <- function(z, idm, sort, sign) {
if( sort && (!sign))
warning("sort is TRUE but sign is FALSE. You sure?")
if((!sort) && sign)
warning("sort is FALSE but sign is TRUE. You sure?")
## Adds the integer, but takes care of sign if sign = TRUE
if(!sign) {
z$NumID <- idm[z$ids]
signs <- grep("^-", z$ids)
if(length(signs))
stop("You have a negative sign and you said sign = FALSE")
} else {
unsigned <- setdiff(unlist(lapply(z$ids,
function(z) strsplit(z, "^-"))), "")
NumID <- idm[unsigned]
signs <- grep("^-", z$ids)
if(length(signs)) {
NumID[signs] <- -1 * NumID[signs]
}
z$NumID <- NumID
}
if(sort) {
## Essential for epistasis, but never do for order effects
o <- order(z$NumID)
z$NumID <- z$NumID[o]
z$ids <- z$ids[o]
}
if(length(unique(z$NumID)) != length(z$NumID))
stop("No node can have the same id multiple times")
if(length(unique(z$ids)) != length(z$ids))
stop("No node can have the same id multiple times")
return(z)
}
checkRT <- function(mdeps) {
if(ncol(mdeps) != 5)
stop("mdeps must be of exactly 5 columns")
if(!identical(colnames(mdeps), c("parent", "child", "s", "sh", "typeDep")))
stop(paste("Column names of mdeps not of appropriate format. ",
"Should be parent, child, s, sh, typeDep"))
}
getNamesID <- function(fp) {
## Return a lookup table for names based on numeric IDs
idname <- c(fp$geneModule$GeneNumID,
fp$long.geneNoInt$GeneNumID,
fp$fitnessLandscape_gene_id$GeneNumID)
names(idname) <- c(fp$geneModule$Gene,
fp$long.geneNoInt$Gene,
fp$fitnessLandscape_gene_id$Gene)
return(idname[-1]) ## remove Root!!
}
getGeneIDNum <- function(geneModule, geneNoInt, fitnessLandscape_gene_id,
drv, sort = TRUE) {
## It returns the genes, as NumID, in the given vector with names drv
## initMutant uses this, for simplicity, without sorting, but noInt
## are always sorted
## Also used for the drivers with sort = TRUE
## Yes, we must do it twice because we do not know before hand which
## is which. This makes sure no NA. Period.
if(any(is.na( match(drv, c(geneModule$Gene, geneNoInt$Gene,
fitnessLandscape_gene_id$Gene))))) {
stop(paste("For driver or initMutant you have passed genes",
"not in the fitness table."))
}
indicesM <- as.vector(na.omit(match( drv, geneModule$Gene)))
indicesI <- as.vector(na.omit(sort(match( drv, geneNoInt$Gene))))
indicesF <- as.vector(na.omit(sort(match( drv, fitnessLandscape_gene_id$Gene))))
if(sort) {
indicesM <- sort(indicesM)
}
return(c(
geneModule$GeneNumID[indicesM],
geneNoInt$GeneNumID[indicesI],
fitnessLandscape_gene_id$GeneNumID[indicesF])
)
}
allFitnessORMutatorEffects <- function(rT = NULL,
epistasis = NULL,
orderEffects = NULL,
noIntGenes = NULL,
geneToModule = NULL,
drvNames = NULL,
keepInput = TRUE,
genotFitness = NULL,
## refFE = NULL,
calledBy = NULL) {
## From allFitnessEffects. Generalized so we deal with Fitness
## and mutator.
## restrictions: the usual rt
## epistasis: as it says, with the ":"
## orderEffects: the ">"
## All of the above can be genes or can be modules (if you pass a
## geneToModule)
## rest: rest of genes, with fitness
## For epistasis and order effects we create the output object but
## missing the numeric ids of genes. With rT we do it in one go, as we
## already know the mapping of genes to numeric ids. We could do the
## same in epistasis and order, but we would be splitting twice
## (whereas for rT extracting the names is very simple).
## called appropriately?
if( !(calledBy %in% c("allFitnessEffects", "allMutatorEffects") ))
stop("How did you call this function?. Bug.")
if(calledBy == "allMutatorEffects") {
## very paranoid check
if( !is.null(rT) || !is.null(orderEffects) ||
!is.null(drvNames) || !is.null(genotFitness))
stop("allMutatorEffects called with forbidden arguments.",
"Is this an attempt to subvert the function?")
}
rtNames <- NULL
epiNames <- NULL
orNames <- NULL
if(!is.null(rT)) {
## This is really ugly, but to prevent the stringsAsFactors I need it here:
rT$parent <- as.character(rT$parent)
rT$child <- as.character(rT$child)
rT$typeDep <- as.character(rT$typeDep)
rtNames <- unique(c(rT$parent, rT$child))
}
if(!is.null(epistasis)) {
long.epistasis <- to.long.epist.order(epistasis, ":")
## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
## deal with the possible negative signs
epiNames <- setdiff(unique(
unlist(lapply(long.epistasis,
function(x) lapply(x$ids,
function(z) strsplit(z, "^-"))))),
"")
} else {
long.epistasis <- list()
}
if(!is.null(orderEffects)) {
long.orderEffects <- to.long.epist.order(orderEffects, ">")
orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
} else {
long.orderEffects <- list()
}
allModuleNames <- unique(c(rtNames, epiNames, orNames))
if(is.null(geneToModule)) {
gMOneToOne <- TRUE
geneToModule <- geneModuleNull(allModuleNames)
} else {
gMOneToOne <- FALSE
if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
stop(paste("Some values in geneToModule not present in any of",
" rT, epistasis, or order effects"))
if(any(is.na(match(allModuleNames, names(geneToModule)))))
stop(paste("Some values in rT, epistasis, ",
"or order effects not in geneToModule"))
}
geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
idm <- unique(geneModule$ModuleNumID)
names(idm) <- unique(geneModule$Module)
if(!is.null(rT)) {
checkRT(rT)
long.rt <- to.long.rt(rT, idm)
} else {
long.rt <- list() ## yes, we want an object of length 0
}
## Append the numeric ids to epistasis and order
if(!is.null(epistasis)) {
long.epistasis <- lapply(long.epistasis,
function(x)
addIntID.epist.order(x, idm,
sort = TRUE,
sign = TRUE))
}
if(!is.null(orderEffects)) {
long.orderEffects <- lapply(long.orderEffects,
function(x)
addIntID.epist.order(x, idm,
sort = FALSE,
sign = FALSE))
}
if(!is.null(noIntGenes)) {
if(inherits(noIntGenes, "character")) {
wm <- paste("noIntGenes is a character vector.",
"This is probably not what you want, and will",
"likely result in an error downstream.",
"You can get messages like",
" 'not compatible with requested type', and others.",
"We are stopping.")
stop(wm)
}
mg <- max(geneModule[, "GeneNumID"])
gnum <- seq_along(noIntGenes) + mg
if(!is.null(names(noIntGenes))) {
ng <- names(noIntGenes)
if( any(grepl(",", ng, fixed = TRUE)) ||
any(grepl(">", ng, fixed = TRUE)) ||
any(grepl(":", ng, fixed = TRUE)) )
stop("The name of some noIntGenes contain a ',' or a '>' or a ':'")
if(any(ng %in% geneModule[, "Gene"] ))
stop("A gene in noIntGenes also present in the other terms")
if(any(duplicated(ng)))
stop("Duplicated gene names in geneNoInt")
if(any(is.na(ng)))
stop("In noIntGenes some genes have names, some don't.",
" Name all of them, or name none of them.")
} else {
ng <- gnum
}
geneNoInt <- data.frame(Gene = as.character(ng),
GeneNumID = gnum,
s = noIntGenes,
stringsAsFactors = FALSE)
} else {
geneNoInt <- data.frame()
}
if(is.null(genotFitness)) {
genotFitness <- matrix(NA, nrow = 0, ncol = 1)
fitnessLandscape_df <- data.frame()
fitnessLandscape_gene_id <- data.frame()
} else {
## Yes, I am duplicating stuff for now.
## This makes life simpler in C++:
## In the map, the key is the genotype name, as
## cnn <- colnames(genotFitness)[-ncol(genotFitness)]
cnn <- 1:(ncol(genotFitness) - 1)
gfn <- apply(genotFitness[, -ncol(genotFitness), drop = FALSE], 1,
function(x) paste(cnn[as.logical(x)],
collapse = ", "))
## rownames(genotFitness) <- gfn
fitnessLandscape_df <-
data.frame(Genotype = gfn,
Fitness = genotFitness[, ncol(genotFitness)],
stringsAsFactors = FALSE)
fitnessLandscape_gene_id <- data.frame(
Gene = colnames(genotFitness)[-ncol(genotFitness)],
GeneNumID = cnn,
stringsAsFactors = FALSE)
}
if( (length(long.rt) + length(long.epistasis) + length(long.orderEffects) +
nrow(geneNoInt) + nrow(genotFitness)) == 0)
stop("You have specified nothing!")
if(calledBy == "allFitnessEffects") {
if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
} else {
graphE <- NULL
}
} else {
graphE <- NULL
}
if(!is.null(drvNames)) {
drv <- unique(getGeneIDNum(geneModule, geneNoInt, fitnessLandscape_gene_id,
drvNames))
## drivers should never be in the geneNoInt; Why!!!???
## Catch the problem. This is an overkill,
## so since we catch the issue, we could leave the geneNoInt. But
## that should not be there in this call.
## if(any(drvNames %in% geneNoInt$Gene)) {
## stop(paste("At least one gene in drvNames is a geneNoInt gene.",
## "That is not allowed.",
## "If that gene is a driver, pass it as gene in the epistasis",
## "component."))
## }
## drv <- getGeneIDNum(geneModule, NULL, drvNames)
} else {
## we used to have this default
## drv <- geneModule$GeneNumID[-1]
drv <- vector(mode = "integer", length = 0L)
}
if(!keepInput) {
rT <- epistasis <- orderEffects <- noIntGenes <- NULL
}
out <- list(long.rt = long.rt,
long.epistasis = long.epistasis,
long.orderEffects = long.orderEffects,
long.geneNoInt = geneNoInt,
geneModule = geneModule,
gMOneToOne = gMOneToOne,
geneToModule = geneToModule,
graph = graphE,
drv = drv,
rT = rT,
epistasis = epistasis,
orderEffects = orderEffects,
noIntGenes = noIntGenes,
fitnessLandscape = genotFitness,
fitnessLandscape_df = fitnessLandscape_df,
fitnessLandscape_gene_id = fitnessLandscape_gene_id
)
if(calledBy == "allFitnessEffects") {
class(out) <- c("fitnessEffects")
} else if(calledBy == "allMutatorEffects") {
class(out) <- c("mutatorEffects")
}
return(out)
}
## Former version, with fitness landscape
## allFitnessORMutatorEffects <- function(rT = NULL,
## epistasis = NULL,
## orderEffects = NULL,
## noIntGenes = NULL,
## geneToModule = NULL,
## drvNames = NULL,
## keepInput = TRUE,
## ## refFE = NULL,
## calledBy = NULL) {
## ## From allFitnessEffects. Generalized so we deal with Fitness
## ## and mutator.
## ## restrictions: the usual rt
## ## epistasis: as it says, with the ":"
## ## orderEffects: the ">"
## ## All of the above can be genes or can be modules (if you pass a
## ## geneToModule)
## ## rest: rest of genes, with fitness
## ## For epistasis and order effects we create the output object but
## ## missing the numeric ids of genes. With rT we do it in one go, as we
## ## already know the mapping of genes to numeric ids. We could do the
## ## same in epistasis and order, but we would be splitting twice
## ## (whereas for rT extracting the names is very simple).
## ## called appropriately?
## if( !(calledBy %in% c("allFitnessEffects", "allMutatorEffects") ))
## stop("How did you call this function?. Bug.")
## if(calledBy == "allMutatorEffects") {
## ## very paranoid check
## if( !is.null(rT) || !is.null(orderEffects) || !is.null(drvNames))
## stop("allMutatorEffects called with forbidden arguments.",
## "Is this an attempt to subvert the function?")
## }
## rtNames <- NULL
## epiNames <- NULL
## orNames <- NULL
## if(!is.null(rT)) {
## ## This is really ugly, but to prevent the stringsAsFactors I need it here:
## rT$parent <- as.character(rT$parent)
## rT$child <- as.character(rT$child)
## rT$typeDep <- as.character(rT$typeDep)
## rtNames <- unique(c(rT$parent, rT$child))
## }
## if(!is.null(epistasis)) {
## long.epistasis <- to.long.epist.order(epistasis, ":")
## ## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
## ## deal with the possible negative signs
## epiNames <- setdiff(unique(
## unlist(lapply(long.epistasis,
## function(x) lapply(x$ids,
## function(z) strsplit(z, "^-"))))),
## "")
## } else {
## long.epistasis <- list()
## }
## if(!is.null(orderEffects)) {
## long.orderEffects <- to.long.epist.order(orderEffects, ">")
## orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
## } else {
## long.orderEffects <- list()
## }
## allModuleNames <- unique(c(rtNames, epiNames, orNames))
## if(is.null(geneToModule)) {
## gMOneToOne <- TRUE
## geneToModule <- geneModuleNull(allModuleNames)
## } else {
## gMOneToOne <- FALSE
## if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
## stop(paste("Some values in geneToModule not present in any of",
## " rT, epistasis, or order effects"))
## if(any(is.na(match(allModuleNames, names(geneToModule)))))
## stop(paste("Some values in rT, epistasis, ",
## "or order effects not in geneToModule"))
## }
## geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
## idm <- unique(geneModule$ModuleNumID)
## names(idm) <- unique(geneModule$Module)
## if(!is.null(rT)) {
## checkRT(rT)
## long.rt <- to.long.rt(rT, idm)
## } else {
## long.rt <- list() ## yes, we want an object of length 0
## }
## ## Append the numeric ids to epistasis and order
## if(!is.null(epistasis)) {
## long.epistasis <- lapply(long.epistasis,
## function(x)
## addIntID.epist.order(x, idm,
## sort = TRUE,
## sign = TRUE))
## }
## if(!is.null(orderEffects)) {
## long.orderEffects <- lapply(long.orderEffects,
## function(x)
## addIntID.epist.order(x, idm,
## sort = FALSE,
## sign = FALSE))
## }
## if(!is.null(noIntGenes)) {
## if(inherits(noIntGenes, "character")) {
## wm <- paste("noIntGenes is a character vector.",
## "This is probably not what you want, and will",
## "likely result in an error downstream.",
## "You can get messages like",
## " 'not compatible with requested type', and others.",
## "We are stopping.")
## stop(wm)
## }
## mg <- max(geneModule[, "GeneNumID"])
## gnum <- seq_along(noIntGenes) + mg
## if(!is.null(names(noIntGenes))) {
## ng <- names(noIntGenes)
## if( grepl(",", ng, fixed = TRUE) || grepl(">", ng, fixed = TRUE)
## || grepl(":", ng, fixed = TRUE))
## stop("The name of some noIntGenes contain a ',' or a '>' or a ':'")
## if(any(ng %in% geneModule[, "Gene"] ))
## stop("A gene in noIntGenes also present in the other terms")
## if(any(duplicated(ng)))
## stop("Duplicated gene names in geneNoInt")
## if(any(is.na(ng)))
## stop("In noIntGenes some genes have names, some don't.",
## " Name all of them, or name none of them.")
## } else {
## ng <- gnum
## }
## geneNoInt <- data.frame(Gene = as.character(ng),
## GeneNumID = gnum,
## s = noIntGenes,
## stringsAsFactors = FALSE)
## } else {
## geneNoInt <- data.frame()
## }
## if( (length(long.rt) + length(long.epistasis) + length(long.orderEffects) +
## nrow(geneNoInt)) == 0)
## stop("You have specified nothing!")
## if(calledBy == "allFitnessEffects") {
## if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
## graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
## } else {
## graphE <- NULL
## }
## } else {
## graphE <- NULL
## }
## if(!is.null(drvNames)) {
## drv <- unique(getGeneIDNum(geneModule, geneNoInt, drvNames))
## ## drivers should never be in the geneNoInt; Why!!!???
## ## Catch the problem. This is an overkill,
## ## so since we catch the issue, we could leave the geneNoInt. But
## ## that should not be there in this call.
## ## if(any(drvNames %in% geneNoInt$Gene)) {
## ## stop(paste("At least one gene in drvNames is a geneNoInt gene.",
## ## "That is not allowed.",
## ## "If that gene is a driver, pass it as gene in the epistasis",
## ## "component."))
## ## }
## ## drv <- getGeneIDNum(geneModule, NULL, drvNames)
## } else {
## ## we used to have this default
## ## drv <- geneModule$GeneNumID[-1]
## drv <- vector(mode = "integer", length = 0L)
## }
## if(!keepInput) {
## rT <- epistasis <- orderEffects <- noIntGenes <- NULL
## }
## out <- list(long.rt = long.rt,
## long.epistasis = long.epistasis,
## long.orderEffects = long.orderEffects,
## long.geneNoInt = geneNoInt,
## geneModule = geneModule,
## gMOneToOne = gMOneToOne,
## geneToModule = geneToModule,
## graph = graphE,
## drv = drv,
## rT = rT,
## epistasis = epistasis,
## orderEffects = orderEffects,
## noIntGenes = noIntGenes
## )
## if(calledBy == "allFitnessEffects") {
## class(out) <- c("fitnessEffects")
## } else if(calledBy == "allMutatorEffects") {
## class(out) <- c("mutatorEffects")
## }
## return(out)
## }
allFitnessEffects <- function(rT = NULL,
epistasis = NULL,
orderEffects = NULL,
noIntGenes = NULL,
geneToModule = NULL,
drvNames = NULL,
genotFitness = NULL,
keepInput = TRUE) {
if(!is.null(genotFitness)) {
if(!is.null(rT) || !is.null(epistasis) ||
!is.null(orderEffects) || !is.null(noIntGenes) ||
!is.null(geneToModule)) {
stop("You have a non-null genotFitness.",
" If you pass the complete genotype to fitness mapping",
" you cannot pass any of rT, epistasis, orderEffects",
" noIntGenes or geneToModule.")
}
genotFitness_std <- to_genotFitness_std(genotFitness, simplify = TRUE)
## epistasis <- from_genotype_fitness(genotFitness)
} else {
genotFitness_std <- NULL
}
allFitnessORMutatorEffects(
rT = rT,
epistasis = epistasis,
orderEffects = orderEffects,
noIntGenes = noIntGenes,
geneToModule = geneToModule,
drvNames = drvNames,
keepInput = keepInput,
genotFitness = genotFitness_std,
calledBy = "allFitnessEffects")
}
## Former version
## allFitnessEffects <- function(rT = NULL,
## epistasis = NULL,
## orderEffects = NULL,
## noIntGenes = NULL,
## geneToModule = NULL,
## drvNames = NULL,
## genotFitness = NULL,
## keepInput = TRUE) {
## if(!is.null(genotFitness)) {
## if(!is.null(rT) || !is.null(epistasis) ||
## !is.null(orderEffects) || !is.null(noIntGenes) ||
## !is.null(geneToModule)) {
## stop("You have a non-null genotFitness.",
## " If you pass the complete genotype to fitness mapping",
## " you cannot pass any of rT, epistasis, orderEffects",
## " noIntGenes or geneToModule.")
## }
## epistasis <- from_genotype_fitness(genotFitness)
## }
## allFitnessORMutatorEffects(
## rT = rT,
## epistasis = epistasis,
## orderEffects = orderEffects,
## noIntGenes = noIntGenes,
## geneToModule = geneToModule,
## drvNames = drvNames,
## keepInput = keepInput,
## calledBy = "allFitnessEffects")
## }
## allFitnessEffects <- function(rT = NULL,
## epistasis = NULL,
## orderEffects = NULL,
## noIntGenes = NULL,
## geneToModule = NULL,
## drvNames = NULL,
## keepInput = TRUE) {
## ## restrictions: the usual rt
## ## epistasis: as it says, with the ":"
## ## orderEffects: the ">"
## ## All of the above can be genes or can be modules (if you pass a
## ## geneToModule)
## ## rest: rest of genes, with fitness
## ## For epistasis and order effects we create the output object but
## ## missing the numeric ids of genes. With rT we do it in one go, as we
## ## already know the mapping of genes to numeric ids. We could do the
## ## same in epistasis and order, but we would be splitting twice
## ## (whereas for rT extracting the names is very simple).
## rtNames <- NULL
## epiNames <- NULL
## orNames <- NULL
## if(!is.null(rT)) {
## ## This is really ugly, but to prevent the stringsAsFactors I need it here:
## rT$parent <- as.character(rT$parent)
## rT$child <- as.character(rT$child)
## rT$typeDep <- as.character(rT$typeDep)
## rtNames <- unique(c(rT$parent, rT$child))
## }
## if(!is.null(epistasis)) {
## long.epistasis <- to.long.epist.order(epistasis, ":")
## ## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
## ## deal with the possible negative signs
## epiNames <- setdiff(unique(
## unlist(lapply(long.epistasis,
## function(x) lapply(x$ids,
## function(z) strsplit(z, "^-"))))),
## "")
## } else {
## long.epistasis <- list()
## }
## if(!is.null(orderEffects)) {
## long.orderEffects <- to.long.epist.order(orderEffects, ">")
## orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
## } else {
## long.orderEffects <- list()
## }
## allModuleNames <- unique(c(rtNames, epiNames, orNames))
## if(is.null(geneToModule)) {
## gMOneToOne <- TRUE
## geneToModule <- geneModuleNull(allModuleNames)
## } else {
## gMOneToOne <- FALSE
## if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
## stop(paste("Some values in geneToModule not present in any of",
## " rT, epistasis, or order effects"))
## if(any(is.na(match(allModuleNames, names(geneToModule)))))
## stop(paste("Some values in rT, epistasis, ",
## "or order effects not in geneToModule"))
## }
## geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
## idm <- unique(geneModule$ModuleNumID)
## names(idm) <- unique(geneModule$Module)
## if(!is.null(rT)) {
## checkRT(rT)
## long.rt <- to.long.rt(rT, idm)
## } else {
## long.rt <- list() ## yes, we want an object of length 0
## }
## ## Append the numeric ids to epistasis and order
## if(!is.null(epistasis)) {
## long.epistasis <- lapply(long.epistasis,
## function(x)
## addIntID.epist.order(x, idm,
## sort = TRUE,
## sign = TRUE))
## }
## if(!is.null(orderEffects)) {
## long.orderEffects <- lapply(long.orderEffects,
## function(x)
## addIntID.epist.order(x, idm,
## sort = FALSE,
## sign = FALSE))
## }
## if(!is.null(noIntGenes)) {
## mg <- max(geneModule[, "GeneNumID"])
## gnum <- seq_along(noIntGenes) + mg
## if(!is.null(names(noIntGenes))) {
## ng <- names(noIntGenes)
## if(any(ng %in% geneModule[, "Gene"] ))
## stop("A gene in noIntGenes also present in the other terms")
## } else {
## ng <- gnum
## }
## geneNoInt <- data.frame(Gene = as.character(ng),
## GeneNumID = gnum,
## s = noIntGenes,
## stringsAsFactors = FALSE)
## } else {
## geneNoInt <- data.frame()
## }
## if( (length(long.rt) + length(long.epistasis) + length(long.orderEffects) +
## nrow(geneNoInt)) == 0)
## stop("You have specified nothing!")
## if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
## graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
## } else {
## graphE <- NULL
## }
## if(!is.null(drvNames)) {
## drv <- getGeneIDNum(geneModule, geneNoInt, drvNames)
## } else {
## drv <- geneModule$GeneNumID[-1]
## }
## if(!keepInput) {
## rT <- epistasis <- orderEffects <- noIntGenes <- NULL
## }
## out <- list(long.rt = long.rt,
## long.epistasis = long.epistasis,
## long.orderEffects = long.orderEffects,
## long.geneNoInt = geneNoInt,
## geneModule = geneModule,
## gMOneToOne = gMOneToOne,
## geneToModule = geneToModule,
## graph = graphE,
## drv = drv,
## rT = rT,
## epistasis = epistasis,
## orderEffects = orderEffects,
## noIntGenes = noIntGenes
## )
## class(out) <- c("fitnessEffects")
## return(out)
## }
## No longer used
## rtAndGeneModule <- function(mdeps, gM = NULL) {
## ## To show a table of restrictions when there are modules. Do not use
## ## for anything else. Maybe as intermediate to plotting.
## ## Specify restriction table of modules and a mapping of modules to
## ## genes. gM is a named vector; names are modules, values are elements
## ## of each module.
## ## We do nothing important if gM is NULL except checks
## ## If there are modules, the table shows the individual genes.
## checkRT(mdeps)
## ## if(ncol(mdeps) != 5)
## ## stop("mdeps must be of exactly 5 columns")
## ## if(!identical(colnames(mdeps), c("parent", "child", "s", "sh", "typeDep")))
## ## stop(paste("Column names of mdeps not of appropriate format. ",
## ## "Should be parent, child, s, sh, typeDep"))
## if(!is.null(gM)) {
## if(any(is.na(match(mdeps[ , 1], names(gM)))))
## stop("Some values in parent not from a known module")
## if(any(is.na(match(mdeps[ , 2], names(gM)))))
## stop("Some values in child not from a known module")
## if(any(is.na(match(names(gM), c(mdeps[, 1], mdeps[, 2])))))
## stop("Some values in module in neither parent or child")
## parent <- gM[mdeps[, 1]]
## child <- gM[mdeps[, 2]]
## df <- data.frame(parent = parent,
## child = child,
## s = mdeps$s,
## sh = mdeps$sh,
## typeDep = mdeps$typeDep,
## stringsAsFactors = FALSE)
## } else {
## df <- mdeps
## }
## rownames(df) <- seq.int(nrow(df))
## return(df)
## }
## wrap.test.rt <- function(rt, gM = NULL) {
## ## FIXME add epistasis and orderEffects
## lrt <- allFitnessEffects(rt, geneToModule = gM)
## ## wrap_test_rt(lrt$long.rt)
## wrap_test_rt(lrt$long.rt, lrt$geneModule)
## }
## No longer used
## wrap.readFitnessEffects <- function(rt, epi, oe, ni, gm, echo = TRUE) {
## tt <- allFitnessEffects(rt, epi, oe, ni, gm)
## readFitnessEffects(tt, echo = echo)
## ## readFitnessEffects(tt$long.rt,
## ## tt$long.epistasis,
## ## tt$long.orderEffects,
## ## tt$long.geneNoInt,
## ## tt$geneModule,
## ## tt$gMOneToOne,
## ## echo = TRUE)
## }
evalGenotypeORMut <- function(genotype,
fmEffects,
verbose = FALSE,
echo = FALSE,
model = "",
calledBy_= NULL) {
## genotype can be a vector of integers, that are the exact same in
## the table of fmEffects or a vector of strings, or a vector (a
## string) with genes separated by "," or ">"
if( !(calledBy_ %in% c("evalGenotype", "evalGenotypeMut") ))
stop("How did you call this function?. Bug.")
## fmEffects could be a mutator effect
if(!exists("fitnessLandscape_gene_id", where = fmEffects)) {
fmEffects$fitnessLandscape_df <- data.frame()
fmEffects$fitnessLandscape_gene_id <- data.frame()
}
if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
(nrow(fmEffects$fitnessLandscape_df) > 0)) {
warning("Bozic model passing a fitness landscape will not work",
" for now.")
}
if(echo)
cat(paste("Genotype: ", genotype))
if(!is.integer(genotype)) {
if(length(grep(">", genotype))) {
genotype <- nice.vector.eo(genotype, ">")
} else if(length(grep(",", genotype))) {
genotype <- nice.vector.eo(genotype, ",")
}
all.g.nums <- c(fmEffects$geneModule$GeneNumID,
fmEffects$long.geneNoInt$GeneNumID,
fmEffects$fitnessLandscape_gene_id$GeneNumID)
all.g.names <- c(fmEffects$geneModule$Gene,
fmEffects$long.geneNoInt$Gene,
fmEffects$fitnessLandscape_gene_id$Gene)
genotype <- all.g.nums[match(genotype, all.g.names)]
}
if(any(is.na(genotype)))
stop("genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
if(!length(genotype))
stop("genotypes must have at least one mutated gene")
if(any(genotype < 0)) {
stop(paste("genotypes cannot contain negative values.",
"If you see this message, you found a bug."))
}
if(model %in% c("Bozic", "bozic1", "bozic2") )
prodNeg <- TRUE
else
prodNeg <- FALSE
ff <- evalRGenotype(rG = genotype,
rFE = fmEffects,
verbose = verbose,
prodNeg = prodNeg,
calledBy_ = calledBy_)
if(echo) {
if(calledBy_ == "evalGenotype") {
if(!prodNeg)
cat(" Fitness: ", ff, "\n")
else
cat(" Death rate: ", ff, "\n")
} else if(calledBy_ == "evalGenotypeMut") {
cat(" Mutation rate product :", ff, "\n")
}
}
return(ff)
}
evalGenotype <- function(genotype, fitnessEffects,
verbose = FALSE,
echo = FALSE,
model = "") {
if(inherits(fitnessEffects, "mutatorEffects"))
stop("You are trying to get the fitness of a mutator specification. ",
"You did not pass an object of class fitnessEffects.")
evalGenotypeORMut(genotype = genotype,
fmEffects = fitnessEffects,
verbose = verbose,
echo = echo,
model = model ,
calledBy_= "evalGenotype"
)
}
evalGenotypeFitAndMut <- function(genotype,
fitnessEffects,
mutatorEffects,
verbose = FALSE,
echo = FALSE,
model = "") {
## Must deal with objects from previous, pre flfast, modifications
if(!exists("fitnessLandscape_gene_id", where = fitnessEffects)) {
fitnessEffects$fitnessLandscape_df <- data.frame()
fitnessEffects$fitnessLandscape_gene_id <- data.frame()
}
if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
(nrow(fitnessEffects$fitnessLandscape_df) > 0)) {
warning("Bozic model passing a fitness landscape will not work",
" for now.")
}
prodNeg <- FALSE
## Next is from evalGenotypeAndMut
if(echo)
cat(paste("Genotype: ", genotype))
if(!is.integer(genotype)) {
if(length(grep(">", genotype))) {
genotype <- nice.vector.eo(genotype, ">")
} else if(length(grep(",", genotype))) {
genotype <- nice.vector.eo(genotype, ",")
}
all.g.nums <- c(fitnessEffects$geneModule$GeneNumID,
fitnessEffects$long.geneNoInt$GeneNumID,
fitnessEffects$fitnessLandscape_gene_id$GeneNumID)
all.g.names <- c(fitnessEffects$geneModule$Gene,
fitnessEffects$long.geneNoInt$Gene,
fitnessEffects$fitnessLandscape_gene_id$Gene)
genotype <- all.g.nums[match(genotype, all.g.names)]
}
if(any(is.na(genotype)))
stop("genotype contains NAs or genes not in fitnessEffects")
if(!length(genotype))
stop("genotypes must have at least one mutated gene")
if(any(genotype < 0)) {
stop(paste("genotypes cannot contain negative values.",
"If you see this message, you found a bug."))
}
full2mutator_ <- matchGeneIDs(mutatorEffects,
fitnessEffects)$Reduced
if(model %in% c("Bozic", "bozic1", "bozic2") )
prodNeg <- TRUE
else
prodNeg <- FALSE
evalRGenotypeAndMut(genotype, fitnessEffects,
mutatorEffects, full2mutator_,
verbose = verbose,
prodNeg = prodNeg)
}
## evalGenotype <- function(genotype, fitnessEffects,
## verbose = FALSE,
## echo = FALSE,
## model = "") {
## ## genotype can be a vector of integers, that are the exact same in
## ## the table of fitnessEffects or a vector of strings, or a vector (a
## ## string) with genes separated by "," or ">"
## if(echo)
## cat(paste("Genotype: ", genotype))
## if(!is.integer(genotype)) {
## if(length(grep(">", genotype))) {
## genotype <- nice.vector.eo(genotype, ">")
## } else if(length(grep(",", genotype))) {
## genotype <- nice.vector.eo(genotype, ",")
## }
## all.g.nums <- c(fitnessEffects$geneModule$GeneNumID,
## fitnessEffects$long.geneNoInt$GeneNumID)
## all.g.names <- c(fitnessEffects$geneModule$Gene,
## fitnessEffects$long.geneNoInt$Gene)
## genotype <- all.g.nums[match(genotype, all.g.names)]
## }
## if(any(is.na(genotype)))
## stop("genotype contains NAs or genes not in fitnessEffects")
## if(!length(genotype))
## stop("genotypes must have at least one mutated gene")
## if(any(genotype < 0)) {
## stop(paste("genotypes cannot contain negative values.",
## "If you see this message, you found a bug."))
## }
## if(model %in% c("Bozic", "bozic1", "bozic2") )
## prodNeg <- TRUE
## else
## prodNeg <- FALSE
## ff <- evalRGenotype(genotype, fitnessEffects, verbose, prodNeg)
## if(echo) {
## if(!prodNeg)
## cat(" Fitness: ", ff, "\n")
## else
## cat(" Death rate: ", ff, "\n")
## } ## else {
## ## return(ff)
## ## }
## return(ff)
## }
## For multiple genotypes, lapply the matching.
## Nope, I think unneeded
## internal.convert_genotypes <- function(genotypes, gm) {
## genotypes <- lapply(lg, function(x) gm$GeneNumID[match(x, gm$Gene)])
## }
## I am here: simplify this
evalAllGenotypesORMut <- function(fmEffects,
order = FALSE, max = 256,
addwt = FALSE,
model = "",
calledBy_ = "") {
## minimal = FALSE) {
if( !(calledBy_ %in% c("evalGenotype", "evalGenotypeMut") ))
stop("How did you call this function?. Bug.")
if( (calledBy_ == "evalGenotype" ) &&
(!inherits(fmEffects, "fitnessEffects")))
stop("You are trying to get the fitness of a mutator specification. ",
"You did not pass an object of class fitnessEffects.")
if( (calledBy_ == "evalGenotypeMut" ) &&
(!inherits(fmEffects, "mutatorEffects")))
stop("You are trying to get the mutator effects of a fitness specification. ",
"You did not pass an object of class mutatorEffects.")
## if(!minimal)
allg <- generateAllGenotypes(fitnessEffects = fmEffects,
order = order, max = max)
## else
## allg <- generateAllGenotypes_minimal(fitnessEffects = fmEffects,
## max = max)
## if(order)
## tot <- function(n) {sum(sapply(seq.int(n),
## function(x) choose(n, x) * factorial(x)))}
## else
## tot <- function(n) {2^n}
## nn <- nrow(fitnessEffects$geneModule) -1 + nrow(fitnessEffects$long.geneNoInt)
## tnn <- tot(nn)
## if(tnn > max) {
## m <- paste("There are", tnn, "genotypes.")
## m <- paste(m, "This is larger than max.")
## m <- paste(m, "Adjust max and rerun if you want")
## stop(m)
## }
## ## With mutator, the ids of genes need not go from 1:n
## vid <- allNamedGenes(fitnessEffects)$GeneNumID
## if(order) {
## f1 <- function(n) {
## lapply(seq.int(n), function(x) permutations(n = n, r = x, v = vid))
## }
## } else {
## f1 <- function(n) {
## lapply(seq.int(n), function(x) combinations(n = n, r = x, v = vid))}
## }
## genotNums <- f1(nn)
## list.of.vectors <- function(y) {
## ## there's got to be a simpler way
## lapply(unlist(lapply(y, function(x) {apply(x, 1, list)}), recursive = FALSE),
## function(m) m[[1]])
## }
## genotNums <- list.of.vectors(genotNums)
## names <- c(fitnessEffects$geneModule$Gene[-1],
## fitnessEffects$long.geneNoInt$Gene)
## genotNames <- unlist(lapply(lapply(genotNums, function(x) names[x]),
## function(z)
## paste(z,
## collapse = if(order){" > "} else {", "} )))
## This ain't the best, as we repeatedly read and convert
## fitnessEffects. If this were slow, prepare C++ function that takes
## vectors of genotypes and uses same fitnessEffects.
if(model %in% c("Bozic", "bozic1", "bozic2") )
prodNeg <- TRUE
else
prodNeg <- FALSE
allf <- vapply(allg$genotNums,
function(x) evalRGenotype(x, fmEffects, FALSE,
prodNeg,
calledBy_),
1.1)
df <- data.frame(Genotype = allg$genotNames, Fitness = allf,
stringsAsFactors = FALSE)
## Why am I doing this? I am not computing the genotype. I test the
## evaluation of the empty genotype in the tests. But its evaluation
## generates warnings that are not needed. And by decree (even in the
## C++) this thing has a fitness of 1, a mutator effect of 1 since
## there are no terms.
if(addwt)
df <- rbind(data.frame(Genotype = "WT", Fitness = 1,
stringsAsFactors = FALSE), df)
if(calledBy_ == "evalGenotype") {
if(prodNeg)
colnames(df)[match("Fitness", colnames(df))] <- "Death_rate"
class(df) <- c("evalAllGenotypes", class(df))
} else if (calledBy_ == "evalGenotypeMut") {
colnames(df)[match("Fitness", colnames(df))] <- "MutatorFactor"
class(df) <- c("evalAllGenotypesMut", class(df))
}
return(df)
}
evalAllGenotypes <- function(fitnessEffects, order = FALSE, max = 256,
addwt = FALSE,
model = "") {
## Must deal with objects from previous, pre flfast, modifications
if(!exists("fitnessLandscape_gene_id", where = fitnessEffects)) {
fitnessEffects$fitnessLandscape_df <- data.frame()
fitnessEffects$fitnessLandscape_gene_id <- data.frame()
}
if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
(nrow(fitnessEffects$fitnessLandscape_df) > 0)) {
warning("Bozic model passing a fitness landscape will not work",
" for now.")
}
evalAllGenotypesORMut(
fmEffects = fitnessEffects,
order = order,
max = max,
addwt = addwt,
model = model,
calledBy_= "evalGenotype"
)
}
generateAllGenotypes <- function(fitnessEffects, order = TRUE, max = 256) {
if(order)
tot <- function(n) {sum(sapply(seq.int(n),
function(x) choose(n, x) * factorial(x)))}
else
tot <- function(n) {2^n}
nn <- nrow(fitnessEffects$geneModule) -1 +
nrow(fitnessEffects$long.geneNoInt) +
nrow(fitnessEffects$fitnessLandscape_gene_id)
tnn <- tot(nn)
if(tnn > max) {
m <- paste("There are", tnn, "genotypes.")
m <- paste(m, "This is larger than max.")
m <- paste(m, "Adjust max and rerun if you want")
stop(m)
}
## With mutator, the ids of genes need not go from 1:n
vid <- allNamedGenes(fitnessEffects)$GeneNumID
if(order) {
f1 <- function(n) {
lapply(seq.int(n), function(x) permutations(n = n, r = x, v = vid))
}
} else {
f1 <- function(n) {
lapply(seq.int(n), function(x) combinations(n = n, r = x, v = vid))}
}
genotNums <- f1(nn)
list.of.vectors <- function(y) {
## there's got to be a simpler way
lapply(unlist(lapply(y, function(x) {apply(x, 1, list)}), recursive = FALSE),
function(m) m[[1]])
}
genotNums <- list.of.vectors(genotNums)
names <- c(fitnessEffects$geneModule$Gene[-1],
fitnessEffects$long.geneNoInt$Gene,
fitnessEffects$fitnessLandscape_gene_id$Gene)
genotNames <- unlist(lapply(lapply(genotNums, function(x) names[x]),
function(z)
paste(z,
collapse = if(order){" > "} else {", "} )))
## This ain't the best, as we repeatedly read and convert
## fitnessEffects. If this were slow, prepare C++ function that takes
## vectors of genotypes and uses same fitnessEffects.
return(list(genotNums = genotNums,
genotNames = genotNames
))
}
evalAllGenotypesFitAndMut <- function(fitnessEffects, mutatorEffects,
order = FALSE, max = 256,
addwt = FALSE,
model = "" ){
## minimal = FALSE) {
## if(!minimal)
allg <- generateAllGenotypes(fitnessEffects = fitnessEffects,
order = order, max = max)
## else
## allg <- generateAllGenotypes_minimal(fitnessEffects = fitnessEffects,
## max = max)
if(model %in% c("Bozic", "bozic1", "bozic2") ) {
prodNeg <- TRUE
} else {
prodNeg <- FALSE
}
## Must deal with objects from previous, pre flfast, modifications
if(!exists("fitnessLandscape_gene_id", where = fitnessEffects)) {
fitnessEffects$fitnessLandscape_df <- data.frame()
fitnessEffects$fitnessLandscape_gene_id <- data.frame()
}
if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
(nrow(fitnessEffects$fitnessLandscape_df) > 0)) {
warning("Bozic model passing a fitness landscape will not work",
" for now.")
}
full2mutator_ <- matchGeneIDs(mutatorEffects,
fitnessEffects)$Reduced
allf <- t(vapply(allg$genotNums,
function(x) evalRGenotypeAndMut(x,
rFE = fitnessEffects,
muEF= mutatorEffects,
full2mutator_ = full2mutator_,
verbose = FALSE,
prodNeg = prodNeg),
c(1.1, 2.2)))
df <- data.frame(Genotype = allg$genotNames, Fitness = allf[, 1],
MutatorFactor = allf[, 2],
stringsAsFactors = FALSE)
if(addwt)
df <- rbind(data.frame(Genotype = "WT", Fitness = 1,
MutatorFactor = 1,
stringsAsFactors = FALSE), df)
if(prodNeg)
colnames(df)[match("Fitness", colnames(df))] <- "Death_rate"
class(df) <- c("evalAllGenotypesFitAndMut", class(df))
return(df)
}
## evalAllGenotypes <- function(fitnessEffects, order = TRUE, max = 256,
## addwt = FALSE,
## model = "") {
## if(order)
## tot <- function(n) {sum(sapply(seq.int(n),
## function(x) choose(n, x) * factorial(x)))}
## else
## tot <- function(n) {2^n}
## nn <- nrow(fitnessEffects$geneModule) -1 + nrow(fitnessEffects$long.geneNoInt)
## tnn <- tot(nn)
## if(tnn > max) {
## m <- paste("There are", tnn, "genotypes.")
## m <- paste(m, "This is larger than max.")
## m <- paste(m, "Adjust max and rerun if you want")
## stop(m)
## }
## if(order) {
## f1 <- function(n) {
## lapply(seq.int(n), function(x) permutations(n = n, r = x))
## }
## } else {
## f1 <- function(n) {
## lapply(seq.int(n), function(x) combinations(n = n, r = x))}
## }
## genotNums <- f1(nn)
## list.of.vectors <- function(y) {
## ## there's got to be a simpler way
## lapply(unlist(lapply(y, function(x) {apply(x, 1, list)}), recursive = FALSE),
## function(m) m[[1]])
## }
## genotNums <- list.of.vectors(genotNums)
## names <- c(fitnessEffects$geneModule$Gene[-1],
## fitnessEffects$long.geneNoInt$Gene)
## genotNames <- unlist(lapply(lapply(genotNums, function(x) names[x]),
## function(z)
## paste(z,
## collapse = if(order){" > "} else {", "} )))
## ## This ain't the best, as we repeatedly read and convert
## ## fitnessEffects. If this were slow, prepare C++ function that takes
## ## vectors of genotypes and uses same fitnessEffects.
## if(model %in% c("Bozic", "bozic1", "bozic2") )
## prodNeg <- TRUE
## else
## prodNeg <- FALSE
## allf <- vapply(genotNums,
## function(x) evalRGenotype(x, fitnessEffects, FALSE, prodNeg),
## 1.1)
## df <- data.frame(Genotype = genotNames, Fitness = allf,
## stringsAsFactors = FALSE)
## if(addwt)
## df <- rbind(data.frame(Genotype = "WT", Fitness = 1,
## stringsAsFactors = FALSE), df)
## if(prodNeg)
## colnames(df)[match("Fitness", colnames(df))] <- "Death_rate"
## return(df)
## }
fitnessEffectsToIgraph <- function(rT, epistasis, orderEffects) {
df0 <- df1 <- df2 <- data.frame()
if(!is.null(rT)) {
df0 <- rT[, c("parent", "child", "typeDep")]
}
if(!is.null(epistasis)) {
df1 <- to.pairs.modules(epistasis, ":")
}
if(!is.null(orderEffects)) {
df2 <- to.pairs.modules(orderEffects, ">")
}
df <- rbind(df0, df1, df2)
## for special case of simple epi effects
if(nrow(df) == 0) return(NULL)
g1 <- graph.data.frame(df)
E(g1)$color <- "black"
E(g1)$color[E(g1)$typeDep == "SM"] <- "blue"
E(g1)$color[E(g1)$typeDep == "XMPN"] <- "red"
E(g1)$color[E(g1)$typeDep == "epistasis"] <- "orange"
E(g1)$color[E(g1)$typeDep == "orderEffect"] <- "orange"
E(g1)$lty <- 1
E(g1)$lty[E(g1)$typeDep == "epistasis"] <- 2
E(g1)$lty[E(g1)$typeDep == "orderEffect"] <- 3
E(g1)$arrow.mode <- 2
E(g1)$arrow.mode[E(g1)$typeDep == "epistasis"] <- 0
E(g1)$arrow.mode[E(g1)$typeDep == "orderEffect"] <- 2
return(g1)
}
plot.fitnessEffects <- function(x, type = "graphNEL",
layout = NULL,
expandModules = FALSE,
autofit = FALSE,
scale_char = ifelse(type == "graphNEL", 1/10, 5),
return_g = FALSE,
lwdf = 1, ## multiplier for lwd for graphNEL
...) {
## some other layouts I find OK
## layout.circle
## layout.reingold.tilford if really a tree
## o.w. it will use the default
g <- x$graph
if(is.null(g))
stop("This fitnessEffects object can not be ploted this way.",
" It is probably one with fitness landscape specification, ",
" so you might want to plot the fitness landscape instead.")
if(type == "igraph") {
if(expandModules && (!x$gMOneToOne)) {
## vlabels <- fe$geneToModule[vertex.attributes(g)$name]
vlabels <- x$geneToModule[V(g)$name]
V(g)$label <- vlabels
}
if(autofit) {
plot(0, type = "n", axes = FALSE, ann = FALSE)
## ideas from http://stackoverflow.com/questions/14472079/match-vertex-size-to-label-size-in-igraph
## vsize <- (strwidth(V(g)$label) + strwidth("oo")) * 200
## but this is a kludge.
vsize <- (nchar(V(g)$label) + 1) * scale_char
plot.igraph(g, vertex.size = vsize, vertex.shape = "rectangle",
layout = layout)
} else {
plot.igraph(g, layout = layout)
}
if(return_g) return(g)
}
else if (type == "graphNEL") {
g1 <- igraph.to.graphNEL(g)
c1 <- unlist(lapply(edgeData(g1), function(x) x$color))
names(c1) <- sub("|", "~", names(c1), fixed = TRUE)
s1 <- unlist(lapply(edgeData(g1), function(x) x$lty))
names(s1) <- names(c1)
a1 <- unlist(lapply(edgeData(g1), function(x) max(x$arrow.mode - 1, 0)))
names(a1) <- names(c1)
lwd <- s1
lwd[lwd == 2] <- 2 ## o.w. too thin
lwd[lwd == 3] <- 2 ## o.w. too thin
lwd <- lwd * lwdf
nAttrs <- list()
if(expandModules && (!x$gMOneToOne)) {
nnodes <- x$geneToModule[nodes(g1)]
names(nnodes) <- nodes(g1)
nAttrs$label <- nnodes
}
if(autofit) {
nAttrs$width <- (nchar(nAttrs$label) + 1) * scale_char
names(nAttrs$width) <- names(nAttrs$label)
plot(g1, edgeAttrs = list(arrowsize = a1, lty = s1, lwd = lwd,
color = c1), attrs=list(node=list(shape = "rectangle")),
nodeAttrs = nAttrs)
} else {
plot(g1, edgeAttrs = list(arrowsize = a1, lty = s1, lwd = lwd,
color = c1),
nodeAttrs = nAttrs)
}
if(return_g) return(g1)
} else {
stop("plot type not recognized")
}
}
## plot.fitnessEffects <- plotFitnessEffects
## FIXME: in the help: say we cannot properly show 3- or higher order
## interactions.
nr_oncoSimul.internal <- function(rFE,
birth,
death,
mu,
initSize,
sampleEvery,
detectionSize,
finalTime,
initSize_species,
initSize_iter,
seed,
verbosity,
speciesFS,
ratioForce,
typeFitness,
max.memory,
mutationPropGrowth,
initMutant,
max.wall.time,
keepEvery,
K,
detectionDrivers,
onlyCancer,
errorHitWallTime,
max.num.tries,
errorHitMaxTries,
minDetectDrvCloneSz,
extraTime,
keepPhylog,
detectionProb,
AND_DrvProbExit,
MMUEF = NULL, ## avoid partial matching, and set default
fixation = NULL ## avoid partial matching
) {
default_min_successive_fixation <- 50 ## yes, set at this for now
if(!inherits(rFE, "fitnessEffects"))
stop(paste("rFE must be an object of class fitnessEffects",
"as created, for instance, with function",
"allFitnessEffects"))
if(countGenesFe(rFE) < 2) {
stop("There must be at least two genes (loci) in the fitness effects.",
"If you only care about a degenerate case with just one,",
"you can enter a second gene (locus)",
"with fitness effect of zero.")
}
if( (length(mu) == 1) && !(is.null(names(mu)))) {
stop("A length 1 mutation, but named. ",
"This is ambiguous. ",
"If you want per-gene mutation rates, each gene",
"must have its entry in the mu vector. ",
"(And regardless of per-gene mutation rates ",
" there must be at least two gene/loci).")
}
## Must deal with objects from previous, pre flfast, modifications
## Could move this way down to the bottom, right before
## .Call
if(!exists("fitnessLandscape_gene_id", where = rFE)) {
rFE$fitnessLandscape_df <- data.frame()
rFE$fitnessLandscape_gene_id <- data.frame()
}
namedGenes <- allNamedGenes(rFE)
if( length(mu) > 1) {
if(is.null(names(mu)))
stop("When using per-gene mutation rates the ",
"mu vector must be named ",
"(and if you have noIntGenes, those must have names).")
if(length(mu) != countGenesFe(rFE))
stop("When using per-gene mutation rates, ",
"there must be the same number of genes in the ",
"mu vector and the fitness effects object.")
if(!identical(sort(namedGenes$Gene),
sort(names(mu))))
stop("When using per-gene mutation rates, ",
"names of genes must match those in the ",
"restriction table.")
mu <- mu[order(match(names(mu), namedGenes$Gene))]
## Hyperparanoid check. Should never, ever, happen.
if(!identical(names(mu), namedGenes$Gene))
stop("Names of re-ordered mu do not match names of genes")
minmu <- 1e-40
if(any(mu <= minmu))
stop(paste("At least one per-gene mutation rate is negative",
"or less than", minmu,". Remember that the per-base",
"mutation rate in the human genome is about 1e-11 to 1e-9."))
}
if(!is.null(initMutant)) {
initMutantString <- initMutant
if(length(grep(">", initMutant))) {
initMutant <- nice.vector.eo(initMutant, ">")
} else if(length(grep(",", initMutant))) {
initMutant <- nice.vector.eo(initMutant, ",")
}
initMutant <- getGeneIDNum(rFE$geneModule,
rFE$long.geneNoInt,
rFE$fitnessLandscape_gene_id,
initMutant,
FALSE)
if(length(initMutant) >= countGenesFe(rFE)) {
stop("For initMutant you passed as many, or more genes, mutated ",
"than the number of genes in the genotype (fitness effects).")
}
} else {
initMutant <- vector(mode = "integer")
initMutantString <- ""
}
## these are never user options
## if(initSize_species < 10) {
## warning("initSize_species too small?")
## }
## if(initSize_iter < 100) {
## warning("initSize_iter too small?")
## }
if(typeFitness %in% c("bozic1", "bozic2")) {
if(nrow(rFE$fitnessLandscape_df) > 0)
warning("Bozic model passing a fitness landscape will not work",
" for now.")
## FIXME: bozic and fitness landscape
## the issue is that in the C++ code we directly do
## s = birth rate - 1
## but we would need something different
## Can be done going through epistasis, etc.
thesh <- unlist(lapply(rFE$long.rt, function(x) x$sh))
## thes <- unlist(lapply(rFE$long.rt, function(x) x$s))
thes <- unlist(c(lapply(rFE$long.rt, function(x) x$s),
lapply(rFE$long.epistasis, function(x) x$s),
lapply(rFE$long.orderEffects, function(x) x$s),
rFE$long.geneNoInt$s))
if(any(thes > 1 )) {
m <- paste("You are using a Bozic model with",
"the new restriction specification, and you have",
"at least one s > 1.",
"But that is the same as setting that to 1:",
"we obviously cannot allow negative death rates,",
"nor problems derived from multiplying odd or even",
"numbers of negative numbers.",
"You can set those > 1 to 1, but then you will",
"eventually get the simulations to abort when the",
"death rate becomes 0.")
stop(m)
}
if(any( (thesh <= -1) & (thesh > -Inf))) {
m <- paste("You are using a Bozic model with",
"the new restriction specification, and you have",
"at least one sh <= -1. Maybe you mean -Inf?")
warning(m)
}
if(any(thes == 1 )) {
m <- paste("You are using a Bozic model with",
"the new restriction specification, and you have",
"at least one s of 1. If that gene is mutated, ",
"this will lead to a death rate of 0",
"and the simulations will abort when you get a non finite",
"value.")
warning(m)
}
}
if(!is.null(MMUEF)) {
if(!inherits(MMUEF, "mutatorEffects"))
stop("muEF must be a mutatorEffects object")
full2mutator_ <- matchGeneIDs(MMUEF, rFE)
if(any(is.na(full2mutator_$Full)))
stop("Genes in mutatorEffects not present in fitnessEffects")
if(any(is.na(full2mutator_)))
stop("full2mutator failed for unknown reasons")
full2mutator_ <- full2mutator_$Reduced
} else {
full2mutator_ <- vector(mode = "numeric", length = 0)
## muEF <- emptyFitnessEffects()
}
dpr <- detectionProbCheckParse(detectionProb, initSize, verbosity)
## if( !is.null(cPDetect) && (sum(!is.null(p2), !is.null(n2)) >= 1 ))
## stop("Specify only cPDetect xor both of p2 and n2")
## if( (is.null(p2) + is.null(n2)) == 1 )
## stop("If you pass one of n2 or p2, you must also pass the other. ",
## "Otherwise, we would not know what to do.")
## stopifnot(PDBaseline >= 0)
## stopifnot(n2 > PDBaseline)
## stopifnot(p2 < 1)
## stopifnot(p2 > 0)
## if( is.null(cPDetect) ) cPDetect <- -9
## if( is.null(p2)) p2 <- 9
## if( is.null(n2)) n2 <- -9
## call <- match.call()
## Process the fixed list, if any
if(!is_null_na(fixation)) {
ng <- namedGenes
## rownames(ng) <- namedGenes[, "Gene"]
## Proposed extension to have exact matching of genotypes
ng <- rbind(
data.frame(Gene = "_", GeneNumID = -9, stringsAsFactors = FALSE),
ng)
rownames(ng) <- ng[, "Gene"]
## FIXME
## Later, accept a last argument, called tolerance.
## If not present, set to 0
## and then, at at the head of fixation_list below
if(is.list(fixation)) {
if(
(is.null(fixation[["fixation_tolerance"]])) ||
(is.na(fixation[["fixation_tolerance"]]))) {
fixation_tolerance <- 0
} else {
fixation_tolerance <- as.numeric(fixation[["fixation_tolerance"]])
fixation <- fixation[-which(names(fixation) == "fixation_tolerance")]
}
if(
(is.null(fixation[["min_successive_fixation"]])) ||
(is.na(fixation[["min_successive_fixation"]]))) {
min_successive_fixation <- default_min_successive_fixation
} else {
min_successive_fixation <- as.integer(fixation[["min_successive_fixation"]])
fixation <- fixation[-which(names(fixation) == "min_successive_fixation")]
}
if(
(is.null(fixation[["fixation_min_size"]])) ||
(is.na(fixation[["fixation_min_size"]]))) {
fixation_min_size <- 0
} else {
fixation_min_size <- as.integer(fixation[["fixation_min_size"]])
fixation <- fixation[-which(names(fixation) == "fixation_min_size")]
}
} else {
if(is_null_na(fixation["fixation_tolerance"])) {
fixation_tolerance <- 0
} else {
fixation_tolerance <- as.numeric(fixation["fixation_tolerance"])
fixation <- fixation[-which(names(fixation) == "fixation_tolerance")]
}
if(is_null_na(fixation["min_successive_fixation"])) {
min_successive_fixation <- default_min_successive_fixation
} else {
min_successive_fixation <- as.integer(fixation["min_successive_fixation"])
fixation <- fixation[-which(names(fixation) == "min_successive_fixation")]
}
if(is_null_na(fixation["fixation_min_size"])) {
fixation_min_size <- 0
} else {
fixation_min_size <- as.integer(fixation["fixation_min_size"])
fixation <- fixation[-which(names(fixation) == "fixation_min_size")]
}
}
if( (fixation_tolerance > 1) || (fixation_tolerance < 0) )
stop("Impossible range for fixation tolerance")
if( (min_successive_fixation < 0) )
stop("Impossible range for min_successive_fixation")
if( (fixation_min_size < 0) )
stop("Impossible range for fixation_min_size")
## Usual genotype specification and might allow ordered vectors
## in the future
fixation_b <- lapply(fixation, nice.vector.eo, sep = ",")
ulf <- unlist(fixation_b)
if(any(ulf == ">"))
stop("Order effects not allowed (yet?) in fixation.")
ulfg <- ng[ulf, 1]
if(any(is.na(ulfg)))
stop(paste("The 'fixation' list contains genes that are not present",
" in the fitness effects."))
## Sorting here is crucial!!
fixation_list <- lapply(fixation_b, function(x) sort(ng[x, 2]))
fixation_list <- list(fixation_list = fixation_list,
fixation_tolerance = fixation_tolerance,
min_successive_fixation = min_successive_fixation,
fixation_min_size = fixation_min_size)
} else {
fixation_list <- list()
}
return(c(
nr_BNB_Algo5(rFE = rFE,
mu_ = mu,
death = death,
initSize = initSize,
sampleEvery = sampleEvery,
detectionSize = detectionSize,
finalTime = finalTime,
initSp = initSize_species,
initIt = initSize_iter,
seed = seed,
verbosity = verbosity,
speciesFS = speciesFS,
ratioForce = ratioForce,
typeFitness_ = typeFitness,
maxram = max.memory,
mutationPropGrowth = as.integer(mutationPropGrowth),
initMutant_ = initMutant,
maxWallTime = max.wall.time,
keepEvery = keepEvery,
K = K,
detectionDrivers = detectionDrivers,
onlyCancer = onlyCancer,
errorHitWallTime = errorHitWallTime,
maxNumTries = max.num.tries,
errorHitMaxTries = errorHitMaxTries,
minDetectDrvCloneSz = minDetectDrvCloneSz,
extraTime = extraTime,
keepPhylog = keepPhylog,
MMUEF = MMUEF,
full2mutator_ = full2mutator_,
n2 = dpr["n2"],
p2 = dpr["p2"],
PDBaseline = dpr["PDBaseline"],
cPDetect_i= dpr["cPDetect"],
checkSizePEvery = dpr["checkSizePEvery"],
AND_DrvProbExit = AND_DrvProbExit,
fixation_list = fixation_list),
Drivers = list(rFE$drv), ## but when doing pops, these will be repeated
geneNames = list(names(getNamesID(rFE))),
InitMutant = initMutantString
))
}
countGenesFe <- function(fe) {
## recall geneModule has Root always
## We want to be able to use objects that did not have
## the fitness landscape component
if(exists("fitnessLandscape_gene_id", where = fe))
return(nrow(fe$geneModule) + nrow(fe$long.geneNoInt) +
nrow(fe$fitnessLandscape_gene_id) - 1)
else
return(nrow(fe$geneModule) + nrow(fe$long.geneNoInt) - 1)
}
allNamedGenes <- function(fe){
## Returns a data frame with genes and their names and verifies all
## genes have names.
## Root should always be first, but just in case
## avoid assuming it
## This does is not a good idea as it assumes the user did not use
## "keepInput = FALSE".
## lni <- length(fe$noIntGenes)
## ## FIXME:test
## if((lni > 0) &&
## (is.null(names(fe$noIntGenes))))
## stop("When using per-gene mutation rates the ",
## "no interaction genes must be named ",
## "(i.e., the noIntGenes vector must have names).")
## accommodate objects w/o fitnessLandscape
if(!is.null(fe$fitnessLandscape) && nrow(fe$fitnessLandscape)) {
gn <-
gtools::mixedsort(
colnames(fe$fitnessLandscape)[-ncol(fe$fitnessLandscape)])
v1 <- data.frame(Gene = gn, GeneNumID = seq.int(length(gn)),
stringsAsFactors = FALSE)
} else {
v1 <- fe$geneModule[, c("Gene", "GeneNumID")]
if(nrow(fe$long.geneNoInt)) {
v1 <- rbind(v1,
fe$long.geneNoInt[, c("Gene", "GeneNumID")])
}
if(any(v1[, "Gene"] == "Root"))
v1 <- v1[-which(v1[, "Gene"] == "Root"), ]
}
rownames(v1) <- NULL
if(any(v1$Gene == "WT")) {
warning("A gene is named WT. You can expect problems ",
"because we use WT to denote the wildtype ",
"genotype. You might want to change it.")
}
return(v1)
}
get.gene.counts <- function(x) {
# , timeSample = "last",
# typeSample = "whole") {
## From get.mut.vector. Used for now just for testing
timeSample <- "last"
the.time <- get.the.time.for.sample(x, timeSample, NULL)
if(the.time < 0) {
counts <- rep(0, length(x$geneNames))
names(counts) <- x$geneNames
freq <- rep(NA, length(x$geneNames))
names(freq) <- x$geneNames
return(list(counts = counts,
freq = freq,
popSize = 0))
}
pop <- x$pops.by.time[the.time, -1]
if(all(pop == 0)) {
stop("You found a bug: this should never happen")
}
## if(typeSample %in% c("wholeTumor", "whole")) {
popSize <- x$PerSampleStats[the.time, 1]
counts <- as.vector(tcrossprod(pop, x$Genotypes))
names(counts) <- x$geneNames
return(list(counts = counts,
freq = counts/popSize,
popSize = popSize))
## return( (tcrossprod(pop,
## x$Genotypes)/popSize) )
## } else if (typeSample %in% c("singleCell", "single")) {
## return(x$Genotypes[, sample(seq_along(pop), 1, prob = pop)])
## } else {
## stop("Unknown typeSample option")
## }
}
geneCounts <- function(x) {
if(inherits(x, "oncosimulpop")) {
return( do.call(rbind,
lapply(x,
function(z)
get.gene.counts(z)$counts))
)
} else {
return( get.gene.counts(x)$counts)
}
}
## geneNumIDReset <- function(x, ref){
## ## Set GeneNumID of a fitnessEffect object, x, using ref as the
## ## reference fitness effect object.
## ## Check also if all in x are in ref.
## gg <- allNamedGenes(ref)
## gnid <- gg$GeneNumID
## names(gnid) <- gg$Gene
## ## FIXME: this later and conditional on what is in thee
## gnid <- c("Root" = 0, gnid)
## if(!all(x$geneModule$Gene %in% names(gnid) ))
## stop("Some genes not in reference fitnessEffects (rebasing geneModule)")
## x$geneModule$GeneNumID <- gnid[geneModule$Gene]
## ## and then go over all the lists in the x object.
## if(nrow(x$long.geneNoInt)) {
## ## now, mapping for the noInt if this is mutator
## if(!all(x$long.geneNoInt$Gene %in% names(gnid) ))
## stop("Some genes not in reference fitnessEffects (rebasing geneNoInt)")
## x$long.geneNoInt$GeneNumID <- gnid[long.geneNoInt$Gene]
## }
## }
matchGeneIDs <- function(x, refFE) {
n1 <- allNamedGenes(x)
n2 <- allNamedGenes(refFE)
colnames(n1)[2] <- "Reduced"
colnames(n2)[2] <- "Full"
## Non standard evaluation thing and a note being generated in check.
## See also http://stackoverflow.com/questions/33299845/when-writing-functions-for-an-r-package-should-i-avoid-non-standard-evaluation
## But does not work well with replace. So use the NULL trick
Reduced <- NULL
dplyr::full_join(n2, n1, by = "Gene") %>%
mutate(Reduced = replace(Reduced, is.na(Reduced), -9))
}
detectionProbCheckParse <- function(x, initSize, verbosity) {
default_p2 <- 0.1
default_n2 <- 2 * initSize
default_PDBaseline <- 1.2 * initSize
default_checkSizePEvery <- 20
## No default cPDetect. That is done from p2 and n2 in C++.
if(is_null_na(x)) {
## if((length(x) == 1) && (is.na(x))) {
y <- vector()
y["cPDetect"] <- -9
y["p2"] <- 9
y["n2"] <- -9
y["PDBaseline"] <- -9
y["checkSizePEvery"] <- Inf
return(y)
} else if((length(x) == 1) && (x == "default")) {
## Populate with defaults
y <- vector()
y["p2"] <- default_p2
y["n2"] <- default_n2
y["PDBaseline"] <- default_PDBaseline
y["checkSizePEvery"] <- default_checkSizePEvery
x <- y
}
if(length(x) >= 1) {
if( !all(names(x) %in% c("cPDetect", "p2", "n2", "PDBaseline", "checkSizePEvery")))
stop("Names of some components of detectionProb are not recognized.",
" Check for possible typos.")
}
## This ain't conditional. If not given, always check
if( !is_null_na(x["cPDetect"]) && (sum(!is_null_na(x["p2"]), !is_null_na(x["n2"])) >= 1 ))
stop("Specify only cPDetect xor both of p2 and n2")
if( (is_null_na(x["p2"]) + is_null_na(x["n2"])) == 1 )
stop("If you pass one of n2 or p2, you must also pass the other. ",
"Otherwise, we would not know what to do.")
if(is_null_na(x["PDBaseline"])) {
x["PDBaseline"] <- default_PDBaseline
if(verbosity > -1)
message("\n PDBaseline set to default value of ", default_PDBaseline, "\n")
}
if(is_null_na(x["checkSizePEvery"])) {
x["checkSizePEvery"] <- default_checkSizePEvery
if(verbosity > -1)
message("\n checkSizePEvery set to default value of ",
default_checkSizePEvery, "\n")
}
## If we get here, we checked consistency of whether cPDetect or p2/n2.
## Fill up with defaults the missing values
if(is_null_na(x["cPDetect"])) {
if(is_null_na(x["p2"])) {
if(!is_null_na(x["n2"])) stop("Eh? no p2 but n2? Bug")
x["n2"] <- default_n2
x["p2"] <- default_p2
if(verbosity > -1)
message("\n n2, p2 set to default values of ",
default_n2, ", ", default_p2, "\n")
}
}
if(x["checkSizePEvery"] <= 0)
stop("checkSizePEvery <= 0")
if(x["PDBaseline"] <= 0)
stop("PDBaseline <= 0")
if(!is_null_na(x["n2"])) {
if(x["n2"] <= x["PDBaseline"])
stop("n2 <= PDBaseline")
if(x["p2"] >= 1) stop("p2 >= 1")
if(x["p2"] <= 0) stop("p2 <= 0")
x["cPDetect"] <- -9
} else { ## only if x["cPDetect"] is not NA
if(is_null_na(x["cPDetect"])) stop("eh? you found a bug")## paranoia
x["n2"] <- -9
x["p2"] <- -9
if(verbosity > -1)
message("\n Using user-specified cPDetect as n2, p2 not given.\n")
}
return(x)
}
sampledGenotypes <- function(y, genes = NULL) {
## From a popSample object, or a matrix for that matter,
## show the sampled genotypes and their frequencies
if(!is.null(genes)) {
cols <- which(colnames(y) %in% genes )
y <- y[, cols]
}
## nn <- colnames(y)
genotlabel <- function(u, nn = colnames(y)) {
if(any(is.na(u))) return(NA)
else {
return(paste(sort(nn[as.logical(u)]), collapse = ", "))
}
}
## df <- data.frame(table(
## apply(y, 1, function(z) paste(sort(nn[as.logical(z)]), collapse = ", ") )
## ))
df <- data.frame(table(
apply(y, 1, genotlabel),
useNA = "ifany"),
stringsAsFactors = TRUE)
gn <- as.character(df[, 1])
gn[gn == ""] <- "WT"
df <- data.frame(Genotype = gn, Freq = df[, 2], stringsAsFactors = FALSE)
attributes(df)$ShannonI <- shannonI(na.omit(df)$Freq)
class(df) <- c("sampledGenotypes", class(df))
return(df)
}
print.sampledGenotypes <- function(x, ...) {
print.data.frame(x, ...)
cat("\n Shannon's diversity (entropy) of sampled genotypes: ")
cat(attributes(x)$ShannonI, "\n")
}
list_g_matches_fixed <- function(x, y) {
## Internal function, for testing the fixation output.
## x and y are vectors
## x is the set of output genotypes, y the set of fixed
## genotypes/subset of genotypes.
## Yes, this function has tests too in test.fixation.R
## All genotypes in x satisfy that they are supersets of at least one
## in y? That is true if, for every element in x, at least one y in
## that x.
if(is.list(y)) y <- unlist(y)
y.nice <- lapply(y, nice.vector.eo, sep = ",")
x.nice <- lapply(x, nice.vector.eo, sep = ",")
fu <- function(u, y.nice)
any(unlist(lapply(y.nice, function(z) all(z %in% u))))
return(all(unlist(lapply(x.nice, fu, y.nice))))
}
## emptyFitnessEffects <- function() {
## list(long.rt = list(),
## long.epistasis = list(),
## long.orderEffects = list(),
## long.geneNoInt = list(),
## geneModule = list(),
## gMOneToOne = TRUE,
## geneToModule = list(),
## graph = NULL,
## drv = vector(mode = "integer", length = 0)
## )
## }
### Later, for all the effects, we will do some kind of dplyr match?
### a, b, c, in fitness, only a, c in mut.
### fitness table for a,b,c
### each row name transformed removing b (so leaving only present)
### each row transformed matched to row in mut table.
## t1 <- data.frame(v1 = c("a,b", "a,c", "b"), v2 = c("b", "c", "b"), v3 = 1:3, stringsAsFactors = FALSE)
## t2 <- data.frame(v2 = c("b", "c"), v4 = c(11, 12), stringsAsFactors = FALSE)
## full_join(t1, t2, by = "v2")
## FIXME
## new exit code
## check:
## baseline < n2
## p2 < 1
## baseline defaults to init size + .1
## use c < -9, n2 < -9, p2 < -9 for no values
## and check one of c or p2 and n2 are valid if using exit
Try the OncoSimulR package in your browser
Any scripts or data that you put into this service are public.
OncoSimulR documentation built on Nov. 8, 2020, 8:31 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 list_g_matches_fixed print.sampledGenotypes sampledGenotypes detectionProbCheckParse matchGeneIDs geneCounts get.gene.counts allNamedGenes countGenesFe nr_oncoSimul.internal plot.fitnessEffects fitnessEffectsToIgraph evalAllGenotypesFitAndMut generateAllGenotypes evalAllGenotypes evalAllGenotypesORMut evalGenotypeFitAndMut evalGenotype evalGenotypeORMut allFitnessEffects allFitnessORMutatorEffects getGeneIDNum getNamesID checkRT addIntID.epist.order to.long.epist.order to.pairs.modules epist.order.to.pairs.modules oe.to.df epist.order.element to.long.rt list.of.deps geneModuleNull gm.to.geneModuleL nice.vector.eo gtm2 check.gm
Documented in allFitnessEffects evalAllGenotypes evalAllGenotypesFitAndMut evalGenotype evalGenotypeFitAndMut plot.fitnessEffects print.sampledGenotypes sampledGenotypes
## Copyright 2013, 2014, 2015, 2016 Ramon Diaz-Uriarte
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
## 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 for more details.
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.
## Say which are drivers: populate the drv vector.
## // In R, the user says which are the drivers. If does not say anuthing,
## // the default (NULL) then drivers are all in poset, epist, restrict. The
## // user can pass a vector with the names of the genes (not modules). Allow
## // also for empty, so this is faster if not needed. And check that if we
## // use a stopping rule based on drivers that drv vectors is not empty.
## - Say that a user can use a "0" as a gene name, but that is BAD idea.
## - Modules and order effects can be kind of funny?
## Gene names can contain no spaces, commas, or ">"
check.gm <- function(gm) {
## Yes, Root: we want no ambiguities.
## Actually, that sucks. So we do not require it, but check for
## consistency.
if(any(gm == "Root") && (gm[1] != "Root") )
stop("If Root is in the module table, it must be the first")
if(any(names(gm) == "Root") && (names(gm)[1] != "Root") )
stop("If the name Root is in the module table, it must be the first")
if( (names(gm)[1] == "Root") && (gm[1] != "Root") )
stop("The name Root is in the module table, but is not of Root")
if( (gm[1] == "Root") && (names(gm)[1] != "Root") )
stop("Root is in the module table, but with a different name")
## if(gm[1] != "Root")
## stop("First value of a module table must be Root")
## if(names(gm)[1] != "Root")
## stop("First name of a module table must be Root")
if(length(unique(names(gm))) != length(gm))
stop("Number of unique module names different from length of vector")
if(gm[1] != "Root")
gm <- c("Root" = "Root", gm)
return(gm)
}
gtm2 <- function(x) {
data.frame(cbind(nice.vector.eo(x, ","), x), stringsAsFactors = TRUE)
}
## nice.vector.eo <- function(z, sep) {
## ## with epistasis, maybe we want sorted?
## setdiff(unlist(lapply(strsplit(z, " "),
## function(u) strsplit(u, sep))), "")
## }
nice.vector.eo <- function(z, sep, rm.sign = FALSE) {
## with epistasis, maybe we want sorted?
if(! rm.sign)
return(setdiff(unlist(lapply(strsplit(z, " "),
function(u) strsplit(u, sep))), ""))
else ## remove the " ", the -, and the sep
return(setdiff(unlist(lapply(strsplit(z, " "), function(u)
strsplit(unlist(strsplit(u, "-")), sep))), ""))
}
gm.to.geneModuleL <- function(gm, one.to.one) {
## the table will be sorted by gene name
gm <- check.gm(gm)
## the named vector with the mapping into the long geneModule df
geneMod <- as.data.frame(rbindlist(lapply(gm, gtm2)),
stringsAsFactors = TRUE) ## this stringsAsFactors is key
geneMod$Module <- names(gm)[geneMod[, 2]] ## reverse lookup table
colnames(geneMod)[1] <- c("Gene")
geneMod <- geneMod[, -2]
geneMod$Gene <- as.character(geneMod$Gene)
## geneMod$Module <- as.character(geneMod$Module) ## already a char
geneMod <- geneMod[c(1, order(geneMod$Gene[-1]) + 1), ]
geneMod$GeneNumID <- 0:(nrow(geneMod) - 1)
## this assumes sorted! and they need not be
## rlem <- rle(geneMod$Module)
## geneMod$ModuleNumID <- rep( 0:(length(rlem$values) - 1), rlem$lengths)
if(one.to.one) {
geneMod$ModuleNumID <- geneMod$GeneNumID
if(!(identical(geneMod$Gene, geneMod$Module)))
stop("Impossible: we are supposed to be in one-to-one for Module-Gene.")
} else {
## Works, but does not give the most ordered module names. But
## keeps implicit order given by user.
idm <- seq.int(length(gm) - 1)
idm <- c("Root" = 0L, idm)
names(idm) <- names(gm)
## This sorts by the names but is not optimal either
## idm1 <- seq.int(length(gm) - 1)
## idm <- c(0L, idm1)
## names(idm) <- c("Root", sort(setdiff(names(gm), "Root")))
geneMod$ModuleNumID <- idm[geneMod[, "Module"]]
}
if(length(unique(geneMod$Gene)) != nrow(geneMod)) {
stop("Are there identical gene names in different modules?")
}
## I think this is unreacheable now. Caught earlier.
if(length(unique(geneMod$GeneNumID)) != nrow(geneMod)) {
stop("Are there identical gene names in different modules?")
}
rownames(geneMod) <- 1:nrow(geneMod)
geneMod
}
geneModuleNull <- function(namesM) {
v <- c("Root", setdiff(namesM, "Root"))
names(v) <- v
return(v)
}
list.of.deps <- function(x) {
## lookupTypeDep <- c("MN" = 1, "monotone" = 1,
## "SM" = 2, "semimonotone" = 2)
lookupTypeDep <- c("MN" = "monotone",
"AND" = "monotone",
"monotone" = "monotone",
"CMPN" = "monotone",
"OR" = "semimonotone",
"SM" = "semimonotone",
"semimonotone" = "semimonotone",
"DMPN" = "semimonotone",
"XOR" = "xmpn",
"xmpn" = "xmpn",
"XMPN" = "xmpn",
"--" = "--",
"-" = "--")
## FIXME: check values of typeDep
if(length(x) > 1) {
if(length(unique(x$s))!= 1)
stop("Not all s identical within a child")
if(length(unique(x$sh))!= 1)
stop("Not all sh identical within a child")
if(length(unique(x$typeDep))!= 1)
stop("Not all typeDep identical within a child")
if(length(unique(x$child))!= 1)
stop("child not unique")
}
typeDep <- lookupTypeDep[unique(x$typeDep)]
if(any(is_null_na(typeDep)))
stop("typeDep value incorrect. Check spelling.")
return(list(
child = unique(x$child),
s = unique(x$s),
sh = unique(x$sh),
typeDep = typeDep,
parents = unlist(x$parent)))
}
to.long.rt <- function(rt, idm) {
## We now do this inconditionally, so that we do not need to use the
## "stringsAsFactors = FALSE". This is now done before
## if(is.numeric(rt$parent))
## rt$parent <- as.character(rt$parent)
## if(is.numeric(rt$child))
## rt$child <- as.character(rt$child)
if(!("Root" %in% rt$parent))
stop("Root must be one parent node")
## rt$parent <- unlist(lapply(rt$parent, nice.string))
## rt$child <- unlist(lapply(rt$child, nice.string))
srt <- rt[order(rt$child), ]
## Not relevant if we allow non-numeric names
## all.child.genes <- as.integer(
## unlist(lapply(rt[, 2],
## function(x) strsplit(x, ","))))
## ## check all childs
## if(!identical(sort(unique(all.child.genes)),
## seq.int(max(all.child.genes))))
## stop("Not all children present")
long.rt <- lapply(split(srt, srt$child), list.of.deps)
## geneModule <- gene.to.module(srt)
## idm <- seq.int(length(names(long.rt)))
## names(idm) <- names(long.rt)
## idm <- c("0" = 0L, idm)
## geneModule$ModuleNumID <- idm[geneModule[, "Module"]]
## idm is just a look up table for the id of the module
## idm <- unique(geneModule$ModuleNumID)
## names(idm) <- unique(geneModule$Module)
## add integer IDs
addIntID <- function(z, idm) {
z$childNumID <- idm[z$child]
z$parentsNumID <- idm[z$parents]
if( any(is.na(z$parentsNumID)) ||
any(is.na(z$childNumID)) ) {
## I think this can no longer be reached ever. Caught before.
stop(paste("An ID is NA:",
"Is a gene part of two different modules?",
"(That includes being by itself and part",
"of a module.)"))
}
## These checks could be somewhere else instead of here.
if(length(unique(z$parentsNumID)) != length(z$parentsNumID))
stop("No node can have the same parent multiple times")
if(length(unique(z$parents)) != length(z$parents))
stop("No node can have the same parent multiple times")
if(length(z$child) > 1)
stop("Child nodes can have one, and only one, member")
## sort the parents, always.
o <- order(z$parentsNumID)
z$parentsNumID <- z$parentsNumID[o]
z$parents <- as.character(z$parents[o])
return(z)
}
long.rt <- lapply(long.rt, function(x) addIntID(x, idm = idm))
## if(verbosity >= 4) {
## message(paste("Number of drivers: ",
## length(unique(geneModule[, "Gene"]))))
## message(paste("Number of modules: ",
## length(unique(geneModule[, "Module"]))))
## }
return(long.rt)
## return(list(long.rt = long.rt, geneModule = geneModule))
}
epist.order.element <- function(x, y, sep, rm.sign = FALSE) {
list(ids = nice.vector.eo(x, sep = sep, rm.sign = rm.sign), s = y)
}
oe.to.df <- function(x) {
ma <- matrix(ncol = 2, nrow = 1 + length(x) - 2)
if(length(x) == 2) {
ma[1, 1] <- x[1]
ma[1, 2] <- x[2]
} else {
ma[, 1] <- x[-length(x)]
ma[, 2] <- x[-1]
}
return(data.frame(ma, stringsAsFactors = FALSE))
}
epist.order.to.pairs.modules <- function(x, sep, rm.sign = TRUE) {
## We discard, do not even accept, the coefficient
tmp <- epist.order.element(x, -99, sep = sep, rm.sign = rm.sign)$ids
if(length(tmp) > 1) {
## if a single gene, as when we specify all genotypes, we do not
## want this
if(sep == ":")
return(data.frame(combinations(n = length(tmp), r = 2, v = tmp),
stringsAsFactors = FALSE))
else if(sep == ">") {
return(oe.to.df(tmp))
}
}
}
to.pairs.modules <- function(x, sep, rm.sign = TRUE) {
df <- data.frame(rbindlist(
lapply(names(x),
function(z) epist.order.to.pairs.modules(z, sep))),
stringsAsFactors = TRUE)
if(nrow(df) == 0L) { ## if only single genes in epist, we get nothing here.
return(df)
} else {
colnames(df) <- c("parent", "child")
if(sep == ":")
df$typeDep <- "epistasis"
else if(sep == ">")
df$typeDep <- "orderEffect"
return(unique(df))
}
}
to.long.epist.order <- function(epor, sep, rm.sign = FALSE) {
## just vectors for now
long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
names(epor), epor)
## if(is.vector(epor))
## long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
## names(epor), epor)
## else if(is.data.frame(epor))
## long <- Map(function(x, y) epist.order.element(x, y, sep, rm.sign),
## as.character(epor$ids),
## epor$s)
names(long) <- NULL
return(long)
}
## addIntID.epist.order <- function(z, idm, sort) {
## z$NumID <- idm[z$ids]
## if(sort) {
## ## essential for epistasis, but never do for order effects
## o <- order(z$NumID)
## z$NumID <- z$NumID[o]
## z$ids <- z$ids[o]
## }
## return(z)
## }
addIntID.epist.order <- function(z, idm, sort, sign) {
if( sort && (!sign))
warning("sort is TRUE but sign is FALSE. You sure?")
if((!sort) && sign)
warning("sort is FALSE but sign is TRUE. You sure?")
## Adds the integer, but takes care of sign if sign = TRUE
if(!sign) {
z$NumID <- idm[z$ids]
signs <- grep("^-", z$ids)
if(length(signs))
stop("You have a negative sign and you said sign = FALSE")
} else {
unsigned <- setdiff(unlist(lapply(z$ids,
function(z) strsplit(z, "^-"))), "")
NumID <- idm[unsigned]
signs <- grep("^-", z$ids)
if(length(signs)) {
NumID[signs] <- -1 * NumID[signs]
}
z$NumID <- NumID
}
if(sort) {
## Essential for epistasis, but never do for order effects
o <- order(z$NumID)
z$NumID <- z$NumID[o]
z$ids <- z$ids[o]
}
if(length(unique(z$NumID)) != length(z$NumID))
stop("No node can have the same id multiple times")
if(length(unique(z$ids)) != length(z$ids))
stop("No node can have the same id multiple times")
return(z)
}
checkRT <- function(mdeps) {
if(ncol(mdeps) != 5)
stop("mdeps must be of exactly 5 columns")
if(!identical(colnames(mdeps), c("parent", "child", "s", "sh", "typeDep")))
stop(paste("Column names of mdeps not of appropriate format. ",
"Should be parent, child, s, sh, typeDep"))
}
getNamesID <- function(fp) {
## Return a lookup table for names based on numeric IDs
idname <- c(fp$geneModule$GeneNumID,
fp$long.geneNoInt$GeneNumID,
fp$fitnessLandscape_gene_id$GeneNumID)
names(idname) <- c(fp$geneModule$Gene,
fp$long.geneNoInt$Gene,
fp$fitnessLandscape_gene_id$Gene)
return(idname[-1]) ## remove Root!!
}
getGeneIDNum <- function(geneModule, geneNoInt, fitnessLandscape_gene_id,
drv, sort = TRUE) {
## It returns the genes, as NumID, in the given vector with names drv
## initMutant uses this, for simplicity, without sorting, but noInt
## are always sorted
## Also used for the drivers with sort = TRUE
## Yes, we must do it twice because we do not know before hand which
## is which. This makes sure no NA. Period.
if(any(is.na( match(drv, c(geneModule$Gene, geneNoInt$Gene,
fitnessLandscape_gene_id$Gene))))) {
stop(paste("For driver or initMutant you have passed genes",
"not in the fitness table."))
}
indicesM <- as.vector(na.omit(match( drv, geneModule$Gene)))
indicesI <- as.vector(na.omit(sort(match( drv, geneNoInt$Gene))))
indicesF <- as.vector(na.omit(sort(match( drv, fitnessLandscape_gene_id$Gene))))
if(sort) {
indicesM <- sort(indicesM)
}
return(c(
geneModule$GeneNumID[indicesM],
geneNoInt$GeneNumID[indicesI],
fitnessLandscape_gene_id$GeneNumID[indicesF])
)
}
allFitnessORMutatorEffects <- function(rT = NULL,
epistasis = NULL,
orderEffects = NULL,
noIntGenes = NULL,
geneToModule = NULL,
drvNames = NULL,
keepInput = TRUE,
genotFitness = NULL,
## refFE = NULL,
calledBy = NULL) {
## From allFitnessEffects. Generalized so we deal with Fitness
## and mutator.
## restrictions: the usual rt
## epistasis: as it says, with the ":"
## orderEffects: the ">"
## All of the above can be genes or can be modules (if you pass a
## geneToModule)
## rest: rest of genes, with fitness
## For epistasis and order effects we create the output object but
## missing the numeric ids of genes. With rT we do it in one go, as we
## already know the mapping of genes to numeric ids. We could do the
## same in epistasis and order, but we would be splitting twice
## (whereas for rT extracting the names is very simple).
## called appropriately?
if( !(calledBy %in% c("allFitnessEffects", "allMutatorEffects") ))
stop("How did you call this function?. Bug.")
if(calledBy == "allMutatorEffects") {
## very paranoid check
if( !is.null(rT) || !is.null(orderEffects) ||
!is.null(drvNames) || !is.null(genotFitness))
stop("allMutatorEffects called with forbidden arguments.",
"Is this an attempt to subvert the function?")
}
rtNames <- NULL
epiNames <- NULL
orNames <- NULL
if(!is.null(rT)) {
## This is really ugly, but to prevent the stringsAsFactors I need it here:
rT$parent <- as.character(rT$parent)
rT$child <- as.character(rT$child)
rT$typeDep <- as.character(rT$typeDep)
rtNames <- unique(c(rT$parent, rT$child))
}
if(!is.null(epistasis)) {
long.epistasis <- to.long.epist.order(epistasis, ":")
## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
## deal with the possible negative signs
epiNames <- setdiff(unique(
unlist(lapply(long.epistasis,
function(x) lapply(x$ids,
function(z) strsplit(z, "^-"))))),
"")
} else {
long.epistasis <- list()
}
if(!is.null(orderEffects)) {
long.orderEffects <- to.long.epist.order(orderEffects, ">")
orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
} else {
long.orderEffects <- list()
}
allModuleNames <- unique(c(rtNames, epiNames, orNames))
if(is.null(geneToModule)) {
gMOneToOne <- TRUE
geneToModule <- geneModuleNull(allModuleNames)
} else {
gMOneToOne <- FALSE
if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
stop(paste("Some values in geneToModule not present in any of",
" rT, epistasis, or order effects"))
if(any(is.na(match(allModuleNames, names(geneToModule)))))
stop(paste("Some values in rT, epistasis, ",
"or order effects not in geneToModule"))
}
geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
idm <- unique(geneModule$ModuleNumID)
names(idm) <- unique(geneModule$Module)
if(!is.null(rT)) {
checkRT(rT)
long.rt <- to.long.rt(rT, idm)
} else {
long.rt <- list() ## yes, we want an object of length 0
}
## Append the numeric ids to epistasis and order
if(!is.null(epistasis)) {
long.epistasis <- lapply(long.epistasis,
function(x)
addIntID.epist.order(x, idm,
sort = TRUE,
sign = TRUE))
}
if(!is.null(orderEffects)) {
long.orderEffects <- lapply(long.orderEffects,
function(x)
addIntID.epist.order(x, idm,
sort = FALSE,
sign = FALSE))
}
if(!is.null(noIntGenes)) {
if(inherits(noIntGenes, "character")) {
wm <- paste("noIntGenes is a character vector.",
"This is probably not what you want, and will",
"likely result in an error downstream.",
"You can get messages like",
" 'not compatible with requested type', and others.",
"We are stopping.")
stop(wm)
}
mg <- max(geneModule[, "GeneNumID"])
gnum <- seq_along(noIntGenes) + mg
if(!is.null(names(noIntGenes))) {
ng <- names(noIntGenes)
if( any(grepl(",", ng, fixed = TRUE)) ||
any(grepl(">", ng, fixed = TRUE)) ||
any(grepl(":", ng, fixed = TRUE)) )
stop("The name of some noIntGenes contain a ',' or a '>' or a ':'")
if(any(ng %in% geneModule[, "Gene"] ))
stop("A gene in noIntGenes also present in the other terms")
if(any(duplicated(ng)))
stop("Duplicated gene names in geneNoInt")
if(any(is.na(ng)))
stop("In noIntGenes some genes have names, some don't.",
" Name all of them, or name none of them.")
} else {
ng <- gnum
}
geneNoInt <- data.frame(Gene = as.character(ng),
GeneNumID = gnum,
s = noIntGenes,
stringsAsFactors = FALSE)
} else {
geneNoInt <- data.frame()
}
if(is.null(genotFitness)) {
genotFitness <- matrix(NA, nrow = 0, ncol = 1)
fitnessLandscape_df <- data.frame()
fitnessLandscape_gene_id <- data.frame()
} else {
## Yes, I am duplicating stuff for now.
## This makes life simpler in C++:
## In the map, the key is the genotype name, as
## cnn <- colnames(genotFitness)[-ncol(genotFitness)]
cnn <- 1:(ncol(genotFitness) - 1)
gfn <- apply(genotFitness[, -ncol(genotFitness), drop = FALSE], 1,
function(x) paste(cnn[as.logical(x)],
collapse = ", "))
## rownames(genotFitness) <- gfn
fitnessLandscape_df <-
data.frame(Genotype = gfn,
Fitness = genotFitness[, ncol(genotFitness)],
stringsAsFactors = FALSE)
fitnessLandscape_gene_id <- data.frame(
Gene = colnames(genotFitness)[-ncol(genotFitness)],
GeneNumID = cnn,
stringsAsFactors = FALSE)
}
if( (length(long.rt) + length(long.epistasis) + length(long.orderEffects) +
nrow(geneNoInt) + nrow(genotFitness)) == 0)
stop("You have specified nothing!")
if(calledBy == "allFitnessEffects") {
if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
} else {
graphE <- NULL
}
} else {
graphE <- NULL
}
if(!is.null(drvNames)) {
drv <- unique(getGeneIDNum(geneModule, geneNoInt, fitnessLandscape_gene_id,
drvNames))
## drivers should never be in the geneNoInt; Why!!!???
## Catch the problem. This is an overkill,
## so since we catch the issue, we could leave the geneNoInt. But
## that should not be there in this call.
## if(any(drvNames %in% geneNoInt$Gene)) {
## stop(paste("At least one gene in drvNames is a geneNoInt gene.",
## "That is not allowed.",
## "If that gene is a driver, pass it as gene in the epistasis",
## "component."))
## }
## drv <- getGeneIDNum(geneModule, NULL, drvNames)
} else {
## we used to have this default
## drv <- geneModule$GeneNumID[-1]
drv <- vector(mode = "integer", length = 0L)
}
if(!keepInput) {
rT <- epistasis <- orderEffects <- noIntGenes <- NULL
}
out <- list(long.rt = long.rt,
long.epistasis = long.epistasis,
long.orderEffects = long.orderEffects,
long.geneNoInt = geneNoInt,
geneModule = geneModule,
gMOneToOne = gMOneToOne,
geneToModule = geneToModule,
graph = graphE,
drv = drv,
rT = rT,
epistasis = epistasis,
orderEffects = orderEffects,
noIntGenes = noIntGenes,
fitnessLandscape = genotFitness,
fitnessLandscape_df = fitnessLandscape_df,
fitnessLandscape_gene_id = fitnessLandscape_gene_id
)
if(calledBy == "allFitnessEffects") {
class(out) <- c("fitnessEffects")
} else if(calledBy == "allMutatorEffects") {
class(out) <- c("mutatorEffects")
}
return(out)
}
## Former version, with fitness landscape
## allFitnessORMutatorEffects <- function(rT = NULL,
## epistasis = NULL,
## orderEffects = NULL,
## noIntGenes = NULL,
## geneToModule = NULL,
## drvNames = NULL,
## keepInput = TRUE,
## ## refFE = NULL,
## calledBy = NULL) {
## ## From allFitnessEffects. Generalized so we deal with Fitness
## ## and mutator.
## ## restrictions: the usual rt
## ## epistasis: as it says, with the ":"
## ## orderEffects: the ">"
## ## All of the above can be genes or can be modules (if you pass a
## ## geneToModule)
## ## rest: rest of genes, with fitness
## ## For epistasis and order effects we create the output object but
## ## missing the numeric ids of genes. With rT we do it in one go, as we
## ## already know the mapping of genes to numeric ids. We could do the
## ## same in epistasis and order, but we would be splitting twice
## ## (whereas for rT extracting the names is very simple).
## ## called appropriately?
## if( !(calledBy %in% c("allFitnessEffects", "allMutatorEffects") ))
## stop("How did you call this function?. Bug.")
## if(calledBy == "allMutatorEffects") {
## ## very paranoid check
## if( !is.null(rT) || !is.null(orderEffects) || !is.null(drvNames))
## stop("allMutatorEffects called with forbidden arguments.",
## "Is this an attempt to subvert the function?")
## }
## rtNames <- NULL
## epiNames <- NULL
## orNames <- NULL
## if(!is.null(rT)) {
## ## This is really ugly, but to prevent the stringsAsFactors I need it here:
## rT$parent <- as.character(rT$parent)
## rT$child <- as.character(rT$child)
## rT$typeDep <- as.character(rT$typeDep)
## rtNames <- unique(c(rT$parent, rT$child))
## }
## if(!is.null(epistasis)) {
## long.epistasis <- to.long.epist.order(epistasis, ":")
## ## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
## ## deal with the possible negative signs
## epiNames <- setdiff(unique(
## unlist(lapply(long.epistasis,
## function(x) lapply(x$ids,
## function(z) strsplit(z, "^-"))))),
## "")
## } else {
## long.epistasis <- list()
## }
## if(!is.null(orderEffects)) {
## long.orderEffects <- to.long.epist.order(orderEffects, ">")
## orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
## } else {
## long.orderEffects <- list()
## }
## allModuleNames <- unique(c(rtNames, epiNames, orNames))
## if(is.null(geneToModule)) {
## gMOneToOne <- TRUE
## geneToModule <- geneModuleNull(allModuleNames)
## } else {
## gMOneToOne <- FALSE
## if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
## stop(paste("Some values in geneToModule not present in any of",
## " rT, epistasis, or order effects"))
## if(any(is.na(match(allModuleNames, names(geneToModule)))))
## stop(paste("Some values in rT, epistasis, ",
## "or order effects not in geneToModule"))
## }
## geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
## idm <- unique(geneModule$ModuleNumID)
## names(idm) <- unique(geneModule$Module)
## if(!is.null(rT)) {
## checkRT(rT)
## long.rt <- to.long.rt(rT, idm)
## } else {
## long.rt <- list() ## yes, we want an object of length 0
## }
## ## Append the numeric ids to epistasis and order
## if(!is.null(epistasis)) {
## long.epistasis <- lapply(long.epistasis,
## function(x)
## addIntID.epist.order(x, idm,
## sort = TRUE,
## sign = TRUE))
## }
## if(!is.null(orderEffects)) {
## long.orderEffects <- lapply(long.orderEffects,
## function(x)
## addIntID.epist.order(x, idm,
## sort = FALSE,
## sign = FALSE))
## }
## if(!is.null(noIntGenes)) {
## if(inherits(noIntGenes, "character")) {
## wm <- paste("noIntGenes is a character vector.",
## "This is probably not what you want, and will",
## "likely result in an error downstream.",
## "You can get messages like",
## " 'not compatible with requested type', and others.",
## "We are stopping.")
## stop(wm)
## }
## mg <- max(geneModule[, "GeneNumID"])
## gnum <- seq_along(noIntGenes) + mg
## if(!is.null(names(noIntGenes))) {
## ng <- names(noIntGenes)
## if( grepl(",", ng, fixed = TRUE) || grepl(">", ng, fixed = TRUE)
## || grepl(":", ng, fixed = TRUE))
## stop("The name of some noIntGenes contain a ',' or a '>' or a ':'")
## if(any(ng %in% geneModule[, "Gene"] ))
## stop("A gene in noIntGenes also present in the other terms")
## if(any(duplicated(ng)))
## stop("Duplicated gene names in geneNoInt")
## if(any(is.na(ng)))
## stop("In noIntGenes some genes have names, some don't.",
## " Name all of them, or name none of them.")
## } else {
## ng <- gnum
## }
## geneNoInt <- data.frame(Gene = as.character(ng),
## GeneNumID = gnum,
## s = noIntGenes,
## stringsAsFactors = FALSE)
## } else {
## geneNoInt <- data.frame()
## }
## if( (length(long.rt) + length(long.epistasis) + length(long.orderEffects) +
## nrow(geneNoInt)) == 0)
## stop("You have specified nothing!")
## if(calledBy == "allFitnessEffects") {
## if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
## graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
## } else {
## graphE <- NULL
## }
## } else {
## graphE <- NULL
## }
## if(!is.null(drvNames)) {
## drv <- unique(getGeneIDNum(geneModule, geneNoInt, drvNames))
## ## drivers should never be in the geneNoInt; Why!!!???
## ## Catch the problem. This is an overkill,
## ## so since we catch the issue, we could leave the geneNoInt. But
## ## that should not be there in this call.
## ## if(any(drvNames %in% geneNoInt$Gene)) {
## ## stop(paste("At least one gene in drvNames is a geneNoInt gene.",
## ## "That is not allowed.",
## ## "If that gene is a driver, pass it as gene in the epistasis",
## ## "component."))
## ## }
## ## drv <- getGeneIDNum(geneModule, NULL, drvNames)
## } else {
## ## we used to have this default
## ## drv <- geneModule$GeneNumID[-1]
## drv <- vector(mode = "integer", length = 0L)
## }
## if(!keepInput) {
## rT <- epistasis <- orderEffects <- noIntGenes <- NULL
## }
## out <- list(long.rt = long.rt,
## long.epistasis = long.epistasis,
## long.orderEffects = long.orderEffects,
## long.geneNoInt = geneNoInt,
## geneModule = geneModule,
## gMOneToOne = gMOneToOne,
## geneToModule = geneToModule,
## graph = graphE,
## drv = drv,
## rT = rT,
## epistasis = epistasis,
## orderEffects = orderEffects,
## noIntGenes = noIntGenes
## )
## if(calledBy == "allFitnessEffects") {
## class(out) <- c("fitnessEffects")
## } else if(calledBy == "allMutatorEffects") {
## class(out) <- c("mutatorEffects")
## }
## return(out)
## }
allFitnessEffects <- function(rT = NULL,
epistasis = NULL,
orderEffects = NULL,
noIntGenes = NULL,
geneToModule = NULL,
drvNames = NULL,
genotFitness = NULL,
keepInput = TRUE) {
if(!is.null(genotFitness)) {
if(!is.null(rT) || !is.null(epistasis) ||
!is.null(orderEffects) || !is.null(noIntGenes) ||
!is.null(geneToModule)) {
stop("You have a non-null genotFitness.",
" If you pass the complete genotype to fitness mapping",
" you cannot pass any of rT, epistasis, orderEffects",
" noIntGenes or geneToModule.")
}
genotFitness_std <- to_genotFitness_std(genotFitness, simplify = TRUE)
## epistasis <- from_genotype_fitness(genotFitness)
} else {
genotFitness_std <- NULL
}
allFitnessORMutatorEffects(
rT = rT,
epistasis = epistasis,
orderEffects = orderEffects,
noIntGenes = noIntGenes,
geneToModule = geneToModule,
drvNames = drvNames,
keepInput = keepInput,
genotFitness = genotFitness_std,
calledBy = "allFitnessEffects")
}
## Former version
## allFitnessEffects <- function(rT = NULL,
## epistasis = NULL,
## orderEffects = NULL,
## noIntGenes = NULL,
## geneToModule = NULL,
## drvNames = NULL,
## genotFitness = NULL,
## keepInput = TRUE) {
## if(!is.null(genotFitness)) {
## if(!is.null(rT) || !is.null(epistasis) ||
## !is.null(orderEffects) || !is.null(noIntGenes) ||
## !is.null(geneToModule)) {
## stop("You have a non-null genotFitness.",
## " If you pass the complete genotype to fitness mapping",
## " you cannot pass any of rT, epistasis, orderEffects",
## " noIntGenes or geneToModule.")
## }
## epistasis <- from_genotype_fitness(genotFitness)
## }
## allFitnessORMutatorEffects(
## rT = rT,
## epistasis = epistasis,
## orderEffects = orderEffects,
## noIntGenes = noIntGenes,
## geneToModule = geneToModule,
## drvNames = drvNames,
## keepInput = keepInput,
## calledBy = "allFitnessEffects")
## }
## allFitnessEffects <- function(rT = NULL,
## epistasis = NULL,
## orderEffects = NULL,
## noIntGenes = NULL,
## geneToModule = NULL,
## drvNames = NULL,
## keepInput = TRUE) {
## ## restrictions: the usual rt
## ## epistasis: as it says, with the ":"
## ## orderEffects: the ">"
## ## All of the above can be genes or can be modules (if you pass a
## ## geneToModule)
## ## rest: rest of genes, with fitness
## ## For epistasis and order effects we create the output object but
## ## missing the numeric ids of genes. With rT we do it in one go, as we
## ## already know the mapping of genes to numeric ids. We could do the
## ## same in epistasis and order, but we would be splitting twice
## ## (whereas for rT extracting the names is very simple).
## rtNames <- NULL
## epiNames <- NULL
## orNames <- NULL
## if(!is.null(rT)) {
## ## This is really ugly, but to prevent the stringsAsFactors I need it here:
## rT$parent <- as.character(rT$parent)
## rT$child <- as.character(rT$child)
## rT$typeDep <- as.character(rT$typeDep)
## rtNames <- unique(c(rT$parent, rT$child))
## }
## if(!is.null(epistasis)) {
## long.epistasis <- to.long.epist.order(epistasis, ":")
## ## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
## ## deal with the possible negative signs
## epiNames <- setdiff(unique(
## unlist(lapply(long.epistasis,
## function(x) lapply(x$ids,
## function(z) strsplit(z, "^-"))))),
## "")
## } else {
## long.epistasis <- list()
## }
## if(!is.null(orderEffects)) {
## long.orderEffects <- to.long.epist.order(orderEffects, ">")
## orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
## } else {
## long.orderEffects <- list()
## }
## allModuleNames <- unique(c(rtNames, epiNames, orNames))
## if(is.null(geneToModule)) {
## gMOneToOne <- TRUE
## geneToModule <- geneModuleNull(allModuleNames)
## } else {
## gMOneToOne <- FALSE
## if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
## stop(paste("Some values in geneToModule not present in any of",
## " rT, epistasis, or order effects"))
## if(any(is.na(match(allModuleNames, names(geneToModule)))))
## stop(paste("Some values in rT, epistasis, ",
## "or order effects not in geneToModule"))
## }
## geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
## idm <- unique(geneModule$ModuleNumID)
## names(idm) <- unique(geneModule$Module)
## if(!is.null(rT)) {
## checkRT(rT)
## long.rt <- to.long.rt(rT, idm)
## } else {
## long.rt <- list() ## yes, we want an object of length 0
## }
## ## Append the numeric ids to epistasis and order
## if(!is.null(epistasis)) {
## long.epistasis <- lapply(long.epistasis,
## function(x)
## addIntID.epist.order(x, idm,
## sort = TRUE,
## sign = TRUE))
## }
## if(!is.null(orderEffects)) {
## long.orderEffects <- lapply(long.orderEffects,
## function(x)
## addIntID.epist.order(x, idm,
## sort = FALSE,
## sign = FALSE))
## }
## if(!is.null(noIntGenes)) {
## mg <- max(geneModule[, "GeneNumID"])
## gnum <- seq_along(noIntGenes) + mg
## if(!is.null(names(noIntGenes))) {
## ng <- names(noIntGenes)
## if(any(ng %in% geneModule[, "Gene"] ))
## stop("A gene in noIntGenes also present in the other terms")
## } else {
## ng <- gnum
## }
## geneNoInt <- data.frame(Gene = as.character(ng),
## GeneNumID = gnum,
## s = noIntGenes,
## stringsAsFactors = FALSE)
## } else {
## geneNoInt <- data.frame()
## }
## if( (length(long.rt) + length(long.epistasis) + length(long.orderEffects) +
## nrow(geneNoInt)) == 0)
## stop("You have specified nothing!")
## if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
## graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
## } else {
## graphE <- NULL
## }
## if(!is.null(drvNames)) {
## drv <- getGeneIDNum(geneModule, geneNoInt, drvNames)
## } else {
## drv <- geneModule$GeneNumID[-1]
## }
## if(!keepInput) {
## rT <- epistasis <- orderEffects <- noIntGenes <- NULL
## }
## out <- list(long.rt = long.rt,
## long.epistasis = long.epistasis,
## long.orderEffects = long.orderEffects,
## long.geneNoInt = geneNoInt,
## geneModule = geneModule,
## gMOneToOne = gMOneToOne,
## geneToModule = geneToModule,
## graph = graphE,
## drv = drv,
## rT = rT,
## epistasis = epistasis,
## orderEffects = orderEffects,
## noIntGenes = noIntGenes
## )
## class(out) <- c("fitnessEffects")
## return(out)
## }
## No longer used
## rtAndGeneModule <- function(mdeps, gM = NULL) {
## ## To show a table of restrictions when there are modules. Do not use
## ## for anything else. Maybe as intermediate to plotting.
## ## Specify restriction table of modules and a mapping of modules to
## ## genes. gM is a named vector; names are modules, values are elements
## ## of each module.
## ## We do nothing important if gM is NULL except checks
## ## If there are modules, the table shows the individual genes.
## checkRT(mdeps)
## ## if(ncol(mdeps) != 5)
## ## stop("mdeps must be of exactly 5 columns")
## ## if(!identical(colnames(mdeps), c("parent", "child", "s", "sh", "typeDep")))
## ## stop(paste("Column names of mdeps not of appropriate format. ",
## ## "Should be parent, child, s, sh, typeDep"))
## if(!is.null(gM)) {
## if(any(is.na(match(mdeps[ , 1], names(gM)))))
## stop("Some values in parent not from a known module")
## if(any(is.na(match(mdeps[ , 2], names(gM)))))
## stop("Some values in child not from a known module")
## if(any(is.na(match(names(gM), c(mdeps[, 1], mdeps[, 2])))))
## stop("Some values in module in neither parent or child")
## parent <- gM[mdeps[, 1]]
## child <- gM[mdeps[, 2]]
## df <- data.frame(parent = parent,
## child = child,
## s = mdeps$s,
## sh = mdeps$sh,
## typeDep = mdeps$typeDep,
## stringsAsFactors = FALSE)
## } else {
## df <- mdeps
## }
## rownames(df) <- seq.int(nrow(df))
## return(df)
## }
## wrap.test.rt <- function(rt, gM = NULL) {
## ## FIXME add epistasis and orderEffects
## lrt <- allFitnessEffects(rt, geneToModule = gM)
## ## wrap_test_rt(lrt$long.rt)
## wrap_test_rt(lrt$long.rt, lrt$geneModule)
## }
## No longer used
## wrap.readFitnessEffects <- function(rt, epi, oe, ni, gm, echo = TRUE) {
## tt <- allFitnessEffects(rt, epi, oe, ni, gm)
## readFitnessEffects(tt, echo = echo)
## ## readFitnessEffects(tt$long.rt,
## ## tt$long.epistasis,
## ## tt$long.orderEffects,
## ## tt$long.geneNoInt,
## ## tt$geneModule,
## ## tt$gMOneToOne,
## ## echo = TRUE)
## }
evalGenotypeORMut <- function(genotype,
fmEffects,
verbose = FALSE,
echo = FALSE,
model = "",
calledBy_= NULL) {
## genotype can be a vector of integers, that are the exact same in
## the table of fmEffects or a vector of strings, or a vector (a
## string) with genes separated by "," or ">"
if( !(calledBy_ %in% c("evalGenotype", "evalGenotypeMut") ))
stop("How did you call this function?. Bug.")
## fmEffects could be a mutator effect
if(!exists("fitnessLandscape_gene_id", where = fmEffects)) {
fmEffects$fitnessLandscape_df <- data.frame()
fmEffects$fitnessLandscape_gene_id <- data.frame()
}
if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
(nrow(fmEffects$fitnessLandscape_df) > 0)) {
warning("Bozic model passing a fitness landscape will not work",
" for now.")
}
if(echo)
cat(paste("Genotype: ", genotype))
if(!is.integer(genotype)) {
if(length(grep(">", genotype))) {
genotype <- nice.vector.eo(genotype, ">")
} else if(length(grep(",", genotype))) {
genotype <- nice.vector.eo(genotype, ",")
}
all.g.nums <- c(fmEffects$geneModule$GeneNumID,
fmEffects$long.geneNoInt$GeneNumID,
fmEffects$fitnessLandscape_gene_id$GeneNumID)
all.g.names <- c(fmEffects$geneModule$Gene,
fmEffects$long.geneNoInt$Gene,
fmEffects$fitnessLandscape_gene_id$Gene)
genotype <- all.g.nums[match(genotype, all.g.names)]
}
if(any(is.na(genotype)))
stop("genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
if(!length(genotype))
stop("genotypes must have at least one mutated gene")
if(any(genotype < 0)) {
stop(paste("genotypes cannot contain negative values.",
"If you see this message, you found a bug."))
}
if(model %in% c("Bozic", "bozic1", "bozic2") )
prodNeg <- TRUE
else
prodNeg <- FALSE
ff <- evalRGenotype(rG = genotype,
rFE = fmEffects,
verbose = verbose,
prodNeg = prodNeg,
calledBy_ = calledBy_)
if(echo) {
if(calledBy_ == "evalGenotype") {
if(!prodNeg)
cat(" Fitness: ", ff, "\n")
else
cat(" Death rate: ", ff, "\n")
} else if(calledBy_ == "evalGenotypeMut") {
cat(" Mutation rate product :", ff, "\n")
}
}
return(ff)
}
evalGenotype <- function(genotype, fitnessEffects,
verbose = FALSE,
echo = FALSE,
model = "") {
if(inherits(fitnessEffects, "mutatorEffects"))
stop("You are trying to get the fitness of a mutator specification. ",
"You did not pass an object of class fitnessEffects.")
evalGenotypeORMut(genotype = genotype,
fmEffects = fitnessEffects,
verbose = verbose,
echo = echo,
model = model ,
calledBy_= "evalGenotype"
)
}
evalGenotypeFitAndMut <- function(genotype,
fitnessEffects,
mutatorEffects,
verbose = FALSE,
echo = FALSE,
model = "") {
## Must deal with objects from previous, pre flfast, modifications
if(!exists("fitnessLandscape_gene_id", where = fitnessEffects)) {
fitnessEffects$fitnessLandscape_df <- data.frame()
fitnessEffects$fitnessLandscape_gene_id <- data.frame()
}
if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
(nrow(fitnessEffects$fitnessLandscape_df) > 0)) {
warning("Bozic model passing a fitness landscape will not work",
" for now.")
}
prodNeg <- FALSE
## Next is from evalGenotypeAndMut
if(echo)
cat(paste("Genotype: ", genotype))
if(!is.integer(genotype)) {
if(length(grep(">", genotype))) {
genotype <- nice.vector.eo(genotype, ">")
} else if(length(grep(",", genotype))) {
genotype <- nice.vector.eo(genotype, ",")
}
all.g.nums <- c(fitnessEffects$geneModule$GeneNumID,
fitnessEffects$long.geneNoInt$GeneNumID,
fitnessEffects$fitnessLandscape_gene_id$GeneNumID)
all.g.names <- c(fitnessEffects$geneModule$Gene,
fitnessEffects$long.geneNoInt$Gene,
fitnessEffects$fitnessLandscape_gene_id$Gene)
genotype <- all.g.nums[match(genotype, all.g.names)]
}
if(any(is.na(genotype)))
stop("genotype contains NAs or genes not in fitnessEffects")
if(!length(genotype))
stop("genotypes must have at least one mutated gene")
if(any(genotype < 0)) {
stop(paste("genotypes cannot contain negative values.",
"If you see this message, you found a bug."))
}
full2mutator_ <- matchGeneIDs(mutatorEffects,
fitnessEffects)$Reduced
if(model %in% c("Bozic", "bozic1", "bozic2") )
prodNeg <- TRUE
else
prodNeg <- FALSE
evalRGenotypeAndMut(genotype, fitnessEffects,
mutatorEffects, full2mutator_,
verbose = verbose,
prodNeg = prodNeg)
}
## evalGenotype <- function(genotype, fitnessEffects,
## verbose = FALSE,
## echo = FALSE,
## model = "") {
## ## genotype can be a vector of integers, that are the exact same in
## ## the table of fitnessEffects or a vector of strings, or a vector (a
## ## string) with genes separated by "," or ">"
## if(echo)
## cat(paste("Genotype: ", genotype))
## if(!is.integer(genotype)) {
## if(length(grep(">", genotype))) {
## genotype <- nice.vector.eo(genotype, ">")
## } else if(length(grep(",", genotype))) {
## genotype <- nice.vector.eo(genotype, ",")
## }
## all.g.nums <- c(fitnessEffects$geneModule$GeneNumID,
## fitnessEffects$long.geneNoInt$GeneNumID)
## all.g.names <- c(fitnessEffects$geneModule$Gene,
## fitnessEffects$long.geneNoInt$Gene)
## genotype <- all.g.nums[match(genotype, all.g.names)]
## }
## if(any(is.na(genotype)))
## stop("genotype contains NAs or genes not in fitnessEffects")
## if(!length(genotype))
## stop("genotypes must have at least one mutated gene")
## if(any(genotype < 0)) {
## stop(paste("genotypes cannot contain negative values.",
## "If you see this message, you found a bug."))
## }
## if(model %in% c("Bozic", "bozic1", "bozic2") )
## prodNeg <- TRUE
## else
## prodNeg <- FALSE
## ff <- evalRGenotype(genotype, fitnessEffects, verbose, prodNeg)
## if(echo) {
## if(!prodNeg)
## cat(" Fitness: ", ff, "\n")
## else
## cat(" Death rate: ", ff, "\n")
## } ## else {
## ## return(ff)
## ## }
## return(ff)
## }
## For multiple genotypes, lapply the matching.
## Nope, I think unneeded
## internal.convert_genotypes <- function(genotypes, gm) {
## genotypes <- lapply(lg, function(x) gm$GeneNumID[match(x, gm$Gene)])
## }
## I am here: simplify this
evalAllGenotypesORMut <- function(fmEffects,
order = FALSE, max = 256,
addwt = FALSE,
model = "",
calledBy_ = "") {
## minimal = FALSE) {
if( !(calledBy_ %in% c("evalGenotype", "evalGenotypeMut") ))
stop("How did you call this function?. Bug.")
if( (calledBy_ == "evalGenotype" ) &&
(!inherits(fmEffects, "fitnessEffects")))
stop("You are trying to get the fitness of a mutator specification. ",
"You did not pass an object of class fitnessEffects.")
if( (calledBy_ == "evalGenotypeMut" ) &&
(!inherits(fmEffects, "mutatorEffects")))
stop("You are trying to get the mutator effects of a fitness specification. ",
"You did not pass an object of class mutatorEffects.")
## if(!minimal)
allg <- generateAllGenotypes(fitnessEffects = fmEffects,
order = order, max = max)
## else
## allg <- generateAllGenotypes_minimal(fitnessEffects = fmEffects,
## max = max)
## if(order)
## tot <- function(n) {sum(sapply(seq.int(n),
## function(x) choose(n, x) * factorial(x)))}
## else
## tot <- function(n) {2^n}
## nn <- nrow(fitnessEffects$geneModule) -1 + nrow(fitnessEffects$long.geneNoInt)
## tnn <- tot(nn)
## if(tnn > max) {
## m <- paste("There are", tnn, "genotypes.")
## m <- paste(m, "This is larger than max.")
## m <- paste(m, "Adjust max and rerun if you want")
## stop(m)
## }
## ## With mutator, the ids of genes need not go from 1:n
## vid <- allNamedGenes(fitnessEffects)$GeneNumID
## if(order) {
## f1 <- function(n) {
## lapply(seq.int(n), function(x) permutations(n = n, r = x, v = vid))
## }
## } else {
## f1 <- function(n) {
## lapply(seq.int(n), function(x) combinations(n = n, r = x, v = vid))}
## }
## genotNums <- f1(nn)
## list.of.vectors <- function(y) {
## ## there's got to be a simpler way
## lapply(unlist(lapply(y, function(x) {apply(x, 1, list)}), recursive = FALSE),
## function(m) m[[1]])
## }
## genotNums <- list.of.vectors(genotNums)
## names <- c(fitnessEffects$geneModule$Gene[-1],
## fitnessEffects$long.geneNoInt$Gene)
## genotNames <- unlist(lapply(lapply(genotNums, function(x) names[x]),
## function(z)
## paste(z,
## collapse = if(order){" > "} else {", "} )))
## This ain't the best, as we repeatedly read and convert
## fitnessEffects. If this were slow, prepare C++ function that takes
## vectors of genotypes and uses same fitnessEffects.
if(model %in% c("Bozic", "bozic1", "bozic2") )
prodNeg <- TRUE
else
prodNeg <- FALSE
allf <- vapply(allg$genotNums,
function(x) evalRGenotype(x, fmEffects, FALSE,
prodNeg,
calledBy_),
1.1)
df <- data.frame(Genotype = allg$genotNames, Fitness = allf,
stringsAsFactors = FALSE)
## Why am I doing this? I am not computing the genotype. I test the
## evaluation of the empty genotype in the tests. But its evaluation
## generates warnings that are not needed. And by decree (even in the
## C++) this thing has a fitness of 1, a mutator effect of 1 since
## there are no terms.
if(addwt)
df <- rbind(data.frame(Genotype = "WT", Fitness = 1,
stringsAsFactors = FALSE), df)
if(calledBy_ == "evalGenotype") {
if(prodNeg)
colnames(df)[match("Fitness", colnames(df))] <- "Death_rate"
class(df) <- c("evalAllGenotypes", class(df))
} else if (calledBy_ == "evalGenotypeMut") {
colnames(df)[match("Fitness", colnames(df))] <- "MutatorFactor"
class(df) <- c("evalAllGenotypesMut", class(df))
}
return(df)
}
evalAllGenotypes <- function(fitnessEffects, order = FALSE, max = 256,
addwt = FALSE,
model = "") {
## Must deal with objects from previous, pre flfast, modifications
if(!exists("fitnessLandscape_gene_id", where = fitnessEffects)) {
fitnessEffects$fitnessLandscape_df <- data.frame()
fitnessEffects$fitnessLandscape_gene_id <- data.frame()
}
if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
(nrow(fitnessEffects$fitnessLandscape_df) > 0)) {
warning("Bozic model passing a fitness landscape will not work",
" for now.")
}
evalAllGenotypesORMut(
fmEffects = fitnessEffects,
order = order,
max = max,
addwt = addwt,
model = model,
calledBy_= "evalGenotype"
)
}
generateAllGenotypes <- function(fitnessEffects, order = TRUE, max = 256) {
if(order)
tot <- function(n) {sum(sapply(seq.int(n),
function(x) choose(n, x) * factorial(x)))}
else
tot <- function(n) {2^n}
nn <- nrow(fitnessEffects$geneModule) -1 +
nrow(fitnessEffects$long.geneNoInt) +
nrow(fitnessEffects$fitnessLandscape_gene_id)
tnn <- tot(nn)
if(tnn > max) {
m <- paste("There are", tnn, "genotypes.")
m <- paste(m, "This is larger than max.")
m <- paste(m, "Adjust max and rerun if you want")
stop(m)
}
## With mutator, the ids of genes need not go from 1:n
vid <- allNamedGenes(fitnessEffects)$GeneNumID
if(order) {
f1 <- function(n) {
lapply(seq.int(n), function(x) permutations(n = n, r = x, v = vid))
}
} else {
f1 <- function(n) {
lapply(seq.int(n), function(x) combinations(n = n, r = x, v = vid))}
}
genotNums <- f1(nn)
list.of.vectors <- function(y) {
## there's got to be a simpler way
lapply(unlist(lapply(y, function(x) {apply(x, 1, list)}), recursive = FALSE),
function(m) m[[1]])
}
genotNums <- list.of.vectors(genotNums)
names <- c(fitnessEffects$geneModule$Gene[-1],
fitnessEffects$long.geneNoInt$Gene,
fitnessEffects$fitnessLandscape_gene_id$Gene)
genotNames <- unlist(lapply(lapply(genotNums, function(x) names[x]),
function(z)
paste(z,
collapse = if(order){" > "} else {", "} )))
## This ain't the best, as we repeatedly read and convert
## fitnessEffects. If this were slow, prepare C++ function that takes
## vectors of genotypes and uses same fitnessEffects.
return(list(genotNums = genotNums,
genotNames = genotNames
))
}
evalAllGenotypesFitAndMut <- function(fitnessEffects, mutatorEffects,
order = FALSE, max = 256,
addwt = FALSE,
model = "" ){
## minimal = FALSE) {
## if(!minimal)
allg <- generateAllGenotypes(fitnessEffects = fitnessEffects,
order = order, max = max)
## else
## allg <- generateAllGenotypes_minimal(fitnessEffects = fitnessEffects,
## max = max)
if(model %in% c("Bozic", "bozic1", "bozic2") ) {
prodNeg <- TRUE
} else {
prodNeg <- FALSE
}
## Must deal with objects from previous, pre flfast, modifications
if(!exists("fitnessLandscape_gene_id", where = fitnessEffects)) {
fitnessEffects$fitnessLandscape_df <- data.frame()
fitnessEffects$fitnessLandscape_gene_id <- data.frame()
}
if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
(nrow(fitnessEffects$fitnessLandscape_df) > 0)) {
warning("Bozic model passing a fitness landscape will not work",
" for now.")
}
full2mutator_ <- matchGeneIDs(mutatorEffects,
fitnessEffects)$Reduced
allf <- t(vapply(allg$genotNums,
function(x) evalRGenotypeAndMut(x,
rFE = fitnessEffects,
muEF= mutatorEffects,
full2mutator_ = full2mutator_,
verbose = FALSE,
prodNeg = prodNeg),
c(1.1, 2.2)))
df <- data.frame(Genotype = allg$genotNames, Fitness = allf[, 1],
MutatorFactor = allf[, 2],
stringsAsFactors = FALSE)
if(addwt)
df <- rbind(data.frame(Genotype = "WT", Fitness = 1,
MutatorFactor = 1,
stringsAsFactors = FALSE), df)
if(prodNeg)
colnames(df)[match("Fitness", colnames(df))] <- "Death_rate"
class(df) <- c("evalAllGenotypesFitAndMut", class(df))
return(df)
}
## evalAllGenotypes <- function(fitnessEffects, order = TRUE, max = 256,
## addwt = FALSE,
## model = "") {
## if(order)
## tot <- function(n) {sum(sapply(seq.int(n),
## function(x) choose(n, x) * factorial(x)))}
## else
## tot <- function(n) {2^n}
## nn <- nrow(fitnessEffects$geneModule) -1 + nrow(fitnessEffects$long.geneNoInt)
## tnn <- tot(nn)
## if(tnn > max) {
## m <- paste("There are", tnn, "genotypes.")
## m <- paste(m, "This is larger than max.")
## m <- paste(m, "Adjust max and rerun if you want")
## stop(m)
## }
## if(order) {
## f1 <- function(n) {
## lapply(seq.int(n), function(x) permutations(n = n, r = x))
## }
## } else {
## f1 <- function(n) {
## lapply(seq.int(n), function(x) combinations(n = n, r = x))}
## }
## genotNums <- f1(nn)
## list.of.vectors <- function(y) {
## ## there's got to be a simpler way
## lapply(unlist(lapply(y, function(x) {apply(x, 1, list)}), recursive = FALSE),
## function(m) m[[1]])
## }
## genotNums <- list.of.vectors(genotNums)
## names <- c(fitnessEffects$geneModule$Gene[-1],
## fitnessEffects$long.geneNoInt$Gene)
## genotNames <- unlist(lapply(lapply(genotNums, function(x) names[x]),
## function(z)
## paste(z,
## collapse = if(order){" > "} else {", "} )))
## ## This ain't the best, as we repeatedly read and convert
## ## fitnessEffects. If this were slow, prepare C++ function that takes
## ## vectors of genotypes and uses same fitnessEffects.
## if(model %in% c("Bozic", "bozic1", "bozic2") )
## prodNeg <- TRUE
## else
## prodNeg <- FALSE
## allf <- vapply(genotNums,
## function(x) evalRGenotype(x, fitnessEffects, FALSE, prodNeg),
## 1.1)
## df <- data.frame(Genotype = genotNames, Fitness = allf,
## stringsAsFactors = FALSE)
## if(addwt)
## df <- rbind(data.frame(Genotype = "WT", Fitness = 1,
## stringsAsFactors = FALSE), df)
## if(prodNeg)
## colnames(df)[match("Fitness", colnames(df))] <- "Death_rate"
## return(df)
## }
fitnessEffectsToIgraph <- function(rT, epistasis, orderEffects) {
df0 <- df1 <- df2 <- data.frame()
if(!is.null(rT)) {
df0 <- rT[, c("parent", "child", "typeDep")]
}
if(!is.null(epistasis)) {
df1 <- to.pairs.modules(epistasis, ":")
}
if(!is.null(orderEffects)) {
df2 <- to.pairs.modules(orderEffects, ">")
}
df <- rbind(df0, df1, df2)
## for special case of simple epi effects
if(nrow(df) == 0) return(NULL)
g1 <- graph.data.frame(df)
E(g1)$color <- "black"
E(g1)$color[E(g1)$typeDep == "SM"] <- "blue"
E(g1)$color[E(g1)$typeDep == "XMPN"] <- "red"
E(g1)$color[E(g1)$typeDep == "epistasis"] <- "orange"
E(g1)$color[E(g1)$typeDep == "orderEffect"] <- "orange"
E(g1)$lty <- 1
E(g1)$lty[E(g1)$typeDep == "epistasis"] <- 2
E(g1)$lty[E(g1)$typeDep == "orderEffect"] <- 3
E(g1)$arrow.mode <- 2
E(g1)$arrow.mode[E(g1)$typeDep == "epistasis"] <- 0
E(g1)$arrow.mode[E(g1)$typeDep == "orderEffect"] <- 2
return(g1)
}
plot.fitnessEffects <- function(x, type = "graphNEL",
layout = NULL,
expandModules = FALSE,
autofit = FALSE,
scale_char = ifelse(type == "graphNEL", 1/10, 5),
return_g = FALSE,
lwdf = 1, ## multiplier for lwd for graphNEL
...) {
## some other layouts I find OK
## layout.circle
## layout.reingold.tilford if really a tree
## o.w. it will use the default
g <- x$graph
if(is.null(g))
stop("This fitnessEffects object can not be ploted this way.",
" It is probably one with fitness landscape specification, ",
" so you might want to plot the fitness landscape instead.")
if(type == "igraph") {
if(expandModules && (!x$gMOneToOne)) {
## vlabels <- fe$geneToModule[vertex.attributes(g)$name]
vlabels <- x$geneToModule[V(g)$name]
V(g)$label <- vlabels
}
if(autofit) {
plot(0, type = "n", axes = FALSE, ann = FALSE)
## ideas from http://stackoverflow.com/questions/14472079/match-vertex-size-to-label-size-in-igraph
## vsize <- (strwidth(V(g)$label) + strwidth("oo")) * 200
## but this is a kludge.
vsize <- (nchar(V(g)$label) + 1) * scale_char
plot.igraph(g, vertex.size = vsize, vertex.shape = "rectangle",
layout = layout)
} else {
plot.igraph(g, layout = layout)
}
if(return_g) return(g)
}
else if (type == "graphNEL") {
g1 <- igraph.to.graphNEL(g)
c1 <- unlist(lapply(edgeData(g1), function(x) x$color))
names(c1) <- sub("|", "~", names(c1), fixed = TRUE)
s1 <- unlist(lapply(edgeData(g1), function(x) x$lty))
names(s1) <- names(c1)
a1 <- unlist(lapply(edgeData(g1), function(x) max(x$arrow.mode - 1, 0)))
names(a1) <- names(c1)
lwd <- s1
lwd[lwd == 2] <- 2 ## o.w. too thin
lwd[lwd == 3] <- 2 ## o.w. too thin
lwd <- lwd * lwdf
nAttrs <- list()
if(expandModules && (!x$gMOneToOne)) {
nnodes <- x$geneToModule[nodes(g1)]
names(nnodes) <- nodes(g1)
nAttrs$label <- nnodes
}
if(autofit) {
nAttrs$width <- (nchar(nAttrs$label) + 1) * scale_char
names(nAttrs$width) <- names(nAttrs$label)
plot(g1, edgeAttrs = list(arrowsize = a1, lty = s1, lwd = lwd,
color = c1), attrs=list(node=list(shape = "rectangle")),
nodeAttrs = nAttrs)
} else {
plot(g1, edgeAttrs = list(arrowsize = a1, lty = s1, lwd = lwd,
color = c1),
nodeAttrs = nAttrs)
}
if(return_g) return(g1)
} else {
stop("plot type not recognized")
}
}
## plot.fitnessEffects <- plotFitnessEffects
## FIXME: in the help: say we cannot properly show 3- or higher order
## interactions.
nr_oncoSimul.internal <- function(rFE,
birth,
death,
mu,
initSize,
sampleEvery,
detectionSize,
finalTime,
initSize_species,
initSize_iter,
seed,
verbosity,
speciesFS,
ratioForce,
typeFitness,
max.memory,
mutationPropGrowth,
initMutant,
max.wall.time,
keepEvery,
K,
detectionDrivers,
onlyCancer,
errorHitWallTime,
max.num.tries,
errorHitMaxTries,
minDetectDrvCloneSz,
extraTime,
keepPhylog,
detectionProb,
AND_DrvProbExit,
MMUEF = NULL, ## avoid partial matching, and set default
fixation = NULL ## avoid partial matching
) {
default_min_successive_fixation <- 50 ## yes, set at this for now
if(!inherits(rFE, "fitnessEffects"))
stop(paste("rFE must be an object of class fitnessEffects",
"as created, for instance, with function",
"allFitnessEffects"))
if(countGenesFe(rFE) < 2) {
stop("There must be at least two genes (loci) in the fitness effects.",
"If you only care about a degenerate case with just one,",
"you can enter a second gene (locus)",
"with fitness effect of zero.")
}
if( (length(mu) == 1) && !(is.null(names(mu)))) {
stop("A length 1 mutation, but named. ",
"This is ambiguous. ",
"If you want per-gene mutation rates, each gene",
"must have its entry in the mu vector. ",
"(And regardless of per-gene mutation rates ",
" there must be at least two gene/loci).")
}
## Must deal with objects from previous, pre flfast, modifications
## Could move this way down to the bottom, right before
## .Call
if(!exists("fitnessLandscape_gene_id", where = rFE)) {
rFE$fitnessLandscape_df <- data.frame()
rFE$fitnessLandscape_gene_id <- data.frame()
}
namedGenes <- allNamedGenes(rFE)
if( length(mu) > 1) {
if(is.null(names(mu)))
stop("When using per-gene mutation rates the ",
"mu vector must be named ",
"(and if you have noIntGenes, those must have names).")
if(length(mu) != countGenesFe(rFE))
stop("When using per-gene mutation rates, ",
"there must be the same number of genes in the ",
"mu vector and the fitness effects object.")
if(!identical(sort(namedGenes$Gene),
sort(names(mu))))
stop("When using per-gene mutation rates, ",
"names of genes must match those in the ",
"restriction table.")
mu <- mu[order(match(names(mu), namedGenes$Gene))]
## Hyperparanoid check. Should never, ever, happen.
if(!identical(names(mu), namedGenes$Gene))
stop("Names of re-ordered mu do not match names of genes")
minmu <- 1e-40
if(any(mu <= minmu))
stop(paste("At least one per-gene mutation rate is negative",
"or less than", minmu,". Remember that the per-base",
"mutation rate in the human genome is about 1e-11 to 1e-9."))
}
if(!is.null(initMutant)) {
initMutantString <- initMutant
if(length(grep(">", initMutant))) {
initMutant <- nice.vector.eo(initMutant, ">")
} else if(length(grep(",", initMutant))) {
initMutant <- nice.vector.eo(initMutant, ",")
}
initMutant <- getGeneIDNum(rFE$geneModule,
rFE$long.geneNoInt,
rFE$fitnessLandscape_gene_id,
initMutant,
FALSE)
if(length(initMutant) >= countGenesFe(rFE)) {
stop("For initMutant you passed as many, or more genes, mutated ",
"than the number of genes in the genotype (fitness effects).")
}
} else {
initMutant <- vector(mode = "integer")
initMutantString <- ""
}
## these are never user options
## if(initSize_species < 10) {
## warning("initSize_species too small?")
## }
## if(initSize_iter < 100) {
## warning("initSize_iter too small?")
## }
if(typeFitness %in% c("bozic1", "bozic2")) {
if(nrow(rFE$fitnessLandscape_df) > 0)
warning("Bozic model passing a fitness landscape will not work",
" for now.")
## FIXME: bozic and fitness landscape
## the issue is that in the C++ code we directly do
## s = birth rate - 1
## but we would need something different
## Can be done going through epistasis, etc.
thesh <- unlist(lapply(rFE$long.rt, function(x) x$sh))
## thes <- unlist(lapply(rFE$long.rt, function(x) x$s))
thes <- unlist(c(lapply(rFE$long.rt, function(x) x$s),
lapply(rFE$long.epistasis, function(x) x$s),
lapply(rFE$long.orderEffects, function(x) x$s),
rFE$long.geneNoInt$s))
if(any(thes > 1 )) {
m <- paste("You are using a Bozic model with",
"the new restriction specification, and you have",
"at least one s > 1.",
"But that is the same as setting that to 1:",
"we obviously cannot allow negative death rates,",
"nor problems derived from multiplying odd or even",
"numbers of negative numbers.",
"You can set those > 1 to 1, but then you will",
"eventually get the simulations to abort when the",
"death rate becomes 0.")
stop(m)
}
if(any( (thesh <= -1) & (thesh > -Inf))) {
m <- paste("You are using a Bozic model with",
"the new restriction specification, and you have",
"at least one sh <= -1. Maybe you mean -Inf?")
warning(m)
}
if(any(thes == 1 )) {
m <- paste("You are using a Bozic model with",
"the new restriction specification, and you have",
"at least one s of 1. If that gene is mutated, ",
"this will lead to a death rate of 0",
"and the simulations will abort when you get a non finite",
"value.")
warning(m)
}
}
if(!is.null(MMUEF)) {
if(!inherits(MMUEF, "mutatorEffects"))
stop("muEF must be a mutatorEffects object")
full2mutator_ <- matchGeneIDs(MMUEF, rFE)
if(any(is.na(full2mutator_$Full)))
stop("Genes in mutatorEffects not present in fitnessEffects")
if(any(is.na(full2mutator_)))
stop("full2mutator failed for unknown reasons")
full2mutator_ <- full2mutator_$Reduced
} else {
full2mutator_ <- vector(mode = "numeric", length = 0)
## muEF <- emptyFitnessEffects()
}
dpr <- detectionProbCheckParse(detectionProb, initSize, verbosity)
## if( !is.null(cPDetect) && (sum(!is.null(p2), !is.null(n2)) >= 1 ))
## stop("Specify only cPDetect xor both of p2 and n2")
## if( (is.null(p2) + is.null(n2)) == 1 )
## stop("If you pass one of n2 or p2, you must also pass the other. ",
## "Otherwise, we would not know what to do.")
## stopifnot(PDBaseline >= 0)
## stopifnot(n2 > PDBaseline)
## stopifnot(p2 < 1)
## stopifnot(p2 > 0)
## if( is.null(cPDetect) ) cPDetect <- -9
## if( is.null(p2)) p2 <- 9
## if( is.null(n2)) n2 <- -9
## call <- match.call()
## Process the fixed list, if any
if(!is_null_na(fixation)) {
ng <- namedGenes
## rownames(ng) <- namedGenes[, "Gene"]
## Proposed extension to have exact matching of genotypes
ng <- rbind(
data.frame(Gene = "_", GeneNumID = -9, stringsAsFactors = FALSE),
ng)
rownames(ng) <- ng[, "Gene"]
## FIXME
## Later, accept a last argument, called tolerance.
## If not present, set to 0
## and then, at at the head of fixation_list below
if(is.list(fixation)) {
if(
(is.null(fixation[["fixation_tolerance"]])) ||
(is.na(fixation[["fixation_tolerance"]]))) {
fixation_tolerance <- 0
} else {
fixation_tolerance <- as.numeric(fixation[["fixation_tolerance"]])
fixation <- fixation[-which(names(fixation) == "fixation_tolerance")]
}
if(
(is.null(fixation[["min_successive_fixation"]])) ||
(is.na(fixation[["min_successive_fixation"]]))) {
min_successive_fixation <- default_min_successive_fixation
} else {
min_successive_fixation <- as.integer(fixation[["min_successive_fixation"]])
fixation <- fixation[-which(names(fixation) == "min_successive_fixation")]
}
if(
(is.null(fixation[["fixation_min_size"]])) ||
(is.na(fixation[["fixation_min_size"]]))) {
fixation_min_size <- 0
} else {
fixation_min_size <- as.integer(fixation[["fixation_min_size"]])
fixation <- fixation[-which(names(fixation) == "fixation_min_size")]
}
} else {
if(is_null_na(fixation["fixation_tolerance"])) {
fixation_tolerance <- 0
} else {
fixation_tolerance <- as.numeric(fixation["fixation_tolerance"])
fixation <- fixation[-which(names(fixation) == "fixation_tolerance")]
}
if(is_null_na(fixation["min_successive_fixation"])) {
min_successive_fixation <- default_min_successive_fixation
} else {
min_successive_fixation <- as.integer(fixation["min_successive_fixation"])
fixation <- fixation[-which(names(fixation) == "min_successive_fixation")]
}
if(is_null_na(fixation["fixation_min_size"])) {
fixation_min_size <- 0
} else {
fixation_min_size <- as.integer(fixation["fixation_min_size"])
fixation <- fixation[-which(names(fixation) == "fixation_min_size")]
}
}
if( (fixation_tolerance > 1) || (fixation_tolerance < 0) )
stop("Impossible range for fixation tolerance")
if( (min_successive_fixation < 0) )
stop("Impossible range for min_successive_fixation")
if( (fixation_min_size < 0) )
stop("Impossible range for fixation_min_size")
## Usual genotype specification and might allow ordered vectors
## in the future
fixation_b <- lapply(fixation, nice.vector.eo, sep = ",")
ulf <- unlist(fixation_b)
if(any(ulf == ">"))
stop("Order effects not allowed (yet?) in fixation.")
ulfg <- ng[ulf, 1]
if(any(is.na(ulfg)))
stop(paste("The 'fixation' list contains genes that are not present",
" in the fitness effects."))
## Sorting here is crucial!!
fixation_list <- lapply(fixation_b, function(x) sort(ng[x, 2]))
fixation_list <- list(fixation_list = fixation_list,
fixation_tolerance = fixation_tolerance,
min_successive_fixation = min_successive_fixation,
fixation_min_size = fixation_min_size)
} else {
fixation_list <- list()
}
return(c(
nr_BNB_Algo5(rFE = rFE,
mu_ = mu,
death = death,
initSize = initSize,
sampleEvery = sampleEvery,
detectionSize = detectionSize,
finalTime = finalTime,
initSp = initSize_species,
initIt = initSize_iter,
seed = seed,
verbosity = verbosity,
speciesFS = speciesFS,
ratioForce = ratioForce,
typeFitness_ = typeFitness,
maxram = max.memory,
mutationPropGrowth = as.integer(mutationPropGrowth),
initMutant_ = initMutant,
maxWallTime = max.wall.time,
keepEvery = keepEvery,
K = K,
detectionDrivers = detectionDrivers,
onlyCancer = onlyCancer,
errorHitWallTime = errorHitWallTime,
maxNumTries = max.num.tries,
errorHitMaxTries = errorHitMaxTries,
minDetectDrvCloneSz = minDetectDrvCloneSz,
extraTime = extraTime,
keepPhylog = keepPhylog,
MMUEF = MMUEF,
full2mutator_ = full2mutator_,
n2 = dpr["n2"],
p2 = dpr["p2"],
PDBaseline = dpr["PDBaseline"],
cPDetect_i= dpr["cPDetect"],
checkSizePEvery = dpr["checkSizePEvery"],
AND_DrvProbExit = AND_DrvProbExit,
fixation_list = fixation_list),
Drivers = list(rFE$drv), ## but when doing pops, these will be repeated
geneNames = list(names(getNamesID(rFE))),
InitMutant = initMutantString
))
}
countGenesFe <- function(fe) {
## recall geneModule has Root always
## We want to be able to use objects that did not have
## the fitness landscape component
if(exists("fitnessLandscape_gene_id", where = fe))
return(nrow(fe$geneModule) + nrow(fe$long.geneNoInt) +
nrow(fe$fitnessLandscape_gene_id) - 1)
else
return(nrow(fe$geneModule) + nrow(fe$long.geneNoInt) - 1)
}
allNamedGenes <- function(fe){
## Returns a data frame with genes and their names and verifies all
## genes have names.
## Root should always be first, but just in case
## avoid assuming it
## This does is not a good idea as it assumes the user did not use
## "keepInput = FALSE".
## lni <- length(fe$noIntGenes)
## ## FIXME:test
## if((lni > 0) &&
## (is.null(names(fe$noIntGenes))))
## stop("When using per-gene mutation rates the ",
## "no interaction genes must be named ",
## "(i.e., the noIntGenes vector must have names).")
## accommodate objects w/o fitnessLandscape
if(!is.null(fe$fitnessLandscape) && nrow(fe$fitnessLandscape)) {
gn <-
gtools::mixedsort(
colnames(fe$fitnessLandscape)[-ncol(fe$fitnessLandscape)])
v1 <- data.frame(Gene = gn, GeneNumID = seq.int(length(gn)),
stringsAsFactors = FALSE)
} else {
v1 <- fe$geneModule[, c("Gene", "GeneNumID")]
if(nrow(fe$long.geneNoInt)) {
v1 <- rbind(v1,
fe$long.geneNoInt[, c("Gene", "GeneNumID")])
}
if(any(v1[, "Gene"] == "Root"))
v1 <- v1[-which(v1[, "Gene"] == "Root"), ]
}
rownames(v1) <- NULL
if(any(v1$Gene == "WT")) {
warning("A gene is named WT. You can expect problems ",
"because we use WT to denote the wildtype ",
"genotype. You might want to change it.")
}
return(v1)
}
get.gene.counts <- function(x) {
# , timeSample = "last",
# typeSample = "whole") {
## From get.mut.vector. Used for now just for testing
timeSample <- "last"
the.time <- get.the.time.for.sample(x, timeSample, NULL)
if(the.time < 0) {
counts <- rep(0, length(x$geneNames))
names(counts) <- x$geneNames
freq <- rep(NA, length(x$geneNames))
names(freq) <- x$geneNames
return(list(counts = counts,
freq = freq,
popSize = 0))
}
pop <- x$pops.by.time[the.time, -1]
if(all(pop == 0)) {
stop("You found a bug: this should never happen")
}
## if(typeSample %in% c("wholeTumor", "whole")) {
popSize <- x$PerSampleStats[the.time, 1]
counts <- as.vector(tcrossprod(pop, x$Genotypes))
names(counts) <- x$geneNames
return(list(counts = counts,
freq = counts/popSize,
popSize = popSize))
## return( (tcrossprod(pop,
## x$Genotypes)/popSize) )
## } else if (typeSample %in% c("singleCell", "single")) {
## return(x$Genotypes[, sample(seq_along(pop), 1, prob = pop)])
## } else {
## stop("Unknown typeSample option")
## }
}
geneCounts <- function(x) {
if(inherits(x, "oncosimulpop")) {
return( do.call(rbind,
lapply(x,
function(z)
get.gene.counts(z)$counts))
)
} else {
return( get.gene.counts(x)$counts)
}
}
## geneNumIDReset <- function(x, ref){
## ## Set GeneNumID of a fitnessEffect object, x, using ref as the
## ## reference fitness effect object.
## ## Check also if all in x are in ref.
## gg <- allNamedGenes(ref)
## gnid <- gg$GeneNumID
## names(gnid) <- gg$Gene
## ## FIXME: this later and conditional on what is in thee
## gnid <- c("Root" = 0, gnid)
## if(!all(x$geneModule$Gene %in% names(gnid) ))
## stop("Some genes not in reference fitnessEffects (rebasing geneModule)")
## x$geneModule$GeneNumID <- gnid[geneModule$Gene]
## ## and then go over all the lists in the x object.
## if(nrow(x$long.geneNoInt)) {
## ## now, mapping for the noInt if this is mutator
## if(!all(x$long.geneNoInt$Gene %in% names(gnid) ))
## stop("Some genes not in reference fitnessEffects (rebasing geneNoInt)")
## x$long.geneNoInt$GeneNumID <- gnid[long.geneNoInt$Gene]
## }
## }
matchGeneIDs <- function(x, refFE) {
n1 <- allNamedGenes(x)
n2 <- allNamedGenes(refFE)
colnames(n1)[2] <- "Reduced"
colnames(n2)[2] <- "Full"
## Non standard evaluation thing and a note being generated in check.
## See also http://stackoverflow.com/questions/33299845/when-writing-functions-for-an-r-package-should-i-avoid-non-standard-evaluation
## But does not work well with replace. So use the NULL trick
Reduced <- NULL
dplyr::full_join(n2, n1, by = "Gene") %>%
mutate(Reduced = replace(Reduced, is.na(Reduced), -9))
}
detectionProbCheckParse <- function(x, initSize, verbosity) {
default_p2 <- 0.1
default_n2 <- 2 * initSize
default_PDBaseline <- 1.2 * initSize
default_checkSizePEvery <- 20
## No default cPDetect. That is done from p2 and n2 in C++.
if(is_null_na(x)) {
## if((length(x) == 1) && (is.na(x))) {
y <- vector()
y["cPDetect"] <- -9
y["p2"] <- 9
y["n2"] <- -9
y["PDBaseline"] <- -9
y["checkSizePEvery"] <- Inf
return(y)
} else if((length(x) == 1) && (x == "default")) {
## Populate with defaults
y <- vector()
y["p2"] <- default_p2
y["n2"] <- default_n2
y["PDBaseline"] <- default_PDBaseline
y["checkSizePEvery"] <- default_checkSizePEvery
x <- y
}
if(length(x) >= 1) {
if( !all(names(x) %in% c("cPDetect", "p2", "n2", "PDBaseline", "checkSizePEvery")))
stop("Names of some components of detectionProb are not recognized.",
" Check for possible typos.")
}
## This ain't conditional. If not given, always check
if( !is_null_na(x["cPDetect"]) && (sum(!is_null_na(x["p2"]), !is_null_na(x["n2"])) >= 1 ))
stop("Specify only cPDetect xor both of p2 and n2")
if( (is_null_na(x["p2"]) + is_null_na(x["n2"])) == 1 )
stop("If you pass one of n2 or p2, you must also pass the other. ",
"Otherwise, we would not know what to do.")
if(is_null_na(x["PDBaseline"])) {
x["PDBaseline"] <- default_PDBaseline
if(verbosity > -1)
message("\n PDBaseline set to default value of ", default_PDBaseline, "\n")
}
if(is_null_na(x["checkSizePEvery"])) {
x["checkSizePEvery"] <- default_checkSizePEvery
if(verbosity > -1)
message("\n checkSizePEvery set to default value of ",
default_checkSizePEvery, "\n")
}
## If we get here, we checked consistency of whether cPDetect or p2/n2.
## Fill up with defaults the missing values
if(is_null_na(x["cPDetect"])) {
if(is_null_na(x["p2"])) {
if(!is_null_na(x["n2"])) stop("Eh? no p2 but n2? Bug")
x["n2"] <- default_n2
x["p2"] <- default_p2
if(verbosity > -1)
message("\n n2, p2 set to default values of ",
default_n2, ", ", default_p2, "\n")
}
}
if(x["checkSizePEvery"] <= 0)
stop("checkSizePEvery <= 0")
if(x["PDBaseline"] <= 0)
stop("PDBaseline <= 0")
if(!is_null_na(x["n2"])) {
if(x["n2"] <= x["PDBaseline"])
stop("n2 <= PDBaseline")
if(x["p2"] >= 1) stop("p2 >= 1")
if(x["p2"] <= 0) stop("p2 <= 0")
x["cPDetect"] <- -9
} else { ## only if x["cPDetect"] is not NA
if(is_null_na(x["cPDetect"])) stop("eh? you found a bug")## paranoia
x["n2"] <- -9
x["p2"] <- -9
if(verbosity > -1)
message("\n Using user-specified cPDetect as n2, p2 not given.\n")
}
return(x)
}
sampledGenotypes <- function(y, genes = NULL) {
## From a popSample object, or a matrix for that matter,
## show the sampled genotypes and their frequencies
if(!is.null(genes)) {
cols <- which(colnames(y) %in% genes )
y <- y[, cols]
}
## nn <- colnames(y)
genotlabel <- function(u, nn = colnames(y)) {
if(any(is.na(u))) return(NA)
else {
return(paste(sort(nn[as.logical(u)]), collapse = ", "))
}
}
## df <- data.frame(table(
## apply(y, 1, function(z) paste(sort(nn[as.logical(z)]), collapse = ", ") )
## ))
df <- data.frame(table(
apply(y, 1, genotlabel),
useNA = "ifany"),
stringsAsFactors = TRUE)
gn <- as.character(df[, 1])
gn[gn == ""] <- "WT"
df <- data.frame(Genotype = gn, Freq = df[, 2], stringsAsFactors = FALSE)
attributes(df)$ShannonI <- shannonI(na.omit(df)$Freq)
class(df) <- c("sampledGenotypes", class(df))
return(df)
}
print.sampledGenotypes <- function(x, ...) {
print.data.frame(x, ...)
cat("\n Shannon's diversity (entropy) of sampled genotypes: ")
cat(attributes(x)$ShannonI, "\n")
}
list_g_matches_fixed <- function(x, y) {
## Internal function, for testing the fixation output.
## x and y are vectors
## x is the set of output genotypes, y the set of fixed
## genotypes/subset of genotypes.
## Yes, this function has tests too in test.fixation.R
## All genotypes in x satisfy that they are supersets of at least one
## in y? That is true if, for every element in x, at least one y in
## that x.
if(is.list(y)) y <- unlist(y)
y.nice <- lapply(y, nice.vector.eo, sep = ",")
x.nice <- lapply(x, nice.vector.eo, sep = ",")
fu <- function(u, y.nice)
any(unlist(lapply(y.nice, function(z) all(z %in% u))))
return(all(unlist(lapply(x.nice, fu, y.nice))))
}
## emptyFitnessEffects <- function() {
## list(long.rt = list(),
## long.epistasis = list(),
## long.orderEffects = list(),
## long.geneNoInt = list(),
## geneModule = list(),
## gMOneToOne = TRUE,
## geneToModule = list(),
## graph = NULL,
## drv = vector(mode = "integer", length = 0)
## )
## }
### Later, for all the effects, we will do some kind of dplyr match?
### a, b, c, in fitness, only a, c in mut.
### fitness table for a,b,c
### each row name transformed removing b (so leaving only present)
### each row transformed matched to row in mut table.
## t1 <- data.frame(v1 = c("a,b", "a,c", "b"), v2 = c("b", "c", "b"), v3 = 1:3, stringsAsFactors = FALSE)
## t2 <- data.frame(v2 = c("b", "c"), v4 = c(11, 12), stringsAsFactors = FALSE)
## full_join(t1, t2, by = "v2")
## FIXME
## new exit code
## check:
## baseline < n2
## p2 < 1
## baseline defaults to init size + .1
## use c < -9, n2 < -9, p2 < -9 for no values
## and check one of c or p2 and n2 are valid if using exit
Try the OncoSimulR package in your browser
Any scripts or data that you put into this service are public.
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.