## ----eval=FALSE---------------------------------------------------------------
#  if (!require("BiocManager", quietly = TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("dce")

## ----message=FALSE------------------------------------------------------------
# fix "object 'guide_edge_colourbar' of mode 'function' was not found"
# when building vignettes
# (see also https://github.com/thomasp85/ggraph/issues/75)
library(ggraph)

library(curatedTCGAData)
library(TCGAutils)
library(SummarizedExperiment)

library(tidyverse)
library(cowplot)
library(graph)
library(dce)

set.seed(42)

## -----------------------------------------------------------------------------
graph_wt <- matrix(c(0, 0, 0, 1, 0, 0, 1, 1, 0), 3, 3)
rownames(graph_wt) <- colnames(graph_wt) <- c("A", "B", "C")
graph_wt

## -----------------------------------------------------------------------------
graph_mt <- graph_wt
graph_mt["A", "B"] <- 2.5 # dysregulation happens here!
graph_mt

cowplot::plot_grid(
  plot_network(graph_wt, edgescale_limits = c(-3, 3)),
  plot_network(graph_mt, edgescale_limits = c(-3, 3)),
  labels = c("WT", "MT")
)

## -----------------------------------------------------------------------------
X_wt <- simulate_data(graph_wt)
X_mt <- simulate_data(graph_mt)

X_wt %>%
  head

## -----------------------------------------------------------------------------
res <- dce(graph_wt, X_wt, X_mt)

res %>%
  as.data.frame %>%
  drop_na

## -----------------------------------------------------------------------------
plot(res) +
  ggtitle("Differential Causal Effects between WT and MT condition")

## -----------------------------------------------------------------------------
set.seed(1337)

# create wild-type and mutant networks
graph_wt <- create_random_DAG(20, 0.3)
graph_mt <- resample_edge_weights(graph_wt)

cowplot::plot_grid(
  plot_network(as(graph_wt, "matrix"), labelsize = 0, arrow_size = 0.01),
  plot_network(as(graph_mt, "matrix"), labelsize = 0, arrow_size = 0.01),
  labels = c("WT", "MT")
)

## -----------------------------------------------------------------------------
# simulate gene expression data for both networks
X_wt <- simulate_data(graph_wt)
X_mt <- simulate_data(graph_mt)

# compute DCEs
res <- dce::dce(graph_wt, X_wt, X_mt)

df_dce <- res %>%
  as.data.frame %>%
  drop_na %>%
  arrange(dce_pvalue)

## -----------------------------------------------------------------------------
# compute ground truth DCEs
dce_gt <- trueEffects(graph_mt) - trueEffects(graph_wt)
dce_gt_ind <- which(dce_gt != 0, arr.ind = TRUE)

# create plot
data.frame(
  source = paste0("n", dce_gt_ind[, "row"]),
  target = paste0("n", dce_gt_ind[, "col"]),
  dce_ground_truth = dce_gt[dce_gt != 0]
) %>%
  inner_join(df_dce, by = c("source", "target")) %>%
  rename(dce_estimate = dce) %>%
ggplot(aes(x = dce_ground_truth, y = dce_estimate)) +
  geom_abline(color = "gray") +
  geom_point() +
  xlab("DCE (ground truth") +
  ylab("DCE (estimate)") +
  theme_minimal()

## -----------------------------------------------------------------------------
brca <- curatedTCGAData(
  diseaseCode = "BRCA",
  assays = c("RNASeq2*"),
  version = "2.0.1",
  dry.run = FALSE
)

## -----------------------------------------------------------------------------
sampleTables(brca)

data(sampleTypes, package = "TCGAutils")
sampleTypes %>%
  dplyr::filter(Code %in% c("01", "06", "11"))

## -----------------------------------------------------------------------------
# split assays
brca_split <- TCGAsplitAssays(brca, c("01", "11"))

# only retain matching samples
brca_matched <- as(brca_split, "MatchedAssayExperiment")

