Probabilities not adding up to 1, and missing doublet-columns in prob_doublet.tsv.gz

[koefoeden]

Hi - I have two questions/issues that might be related, which is why I post them together here! And first of all, thanks for such a cool tool!

My first issue is that when plotting prob_max vs prob_doublet as shown below, there is a fraction of cells that have quite low doublet probabilities, while also having a near-zero prob_max (visualized in the orange circle). From my understanding of the calculation of prob_max and prob_doublet, this shouldn't be possible with only 8 donors as used here - since 8 times ~0 is of course not enough to account for the missing probability (i.e. (1-prob_doublet)). But perhaps the 8 single-probabilities and 64-doublet-pair probabilities does not need to add up to 1?

Secondly, and perhaps related, I looked into the prob_doublet.tsv.gz table as seen in the image below (I have obscurred column names for privacy-reasons). Here, it looks like there are many columns with only missing data (comb29 to comb56) and even columns missing, as seen from the total column count of 57, which should be 1+8^2 = 65. Do you know why this might be the case?

[koefoeden]

I did some further digging with a specific one of the problematic barcodes with prob_doublet = 0.25, and prob_max = ~0, and it turns out the summarized doublet-probability in the matrix row (even with missing cells) is correct, i.e. adding to ~1.

However, suspiciously, the max probability of these many doublet-combinations is 0.25 - exactly what was reported as prob_doublet in donor_ids.tsv. Is there a chance that the prob_doublet, as reported in donor_ids.tsv is in fact not a summarized doublet probability (as mentioned in the paper), but instead a max_doublet_prob? Maybe this is already mentioned elsewhere, and I simply missed this!

[xueyazhou2017]

@koefoeden I think the prob_doublet reported in the donor_ids.tsv file is the max. doublet_prob based on line 98 in utils/io_utils.py.

However, the second issue you mentioned in the earlier thread is most likely due to a bug. The total number of doublet pairs should be n*(n-1)/2, instead of n*(n-1) where n is the number of donors. In prob_doublet.tsv.gz file contains n*(n-1) ID pairs, but subsequent lines contains n*(n-1)/2 probabilities. This bug also makes IDs of best_doublet incorrect in the donor_ids.tsv file (although prob_doublet is still valid), so be cautious when you use that column. The I have submitted a PR on this issue.

[huangyh09]

Thanks @koefoeden, for sharing the issue. For issue 1, you are right: the max_doublet_prob is only the max but the sum of all doublet combinations. doublet_prob vector and singlet_prob vector should sum to 1, as shown these lines

For issue 2, @xueyazhou2017 has helped fix the bug by this commit - thanks. Now, it has been upgraded to 0.5.9, so you can upgrade via pip -U vireoSNP and re-run your experiment. This will only affect the best doublet combination, and no other results will change.

[koefoeden]

Hi @huangyh09 - regarding issue 1, the problem is exactly that it's not called max_doublet_prob in the donor_ids.tsv file, but just doublet_prob, even though it's using a max.

I went ahead and submitted a pull request #116 using 'sum' instead of 'max' so the naming is more consistent - feel free to decline it if there is something I've misunderstood!