MultiSEp: Predicting Gene Dependency Relationships (GDRs) from Multiomics Data

In MultiSEp: Predict Synthetic Lethality and Other Gene Dependency Relationships from Multiomics Data

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)

library(MultiSEp)
library(kableExtra)
library(dplyr)
library(knitr)
options(ignore.interactive = FALSE)

1. Introduction

The MultiSEp R package integrates 'omics data to predict Gene Dependency Relationships (GDRs) such as synthetic lethality or induced dependency. Functionality is also available to analyse and visualise the results.

MultiSEp (Multimodal Subsets of Expression) takes a matrix of continuous data in a gene by sample format (normally gene expression but possibly other gene-linked 'omics data such as DNA methylation), determines the optimal number of clusters for each gene, and produces cluster assignments for each sample. These clusters are integrated with different 'omics data types, such as gene effect score data from RNAi or CRISPR screens, or mutation calls from genome sequencing data to generate candidate GDRs. The MultiSEp package builds upon the web server SynLeGG (https://www.overton-lab.uk/synlegg/, Wappett \textit{et al.} 2021), with significant additional functionality including new analytics for drug discovery applications with patient data.

This vignette demonstrates how MultiSEp may be applied to predict, analyse and visualise GDRs. We will follow the flow of the package structure outlined in Figure 1, below, which summarises key functions. Wrapper functions are shown with dashed lines in Figure 1 and will be considered in the 'Quickstart' section, below.

Some of the example data used in this package was provided by the Cancer Dependency Map 2024Q2 (https://doi.org/10.25452/figshare.plus.25880521.v1, Arafeh \textit{et al.} 2025).

Citing MultiSEp

If you use MultiSEp in your work, please cite: AS McKie, M Wappett, H Vandierendonck, IM Overton (2026). The MultiSEp R package. https://cran.r-project.org/package=MultiSEp Please note: a preprint article will replace the above citation in the near future.

\clearpage

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 = "/"
)

\includepdf[pages=-, fitpaper=true, pagecommand={\thispagestyle{fancy}\fancyfoot[C]{\thepage}}]{r fig1_path} \clearpage

2. Quickstart

Wrapper functions are available for key MultiSEp workflows, summarised below. Broadly, two analysis strategies are available for Gene Dependency Relationship (GDR) discovery. Both depend upon Gaussian mixture modelling (GMM). One approach is suitable for analysis of various samples (tissue, organoids or cell lines \textit{etc.}) and deploys the binomial test to investigate GDRs in multiple data types (left-hand side, Figure 1). For example, by analysis of the GMM-derived gene expression clusters for two genes; analysing gene expression clusters with a categorical data type such as mutations or methylation status; and comparing categorical data types against themselves or each other (for example mutation vs mutation). A further approach was primarily developed for analysis of \textit{in vitro} data such as from cell lines or organoids (SynLeGG methods) and a key application evaluates differences in CRISPR scores for one gene across GMM-derived expression clusters for a second gene; with statistical significance estimated using the t-test (right-hand side, Figure 1).

Predicting Synthetic Lethality (SL) and other gene dependencies with mts_omics()

mts_omics() can be applied to various 'omics data types to predict pairwise SL and other Gene Dependency Relationships (GDRs), including when only gene expression data is available (\textit{e.g.} RNA-seq). Therefore, mts_omics() may be used for analysis of data from many different kinds of samples, such as tumour biopsies. Gaussian mixture modelling (GMM) defines clusters of samples where gene function is diminished, or lost (LoF); for example corresponding to the lowest expression values for a gene. Depletion of samples in the LoF clusters combined for two genes is characteristic of SL for most data types. An exception is that enrichment in the LoF clusters is expected for SL with CRISPR data, for example correlating CRISPR and gene expression (XPR). Additionally, mts_omics() can be used to discover arbitrary GDR patterns by setting \textit{SyntheticLethalityPrediction=FALSE}. Examples of mts_omics() are given below for analysis of XPR data, as well as XPR with CRISPR data. When the same data-type is used, GDRs are symmetrical. However, GDRs arising from two data types have two orientations. For example, the analysis of gene A's XPR against gene B's CRISPR is a different orientation to analysing gene B's XPR against gene A's CRISPR. Both orientations are evaluated by default, if the data are available.

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
)

Analysis with mts_omics() using categorical data, such as mutations and (continuous) gene expression, requires binary values for the categorical data type:

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)

