RuntimeError when creating custom reference

[KaczorowskiLab]

Hello @danielsf,

Thank you for putting together this wonderful tool. For our analysis workflow, we would like to be able to create region specific reference based on the ABC WMB taxonomy. I followed the instructions to create the precomputed_stats.h5 file from a ABC release as documented here.

The inputs for the function include all of the original taxonomy CSV files, the region specific cell metadata file and region specific list of h5ad files. However, I am running into a RuntimeError raised by the validate_taxonomy_tree function.

RuntimeError: node CS20230722_CLUS_0001 (child of CCN20230722_SUPT:CS20230722_SUPT_0001 -- type <class 'str'>) is not present in the keys at level CCN20230722_CLUS e.g. of keys ['CS20230722_CLUS_0326', 'CS20230722_CLUS_0328', 'CS20230722_CLUS_0033', 'CS20230722_CLUS_0038', 'CS20230722_CLUS_0044']

I separately confirmed that the CSV file includes all the parent-child relationships and that level CCN20230722_CLUS does include the node CS20230722_CLUS_0001. I am a bit confused about the setup since the above function is getting the taxonomy tree input from the h5ad file which for the ABC releases are devoid of taxonomy information. Is there a work around for this error ? Thank you for your help!

[danielsf]

Hi @KaczorowskiLab

This is a use case I did not anticipate. I am assuming that the problem is that CS20230722_CLUS_001 does not appear in the region for which you are constructing a taxonomy. The way the code is set up, the nodes on leaf level of the taxonomy tree (in our case the CLUS_NNNN nodes) end up pointing to lists of the cells associated with that leaf node. If you create the taxonomy tree with cell_metadata_path = None, this step gets skipped, and each cluster points to an empty list. If you create the taxonomy tree with non-NULL cell_metadata_path, it is expected that every cluster will have some cells associated with it (if not, you would probably end up with trouble when you go to compute the mean expression values of cell type clusters, since there will be some divide-by-zero errors running around).

Unless you want to re-engineer cluster_annotation_term.csv and cluster_to_cluster_annotation_membership.csv to only include cell types that occur in the brain region you are working with, I think this is going to require a fix on our end.

What I don't know for sure (and would appreciate your feedback on) is: what is the right way to handle this? My inclination is to add a step to the taxonomy construction code that trims the tree, removing all cell types that have zero cells assigned to them (effectively pretending that cell types which occur outside of your region of interest do not exist). Is this acceptable to you, or do you need/want those cell types to persist in the taxonomy*?

*if they did persist in the taxonomy, I don't think that anything would ever get mapped to them, as they would all have NULL mean expression profiles.

Let me know what you think of this rough sketch of a solution.

[danielsf]

By the way: I am excited to see what comes of this analysis. Though I didn't anticipate this particular use case, that was just a deficiency in my imagination. What you are proposing ought to be possible in a perfect world.

[KaczorowskiLab]

Thank you for the detailed explanation and proposed solutions @danielsf ! I appreciate the clarification regarding the root cause of the issue. I will give re-engineering the cluster specific CSVs a shot.

Since the re-engineering the CSV files has constraints and potential problems, I find the second solution to trim trees of cell types with zero cells to be practical and efficient.

However, I have some concerns about the potential implications of trimming the tree, particularly in scenarios such as:

• Comparing different disease states or conditions where certain cell types might be present in one state but absent in another.
• Integrating our current dataset with future datasets that might include these missing cell types.
• Maintaining consistency across different analyses and regions.
• Biological relevance of some of the cell types even when not detected in the specific region or condition.

To mitigate these concerns, is a hybrid approach where the cell types with zero counts are flagged and excluded rather than permanently removed? This would provide the ability to retain the complete taxonomy structure for future use while avoiding immediate computational issues.

Please let me know your thoughts on this approach. Thank you for your assistance!

[danielsf]

Flagging and excluding empty cell types would take a lot more coding to implement, since several steps in the pipeline involve iterating over the taxonomy and all of its nodes. Each of those steps would have to be altered to check for the validity of the cell type.

In principle, we could just fill in the leaf level of the taxonomy with empty lists of cells wherever necessary. I think that will bog down the code, though, as the marker gene finder tries to find markers to distinguish between cell types with no data assigned to them. It would also undermine the advantage of doing a region-specific taxonomy (which, I assume, is finding marker genes that are targeted for discriminating only between cell types that exist in the chosen region), since the taxonomy would still appear to have all ~ 5,000 cell type clusters from the full Whole Mouse Brain data release.

