Delete trace_plots.R

source_code/trace_plots.R

-# edge_freq_thres is the edge weights quantile, which is the minimal edge weight accepted in the resulting network
-#' @export
-trace_plots <- function(mcmc_res, burn_in, thin, figures_dir, gene_annot, PK, OMICS_module_res, edge_weights = "MCMC_freq", edge_freq_thres = NULL, gene_ID, TFtargs)
-{
-  if(!(edge_weights %in% c("MCMC_freq","empB")))
-  {
-    print('Warning: edge_weights argument must be either "MCMC_freq" or "empB".')  
-  }
-  
-  if(!dir.exists(figures_dir)){dir.create(figures_dir)}
-
-  # Trace plot of beta values
-  df1 <- data.frame(beta = mapply(mcmc_res$beta_tuning,FUN=function(x) x$value), k=1:length(mapply(mcmc_res$beta_tuning,FUN=function(x) x$value)), accept = 1)
-
-  # RMS strength (threshold)
-  rms_strength <- abs(diff(mcmc_res$sampling.phase_res$rms))
-  strength_threshold <- quantile(rms_strength, 0.75, na.rm = TRUE)
-
-  # custom.strength estimates the strength of each arc as its empirical frequency over a set of networks:
-  cpdags1 <- unique(mcmc_res$sampling.phase_res$mcmc_sim_part_res$seed1$cpdags[(burn_in/thin+1):length(mcmc_res$sampling.phase_res$mcmc_sim_part_res$seed1$cpdags)])
-  cpdags2 <- unique(mcmc_res$sampling.phase_res$mcmc_sim_part_res$seed2$cpdags[(burn_in/thin+1):length(mcmc_res$sampling.phase_res$mcmc_sim_part_res$seed2$cpdags)])
-
-  cpdag_weights1 <- custom.strength(cpdags1, nodes = bnlearn::nodes(cpdags1[[1]]), weights = NULL)
-  cpdag_weights2 <- custom.strength(cpdags2, nodes = bnlearn::nodes(cpdags2[[1]]), weights = NULL)
-  cpdag_weights1 <- cpdag_weights1[cpdag_weights1$direction>=0.5,]
-  cpdag_weights2 <- cpdag_weights2[cpdag_weights2$direction>=0.5,]
-
-  cpdag_weights1$edge <- paste(cpdag_weights1$from, cpdag_weights1$to, sep="_")
-  cpdag_weights2$edge <- paste(cpdag_weights2$from, cpdag_weights2$to, sep="_")
-  cpdag_weights <- merge(cpdag_weights1, cpdag_weights2, by = "edge")
-  cpdag_weights$strength <- round(rowMeans(cbind(cpdag_weights$strength.x, cpdag_weights$strength.y)),2)
-  
-  if(!is.null(edge_freq_thres))
-  {
-    strength_quant <- quantile(x = cpdag_weights$strength, probs = edge_freq_thres)
-    cpdag_weights <- cpdag_weights[cpdag_weights$strength >= strength_quant,]
-  }
-
-  total <- merge(cpdag_weights1, cpdag_weights2, by = c("from","to"))
-
-  svg(paste(figures_dir,"beta_values.svg",sep="/"))
-  plot(df1$beta ~ df1$k, type = "l", col= "darkblue",
-       main = "Beta values of adaptive MCMC",
-       xlab = "iteration",
-       ylab = "beta")
-  dev.off()
-
-  svg(paste(figures_dir,"post_prob_edges.svg",sep="/"))
-  plot(total$strength.x ~ total$strength.y,
-       main = "Consistency of edges posterior probabilities",
-       xlab="MCMC run 2",
-       ylab = "MCMC run 1")
-  abline(0,1, col="orange")
-  dev.off()
-
-  svg(paste(figures_dir,"convergence_RMS.svg",sep="/"))
-  plot(rms_strength, main="Convergence RMS strength (C.RMS.str)", pch = 18, col="gray30")
-  abline(h=strength_threshold, col="#E69F00", lwd = 1.5)
-  text(label = paste("3rd quartile of C.RMS.str = ", round(strength_threshold,3),sep=""), x = 100, y = strength_threshold+0.015, col="#E69F00")
-  dev.off()
-
-  PK <- PK[PK$src_entrez %in% unlist(lapply(OMICS_module_res$omics,colnames)),]
-  PK <- PK[PK$dest_entrez %in% unlist(lapply(OMICS_module_res$omics,colnames)),]
-  
-  if(gene_ID=="entrezID")
-  {
-    edge_list <- matrix(data = c(cpdag_weights$from.x,
-                                 cpdag_weights$to.x,
-                                 cpdag_weights$strength,
-                                 rep(NA,length(cpdag_weights$strength)),
-                                 rep(NA,length(cpdag_weights$strength))),
-                        nrow = length(cpdag_weights$strength),
-                        dimnames = list(c(), c("from", "to", "weight", "edge_type", "edge")))
-    # needs to be sorted because of colors in the final figure
-    node_list <- unique(c(edge_list[,"from"], edge_list[,"to"]))
-    edge_list[,"edge"] <- paste(edge_list[,"from"], edge_list[,"to"], sep="_")
-    return_list <- edge_types(mcmc_res = mcmc_res, PK = PK, gene_annot = gene_annot, edge_list = edge_list, node_list = node_list, OMICS_module_res = OMICS_module_res, edge_weights = edge_weights, TFtargs = TFtargs)
-  } else if (gene_ID=="gene_symbol") {
-    from <- as.character(gene_annot$gene_symbol[match(cpdag_weights$from.x, gene_annot$entrezID)])
-    from[is.na(from)] <- cpdag_weights$from.x[is.na(from)]
-    from[regexpr("entrezid",from)>0] <- tolower(as.character(gene_annot$gene_symbol[match(toupper(from[regexpr("entrezid",from)>0]), gene_annot$entrezID)]))
-    to <- as.character(gene_annot$gene_symbol[match(cpdag_weights$to.x, gene_annot$entrezID)])
-    edge_list <- matrix(data = c(from,
-                                 to,
-                                 cpdag_weights$strength,
-                                 rep(NA,length(cpdag_weights$strength)),
-                                 rep(NA,length(cpdag_weights$strength))),
-                        nrow = length(cpdag_weights$strength),
-                        dimnames = list(c(), c("from", "to", "weight", "edge_type", "edge")))
-    # needs to be sorted because of colors in the final figure
-    node_list <- unique(c(edge_list[,"from"], edge_list[,"to"]))
-    edge_list[,"edge"] <- paste(edge_list[,"from"], edge_list[,"to"], sep="_")
-    return_list <- edge_types(mcmc_res = mcmc_res, PK = PK, gene_annot = gene_annot, edge_list = edge_list, node_list = node_list, OMICS_module_res = OMICS_module_res, edge_weights = edge_weights, TFtargs = TFtargs)
-  } else {
-    print('gene_ID argument must be either "entrezID" or "gene_symbol"')
-  } # end if else if else (gene_ID=="entrezID")
-
-  net_weighted <- graph_from_edgelist(return_list$edge_list[,c("from","to")])
-  V(net_weighted)$color <- return_list$node_list[match(as_ids(V(net_weighted)),return_list$node_list[,"label"]),"color"]
-  palette <- return_list$node_palette
-  names(palette) <- 1:length(palette)
-  palette <- palette[unique(V(net_weighted)$color)]
-  V(net_weighted)$label <- return_list$node_list[match(as_ids(V(net_weighted)),return_list$node_list[,"label"]),"label"]
-  E(net_weighted)$edge <- return_list$edge_list[match(sub("|","_",as_ids(E(net_weighted)), fixed = TRUE),return_list$edge_list[,"edge"]),"edge_type"]
-  E(net_weighted)$weight <- return_list$edge_list[match(sub("|","_",as_ids(E(net_weighted)), fixed = TRUE),return_list$edge_list[,"edge"]),"weight"]
-  
-  # arrow in the network plot
-  V(net_weighted)$degree <- igraph::degree(net_weighted, mode = "in")
-  V(net_weighted)$degree <- normalise(V(net_weighted)$degree, to = c(3, 11))
-  
-  return(list(edge_list = return_list$edge_list,
-              node_list = return_list$node_list,
-              borders_GE = return_list$borders_GE,
-              borders_CNV = return_list$borders_CNV,
-              borders_METH = return_list$borders_METH,
-              node_palette = palette,
-              net_weighted = net_weighted))
-}