SL prediction with only categorical data is also possible, for example with mutations:

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
  )

Predicting SL with CRISPR and gene expression (SynLeGG)

mts_GeneDepCrispr() runs SL discovery with Gaussian mixture modelling (GMM) and t-test; expected inputs are a gene expression matrix and a CRISPR (or RNAi) gene effect score matrix in gene by sample format. An earlier version of this function generated the 'CRISPR' results for the SynLeGG resource (https://www.overton-lab.uk/synlegg/, Wappett \textit{et al.} 2021). Data should be continuous and in log2 format. We recommend running in Linux with multiple cores for timely production of results. The standard analysis with mts_GeneDepCrispr() predicts SL pairs. An example is given below:

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)

Induced dependency prediction with CRISPR and gene expression (SynLeGG)

Predicted induced dependency pairs (synthetic dosage lethal relationships) are output from mts_GeneDepID(), which runs the induced dependency t-test pipeline using a gene expression matrix and a CRISPR gene effect score matrix in gene by sample format. An example is given below:

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)

Predicting SL with mutation and gene expression (SynLeGG)

mts_GeneDepMutation() runs the chi-squared test for SL discovery using a gene expression matrix and a mutation summary matrix in gene by sample format. A previous version of this function generated results in the 'Mutation' section of the SynLeGG resource (https://www.overton-lab.uk/synlegg/, Wappett \textit{et al.} 2021). The output from mts_GeneDepMutation() is predicted SL gene pairs. Please see the example below:

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

3. Prediction of Gene Dependency Relationships (GDRs) with Functional Genomics Data

The availability of large-scale 'omics data enables analysis of the consequences of gene loss of function (LoF) across many different genetic, transcriptional and disease backgrounds. These context-specific differences in the impact of gene LoF upon cell viability can reveal Gene Dependency Relationships (GDRs), including Synthetic Lethality (SL). Data from CRISPR and RNAi screens coupled with gene expression and other 'omics data across many different cell lines can sample the co-occurrence of gene pair LoF and so may be leveraged to discover depletion or enrichment patterns that are characteristic of SL or other GDRs. Tumours have endogenous perturbations, for example arising from mutations, which influence gene functional status. The variability in endogenous perturbations across different tumours may also be analysed to discover SL vulnerabilities by investigating mutually exclusive LoF patterns for pairs of genes.

Data partitioning to map gene function states

In order to discover mutually exclusive LoF patterns, MultiSEp requires a map of gene function states. This is relatively straightforward for categorical data, for example high impact mutations are expected to disrupt gene function. One of many possible mechanisms is where a mutation creates a premature STOP codon that leads to nonsense-mediated decay or a truncated protein, resulting in LoF due to the absence of the protein (or part of the protein) in the cell. On the other hand, more sophisticated methods are required to partition continuous data such as gene expression. The MultiSEp functionality for assigning gene function states using categorical and continuous data is outlined in the sections below.

Unsupervised clustering to reveal gene states for continuous data

MultiSEp runs Gaussian mixture modelling (GMM) with Expectation-Maximisation to discover clusters (sometimes called modes), that represent different gene function states. Cardinality (#clusters) is determined by Bayesian Information Criterion regularisation (Lubbock \textit{et al.} 2013). A typical use case derives clusters of samples from gene expression (XPR) data.

The following example demonstrates GMM for gene expression data with mts_mixModelCluster_XPR() to derive cluster assignments for genes based on their expression patterns across the samples. The output of mts_mixModelCluster_XPR() is a list of data frames which correspond to a gene (row name) from the input gene expression object. The 'cores' argument may be increased to enable multithreading, resulting in faster execution time (not available in MS Windows at present). In addition to GMM, mts_mixModelCluster_XPR() applies filters that ensure the genes taken forwards for analysis have sufficient data to produce meaningful results. In particular, these filters ensure that the gene is expressed above background levels in a sufficient number of samples:

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

The ExpressionClusters object (above) is a list of dataframes where each gene corresponds to a separate dataframe containing the samples, sample values (expression) and cluster assignment:

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

A data-agnostic function is available, mts_mixModelCluster(); which performs GMM without the above filters. Quality metrics vary between data types and should be considered before running this function. A value of at least 20 samples with values above the background value is recommended. The input is a data matrix and therefore many different data types may be analysed; such as drug sensitivity, transcriptomics, proteomics, CRISPR, methylation etc. The output of mts_mixModelCluster() is a list of dataframes, as shown above.

Threshold-based partitioning of continuous data

Data may also be partitioned using a fixed threshold value. The mts_crisprPartition() function is designed to separate scores (e.g. from Chronos) into two classes that approximate to (1) 'cell death' or (2) 'little or no cell death'. This function could also be applied to partition other data types into classes representing different gene functional states or responses, assuming that a reasoned threshold value can be chosen. This approach is motivated by the observation that unsupervised clustering may sometimes define cluster boundaries that do not align well with functional states of interest - because cluster assignment is driven by the distribution of data points. For example, CRISPR disruption of some genes may kill most cell lines (or leave most cell lines unaffected). In these cases, GMM might define cluster boundaries that do not distinguish between samples that undergo elevated cell death and those that are relatively unaffected by the gene LoF. Accordingly, setting a cluster boundary based upon prior knowledge can be beneficial. The t-test (SynLeGG) methods such as mts_GeneDepCrispr() do not require boundaries to be set in the CRISPR data, which could have advantages especially when there are small sample sizes. The output of mts_crisprPartition() is a list of dataframes, as shown in the section above.

Partitioning categorical data

Categorical data, for example mutational status, does not typically require clustering. The mts_formatMatrix() function may be used for preprocessing if the input is a binary matrix. For example a matrix representing mutation or methylation data may encode gene functional status with numerical values. A value of 1 within the input binary matrix is assigned by the mts_formatMatrix() to 'HIGH' impact (loss of function) and a value of 2 is 'LOW' impact. The output follows a similar format to that produced by mts_mixModelCluster_XPR() where each row name (e.g. gene) corresponds to a separate dataframe that contains the impact groups for the samples evaluated. The following example predicts pairwise genetic dependencies using a mutational matrix formatted by mts_formatMatrix():

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

Analysing canonical dependency patterns for gene pairs

The mts_patternDetection() function can identify many different GDR patterns by evaluation of sample depletion (or enrichment) in pairwise combinations of gene clusters - which form cells in a contingency table. The mts_omics() function is a wrapper for mts_patternDetection() and therefore may provide convenient access to the same functionality. Statistical significance is estimated with the binomial test. Depletion of samples in the bottom-left table cell is evaluated by default; where both genes are assumed to have loss of, or low, function. However, all of the contingency table cells may be evaluated for depletion or enrichment patterns. The example below predicts synthetic lethality (SL) for gene expression data by evaluating the canonical mutually exclusive LoF pattern:

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)

The output from mts_patternDetection() characterises the table cell evaluated ('Cluster_Combination'); including the actual and expected count of samples, as well as the statistical significance of the difference in the counts (q-value). The total number of cells in the contingency table is also given in the 'TotalCluster_Combinations' column:

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

The pattern of samples for a gene pair can be visualised with mts_plotClusterDistribution(). The resulting 'MultiSEp plot' shows boundaries between the expression clusters; the clusters for gene1 and gene2 are respectively shown on the x-axis and y-axis. Figure 2 shows depletion of samples in the bottom-left contingency table cell for PSMB8 and TTC7B, consistent with an SL relationship:

\addtocounter{figure}{1}

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

The mts_patternDetection() function may also be applied to evaluate all of the cluster combinations in the contingency table (i.e. across all table cells), by altering the SyntheticLethalityPrediction flag:

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

Results are shown on the next page:

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

\newpage For example, SNAP25 and TUBA1A have depletion of samples in the bottom right table cell (4x1 Cluster_Combination), shown in Figure 3:

```r. 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")

### Mapping GDRs with sample enrichment: CRISPR data

We expect synthetic lethal (SL) relationships to involve enrichment of
samples in the bottom left contingency table cell if we are analysing CRISPR data,
as long as the second data type (e.g. gene expression) codes loss of function (LoF)
with low values. As before, the patterns formed by the samples may also be explored
across all cluster combinations in the contingency table by setting the
SyntheticLethalityPrediction argument to FALSE. The code on the next page analyses CRISPR and
gene expression data:

```r
 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)