This branch implements the solution I originally proposed (it had the virtue of only taking a few hours to code up). If you add --do_pruning True to your command line call to cell_type_mapper.cli.precompute_stats_abc, it will add a step before validating that taxonomy that removes any nodes that are not directly connected to lists of cells. Nothing about the inheritance structure of the existing taxonomy will be changed. It's just that cell types you don't "need" (because they aren't in your reference dataset) will be removed from the tree.

If you want a reference taxonomy that includes all of the nodes, you could do something like

from cell_type_mapper.taxonomy.taxonomy_tree import TaxonomyTree

tree = TaxonomyTree.from_data_release(
    cell_metadata_path='/path/to/full/cell_metadata.csv',
    cluster_annotation_path='/path/to/cluster_annotation_term.csv',
    cluster_membership_path='/path/to/cluster_to_cluster_annotation_membership.csv',
    hierarchy=['CCN20230722_CLAS', 'CCN20230722_SUBC', 'CCN20230722_SUPT', 'CCN20230722_CLUS'],
    do_pruning=False
)

# save the tree to a text file in JSON serialized form
with open('reference_taxonomy.json', 'w') as src:
    src.write(tree.to_str(indent=2, drop_cells=True))

The form of the taxonomy tree as serialized to reference_taxonomy.json is documented on this page.

Does this solution work for you?

[KaczorowskiLab]

Thanks @danielsf for the quick implementation! I will give this solution a try. I am curious to know if one would expect to see major differences going this route when compared to building a single region specific h5ad file and using the python -m cell_type_mapper.cli.precompute_stats_scrattch instead ?

[danielsf]

Building a region-specific h5ad file and calling precompute_stats_scrattch should achieve exactly the same result as the solution in the branch I pointed you towards.

precompute_stats_scrattch infers the cell type taxonomy tree from the annotations applied to each individual cell in the h5ad file's obs dataframe. Therefore, there cannot be any cell types with no cells associated with them. If no cells are annotated to belong to a cell type, that cell type will not appear in the inferred taxonomy tree.

The problem you've run into with precompute_stats_abc is that the inheritance structure of the taxonomy tree is defined in cluster_annotation_term.csv and cluster_to_cluster_annotation_membership.csv, while the assignment of cells to cell types is defined in cell_metadata.csv. There is technically nothing keeping those files in sync. At the point where you replaced the full cell_metadata.csv with single-region cell_metadata.csv, you brought the cell-to-cluster assignment out of sync with the cluster-to-supertype-to-subclass-to-class assignment. There were now clusters in cluster_annotation_term.csv and cluster_to_cluster_annotation_membership.csv that had not associated cells.

So: if you feel comfortable assembling a region-specific h5ad file in which the cell type annotations are in the obs dataframe, then: yes. That would be another way to work around this problem.

[danielsf]

Incidentally, this Jupyter notebook constructs an h5ad file with cartoon data and a cartoon taxonomy and then walks through all of the steps in the process of mapping data onto that cartoon taxonomy, starting with precompute_stats_scrattch, in case it is helpful.

[KaczorowskiLab]

Thank you so much for the quick clarification. I was able to complete the single h5ad route to use python -m cell_type_mapper.cli.precompute_stats_scrattch. However, since the WMB data is distributed across chemistries, in the cell_type_mapper.cli.map_to_on_the_fly_markers I see an option to pass multiple precomputed files (--reference_markers.precomputed_path_list). My current approach collapses all of the data into a single file, but it is possible to modify the code to create multiple precompute_stats files. Does the cell_type_mapper.cli.precompute_stats_abc with the newly implemented --do_pruning True also keep each region and chemistry files separate ? The 10Xv2 and 10Xv3 files are easier to manage but the 10XMulti files need an additional subsetting step since all regions are collapsed into one file.

[danielsf]

The new branch (precomputed_stats_abc with --do_pruning True) ought to keep the chemistry files separate.

[KaczorowskiLab]

Hi @danielsf,

After successful creation of region specific multiple chemistry files, I am running into an error that I am a bit confused about. I am running python -m cell_type_mapper.cli.map_to_on_the_fly_markers using the
--reference_markers.precomputed_path_list option to supply the list of all of the reference precomputed stat files.

Not supplying the --precomputed_stats.path in the above command gives me a ValidationError. When I do supply the same reference list to this argument, it complains that it is not a file. So based on this, can I only run pairwise mapping between the query and the precomputed files ?

[danielsf]

There is a subtlety that I obviously failed to capture in the documentation.

--reference_markers.precomputed_path_list and --precomputed_stats.path serve two different purposes. --reference_markers.precomputed_path_list points to the set of precomputed stats files which are used for finding marker genes. This is the step at which different chemistries can be kept separate. --precomputed_stats.path is used as the source for the mean gene expression profile against which cells are compared during the actual mapping process. At this mapping stage, there is no point in keeping the chemistries separate. The way the algorithm works, a cell type can only have one mean gene expression profile for purposes of mapping.

When you created the multiple precomputed_stats files, you should have seen files like

ouptut/directory/precomputed_stats.combined.h5
output/directory/precomputed_stats.WMB-10XMulti.h5
output/directory/precomputed_stats.WMB-10Xv2.h5
output/directory/precomputed_stats.WMB-10Xv3.h5

The WMB-10X* files are the files that should be passed to --reference_markers.precomputed_path_list. The *combined.h5 file is the file that should be passed to --precomputed_stats.path. It is a single file in which each cell type is represented by the chemistry in which it was most prevalent.

Sorry for the confusion. I'll try to put together an example notebook demonstrating this process (probably after SfN is past).

Quick note: the example notebook I pointed you towards intentionally sets n_processors to a very low value. This is because the data being processed is very small, and I want users to be able to run the example on their laptop. Marker gene selection can be very expensive, so I would recommend setting n_processors to as high as you can on whatever machine you are using. The back-of-the-envelope estimate is that, for the full Whole Mouse Brain taxonomy, marker gene selection took 3 hours with 32 cores and ~ 64 GB of memory. The process scales as N^2 where N is the number of leaf nodes in the taxonomy. I know you have drastically reduced this number from the full 5322, but I'm not sure by how much.

[KaczorowskiLab]

Thank you for the clarification! The errors and the description had me stumped. I am hopeful about the process running without any more errors! If you are open to it, I am also happy to contribute to the documentation to showcase the region specific reference file setup as an example.

[KaczorowskiLab]

@danielsf - Is there a threshold for minimum cells to be included in the reference files to lead to meangingful analysis? Another error I ran into which took a long time to debug was due to the sparse amount of cells in one of the chemistry files contributing to the region specific reference file. The precompute_stats_abc step ran without issues.

But when I move on to the mapping step using map_to_on_the_fly_markers, I ran into ValueError: All chunk dimensions must be positive. The initial error log pointed towards a problem with the query dataset. But I ran into the same problem using the reference_markers pipeline. This error stemmed from the _merge_sparse_by_pair_files function in markers.py which has a step to calculate chunk dimensions based on the n_up_indices or n_down_indices which ends up being 0 for the reference file with the sparse data. The problem was solved by excluding this file and re-running the pipeline. In this case, the extreme sparse data in this reference file made it easier to identify the problem, but in cases when the distinction is unclear, it would be helpful to know some metric or way to validate the reference data.

[danielsf]

Good debugging.

I don't have an off-the-cuff answer to the question that you have actually asked. I've run this code on taxonomies that contain leaf nodes with a few dozen cells in them and the code gave somewhat reasonable results. What ultimately matters is the statistical distribution of genes in the cells you do have. Can you quickly assess the number of cells in the least populous cluster? Since a variance needs to be calculated, a cluster with 2 or fewer cells in it will put some zeros in unfortunate (but not necessarily fatal) places. Note: as long as each cluster has a sensible (again, ~ a dozen) cell in at least one of the chemistry files, your results should be fine. The final step in mapping only considers the most populous chemistry file for each cluster.

On to the question you didn't ask: I've made a quick bugfix to the code in the trim/tree/240918 branch such that chunks is set to None if there is nothing to chunk in that code block. That ought to allow the code to run to completion, including the offending chemistry file. To a pull and try again. Let me know if I am wrong and the code still crashes.

[danielsf]

And thanks for the catch.