S3. Score Consistency

Here, we briefly demonstrate the similarities and discrepancies between LIANA’s output in python and R.

[1]:

library(liana, quietly = TRUE)
library(plyr, quietly = TRUE)
library(parallel, quietly = TRUE)

data.path <- file.path('..', '..', 'data')
output_folder <- file.path(data.path, 'liana-outputs/')

[2]:

# parallelization
n.cores<-parallel::detectCores() # default to all
if (n.cores<=1){
    par = F
}else{
    par = T
}

Comparison with Python output:

There are minor differences with the LIANA implementation in python that lead to outputs not being identical

SingleCellSignalR Magnitude (LRscore): precision - slightly different after 3rd decimal place
LogFC Specificity (logfc_comb): similar relative differences but different exact values
CellPhoneDB Specificity (pvalue): similar relative differences but different exact values
CellChat: not run by default in R

Let’s check the consistency in the magnitude aggregate rank score when running the different methods that report magnitude (excluding CellChat, which is not present by default in R).

[3]:

covid_data <- readRDS(file.path(data.path, 'covid_balf_norm.rds'))
md <- covid_data@colData
barcodes <- rownames(md[md$sample == 'C100', ])

sdata <- covid_data[, barcodes]

[4]:

liana_res_partial <- liana_wrap(sce = sdata,
                                idents_col='cell.type',
                                assay.type = 'logcounts',
                                verbose=T,
                                method = c("cellphonedb", "natmi", 'sca'),
                                parallelize = par, workers = n.cores)
liana_aggregate_partial <- liana_aggregate(liana_res = liana_res_partial,
                                           aggregate_how='magnitude', verbose = F)
liana_aggregate_partial<-liana_aggregate_partial[c('source', 'target', 'ligand.complex', 'receptor.complex', 'aggregate_rank')]

Running LIANA with `cell.type` as labels!

`Idents` were converted to factor

Cell identities with less than 5 cells: Plasma were removed!

Warning message in exec(output, ...):
“5591 genes and/or 0 cells were removed as they had no counts!”
Warning message:
“`invoke()` is deprecated as of rlang 0.4.0.
Please use `exec()` or `inject()` instead.
This warning is displayed once every 8 hours.”
LIANA: LR summary stats calculated!

Now Running: Cellphonedb

Warning message:
“`progress_estimated()` was deprecated in dplyr 1.0.0.
ℹ The deprecated feature was likely used in the liana package.
  Please report the issue at <https://github.com/saezlab/liana/issues>.”
Now Running: Natmi

Now Running: Sca

[5]:

liana_aggregate_partial <- plyr::arrange(liana_aggregate_partial, source, target, ligand.complex, receptor.complex)
write.csv(liana_aggregate_partial, paste0(output_folder, 'magnitude_ranks_R.csv'))
head(liana_aggregate_partial)

A tibble: 6 × 5
source	target	ligand.complex	receptor.complex	aggregate_rank
<chr>	<chr>	<chr>	<chr>	<dbl>
B	B	ACTR2	ADRB2	0.44463670
B	B	ADAM17	ITGB1	0.08916285
B	B	ADAM17	RHBDF2	0.36320725
B	B	ADAM28	ITGA4	0.43629682
B	B	APOC2	LRP1	1.00000000
B	B	APOE	LRP1	0.72715808

Note, to run the correlation, make sure to have run the companion Python tutorial up to the point where you save the csv named “magnitude_ranks_python.csv”.

[6]:

# read and format python aggregate rank
lap_python<-read.csv(paste0(output_folder, 'magnitude_ranks_python.csv'))
clp<-colnames(lap_python)
clp[4]<-'ligand.complex'
clp[5]<-'receptor.complex'
colnames(lap_python)<-clp
lap_python<-lap_python[colnames(lap_python) != 'X']

# merge the two scores
la<-dplyr::inner_join(liana_aggregate_partial, lap_python)

spearmanr<-cor.test(la$aggregate_rank,
         la$magnitude_rank, method = c('spearman'))$estimate[['rho']]
print(paste0('The spearman correlation bewteen R and python aggregate magnitude scores is: ',
             format(round(spearmanr, 2), nsmall = 2)))

Joining with `by = join_by(source, target, ligand.complex, receptor.complex)`
Warning message in cor.test.default(la$aggregate_rank, la$magnitude_rank, method = c("spearman")):
“Cannot compute exact p-value with ties”

[1] "The spearman correlation bewteen R and python aggregate magnitude scores is: 0.98"

[ ]: