Error when running on a subset of 10x Xenium Human Ovarian Cancer

[VPetukhov]

I'm running ComSeg on a Xenium Human Ovarian Cancer dataset. Running on a full dataset (~120M molecules) takes forever, so as a test I made a subset of 290k molecules. Below is my code:

dict_scale = {"x": 0.2125, 'y': 0.2125, "z": 0.2125}

### create the dataset object
dataset = comseg.dataset.ComSegDataset(
    path_dataset_folder=path_dataset_folder,
    path_to_mask_prior=path_to_mask_prior,
    dict_scale=dict_scale,
    mask_file_extension=".tif",
    mean_cell_diameter=mean_cell_diameter,
    prior_name='in_nucleus'
)

dataset.add_prior_from_mask(overwrite=True)

dataset.compute_edge_weight( 
    images_subset=None,
    n_neighbors=40,
    sampling=True,
    sampling_size=10000
)

Comsegdict = dictionary.ComSegDict(
    dataset=dataset,
    mean_cell_diameter=mean_cell_diameter,
    community_detection="with_prior"
)

Comsegdict.compute_community_vector()

Comsegdict.compute_insitu_clustering(
    size_commu_min=3,
    norm_vector=True,
    ### parameter clustering
    n_pcs=3,
    n_comps=3,
    clustering_method="leiden",
    n_neighbors=20,
    resolution=1,
    n_clusters_kmeans=4,
    palette=None,
    nb_min_cluster=0,
    min_merge_correlation=0.8,
)

It works fine till the last function call (Comsegdict.compute_insitu_clustering), which fails with the following logs:

0%|                                                                                                                                                        | 0/1 [00:00<?, ?it/s]
improvement of modularity 0.04263385318290602
improvement of modularity 0.005007403551788303
improvement of modularity 0.0006468670214521133
improvement of modularity 2.5626005440537725e-05
[/home/vpetukhov/seg-errors/segmentation/env_sopa/lib/python3.12/site-packages/comseg/model.py:352](https://vscode-remote+ssh-002dremote-002b7b22686f73744e616d65223a22706c616e636b227d.vscode-resource.vscode-cdn.net/home/vpetukhov/seg-errors/segmentation/env_sopa/lib/python3.12/site-packages/comseg/model.py:352): FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  for node in self.community_anndata.obs["node_index"][comm_index]:
[/home/vpetukhov/seg-errors/segmentation/env_sopa/lib/python3.12/site-packages/comseg/model.py:354](https://vscode-remote+ssh-002dremote-002b7b22686f73744e616d65223a22706c616e636b227d.vscode-resource.vscode-cdn.net/home/vpetukhov/seg-errors/segmentation/env_sopa/lib/python3.12/site-packages/comseg/model.py:354): FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  nn_expression_vector = nn_expression_vector [/](https://vscode-remote+ssh-002dremote-002b7b22686f73744e616d65223a22706c616e636b227d.vscode-resource.vscode-cdn.net/) len(self.community_anndata.obs["node_index"][comm_index])
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [03:30<00:00, 210.56s/it]
Writing temporary files...
Running scTransform via Rscript...

Attaching package: arrow

The following objects are masked from package:feather:

    read_feather, write_feather

The following object is masked from package:utils:

    timestamp

Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 3784 by 18972
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 18972 cells
  |                                                                      |   0%Error in getGlobalsAndPackages(expr, envir = envir, globals = globals) : 
  The total size of the 19 globals exported for future expression (FUN()) is 1.21 GiB.. This exceeds the maximum allowed size of 500.00 MiB (option 'future.globals.maxSize'). The three largest globals are FUN (1.14 GiB of class function), umi_bin (72.60 MiB of class numeric) and data_step1 (526.63 KiB of class list)
 ... getGlobalsAndPackagesXApply -> getGlobalsAndPackages
Execution halted
Reading output files...

[tdefa]

Hello,
sorry for the late answer.

I think the problem comes from the scTransform R package. It may take too much memory when computing the normalisation parameters.
I added a default samling parameter sample_size_for_scTranform=5000, in compute_insitu_clustering so that the normalisation is computed on a smaller count matrix.

def compute_insitu_clustering(self, size_commu_min=3, norm_vector=False, # parameter clustering n_pcs=3, n_comps=3, clustering_method="leiden", n_neighbors=20, resolution=1, n_clusters_kmeans=4, palette=None, nb_min_cluster=0, min_merge_correlation=0.8, merge_cluster=True, sample_size_for_scTranform=5000, ):

Tell me if solve or not the problem

Maybe for large scale data it is more efficient to use sopa+comseg ?

[VPetukhov]

Thank you, I'll test it!

Maybe for large scale data it is more efficient to use sopa+comseg ?

I tried. :( It worked for ~20h and then failed with some incomprehensible errors. So I switched to a subset. And it's only 290k molecules, which is nothing. When I tested ComSeg on smaller gene panels it worked well though.

A bit off-topic, do you think scTransform is actually necessary here? :) scanpy has a similar method if a bit simpler. And it would make the package much more usable!

[tdefa]

Thanks for the suggestion! I wasn’t aware that it’s possible to compute Pearson residuals with Scanpy.

I tried. :( It worked for ~20h and then failed with some incomprehensible errors. So I switched to a subset. And it's only 290k molecules, which is nothing. When I tested ComSeg on smaller gene panels it worked well though.

Ok strange as I used the Xenium Human Ovarian Cancer datasets to test comseg+sopa+snakemake.

I will retry this week to run the code with the new versions of sopa and comseg. I will keep you updated if I find the problem or send you my set-up if I make it works.

[VPetukhov]

Thank you! Would you also be able to provide your code? We're running benchmarks, and that would help a ton.

[tdefa]

Hello,

I have indeed identified a non-deterministic bug in comseg+sopa+snakemake. Sometimes it works as expected, but other times it gets stuck at the resolution step. I am still investigating the issue and will keep you updated.

Additionally, I am preparing a new version of the package that includes the "scanpy residual normalization" you mentioned. I plan to release it in a few days after performing additional tests.

[tdefa]

In case you're interested, we have just developed a new segmentation method called RNA2seg. It takes both RNA point cloud data and staining as input, and it runs directly on spatial data. I hope it make RNA2seg much easier to use.

You can find the code here: GitHub - RNA2seg
Pre-print: bioRxiv Pre-print

It is not yet integrated in sopa