## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  out.width = "100%",
  fig.align = "center",
  message = FALSE,
  dev = "png", dpi = 150
)
options(pillar.width = 85)
options(pillar.max_dec_width = 3)
options(pillar.sigfig = 2)

## ----setup, include=FALSE-----------------------------------------------------
library(MultiSEp)
library(kableExtra)
library(dplyr)
library(knitr)
options(ignore.interactive = FALSE)

## ----include=FALSE------------------------------------------------------------
input_file <- knitr::current_input(dir = TRUE)
input_dir <- if (is.null(input_file)) getwd() else dirname(input_file)
fig1_path <- normalizePath(
  file.path(input_dir, "figures", "Figure1.pdf"),
  winslash = "/"
)

## ----eval=FALSE---------------------------------------------------------------
# data("depMapXPR_subset")
# SL_XPRvsXPR <- mts_omics(
#   dataMatrix=depMapXPR_subset[1:5,],
#   qVal=1, # 0.05 is recommended (default)
#   cores=1) # increase if possible
# 
# SL_XPRvsCRISPR = mts_omics(
#   dataMatrix = depMapCRISPRscores_subset[1:5,],
#   dataMatrix2 = depMapXPR_subset[1:5,],
#   directionality = "enrichment",
#   effectsize = FALSE, # TRUE is recommended for CRISPR/XPR
#   qVal = 1, # 0.05 is recommended (default)
#   cores = 1 # increase if possible
# )

## ----eval=FALSE---------------------------------------------------------------
# data("depMapMUT_small")
# data("depMapXPR_small")
# 
# mutImpact <- as.data.frame( # convert to binary values
#     ifelse(as.matrix(depMapMUT_small) == "WT", 2, 1),
#     stringsAsFactors = FALSE, check.names = FALSE)
#   rownames(mutImpact) <- rownames(depMapMUT_small)
# 
# SL_mut_vs_XPR <- mts_omics(
#     dataMatrix   = mutImpact,
#     dataMatrix2  = depMapXPR_small,
#     categorical1 = TRUE, # mutation data is categorical
#     directionality = "depletion",
#     qVal         = 0.05,
#     effectsize   = FALSE, # recommended for categorical data
#     cores        = 1)
# 

## ----eval=FALSE---------------------------------------------------------------
# mutClusters <- mts_formatMatrix(matrix = mutImpact, cores = 1)
# 
# SL_mut_only = mts_omics(
#   mixModelClusters1 = mutClusters,
#   directionality    = "depletion",
#   qVal              = 1,
#   effectsize        = FALSE, # recommended for categorical data
#   cores             = 1
#   )
# 

## ----eval=FALSE---------------------------------------------------------------
# data("depMapXPR_subset")
# data("depMapCRISPRscores_subset")
# mtsGeneDepResults <- mts_GeneDepCrispr(exprsMatrix=depMapXPR_subset[1:5,],
#                                        crisprMatrix=depMapCRISPRscores_subset,
#                                        cores=1, fcVal=-0.1, pVal=0.1)

## ----eval=FALSE---------------------------------------------------------------
# data("depMapXPR_subset")
# data("depMapCRISPRscores_subset")
# mtsGeneDepIDResults <- mts_GeneDepID(exprsMatrix=depMapXPR_subset[1:5,],
#                                      crisprMatrix=depMapCRISPRscores_subset,
#                                      cores=1, fcVal=0.1, pVal=0.1)

## ----eval=FALSE---------------------------------------------------------------
# data("depMapMUT_subset")
# data("depMapTissue_subset")
# mtsGeneDepMutResults <- mts_GeneDepMutation(exprsMatrix=depMapXPR_subset[1:5,],
#                         tissueMatrix = depMapTissue_subset,
#                         mutMatrix = depMapMUT_subset, pVal=0.1)

## ----results='hide'-----------------------------------------------------------
data("depMapXPR_subset")
ExpressionClusters <- mts_mixModelCluster_XPR(dataMatrix = depMapXPR_subset[1:5,],
                          GeneXPRthresh = 3.321928, # equal to log2(10) (default)
                          NumSampleThresh = 20)  # (default)

## -----------------------------------------------------------------------------
knitr::kable(head(ExpressionClusters[[1]], 10))
names(ExpressionClusters)

