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)
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"
[ ]: