Why remove cells with 0 reads during LSI projection?

[dagarfield]

I think ArchR's implementation of iterative LSI is the greatest thing since sliced bread for feature selection in an ATAC dataset.
Looking through the code, I found something I don't quite understand in the ProjectLSI part of the code.

The LSI itself is constructed using top features, which leads to the (probably rare) situation that a cell could be a valid cell but still not have any reads in the top, say, 25,000 most variable features.

ArchR/R/IterativeLSI.R

Lines 1269 to 1273 in 968e442

	if(any(colSm == 0)){
	exclude <- which(colSm==0)
	mat <- mat[,-exclude]
	colSm <- colSm[-exclude]
	}

This is probably sensible, but does it lead downstream to the possibility that valid cells don't get incorporated into the LSI projection (and thus fall out of downstream steps like UMAP and clustering)?

This came up because I modified the code to better run on peak matrices. Support for my hackery is well outside the range of support, but curious if I am correct in this observation. Or if I've managed to fub up something.

[rcorces]

I think ArchR's implementation of iterative LSI is the greatest thing since sliced bread for feature selection in an ATAC dataset.

agreed

This is probably sensible, but does it lead downstream to the possibility that valid cells don't get incorporated into the LSI projection (and thus fall out of downstream steps like UMAP and clustering)?

This seems totally possible but I'm not sure how it would best be addressed. We could throw a warning saying how many cells are excluded and provide a recommendation to use more features?

This came up because I modified the code to better run on peak matrices.

If you ever want to contribute to making ArchR better, just let me know.

[dagarfield]

In a way, I think you did already deal with it (in a non-obvious way) -- there is a slot in the LSI object for excluded that lists all the cells left out. I'd just never seen it filled before I tried running LSI on a peak matrix (with tweaking), so wanted to make sure I hadn't missed a step (e.g. the program rejects "informative" features that don't have a read in all cells).

As for contributing, we should probably have that talk...probably when we invite you to give a talk as most of us here at GNE are now using ArchR for our scATAC....more soon.

[Xiaxhan]

I am curious about this.

If we use less varFeatures and some cells are left out in LSI projection. Although we even can list these cells, we need to handle them for downstream analysis - and I am not sure whether it would make some errors.

So maybe varFeatures should be recommended to set that each cell have reads in at least one top feature?

I tried this parameter to see how it affect batch effects - I should try other parameters and addHarmony if it cause cell lost. Although decreasing varFeatures is really helpful to reduce batch effects.

[rcorces]

I am not sure whether it would make some errors

I'm not sure either. would be good for you to report back on how this affects your downstream analysis. I think that this is handled correctly.

So maybe varFeatures should be recommended to set that each cell have reads in at least one top feature?

I think it is worth adding a warning during LSI for how many cells are removed. I think it is too analytically complex to try to automatically adjust the number of var features used to ensure that no cells are lost. I'm mostly concerned about edge cases where people have very lowly-sequenced or low quality data.

[dagarfield]

I have yet to actually encounter a case in which I lost cells during "normal" use of LSI (e.g. on the tile matrix). It was only messing about with edge cases and peak-call level signals

My own suggestions would be:

1. Probably a warning if/when it happens. I don't think it is a problem in most cases, and if/when it is, the user can recover them via projection (the LSI process actually does this internally). I would also guess that cells left out (esp. in the case of TileMatrix) are likely to be enriched for bad cells.

2. The option to call the LSI projection sub-function without discarding cells.

For context for anyone reading this, my own messing about with an edge case led to the removal of 321 cells out of ~250,000 cells. This is very much in the category of no big deal.

[Zepeng-Mu]

I had the same problem when running on 211571 cells and only got 207209 in the LSI when I used top 35,000 variable tiles. Those removed cells could be low quality cells, but it could be problematic for downstream analysis when the dimension of reducedDim and matrix do not match. Maybe one idea is to automatically use more tiles so all the cells are covered?

[Zepeng-Mu]

increasing varFeatures from 35000 to 50000 does not change the number of cells removed. So I'm not sure whether this is due to varFeatures.