Results are shown below:

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

For example, there is a significant enrichment of samples in the bottom-left contingency table cell for DNAJC19 CRISPR scores and DNAJC15 gene expression; this pattern is consistent with a synthetic lethal relationship (Figure 4).

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

A table of GDRs for gene expression with CRISPR effect scores is shown below:

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

\newpage For example ATP1B1 and ATP1B3 have significant enrichment in the 1_x_1 table cell (death occurs in cells with ATP1B3 loss and low ATP1B1 gene expression, consistent with an SL relationship) and in the 2_x_2 table cell (cells with high ATP1B1 gene expression are viable when ATP1B3 is lost), shown in Figure 5:

```r; 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")

\newpage
## Tissue-specific GDR detection
Context-specific gene dependencies are enriched within particular tissues or lineages (\textit{i.e.} contexts). These relationships may reveal selective vulnerabilities with a potentially wider therapeutic window, due to sparing healthy tissue outside of the biological context where SL occurs. The mts_patternDetection() function can be used to detect these context-specific gene dependency patterns.

### Predicting tissue-specific GDRs with transcriptome data
Transcriptomic data can be analysed to predict context-specific GDRs, including synthetic lethal (SL) relationships. The example below deploys mts_patternDetection() to evaluate candidate lung cancer SL relationships.
```r
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))

