## ----include = TRUE, echo = TRUE, eval=FALSE----------------------------------
#  
#  # bioconductor install
#  if (!requireNamespace("BiocManager", quietly = TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("IntOMICS")
#  
#  # install the newest (development) version from GitHub
#  # install.packages("remotes")
#  remotes::install_github("anna-pacinkova/IntOMICS")
#  

## ----include = TRUE, echo = TRUE----------------------------------------------

library(IntOMICS)

# load required libraries
suppressPackageStartupMessages({
library(curatedTCGAData)
library(TCGAutils)
library(bnlearn)
library(bnstruct)
library(matrixStats)
library(parallel)
library(RColorBrewer)
library(bestNormalize)
library(igraph)
library(gplots)
library(methods)
library(ggraph)
library(ggplot2)})


## ----echo=TRUE, include=TRUE, message = FALSE---------------------------------

coad_SE <- curatedTCGAData("COAD",
                           c("RNASeq2GeneNorm_illuminahiseq",
                             "Methylation_methyl450","GISTIC_AllByGene"), 
                           version = "2.0.1", dry.run = FALSE)
coad_SE <- TCGAprimaryTumors(coad_SE)
 
## keep NA for Methylation data
coad_ma <- subsetByColData(coad_SE, coad_SE$MSI_status == "MSI-H" | 
                             is.na(coad_SE$MSI_status))
# choose random 50 samples for illustration
random_samples <- sample(names(which(table(sampleMap(coad_ma)$primary)==3)),50)
coad_ma <- subsetByColData(coad_ma, random_samples)

rownames(coad_ma[["COAD_GISTIC_AllByGene-20160128"]]) <- rowData(coad_ma[["COAD_GISTIC_AllByGene-20160128"]])[["Gene.Symbol"]]
 
data(list=c("gene_annot","annot"))
rowselect <- list(gene_annot$gene_symbol, gene_annot$gene_symbol, 
                  unlist(annot))
names(rowselect) <- names(coad_ma)
omics <- coad_ma[rowselect, , ]
names(omics) <- c("cnv","ge","meth")


## ----echo=TRUE, include=TRUE--------------------------------------------------

t(assay(omics[["ge"]]))[1:5,1:5]


## ----echo=TRUE, include=TRUE--------------------------------------------------

t(assay(omics[["cnv"]]))[1:5,1:5]


## ----echo=TRUE, include=TRUE--------------------------------------------------

t(assay(omics[["meth"]]))[1:5,1:5]


## ----echo=TRUE, include=TRUE--------------------------------------------------

str(annot)


## ----echo=TRUE, include=TRUE--------------------------------------------------

gene_annot


## ----echo=TRUE, include=TRUE--------------------------------------------------

data("PK")
PK


## ----echo=TRUE, include=TRUE--------------------------------------------------

data("TFtarg_mat")
TFtarg_mat[1:5,1:5]


## ----echo=TRUE, include=TRUE--------------------------------------------------

data("layers_def")
layers_def


## ----echo=TRUE, include=TRUE--------------------------------------------------

OMICS_mod_res <- omics_module(omics = omics, 
                              PK = PK, 
                              layers_def = layers_def, 
                              TFtargs = TFtarg_mat,
                              annot = annot, 
                              gene_annot = gene_annot,
                              lm_METH = TRUE,
                              r_squared_thres = 0.3,
                              p_val_thres = 0.1)


## ----echo=TRUE, include=TRUE--------------------------------------------------

names(OMICS_mod_res)


## ----echo=TRUE, include=TRUE--------------------------------------------------

if(interactive())
{
  BN_mod_res_sparse <- bn_module(burn_in = 500, 
                               thin = 50, 
                               OMICS_mod_res = OMICS_mod_res,
                               minseglen = 5)
}


## ----echo=TRUE, include=TRUE--------------------------------------------------

data("BN_mod_res")
getSlots(class(BN_mod_res))


## ----echo=TRUE, include=TRUE--------------------------------------------------

trace_plots(mcmc_res = BN_mod_res,
            burn_in = 10000,
            thin = 500, 
            edge_freq_thres = 0.5)


## ----echo=TRUE, include=TRUE--------------------------------------------------

res_weighted <- edge_weights(mcmc_res = BN_mod_res, 
                            burn_in = 10000, 
                            thin = 500, 
                            edge_freq_thres = 0.5)

weighted_net_res <- weighted_net(cpdag_weights = res_weighted,
                                 gene_annot = gene_annot, 
                                 PK = PK, 
                                 OMICS_mod_res = OMICS_mod_res, 
                                 gene_ID = "gene_symbol", 
                                 TFtargs = TFtarg_mat, 
                                 B_prior_mat_weighted = B_prior_mat_weighted(BN_mod_res))


## ----fig_weights, fig.height = 7, fig.width = 7-------------------------------

ggraph_weighted_net(net = weighted_net_res, 
                    node_size = 10, 
                    node_label_size = 4, 
                    edge_label_size = 4)


## ----echo=TRUE, include=TRUE--------------------------------------------------

weighted_net_res <- weighted_net(cpdag_weights = res_weighted,
                                 gene_annot = gene_annot, 
                                 PK = PK, 
                                 OMICS_mod_res = OMICS_mod_res, 
                                 gene_ID = "gene_symbol", 
                                 edge_weights = "empB",
                                 TFtargs = TFtarg_mat, 
                                 B_prior_mat_weighted = B_prior_mat_weighted(BN_mod_res))


## ----fig_empB, fig.height = 7, fig.width = 7----------------------------------

ggraph_weighted_net(net = weighted_net_res)


## ----fig_heat, fig.height = 7, fig.width = 7----------------------------------

emp_b_heatmap(mcmc_res = BN_mod_res, 
             OMICS_mod_res = OMICS_mod_res, 
             gene_annot = gene_annot, 
             TFtargs = TFtarg_mat)


## ----fig_ew, fig.height = 4.5, fig.width = 5----------------------------------

res_weighted <- edge_weights(mcmc_res = BN_mod_res, 
                            burn_in = 10000, 
                            thin = 500,
                            edge_freq_thres = NULL)
                            
weighted_net_res <- weighted_net(cpdag_weights = res_weighted,
                                 gene_annot = gene_annot, 
                                 PK = PK, 
                                 OMICS_mod_res = OMICS_mod_res, 
                                 gene_ID = "gene_symbol", 
                                 TFtargs = TFtarg_mat, 
                                 B_prior_mat_weighted = B_prior_mat_weighted(BN_mod_res))

dens_edge_weights(weighted_net_res)


## ----echo=TRUE, include=TRUE--------------------------------------------------

sessionInfo()