Implement the the order of magnitude algorithm from Cell Ranger to emptyDropsCellRanger

[an-altosian]

Hi @LTLA,

Happy New Year!! Here is Dongze using my Altos's GitHub account. Finally I got some cycles to address #88 (and #109 as a bonus ;P).

According to CellRanger's GitHub repo, I implemented the core functions for the order of magnitude algorithm to estimate the n.expected.cells parameter if it is set as NULL (the default). I have tested that the R implementation resulted in the same number as the python implementation in CellRanger if ignoring the effects of random seeds.

The only concern I have is that the current implementation uses population variance in this line (or here in Cell Ranger). Please feel free to change this to sample var (var(top_n_boot)) if you think it is more appropriate.

Best,
Dongze

[LTLA]

After reading this for half an hour, I have no idea what this is doing. Why is there a need for bootstrapping? Why can't the threshold be found exactly?

[an-altosian]

I guess the bootstrapping is saying those barcodes represent a tiny fraction of the population, so they want to use bootstrapping to get a "robust" threshold that is a "guaranteed approximation"?

Haha, I guess you found the reason why I did not copy and paste the function you wrote in #119 but spent a day translating their python functions into R. 😉

I have tested that the numbers I got using the function you wrote were only tens of cells' different than their implementation. It is totally fine if we just apply your succinct and straightforward function, but I guess there will be issues popping up requesting "Cell Ranger"'s algorithms.

[LTLA]

Hm. I'm caught between two concerns:

1. Reproducing the CellRanger calls is desirable.
2. But I don't want to maintain that code or keep it up to date with new CellRanger versions.

So how about this. We'll spin off the existing existing emptyDropsCellRanger function in a new Bioconductor package that you (and possibly other interested parties) can maintain. Then I don't have to worry about staying up to date with CellRanger, and you can achieve perfect parity at your own discretion.

Depending on the scope of this new package, I am also willing to offer the other 10X-specific functions, e.g., get10xMolInfoStats, read10xMolInfo. These functions also cause maintenance headaches for me whenever 10X decides to change the CellRanger format, so if someone is actively interested in staying aligned with CellRanger, they would be better suited to as maintainers. (In fact, I do have https://github.com/LTLA/crio, which was my plan for isolating these headaches in a separate package, but I ran out of motivation.)

[an-altosian]

I am totally fine with the plan. I will try my best to help implement any future changes in these functions!

[LTLA]

Alright, let me ping BioC slack and see if there are any other interested parties, then I will start the migration process.

[LTLA]

Well, I guess there aren't any more interested parties. Can you set up a repo and put me as a collaborator, then I'll move over some code.

[DongzeHE]

Haha it's OK. I will work on this soon!

[LTLA]

Sorry for the lateness @DongzeHE, I just remembered to deal with this.

I set up https://github.com/LTLA/cram (placeholder name, actual name is up to you). I've already moved your emptyDropsCellRanger code there and it passes CHECK. Currently it still has a few ::: dependencies on DropletUtils but we can get rid of these later.

Let me know where you want me to transfer the repo. Then you add your new stuff, submit it to BioC, etc. Happy to help in that regard. If we are quick, you might be able to get it into the next release (3.22).

[DongzeHE]

Hi @LTLA , wonderful! I will find some time ASAP to move things!

BTW, we also used what were discussed our tool called QCatch for QCing alevin-fry's quantification results. Your name is in the acknowledgment section of the preprint 🫶🫶🫶