When MYC and CCDC88A have low expression, there is a depletion of lung cancer cell lines (red) in the bottom left cell of the MultiSEp plot, consistent with an SL relationship and visualised in Figure 6 with mts_plotClusterDistribution():

mts_plotClusterDistribution(mixModelClusters = mixturemodelClusters,
                            TScluster1 = LungXPR_MCC,
                            gene1 = "MYC", gene2 = "CCDC88A")

Predicting tissue-specific GDRs with transcriptome and CRISPR data

This toy example shows how mts_patternDetection() can investigate SL relationships in a tissue-specific context with CRISPR and gene expression (XPR) data; specifically for lung cancer cell lines. Please note that the effect size threshold should not normally be specified and the q-value may be set to 0.05 for standard usage.

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

We can visualise the candidate lung cancer SL relationship between ACSL1 and PSMB5 with the mts_plotClusterDistribution() function (Figure 7). The lung cancer samples (red) are largely located in the bottom left cell of the contingency table, which is consistent with a candidate SL relationship.

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

Predicting SL with CRISPR and gene expression (SynLeGG)

In the quickstart section we outlined mts_GeneDepCrispr() which is the simplest way to invoke the SynLeGG methods for predicting SL. Here we consider individual steps in the pipeline in more detail. mts_clusterAvg() takes as input a log2 gene expression matrix, and a gene effect score matrix in gene by sample format, performs cluster assignment and calculates cluster arithmetic means for each CRISPR gene.

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

The object clusterAssign is a list where each item corresponds to a gene result. Each gene result is a list of two data frames. Data frame 1 is a table of mean gene effect scores for every CRISPR gene:

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

Data frame 2 returns the cluster assignment and log2 gene expression value for each sample:

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

The clusterAssign object is input for the mts_Crispr() function which performs statistical evaluation of dependency relationships by calculating a log2 fold-change and two-tailed t-test p-value between the CRISPR scores for each consecutive pair of mRNA expression clusters. The results are a ranked table of candidate synthetic lethal gene pairs, ordered by qvalue.

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

The ttestCrisprResults output table contains lots of useful information including ranked gene pair names and the number of gene expression clusters (the CRISPR data is not clustered for this method).

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

