min-reads seems to have no effect in some cases

[LinearParadox]

Hi Simpleaf team,

Maybe I am misunderstanding something here, but it seems that min-reads is not actually filtering cells. For example, running a sample with the following config:

simpleaf quant --index GRCH38 \
        -1 r1.fastq.gz \
        -2 r2.fastq.gz \
        --threads 8 \
        -c 10xv2-5p \
        -o test_sample \
        -m GRCH38/ref/t2g_3col.tsv  \
        -r parsimony-em  \
        --anndata-out \
       -u --min-reads 10

However, when I open the resultant anndata file, there are cells with 0 counts. Confirmed using the following:

import scanpy as sc

adat = sc.read_h5ad("quants.h5ad")
adat.X.sum(axis=1).min()

outputs np.float64(0.0)

Looking at the docs, this iterates over mapped reads and selects barcodes occurring at least 10 times. Shouldn't this retain cells with at least one count, given the barcode is mapped?

[rob-p]

Hi @LinearParadox,

Great question, and I think there should probably be a better description of the effect of min-reads. It is not intended as e.g. a QC filter on cells, but rather as a way to determine what subset of the barcode universe should be considered as potentially valid (e.g. for the purpose of correction). Specifically, the min-read threshold applies to all reads having a barcode exactly, not just mapped reads. So, at least one way this might happen is if there are >= min-read unmapped reads exhibiting this barcode. In that case, the barcode is added to the list of barcodes to be quantified, but it may result in no assigned reads after the fact (which, I admit, is somewhat confusing).

This could be tested by looking at the unmapped_bc_count.bin in the mapping directory, which is just a binary file with (barcode, count) pairs written back to back. It's possible the offending barcode(s) appear in that list with a count >= min-count.

Best,
Rob