Loess error in compute-score after subsampling

[twoneu]

Hi scDRS team,

I am encountering a loess error when trying to run a subsampled, lognormalized h5ad file through scdrds compute-score:

Call: scdrs compute-score
--h5ad-file ordovas_immune_downsampled.h5ad
--h5ad-species human
--cov-file None
--gs-file all.gs
--gs-species human
--ctrl-match-opt mean_var
--weight-opt vs
--adj-prop None
--flag-filter-data False
--flag-raw-count False
--n-ctrl 1000
--flag-return-ctrl-raw-score False
--flag-return-ctrl-norm-score True
--out-folder scDRS/immune

Loading data:
--h5ad-file loaded: n_cell=10727, n_gene=5000 (sys_time=0.2s)
First 3 cells: ['N51.LPB.TCAGCAACATACGCTA-0', 'N111.LPA2.CTTCTCTGTGGCAAAC-0', 'N20.LPA.TGGATGTGCACTTT-0']
First 5 genes: ['A2M', 'A2M-AS1', 'A2ML1', 'AAED1', 'AAK1']
--gs-file loaded: n_trait=6 (sys_time=0.2s)
Print info for first 3 traits:
First 3 elements for 'ukbb_CRC': ['CLTCL1', 'DAB2IP', 'NBPF1'], [4.0019, 3.5014, 3.4717]
First 3 elements for 'ukbb_IBD': ['ZNF236', 'C3orf38', 'LAMB2'], [3.17, 3.1623, 3.1358]
First 3 elements for 'alkes_UC': ['FCGR2A', 'IL23R', 'DLD'], [7.964, 7.0309, 6.897]

Preprocessing:
/Users/user/opt/anaconda3/lib/python3.9/site-packages/scdrs/pp.py:323: RuntimeWarning: invalid value encountered in log10
x = np.log10(df_gene["ct_mean"].values[not_const])
Traceback (most recent call last):
File "/Users/user/opt/anaconda3/bin/scdrs", line 723, in
fire.Fire()
File "/Users/user/opt/anaconda3/lib/python3.9/site-packages/fire/core.py", line 141, in Fire
component_trace = _Fire(component, args, parsed_flag_args, context, name)
File "/Users/user/opt/anaconda3/lib/python3.9/site-packages/fire/core.py", line 466, in _Fire
component, remaining_args = _CallAndUpdateTrace(
File "/Users/user/opt/anaconda3/lib/python3.9/site-packages/fire/core.py", line 681, in _CallAndUpdateTrace
component = fn(*varargs, **kwargs)
File "/Users/user/opt/anaconda3/bin/scdrs", line 211, in compute_score
scdrs.preprocess(adata, cov=df_cov, adj_prop=ADJ_PROP, n_mean_bin=20, n_var_bin=20, copy=False)
File "/Users/user/opt/anaconda3/lib/python3.9/site-packages/scdrs/pp.py", line 205, in preprocess
df_gene, df_cell = compute_stats(
File "/Users/user/opt/anaconda3/lib/python3.9/site-packages/scdrs/pp.py", line 325, in compute_stats
model.fit()
File "_loess.pyx", line 899, in _loess.loess.fit
ValueError: b'Extrapolation not allowed with blending'

I used the following code to generate the subsampled dataset from the original data:

crc_immune_adata = sc.read_h5ad("ordovas_immune_processed.h5ad")
target_cells = 250

crc_immune_adata_ds = [crc_immune_adata[crc_immune_adata.obs.Cluster == clust] for clust in crc_immune_adata.obs.Cluster.cat.categories]

for dat in crc_immune_adata_ds:
    if dat.n_obs > target_cells:
         sc.pp.subsample(dat, n_obs=target_cells)

crc_immune_adata_downsampled = crc_immune_adata_ds[0].concatenate(*crc_immune_adata_ds[1:])
crc_immune_adata_downsampled.write("ordovas_immune_downsampled_250.h5ad")

The subsampled dataset is attached here:
ordovas_immune_downsampled_250.h5ad.zip

Thanks for all your time and help!

[KangchengHou]

/Users/user/opt/anaconda3/lib/python3.9/site-packages/scdrs/pp.py:323: RuntimeWarning: invalid value encountered in log10
x = np.log10(df_gene["ct_mean"].values[not_const])

looks like there is a log(0) problem (which is quite strange). Somehow I am not able to read h5ad file. Could you send a raw text file (perhaps .csv)?

[twoneu]

ordovas_immune_downsampled.zip
Hopefully this works! I'm very sorry if it ends up being a data format issue - I really appreciate your time and patience.

[KangchengHou]

@twoneu I was being confusing. I was referring to the expression matrix in .csv format (rows are cell, columns are genes). I can not open the h5ad matrix somehow.

Using df = pd.DataFrame(data=adata.X, index=adata.obs_names, column=adata.var_names) should work.

PS. within scdrs, there is a step for calculating the ct_mean (gene-wise count matrix average) where scdrs transform your data in log-transformed space to count space using np.expm1. One thing to check is to make sure your log-transformed expression values are all >= 0.

[twoneu]

Thank you so much - it turns out there was one NaN within my matrix that was throwing the error. This package is great and so are you!

[KangchengHou]

that's great! we will do some more checking on the input matrix in next version (including NaN, expression values >= 0 for log-transformed matrix).