\newpage The number of samples per gene expression cluster (nMode) and the mean CRISPR score (mMode) per gene expression cluster are also available:

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

Finally, ttestCrisprResults contains some important statistics for each gene pair. Considering consecutive clusters, the gene effect score log2 fold change values are given (in the 'shift' columns) - this is the difference in average gene effect score between neighbouring clusters; p-values and q-values are also calculated for these changes. Negative fold-changes indicate a lower gene effect score in the lower expression cluster.

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 shows the CRISPR gene effect score for the gene NMT1, split according to the NMT2 gene expression clusters. NMT1 CRISPR scores are significantly more negative in the low NMT2 expression cluster, predicting an SL relationship.

data("depMapTissue_subset")
mts_plotCRISPRGeneCluster(mrna_gene = "NMT2", crispr_gene = "NMT1",
                             crisprMatrix = depMapCRISPRscores_subset[c(4,27),],
                             tissueMatrix = depMapTissue_subset, 
                             resultList = clusterAssign, plotType = "Integrated")

\newpage It is also possible to view the gene expression distribution in the MultiSEp clusters (Figure 9):

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")

\newpage Or a waterfall plot of the CRISPR data and gene expression clusters (Figure 10):

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")

\newpage

Tissue-specific GDR detection (SynLeGG)

GDRs are often restricted to a single tissue type, or a subset of tissues. Many genes are expressed in a subset of tissues, which is one factor that leads to tissue-specific dependencies and is sometimes referred to as tissue penetrance. MultiSEp includes functionality to evaluate tissue-specific GDRs:

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

As with the pan-tissue results, the output (i.e. ttestCrisprResultsTS) contains lots of useful information including ranked gene pair names, tissue and cluster number:

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

The number of samples per gene expression cluster and the average CRISPR score per gene expression cluster are shown below. This table could be quite sparse due to low tissue-specific sample number, or non-represented modes for some tissues.

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

The results also include some important statistics for each gene pair. A set of tissue-specific consecutive cluster gene effect shift values (log2 fold change) are returned, which indicate the difference in average gene effect score between ascending, neighbouring clusters. The p-values and q-values for these differences are also available:

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

For the functional paralogues NMT1 and NMT2, the mutually exclusive loss relationship is most statistically significant in Ovarian Cancer (Figure 11):

mts_plotCRISPRGeneCluster(mrna_gene = "NMT2", crispr_gene = "NMT1",
                          crisprMatrix = depMapCRISPRscores_subset[c(4,27),],
                          tissueMatrix = depMapTissue_subset, 
                          diseaseFilter = "Ovarian Cancer",
                          resultList = clusterAssign, plotType = "Integrated")

Induced dependency prediction (t-test)

Prediction of Synthetic Lethality (SL) is canonically based upon a pattern of mutually exclusive loss. However it is possible that increased expression may be a signature of dependency, sometimes called synthetic dosage lethality. The mts_InducedDependency() function defaults to return GDRs where cell death is enhanced for loss of gene 1 and high expression of gene 2. Here, the clusterAssign object is input into the mts_InducedDependency() function which performs statistical evaluation of dependency relationships by calculating a log2 fold-change and two-tailed t-test p-value between the CRISPR scores for each consecutive pair of mRNA expression clusters. For induced dependency, a positive fold-change threshold is needed (default):

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

The mts_InducedDependency() output is structured in the same way as for mts_Crispr(). Positive shift values (fold-changes) indicate lower gene 1 effect scores (implying cell death/growth inhibition) in the higher gene 2 expression cluster:

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

For example, CHMP4B CRISPR gene effect scores are lower for high STX2 expression, predicting an induced dependency (synthetic dosage lethal) relationship (Figure 12):

data("depMapTissue_subset")
mts_plotCRISPRGeneCluster(mrna_gene = "STX2", crispr_gene = "CHMP4B", 
                             crisprMatrix = depMapCRISPRscores_subset, 
                             tissueMatrix = depMapTissue_subset, 
                             resultList = clusterAssign, plotType = "Integrated")

We can again run the mts_CrisprTS() pipeline to assess tissue-specific penetrance for induced dependency:

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)

