LIANA++ for Benchmarks and Developers
Daniel Dimitrov
2023-02-24
Source:vignettes/liana_devel.Rmd
liana_devel.RmdLIANA ++
To facilitate further method development and the benchmarks of current and future methods, in this vignette we provide further instructions relevant for method developers.
liana_pipe and Contributions
To facilitate method development, we provide a guide how to make use of the architecture already implemented in LIANA.
First, we start by converting the object to SCE (if Seurat) and basic
checks via liana_rep. Then we generate a LR_res table with
many relevant columns, such as stats for ligands and receptors per cell
type. Note that liana_pipe also takes complexes into account by
dissociating them into subunits (see
?liana:::decomplexify). These are later re-assembled back
into complexes using the recomplexify function.
sce <- liana_prep(testdata)
# Obtain LR summary res with CellPhoneDB and dissociate complexes
lr_res <- liana_pipe(sce,
op_resource = select_resource("CellPhoneDB")[[1]] %>%
liana:::decomplexify(),
base = exp(1) # here we assume log-transformation
# of the library-normalized counts
) %>%
glimpse()## Rows: 1,818
## Columns: 25
## $ ligand <chr> "DLL1", "CCL5", "CCL4", "TNFSF14", "CD40LG", "LTA", "…
## $ receptor <chr> "NOTCH1", "CCR5", "CCR5", "LTBR", "CD40", "TNFRSF1A",…
## $ ligand.pval <dbl> 1.0000000000, 0.0000114871, 0.0279585592, 1.000000000…
## $ ligand.FDR <dbl> 1.0000000000, 0.0003948689, 0.2078001025, 1.000000000…
## $ ligand.stat <dbl> 0.5000000, 0.1888889, 0.3922222, 0.4994444, 0.4994444…
## $ receptor.pval <dbl> 1.0000000000, 0.5440267057, 0.5440267057, 0.333710695…
## $ receptor.FDR <dbl> 1.000000000, 0.960392173, 0.960392173, 0.679306779, 0…
## $ receptor.stat <dbl> 0.4994444, 0.4822222, 0.4822222, 0.4833333, 0.7272222…
## $ source <chr> "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B"…
## $ target <chr> "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B"…
## $ ligand.expr <dbl> 0.00000000, 0.49167104, 0.08073624, 0.03217220, 0.039…
## $ receptor.expr <dbl> 0.04998350, 0.04227786, 0.04227786, 0.00000000, 0.685…
## $ ligand.scaled <dbl> -0.149312933, -0.945434584, -0.548615303, 0.017698966…
## $ receptor.prop <dbl> 0.03333333, 0.03333333, 0.03333333, 0.00000000, 0.500…
## $ ligand.prop <dbl> 0.00000000, 0.30000000, 0.06666667, 0.03333333, 0.033…
## $ receptor.scaled <dbl> -0.08950937, -0.02089337, -0.02089337, -0.14931046, 0…
## $ ligand.sum <dbl> 0.08632215, 6.59848126, 2.18621009, 0.08603449, 0.214…
## $ receptor.sum <dbl> 0.23999457, 0.14310802, 0.14310802, 0.09812520, 0.780…
## $ receptor.pem <dbl> 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.490…
## $ ligand.pem <dbl> 0.000000000, 0.000000000, 0.000000000, 0.006821435, 0…
## $ ligand.log2FC <dbl> -0.12768047, -4.01952187, -2.98806898, -0.02103795, -…
## $ receptor.log2FC <dbl> -0.071214715, -0.042705600, -0.042705600, -0.17490143…
## $ global_mean <dbl> 0.1990563, 0.1990563, 0.1990563, 0.1990563, 0.1990563…
## $ ligand.complex <chr> "DLL1", "CCL5", "CCL4", "TNFSF14", "CD40LG", "LTA", "…
## $ receptor.complex <chr> "NOTCH1", "CCR5", "CCR5", "LTBR", "CD40", "TNFRSF1A",…``
The LR result summary is then passed to the get_*
functions then provide ways to summarize this information into a score
(LRscore in this case) that can be used to the prioritize
interactions.
We also need to pass the function that we wish to use to account for heteromeric complexes
get_sca(lr_res, complex_policy="mean",
expr_prop=liana_defaults()[["expr_prop"]],
return_all=FALSE) %>% glimpse## Rows: 172
## Columns: 12
## $ source <chr> "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B"…
## $ target <chr> "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B"…
## $ ligand.complex <chr> "ALOX5", "MIF", "ANXA1", "CD72", "GRN", "SELL", "TNFS…
## $ ligand <chr> "ALOX5", "MIF", "ANXA1", "CD72", "GRN", "SELL", "TNFS…
## $ receptor.complex <chr> "ALOX5AP", "CD74", "FPR1", "SEMA4D", "TNFRSF1B", "SEL…
## $ receptor <chr> "ALOX5AP", "CD74", "FPR1", "SEMA4D", "TNFRSF1B", "SEL…
## $ receptor.prop <dbl> 0.3666667, 1.0000000, 0.1000000, 0.1666667, 0.1666667…
## $ ligand.prop <dbl> 0.3000000, 0.6000000, 0.2666667, 0.2333333, 0.2666667…
## $ ligand.expr <dbl> 0.3937931, 0.9104887, 0.3938854, 0.3774934, 0.3739728…
## $ receptor.expr <dbl> 0.5182538, 4.5311493, 0.1062514, 0.2230945, 0.2157197…
## $ global_mean <dbl> 0.1990563, 0.1990563, 0.1990563, 0.1990563, 0.1990563…
## $ LRscore <dbl> 0.6941424, 0.9107452, 0.5068361, 0.5931461, 0.5879491…get_* functions refer to liana_call and
liana_score functions behind the scenes, so please refer to
those for more info.
Install LIANA++
Install a LIANA++ environment with all original software as a way to minimize bias for benchmark studies.
LIANA++ uses CCC methods from both R and Python, as such dependencies for both need to be installed.
To install the LIANA++ framework run the following code in R:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!requireNamespace("devtools", quietly = TRUE))
install.packages("devtools")
library(devtools)
Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS = TRUE) # ignore warning from iTALK
BiocManager::install("ComplexHeatmap") # required for Connectome
devtools::install_github('saezlab/OmnipathR@ff3ad88e3915747e1b557bf44ac5396f9525dd7e') # install 4.0 version of OmnipathR
# install tools
devtools::install_github("sqjin/CellChat")
devtools::install_github('msraredon/Connectome', ref = 'master')
devtools::install_github("Coolgenome/iTALK", build_vignettes = FALSE)
# A modified version of SingleCellSignalR (SCA) that enables external resources
devtools::install_github(repo = "saezlab/SingleCellSignalR_v1",
subdir = "SingleCellSignalR")
# Finally, install LIANA
devtools::install_github('saezlab/liana')Squidpy and NATMI are written in Python, as such a python environment with the prerequisites of these methods needs to be set up.
If you also wish to run the CellPhoneDB re-implementation from Squidpy, please set up a conda environment by running the following lines in the terminal:
conda create -n liana_env
conda activate liana_env
conda install -c anaconda python=3.8.5
pip install squidpyPlease use the .yml file to set up a conda environment:
conda env create -f liana_env.ymlNote! NATMI and Squidpy are set by default to look for this conda environment (“liana_env”).
To use NATMI with LIANA, clone the modified NATMI repo into liana path by running the following in the terminal:
cd *insert fullpath*/liana
git clone https://github.com/saezlab/NATMINote that, NATMI was forked from asrhou/NATMI and changed to be agnostic of the working directory when loading the resources.
If needed, use the following to locate the liana package
system.file(package = "liana")
# "/home/user/R/x86_64.../4.0/liana"It might be easier to just also clone the liana repo.
call_* functions
Use the call_ functions to run the pipelines of the original methods, rather than their re-implementations. These functions work with Seurat objects and they are instead converted to SCE. These are largely getting deprecated and will be moved to a seperate repository in the future.
# RUN cellchat alone with OmniPath's Consensus
sca_res <- call_sca(
testdata,
op_resource = select_resource('Consensus')[[1]]
)
# Show CellChat Results
sca_resSimilarly one could call all call_* functions via
liana_wrap
liana_calls <- liana_wrap(testdata,
method = c('call_connectome',
'call_sca',
'call_italk',
'call_natmi'
),
resource = c("Default")) # Run with Default Resources
liana_calls %>% liana_aggregate() # Connectome finds no significant hits here
# Note that for italk we use mean(ligand.logFC, receptor.logFC) as a way to compare
# to the other methodsSession information
options(width = 120)
sessioninfo::session_info()