dbr.per1K and equivalent dbr gave different number of doublets removed

[MengjunWu]

Hi,

I am running scDblFinder using dbr.per1K (0.004) as well as manually calculated dbr (ncol(sre)/1000 * 0.004), they two should be equivalent from my understanding, however they resulted in different number of doublets removed. More doublets removed using dbr.per1K compared to dbr. I plotted histogram of scDblFinder.score, they look the same, just the thresholding is different -- dbr.per1K threshold is more strict than dbr. Any idea why it is? Thanks!

[plger]

Hi,

1. are you running this on a single sample? (if not, please do so for testing purposes)
2. how different are the thresholds?
3. are you running this with the same random seed?

[MengjunWu]

Thanks for the reply!

1. Yes, I am running on a single sample, with/without cluster mode
   2&3. Yes, I am running the same random seed, here are the code and results:

Use dbr.per1K without cluster:

set.seed(137)
sce.filt <- suppressWarnings(scDblFinder(sce, dbr.per1k =  0.004))
hist(sce.filt$scDblFinder.score)

The results:

Creating ~13064 artificial doublets...

Dimensional reduction

Evaluating kNN...

Training model...

iter=0, 3238 cells excluded from training.

iter=1, 3238 cells excluded from training.

iter=2, 3238 cells excluded from training.

Threshold found:0.363

2256 (13.8%) doublets called

Use dbr without cluster:

set.seed(137)
doublet_rate <- ncol(sre.full)/1000*0.004
doublet_rate <- ceiling(doublet_rate*100)/100
sce.filt2 <- suppressWarnings(scDblFinder(sce, dbr =  doublet_rate))
hist(sce.filt2$scDblFinder.score)

The results:

Creating ~13064 artificial doublets...

Dimensional reduction

Evaluating kNN...

Training model...

iter=0, 2805 cells excluded from training.

iter=1, 2805 cells excluded from training.

iter=2, 2805 cells excluded from training.

Threshold found:0.413

1872 (11.5%) doublets called

Use dbr.per1K with clusters:

set.seed(137)
sce.filt5 <- suppressWarnings(scDblFinder(sce, clusters=clusters, dbr.per1k =  0.004))
hist(sce.filt5$scDblFinder.score)

The results:

18 clusters

Creating ~13064 artificial doublets...

Dimensional reduction

Evaluating kNN...

Training model...

iter=0, 2551 cells excluded from training.

iter=1, 2551 cells excluded from training.

iter=2, 2551 cells excluded from training.

Threshold found:0.331

2794 (17.1%) doublets called

Use dbr with clusters:

set.seed(137)
sce.filt3 <- suppressWarnings(scDblFinder(sce, clusters=clusters, dbr =  doublet_rate))
hist(sce.filt3$scDblFinder.score)

The results:

18 clusters

Creating ~13064 artificial doublets...

Dimensional reduction

Evaluating kNN...

Training model...

iter=0, 2067 cells excluded from training.

iter=1, 2067 cells excluded from training.

iter=2, 2067 cells excluded from training.

Threshold found:0.424

2299 (14.1%) doublets called

The score distributions are pretty similar of the four cases:

The version of scDblFinder I am using is scDblFinder_1.20.2

Many thanks in advance!

[MengjunWu]

Another thing I noticed is dbr.per1K is to some extent insensitive to the value passed to it. I tried both 0.004 and default 0.008, the percentage of doublet removed are similar.

set.seed(137)
sce.filt4 <- suppressWarnings(scDblFinder(sce, clusters=clusters))
hist(sce.filt4$scDblFinder.score)

The results:

18 clusters

Creating ~13064 artificial doublets...

Dimensional reduction

Evaluating kNN...

Training model...

iter=0, 2551 cells excluded from training.

iter=1, 2551 cells excluded from training.

iter=2, 2551 cells excluded from training.

Threshold found:0.307

2912 (17.8%) doublets called

The distribution is the same as the case when using 0.004

[plger]

Not sure why you're doing ceiling(doublet_rate*100)/100, but given your large number of cells this isn't likely to be a problem here.

I quickly tried on dummy data, and setting it manually gave the same results:

> sce <- mockDoubletSCE(ncells=c(1000,2000,200,500,3000))
> set.seed(123); sce <- scDblFinder(sce)
Creating ~5754 artificial doublets...
Dimensional reduction
Evaluating kNN...
Training model...
iter=0, 961 cells excluded from training.
iter=1, 887 cells excluded from training.
iter=2, 848 cells excluded from training.
Threshold found:0.213
536 (7.5%) doublets called
> set.seed(123); sce <- scDblFinder(sce, dbr=ncol(sce)/1000*0.008)
Creating ~5754 artificial doublets...
Dimensional reduction
Evaluating kNN...
Training model...
iter=0, 961 cells excluded from training.
iter=1, 887 cells excluded from training.
iter=2, 848 cells excluded from training.
Threshold found:0.213
536 (7.5%) doublets called

I also had a quick look in a real data and same the same thing. Also changing changing the dbr.per1k has a reproducible effect (note however that the expected doublet rate is only part of what is used in the thresholding: the classification accuracy is typically more important).

So it looks like it's very dataset-specific. Your distributions are actually odd, I normally have much more of a bimodal distribution, i.e. with a peak close to one, and fewer values in-between. I can see why this makes thresholding much harder.

From your outputs one can see that the difference happens right at the beginning: the number of real cells excluded as doublets at iteration 0 is already different, and (unexpectedly) it doesn't update at all, which is weird. I will investigate why the difference at iteration 0 comes up, but I suspect you have a more fundamental problem. I noticed you're suppressing warnings, what is it that comes up? Do you have some very weird or low quality cells in there?

I would also suggest that you once run scDblFinder with dbr.sd=1, so that the threshold is not influenced by the expected rate, and see what comes up and how that compares to your expectation.

[MengjunWu]

Many thanks for the detailed reply! The idea of ceiling is to always take slightly larger doublet rate compared to the exactly doublet rate to be more conservative. (Also as I observed the dbr.per1k is slightly more strict so I was wondering if there is any rounding thing there)

I have already did QC to remove low quality cells. One thing specific to the data is that it is from a well staged developmental process with dramatic RNA content changes , sometimes there can be two-to-three times difference in the median RNA umi count between cell types in different stages (the highest has median about 18000 UMI, 6000 genes detected, while the lowest has median about 6000 UMI, 3000 gene detected). So would this affect the performance? I have also tried dbr.sd=1, I think it is a bit too harsh in removing things, so decided to go with the default setting.

The warning message is:

Warning message in check.deprecation(deprecated_cv_params, match.call(), ...):
“Passed invalid function arguments: max_depth, eval_metric, tree_method, subsample, nthread, eta. These should be passed as a list to argument 'params'. Conversion from argument to 'params' entry will be done automatically, but this behavior will become an error in a future version.”
Warning message in throw_err_or_depr_msg("Parameter(s) have been removed from this function: ", :
“Parameter(s) have been removed from this function: label. This warning will become an error in a future version.”

[plger]

Hi,
Sorry I didn't get back to you on your last question.
Just wanted to let you know that there was some trouble with being compatible with the different xgboost releases (see #132 ), which may have affected your results.