## ----eval=FALSE---------------------------------------------------------------
# data("MUT_impactMatrix")
# Formatted_MUTimpactMatrix <- mts_formatMatrix(matrix = MUT_impactMatrix, cores = 1)
# mutation_vs_XPR_SL <- mts_patternDetection(
#                                 mixModelClusters1 = Formatted_MUTimpactMatrix,
#                                 mixModelClusters2 = ExpressionClusters,
#                                 effectsize = TRUE,
#                                 effectsize_threshold = 0, # not recommended, just for this example
#                                 directionality = "depletion",
#                                 qVal = 1) # 0.05 is recommended and the default

## ----results='hide'-----------------------------------------------------------
data(mixModelClusters_vignette) # precomputed for speed
mixturemodelClusters = mixModelClusters_vignette
XPRvsXPR_SL <- mts_patternDetection(mixModelClusters1=mixturemodelClusters,
                                                   directionality = "depletion",
                                                   p_adjustMethod = "BY",
                                                   qVal = 0.01, effectsize = TRUE,
                                                   SyntheticLethalityPrediction = TRUE)

## -----------------------------------------------------------------------------
XPRvsXPR_SL %>%
  arrange(q_value) %>%
  arrange(Actual_Count) %>%
  kable(row.names = FALSE, format = "latex", booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped", "scale_down", "hold_position")) %>%
  add_header_above(c("Results" = 12))

## ----Figure-2, fig.width=10.5, fig.height=7, fig.cap="MultiSEp plot of PSMB8 and TTC7B expression. Clusters from Gaussian mixture modelling of the gene expression data are shown (top for PSMB8, left for TTC7B) and the cluster boundaries form a contingency table. Each circle in the table represents a sample, coloured by the pairwise cluster combination. The bottom-left table cell corresponds to low expression of both genes, and may be assumed to represent low or loss of function. The absence of samples in the bottom-left cell is consistent with an SL relationship between TTC7B and PSMB8."----

mts_plotClusterDistribution(mixModelClusters = mixturemodelClusters, 
                            gene1 = "PSMB8", gene2 = "TTC7B")

## ----results='hide'-----------------------------------------------------------
GDRdepletionPattern <- mts_patternDetection(mixModelClusters1=mixturemodelClusters,
                                     directionality = "depletion",
                                     qVal=0.01,
                                     SyntheticLethalityPrediction = FALSE)

## -----------------------------------------------------------------------------
GDRdepletionPattern %>%
  arrange(Actual_Count / Expected_Count) %>%
  kable(row.names = FALSE, format = "latex", booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped", "scale_down", "hold_position")) %>%
  add_header_above(c("Results" = 12))

## ----Figure-3, fig.dim = c(10.5, 7), fig.cap="MultiSEp plot of SNAP25 and TUBA1A expression. Depletion of samples in the bottom-right cell was identified by \\textit{mts\\_patternDetection()}. Clusters from Gaussian mixture modelling of the gene expression data are shown (top for SNAP25, left for TUBA1A) and the cluster boundaries form a contingency table. Each circle in the table represents a sample, coloured by the pairwise cluster combination."----

mts_plotClusterDistribution(mixModelClusters = mixturemodelClusters, 
                            gene1 = "SNAP25", gene2 = "TUBA1A")

## ----results='hide'-----------------------------------------------------------
 CRISPR_clusters = mts_crisprPartition(dataMatrix = depMapCRISPRscores_subset[c(3,9,17),])
CRISPR_XPR_enrichment_SL <- mts_patternDetection(
                                        mixModelClusters1=mixturemodelClusters,
                                        mixModelClusters2 = CRISPR_clusters,
                                        p_adjustMethod = "BH",
                                        qVal = 0.01, effectsize = FALSE, #  effectsize = TRUE is generally recommended, is set to FALSE just for this example
                                        directionality = "enrichment",
                                        SyntheticLethalityPrediction = TRUE)

CRISPR_XPR_enrichment_GDRs <- mts_patternDetection(
                                      mixModelClusters1=mixturemodelClusters,
                                      mixModelClusters2 = CRISPR_clusters,
                                      p_adjustMethod = "BH",
                                      qVal = 0.01, effectsize = FALSE,  #  effectsize = TRUE is generally recommended, is set to FALSE just for this example
                                      directionality = "enrichment",
                                      SyntheticLethalityPrediction = FALSE)

## -----------------------------------------------------------------------------
CRISPR_XPR_enrichment_SL %>%
  arrange(q_value) %>%
  kable(row.names = FALSE, format = "latex", booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped", "scale_down", "hold_position")) %>%
  add_header_above(c("Results" = 11))

## ----Figure-4, fig.dim = c(10.5, 7), fig.cap="MultiSEp plot of DNAJC19 CRISPR scores and DNAJC15 gene expression. Clusters from Gaussian mixture modelling of the gene expression data are shown at the top for DNAJC15, and the results of partitioning the DNAJC19 CRISPR scores at a value of -0.5 are shown in the density on the left-hand side of the Figure. The cluster boundaries form a contingency table and each circle in the table represents a sample, coloured by the pairwise cluster combination. Enrichment of samples in the bottom-left contingency table cell is consistent with an SL relationship."----
mts_plotClusterDistribution(mixModelClusters = mixturemodelClusters,
                            mixModelClusters2 = CRISPR_clusters,
                            gene1 = "DNAJC15", gene2 = "DNAJC19",
                            gene1_datatype="log2 gene expression",
                            gene2_datatype = "CRISPR gene effect score")

## -----------------------------------------------------------------------------
CRISPR_XPR_enrichment_GDRs %>%
  arrange(q_value) %>%
  kable(row.names = FALSE, format = "latex", booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped", "scale_down", "hold_position")) %>%
  add_header_above(c("Results" = 11))

## ----Figure-5, fig.dim = c(10.5, 7), fig.cap="MultiSEp plot of ATP1B3 CRISPR scores and ATP1B1 gene expression. Clusters from Gaussian mixture modelling of the gene expression data are shown at the top for ATP1B1, and the results of partitioning the ATP1B3 CRISPR scores at a value of -0.5 are shown by the distributions at the left-hand side. Each circle in the plot represents a sample, coloured by the pairwise cluster combination. This example highlights how enrichment of samples may be evaluated across the whole contingency table with \\textit{mts\\_patternDetection()}; the bottom-left (1,1) and top-right (2,2) table cells have statistically significant enrichment."----

mts_plotClusterDistribution(mixModelClusters = mixturemodelClusters,
                            mixModelClusters2 = CRISPR_clusters,
                            gene1 = "ATP1B1", gene2 = "ATP1B3",
                            gene1_datatype="log2 gene expression",
                            gene2_datatype = "CRISPR gene effect score")

## ----results='hide'-----------------------------------------------------------
cellLines <- depMapTissue_subset$cell_line[
  depMapTissue_subset$tissue == "Lung Cancer"]

LungXPR_MCC <- lapply(mixturemodelClusters, function(GMM) {
  GMM[GMM$Sample %in% cellLines, ]
})

Lung_SL_depletionPattern <- mts_patternDetection(
  mixModelClusters1 = LungXPR_MCC,
  directionality = "depletion",
  qVal = 1,
  SyntheticLethalityPrediction = TRUE)

## -----------------------------------------------------------------------------
Lung_SL_depletionPattern %>% 
  arrange(Actual_Count / Expected_Count) %>%
  head(10) %>% 
  kable(row.names = FALSE, format = "latex", booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped", "scale_down", "hold_position")) %>%
  add_header_above(c("Results" = 12))

## ----Figure-6, fig.dim = c(10.5, 7), fig.cap="MultiSEp plot for CCDC88A and MYC gene expression, highlighting lung cancer cell lines. Clusters from Gaussian mixture modelling of the gene expression data are shown at the top and left-hand-side for MYC, CCDC88A respectively. The cluster boundaries form a contingency table and each circle in the table represents a cell line. Lung cancer cell lines are coloured red and the remaining cell lines, from various tissues, are white. The lung cell lines are depleted from the bottom-left contingency table cell, consistent with a tissue-specific SL relationship between CCDC88A and MYC."----
mts_plotClusterDistribution(mixModelClusters = mixturemodelClusters,
                            TScluster1 = LungXPR_MCC,
                            gene1 = "MYC", gene2 = "CCDC88A")

## ----results='hide'-----------------------------------------------------------
LungCRISPR_MCC <- lapply(CRISPR_clusters, function(GMM) {
  GMM[GMM$Sample %in% cellLines, ]
})

Lung_SL_enrichmentPattern <- mts_patternDetection(
  mixModelClusters1 = LungXPR_MCC,
  mixModelClusters2 = LungCRISPR_MCC, 
  include_reverse_pairs = TRUE,
  SyntheticLethalityPrediction = TRUE,
  directionality = "enrichment",
  qVal = 1) # the default value is 0.05


## -----------------------------------------------------------------------------
Lung_SL_enrichmentPattern[Lung_SL_enrichmentPattern$Gene1=='ACSL1',] %>%
  kable(row.names = FALSE, format = "latex", booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped", "scale_down", "hold_position")) %>%
  add_header_above(c("Results" = 12))

## ----Figure-7, fig.dim = c(10.5, 7), fig.cap="MultiSEp plot for PSMB5 CRISPR scores and ACSL1 gene expression, highlighting lung cancer cell lines. Clusters from Gaussian mixture modelling of ACSL1 gene expression are shown at the top and the results of partitioning the PSMB5 CRISPR scores at a value of -0.5 are shown by the distributions at the left-hand side. The cluster boundaries form a contingency table and each circle in the table represents a cell line. Lung cancer cell lines are coloured red and the remaining cell lines, from various tissues, are white. The lung cell lines are largely located in the bottom-left contingency table cell, consistent with a tissue-specific SL relationship between PSMB5 and ACSL1."----

mts_plotClusterDistribution(mixModelClusters = mixturemodelClusters,
                            mixModelClusters2 = CRISPR_clusters,
                            TScluster1 = LungXPR_MCC,
                            TScluster2 = LungCRISPR_MCC,
                            gene1 = "ACSL1", gene2 = "PSMB5")

## ----results='hide'-----------------------------------------------------------
clusterAssign <- mts_clusterAvg(exprsMatrix = depMapXPR_subset[1:5,], 
                            depMapCRISPRscores_subset[c(4,27),])

## -----------------------------------------------------------------------------
knitr::kable(clusterAssign[[1]][[1]], format = "latex", longtable = FALSE)%>%
  kable_styling(latex_options = c("striped"))

## -----------------------------------------------------------------------------
knitr::kable(head(clusterAssign[[1]][[2]], 10))%>%
  kable_styling(latex_options = c("striped"))

## ----results='hide'-----------------------------------------------------------
ttestCrisprResults <- mts_Crispr(resultList = clusterAssign, 
                               exprsMatrix = depMapXPR_subset[1:5,],
                               depMapCRISPRscores_subset[c(4,27),])

## -----------------------------------------------------------------------------
knitr::kable(head(ttestCrisprResults[,1:3]), row.names = FALSE) %>%
  kable_styling(latex_options = c("striped", "hold_position")) %>%
  add_header_above(c("Results: Gene_Summary"= 3))

## -----------------------------------------------------------------------------
knitr::kable(head(ttestCrisprResults[,4:13]), row.names = FALSE, format = "latex", booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped", "scale_down", "hold_position")) %>%
  add_header_above(c("Results: Mode_Information"= 10))

## -----------------------------------------------------------------------------
knitr::kable(head(ttestCrisprResults[,14:25]), row.names = FALSE, format = "latex", booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped", "scale_down", "hold_position")) %>%
  add_header_above(c("Results: Statistics"= 12))

## ----Figure-8, fig.dim = c(10.5, 7), fig.cap="NMT1 CRISPR Scores are more negative with low NMT2 gene expression. The y axis shows NMT1 CRISPR scores and each box corresponds to a NMT2 gene expression cluster. The cell lines (circles) are coloured according to tissue type, please see the key for details."----
data("depMapTissue_subset")
mts_plotCRISPRGeneCluster(mrna_gene = "NMT2", crispr_gene = "NMT1",
                             crisprMatrix = depMapCRISPRscores_subset[c(4,27),],
                             tissueMatrix = depMapTissue_subset, 
                             resultList = clusterAssign, plotType = "Integrated")

## ----Figure-9, fig.dim = c(10.5, 7), fig.cap="Visualising MultiSEp gene expression clusters. The x-axis corresponds to NMT2 log2 mRNA expression and each curve represents a gene expression cluster, the allocation of colourings to cluster identity (mode) is shown in the key."----
data("depMapTissue_subset")
mts_plotCRISPRGeneCluster(mrna_gene = "NMT2", crispr_gene = "NMT1",
                             crisprMatrix = depMapCRISPRscores_subset[c(4,27),],
                             tissueMatrix = depMapTissue_subset, 
                             resultList = clusterAssign, plotType = "mrna_only")

## ----Figure-10, fig.dim = c(10.5, 7), fig.cap="Waterfall plot of NMT1 CRISPR score and NMT2 gene expression cluster. The bars are ordered from low to high CRISPR score and coloured according to the MultiSEp gene expression cluster, allocation of colourings to cluster identity is given in the key. Mode (cluster) 1 is enriched for low CRISPR scores, indicating a greater effect on cell viability/growth with loss of NMT1 when NMT2 has low expression."----
data("depMapTissue_subset")
mts_plotCRISPRGeneCluster(mrna_gene = "NMT2", crispr_gene = "NMT1",
                             crisprMatrix = depMapCRISPRscores_subset[c(4,27),],
                             tissueMatrix = depMapTissue_subset, 
                             resultList = clusterAssign, plotType = "crispr_only")

## ----results='hide'-----------------------------------------------------------
ttestCrisprResultsTS <- mts_CrisprTS(resultList = clusterAssign, 
                                     crisprMatrix = depMapCRISPRscores_subset[c(4,27),],
                                     fcVal = -0.1, pVal = 0.25, 
                                     tissueMatrix = depMapTissue_subset, 
                                     allDisRes = ttestCrisprResults, cores=1)

## -----------------------------------------------------------------------------
knitr::kable(head(ttestCrisprResultsTS[,1:4]), row.names = FALSE) %>%
  kable_styling(latex_options = c("striped", "hold_position")) %>%
  add_header_above(c("Results: Gene_Summary"= 4))

## -----------------------------------------------------------------------------
knitr::kable(head(ttestCrisprResultsTS[,5:14]), row.names = FALSE, format = "latex", booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped", "scale_down", "hold_position")) %>%
  add_header_above(c("Results: Mode_Information"= 10))

## -----------------------------------------------------------------------------
knitr::kable(head(ttestCrisprResultsTS[,15:26]), row.names = FALSE, format = "latex", booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped", "scale_down", "hold_position")) %>%
  add_header_above(c("Results: Statistics"= 12))

## ----Figure-11, fig.dim = c(10.5, 7), fig.cap="NMT1 CRISPR scores for NMT2 Gene expression clusters in Ovarian Cancer. The y axis shows NMT1 CRISPR scores and each box corresponds to a NMT2 gene expression cluster. Only Ovarian cancer cell lines are shown (circles)."----
mts_plotCRISPRGeneCluster(mrna_gene = "NMT2", crispr_gene = "NMT1",
                          crisprMatrix = depMapCRISPRscores_subset[c(4,27),],
                          tissueMatrix = depMapTissue_subset, 
                          diseaseFilter = "Ovarian Cancer",
                          resultList = clusterAssign, plotType = "Integrated")

## ----results='hide'-----------------------------------------------------------
ttestIDResults <- mts_InducedDependency(resultList = clusterAssign, 
                               exprsMatrix = depMapXPR_subset[1:5,],
                               crisprMatrix = depMapCRISPRscores_subset[c(4,27),], fcVal = 0.1)

## -----------------------------------------------------------------------------
knitr::kable(head(ttestIDResults[,14:25]), row.names = FALSE, format = "latex", booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped", "scale_down", "hold_position")) %>%
  add_header_above(c("Statistics"= 12))

## ----Figure-12, fig.dim = c(10.5, 7), fig.cap="Candidate induced dependency relationship between CHMP4B and STX2. CHMP4B CRISPR Scores are more negative with high STX2 gene expression, consistent with an induced dependency relationship. The y axis shows CHMP4B CRISPR scores and each box corresponds to a STX2 gene expression cluster. The cell lines (circles) are coloured according to tissue type, please see the key for details."----
data("depMapTissue_subset")
mts_plotCRISPRGeneCluster(mrna_gene = "STX2", crispr_gene = "CHMP4B", 
                             crisprMatrix = depMapCRISPRscores_subset, 
                             tissueMatrix = depMapTissue_subset, 
                             resultList = clusterAssign, plotType = "Integrated")

## ----eval=FALSE---------------------------------------------------------------
# ttestIDResultsTS <- mts_CrisprTS(resultList = clusterAssign,
#                                  crisprMatrix = depMapCRISPRscores_subset[c(4,27),],
#                                  fcVal = 0.1, # set to 0.1 to predict induced dependency relationships
#                                  pVal = 0.25,
#                                  tissueMatrix = depMapTissue_subset,
#                                  allDisRes = ttestIDResults, cores=1)

## ----Figure-13, fig.dim = c(10.5, 7), fig.cap="Candidate induced dependency relationship between CHMP4B and STX2 in breast cancer. CHMP4B CRISPR Scores are more negative with high STX2 gene expression, consistent with an induced dependency relationship. The y axis shows CHMP4B CRISPR scores and each box corresponds to a STX2 gene expression cluster. Only breast cancer cell lines are shown, where this predicted induced dependency relationship has the strongest p-value."----
mts_plotCRISPRGeneCluster(mrna_gene = "STX2", crispr_gene = "CHMP4B", 
                          crisprMatrix = depMapCRISPRscores_subset, 
                          tissueMatrix = depMapTissue_subset, 
                          diseaseFilter = "Breast Cancer",
                          resultList = clusterAssign, plotType = "Integrated")

## ----results='hide'-----------------------------------------------------------

mixModelClusters <- mts_mixModelCluster(dataMatrix = depMapXPR_subset[c(1,12,15,20,25),])

multisepMutAll <- mts_Mutation(resultList = mixModelClusters, mutMatrix = depMapMUT_subset,
                               pVal = 0.01)

## -----------------------------------------------------------------------------
knitr::kable(head(multisepMutAll), row.names = FALSE) %>%
  kable_styling(latex_options = c("striped", "hold_position")) %>%
  add_header_above(c("Mutation Cluster Enrichment"= 5))

## ----Figure-14, fig.dim = c(10.5, 7), fig.cap="TP53 mutations associate with lower CDKN1A expression. This example shows how mutations may be visualised across gene expression clusters. The y axis shows expression of CDKN1A and each box corresponds to the CDKN1A expression cluster. The samples (cell lines) are coloured according to the TP53 mutation class, which is detailed in the key."----
mts_plotMutation(resultList = mixModelClusters, mrna_gene = "CDKN1A", mut_gene = "TP53", 
                 mutMatrix = depMapMUT_subset, tissueMatrix = depMapTissue_subset)

## ----results='hide'-----------------------------------------------------------
multisepMutTS <- mts_Mutation(resultList = mixModelClusters, 
                              mutMatrix = depMapMUT_subset, pVal = 0.01, 
                              tissueMatrix = depMapTissue_subset)

## -----------------------------------------------------------------------------
knitr::kable(head(multisepMutTS, 5), row.names = FALSE) %>%
  kable_styling(latex_options = c("striped", "hold_position")) %>%
  add_header_above(c("Mutation Cluster Enrichment"= 5))

## ----Figure-15, fig.dim = c(10.5, 7), fig.cap="Tissue-specific analysis of TP53 mutations and CDKN1A expression. The enrichment of TP53 mutations in the low CDKN1A expression cluster is shown for Ovarian cancer cell lines. The y axis shows expression of CDKN1A and each box corresponds to the CDKN1A expression cluster. Samples are coloured according to the TP53 mutation class, which is detailed in the key."----
mts_plotMutation(resultList = mixModelClusters, mrna_gene = "CDKN1A", 
                 mut_gene = "TP53", mutMatrix = depMapMUT_subset, 
                 tissueMatrix = depMapTissue_subset, 
                 diseaseFilter = "Ovarian Cancer")

## -----------------------------------------------------------------------------
data("depMapMUT_subset")
netN <- mts_targetCoverageFromMutations(mutData=depMapMUT_subset, 
                                  edgeData=multisepMutAll, cores=1)

## -----------------------------------------------------------------------------
knitr::kable(head(netN[1:4]), row.names=FALSE) %>%
  kable_styling(latex_options = c("striped", "hold_position"))


## ----eval=FALSE---------------------------------------------------------------
# CoMMpassXPRClusters <- mts_mixModelCluster_XPR(
#   dataMatrix = CoMMpass_transcriptomicData,
#   GeneXPRthresh = 10, # data is not logged
#   NumSampleThresh = 20,
#   cores = 50)
# 
# save(CoMMpassXPRClusters, file = "CoMMpassXPRClusters.RData")

## ----eval=FALSE---------------------------------------------------------------
# result <- mts_genepairsChunkGeneration(mixModelClusters1 = CoMMpassXPRClusters,
#                                       num_tasks = 1000,
#                                       output_dir = "/example/outputdirectory/array_chunks/")

## ----eval=FALSE---------------------------------------------------------------
# library(MultiSEp)
# # Identify the SLURM task ID
# args <- commandArgs(trailingOnly=TRUE)
# SlurmTaskID  <- args[1]
# chunkdirectory  <- "/example/outputdirectory/array_chunks/"
# filename <- paste0(chunkdirectory,"genepairs_chunk_", SlurmTaskID, ".RData")
# 
# # Run binomial test
# load(filename)
# load("/example/outputdirectory/CoMMpassXPRClusters.RData")
# SL_CoMMpass <- mts_omics(
#     genepairs = subset_genepairs,
#     mixModelClusters1 = CoMMpassXPRClusters, # the GMM clustering object
#     SyntheticLethalityPrediction = TRUE,
#     p_adjustMethod = "BY",
#     qVal = 2,
#     effectsize = TRUE,
#     effectsize_threshold = 0,
#     directionality = "depletion",
#     cores = 10 # this number matches the cpus-per-task in the bash script below
# )
# # Save the results to file
# output_directory2 <- "/example/outputdirectory/"
# text_file <- paste0(output_directory2, "SL_CoMMpass_", SlurmTaskID, ".txt")
# write.table(SL_CoMMpass, file = text_file, quote = FALSE, sep = "\t",
#             row.names = FALSE)

## ----eval=FALSE---------------------------------------------------------------
# networkBinomial = read.table("combined_SL_CoMMpass.txt",
#                              header = TRUE)
# # Adjust p-values for multiple comparisons with the Benjamini & Yekutieli (BY) method
# p_adjusted = p.adjust(networkBinomial$p_value, method = "BY")
# networkBinomial$New_q_value = p_adjusted
# 
# # Filtering by recomputed q-value < 0.05 and effect size >= 0.9621389
# networkBinomialFiltered = subset(networkBinomial,
#                                  New_q_value < 0.05 & Effect_Size >= 0.9621389)
# 
# # Save to a file
# write.table(networkBinomialFiltered,
#             "multiSEp_SL_XPR_qval05.txt",
#             quote = FALSE, row.names=FALSE)

## ----eval=FALSE---------------------------------------------------------------
# SLedges <- read.table("multiSEp_SL_XPR_qval05.txt", header=TRUE)
# MutationMatrix <- read.table("MutationalMatrix_CoMMpassSLnet.txt",
#                             sep ="\t", header=TRUE, row.names=1)
# 
# PopulationCoverage <- mts_targetCoverageFromMutations(
#   mutData = MutationMatrix, edgeData = SLedges, cores = 10)
# 

## ----Figure-16, fig.dim = c(10.5, 7), fig.cap="Predicted SL network for Multiple Myeloma. Edges (connections) represent predicted SL relationships from MultiSEp analysis of the MMRF CoMMpass transcriptome data. The network (4,199 genes, 36,595 edges, \textit{q}<0.05) contains 0.035% of the 105,829,426 gene pairs evaluated by MultiSEp. Node (gene) colouring and size indicates the proportion of patients that are predicted to respond to therapeutic inhibition of the candidate target gene. Larger genes and warmer colours correspond to a higher predicted response rate. The smallest genes, with lowest predicted population coverage are shown at the same width and colour as the edges."----
knitr::include_graphics("figures/CoMMpassSLnetwork.png")

