Upload New File

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))
+}