## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----eval=FALSE---------------------------------------------------------------
# if (!requireNamespace("BiocManager"))
#     install.packages("BiocManager")
# BiocManager::install("CBN2Path")

## -----------------------------------------------------------------------------
library(CBN2Path)

## ----include=FALSE------------------------------------------------------------
library(TCGAbiolinks)
library(BiocParallel)

## ----message=FALSE, include=FALSE---------------------------------------------
# Obtaining the "TCGA-BLCA" data
rawData <- getRawTCGAData("TCGA-BLCA")

# Specifying the genes
genes <- c("TP53","ARID1A","KDM6A","PIK3CA")

# Generating the genotype matrix
gMat <- generateTCGAMatrix(rawData, genes)

## -----------------------------------------------------------------------------
# The poset
dag <- matrix(c(1, 1, 1, 4, 2, 3, 4, 3), 4, 2)

# The genotype matrix: 
# gMat, which was generated using the generateTCGAMatrix function (see above). 

# Preparing input of the ct-cbn/h-cbn methods
bc <- Spock$new(
     poset = dag,
     numMutations = 4,
     genotypeMatrix = gMat
)

# Running the ct-cbn model
resultsC <- ctcbnSingle(bc)

## ----fig.width=4.25, fig.height=4.25------------------------------------------
visualizeCBNModel(dag)

## -----------------------------------------------------------------------------
mlLambda <- resultsC$lambda
logLikelihood <- resultsC$summary["Loglike"]

## -----------------------------------------------------------------------------
gMatMut <- genotypeMatrixMutator(gMat, 0.3, 0.2)

## -----------------------------------------------------------------------------
# The poset
dag <- matrix(c(1, 1, 1, 4, 2, 3, 4, 3), 4, 2)
# Preparing the inputs of the ct-cbn method
bc <- Spock$new(
     poset = dag,
     numMutations = 4,
     genotypeMatrix = gMatMut
)
# Running the ct-cbn model
resultsCMut <- ctcbnSingle(bc)

## -----------------------------------------------------------------------------
examplePath <- getExamples()[1]
bc <- Spock$new(
     poset = readPoset(examplePath)$sets,
     numMutations = readPoset(examplePath)$mutations,
     genotypeMatrix = readPattern(examplePath)
)
resultsC2 <- ctcbnSingle(bc)

## ----message=FALSE------------------------------------------------------------
posets <- readRDS(system.file("extdata", "Posets.rds", package = "CBN2Path"))

bc <-  Spock$new(
     poset = posets,
     numMutations = 4,
     genotypeMatrix = gMat
)
resultsC3 <- ctcbn(bc, nCores = 3)

## -----------------------------------------------------------------------------
logLik <- sapply(resultsC3, \(x) x$summary["Loglike"])

## -----------------------------------------------------------------------------
indx <- which.max(logLik)
mlPosetC <- posets[[indx]]

## ----fig.width=4.25, fig.height=4.25------------------------------------------
visualizeCBNModel(mlPosetC)

## -----------------------------------------------------------------------------
# The poset
dag <- matrix(c(1, 1, 1, 4, 2, 3, 4, 3), 4, 2)

# Preparing the inputs of the h-cbn method
bc <- Spock$new(
     poset = dag,
     numMutations = 4,
     genotypeMatrix = gMat
)
# Running the h-cbn model
resultsH <- hcbnSingle(bc)

## -----------------------------------------------------------------------------
mlLambdaH <- resultsH$lambda
logLikelihoodH <- resultsH$summary["Loglike"]

## -----------------------------------------------------------------------------
# The poset
dag <- matrix(c(3, 3, 4, 4, 1, 2, 1, 2), 4, 2)
# Preparing the inputs of the h-cbn method
bc <- Spock$new(
     poset = dag,
     numMutations = 4,
     genotypeMatrix = gMatMut
)
# Running the h-cbn model
resultsHMut <- hcbnSingle(bc)

## -----------------------------------------------------------------------------
examplePath <- getExamples()[1]
bc <- Spock$new(
     poset = readPoset(examplePath)$sets,
     numMutations = readPoset(examplePath)$mutations,
     genotypeMatrix = readPattern(examplePath)
)
resultsH2 <- hcbnSingle(bc)

## ----message=FALSE------------------------------------------------------------
posets <- readRDS(system.file("extdata", "Posets.rds", package = "CBN2Path"))

bc <- Spock$new(
     poset = posets,
     numMutations = 4,
     genotypeMatrix = gMat
)
resultsH3 <- hcbn(bc, nCores = 3)

## -----------------------------------------------------------------------------
logLikH <- sapply(resultsH3, \(x) x$summary["Loglike"])

## -----------------------------------------------------------------------------
indx <- which.max(logLikH)
mlPosetH <- posets[[indx]]

