Pseudotime computation fails when cells_back has fewer columns than boo

[tianxin-wang]

Hi Louis,

When using scFates.tl.pseudotime() with n_map > 100 on complex branching datasets, I consistently encountered the following error:

It appears that cells_back has fewer columns than boo, due to the fact that some segments are never assigned in the multiple probabilistic mappings (n_map > 1). This shape mismatch causes the DataFrame multiplication to fail. However, the tree learned by scf.tl.tree() isn't perfect no matter how I tune the hyperparameters; not all cells get assigned to all possible segments.

Would it be possible to enforce the full segment space when building cells_back, e.g., cells_back = cells_back.reindex(columns=full_segment_list, fill_value=0), to avoid this shape mismatch??

Thank you,
Cynthia

[tianxin-wang]

Similarly, I've run into an issue where adata.obs['milestones'] contains NaN values and cannot be directly cast to int.

The root cause seems to be that not all segments (pp_seg.n) end up assigned to any cells — which leads to NaN values in the milestones Series. These NaNs then cause the .astype(int) conversion to fail.

[LouisFaure]

Hi Cynthia,

Following your remarks and suggestions, I have prepared a new version that addresses your issues. Could you to run with the following version:

pip install git+https://github.com/LouisFaure/scFates@v1.1.2dev

The first issue is indeed due to a segment not having cells, when taking into account all the mappings, now allsegs contains all segments as you suggested. The second error is more likely due to when a segment is composed of a single cell, rendering milestone estimation impossible (as it requires at least two cells). In that case, these cells are assigned a NaN value in the final obs.milestones.

In both cases, I have added warnings to let the user know of these occurrences. Let me know if that solves your issue; in such case, I will release these fixes into a new version!

[tianxin-wang]

Hi Louis,

Thank you so much for your prompt response! I have tried again which ran into the following problem. Would it be solved by first importing log, then log.info?

[LouisFaure]

It seems that my tests never caught any situation requiring to call logging in palette management, I have added the import in the latest commit of v1.1.2dev, let me know if this works now

[tianxin-wang]

Hi Louis, it now works without issues - thank you so much!

[tianxin-wang]

I encountered another maybe related issue while using scf.tl.fit():

The error originates from this line in fit():
Xgenes = get_X(adata, subtree.index, features, layer, togenelist=True), and the suspected cause is the subtree.index. subtree is generated from getpath(img, root, tips, tip, graph, df). Looking further, the issue seems to originate from getpath(), where gatherpath() may return None if a tip has no path or the segment has no mapped cells. This leads to a list containing None, which causes pd.concat() in fit() to fail downstream, raising obscure errors like "negative dimensions are not allowed."

Would it be possible to return empty DataFrames instead of None, or explicitly filter out None values before concatenation?

Let me know what you think, and thanks again for all your work on scFates!

[LouisFaure]

I have implemented what you asked on the latest commit, but I am not sure this will ensure the rest of downstream task to run successfully, as this is quite an edge case that was not initially encountered

[NicolasSompairac]

Hi @LouisFaure !

I've run into a similar issue mentioned in this thread. I've tried installing the version your provided but while fixing the pandas cohersion problem it seems to have created another.

Here is the error log I got:

Any idea how to get around that?

[LouisFaure]

Hi @NicolasSompairac
Does your tree contains single edge connecting two forks? a quick way to check this it to run the following:

B = adata.uns["graph"]["B"]

g = igraph.Graph.Adjacency((B > 0).tolist(), mode="undirected")
degrees = g.degree()
high_degree_nodes = [i for i, d in enumerate(degrees) if d >= 3]

for edge in g.es:
    u, v = edge.tuple
    if u in high_degree_nodes and v in high_degree_nodes:
        print((u, v))

[janitorhumor]

I encountered same problem and then updated to the dev verson 1.1.2. Then I got a similar error:

~/miniconda3/envs/scfates/lib/python3.11/site-packages/scFates/tools/pseudotime.py:161: UserWarning: Some segs have no cell assigned: ['1', '39', '43', '57', '58']
  warnings.warn(message)