In the case of the functional paralogues STX2 and CHMP4B, the induced dependency relationship is most statistically significant in Breast Cancer (Figure 13):

mts_plotCRISPRGeneCluster(mrna_gene = "STX2", crispr_gene = "CHMP4B", 
                          crisprMatrix = depMapCRISPRscores_subset, 
                          tissueMatrix = depMapTissue_subset, 
                          diseaseFilter = "Breast Cancer",
                          resultList = clusterAssign, plotType = "Integrated")

Mutation and expression GDRs (SynLeGG)

The mts_Mutation() function takes the clusters from mts_mixModelCluster(), a gene by sample mutation matrix containing 'WT' or a mutation call, and optionally a tissue to sample tissue map. mts_Mutation() detects dependencies between mutations and gene expression clusters by performing the chi-squared test, which produces p-values. Analysis may be performed across all tissues, or in a tissue-specific manner.

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

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

The results are a data frame containing data for the dependency of mutations with one or more gene expression modes (clusters), quantified by chi-squared p-values, and the number of mutations in each mode. Both classical SL and induced dependency GDR patterns are evaluated. The mutation_per_mode attribute gives details of the number of mutations identified in the different gene expression modes (clusters):

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

The top hit in the example above is an enrichment of TP53 mutation in the low expression cluster of CDKN1A (Figure 14). However, wild-type TP53 is expected to activate and lead to higher expression of CDKN1A. Therefore this example is not a true induced dependency (synthetic dosage lethal) relationship and likely represents a false positive.

mts_plotMutation(resultList = mixModelClusters, mrna_gene = "CDKN1A", mut_gene = "TP53", 
                 mutMatrix = depMapMUT_subset, tissueMatrix = depMapTissue_subset)

\newpage Addition of a tissue matrix allows mts_Mutation() to be run in a tissue-specific manner.

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

Tissue-specific results are now available:

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

The enrichment of TP53 mutations in the CDKN1A low expression cluster is particularly clear in Ovarian cancer (Figure 15):

mts_plotMutation(resultList = mixModelClusters, mrna_gene = "CDKN1A", 
                 mut_gene = "TP53", mutMatrix = depMapMUT_subset, 
                 tissueMatrix = depMapTissue_subset, 
                 diseaseFilter = "Ovarian Cancer")

Analysis of mutational neighbourhoods in SL networks for therapeutic target prioritisation

The mts_targetCoverageFromMutations() function can integrate mutational data with a network of predicted SL relationships in order to estimate the proportion of a population that could benefit from inhibiting any one of the candidate target genes. For example, the population could be cancer patients. We reason that an inactivating mutation in a synthetic lethal partner gene exposes (cancer) cells to killing when the functional gene in the pair is inhibited by a drug. mts_targetCoverageFromMutations() takes as input an object with gene pairs (two columns, corresponding to each gene in the pair) and a gene by sample mutation matrix containing 'WT' or a mutation call (e.g. details of amino acid changes produced by point mutations or INDELs). Gene pairs of interest may be provided through the \textit{edgeData} argument; alternatively, results generated with mts_Mutation() may be assigned to \textit{edgeData}. Increasing the value of \textit{cores} is recommended if possible:

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

The mts_targetCoverageFromMutations() results are a data frame containing: the candidate target gene ('Gene'), neighbours ('NeighbourGenes'), the samples that have a mutation in a neighbouring gene ('MutatedSamples'), and the population frequency of mutations in any of the genes neighbouring the candidate target gene ('NeighbourMutationFrequency').

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

Ordering by the neighbour mutation frequency reveals the candidate target genes with highest predicted population coverage. In this toy example, AP1M1 has 44 neighbours and 1721 samples have mutations in one of the neighbouring genes. Please note that this function is intended for use with patient data to estimate population coverage, whereas the example here is with cell line data; additionally the data were not filtered to ensure that the mutation produces loss of function (e.g. selecting 'high impact' mutations) - which is an important step.

4. Running MultiSEp in a High-Performance Computing (HPC) Environment for Evaluation of Context-specific SL at Genome Scale

