Allow only high confidence filtering when reading Cell Ranger output

[ktpolanski]

Description of feature

There are situations when it's desirable to forego Cell Ranger's cell calling, and the VDJ equivalent of that is to go back to all_contig_annotations and subset accordingly. Scirpy supports ingesting the file, but its subsequent filtering option enforces both is_cell and high_confidence. It would be convenient to be able to look at just high_confidence, leaving any errant empty droplets to get filtered out when calling pp.merge_with_ir() later.

[grst]

Hi @ktpolanski,

can't you achieve this by using read_10x_vdj(..., filtered=False) and then subsetting by the high_confidence column in .obs?

[ktpolanski]

I did not spot that column in .obs before, sorry about that. However, high_confidence is a per-contig call, the values are not guaranteed to be consistent across a cell's contig pool. Here's an example:

$ grep "TTTGTCATCAGAAATG-1" all_contig_annotations.csv  | cut -f 1-4 -d ,
TTTGTCATCAGAAATG-1,false,TTTGTCATCAGAAATG-1_contig_1,false
TTTGTCATCAGAAATG-1,false,TTTGTCATCAGAAATG-1_contig_2,false
TTTGTCATCAGAAATG-1,false,TTTGTCATCAGAAATG-1_contig_3,true

Attempting to read the file with filtered=False ends up throwing an understandable error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/jovyan/my-conda-envs/d2c/lib/python3.8/site-packages/scirpy/io/_io.py", line 290, in read_10x_vdj
    return _read_10x_vdj_csv(path, filtered, include_fields)
  File "/home/jovyan/my-conda-envs/d2c/lib/python3.8/site-packages/scirpy/io/_io.py", line 239, in _read_10x_vdj_csv
    ir_obj.add_chain(chain_dict)
  File "/home/jovyan/my-conda-envs/d2c/lib/python3.8/site-packages/scirpy/io/_datastructures.py", line 139, in add_chain
    self[tmp_field] = chain.pop(tmp_field)
  File "/home/jovyan/my-conda-envs/d2c/lib/python3.8/site-packages/scirpy/io/_datastructures.py", line 113, in __setitem__
    raise ValueError(
ValueError: Cell-level attributes differ between different chains. Already present: `False`. Tried to add `True`.

[grst]

I see, so I had wrong assumptions about the high_confidence field when implementing that. At least the datastructure caught it.

I am working on a new datastructure for scirpy (#327 / #356). This will make it easier to read in all chains and use custom filtering strategies (e.g. keep non-productive chains etc).

[ktpolanski]

Sure, so hopefully this works itself nicely into that. Thanks!

[ktpolanski]

By the way, I took this as an excuse to acquaint myself a little more with awk to make a workaround before this is rolled out:

awk 'BEGIN {IGNORECASE=1; FS=","}; NR==1 {print}; $4=="true" {print}' all_contig_annotations.csv > hiconf.csv

This can then be read in via scirpy.io.read_10x_vdj("hiconf.csv", filtered=False) for the desired outcome.

[grst]

Hi @ktpolanski,

I implemented a solution using the new data structure which is now available from the latest development version. Would be glad if you could test it!

1. Install the dev version:

   pip install git+https://github.com/scverse/scirpy.git

2. You should now be able to read in the 10x data with filtered=False. This will keep all cells and store all chains in adata.obsm["airr"].

   ir.io.read_10x_vdj("all_contig_annotations.csv", filtered=False)

3. Then, there is now the new index_chains step (explanation | function doc) that will create a reference to the chains that are actually used in the analysis (while still keeping all chains around). By default, it filters out chains that are marked as non-productive and don't have a CDR3 amino acid sequence. But it is possible to customize that using the filter attribute and a callback function:

   ir.pp.index_chains(adata, filter = lambda x: x["high_confidence"] is True)

   or, if you want to keep the "productive" and "CDR3" filter

   ir.pp.index_chains(
       adata,
       filter = [
           "productive", 
           "require_junction_aa",
           "lambda x: x["high_confidence"] is True
       ]
   )

I hope it is clear how this works. Let me know if not or if you have any suggestions!

Additionally, I was wondering if in the read_10x_vdj, filtered=False should become the new default. It would be a breaking change though.