p-value combination

[Captain-Pam]

Hi @plger

Thank you for creating such a great tool!

I have a question: in your paper and documentation, you used combined p-values from two methods (scDblFinder and Amulet). I'm planning to do the same, but I'm not sure how the final threshold is defined.

Should I filter by combined p-value < 0.05, or use FDR < 0.05?

Thanks a lot for your help!

[plger]

Hi,

I'm afraid I don't have a great answer beyond educated guess.

The aim of FDR correction is to control for FDR (disregarding type II errors), but with doublet identification the aim should arguably be more towards not missing too many doublets, or at the very least to strike a balance between the two. I would therefore not do FDR correction, or if you want simply use the adjusted Amulet p-value.

Regarding the exact threshold, I'm afraid the procedure implemented in the scDblFinder function cannot be applied for this case. In most datasets, if you do a histogram of the scores you'll see that most values are either very close to 0, or very close to 1. You want to ensure that your threshold is somewhere between the two modes, and as long as it is,, where exactly you set it will typically change the call for only some of the cells.

This being said, I understand that you do have to set it somewhere. I've seen three decent options:

1. set it to 0.5, the rational being that if a cell is 'more likely' to be a doublet than a real cell, you want to exclude it (although this isn't actually a correct interpretation of the amulet score)
2. use the expected proportion of doublets to set the threshold
3. use posterior >0.5 on a mixture model of the (logit) distributions at the two modes. That sounds fancier than it is, and in practice will result to an equivalent of setting the threshold at well over 0.5.
   I think there are arguments for an against each approach...

Pierre-Luc

[Captain-Pam]

Hi @plger
Thank you for your reply.
I agree that using FDR for doublet detection may not be the most appropriate approach.

I’m currently working with multiome data and am unsure what strategy is most suitable for removing doublets. Below is the approach I’m currently using — I would greatly appreciate it if you could take a look and let me know whether it seems reasonable:

1. RNA modality:
   I use both scDblFinder and DoubletFinder to obtain scDblFinder.score and pANN, respectively. Then I integrate the scores using Fisher’s method:
   sample$combn.fisher <- 1 - apply(cbind(sample$scDblFinder.score, sample$pANN), 1, function(x) {
   suppressWarnings(aggregation::fisher(1 - x))
   })
   Cells with combn.fisher > 0.5 are removed as predicted doublets.

2. ATAC modality:
   I use scDblFinder and scDblFinder::Amulet to obtain scDblFinder.p and p.value, then integrate the p-values as follows:
   res$combined <- apply(res[, c("scDblFinder.p", "p.value")], 1, function(x) {
   x[x < 0.001] <- 0.001 # to avoid over-skewing from very small values
   suppressWarnings(aggregation::fisher(x))
   })
   Cells with combined < 0.01 are considered doublets and removed.

Finally, I take the intersection of the remaining RNA and ATAC cells for downstream analysis.
Does this approach sound reasonable to you? I'm open to any suggestions you might have for improving it.

Pam

[plger]

This isn't unreasonable, if quite strict.

With multimodal data, I think there's a value in combining scores from the two modalities before calling doublets, but how exactly to best do it I'm not sure.

In particular, since AMULET calls also homotypic doublets (and we typically don't care much about those), I personally don't want to rely only on it alone, so I sometimes end up doing something like this:
(b>0.3 & a<0.01) | (b>0.6), where a is the ATAC/amulet-like score, and b is the RNA-based scDblFinder score.

But in your case I suppose you could also simply further combine the two scores. I really can't say what's best.

[Captain-Pam]

Thank you for your answer.

How to choose a better strategy is really difficult, I have been struggling for a long time. Maybe using your threshold gives me a better choice (b>0.3 & a<0.01) | (b>0.6).

1. After the last time, I noticed that it may be too strict. Also I read your document, I also intend to run scDblFinder twice on the RNA data and calculate the average score or votes.

The results will provide the category of prediction (scDblFinder.class), but also the score. As you mentioned earlier, the strategy of using a score of 0.3 or 0.6 as the threshold confused me whether to use the score or directly provide the predicted category. Perhaps there isn't much difference between them?

2. For ATAC part, I used the combination of scDblFinder.p and p.value (scDblFinder::Amulet) to calculate the combined combined.p value, with a threshold of combined.p<0.01. We typically don't care much about homotypic doublets, this seems cumbersome to use. Maybe just using scDblFinder::Amulet < 0.01 is enough and do not need scDblFinder in ATAC?

3. So in the end, we're using the intersection cells of two parts for integrating like WNN or MultiVI? Then we estimate scores in each cluster before deciding whether to remove some high-scoring clusters?

Many Thanks!
Pam

[plger]

Hi,

I agree. The main drawback of using amulet on its own is that it doesn't work well with cells that have a low library size, so that's why we use a combination. However, in multiome you're getting information from the RNA side, so I don't think it's necessary in such a context to use two methods on the ATAC modality. On such data I personally use just a combination of ATAC-amulet + RNA-scDblFinder.

And yes, I think it's wise afterwards to compute median scores per cluster to screen for potential high-scoring clusters.

Pierre-Luc

[Captain-Pam]

Thank you for your patient answer.

Finally, still have a few final confusing questions, I hope you can help me.

1. As you mentioned earlier, the strategy of using a score of 0.3 or 0.6 as the threshold, but the results also provide the predicted category (scDblFinder.class). Why don't you just delete the cells that are both double cells twice. Then, after calculating the score of each cluster (resolution = 1), decide to remove these subclusters.

2. If applying (b>0.3 & a<0.01) | (b>0.6) to remove doublets, that means there are no high-score cells in subclusters, because all cells are b < 0.6. Or would you delete subclusters with highest scores?

3. It is my pipeline:

1. I use cellbender to remove ambintRNA
2. simple QC to filter possible empty droplets (min.cells = 5, min.feature = 200, nCounts > 200)
3. scDblfinder and scDblfinder::Amulet. scDblFinder > 0.6 | (scDblFinder > 0.3 & Amulet.p < 0.01)
4. QC again like nCounts < +- 5MAD, nFeature < +- 5MAD, mt < 10

I really don't know how to deal with doublets, it's just my result (the figure below, scVI cluster with resolution = 1). Should I calculate the score for each cluster next, like removing a cluster with highest scores.

Sorry, I have a lot of questions, but I have been struggling for a long time. Any suggestions are worth it.

Million thanks.

[plger]

Hi,

1. Please clarify what you mean by "delete the cells that are both double cells twice."
2. It's common (especially if the artificial doublets are random) for some singlet cells to have a high-ish scDblFinder.score (e.g. because it happens to be very similar to homotypic artificial doublets), but the vast majority of cells should have a very low score. For this reason, when we compute the median scDblFinder.score per cluster, one expects that most clusters have a very low median score (e.g. <0.05). However it often happens that you have small clusters that have something much higher, even though it's obviously still below whatever threshold you used for filtering. When the median score in a cluster is 0.3, this is very suspicious, even though a single cell with 0.3 might not be. I'm not saying you should necessarily discard those clusters (although I normally do if they are very small), but you should at least treat them as highly suspicious.

One thing that worries me about your figure is that the vast majority of your nuclei seem to have a high doublet score (and a high % mito), although I'm not sure if this is simply due to overplotting. Could you maybe share a histogram of scDblFinder.score ?

[Captain-Pam]

Hi,

Thank you for your quick reply.

1. "delete the cells that are both double cells twice." I ran it twice as you mentioned, and one strategy was to remove cells that were doublets both times ((scDblFinder.class_1 == "doublets" & scDblFinder.class_2 == "doublets"). Then I found that the remaining some cells still have high scores, so wo applies another strategy (b>0.3 & a<0.01) | (b>0.6) .

2.The current threshold for mitochondria is 10%, but there are still many outlier cells on vlnplot. Should we choose a stricter threshold (8%)?

3. I ran two times and computed the mean scores. Here is their score histogram. Maybe I should choose a stricter threshold, like 0.2.

Thank you again.

[plger]

1. Sounds good.

So indeed the UMAP was misleading, and your doublet scores are very bimodal, as they should be.

[Captain-Pam]

Yes, you are right. umap was misleading.

1. I used the threshold scDblFinder.mean.score > 0.6 | (scDblFinder.mean.score > 0.3
   & scDblFinder.amulet.q.value < 0.01) . Is average score a good choice?I found that there are still doublets tags in the scDblFinder.class

2. Some of my samples are combined in a single library, but I'm currently using the default doublet rate (dbr.per1k.), which I believe is proportional to the total number of cells. Given this, is it still necessary to perform demultiplexing?

3. In addition, my siVI results are based on resolution 1, with the median nCounts in the middle of umap around 1000. Is this possible double cells? I see that its doublets score is not very high.

Thanks a Million