This section presents an approach for predicting synthetic lethal (SL) relationships at genome scale with high-performance computing (HPC). We investigate candidate SL relationships for Multiple Myeloma gene expression from the CoMMpass Study (Skerget \textit{et al.} 2024). This 'all-vs-all' comparison involves running the SL discovery workflow for approximately 105 million gene pairs, therefore use of a HPC cluster is required for timely production of results.

Step 1: Produce gene expression clusters

The first step in predicting the synthetic lethal network is to obtain gene expression clusters with Gaussian mixture modelling (GMM) using the mts_mixModelCluster_XPR() function to analyse the MMRF CoMMpass transcriptomic dataset (n=718 samples). When run in a HPC environment, where an individual cluster node may have a large number of cores, you may increase the value of the cores argument (for example, set to 50 below). We recommend saving the output of mts_mixModelCluster_XPR() as an .RData file for downstream analysis, for example:

CoMMpassXPRClusters <- mts_mixModelCluster_XPR(
  dataMatrix = CoMMpass_transcriptomicData,
  GeneXPRthresh = 10, # data is not logged
  NumSampleThresh = 20, 
  cores = 50)

save(CoMMpassXPRClusters, file = "CoMMpassXPRClusters.RData")

Step 2: Preparing the clustered genes for parallel processing

The second step runs the mts_genepairsChunkGeneration() function to prepare the data for parallel processing and involves splitting the CoMMpassXPRClusters object into multiple .RData objects, each containing a subset of the genes to be analysed. mts_genepairsChunkGeneration() also removes any genes where expression clustering failed, for example if there is only one expression value for the gene. The num_tasks argument stores the number of 'chunks' (tasks) to be run in the array job on the HPC cluster and determines the number of files to be saved to the current working directory. For this exemplar, involving a genome scale multiple myeloma GDR network, mts_mixModelCluster_XPR() found clusters for 14,549 genes. Accordingly, all vs all GDR discovery requires binomial tests for 105,829,426 gene pairs; when num_tasks is set to 1000, there will be 105,830 gene pairs in 999 'chunks' (tasks) and 105,256 pairs analysed in the 1000th 'chunk' (task):

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

Step 3: Predicting SL with mts_omics()

To calculate 105,829,426 binomial test depletion p-values (undirected pairs for the 14,549 genes) a script is required to run mts_omics() in the HPC environment. The code below is run as a Rscript within a SLURM array job (with 1000 tasks), and where the SLURM task ID is given as a command-line argument. Please note that the cluster node cores requested via SLURM should match the cores argument - set to 10 in the code below; we recommend 4 Gb RAM per CPU for this application (specified in the bash script below).

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)

In this example, we run a bash script to call the above Rscript:

```{bash, eval=FALSE}

!/bin/bash

SBATCH --job-name=mts_omics_array

SBATCH --mail-type=END

SBATCH --mail-user=[e-mail]

SBATCH --ntasks=1

SBATCH --cpus-per-task=10

SBATCH --mem-per-cpu=4G

SBATCH --partition=optional_queue_name_on_your_cluster

SBATCH --chdir=/path/to/your/directory/for/logfiles/log/

SBATCH --array=1-1000

other commands that you might need to set up the environment, for example "module load R"

Rscript /path/to/directory/where/you/have/the/above/Rscript/rscriptName.R $SLURM_ARRAY_TASK_ID

When run as a SLURM array job with 1000 tasks there are 1000 files which
may be combined as follows:

```{bash, eval=FALSE}
awk 'FNR == 1 && NR != 1 { next } { print }' SL_CoMMpass_{1..1000}.txt > combined_SL_CoMMpass.txt

The analysis for predicting synthetic lethal depletion patterns was performed in parallel for each subset of gene pairs, therefore adjusting p-values for multiple comparisons and q-value filtering should be calculated on the combined set of gene pairs:

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)

In the example above, the predicted multiple myeloma SL network is saved to the file multiSEp_SL_XPR_qval05.txt

Step 4: Identifying deleterious mutations to inform drug discovery

