| 6fe809cf |
################################################################################
## FUNCTIONS FOR STABILITY ANALYSIS OF THE ONCOGENETIC TREES MIXTURE MODEL ##
################################################################################
## Function for calculating the P value of a given similarity measure.
## What is the probability for obtaining the same or smaller
## value in a random vector of similarity values.
Pval.dist <- function(dist.val, random.vals) {
return((sum(random.vals <= dist.val) + 1) /(length(random.vals) + 1))
}
######################################################################
## Function for calculating the Kullback-Leibler divergence between
## two discrete probability distributions p and q.
kullback.leibler <- function(p, q) {
if(length(p) != length(q))
stop("Error: The distribution vectors have different lengths!")
return(sum(p * log(p / q)))
}
######################################################################
## Function for calculating the L1 distance between
## two discrete probability distributions p and q.
L1.dist <- function(p, q) {
if(length(p) != length(q))
stop("Error: The distribution vectors have different lengths!")
return(sum(abs(p - q)))
}
######################################################################
## Function for calculating the cosine distance between
## two discrete probability distributions p and q.
cosin.dist <- function(p, q) {
if(length(p) != length(q))
stop("Error: The distribution vectors have different lengths!")
return(1 - (sum (p * q) / (L2.norm(p) * L2.norm(q)) ))
}
######################################################################
## Function for calculating the L2 norm of a given vector x.
L2.norm <- function(x) {
return(sqrt(sum(x * x)))
}
######################################################################
## Function for calculating the Euclidian distance between
## two vectors x and y.
euclidian.dist <- function(x, y) {
if(length(x) != length(y))
stop("Error: The vectors have different lengths!")
return(sqrt(drop((x - y) %*% (x - y))))
}
######################################################################
## Function for calculating the rank-correlation distance between
## two vectors x and y.
rank.cor.dist <- function(x, y) {
if(length(x) != length(y))
stop("Error: The vectors have different length!")
return(1 - rcorr(x, y, type = "spearman")[[1]][1, 2])
}
######################################################################
## Function that implements a similarity measure for comparing the topologies
## of the trees of two mixture models mixture1 and mixture2.
## It returns a value that uses the number of different edges for caracterizing
## the difference between the models.
## For the comparisson it is necessary that the two models have the same
## number of tree components build on the same number of genetic events.
## It is assumed that the mixtures have at least two components (the first one is the noise).
comp.models <- function(mixture1, mixture2) {
if((class(mixture1) != "RtreemixModel") || (class(mixture2) != "RtreemixModel"))
stop("The specified mixture models should be of type RtreemixModel!")
if(numTrees(mixture1) != numTrees(mixture2))
stop("The two specified mixture models don't have the same number of trees!")
if(eventsNum(mixture1) != eventsNum(mixture2))
stop("The two specified mixture models don't have the same number of events!")
no.trees <- numTrees(mixture1)
if (no.trees == 1)
stop("The mixture models should have at least two tree components, where the first one is the noise component!")
no.events <- eventsNum(mixture1)
true.order <- list() ## list specifying the edges in the components of mixture1
fit.order <- list() ## list specifying the edges in the components of mixture2
## Build the vectors that characterize the edges of the components of the mixture models.
## The noise components are ignored since they always have the same (star) topology.
for(k in 2:no.trees) {
true.vec <- character(0)
fit.vec <- character(0)
for(l in 2:no.events) {
ch <- names(which(sapply(edges(getTree(mixture1, k)), function(x, el) {el %in% x}, nodes(getTree(mixture1, k))[l])))
true.vec <- c(true.vec, ifelse((identical(ch, character(0)) || is.null(ch)), "e", ch))
ch <- names(which(sapply(edges(getTree(mixture2, k)), function(x, el) {el %in% x}, nodes(getTree(mixture2, k))[l])))
fit.vec <- c(fit.vec, ifelse((identical(ch, character(0)) || is.null(ch)), "e", ch))
}
names(true.vec) <- nodes(getTree(mixture1, k))[-1]
true.order <- c(true.order, list(true.vec))
names(fit.vec) <- nodes(getTree(mixture2, k))[-1]
fit.order <- c(fit.order, list(fit.vec))
}
## Build the comparisson matrix.
## The (i, j) element is the number of distinct edges of the
## (i + 1)-th and (j + 1)-th component of the models
## mixture1 and mixture2 respectively.
comp.mat <- matrix(nrow = no.trees - 1, ncol = no.trees - 1)
for(i in 1:(no.trees - 1))
for(j in 1:(no.trees - 1)) {
comp.mat[i, j] <- sum(true.order[[i]] != fit.order[[j]])
}
## Form the similarity pairs and calculate the similarity of the models
## as a sum of the number of different edges of the trees in the similarity pairs:
if(no.trees > 2) {
rc <- which(comp.mat == min(comp.mat), arr.ind = TRUE)
diff.sum <- min(comp.mat)
row.index <- c(rc[1, 1]) ## get the row index of the minimum
col.index <- c(rc[1, 2]) ## get the column index of the minimum
for(m in 1:(nrow(comp.mat) - 1)) {
rc <- which(comp.mat == min(comp.mat[-(row.index), -(col.index)]), arr.ind = TRUE)
diff.sum <- diff.sum + min(comp.mat[-(row.index), -(col.index)])
row.index <- c(row.index, rc[1, 1])
col.index <- c(col.index, rc[1, 2])
}
} else {
diff.sum <- comp.mat[1, 1]
}
## Return a result between 0 and 1.
return(diff.sum/((no.trees - 1) * (no.events - 1)))
}
######################################################################
## Function that assignes to each node the level at which that node is
## in a specific tree (tree.num) of the trees mixture model.
## The start.val is the maximum depth of the tree with which the tree
## specified with tree.num will be compared.
get.tree.levels <- function(mixture, tree.num, start.val) {
root <- nodes(getTree(mixture, tree.num))[1]
levels <- acc(getTree(mixture, tree.num), root)
## vec <- rep(max(levels[[1]]), eventsNum(mixture) - 1)
vec <- rep(start.val, eventsNum(mixture) - 1)
names(vec) <- nodes(getTree(mixture, tree.num))[-1]
vec[names(levels[[1]])] <- levels[[1]]
return(vec)
}
######################################################################
## Function that implements a similarity measure for comparing the topologies
## of the trees of two mixture models mixture1 and mixture2.
## It returns a value that uses the number of different edges and the depth at
## which the events are, for caracterizing the difference between the models.
## This measure is more application specific, since the depth at which the
## events are in a tree is very important for disease progression.
## For the comparisson it is necessary that the two models have the same
## number of tree components build on the same number of genetic events.
## It is assumed that the mixtures have at least two components (the first one is the noise).
comp.models.levels <- function(mixture1, mixture2) {
if((class(mixture1) != "RtreemixModel") || (class(mixture2) != "RtreemixModel"))
stop("The specified mixture models should be of type RtreemixModel!")
if(numTrees(mixture1) != numTrees(mixture2))
stop("The two specified mixture models don't have the same number of trees!")
if(eventsNum(mixture1) != eventsNum(mixture2))
stop("The two specified mixture models don't have the same number of events!")
no.trees <- numTrees(mixture1)
if (no.trees == 1)
stop("The mixture models should have at least two tree components, where the first one is the noise component!")
no.events <- eventsNum(mixture1)
true.order <- list() ## list specifying the edges in the components of mixture1
fit.order <- list() ## list specifying the edges in the components of mixture2
start.vals1 <- vector(mode = "numeric", length = (no.trees - 1)) ## the depth of the tree components in mixture1 (without the noise component)
start.vals2 <- vector(mode = "numeric", length = (no.trees - 1)) ## the depth of the tree components in mixture2 (without the noise component)
## Build the vectors that characterize the edges of the components of the mixture models.
## The noise components are ignored since they always have the same (star) topology.
for(k in 2:no.trees) {
true.vec <- character(0)
fit.vec <- character(0)
for(l in 2:no.events) {
ch <- names(which(sapply(edges(getTree(mixture1, k)), function(x, el) {el %in% x}, nodes(getTree(mixture1, k))[l])))
true.vec <- c(true.vec, ifelse((identical(ch, character(0)) || is.null(ch)), "e", ch))
ch <- names(which(sapply(edges(getTree(mixture2, k)), function(x, el) {el %in% x}, nodes(getTree(mixture2, k))[l])))
fit.vec <- c(fit.vec, ifelse((identical(ch, character(0)) || is.null(ch)), "e", ch))
}
names(true.vec) <- nodes(getTree(mixture1, k))[-1]
true.order <- c(true.order, list(true.vec))
names(fit.vec) <- nodes(getTree(mixture2, k))[-1]
fit.order <- c(fit.order, list(fit.vec))
## Assign the (maximal depths + 1) for each tree in the mixtures
start.vals1[k - 1] <- max(acc(getTree(mixture1, k), nodes(getTree(mixture1, k))[1])[[1]] + 1, na.remove = TRUE)
start.vals2[k - 1] <- max(acc(getTree(mixture2, k), nodes(getTree(mixture2, k))[1])[[1]] + 1, na.remove = TRUE)
}
## Build the comparisson matrix.
## The (i, j) element is the number of distinct edges of the
## (i + 1)-th and (j + 1)-th component of the models
## mixture1 and mixture2 respectively.
comp.mat <- matrix(nrow = no.trees - 1, ncol = no.trees - 1)
for(i in 1:(no.trees - 1))
for(j in 1:(no.trees - 1)) {
comp.mat[i, j] <- sum(true.order[[i]] != fit.order[[j]])
}
## Form the similarity pairs and calculate the similarity of the models
## as a sum of the number of different edges of the trees in the similarity pairs
## and their corresponding L1-distance of the levels of the events:
if(no.trees > 2) {
rc <- which(comp.mat == min(comp.mat), arr.ind = TRUE)
diff.sum <- min(comp.mat) +
L1.dist(get.tree.levels(mixture1, rc[1, 1] + 1, start.vals2[rc[1, 2]]),
get.tree.levels(mixture2, rc[1, 2] + 1, start.vals1[rc[1, 1]]))
row.index <- c(rc[1, 1]) ## get the row index of the minimum
col.index <- c(rc[1, 2]) ## get the column index of the minimum
for(m in 1:(nrow(comp.mat) - 1)) {
rc <- which(comp.mat == min(comp.mat[-(row.index), -(col.index)]), arr.ind = TRUE)
diff.sum <- diff.sum + min(comp.mat[-(row.index), -(col.index)]) +
L1.dist(get.tree.levels(mixture1, rc[1, 1] + 1, start.vals2[rc[1, 2]]),
get.tree.levels(mixture2, rc[1, 2] + 1, start.vals1[rc[1, 1]]))
row.index <- c(row.index, rc[1, 1])
col.index <- c(col.index, rc[1, 2])
}
} else {
if(no.trees == 2) {
diff.sum <- comp.mat[1, 1] +
L1.dist(get.tree.levels(mixture1, 2, start.vals2[1]),
get.tree.levels(mixture2, 2, start.vals1[1]))
} else {
stop("The specified mixture models must have at least two tree components, where the first one is the noise component!")
}
}
return(diff.sum)
}
######################################################################
## Function that implements a similarity measure for comparing the topologies
## of the nontrivial tree components of a specified mixture model.
## It returns a value that uses the number of different edges for caracterizing
## the difference of the trees in the model.
## For the comparisson it is necessary that the model has at least two
## nontrivial components.
comp.trees <- function(model) {
no.trees <- numTrees(model)
if(no.trees <= 2)
stop("The specified mixture model should have at least two nontrivial tree components.")
no.events <- eventsNum(model)
true.order <- list() ## list specifying the edges in the nontrivial components of the model
## Build the vectors that characterize the
## edges of the nontrivial components of the mixture model.
for(k in 2:no.trees) {
true.vec <- character(0)
for(l in 2:no.events) {
ch <- names(which(sapply(edges(getTree(model, k)), function(x, el) {el %in% x}, nodes(getTree(model, k))[l])))
true.vec <- c(true.vec, ifelse((identical(ch, character(0)) || is.null(ch)), "e", ch))
}
names(true.vec) <- nodes(getTree(model, k))[-1]
true.order <- c(true.order, list(true.vec))
}
## Calculate the similarity of the nontrivial components of the model
## as a sum of the number of different edges of all combinations of
## two different nontrivial trees:
diff.sum <- 0
for(k1 in 2:(no.trees - 1))
for(k2 in (k1 + 1):no.trees)
diff.sum <- diff.sum + sum(true.order[[k1 - 1]] != true.order[[k2 - 1]])
## Return a result between 0 and 1.
return(diff.sum / (choose((no.trees - 1), 2) * (no.events - 1)))
}
######################################################################
## Function that implements a similarity measure for comparing the topologies
## of the nontrivial tree components of a specified mixture model.
## It returns a value that uses the number of different edges and the depth at
## which the events are, for caracterizing the difference of the trees in the model.
## For the comparisson it is necessary that the model has at least two
## nontrivial components.
comp.trees.levels <- function(model) {
no.trees <- numTrees(model)
if(no.trees <= 2)
stop("The specified mixture model should have at least two nontrivial tree components.")
no.events <- eventsNum(model)
start.vals <- vector(mode = "numeric", length = (no.trees - 1)) ## the depth of the nontrivial tree components in the model
true.order <- list() ## list specifying the edges in the nontrivial components of the model
## Build the vectors that characterize the
## edges of the nontrivial components of the mixture model.
for(k in 2:no.trees) {
true.vec <- character(0)
for(l in 2:no.events) {
ch <- names(which(sapply(edges(getTree(model, k)), function(x, el) {el %in% x}, nodes(getTree(model, k))[l])))
true.vec <- c(true.vec, ifelse((identical(ch, character(0)) || is.null(ch)), "e", ch))
}
names(true.vec) <- nodes(getTree(model, k))[-1]
true.order <- c(true.order, list(true.vec))
## Assign the maximal depths + 1 for each tree in the mixtures
start.vals[k - 1] <- max(acc(getTree(model, k), nodes(getTree(model, k))[1])[[1]] + 1, na.remove = TRUE)
}
## Calculate the similarity of the models as a sum of the number of
## different edges of all possible combinations of two different
## nontrivial trees in the model and their corresponding L1-distance
## of the levels of the events:
diff.sum <- 0
for(k1 in 2:(no.trees - 1)) {
for(k2 in (k1 + 1):no.trees) {
diff.sum <- diff.sum + (sum(true.order[[k1 - 1]] != true.order[[k2 - 1]])
+ L1.dist(get.tree.levels(model, k1, start.vals[k2 - 1]),
get.tree.levels(model, k2, start.vals[k1 - 1])))
}
}
return(diff.sum)
}
######################################################################
## Stability analysis of the oncogenetic trees mixture model.
## This includes stability analysis on different levels of the mixture
## model: GPS values, encoded probability distribution, tree topologies.
## The function outputs a list of analysis for the mentioned features.
## Each analysis contains the values of different similarity measures
## with their corresponding p-values.
## The models should have at least two components (the first one is the noise).
stability.sim <- function(no.trees = 3, ## number of tree components
no.events = 9, ## number of genetic events
prob = c(0.2, 0.8), ## interval for the edgeweights of the random mixture models
no.draws = 300, ## number of samples drawn from the random models used for learning a mixture model
no.rands = 100, ## number of rands for calculatin the p-values
no.sim = 1 ## number of simulation iterations
) {
## Set the true types.
no.trees <- as.integer(no.trees)
no.events <- as.integer(no.events)
no.draws <- as.integer(no.draws)
no.rands <- as.integer(no.rands)
no.sim <- as.integer(no.sim)
## Check if the necessary parameters are provided and have the correct form.
if(no.trees < 2)
stop("The specified mixture model should have at least two
tree components (where the first one is the noise).")
if (no.events < 1)
stop("The number of events must be an integer greater than zero!")
if(no.draws < 1)
stop("The number of draws (number of samples) must be an integer greater than zero!")
if (no.rands < 1)
stop("The number of random models for the p-values must be an integer greater than zero!")
if (no.sim < 1)
stop("The number of iterations for the waiting time simulation
must be an integer greater than zero!")
if(!missing(prob)) {
if(!is.numeric(prob) || (length(prob) != 2))
stop("Specify the boundaries of the conditional probabilities as a numeric vector
of length two = c(min, max)!")
if(prob[2] < prob[1])
stop("In the probability vector the lower boundary must be smaller than the
upper boundary!")
}
simGPS <- numeric(0) ## results from the stability analysis of the GPS values
topo.dif <- numeric(0) ## results from the stability analysis of the topologies of the tree components of mixture models based on comp.models
topo.levels.dif <- numeric(0) ## results from the stability analysis of the topologies of the tree components of mixture models based on comp.models.levels
result.distr <- numeric(0) ## results from the stability analysis of the probability distribution encoded by the model
mat.true.gps <- numeric(0) ## matrix containing the true GPS values from each simulation iteration
mat.fit.gps <- numeric(0) ## matrix containing the corresponding fitted GPS values from each simulation iteration
list.true.models <- list() ## list containing the true models from each simulation iteration
list.fit.models <- list() ## list containing the corresponding fitted models from each simulation iteration
## Simulation iterations:
for(i in 1:as.integer(no.sim)) {
print(i) ## simulation iteration
## Pick a true model from the space of random models and draw data from it.
true.m <- generate(K = no.trees, no.events = no.events, prob = prob)
Weights(true.m) <- c(0.05, rep(0.95/(no.trees - 1), (no.trees - 1)))
sim.data <- sim(model = true.m, no.draws = no.draws)
list.true.models <- c(list.true.models, true.m)
## Calculate the GPS and the distribution encoded by true.m
true.gps <- GPS(gps(model = true.m, data = sim.data, no.sim = 10000))
mat.true.gps <- cbind(mat.true.gps, true.gps)
true.distr <- distribution(model = true.m)$probability
## Fit a mixture model to the data sim.data
fm <- fit(data = sim.data, K = no.trees, noise = TRUE, equal.edgeweights = TRUE, eps = 0.01)
list.fit.models <- c(list.fit.models, fm)
## Calculate the GPS and the distribution encoded by fm
fit.gps <- GPS(gps(model = fm, data = sim.data, no.sim = 10000))
mat.fit.gps <- cbind(mat.fit.gps, fit.gps)
fit.distr <- distribution(model = fm)$probability
## Compute different distances between the distributions induced
## by the true and fitted models:
true.cosin <- cosin.dist(true.distr, fit.distr)
true.l1 <- L1.dist(true.distr, fit.distr)
true.kull <- kullback.leibler(true.distr, fit.distr)
## Compute different distances between the GPS values
## for the true and fitted models:
true.euclGPS <- euclidian.dist(true.gps, fit.gps)
true.rank.dist <- rank.cor.dist(true.gps, fit.gps)
## Compute different similarity measures for comparing the
## topologies of the nontrivial components of the true and fitted models:
true.topo.dif <- comp.models(true.m, fm)
true.topo.levels.dif <- comp.models.levels(true.m, fm)
## Vectors for keeping the values for the different features
## of the random mixture models needed for the p-values calculation:
rand.euclGPS <- numeric(0)
rand.rankGPS <- numeric(0)
rand.cos.distr <- numeric(0)
rand.l1.distr <- numeric(0)
rand.kull.distr <- numeric(0)
rand.topo <- numeric(0)
rand.topo.levels <- numeric(0)
## Create the vectors with random values for the p-values calculation:
for(j in 1:(as.integer(no.rands) - 1)) {
## Generate random model and calculate the GPS and distribution:
model <- generate(K = no.trees, no.events = no.events, prob = prob)
Weights(model) <- c(0.05, rep(0.95/(no.trees - 1), (no.trees - 1)))
random.gps <- GPS(gps(model = model, data = sim.data, no.sim = 10000))
random.distr <- distribution(model = model)$probability
## GPS:
rand.euclGPS <- c(rand.euclGPS, euclidian.dist(true.gps, random.gps))
rand.rankGPS <- c(rand.rankGPS, rank.cor.dist(true.gps, random.gps))
## Distribution:
rand.cos.distr <- c(rand.cos.distr, cosin.dist(true.distr, random.distr))
rand.l1.distr <- c(rand.l1.distr, L1.dist(true.distr, random.distr))
rand.kull.distr <- c(rand.kull.distr, kullback.leibler(true.distr, random.distr))
## Tree topologies:
rand.topo <- c(rand.topo, comp.models(true.m, model))
rand.topo.levels <- c(rand.topo.levels, comp.models.levels(true.m, model))
}
## Stability analysis of GPS values
## (using Euclidian distance and rank correlation distance as similarity measures):
simGPS <- rbind(simGPS,
c(true.euclGPS, Pval.dist(true.euclGPS, rand.euclGPS),
true.rank.dist, Pval.dist(true.rank.dist, rand.rankGPS)))
## Stability analysis of induced distributions
## (using cosine distance, L1 distance, Kullback-Leibler divergence as similarity measures):
result.distr <- rbind(result.distr,
c(true.cosin, Pval.dist(true.cosin, rand.cos.distr),
true.l1, Pval.dist(true.l1, rand.l1.distr),
true.kull, Pval.dist(true.kull, rand.kull.distr)))
## Stability analysis of tree topologies
## (using comp.models and comp.models.levels as similarity measures):
topo.dif <- rbind(topo.dif,
c(true.topo.dif, Pval.dist(true.topo.dif, rand.topo)))
topo.levels.dif <- rbind(topo.levels.dif,
c(true.topo.levels.dif, Pval.dist(true.topo.levels.dif, rand.topo.levels)))
}
## Organize the output:
colnames(simGPS) <- c("eucl.dist", "p.val.eucl", "rank.cor.dist", "p.val.rank")
colnames(result.distr) <- c("cos", "p.val.cos", "L1", "p.val.L1", "KL", "p.val.KL")
colnames(topo.dif) <- c("topo.dif", "p.value")
colnames(topo.levels.dif) <- c("topo.levels.dif", "p.value")
colnames(mat.true.gps) <- c(1:no.sim)
colnames(mat.fit.gps) <- c(1:no.sim)
output <- list(simGPS, result.distr,
topo.dif, topo.levels.dif,
mat.true.gps, mat.fit.gps,
list.true.models, list.fit.models)
names(output) <- c("GPS", "Distribution", "Tree topologies (distinct edges)",
"Tree topologies (distinct edges + levelsL1dist)", "True GPS",
"Fitted GPS", "True models", "Fitted models")
return(output)
}
|
