Nothing
TNBC
In funkycells: Functional Data Analysis for Multiplexed Cell Images
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
# knitr::opts_chunk$set(eval = !is_check)
library(funkycells)
full_run <- FALSE # This runs only computations that are feasibly quick
paperRun <- FALSE * full_run # This creates and saves figures as in paper
save_path <- tempdir() # Location to save files
specify_decimal <- function(x, k) {
trimws(format(round(x, k), nsmall = k))
}
This vignette documents a sample analysis of Triple Negative Breast Cancer (TNBC) data, retrieved from \url{https://www.angelolab.com/mibi-data}. We consider the classified phenotypes rather than the direct protein data, pre-loaded in this package.
TNBC_pheno <- TNBC_pheno
TNBC_meta <- TNBC_meta
We begin by defining the phenotypes of interest as well as the related interactions.
phenos <- unique(TNBC_pheno$Phenotype)
pheno_interactions <- rbind(
data.frame(t(combn(phenos, 2))),
data.frame("X1" = phenos, "X2" = phenos)
)
phenos_subset <- c(
"Tumor", "CD3T", "CD4T", "CD8T", "B", "DC",
"DCMono", "Macrophage", "MonoNeu", "NK", "Treg"
)
pheno_interactions_subset <- data.frame(
Var1 = rep("Tumor", 11),
Var2 = c(
"Tumor", "CD3T", "CD4T", "CD8T", "B", "DC",
"DCMono", "Macrophage", "MonoNeu", "NK", "Treg"
)
)
We notice there are r length(phenos) total phenotypes with
r nrow(pheno_interactions) total possible interactions. Alternatively, we can
also focus on a subset of phenotypes (r length(phenos_subset)) and
interactions (r nrow(pheno_interactions_subset)) of particular interest. Another analysis on the subset interaction set has been considered but was suppressed for brevity.
In order to build confidence in our approach, we begin by simulating data with a similar structure to the full phenotype data and evaluating the effectiveness of the approach.
To start, let us determine the agent information in a typical image.
data_pheno_stat <- data.frame("pheno" = phenos)
for (person in unique(TNBC_pheno$Person)) {
tmp <- as.data.frame(table(TNBC_pheno[TNBC_pheno$Person == person, "Phenotype"]))
colnames(tmp) <- c("pheno", person)
data_pheno_stat <- merge(x = data_pheno_stat, y = tmp, all.x = TRUE)
}
agent_info_subset <- data.frame(data_pheno_stat[1],
"mean" = rowMeans(data_pheno_stat[-1], na.rm = TRUE),
"med" = apply(data_pheno_stat[-1], 1, median, na.rm = TRUE)
)
tnbc_fig <- plotPP(TNBC_pheno[TNBC_pheno$Person == 2, c("cellx", "celly", "Phenotype")],
ptSize = 3, colorGuide = "none", dropAxes = T
)
png(paste0(save_path,"/paper/tnbc_person2_phenotypes.png"),
width = 800, height = 800)
print(tnbc_fig)
dev.off()
kf <- getKFunction(
data = TNBC_pheno[TNBC_pheno$Person == 2, -1],
agents = c("Tumor", "MonoNeu"),
unit = "Person",
rCheckVals = seq(0, 50, 1)
)
k_fig <- ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = r, y = K1), data = kf, linewidth = 2) +
ggplot2::geom_line(ggplot2::aes(x = seq(0, 50, 1), y = pi * seq(0, 50, 1)^2),
linewidth = 2, col = "red", linetype = "dashed"
) +
ggplot2::theme_bw() +
ggplot2::theme(
axis.title = ggplot2::element_blank(),
axis.text = ggplot2::element_blank()
)
png(paste0(save_path,"/paper/tnbc_person2_tumor_mononeu_k.png"),
width = 800, height = 800)
print(k_fig)
dev.off()
Now let us generate two data sets: (1) where no informative interactions are
present and (2) where some informative interactions are present. We generate
$34$ patients, comparative to the r length(unique(TNBC_pheno$Person)) TNBC
patients. We generate $15$ for each class with possible interactions, if
present, and $2$ for each class that are indistinguishable to add some noise to
the data.
The non-interactive data is simulated as follows.
set.seed(123)
dat0 <- simulatePP(
agentVarData =
data.frame(
"outcome" = c(0, 1),
"c1" = c(0, 0),
"c2" = c(1 / 25, 1 / 25), "c3" = c(1 / 50, 1 / 50),
"c4" = c(0, 0),
"c5" = c(0, 0), "c6" = c(0, 0),
"c7" = c(0, 0), "c8" = c(1 / 100, 1 / 100),
"c9" = c(1 / 20, 1 / 20), "c10" = c(1 / 250, 1 / 250),
"c11" = c(1 / 100, 1 / 100), "c12" = c(1 / 80, 1 / 80),
"c13" = c(0, 0), "c14" = c(0, 0),
"c15" = c(0, 0), "c16" = c(1 / 10, 1 / 10)
),
agentKappaData = data.frame(
"agent" = paste0("c", 1:16),
"clusterAgent" = c(NA, "c1", "c1", rep(NA, 12), "c1"),
"kappa" = c(
rbinom(1, 100, 0.5),
rbinom(1, 100, 0.5),
rbinom(1, 30, 0.5),
rbinom(1, 80, 0.5),
rbinom(1, 350, 0.5),
rbinom(1, 100, 0.5),
rbinom(1, 120, 0.5),
rbinom(1, 150, 0.5),
rbinom(2, 250, 0.5),
rbinom(1, 600, 0.5),
rbinom(1, 60, 0.5),
rbinom(2, 140, 0.5),
rbinom(1, 20, 0.5),
rbinom(1, 5, 0.5)
)
),
unitsPerOutcome = 15,
replicatesPerUnit = 1,
silent = FALSE
)
dat1 <- simulatePP(
agentVarData =
data.frame(
"outcome" = c(0, 1),
"c1" = c(0, 0),
"c2" = c(1 / 25, 1 / 25), "c3" = c(1 / 50, 1 / 50),
"c4" = c(0, 0),
"c5" = c(0, 0), "c6" = c(0, 0),
"c7" = c(0, 0), "c8" = c(1 / 100, 1 / 100),
"c9" = c(1 / 20, 1 / 20), "c10" = c(1 / 250, 1 / 250),
"c11" = c(1 / 100, 1 / 100), "c12" = c(1 / 80, 1 / 80),
"c13" = c(0, 0), "c14" = c(0, 0),
"c15" = c(0, 0), "c16" = c(1 / 10, 1 / 10)
),
agentKappaData = data.frame(
"agent" = paste0("c", 1:16),
"clusterAgent" = c(NA, "c1", "c1", rep(NA, 12), "c1"),
"kappa" = c(
rbinom(1, 100, 0.5),
rbinom(1, 100, 0.5),
rbinom(1, 30, 0.5),
rbinom(1, 80, 0.5),
rbinom(1, 350, 0.5),
rbinom(1, 100, 0.5),
rbinom(1, 120, 0.5),
rbinom(1, 150, 0.5),
rbinom(2, 250, 0.5),
rbinom(1, 600, 0.5),
rbinom(1, 60, 0.5),
rbinom(2, 140, 0.5),
rbinom(1, 20, 0.5),
rbinom(1, 5, 0.5)
)
),
unitsPerOutcome = 2,
replicatesPerUnit = 1,
silent = F
)
dat1$unit <- ifelse(dat1$unit == "u1", "u31",
ifelse(dat1$unit == "u2", "u32",
ifelse(dat1$unit == "u3", "u33",
ifelse(dat1$unit == "u4", "u34", NA)
)
)
)
dat1$replicate <- ifelse(dat1$replicate == "1", "31",
ifelse(dat1$replicate == "2", "32",
ifelse(dat1$replicate == "3", "33",
ifelse(dat1$replicate == "4", "34", NA)
)
)
)
dat <- rbind(dat0, dat1)
Compare images of the data.
plotPP(
data = TNBC_pheno[TNBC_pheno$Person == 2, c("cellx", "celly", "Phenotype")],
ptSize = 1, dropAxes = T, colorGuide = "none"
)
plotPP(dat[dat$replicate == 18, c("x", "y", "type")],
ptSize = 1, dropAxes = T, colorGuide = "none"
)
tmp <- length(table(dat[dat$replicate == 1, c("type")]))
clrs <- scales::hue_pal()(tmp)
tnbc_image <- plotPP(
data = TNBC_pheno[TNBC_pheno$Person == 2, c("cellx", "celly", "Phenotype")],
ptSize = 3, dropAxes = T, colorGuide = "none",
colors = clrs[-length(clrs)]
)
png(paste0(save_path,"/paper/tnbc_image.png"),
width = 800, height = 800)
print(tnbc_image)
dev.off()
sim_image <- plotPP(
data = dat[dat$replicate == 18, c("x", "y", "type")],
ptSize = 3, dropAxes = T, colorGuide = "none",
colors = c(clrs[1:2], clrs[16], clrs[3:7], clrs[15], clrs[8:14])
)
png(paste0(save_path,"/paper/sim_image.png"),
width = 800, height = 800)
print(sim_image)
dev.off()
Now we compute PCA and attach an age meta-variable (with no effect).
cells <- paste0("c", 1:16)
cells_interactions <- rbind(
data.frame(t(combn(cells, 2))),
data.frame("X1" = cells, "X2" = cells)
)
set.seed(12345)
pcaData <- getKsPCAData(
data = dat, replicate = "replicate",
agents_df = cells_interactions,
xRange = c(0, 1), yRange = c(0, 1),
silent = F
)
pcaMeta <- simulateMeta(pcaData,
metaInfo = data.frame(
"var" = c("age"),
"rdist" = c("rnorm"),
"outcome_0" = c("25"),
"outcome_1" = c("25")
)
)
We analyze the data.
set.seed(123)
rfcv <- funkyModel(
data = pcaMeta,
outcome = "outcome",
unit = "unit",
metaNames = c("age")
)
Resulting in the variable importance figures. The code for these is given below, but may not be shown due to computational time in creating this vignette.
rfcv$viPlot
rfcv$subset_viPlot
# Get Vars
viData <- rfcv$VariableImportance
accData <- rfcv$AccuracyEstimate
NoiseCutoff <- rfcv$NoiseCutoff
InterpolationCutoff <- rfcv$InterpolationCutoff
subsetPlotSize <- rfcv$AdditionalParams$subsetPlotSize
# Plot Full
maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est)
plot_full <- ggplot2::ggplot(
data = viData,
mapping = ggplot2::aes(
x = factor(stats::reorder(var, est)),
y = ifelse(est / maxVal > 1, 1,
ifelse(est / maxVal < 0, 0,
est / maxVal
)
),
group = 1
)
) +
ggplot2::geom_errorbar(
ggplot2::aes(
ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal),
ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal)
),
color = "black", width = 0.2
) +
ggplot2::geom_point(color = "black", size = 3) +
ggplot2::geom_hline(
ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))),
color = "red", linetype = "dotted", linewidth = 2
) +
ggplot2::coord_flip(ylim = c(0, 1)) +
ggplot2::xlab(NULL) +
ggplot2::ylim(c(0, 1)) +
ggplot2::ylab(NULL) +
ggplot2::theme_bw() +
ggplot2::theme(
axis.text = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_blank(),
panel.grid.major.y = ggplot2::element_blank(),
axis.title = ggplot2::element_text(size = 20)
) +
ggplot2::geom_line(ggplot2::aes(
x = ordered(viData[order(-est), "var"]),
y = InterpolationCutoff / maxVal
), color = "orange", linetype = "dashed", linewidth = 2) +
ggplot2::ylab(paste0(
"OOB (",
specify_decimal(accData$OOB, 2),
"), Guess (",
specify_decimal(accData$guess, 2),
"), Bias (",
specify_decimal(accData$bias, 2), ")"
))
png(paste0(save_path,"/paper/noeffect_variableimportance.png"))
print(plot_full)
dev.off()
# Plot Subset
InterpolationCutoff <- InterpolationCutoff[1:subsetPlotSize]
viData <- viData[order(-viData$est), ]
viData <- viData[1:subsetPlotSize, ]
maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est)
plot_top25 <- ggplot2::ggplot(
data = viData,
mapping = ggplot2::aes(
x = factor(stats::reorder(var, est)),
y = ifelse(est / maxVal > 1, 1,
ifelse(est / maxVal < 0, 0,
est / maxVal
)
),
group = 1
)
) +
ggplot2::geom_errorbar(
ggplot2::aes(
ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal),
ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal)
),
color = "black", width = 0.2
) +
ggplot2::geom_point(color = "black", size = 3) +
ggplot2::geom_hline(
ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))),
color = "red", linetype = "dotted", linewidth = 2
) +
ggplot2::coord_flip(ylim = c(0, 1)) +
ggplot2::xlab(NULL) +
ggplot2::ylim(c(0, 1)) +
ggplot2::ylab(NULL) +
ggplot2::theme_bw() +
ggplot2::theme(
axis.text = ggplot2::element_text(size = 18),
axis.text.x = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_blank(),
panel.grid.major.y = ggplot2::element_blank(),
axis.title = ggplot2::element_text(size = 20)
) +
ggplot2::geom_line(ggplot2::aes(
x = ordered(viData[order(-est), "var"]),
y = InterpolationCutoff / maxVal
), color = "orange", linetype = "dashed", linewidth = 2) +
ggplot2::ylab(paste0(
"OOB (",
specify_decimal(accData$OOB, 2),
"), Guess (",
specify_decimal(accData$guess, 2),
"), Bias (",
specify_decimal(accData$bias, 2), ")"
))
png(paste0(save_path,"/paper/noeffect_variableimportance_top25.png"))
print(plot_top25)
dev.off()
Now we consider the simulated model with effects.
set.seed(123456)
dat0 <- simulatePP(
agentVarData =
data.frame(
"outcome" = c(0, 1),
"c1" = c(0, 0),
"c2" = c(1 / 25, 1 / 60), "c3" = c(1 / 50, 1 / 10),
"c4" = c(0, 0),
"c5" = c(0, 0), "c6" = c(0, 0),
"c7" = c(0, 0), "c8" = c(1 / 100, 1 / 100),
"c9" = c(1 / 20, 1 / 20), "c10" = c(1 / 250, 1 / 250),
"c11" = c(1 / 100, 1 / 100), "c12" = c(1 / 80, 1 / 80),
"c13" = c(0, 0), "c14" = c(0, 0),
"c15" = c(0, 0), "c16" = c(1 / 10, 1 / 10)
),
agentKappaData = data.frame(
"agent" = paste0("c", 1:16),
"clusterAgent" = c(NA, "c1", "c1", rep(NA, 12), "c1"),
"kappa" = c(
rbinom(1, 100, 0.5),
rbinom(1, 100, 0.5),
rbinom(1, 30, 0.5),
rbinom(1, 80, 0.5),
rbinom(1, 350, 0.5),
rbinom(1, 100, 0.5),
rbinom(1, 120, 0.5),
rbinom(1, 150, 0.5),
rbinom(2, 250, 0.5),
rbinom(1, 600, 0.5),
rbinom(1, 60, 0.5),
rbinom(2, 140, 0.5),
rbinom(1, 20, 0.5),
rbinom(1, 5, 0.5)
)
),
unitsPerOutcome = 15,
replicatesPerUnit = 1,
silent = F
)
dat1 <- simulatePP(
agentVarData =
data.frame(
"outcome" = c(0, 1),
"c1" = c(0, 0),
"c2" = c(1 / 25, 1 / 25), "c3" = c(1 / 50, 1 / 50),
"c4" = c(0, 0),
"c5" = c(0, 0), "c6" = c(0, 0),
"c7" = c(0, 0), "c8" = c(1 / 100, 1 / 100),
"c9" = c(1 / 20, 1 / 20), "c10" = c(1 / 250, 1 / 250),
"c11" = c(1 / 100, 1 / 100), "c12" = c(1 / 80, 1 / 80),
"c13" = c(0, 0), "c14" = c(0, 0),
"c15" = c(0, 0), "c16" = c(1 / 10, 1 / 10)
),
agentKappaData = data.frame(
"agent" = paste0("c", 1:16),
"clusterAgent" = c(NA, "c1", "c1", rep(NA, 12), "c1"),
"kappa" = c(
rbinom(1, 100, 0.5),
rbinom(1, 100, 0.5),
rbinom(1, 30, 0.5),
rbinom(1, 80, 0.5),
rbinom(1, 350, 0.5),
rbinom(1, 100, 0.5),
rbinom(1, 120, 0.5),
rbinom(1, 150, 0.5),
rbinom(2, 250, 0.5),
rbinom(1, 600, 0.5),
rbinom(1, 60, 0.5),
rbinom(2, 140, 0.5),
rbinom(1, 20, 0.5),
rbinom(1, 5, 0.5)
)
),
unitsPerOutcome = 2,
replicatesPerUnit = 1,
silent = F
)
dat1$unit <- ifelse(dat1$unit == "u1", "u31",
ifelse(dat1$unit == "u2", "u32",
ifelse(dat1$unit == "u3", "u33",
ifelse(dat1$unit == "u4", "u34", NA)
)
)
)
dat1$replicate <- ifelse(dat1$replicate == "1", "31",
ifelse(dat1$replicate == "2", "32",
ifelse(dat1$replicate == "3", "33",
ifelse(dat1$replicate == "4", "34", NA)
)
)
)
dat <- rbind(dat0, dat1)
cells <- paste0("c", 1:16)
cells_interactions <- rbind(
data.frame(t(combn(cells, 2))),
data.frame("X1" = cells, "X2" = cells)
)
pcaData <- getKsPCAData(
data = dat, replicate = "replicate",
agents_df = cells_interactions,
xRange = c(0, 1), yRange = c(0, 1),
silent = F
)
pcaMeta <- simulateMeta(pcaData,
metaInfo = data.frame(
"var" = c("age"),
"rdist" = c("rnorm"),
"outcome_0" = c("25"),
"outcome_1" = c("27")
)
)
rfcv <- funkyModel(
data = pcaMeta,
outcome = "outcome",
unit = "unit",
metaNames = c("age")
)
Creating the variable importance plots as before (code below, perhaps figure not given for computational reasons)
rfcv$viPlot
rfcv$subset_viPlot
# Get Vars
viData <- rfcv$VariableImportance
accData <- rfcv$AccuracyEstimate
NoiseCutoff <- rfcv$NoiseCutoff
InterpolationCutoff <- rfcv$InterpolationCutoff
subsetPlotSize <- rfcv$AdditionalParams$subsetPlotSize
# Plot Full
maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est)
plot_full <- ggplot2::ggplot(
data = viData,
mapping = ggplot2::aes(
x = factor(stats::reorder(var, est)),
y = ifelse(est / maxVal > 1, 1,
ifelse(est / maxVal < 0, 0,
est / maxVal
)
),
group = 1
)
) +
ggplot2::geom_errorbar(
ggplot2::aes(
ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal),
ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal)
),
color = "black", width = 0.2
) +
ggplot2::geom_point(color = "black", size = 3) +
ggplot2::geom_hline(
ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))),
color = "red", linetype = "dotted", linewidth = 2
) +
ggplot2::coord_flip(ylim = c(0, 1)) +
ggplot2::xlab(NULL) +
ggplot2::ylim(c(0, 1)) +
ggplot2::ylab(NULL) +
ggplot2::theme_bw() +
ggplot2::theme(
axis.text = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_blank(),
panel.grid.major.y = ggplot2::element_blank(),
axis.title = ggplot2::element_text(size = 20)
) +
ggplot2::geom_line(ggplot2::aes(
x = ordered(viData[order(-est), "var"]),
y = InterpolationCutoff / maxVal
), color = "orange", linetype = "dashed", linewidth = 2) +
ggplot2::ylab(paste0(
"OOB (",
specify_decimal(accData$OOB, 2),
"), Guess (",
specify_decimal(accData$guess, 2),
"), Bias (",
specify_decimal(accData$bias, 2), ")"
))
png(paste0(save_path,"/paper/effect_variableimportance.png"))
print(plot_full)
dev.off()
# Plot Subset
InterpolationCutoff <- InterpolationCutoff[1:subsetPlotSize]
viData <- viData[order(-viData$est), ]
viData <- viData[1:subsetPlotSize, ]
maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est)
plot_top25 <- ggplot2::ggplot(
data = viData,
mapping = ggplot2::aes(
x = factor(stats::reorder(var, est)),
y = ifelse(est / maxVal > 1, 1,
ifelse(est / maxVal < 0, 0,
est / maxVal
)
),
group = 1
)
) +
ggplot2::geom_errorbar(
ggplot2::aes(
ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal),
ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal)
),
color = "black", width = 0.2
) +
ggplot2::geom_point(color = "black", size = 3) +
ggplot2::geom_hline(
ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))),
color = "red", linetype = "dotted", linewidth = 2
) +
ggplot2::coord_flip(ylim = c(0, 1)) +
ggplot2::xlab(NULL) +
ggplot2::ylim(c(0, 1)) +
ggplot2::ylab(NULL) +
ggplot2::theme_bw() +
ggplot2::theme(
axis.text = ggplot2::element_text(size = 18),
axis.text.x = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_blank(),
panel.grid.major.y = ggplot2::element_blank(),
axis.title = ggplot2::element_text(size = 20)
) +
ggplot2::geom_line(ggplot2::aes(
x = ordered(viData[order(-est), "var"]),
y = InterpolationCutoff / maxVal
), color = "orange", linetype = "dashed", linewidth = 2) +
ggplot2::ylab(paste0(
"OOB (",
specify_decimal(accData$OOB, 2),
"), Guess (",
specify_decimal(accData$guess, 2),
"), Bias (",
specify_decimal(accData$bias, 2), ")"
))
png(paste0(save_path,"/paper/effect_variableimportance_top25.png"))
print(plot_top25)
dev.off()
With this confidence, we now consider the TNBC phenotype data.
dataPCA_pheno <- getKsPCAData(
data = TNBC_pheno, unit = "Person",
agents_df = pheno_interactions,
rCheckVals = seq(0, 50, 1)
)
dataPCAAge_pheno <- merge(dataPCA_pheno, TNBC_meta)
set.seed(123456)
rfcv <- funkyModel(
data = dataPCAAge_pheno, K = 10,
outcome = "Class", unit = "Person",
metaNames = c("Age"), synthetics = 100,
alpha = 0.05, silent = FALSE,
subsetPlotSize = 25
)
And the related variable importance plots (may not be shown for computational reasons).
rfcv$viPlot
rfcv$subset_viPlot
# Get Vars
viData <- rfcv$VariableImportance
accData <- rfcv$AccuracyEstimate
NoiseCutoff <- rfcv$NoiseCutoff
InterpolationCutoff <- rfcv$InterpolationCutoff
subsetPlotSize <- rfcv$AdditionalParams$subsetPlotSize
# Plot Full
maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est)
plot_full <- ggplot2::ggplot(
data = viData,
mapping = ggplot2::aes(
x = factor(stats::reorder(var, est)),
y = ifelse(est / maxVal > 1, 1,
ifelse(est / maxVal < 0, 0,
est / maxVal
)
),
group = 1
)
) +
ggplot2::geom_errorbar(
ggplot2::aes(
ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal),
ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal)
),
color = "black", width = 0.2
) +
ggplot2::geom_point(color = "black", size = 3) +
ggplot2::geom_hline(
ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))),
color = "red", linetype = "dotted", linewidth = 2
) +
ggplot2::coord_flip(ylim = c(0, 1)) +
ggplot2::xlab(NULL) +
ggplot2::ylim(c(0, 1)) +
ggplot2::ylab(NULL) +
ggplot2::theme_bw() +
ggplot2::theme(
# axis.text = ggplot2::element_text(size = 18),
axis.text = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_blank(),
panel.grid.major.y = ggplot2::element_blank(),
axis.title = ggplot2::element_text(size = 20)
) +
ggplot2::geom_line(ggplot2::aes(
x = ordered(viData[order(-est), "var"]),
y = InterpolationCutoff / maxVal
), color = "orange", linetype = "dashed", linewidth = 2) +
ggplot2::ylab(paste0(
"OOB (",
specify_decimal(accData$OOB, 2),
"), Guess (",
specify_decimal(accData$guess, 2),
"), Bias (",
specify_decimal(accData$bias, 2), ")"
))
png(paste0(save_path,"/paper/tnbc_variableimportance.png"))
print(plot_full)
dev.off()
# Plot Subset
InterpolationCutoff <- InterpolationCutoff[1:subsetPlotSize]
viData <- viData[order(-viData$est), ]
viData <- viData[1:subsetPlotSize, ]
maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est)
plot_top25 <- ggplot2::ggplot(
data = viData,
mapping = ggplot2::aes(
x = factor(stats::reorder(var, est)),
y = ifelse(est / maxVal > 1, 1,
ifelse(est / maxVal < 0, 0,
est / maxVal
)
),
group = 1
)
) +
ggplot2::geom_errorbar(
ggplot2::aes(
ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal),
ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal)
),
color = "black", width = 0.2
) +
ggplot2::geom_point(color = "black", size = 3) +
ggplot2::geom_hline(
ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))),
color = "red", linetype = "dotted", linewidth = 2
) +
ggplot2::coord_flip(ylim = c(0, 1)) +
ggplot2::xlab(NULL) +
ggplot2::ylim(c(0, 1)) +
ggplot2::ylab(NULL) +
ggplot2::theme_bw() +
ggplot2::theme(
axis.text.x = ggplot2::element_blank(),
axis.text.y = ggplot2::element_text(size = 14),
axis.ticks.y = ggplot2::element_blank(),
panel.grid.major.y = ggplot2::element_blank(),
axis.title = ggplot2::element_text(size = 17)
) +
ggplot2::geom_line(ggplot2::aes(
x = ordered(viData[order(-est), "var"]),
y = InterpolationCutoff / maxVal
), color = "orange", linetype = "dashed", linewidth = 2) +
ggplot2::ylab(paste0(
"OOB (",
specify_decimal(accData$OOB, 2),
"), Guess (",
specify_decimal(accData$guess, 2),
"), Bias (",
specify_decimal(accData$bias, 2), ")"
))
png(paste0(save_path,"/paper/tnbc_variableimportance_top25.png"))
print(plot_top25)
dev.off()
Also consider $K$ functions from significant and insignificant interactions.
tmp <- getKFunction(TNBC_pheno[TNBC_pheno$Class == 0, -1],
c("Tumor", "Tumor"),
unit = "Person",
rCheckVals = seq(0, 50, 1)
)
tmp1 <- getKFunction(TNBC_pheno[TNBC_pheno$Class == 1, -1],
c("Tumor", "Tumor"),
unit = "Person",
rCheckVals = seq(0, 50, 1)
)
tmp_1 <- tidyr::pivot_longer(data = tmp, cols = -r)
tmp1_1 <- tidyr::pivot_longer(data = tmp1, cols = -r)
data_plot <- rbind(
data.frame(
"r" = tmp_1$r,
"K" = tmp_1$value,
"unit" = tmp_1$name,
"outcome" = "0"
),
data.frame(
"r" = tmp1_1$r,
"K" = tmp1_1$value,
"unit" = paste0(tmp1_1$name, "_1"),
"outcome" = "1"
)
)
Creating the following figure.
plot_K_functions(data_plot)
data <- data_plot
# Prep Legend
info <- unique(data[, c("unit", "outcome")])
info$Missing <- NA
for (i in 1:nrow(info)) {
info[i, 3] <- nrow(data[data$unit == info$unit[i] &
!stats::complete.cases(data), ]) > 0
}
info_labels <- data.frame(
"outcome" = unique(data$outcome),
"units" = NA,
"Missing" = NA,
"Label" = NA
)
for (i in 1:nrow(info_labels)) {
info_labels[i, 2] <- nrow(info[info$outcome == info_labels[i, 1], ])
info_labels[i, 3] <- nrow(info[info$outcome == info_labels[i, 1] &
info$Missing, ])
# Only label outcomes with numbers if at least one is missing in whole set
if (sum(info$Missing) != 0) {
info_labels[i, 4] <- paste0(
info_labels[i, 1], " (",
info_labels[i, 2] - info_labels[i, 3], "/",
info_labels[i, 2], ")"
)
} else {
info_labels[i, 4] <- info_labels[i, 1]
}
}
# Build Averages
data_wide <- tidyr::pivot_wider(data,
names_from = "unit",
values_from = "K"
)
data_avg <- data.frame("r" = unique(data_wide$r))
for (i in 1:nrow(info_labels)) {
data_avg[[info_labels[i, "outcome"]]] <-
rowMeans(data_wide[data_wide$outcome == info_labels[i, "outcome"], -c(1:2)],
na.rm = TRUE
)
}
data_avg <- tidyr::pivot_longer(data_avg, cols = -r)
colnames(data_avg) <- c("r", "outcome", "Value")
# Plot
return_plot <-
ggplot2::ggplot(
data = stats::na.omit(data),
ggplot2::aes(
x = r, y = K,
group = unit, color = outcome
)
) +
ggplot2::geom_line(alpha = 0.2, linewidth = 1.25) +
ggplot2::geom_line(
ggplot2::aes(
x = r, y = Value,
group = outcome, color = outcome
),
data = data_avg, linewidth = 2
) +
ggplot2::theme_bw() +
ggplot2::theme(
axis.text = ggplot2::element_text(size = 18),
axis.text.y = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_blank(),
panel.grid.major.y = ggplot2::element_blank(),
# axis.title = ggplot2::element_text(size = 22),
axis.title = ggplot2::element_blank(),
legend.position = "none"
) # +
# ggplot2::geom_line(ggplot2::aes(x = r, y = pi * r^2,
# group = outcome),
# data = data_avg, linewidth = 1.25,
# color = "black", linetype = "dashed"
# )
png(paste0(save_path,"/paper/TNBC_sig_KFunctions.png"))
print(return_plot)
dev.off()
And for a insignificant interaction.
tmp <- getKFunction(TNBC_pheno[TNBC_pheno$Class == 0, -1],
c("CD4T", "Endothelial"),
unit = "Person",
rCheckVals = seq(0, 50, 1)
)
tmp1 <- getKFunction(TNBC_pheno[TNBC_pheno$Class == 1, -1],
c("CD4T", "Endothelial"),
unit = "Person",
rCheckVals = seq(0, 50, 1)
)
tmp_1 <- tidyr::pivot_longer(data = tmp, cols = -r)
tmp1_1 <- tidyr::pivot_longer(data = tmp1, cols = -r)
data_plot <- rbind(
data.frame(
"r" = tmp_1$r,
"K" = tmp_1$value,
"unit" = tmp_1$name,
"outcome" = "0"
),
data.frame(
"r" = tmp1_1$r,
"K" = tmp1_1$value,
"unit" = paste0(tmp1_1$name, "_1"),
"outcome" = "1"
)
)
Creating the following figure.
plot_K_functions(data_plot)
data <- data_plot
# Prep Legend
info <- unique(data[, c("unit", "outcome")])
info$Missing <- NA
for (i in 1:nrow(info)) {
info[i, 3] <- nrow(data[data$unit == info$unit[i] &
!stats::complete.cases(data), ]) > 0
}
info_labels <- data.frame(
"outcome" = unique(data$outcome),
"units" = NA,
"Missing" = NA,
"Label" = NA
)
for (i in 1:nrow(info_labels)) {
info_labels[i, 2] <- nrow(info[info$outcome == info_labels[i, 1], ])
info_labels[i, 3] <- nrow(info[info$outcome == info_labels[i, 1] &
info$Missing, ])
# Only label outcomes with numbers if at least one is missing in whole set
if (sum(info$Missing) != 0) {
info_labels[i, 4] <- paste0(
info_labels[i, 1], " (",
info_labels[i, 2] - info_labels[i, 3], "/",
info_labels[i, 2], ")"
)
} else {
info_labels[i, 4] <- info_labels[i, 1]
}
}
# Build Averages
data_wide <- tidyr::pivot_wider(data,
names_from = "unit",
values_from = "K"
)
data_avg <- data.frame("r" = unique(data_wide$r))
for (i in 1:nrow(info_labels)) {
data_avg[[info_labels[i, "outcome"]]] <-
rowMeans(data_wide[data_wide$outcome == info_labels[i, "outcome"], -c(1:2)],
na.rm = TRUE
)
}
data_avg <- tidyr::pivot_longer(data_avg, cols = -r)
colnames(data_avg) <- c("r", "outcome", "Value")
# Plot
return_plot <-
ggplot2::ggplot(
data = stats::na.omit(data),
ggplot2::aes(
x = r, y = K,
group = unit, color = outcome
)
) +
ggplot2::geom_line(alpha = 0.2, linewidth = 1.25) +
ggplot2::geom_line(
ggplot2::aes(
x = r, y = Value,
group = outcome, color = outcome
),
data = data_avg, linewidth = 2
) +
ggplot2::theme_bw() +
ggplot2::theme(
axis.text = ggplot2::element_text(size = 18),
axis.text.y = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_blank(),
panel.grid.major.y = ggplot2::element_blank(),
# axis.title = ggplot2::element_text(size = 22),
axis.title = ggplot2::element_blank(),
legend.position = "none"
) # +
# ggplot2::geom_line(ggplot2::aes(x = r, y = pi * r^2,
# group = outcome),
# data = data_avg, linewidth = 1.25,
# color = "black", linetype = "dashed"
# )
png(paste0(save_path,"/paper/TNBC_insig_KFunctions.png"))
print(return_plot)
dev.off()
Further analysis of this data can be found in our paper.
Try the funkycells package in your browser
Any scripts or data that you put into this service are public.
funkycells documentation built on Aug. 9, 2023, 5:10 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.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) # knitr::opts_chunk$set(eval = !is_check)
library(funkycells)
full_run <- FALSE # This runs only computations that are feasibly quick paperRun <- FALSE * full_run # This creates and saves figures as in paper save_path <- tempdir() # Location to save files specify_decimal <- function(x, k) { trimws(format(round(x, k), nsmall = k)) }
This vignette documents a sample analysis of Triple Negative Breast Cancer (TNBC) data, retrieved from \url{https://www.angelolab.com/mibi-data}. We consider the classified phenotypes rather than the direct protein data, pre-loaded in this package.
TNBC_pheno <- TNBC_pheno TNBC_meta <- TNBC_meta
We begin by defining the phenotypes of interest as well as the related interactions.
phenos <- unique(TNBC_pheno$Phenotype) pheno_interactions <- rbind( data.frame(t(combn(phenos, 2))), data.frame("X1" = phenos, "X2" = phenos) ) phenos_subset <- c( "Tumor", "CD3T", "CD4T", "CD8T", "B", "DC", "DCMono", "Macrophage", "MonoNeu", "NK", "Treg" ) pheno_interactions_subset <- data.frame( Var1 = rep("Tumor", 11), Var2 = c( "Tumor", "CD3T", "CD4T", "CD8T", "B", "DC", "DCMono", "Macrophage", "MonoNeu", "NK", "Treg" ) )
We notice there are r length(phenos) total phenotypes with
r nrow(pheno_interactions) total possible interactions. Alternatively, we can
also focus on a subset of phenotypes (r length(phenos_subset)) and
interactions (r nrow(pheno_interactions_subset)) of particular interest. Another analysis on the subset interaction set has been considered but was suppressed for brevity.
In order to build confidence in our approach, we begin by simulating data with a similar structure to the full phenotype data and evaluating the effectiveness of the approach.
To start, let us determine the agent information in a typical image.
data_pheno_stat <- data.frame("pheno" = phenos) for (person in unique(TNBC_pheno$Person)) { tmp <- as.data.frame(table(TNBC_pheno[TNBC_pheno$Person == person, "Phenotype"])) colnames(tmp) <- c("pheno", person) data_pheno_stat <- merge(x = data_pheno_stat, y = tmp, all.x = TRUE) } agent_info_subset <- data.frame(data_pheno_stat[1], "mean" = rowMeans(data_pheno_stat[-1], na.rm = TRUE), "med" = apply(data_pheno_stat[-1], 1, median, na.rm = TRUE) )
tnbc_fig <- plotPP(TNBC_pheno[TNBC_pheno$Person == 2, c("cellx", "celly", "Phenotype")], ptSize = 3, colorGuide = "none", dropAxes = T ) png(paste0(save_path,"/paper/tnbc_person2_phenotypes.png"), width = 800, height = 800) print(tnbc_fig) dev.off() kf <- getKFunction( data = TNBC_pheno[TNBC_pheno$Person == 2, -1], agents = c("Tumor", "MonoNeu"), unit = "Person", rCheckVals = seq(0, 50, 1) ) k_fig <- ggplot2::ggplot() + ggplot2::geom_line(ggplot2::aes(x = r, y = K1), data = kf, linewidth = 2) + ggplot2::geom_line(ggplot2::aes(x = seq(0, 50, 1), y = pi * seq(0, 50, 1)^2), linewidth = 2, col = "red", linetype = "dashed" ) + ggplot2::theme_bw() + ggplot2::theme( axis.title = ggplot2::element_blank(), axis.text = ggplot2::element_blank() ) png(paste0(save_path,"/paper/tnbc_person2_tumor_mononeu_k.png"), width = 800, height = 800) print(k_fig) dev.off()
Now let us generate two data sets: (1) where no informative interactions are
present and (2) where some informative interactions are present. We generate
$34$ patients, comparative to the r length(unique(TNBC_pheno$Person)) TNBC
patients. We generate $15$ for each class with possible interactions, if
present, and $2$ for each class that are indistinguishable to add some noise to
the data.
The non-interactive data is simulated as follows.
set.seed(123) dat0 <- simulatePP( agentVarData = data.frame( "outcome" = c(0, 1), "c1" = c(0, 0), "c2" = c(1 / 25, 1 / 25), "c3" = c(1 / 50, 1 / 50), "c4" = c(0, 0), "c5" = c(0, 0), "c6" = c(0, 0), "c7" = c(0, 0), "c8" = c(1 / 100, 1 / 100), "c9" = c(1 / 20, 1 / 20), "c10" = c(1 / 250, 1 / 250), "c11" = c(1 / 100, 1 / 100), "c12" = c(1 / 80, 1 / 80), "c13" = c(0, 0), "c14" = c(0, 0), "c15" = c(0, 0), "c16" = c(1 / 10, 1 / 10) ), agentKappaData = data.frame( "agent" = paste0("c", 1:16), "clusterAgent" = c(NA, "c1", "c1", rep(NA, 12), "c1"), "kappa" = c( rbinom(1, 100, 0.5), rbinom(1, 100, 0.5), rbinom(1, 30, 0.5), rbinom(1, 80, 0.5), rbinom(1, 350, 0.5), rbinom(1, 100, 0.5), rbinom(1, 120, 0.5), rbinom(1, 150, 0.5), rbinom(2, 250, 0.5), rbinom(1, 600, 0.5), rbinom(1, 60, 0.5), rbinom(2, 140, 0.5), rbinom(1, 20, 0.5), rbinom(1, 5, 0.5) ) ), unitsPerOutcome = 15, replicatesPerUnit = 1, silent = FALSE ) dat1 <- simulatePP( agentVarData = data.frame( "outcome" = c(0, 1), "c1" = c(0, 0), "c2" = c(1 / 25, 1 / 25), "c3" = c(1 / 50, 1 / 50), "c4" = c(0, 0), "c5" = c(0, 0), "c6" = c(0, 0), "c7" = c(0, 0), "c8" = c(1 / 100, 1 / 100), "c9" = c(1 / 20, 1 / 20), "c10" = c(1 / 250, 1 / 250), "c11" = c(1 / 100, 1 / 100), "c12" = c(1 / 80, 1 / 80), "c13" = c(0, 0), "c14" = c(0, 0), "c15" = c(0, 0), "c16" = c(1 / 10, 1 / 10) ), agentKappaData = data.frame( "agent" = paste0("c", 1:16), "clusterAgent" = c(NA, "c1", "c1", rep(NA, 12), "c1"), "kappa" = c( rbinom(1, 100, 0.5), rbinom(1, 100, 0.5), rbinom(1, 30, 0.5), rbinom(1, 80, 0.5), rbinom(1, 350, 0.5), rbinom(1, 100, 0.5), rbinom(1, 120, 0.5), rbinom(1, 150, 0.5), rbinom(2, 250, 0.5), rbinom(1, 600, 0.5), rbinom(1, 60, 0.5), rbinom(2, 140, 0.5), rbinom(1, 20, 0.5), rbinom(1, 5, 0.5) ) ), unitsPerOutcome = 2, replicatesPerUnit = 1, silent = F ) dat1$unit <- ifelse(dat1$unit == "u1", "u31", ifelse(dat1$unit == "u2", "u32", ifelse(dat1$unit == "u3", "u33", ifelse(dat1$unit == "u4", "u34", NA) ) ) ) dat1$replicate <- ifelse(dat1$replicate == "1", "31", ifelse(dat1$replicate == "2", "32", ifelse(dat1$replicate == "3", "33", ifelse(dat1$replicate == "4", "34", NA) ) ) ) dat <- rbind(dat0, dat1)
Compare images of the data.
plotPP( data = TNBC_pheno[TNBC_pheno$Person == 2, c("cellx", "celly", "Phenotype")], ptSize = 1, dropAxes = T, colorGuide = "none" )
plotPP(dat[dat$replicate == 18, c("x", "y", "type")], ptSize = 1, dropAxes = T, colorGuide = "none" )
tmp <- length(table(dat[dat$replicate == 1, c("type")])) clrs <- scales::hue_pal()(tmp) tnbc_image <- plotPP( data = TNBC_pheno[TNBC_pheno$Person == 2, c("cellx", "celly", "Phenotype")], ptSize = 3, dropAxes = T, colorGuide = "none", colors = clrs[-length(clrs)] ) png(paste0(save_path,"/paper/tnbc_image.png"), width = 800, height = 800) print(tnbc_image) dev.off() sim_image <- plotPP( data = dat[dat$replicate == 18, c("x", "y", "type")], ptSize = 3, dropAxes = T, colorGuide = "none", colors = c(clrs[1:2], clrs[16], clrs[3:7], clrs[15], clrs[8:14]) ) png(paste0(save_path,"/paper/sim_image.png"), width = 800, height = 800) print(sim_image) dev.off()
Now we compute PCA and attach an age meta-variable (with no effect).
cells <- paste0("c", 1:16) cells_interactions <- rbind( data.frame(t(combn(cells, 2))), data.frame("X1" = cells, "X2" = cells) ) set.seed(12345) pcaData <- getKsPCAData( data = dat, replicate = "replicate", agents_df = cells_interactions, xRange = c(0, 1), yRange = c(0, 1), silent = F ) pcaMeta <- simulateMeta(pcaData, metaInfo = data.frame( "var" = c("age"), "rdist" = c("rnorm"), "outcome_0" = c("25"), "outcome_1" = c("25") ) )
We analyze the data.
set.seed(123) rfcv <- funkyModel( data = pcaMeta, outcome = "outcome", unit = "unit", metaNames = c("age") )
Resulting in the variable importance figures. The code for these is given below, but may not be shown due to computational time in creating this vignette.
rfcv$viPlot
rfcv$subset_viPlot
# Get Vars viData <- rfcv$VariableImportance accData <- rfcv$AccuracyEstimate NoiseCutoff <- rfcv$NoiseCutoff InterpolationCutoff <- rfcv$InterpolationCutoff subsetPlotSize <- rfcv$AdditionalParams$subsetPlotSize # Plot Full maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est) plot_full <- ggplot2::ggplot( data = viData, mapping = ggplot2::aes( x = factor(stats::reorder(var, est)), y = ifelse(est / maxVal > 1, 1, ifelse(est / maxVal < 0, 0, est / maxVal ) ), group = 1 ) ) + ggplot2::geom_errorbar( ggplot2::aes( ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal), ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal) ), color = "black", width = 0.2 ) + ggplot2::geom_point(color = "black", size = 3) + ggplot2::geom_hline( ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))), color = "red", linetype = "dotted", linewidth = 2 ) + ggplot2::coord_flip(ylim = c(0, 1)) + ggplot2::xlab(NULL) + ggplot2::ylim(c(0, 1)) + ggplot2::ylab(NULL) + ggplot2::theme_bw() + ggplot2::theme( axis.text = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), panel.grid.major.y = ggplot2::element_blank(), axis.title = ggplot2::element_text(size = 20) ) + ggplot2::geom_line(ggplot2::aes( x = ordered(viData[order(-est), "var"]), y = InterpolationCutoff / maxVal ), color = "orange", linetype = "dashed", linewidth = 2) + ggplot2::ylab(paste0( "OOB (", specify_decimal(accData$OOB, 2), "), Guess (", specify_decimal(accData$guess, 2), "), Bias (", specify_decimal(accData$bias, 2), ")" )) png(paste0(save_path,"/paper/noeffect_variableimportance.png")) print(plot_full) dev.off() # Plot Subset InterpolationCutoff <- InterpolationCutoff[1:subsetPlotSize] viData <- viData[order(-viData$est), ] viData <- viData[1:subsetPlotSize, ] maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est) plot_top25 <- ggplot2::ggplot( data = viData, mapping = ggplot2::aes( x = factor(stats::reorder(var, est)), y = ifelse(est / maxVal > 1, 1, ifelse(est / maxVal < 0, 0, est / maxVal ) ), group = 1 ) ) + ggplot2::geom_errorbar( ggplot2::aes( ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal), ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal) ), color = "black", width = 0.2 ) + ggplot2::geom_point(color = "black", size = 3) + ggplot2::geom_hline( ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))), color = "red", linetype = "dotted", linewidth = 2 ) + ggplot2::coord_flip(ylim = c(0, 1)) + ggplot2::xlab(NULL) + ggplot2::ylim(c(0, 1)) + ggplot2::ylab(NULL) + ggplot2::theme_bw() + ggplot2::theme( axis.text = ggplot2::element_text(size = 18), axis.text.x = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), panel.grid.major.y = ggplot2::element_blank(), axis.title = ggplot2::element_text(size = 20) ) + ggplot2::geom_line(ggplot2::aes( x = ordered(viData[order(-est), "var"]), y = InterpolationCutoff / maxVal ), color = "orange", linetype = "dashed", linewidth = 2) + ggplot2::ylab(paste0( "OOB (", specify_decimal(accData$OOB, 2), "), Guess (", specify_decimal(accData$guess, 2), "), Bias (", specify_decimal(accData$bias, 2), ")" )) png(paste0(save_path,"/paper/noeffect_variableimportance_top25.png")) print(plot_top25) dev.off()
Now we consider the simulated model with effects.
set.seed(123456) dat0 <- simulatePP( agentVarData = data.frame( "outcome" = c(0, 1), "c1" = c(0, 0), "c2" = c(1 / 25, 1 / 60), "c3" = c(1 / 50, 1 / 10), "c4" = c(0, 0), "c5" = c(0, 0), "c6" = c(0, 0), "c7" = c(0, 0), "c8" = c(1 / 100, 1 / 100), "c9" = c(1 / 20, 1 / 20), "c10" = c(1 / 250, 1 / 250), "c11" = c(1 / 100, 1 / 100), "c12" = c(1 / 80, 1 / 80), "c13" = c(0, 0), "c14" = c(0, 0), "c15" = c(0, 0), "c16" = c(1 / 10, 1 / 10) ), agentKappaData = data.frame( "agent" = paste0("c", 1:16), "clusterAgent" = c(NA, "c1", "c1", rep(NA, 12), "c1"), "kappa" = c( rbinom(1, 100, 0.5), rbinom(1, 100, 0.5), rbinom(1, 30, 0.5), rbinom(1, 80, 0.5), rbinom(1, 350, 0.5), rbinom(1, 100, 0.5), rbinom(1, 120, 0.5), rbinom(1, 150, 0.5), rbinom(2, 250, 0.5), rbinom(1, 600, 0.5), rbinom(1, 60, 0.5), rbinom(2, 140, 0.5), rbinom(1, 20, 0.5), rbinom(1, 5, 0.5) ) ), unitsPerOutcome = 15, replicatesPerUnit = 1, silent = F ) dat1 <- simulatePP( agentVarData = data.frame( "outcome" = c(0, 1), "c1" = c(0, 0), "c2" = c(1 / 25, 1 / 25), "c3" = c(1 / 50, 1 / 50), "c4" = c(0, 0), "c5" = c(0, 0), "c6" = c(0, 0), "c7" = c(0, 0), "c8" = c(1 / 100, 1 / 100), "c9" = c(1 / 20, 1 / 20), "c10" = c(1 / 250, 1 / 250), "c11" = c(1 / 100, 1 / 100), "c12" = c(1 / 80, 1 / 80), "c13" = c(0, 0), "c14" = c(0, 0), "c15" = c(0, 0), "c16" = c(1 / 10, 1 / 10) ), agentKappaData = data.frame( "agent" = paste0("c", 1:16), "clusterAgent" = c(NA, "c1", "c1", rep(NA, 12), "c1"), "kappa" = c( rbinom(1, 100, 0.5), rbinom(1, 100, 0.5), rbinom(1, 30, 0.5), rbinom(1, 80, 0.5), rbinom(1, 350, 0.5), rbinom(1, 100, 0.5), rbinom(1, 120, 0.5), rbinom(1, 150, 0.5), rbinom(2, 250, 0.5), rbinom(1, 600, 0.5), rbinom(1, 60, 0.5), rbinom(2, 140, 0.5), rbinom(1, 20, 0.5), rbinom(1, 5, 0.5) ) ), unitsPerOutcome = 2, replicatesPerUnit = 1, silent = F ) dat1$unit <- ifelse(dat1$unit == "u1", "u31", ifelse(dat1$unit == "u2", "u32", ifelse(dat1$unit == "u3", "u33", ifelse(dat1$unit == "u4", "u34", NA) ) ) ) dat1$replicate <- ifelse(dat1$replicate == "1", "31", ifelse(dat1$replicate == "2", "32", ifelse(dat1$replicate == "3", "33", ifelse(dat1$replicate == "4", "34", NA) ) ) ) dat <- rbind(dat0, dat1) cells <- paste0("c", 1:16) cells_interactions <- rbind( data.frame(t(combn(cells, 2))), data.frame("X1" = cells, "X2" = cells) ) pcaData <- getKsPCAData( data = dat, replicate = "replicate", agents_df = cells_interactions, xRange = c(0, 1), yRange = c(0, 1), silent = F ) pcaMeta <- simulateMeta(pcaData, metaInfo = data.frame( "var" = c("age"), "rdist" = c("rnorm"), "outcome_0" = c("25"), "outcome_1" = c("27") ) ) rfcv <- funkyModel( data = pcaMeta, outcome = "outcome", unit = "unit", metaNames = c("age") )
Creating the variable importance plots as before (code below, perhaps figure not given for computational reasons)
rfcv$viPlot
rfcv$subset_viPlot
# Get Vars viData <- rfcv$VariableImportance accData <- rfcv$AccuracyEstimate NoiseCutoff <- rfcv$NoiseCutoff InterpolationCutoff <- rfcv$InterpolationCutoff subsetPlotSize <- rfcv$AdditionalParams$subsetPlotSize # Plot Full maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est) plot_full <- ggplot2::ggplot( data = viData, mapping = ggplot2::aes( x = factor(stats::reorder(var, est)), y = ifelse(est / maxVal > 1, 1, ifelse(est / maxVal < 0, 0, est / maxVal ) ), group = 1 ) ) + ggplot2::geom_errorbar( ggplot2::aes( ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal), ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal) ), color = "black", width = 0.2 ) + ggplot2::geom_point(color = "black", size = 3) + ggplot2::geom_hline( ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))), color = "red", linetype = "dotted", linewidth = 2 ) + ggplot2::coord_flip(ylim = c(0, 1)) + ggplot2::xlab(NULL) + ggplot2::ylim(c(0, 1)) + ggplot2::ylab(NULL) + ggplot2::theme_bw() + ggplot2::theme( axis.text = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), panel.grid.major.y = ggplot2::element_blank(), axis.title = ggplot2::element_text(size = 20) ) + ggplot2::geom_line(ggplot2::aes( x = ordered(viData[order(-est), "var"]), y = InterpolationCutoff / maxVal ), color = "orange", linetype = "dashed", linewidth = 2) + ggplot2::ylab(paste0( "OOB (", specify_decimal(accData$OOB, 2), "), Guess (", specify_decimal(accData$guess, 2), "), Bias (", specify_decimal(accData$bias, 2), ")" )) png(paste0(save_path,"/paper/effect_variableimportance.png")) print(plot_full) dev.off() # Plot Subset InterpolationCutoff <- InterpolationCutoff[1:subsetPlotSize] viData <- viData[order(-viData$est), ] viData <- viData[1:subsetPlotSize, ] maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est) plot_top25 <- ggplot2::ggplot( data = viData, mapping = ggplot2::aes( x = factor(stats::reorder(var, est)), y = ifelse(est / maxVal > 1, 1, ifelse(est / maxVal < 0, 0, est / maxVal ) ), group = 1 ) ) + ggplot2::geom_errorbar( ggplot2::aes( ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal), ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal) ), color = "black", width = 0.2 ) + ggplot2::geom_point(color = "black", size = 3) + ggplot2::geom_hline( ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))), color = "red", linetype = "dotted", linewidth = 2 ) + ggplot2::coord_flip(ylim = c(0, 1)) + ggplot2::xlab(NULL) + ggplot2::ylim(c(0, 1)) + ggplot2::ylab(NULL) + ggplot2::theme_bw() + ggplot2::theme( axis.text = ggplot2::element_text(size = 18), axis.text.x = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), panel.grid.major.y = ggplot2::element_blank(), axis.title = ggplot2::element_text(size = 20) ) + ggplot2::geom_line(ggplot2::aes( x = ordered(viData[order(-est), "var"]), y = InterpolationCutoff / maxVal ), color = "orange", linetype = "dashed", linewidth = 2) + ggplot2::ylab(paste0( "OOB (", specify_decimal(accData$OOB, 2), "), Guess (", specify_decimal(accData$guess, 2), "), Bias (", specify_decimal(accData$bias, 2), ")" )) png(paste0(save_path,"/paper/effect_variableimportance_top25.png")) print(plot_top25) dev.off()
With this confidence, we now consider the TNBC phenotype data.
dataPCA_pheno <- getKsPCAData( data = TNBC_pheno, unit = "Person", agents_df = pheno_interactions, rCheckVals = seq(0, 50, 1) ) dataPCAAge_pheno <- merge(dataPCA_pheno, TNBC_meta) set.seed(123456) rfcv <- funkyModel( data = dataPCAAge_pheno, K = 10, outcome = "Class", unit = "Person", metaNames = c("Age"), synthetics = 100, alpha = 0.05, silent = FALSE, subsetPlotSize = 25 )
And the related variable importance plots (may not be shown for computational reasons).
rfcv$viPlot
rfcv$subset_viPlot
# Get Vars viData <- rfcv$VariableImportance accData <- rfcv$AccuracyEstimate NoiseCutoff <- rfcv$NoiseCutoff InterpolationCutoff <- rfcv$InterpolationCutoff subsetPlotSize <- rfcv$AdditionalParams$subsetPlotSize # Plot Full maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est) plot_full <- ggplot2::ggplot( data = viData, mapping = ggplot2::aes( x = factor(stats::reorder(var, est)), y = ifelse(est / maxVal > 1, 1, ifelse(est / maxVal < 0, 0, est / maxVal ) ), group = 1 ) ) + ggplot2::geom_errorbar( ggplot2::aes( ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal), ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal) ), color = "black", width = 0.2 ) + ggplot2::geom_point(color = "black", size = 3) + ggplot2::geom_hline( ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))), color = "red", linetype = "dotted", linewidth = 2 ) + ggplot2::coord_flip(ylim = c(0, 1)) + ggplot2::xlab(NULL) + ggplot2::ylim(c(0, 1)) + ggplot2::ylab(NULL) + ggplot2::theme_bw() + ggplot2::theme( # axis.text = ggplot2::element_text(size = 18), axis.text = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), panel.grid.major.y = ggplot2::element_blank(), axis.title = ggplot2::element_text(size = 20) ) + ggplot2::geom_line(ggplot2::aes( x = ordered(viData[order(-est), "var"]), y = InterpolationCutoff / maxVal ), color = "orange", linetype = "dashed", linewidth = 2) + ggplot2::ylab(paste0( "OOB (", specify_decimal(accData$OOB, 2), "), Guess (", specify_decimal(accData$guess, 2), "), Bias (", specify_decimal(accData$bias, 2), ")" )) png(paste0(save_path,"/paper/tnbc_variableimportance.png")) print(plot_full) dev.off() # Plot Subset InterpolationCutoff <- InterpolationCutoff[1:subsetPlotSize] viData <- viData[order(-viData$est), ] viData <- viData[1:subsetPlotSize, ] maxVal <- max(InterpolationCutoff, NoiseCutoff, viData$est) plot_top25 <- ggplot2::ggplot( data = viData, mapping = ggplot2::aes( x = factor(stats::reorder(var, est)), y = ifelse(est / maxVal > 1, 1, ifelse(est / maxVal < 0, 0, est / maxVal ) ), group = 1 ) ) + ggplot2::geom_errorbar( ggplot2::aes( ymin = ifelse((est - sd) / maxVal < 0, 0, (est - sd) / maxVal), ymax = ifelse((est + sd) / maxVal > 1, 1, (est + sd) / maxVal) ), color = "black", width = 0.2 ) + ggplot2::geom_point(color = "black", size = 3) + ggplot2::geom_hline( ggplot2::aes(yintercept = max(0, min(1, NoiseCutoff / maxVal))), color = "red", linetype = "dotted", linewidth = 2 ) + ggplot2::coord_flip(ylim = c(0, 1)) + ggplot2::xlab(NULL) + ggplot2::ylim(c(0, 1)) + ggplot2::ylab(NULL) + ggplot2::theme_bw() + ggplot2::theme( axis.text.x = ggplot2::element_blank(), axis.text.y = ggplot2::element_text(size = 14), axis.ticks.y = ggplot2::element_blank(), panel.grid.major.y = ggplot2::element_blank(), axis.title = ggplot2::element_text(size = 17) ) + ggplot2::geom_line(ggplot2::aes( x = ordered(viData[order(-est), "var"]), y = InterpolationCutoff / maxVal ), color = "orange", linetype = "dashed", linewidth = 2) + ggplot2::ylab(paste0( "OOB (", specify_decimal(accData$OOB, 2), "), Guess (", specify_decimal(accData$guess, 2), "), Bias (", specify_decimal(accData$bias, 2), ")" )) png(paste0(save_path,"/paper/tnbc_variableimportance_top25.png")) print(plot_top25) dev.off()
Also consider $K$ functions from significant and insignificant interactions.
tmp <- getKFunction(TNBC_pheno[TNBC_pheno$Class == 0, -1], c("Tumor", "Tumor"), unit = "Person", rCheckVals = seq(0, 50, 1) ) tmp1 <- getKFunction(TNBC_pheno[TNBC_pheno$Class == 1, -1], c("Tumor", "Tumor"), unit = "Person", rCheckVals = seq(0, 50, 1) ) tmp_1 <- tidyr::pivot_longer(data = tmp, cols = -r) tmp1_1 <- tidyr::pivot_longer(data = tmp1, cols = -r) data_plot <- rbind( data.frame( "r" = tmp_1$r, "K" = tmp_1$value, "unit" = tmp_1$name, "outcome" = "0" ), data.frame( "r" = tmp1_1$r, "K" = tmp1_1$value, "unit" = paste0(tmp1_1$name, "_1"), "outcome" = "1" ) )
Creating the following figure.
plot_K_functions(data_plot)
data <- data_plot # Prep Legend info <- unique(data[, c("unit", "outcome")]) info$Missing <- NA for (i in 1:nrow(info)) { info[i, 3] <- nrow(data[data$unit == info$unit[i] & !stats::complete.cases(data), ]) > 0 } info_labels <- data.frame( "outcome" = unique(data$outcome), "units" = NA, "Missing" = NA, "Label" = NA ) for (i in 1:nrow(info_labels)) { info_labels[i, 2] <- nrow(info[info$outcome == info_labels[i, 1], ]) info_labels[i, 3] <- nrow(info[info$outcome == info_labels[i, 1] & info$Missing, ]) # Only label outcomes with numbers if at least one is missing in whole set if (sum(info$Missing) != 0) { info_labels[i, 4] <- paste0( info_labels[i, 1], " (", info_labels[i, 2] - info_labels[i, 3], "/", info_labels[i, 2], ")" ) } else { info_labels[i, 4] <- info_labels[i, 1] } } # Build Averages data_wide <- tidyr::pivot_wider(data, names_from = "unit", values_from = "K" ) data_avg <- data.frame("r" = unique(data_wide$r)) for (i in 1:nrow(info_labels)) { data_avg[[info_labels[i, "outcome"]]] <- rowMeans(data_wide[data_wide$outcome == info_labels[i, "outcome"], -c(1:2)], na.rm = TRUE ) } data_avg <- tidyr::pivot_longer(data_avg, cols = -r) colnames(data_avg) <- c("r", "outcome", "Value") # Plot return_plot <- ggplot2::ggplot( data = stats::na.omit(data), ggplot2::aes( x = r, y = K, group = unit, color = outcome ) ) + ggplot2::geom_line(alpha = 0.2, linewidth = 1.25) + ggplot2::geom_line( ggplot2::aes( x = r, y = Value, group = outcome, color = outcome ), data = data_avg, linewidth = 2 ) + ggplot2::theme_bw() + ggplot2::theme( axis.text = ggplot2::element_text(size = 18), axis.text.y = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), panel.grid.major.y = ggplot2::element_blank(), # axis.title = ggplot2::element_text(size = 22), axis.title = ggplot2::element_blank(), legend.position = "none" ) # + # ggplot2::geom_line(ggplot2::aes(x = r, y = pi * r^2, # group = outcome), # data = data_avg, linewidth = 1.25, # color = "black", linetype = "dashed" # ) png(paste0(save_path,"/paper/TNBC_sig_KFunctions.png")) print(return_plot) dev.off()
And for a insignificant interaction.
tmp <- getKFunction(TNBC_pheno[TNBC_pheno$Class == 0, -1], c("CD4T", "Endothelial"), unit = "Person", rCheckVals = seq(0, 50, 1) ) tmp1 <- getKFunction(TNBC_pheno[TNBC_pheno$Class == 1, -1], c("CD4T", "Endothelial"), unit = "Person", rCheckVals = seq(0, 50, 1) ) tmp_1 <- tidyr::pivot_longer(data = tmp, cols = -r) tmp1_1 <- tidyr::pivot_longer(data = tmp1, cols = -r) data_plot <- rbind( data.frame( "r" = tmp_1$r, "K" = tmp_1$value, "unit" = tmp_1$name, "outcome" = "0" ), data.frame( "r" = tmp1_1$r, "K" = tmp1_1$value, "unit" = paste0(tmp1_1$name, "_1"), "outcome" = "1" ) )
Creating the following figure.
plot_K_functions(data_plot)
data <- data_plot # Prep Legend info <- unique(data[, c("unit", "outcome")]) info$Missing <- NA for (i in 1:nrow(info)) { info[i, 3] <- nrow(data[data$unit == info$unit[i] & !stats::complete.cases(data), ]) > 0 } info_labels <- data.frame( "outcome" = unique(data$outcome), "units" = NA, "Missing" = NA, "Label" = NA ) for (i in 1:nrow(info_labels)) { info_labels[i, 2] <- nrow(info[info$outcome == info_labels[i, 1], ]) info_labels[i, 3] <- nrow(info[info$outcome == info_labels[i, 1] & info$Missing, ]) # Only label outcomes with numbers if at least one is missing in whole set if (sum(info$Missing) != 0) { info_labels[i, 4] <- paste0( info_labels[i, 1], " (", info_labels[i, 2] - info_labels[i, 3], "/", info_labels[i, 2], ")" ) } else { info_labels[i, 4] <- info_labels[i, 1] } } # Build Averages data_wide <- tidyr::pivot_wider(data, names_from = "unit", values_from = "K" ) data_avg <- data.frame("r" = unique(data_wide$r)) for (i in 1:nrow(info_labels)) { data_avg[[info_labels[i, "outcome"]]] <- rowMeans(data_wide[data_wide$outcome == info_labels[i, "outcome"], -c(1:2)], na.rm = TRUE ) } data_avg <- tidyr::pivot_longer(data_avg, cols = -r) colnames(data_avg) <- c("r", "outcome", "Value") # Plot return_plot <- ggplot2::ggplot( data = stats::na.omit(data), ggplot2::aes( x = r, y = K, group = unit, color = outcome ) ) + ggplot2::geom_line(alpha = 0.2, linewidth = 1.25) + ggplot2::geom_line( ggplot2::aes( x = r, y = Value, group = outcome, color = outcome ), data = data_avg, linewidth = 2 ) + ggplot2::theme_bw() + ggplot2::theme( axis.text = ggplot2::element_text(size = 18), axis.text.y = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), panel.grid.major.y = ggplot2::element_blank(), # axis.title = ggplot2::element_text(size = 22), axis.title = ggplot2::element_blank(), legend.position = "none" ) # + # ggplot2::geom_line(ggplot2::aes(x = r, y = pi * r^2, # group = outcome), # data = data_avg, linewidth = 1.25, # color = "black", linetype = "dashed" # ) png(paste0(save_path,"/paper/TNBC_insig_KFunctions.png")) print(return_plot) dev.off()
Further analysis of this data can be found in our paper.
Try the funkycells 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.