A gene that has an SL relationship with one or more genes that are mutated may be an attractive drug target for precision oncology. In this case study, mutational data was sourced from the CoMMpass Study (Skerget \textit{et al.} 2024). Our analysis only utilised deleterious mutations predicted to be high impact by the variant effect prediction tools, SnpEff and SnpSift. The bash command below can be run from the command line to extract high impact mutations:

```{bash, eval=FALSE} (head -n 1 Somatic_Observed_SNV_INDEL_mutations.txt; grep HIGH Somatic_Observed_SNV_INDEL_mutations.txt) > HighImpact_Somatic_Observed_SNV_INDEL_mutations.txt

The following python script may be used to process the data into a
format suitable for analysis by MultiSEp - creating a mutational score
matrix in gene by patient format with values either "WT" or
"MUT":

```{python, eval=FALSE, python.reticulate = FALSE}
"""Build a matrix of genes (rows) by samples (columns) from a list of
high-impact variant calls. Each cell is 'MUT' if the gene carries >=1 variant
in that sample, otherwise 'WT'.
"""
import pandas as pd
import os
# directory of the script
script_dir = os.path.dirname(os.path.realpath(__file__))
data_path = "HighImpact_Somatic_Observed_SNV_INDEL_mutations.txt"
data = pd.read_csv(data_path, sep='\t')

matrix = (pd.crosstab(data['GENEID'], data['sample'])
            .gt(0)
            .replace({True: 'MUT', False: 'WT'}))
# Name the row and column headers used in the saved file.
matrix.index.name = 'Gene'
matrix.columns.name = 'Sample'

print(matrix)
print(matrix.isnull().values.any())   # expect False: every cell is filled
print(matrix.isnull().sum().sum())     # expect 0

# save the matrix to a tab-separated file in the output folder
output_filename = "MutationalMatrix_CoMMpassSLnet.txt"
output_path = os.path.join(script_dir, "..", "output", output_filename)
matrix.to_csv(output_path, sep='\t', index=True)

Step 5: Predicting the population coverage achieved by inhibition of a candidate target gene

The mutational score matrix and predicted SL network are analysed by the mts_targetCoverageFromMutations() function in order to predict the proportion of patients that may benefit from targeting a candidate gene.

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)

Visualisation of predicted SL network for Multiple Myeloma

The network resulting from the above analysis may be visualised with applications such as Cytoscape (Shannon \textit{et al.} 2003). Figure 16 visualises MultiSEp results for the Myeloma transcriptome data (MMRF CoMMpass, Skerget \textit{et al.} 2024):

r<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")

\newpage

5. References

Arafeh, R., Shibue, T., Dempster, J. M., Hahn, C. W., & Vazquez, F. (2025). The present and future of the Cancer Dependency Map. \textit{Nature Reviews Cancer} 25, 59–73. https://doi.org/10.1038/s41568-024-00763-x

Lubbock, A. L. R., Katz, E., Harrison, D. J., & Overton, I. M. (2013). TMA Navigator: network inference, patient stratification and survival analysis with tissue microarray data. \textit{Nucleic Acids Research} 41(W1), W562–W568. https://doi.org/10.1093/nar/gkt529

Skerget, S., Penaherrera, D., Chari, A., Jagannath, S., Siegel, D. S., Vij, R., ... & Keats, J. J. (2024). Comprehensive molecular profiling of multiple myeloma identifies refined copy number and expression subtypes. \textit{Nature Genetics} 56, 1878–1889. https://doi.org/10.1038/s41588-024-01853-0

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., Amin, N., Schwikowski, B., & Ideker, T. (2003). Cytoscape: A software environment for integrated models of biomolecular interaction networks. \textit{Genome Research} 13(11), 2498–2504. https://doi.org/10.1101/gr.1239303

Wappett, M., Harris, A., Lubbock, A. L. R., Lobb, I., McDade, S., & Overton, I. M. (2021). SynLeGG: analysis and visualization of multiomics data for discovery of cancer ‘Achilles Heels’ and gene function relationships. \textit{Nucleic Acids Research} 49(W1), W613–W618. https://doi.org/10.1093/nar/gkab338

MultiSEp documentation built on Aug. 27, 2026, 5:07 p.m.