## ----fig.width=4.25, fig.height=4.25------------------------------------------
visualizeCBNModel(mlPosetH)

## ----eval=FALSE---------------------------------------------------------------
# lambdaC <- as.numeric(resultsC2$lambda)
# lambdaH <- as.numeric(resultsH2$lambda)
# dagC <- resultsC2$poset$sets
# dagH <- resultsH2$poset$sets
# 
# probC <- pathProbCBN(dagC, lambdaC, 10)
# probH <- pathProbCBN(dagH, lambdaH, 10)

## ----message=FALSE------------------------------------------------------------
probC1 <- pathProbQuartetCTCBN(gMat)
probC2 <- pathProbQuartetCTCBN(gMatMut)

probH1 <- pathProbQuartetHCBN(gMat)
probH2 <- pathProbQuartetHCBN(gMatMut)

## -----------------------------------------------------------------------------
visualizeProbabilities(probC1)
visualizeProbabilities(probC2)

## -----------------------------------------------------------------------------
visualizeProbabilities(probC2, geneNames = genes)

## -----------------------------------------------------------------------------
visualizeProbabilities(probH1)
visualizeProbabilities(probH2)

## -----------------------------------------------------------------------------
jsdC <- jensenShannonDivergence(probC1, probC2)
jsdH <- jensenShannonDivergence(probH1, probH2)
jsdC
jsdH

## -----------------------------------------------------------------------------
predC1 <- predictability(probC1, 4)
predC2 <- predictability(probC2, 4)
predC1
predC2
predC1 - predC2

## -----------------------------------------------------------------------------
predH1 <- predictability(probH1, 4)
predH2 <- predictability(probH2, 4)
predH1
predH2
predH1 - predH2

## -----------------------------------------------------------------------------
pathwayC1 <- pathwayCompatibilityQuartet(gMat)
pathwayC2 <- pathwayCompatibilityQuartet(gMatMut)

## -----------------------------------------------------------------------------
rhoC1 <- cor(pathwayC1, probC1, method = "spearman")
rhoC2 <- cor(pathwayC2, probC2, method = "spearman")
rhoC1
rhoC2
rhoH1 <- cor(pathwayC1, probH1, method = "spearman")
rhoH2 <- cor(pathwayC2, probH2, method = "spearman")
rhoH1
rhoH2

## ----message=FALSE------------------------------------------------------------
probR1 <- pathProbQuartetRCBN(gMat)
probR2 <- pathProbQuartetRCBN(gMatMut)

## -----------------------------------------------------------------------------
visualizeProbabilities(probR1)
visualizeProbabilities(probR2)

## -----------------------------------------------------------------------------
jsdR <- jensenShannonDivergence(probR1, probR2)
jsdR

## -----------------------------------------------------------------------------
predR1 <- predictability(probR1, 4)
predR2 <- predictability(probR2, 4)
predR1
predR2
predR1 - predR2

## -----------------------------------------------------------------------------
rhoR1 <- cor(pathwayC1, probR1, method = "spearman")
rhoR2 <- cor(pathwayC2, probR2, method = "spearman")
rhoR1
rhoR2

## ----warning=FALSE, results='hide',eval=FALSE---------------------------------
# probB1 <- pathProbQuartetBCBN(gMat)
# probB2 <- pathProbQuartetBCBN(gMatMut)

## ----eval=FALSE---------------------------------------------------------------
# visualizeProbabilities(probB1)
# visualizeProbabilities(probB2)

## ----eval=FALSE---------------------------------------------------------------
# jsdB <- jensenShannonDivergence(probB1, probB2)
# jsdB

## ----eval=FALSE---------------------------------------------------------------
# predB1 <- predictability(probB1, 4)
# predB2 <- predictability(probB2, 4)
# predB1
# predB2
# predB1 - predB2

## ----eval=FALSE---------------------------------------------------------------
# rhoB1 <- cor(pathwayC1, probB1, method = "spearman")
# rhoB2 <- cor(pathwayC2, probB2, method = "spearman")
# rhoB1
# rhoB2

## -----------------------------------------------------------------------------
fitnessVector <- c(0, 0.1, 0.2, 0.1, 0.2, 0.4, 0.3, 0.2, 0.2, 0.1, 0, 0.6, 0.4, 0.3, 0.2, 1)
g <- generateMatrixGenotypes(4)

## -----------------------------------------------------------------------------
visualizeFitnessLandscape(fitnessVector)

## -----------------------------------------------------------------------------
probW <- pathProbSSWM(fitnessVector,4)

## -----------------------------------------------------------------------------
visualizeProbabilities(probW)

## ----fig.width=4.25, fig.height=4.25------------------------------------------
predW <- predictability(probW, 4)

## -----------------------------------------------------------------------------
sessionInfo()