-
Notifications
You must be signed in to change notification settings - Fork 4
/
apply_ggm.R
executable file
·95 lines (77 loc) · 3.29 KB
/
apply_ggm.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
#!/usr/bin/env Rscript
# ------------------------------------------------------------------------------
#' Apply graphical model inference on data collected for a specific sentinel.
#'
#' @author Johann Hawe <[email protected]>
# ------------------------------------------------------------------------------
log <- file(snakemake@log[[1]], open="wt")
sink(log)
sink(log, type="message")
# ------------------------------------------------------------------------------
print("Prepare libraries and source scripts.")
# ------------------------------------------------------------------------------
library(pheatmap)
suppressPackageStartupMessages(library(GenomicRanges))
library(igraph)
library(graph)
library(reshape2)
source("scripts/lib.R")
source("scripts/reg_net.R")
source("scripts/reg_net_utils.R")
# ------------------------------------------------------------------------------
print("Get snakemake parameters.")
# ------------------------------------------------------------------------------
#input
franges <- snakemake@input[["ranges"]]
fdata <- snakemake@input[["data"]]
fppi_db <- snakemake@input[["ppi_db"]]
fpriors <- snakemake@input[["priors"]]
fcpg_context <- snakemake@input[["cpg_context"]]
ftss_context <- snakemake@input[["tss_context"]]
# output
fout <- snakemake@output$fit
fsummary_plot <- snakemake@output$summary_file
# params
threads <- snakemake@threads
# ------------------------------------------------------------------------------
print("Load and prepare data.")
# ------------------------------------------------------------------------------
data <- readRDS(fdata)
# remove (rare) all-NA cases. This can happen due to scaling of all-zero entities,
# which can arise due to a very large number of cis-meQTLs which effects get
# removed from the CpGs during data preprocessing.
# NOTE: we could possibly handle this differently? Seems that these particular
# cpgs are highly influenced by genetic effects?
use <- apply(data,2,function(x) (sum(is.na(x)) / length(x)) < 1)
data <- data[,use]
print("Dimensions of data:")
print(dim(data))
priors <- readRDS(fpriors)
# filter for available data in priros
priors <- priors[colnames(data), colnames(data)]
ranges <- readRDS(franges)
# load PPI DB
ppi_db <- readRDS(fppi_db)
# ------------------------------------------------------------------------------
# for catching the genenet summary plots, we open the pdf connection here
# ------------------------------------------------------------------------------
pdf(fsummary_plot)
# ------------------------------------------------------------------------------
print("Infer regulatory networks.")
# ------------------------------------------------------------------------------
if(ranges$seed == "meqtl") {
fcontext <- fcpg_context
} else {
fcontext <- ftss_context
}
result <- infer_all_graphs(data, priors, ranges, fcontext, ppi_db,
threads)
dev.off()
# ------------------------------------------------------------------------------
print("All done. Saving results.")
# ------------------------------------------------------------------------------
saveRDS(file=fout, result)
# ------------------------------------------------------------------------------
print("Session info:")
# ------------------------------------------------------------------------------
sessionInfo()