~/miniconda3/envs/scfates/lib/python3.11/site-packages/scFates/tools/pseudotime.py:248: UserWarning: Some cells have no milestones assigned. This is likely due to the fact that these uniquely compose a segment.
  warnings.warn(message)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "~/miniconda3/envs/scfates/lib/python3.11/site-packages/scFates/tools/pseudotime.py", line 270, in pseudotime
    adata.uns["seg_colors"] = [
                              ^
  File "~/miniconda3/envs/scfates/lib/python3.11/site-packages/scFates/tools/pseudotime.py", line 271, in <listcomp>
    np.array(adata.uns["milestones_colors"])[
IndexError: index 0 is out of bounds for axis 0 with size 0

I checked the connection between highly branched nodes, as you mentioned, and the answer is yes.

(0, 1)
(1, 6)
(4, 51)
(11, 14)
(11, 42)

I guess a graph-like structure would surely bring some problem to the pseudotime procedure, so copy the parameters of scf.tl.tree as provided in Tree_Analysis_Bone_marrow_fates can be a bad idea. I tried to lessen node, turn up regularization parameter, etc. However, I cannot solve the problem, and there are always fork-to-fork egdes.

Updated:
Basically the line 271 warning is triggered when adata.uns["graph"]["pp_seg"].to doesn't match adata.uns["graph"]["milestones"], but I don't see the influence on downstream analysis. So I tentatiely rewrite scf.tl.pseudotime code to:

try:
    print("Trying to assign color to milestone.")
    mst_dict = adata.uns["graph"]["milestones"]  # dict
    milestone_colors = np.array(adata.uns["milestones_colors"])
    
    # 先构建 value->index 的字典，避免每次都调用 list.index()
    val_to_idx = {v: i for i, v in enumerate(mst_dict.values())}
    
    adata.uns["seg_colors"] = [
        milestone_colors[val_to_idx[t]] if t in val_to_idx else "#808080"
        for t in adata.uns["graph"]["pp_seg"].to
    ]
    grey_count = sum(c == "#808080" for c in adata.uns["seg_colors"])
    print(f"Successfully assigned colors. Grey fallback used {grey_count} times.")
    
except Exception as e:
    print(f"An error occurred: {e}, please assign color manually at adata.uns['seg_colors']")

... making it easier at least to manipulate with the adata manually.

In case anyone would be interested in a look-alike error in scf.pl.milestones (see following), which intrinsically related to line 248 warning, I recomend to run a scf.tl.cleanup first, that would prune spurious branches, such as a segment composed by only few cells.

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/data/HeLab/bio/IBD_analysis/src/ScanpyTools/Scanpy_Plot.py", line 46, in __call__
    self.func(*args, **kwargs)
  File "/data/HeLab/miniconda3/envs/scfates/lib/python3.11/site-packages/scFates/plot/milestones.py", line 86, in milestones
    c=gen_milestones_gradients(adata)[order].values,
      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/HeLab/miniconda3/envs/scfates/lib/python3.11/site-packages/scFates/plot/utils.py", line 208, in gen_milestones_gradients
    return pd.concat(list(map(milestones_prog, seg_order)))
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/HeLab/miniconda3/envs/scfates/lib/python3.11/site-packages/scFates/plot/utils.py", line 193, in milestones_prog
    np.argwhere(adata.obs.milestones.cat.categories == mfrom)[0][0]
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^
IndexError: index 0 is out of bounds for axis 0 with size 0

Otherwise you can still dodge this by simply remove these cells, normally there are not too many of them:
adata_notna = adata_pseudo[adata_pseudo.obs['milestones'].isna()==False]

Segment with too little cells also triggers another error during scf.tl.test_association:

...
  File "/data/HeLab/miniconda3/envs/scfates/lib/python3.11/site-packages/scFates/tools/test_association.py", line 362, in test_assoc
    mdl = list(map(gamfit, sdf.seg.unique()))
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/HeLab/miniconda3/envs/scfates/lib/python3.11/site-packages/scFates/tools/test_association.py", line 357, in gamfit
    m = rmgcv.gam(
        ^^^^^^^^^^
  File "/data/HeLab/miniconda3/envs/scfates/lib/python3.11/site-packages/rpy2/robjects/functions.py", line 208, in __call__
    return (super(SignatureTranslatedFunction, self)
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/HeLab/miniconda3/envs/scfates/lib/python3.11/site-packages/rpy2/robjects/functions.py", line 131, in __call__
    res = super(Function, self).__call__(*new_args, **new_kwargs)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/HeLab/miniconda3/envs/scfates/lib/python3.11/site-packages/rpy2/rinterface_lib/conversion.py", line 45, in _
    cdata = function(*args, **kwargs)
            ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/HeLab/miniconda3/envs/scfates/lib/python3.11/site-packages/rpy2/rinterface.py", line 817, in __call__
    raise embedded.RRuntimeError(_rinterface._geterrmessage())
rpy2.rinterface_lib.embedded.RRuntimeError: Error in (function (formula, family = gaussian(), data = list(), weights = NULL,  : 
  Not enough (non-NA) data to do anything meaningful
...

Although I am not very familiar with GAM, but in my cases the easiest way to avoid this error is again to remove segments, but this time more strict than above, for the minimal count of cells in each branch should not be less than spline_df.

def filter_adata_by_seg(adata, spline_df=5, seg_col="seg"):
    """
    根据每个 seg 的细胞数过滤 AnnData 对象，并打印筛选前后细胞数量变化。

    Parameters
    ----------
    adata : AnnData
        输入的 AnnData 对象，要求 .obs 里有 seg 列。
    min_cells_per_seg : int
        每个 seg 至少保留多少细胞，否则丢弃。
    seg_col : str
        .obs 中存放 seg 信息的列名，默认 "seg"。

    Returns
    -------
    AnnData
        过滤后的 AnnData。
    """
    import pandas as pd
    
    # pro-filter info
    total_before = adata.n_obs
    value_df = adata.obs[seg_col].value_counts()
    print(f"Before filtering: total cells = {total_before}")
    print(f"Segment counts:\n{value_df}")
    
    # filter
    keep_segs = value_df.index[value_df.values >= spline_df]
    adata_filtered = adata[adata.obs[seg_col].isin(keep_segs)].copy()
    
    # after-filter info
    total_after = adata_filtered.n_obs
    print(f"\nAfter filtering: total cells = {total_after} (removed {total_before - total_after})")
    value_df_after = adata_filtered.obs[seg_col].value_counts()
    print(f"Remaining segment counts:\n{value_df_after}")
    
    return adata_filtered

→ An extra tip here: scf.tl.cleanup may cause nan line in adata.obsm["X_R"], when some cells are isolated (only have non-zero soft assignment value on principal points). Deleting those principal points results in full zero rows in the matrix, thus a safer normalisation is recomended. Instead of:
(Line 122, scf.tl.cleanup)

R = (R.T / R.sum(axis=1)).T

Try:

row_sums = R.sum(axis=1, keepdims=True)
row_sums[row_sums == 0] = 1  # 防止除零
R = R / row_sums

Nan lines don't trigger any warning but impede scf.pl.graph generating nodes properly.

[MariaRosariaNucera]

It seems that my tests never caught any situation requiring to call logging in palette management, I have added the import in the latest commit of v1.1.2dev, let me know if this works now

is this scFates-1.1.1.post4+git.6a352310 the same version?