search

cistrome / MIRA

Python package for analysis of multiomic single cell RNA-seq and ATAC-seq.
60 stars 8 forks source link

LinAlgError: Singular matrix for some of the data #42

Open Yansr3 opened 9 months ago

Yansr3 commented 9 months ago

I was trying to carte lineage tree and calculate branch probabilities with mira.time.get_branch_probabilities using terminal cells detected from mira.time.find_terminal_cells. I got the error of "LinAlgError: Singular matrix" for some of my data. May I ask for some advice for creating the tree?

terminal_cells = mira.time.find_terminal_cells(rna_main, threshold = 1e-3)
mira.time.get_branch_probabilities(rna_main, 
                                   terminal_cells= {'branch1' : terminal_cells[1],
                                                    'branch2' : terminal_cells[2],
                                                    'branch3' : terminal_cells[0]})
---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
Cell In[176], [line 1](vscode-notebook-cell:?execution_count=176&line=1)
----> [1](vscode-notebook-cell:?execution_count=176&line=1) mira.time.get_branch_probabilities(rna_main, 
      [2](vscode-notebook-cell:?execution_count=176&line=2)                                    terminal_cells= {'branch1' : terminal_cells[1],
      [3](vscode-notebook-cell:?execution_count=176&line=3)                                                     'branch2' : terminal_cells[2],
      [4](vscode-notebook-cell:?execution_count=176&line=4)                                                     'branch3' : terminal_cells[0]})

File [~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/adata_interface/core.py:73](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/adata_interface/core.py:73), in wraps_functional.<locals>.run.<locals>._run(adata, *args, **kwargs)
     [68](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/adata_interface/core.py:68)         raise TypeError('{} is not a valid keyword arg for this function.'.format(kwarg))
     [70](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/adata_interface/core.py:70) #print(function_kwargs)
     [71](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/adata_interface/core.py:71) #print(fetch(None, adata, **getter_kwargs))
---> [73](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/adata_interface/core.py:73) output = func(**fetch(None, adata, **getter_kwargs), **function_kwargs)
     [74](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/adata_interface/core.py:74) #print(output, adata, adder_kwargs)
     [75](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/adata_interface/core.py:75) return add(adata, output, **adder_kwargs)

File [~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/pseudotime/pseudotime.py:560](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/pseudotime/pseudotime.py:560), in get_branch_probabilities(transport_map, terminal_cells)
    [557](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/pseudotime/pseudotime.py:557) # Fundamental matrix
    [558](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/pseudotime/pseudotime.py:558) mat = np.eye(Q.shape[0]) - Q
--> [560](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/pseudotime/pseudotime.py:560) N = inv(mat)
    [562](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/pseudotime/pseudotime.py:562) # Absorption probabilities
    [563](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/mira/pseudotime/pseudotime.py:563) branch_probs = np.dot(N, transport_map[trans_states, :][:, absorbing_states_idx])

File [~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/numpy/linalg/linalg.py:561](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/numpy/linalg/linalg.py:561), in inv(a)
    [559](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/numpy/linalg/linalg.py:559) signature = 'D->D' if isComplexType(t) else 'd->d'
    [560](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/numpy/linalg/linalg.py:560) extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> [561](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/numpy/linalg/linalg.py:561) ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    [562](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/numpy/linalg/linalg.py:562) return wrap(ainv.astype(result_t, copy=False))

File [~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/numpy/linalg/linalg.py:112](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/numpy/linalg/linalg.py:112), in _raise_linalgerror_singular(err, flag)
    [111](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/numpy/linalg/linalg.py:111) def _raise_linalgerror_singular(err, flag):
--> [112](https://file+.vscode-resource.vscode-cdn.net/Users/siruiyan/Workspace/EryProject/mira/trimod_MIRA_rna_atac_totalVI_mito20/lineage_tree/~/miniconda3/envs/Mira-2.1.1a4-fa2/lib/python3.11/site-packages/numpy/linalg/linalg.py:112)     raise LinAlgError("Singular matrix")

LinAlgError: Singular matrix
craniolab commented 7 months ago

I'm currently experiencing the same error and I have not yet been able to determine what is causing it.

AllenWLynch commented 7 months ago

Hi all, can you provide the scipy and mira versions which you are using?

Sometimes this issue arises because of scipy, sometimes because of the underlying topology of the diffusion graph. If you try selecting different start and end cells (by choosing a nearest neighbor, perhaps), does that work?

craniolab commented 7 months ago

This happened with scipy 1.10.1 and mira 2.1.1

Changing the start and end cells didn't resolve the issue. In the end I tried switching to the pseudo-inverse of the matrix (as palantir has a similar exception), I assumed it was Hermitian so I used the pinvh function from scipy (numpy pinv kept crashing my pc) , which produced some results.

from numpy.linalg import inv, LinAlgError from scipy.linalg import pinvh try: N = inv(mat) except LinAlgError: N = pinvh(mat)

Since this part is based on palantir I also used ran it with the mira topic model results, and the palantir results were quite different from what I got with mira, which makes me wonder what I'm doing wrong.

Honestly, none of this math stuff is my forte, so I don't understand it much.

Yansr3 commented 7 months ago

I'm also using scipy 1.10.1 and mira 2.1.1. It worked for all my cells, and some subsets of the cells. I encountered this problem for some subsets of cells.