brca_wt <- assay(brca_matched, "01_BRCA_RNASeq2GeneNorm-20160128")
brca_mt <- assay(brca_matched, "11_BRCA_RNASeq2GeneNorm-20160128")

## -----------------------------------------------------------------------------
pathways <- get_pathways(pathway_list = list(kegg = c("Breast cancer")))
brca_pathway <- pathways[[1]]$graph

## -----------------------------------------------------------------------------
shared_genes <- intersect(nodes(brca_pathway), rownames(brca_wt))
glue::glue(
  "Covered nodes: {length(shared_genes)}/{length(nodes(brca_pathway))}"
)

## ----warning=FALSE------------------------------------------------------------
res <- dce::dce(brca_pathway, t(brca_wt), t(brca_mt))

## -----------------------------------------------------------------------------
res %>%
  as.data.frame %>%
  drop_na %>%
  arrange(desc(abs(dce))) %>%
  head

plot(
  res,
  nodesize = 20, labelsize = 1,
  node_border_size = 0.05, arrow_size = 0.01,
  use_symlog = TRUE,
  shadowtext = TRUE
)

## -----------------------------------------------------------------------------
set.seed(1)
epsilon <- 1e-100
network_size <- 50
graph_wt <- as(create_random_DAG(network_size, prob = .2), "matrix")
graph_wt["n1", "n2"] <- epsilon
graph_mt <- graph_wt
graph_mt["n1", "n2"] <- 2

cowplot::plot_grid(
  plot_network(graph_wt, edgescale_limits = c(-2, 2)),
  plot_network(graph_mt, edgescale_limits = c(-2, 2)),
  labels = c("WT", "MT")
)

truth <- trueEffects(graph_mt) - trueEffects(graph_wt)
plot_network(graph_wt, value_matrix = truth, edgescale_limits = c(-2, 2))

X_wt <- simulate_data(n = 100, graph_wt)
X_mt <- simulate_data(n = 100, graph_mt)

## -----------------------------------------------------------------------------
res <- dce(graph_mt, X_wt, X_mt, deconfounding = FALSE)

qplot(truth[truth != 0], res$dce[truth != 0]) +
  geom_abline(color = "red", linetype = "dashed") +
  xlab("true DCE") +
  ylab("estimated DCE")

plot_network(graph_wt, value_matrix = -log(res$dce_pvalue)) +
  ggtitle("-log(p-values) for DCEs between WT and MT condition")

## -----------------------------------------------------------------------------
batch <- sample(c(0, 1), replace = TRUE, nrow(X_wt))
bX_wt <- apply(X_wt, 2, function(x) x + max(x) * runif(1) * batch)
bX_mt <- apply(X_mt, 2, function(x) x + max(x) * runif(1) * batch)

res_without_deconf <- dce(graph_mt, bX_wt, bX_mt, deconfounding = FALSE)

cowplot::plot_grid(
  plot_network(
    graph_wt,
    value_matrix = -log(res_without_deconf$dce_pvalue + epsilon)
  ) +
    ggtitle("-log(p-values) without deconfounding"),
  qplot(truth[truth != 0], res_without_deconf$dce[truth != 0]) +
    geom_abline(color = "red", linetype = "dashed") +
    xlab("true DCE") +
    ylab("estimated DCE"),
  nrow = 1
)

## -----------------------------------------------------------------------------
res_with_deconf <- dce(graph_mt, bX_wt, bX_mt, deconfounding = 1)

cowplot::plot_grid(
  plot_network(
    graph_wt,
    value_matrix = -log(res_with_deconf$dce_pvalue + epsilon)
  ) +
    ggtitle("-log(p-values) with deconfounding"),
  qplot(truth[truth != 0], res_with_deconf$dce[truth != 0]) +
    geom_abline(color = "red", linetype = "dashed") +
    xlab("true DCE") +
    ylab("estimated DCE"),
  nrow = 1
)

## -----------------------------------------------------------------------------
sessionInfo()

