## ----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()