wangjiawen2013 commented 4 years ago

Hi, how to use binarization.py in jupyter notebook and shell ? I ran auc_mtx_binary = binarize(auc_mtx,num_workers=15) in jupyter notebook Some warnings and errors occured: /home/wangjw/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/diptest.py:30: RuntimeWarning: divide by zero encountered in true_divide slopes = (work_cdf[1:] - work_cdf[0]) / distances /home/wangjw/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/diptest.py:30: RuntimeWarning: divide by zero encountered in true_divide slopes = (work_cdf[1:] - work_cdf[0]) / distances /home/wangjw/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/diptest.py:30: RuntimeWarning: divide by zero encountered in true_divide slopes = (work_cdf[1:] - work_cdf[0]) / distances /home/wangjw/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/diptest.py:30: RuntimeWarning: divide by zero encountered in true_divide slopes = (work_cdf[1:] - work_cdf[0]) / distances RemoteTraceback Traceback (most recent call last) RemoteTraceback: """ Traceback (most recent call last): File "/home/wangjw/bin/miniconda3/envs/scenic/lib/python3.6/multiprocessing/pool.py", line 119, in worker result = (True, func(*args, **kwds)) File "/home/wangjw/bin/miniconda3/envs/scenic/lib/python3.6/multiprocessing/pool.py", line 47, in starmapstar return list(itertools.starmap(args[0], args[1])) File "/home/wangjw/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/binarization.py", line 52, in derivethreshold if not isbimodal(data, method): File "/home/wangjw/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/binarization.py", line 43, in isbimodal , pval, = diptst(np.msort(data)) File "/home/wangjw/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/diptest.py", line 53, in diptst d, (, idxs, left, , right, ) = dip_fn(dat, is_hist) File "/home/wangjw/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/diptest.py", line 115, in dip_fn xl = left_touchpoints[d_left == left_diffs][0] IndexError: index 0 is out of bounds for axis 0 with size 0 """

The above exception was the direct cause of the following exception:

IndexError Traceback (most recent call last) in ----> 1 auc_mtx_binary = binarize(auc_mtx,num_workers=15)

~/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/binarization.py in binarize(auc_mtx, threshold_overides, seed, num_workers) 77 thrs = p.starmap( derive_threshold, [ (auc_mtx, c, seed) for c in auc_mtx.columns ] ) 78 return pd.Series(index=auc_mtx.columns, data=thrs) ---> 79 thresholds = derive_thresholds(auc_mtx) 80 if threshold_overides is not None: 81 thresholds[list(threshold_overides.keys())] = list(threshold_overides.values())

~/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/binarization.py in derive_thresholds(auc_mtx, seed) 75 def derive_thresholds(auc_mtx, seed=seed): 76 with Pool( processes=num_workers ) as p: ---> 77 thrs = p.starmap( derive_threshold, [ (auc_mtx, c, seed) for c in auc_mtx.columns ] ) 78 return pd.Series(index=auc_mtx.columns, data=thrs) 79 thresholds = derive_thresholds(auc_mtx)

~/bin/miniconda3/envs/scenic/lib/python3.6/multiprocessing/pool.py in starmap(self, func, iterable, chunksize) 272 func and (a, b) becomes func(a, b). 273 ''' --> 274 return self._map_async(func, iterable, starmapstar, chunksize).get() 275 276 def starmap_async(self, func, iterable, chunksize=None, callback=None,

~/bin/miniconda3/envs/scenic/lib/python3.6/multiprocessing/pool.py in get(self, timeout) 642 return self._value 643 else: --> 644 raise self._value 645 646 def _set(self, i, obj):

~/bin/miniconda3/envs/scenic/lib/python3.6/multiprocessing/pool.py in worker() 117 job, i, func, args, kwds = task 118 try: --> 119 result = (True, func(*args, **kwds)) 120 except Exception as e: 121 if wrap_exception and func is not _helper_reraises_exception:

~/bin/miniconda3/envs/scenic/lib/python3.6/multiprocessing/pool.py in starmapstar() 45 46 def starmapstar(args): ---> 47 return list(itertools.starmap(args[0], args[1])) 48 49 #

~/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/binarization.py in derive_threshold() 50 return gmm2.bic(X) <= gmm1.bic(X) 51 ---> 52 if not isbimodal(data, method): 53 # For a unimodal distribution the threshold is set as mean plus two standard deviations. 54 return data.mean() + 2.0*data.std()

~/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/binarization.py in isbimodal() 41 if method == 'hdt': 42 # Use Hartigan's dip statistic to decide if distribution deviates from unimodality. ---> 43 , pval, = diptst(np.msort(data)) 44 return (pval is not None) and (pval <= 0.05) 45 else:

~/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/diptest.py in diptst() 51 """ diptest with pval """ 52 # sample dip ---> 53 d, (, idxs, left, , right, _) = dip_fn(dat, is_hist) 54 55 # simulate from null uniform

~/bin/miniconda3/envs/scenic/lib/python3.6/site-packages/pyscenic/diptest.py in dip_fn() 113 d = d_right 114 else: --> 115 xl = left_touchpoints[d_left == left_diffs][0] 116 xr = right_touchpoints[right_touchpoints >= xl][0] 117 d = d_left

IndexError: index 0 is out of bounds for axis 0 with size 0

wangjiawen2013 commented 4 years ago

My dataset have 59847 rows and 401 columns, I can run binarization with top 10000 rows and 401 columns successfully, and top 40000 rows or random sample 40000 rows from the dataset worked well too. However, the full dataset failed to run. May be there are some potential bugs in diptest, or this is only a occasional case specific to the dataset. Finally, I sampled 40000 cells randomly to generated the thresholds, then used (auc_mtx > thresholds).astype(int) to get the binarized auc matrix for all 59847 cells approximately

wangjiawen2013 commented 4 years ago

Besides, theere are no good documentation on binarization, for example, the "seed" parameter in derive_threshold() are not well explained, and some tutorials are needed.

wangjiawen2013 commented 4 years ago

here is my conda env:

Name Version Build Channel

_libgcc_mutex 0.1 conda_forge conda-forge _openmp_mutex 4.5 0_gnu conda-forge arboreto 0.1.5 pypi_0 pypi attrs 19.3.0 pypi_0 pypi backcall 0.1.0 py_0 conda-forge bokeh 2.0.1 py36h9f0ad1d_0 conda-forge boltons 20.1.0 pypi_0 pypi ca-certificates 2020.4.5.1 hecc5488_0 conda-forge certifi 2020.4.5.1 py36h9f0ad1d_0 conda-forge click 7.1.2 pyh9f0ad1d_0 conda-forge cloudpickle 1.4.1 py_0 conda-forge cycler 0.10.0 py_2 conda-forge cytoolz 0.10.1 py36h516909a_0 conda-forge dask 1.0.0 py_1 conda-forge dask-core 1.0.0 py_0 conda-forge dbus 1.13.6 he372182_0 conda-forge decorator 4.4.2 py_0 conda-forge dill 0.3.1.1 pypi_0 pypi distributed 1.28.1 py36_0 conda-forge expat 2.2.9 he1b5a44_2 conda-forge fontconfig 2.13.1 h86ecdb6_1001 conda-forge freetype 2.10.2 he06d7ca_0 conda-forge frozendict 1.2 pypi_0 pypi gettext 0.19.8.1 hc5be6a0_1002 conda-forge glib 2.64.2 h6f030ca_1 conda-forge gst-plugins-base 1.14.5 h0935bb2_2 conda-forge gstreamer 1.14.5 h36ae1b5_2 conda-forge h5py 2.10.0 pypi_0 pypi heapdict 1.0.1 py_0 conda-forge icu 64.2 he1b5a44_1 conda-forge interlap 0.2.6 pypi_0 pypi ipython 7.14.0 py36h9f0ad1d_0 conda-forge ipython_genutils 0.2.0 py_1 conda-forge jedi 0.17.0 py36h9f0ad1d_0 conda-forge jinja2 2.11.2 pyh9f0ad1d_0 conda-forge joblib 0.15.1 pypi_0 pypi jpeg 9c h14c3975_1001 conda-forge kiwisolver 1.2.0 py36hdb11119_0 conda-forge ld_impl_linux-64 2.34 h53a641e_0 conda-forge libblas 3.8.0 14_openblas conda-forge libcblas 3.8.0 14_openblas conda-forge libclang 9.0.1 default_hde54327_0 conda-forge libffi 3.2.1 he1b5a44_1007 conda-forge libgcc-ng 9.2.0 h24d8f2e_2 conda-forge libgfortran-ng 7.5.0 hdf63c60_6 conda-forge libgomp 9.2.0 h24d8f2e_2 conda-forge libiconv 1.15 h516909a_1006 conda-forge liblapack 3.8.0 14_openblas conda-forge libllvm9 9.0.1 he513fc3_1 conda-forge libopenblas 0.3.7 h5ec1e0e_6 conda-forge libpng 1.6.37 hed695b0_1 conda-forge libstdcxx-ng 9.2.0 hdf63c60_2 conda-forge libtiff 4.1.0 hc7e4089_6 conda-forge libuuid 2.32.1 h14c3975_1000 conda-forge libwebp-base 1.1.0 h516909a_3 conda-forge libxcb 1.13 h14c3975_1002 conda-forge libxkbcommon 0.10.0 he1b5a44_0 conda-forge libxml2 2.9.10 hee79883_0 conda-forge llvmlite 0.32.1 pypi_0 pypi locket 0.2.0 py_2 conda-forge loompy 3.0.6 pypi_0 pypi lz4-c 1.9.2 he1b5a44_1 conda-forge markupsafe 1.1.1 py36h8c4c3a4_1 conda-forge matplotlib 3.2.1 0 conda-forge matplotlib-base 3.2.1 py36hb8e4980_0 conda-forge msgpack-python 0.6.2 py36hc9558a2_0 conda-forge multiprocessing-on-dill 3.5.0a4 pypi_0 pypi ncurses 6.1 hf484d3e_1002 conda-forge networkx 2.4 pypi_0 pypi nspr 4.25 he1b5a44_0 conda-forge nss 3.47 he751ad9_0 conda-forge numba 0.49.1 pypi_0 pypi numpy 1.18.4 py36h7314795_0 conda-forge numpy-groupies 0+unknown pypi_0 pypi olefile 0.46 py_0 conda-forge openssl 1.1.1g h516909a_0 conda-forge packaging 20.1 py_0 conda-forge pandas 0.25.3 py36hb3f55d8_0 conda-forge parso 0.7.0 pyh9f0ad1d_0 conda-forge partd 1.1.0 py_0 conda-forge patsy 0.5.1 py_0 conda-forge pcre 8.44 he1b5a44_0 conda-forge pexpect 4.8.0 py36h9f0ad1d_1 conda-forge pickleshare 0.7.5 py36h9f0ad1d_1001 conda-forge pillow 7.1.2 py36h8328e55_0 conda-forge pip 20.1 pyh9f0ad1d_0 conda-forge prompt-toolkit 3.0.5 py_0 conda-forge psutil 5.7.0 py36h8c4c3a4_1 conda-forge pthread-stubs 0.4 h14c3975_1001 conda-forge ptyprocess 0.6.0 py_1001 conda-forge pyarrow 0.16.0 pypi_0 pypi pygments 2.6.1 py_0 conda-forge pyparsing 2.4.7 pyh9f0ad1d_0 conda-forge pyqt 5.12.3 py36haa643ae_3 conda-forge pyqt5-sip 4.19.18 pypi_0 pypi pyqtchart 5.12 pypi_0 pypi pyqtwebengine 5.12.1 pypi_0 pypi pyscenic 0.10.1 pypi_0 pypi python 3.6.10 h8356626_1011_cpython conda-forge python-dateutil 2.8.1 py_0 conda-forge python_abi 3.6 1_cp36m conda-forge pytz 2020.1 pyh9f0ad1d_0 conda-forge pyyaml 5.3.1 py36h8c4c3a4_0 conda-forge qt 5.12.5 hd8c4c69_1 conda-forge readline 8.0 hf8c457e_0 conda-forge scikit-learn 0.23.0 pypi_0 pypi scipy 1.4.1 py36h2d22cac_3 conda-forge seaborn 0.10.1 py_0 conda-forge setuptools 46.3.1 py36h9f0ad1d_0 conda-forge six 1.14.0 py_1 conda-forge sortedcontainers 2.1.0 py_0 conda-forge sqlite 3.30.1 hcee41ef_0 conda-forge statsmodels 0.11.1 py36h8c4c3a4_1 conda-forge tbb 2020.0.133 pypi_0 pypi tblib 1.6.0 py_0 conda-forge threadpoolctl 2.0.0 pypi_0 pypi tk 8.6.10 hed695b0_0 conda-forge toolz 0.10.0 py_0 conda-forge tornado 6.0.4 py36h8c4c3a4_1 conda-forge tqdm 4.46.0 pypi_0 pypi traitlets 4.3.3 py36h9f0ad1d_1 conda-forge typing_extensions 3.7.4.2 py_0 conda-forge umap-learn 0.4.3 pypi_0 pypi wcwidth 0.1.9 pyh9f0ad1d_0 conda-forge wheel 0.34.2 py_1 conda-forge xorg-libxau 1.0.9 h14c3975_0 conda-forge xorg-libxdmcp 1.1.3 h516909a_0 conda-forge xz 5.2.5 h516909a_0 conda-forge yaml 0.2.4 h516909a_0 conda-forge zict 2.0.0 py_0 conda-forge zlib 1.2.11 h516909a_1006 conda-forge zstd 1.4.4 h6597ccf_3 conda-forge

wangjiawen2013 commented 4 years ago

And pip list: Package Version

arboreto 0.1.5 attrs 19.3.0 backcall 0.1.0 bokeh 2.0.1 boltons 20.1.0 certifi 2020.4.5.1 click 7.1.2 cloudpickle 1.4.1 cycler 0.10.0 cytoolz 0.10.1 dask 1.0.0 decorator 4.4.2 dill 0.3.1.1 distributed 1.28.1 fbpca 1.0 frozendict 1.2 h5py 2.10.0 HeapDict 1.0.1 interlap 0.2.6 intervaltree 2.1.0 ipython 7.14.0 ipython-genutils 0.2.0 jedi 0.17.0 Jinja2 2.11.2 joblib 0.15.1 kiwisolver 1.2.0 llvmlite 0.32.1 locket 0.2.0 loompy 3.0.6 MarkupSafe 1.1.1 matplotlib 3.2.1 msgpack 0.6.2 multiprocessing-on-dill 3.5.0a4 networkx 2.4 numba 0.49.1 numpy 1.18.4 numpy-groupies 0+unknown olefile 0.46 packaging 20.1 pandas 0.25.3 parso 0.7.0 partd 1.1.0 patsy 0.5.1 pexpect 4.8.0 pickleshare 0.7.5 Pillow 7.1.2 pip 20.1 prompt-toolkit 3.0.5 psutil 5.7.0 ptyprocess 0.6.0 pyarrow 0.16.0 Pygments 2.6.1 pyparsing 2.4.7 PyQt5 5.12.3 PyQt5-sip 4.19.18 PyQtChart 5.12 PyQtWebEngine 5.12.1 pyscenic 0.10.1 python-dateutil 2.8.1 pytz 2020.1 PyYAML 5.3.1 scikit-learn 0.23.0 scipy 1.4.1 seaborn 0.10.1 setuptools 46.3.1.post20200515 six 1.14.0 sortedcontainers 2.1.0 statsmodels 0.11.1 tbb 2020.0.133 tblib 1.6.0 threadpoolctl 2.0.0 toolz 0.10.0 tornado 6.0.4 tqdm 4.46.0 traitlets 4.3.3 typing-extensions 3.7.4.2 umap-learn 0.4.3 wcwidth 0.1.9 wheel 0.34.2 zict 2.0.0

wangjiawen2013 commented 4 years ago

If num_workers was set default, binarize can run, while failed when set to more than 1. I am working on Centos 7 server.

cflerin commented 4 years ago

Hi @wangjiawen2013 ,

I haven't seen this error before. Does it work when setting num_workers to 1 for the full dataset?

I wonder if it could be something specific to your dataset. Have you done any cell/gene filtering? Is it possible that there are cells or genes with zero expression across the whole matrix?

wangjiawen2013 commented 4 years ago

No, it didn't work for the full dataset, no matter what num_workers is. I agree with you, this could be something specific to my dataset. I filtered cells and genes , and there are no cells and genes with zero expression. Perhaps it will never happen again for other datasets.

cflerin commented 4 years ago

Can you post the output of:

np.sum(auc_mtx, axis=0).sort_values()

and

np.sum(auc_mtx, axis=1).sort_values()

For comparison, I've successfully run the binarization step with 100s of thousands of cells with no problems (although it does use a lot of memory).

wangjiawen2013 commented 4 years ago

It took several hours to run 40000 cells and 401 regulons (num_workers=1). Will the program be interrupted when two cells have the same expression pattern or some regulons have the same distribution across cells ?

wangjiawen2013 commented 4 years ago

np.sum(auc_mtx, axis=0).sort_values():

SRY(+) | 12.34071 ZIC5(+) | 45.86162 NEUROG1(+) | 70.31318 KLF14(+) | 87.62278 CUX2(+) | 93.68485 ZNF302(+) | 114.3998 PKNOX2(+) | 121.1167 FOXA2(+) | 129.1522 BHLHE22(+) | 147.101 DLX6(+) | 168.3566 NKX6-2(+) | 170.608 HSFY2(+) | 175.1683 HOXA3(+) | 232.2691 SOX21(+) | 233.4766 POU3F2(+) | 238.2204 ZNF768(+) | 245.4822 ATF6B(+) | 249.4129 NKX6-3(+) | 252.8736 MSX2(+) | 275.4885 ZNF577(+) | 305.9244 HOXD3(+) | 320.7293 LMX1B(+) | 355.3763 NR1I3(+) | 363.1422 HOXD10(+) | 365.188 RXRA(+) | 383.9088 POU6F1(+) | 488.0977 MYPOP(+) | 529.4444 ZNF281(+) | 536.2759 EMX2(+) | 545.0218 OTP(+) | 548.5062 HOXC9(+) | 558.2894 SMAD3(+) | 563.0815 ZFP90(+) | 588.0395 SOX8(+) | 596.8208 HOXA7(+) | 600.9124 SOX2(+) | 616.2755 ZNF605(+) | 619.1444 NR2F1(+) | 640.8229 SOX10(+) | 642.2514 HES5(+) | 653.1455 ATF2(+) | 653.8665 ZNF658(+) | 663.552 HOXD9(+) | 663.7787 HOXC6(+) | 674.6403 HOXB6(+) | 693.6248 MNT(+) | 710.7614 CIC(+) | 725.0597 RORC(+) | 725.1168 PBX1(+) | 740.2425 RARB(+) | 748.1697 ZNF471(+) | 759.2646 HOXC8(+) | 785.8955 NFATC3(+) | 786.6783 HOXD1(+) | 832.5244 GATA1(+) | 838.0994 NFE2(+) | 845.9267 NR2C1(+) | 882.7984 DEAF1(+) | 895.9629 MECOM(+) | 897.5863 VEZF1(+) | 905.6005 PPARA(+) | 914.7326 ZFP1(+) | 915.237 MXD4(+) | 930.506 MGA(+) | 931.918 ZNF746(+) | 942.2529 KLF11(+) | 944.3009 PRRX2(+) | 958.0921 PATZ1(+) | 962.8519 LHX9(+) | 972.1979 HESX1(+) | 1025.846 NR2F6(+) | 1032.104 HOXD11(+) | 1057.519 ZNF48(+) | 1059.075 MZF1(+) | 1100.852 ZIC1(+) | 1107.421 ZBTB7B(+) | 1123.188 SOX11(+) | 1137.478 ZNF740(+) | 1153.313 TFAP4(+) | 1180.058 TCF7(+) | 1180.396 HOXA2(+) | 1188.649 RXRB(+) | 1211.613 SMAD1(+) | 1226.662 SOX18(+) | 1267.201 BCL6B(+) | 1276.373 POU1F1(+) | 1288.005 NR6A1(+) | 1298.689 HOXB7(+) | 1341.541 MLXIP(+) | 1368.588 TLX1(+) | 1371.56 HOXA10(+) | 1373.19 TP73(+) | 1386.458 ZNF568(+) | 1396.57 TCF7L1(+) | 1426.859 PPARD(+) | 1447.833 TFDP2(+) | 1464.361 ZNF516(+) | 1470.786 ELK4(+) | 1504.033 CREBL2(+) | 1521.356 PGR(+) | 1521.69 MYB(+) | 1534.124 GATA5(+) | 1543.824 E2F3(+) | 1551.863 ZFP64(+) | 1554.739 SOX13(+) | 1563.672 NEUROG3(+) | 1614.468 ZNF362(+) | 1623.706 TEF(+) | 1627.305 SP1(+) | 1657.611 PAX5(+) | 1676.42 SETDB1(+) | 1713.202 CLOCK(+) | 1739.537 SOX7(+) | 1755.179 LYL1(+) | 1761.088 VDR(+) | 1761.833 FOXN1(+) | 1777.841 E2F6(+) | 1800.28 FOXN2(+) | 1805.531 SHOX2(+) | 1848.543 FOXK1(+) | 1876.854 HOXA5(+) | 1887.292 MEF2D(+) | 1887.824 GTF2IRD1(+) | 1894.581 GATA4(+) | 1895.192 CEBPZ(+) | 1920.595 ISL1(+) | 1934.152 GATA2(+) | 1953.129 KLF8(+) | 1954.389 RFX5(+) | 1956.885 THAP11(+) | 1965.183 HEY1(+) | 1966.461 MEIS1(+) | 1966.89 KLF16(+) | 2002.179 ZBTB7A(+) | 2007.223 CTCF(+) | 2027.234 SP2(+) | 2034.702 FOXL1(+) | 2039.262 NR1H3(+) | 2056.085 HOXA9(+) | 2078.037 FOXJ3(+) | 2087.617 EOMES(+) | 2156.577 T(+) | 2166.332 FOXN3(+) | 2168.147 FOXO1(+) | 2170.557 ZNF143(+) | 2182.869 ZNF250(+) | 2202.808 STAT6(+) | 2208.022 ZNF274(+) | 2217.088 NR1D2(+) | 2245.602 TFCP2(+) | 2246.678 SOX17(+) | 2246.866 SIX3(+) | 2251.859 KLF1(+) | 2252.562 TBX1(+) | 2256.483 RXRG(+) | 2268.295 PBX3(+) | 2269.903 TFAP2A(+) | 2310.046 ZNF569(+) | 2322.319 NR2C2(+) | 2325.993 CREB1(+) | 2339.487 CEBPG(+) | 2345.383 HLF(+) | 2363.906 FOXF1(+) | 2364.763 SOX9(+) | 2364.973 KDM5B(+) | 2398.09 VAX2(+) | 2409.006 ELF4(+) | 2409.016 SIX5(+) | 2412.478 JDP2(+) | 2470.795 IRF2(+) | 2480.604 ZNF76(+) | 2531.21 OLIG1(+) | 2533.715 KLF9(+) | 2543.128 GFI1(+) | 2568.032 NR1D1(+) | 2575.266 GBX1(+) | 2585.541 RFX1(+) | 2592.207 POU3F1(+) | 2631.142 GLI1(+) | 2638.277 RARA(+) | 2659.474 JUND(+) | 2677.523 ARNT(+) | 2698.744 TBX2(+) | 2700.992 FOXL2(+) | 2702.008 MYBL1(+) | 2706.28 SPIC(+) | 2746.357 KLF12(+) | 2763.05 ATF1(+) | 2786.704 MAFK(+) | 2793.21 ZIC4(+) | 2812.191 E2F2(+) | 2818.49 DBP(+) | 2826.448 HINFP(+) | 2831.727 ELF2(+) | 2836.868 TFEB(+) | 2841.871 E2F8(+) | 2850.81 ZFY(+) | 2878.679 STAT2(+) | 2894.998 SP3(+) | 2896.092 ELK3(+) | 2899.194 ZNF680(+) | 2905.019 BRF2(+) | 2911.527 PPARG(+) | 2919.174 TCF12(+) | 2931.235 GABPA(+) | 2942.645 SETBP1(+) | 2962.109 HOXA6(+) | 2980.017 ZNF91(+) | 2989.153 GATA3(+) | 3002.839 TBX21(+) | 3023.608 ZNF502(+) | 3026.903 BATF(+) | 3036.242 ZNF354B(+) | 3062.293 SIX4(+) | 3074.069 KLF2(+) | 3093.305 RUNX3(+) | 3131.734 ZNF284(+) | 3161.727 ZNF23(+) | 3161.905 REST(+) | 3162.315 ETV5(+) | 3192.456 USF1(+) | 3207.169 IRF5(+) | 3207.915 HIC1(+) | 3218.266 GLIS2(+) | 3219.107 SREBF2(+) | 3221.748 RUNX2(+) | 3251.98 WT1(+) | 3253.505 ZNF254(+) | 3264.203 RFX2(+) | 3279.211 NFYB(+) | 3279.299 HIVEP3(+) | 3283.56 NFYA(+) | 3288.713 NFYC(+) | 3294.191 KLF7(+) | 3300.192 STAT1(+) | 3303.049 GATA6(+) | 3329.976 NFATC4(+) | 3353.353 TRPS1(+) | 3361.164 TFCP2L1(+) | 3389.718 KLF13(+) | 3435.941 TEAD1(+) | 3491.768 TBP(+) | 3499.368 MEF2B(+) | 3554.533 IKZF1(+) | 3556.223 E2F1(+) | 3578.852 RBPJ(+) | 3581.017 RUNX1(+) | 3614.236 TBX5(+) | 3615.258 FOXO4(+) | 3626.073 USF2(+) | 3645.442 KLF3(+) | 3645.573 GLI2(+) | 3667.539 ALX1(+) | 3669.777 HIVEP1(+) | 3672.564 FLI1(+) | 3687.006 TFE3(+) | 3706.423 MAF(+) | 3720.33 YBX1(+) | 3732.038 ZNF597(+) | 3804.809 ELF1(+) | 3804.988 ZNF148(+) | 3818.276 ETS1(+) | 3839.075 BCL6(+) | 3841.861 TBX15(+) | 3864.038 CREB3L2(+) | 3903.628 CREB3L1(+) | 3911.438 TFEC(+) | 3920.081 ZNF189(+) | 3930.29 ETV3(+) | 3954.537 FOXC2(+) | 3956.751 MITF(+) | 3964.059 STAT5A(+) | 3975.079 TEAD2(+) | 4005.024 SPI1(+) | 4037.984 THAP1(+) | 4048.62 IRF3(+) | 4085.447 EGR4(+) | 4091.59 BACH1(+) | 4091.763 GMEB1(+) | 4109.809 IRF8(+) | 4118.857 PRDM1(+) | 4122.357 TWIST1(+) | 4127.185 IRF9(+) | 4194.815 TBX20(+) | 4201.922 REL(+) | 4209.016 MAZ(+) | 4242.032 MSX1(+) | 4244.862 CEBPE(+) | 4273.872 FOXP2(+) | 4276.078 HIVEP2(+) | 4278.584 ZNF398(+) | 4296.374 RARG(+) | 4306.508 LEF1(+) | 4331.653 NRF1(+) | 4332.303 E2F7(+) | 4348.817 EGR1(+) | 4354.374 POU2F1(+) | 4366.101 BHLHE41(+) | 4378.632 ELF5(+) | 4386.305 EGR2(+) | 4421.862 FOS(+) | 4442.377 NFKB2(+) | 4511.846 SOX4(+) | 4526.192 IRF4(+) | 4563.289 NFKB1(+) | 4580.035 ZBTB41(+) | 4585.923 HOXB5(+) | 4588.33 MYC(+) | 4589.097 ERG(+) | 4618.471 BHLHE40(+) | 4619.495 ZNF467(+) | 4662.942 FOXC1(+) | 4692.491 ETS2(+) | 4708.851 TEAD3(+) | 4748.293 GBX2(+) | 4748.883 ZNF264(+) | 4774.975 SREBF1(+) | 4777.982 MEF2A(+) | 4784.573 LHX8(+) | 4800.355 NFE2L3(+) | 4808.534 DBX2(+) | 4830.416 YY1(+) | 4861.617 RELB(+) | 4880.234 TP53(+) | 4922.82 NFATC1(+) | 4944.212 ATF6(+) | 4947.119 ESRRA(+) | 5000.332 NR2F2(+) | 5017.133 ATF5(+) | 5034.89 SP4(+) | 5046.925 MAX(+) | 5049.23 IRF7(+) | 5064.14 TCF21(+) | 5127.768 TWIST2(+) | 5167.152 NR3C1(+) | 5188.891 CREB5(+) | 5212.122 EGR3(+) | 5220.825 FOXO3(+) | 5266.168 SOX15(+) | 5297.587 CEBPD(+) | 5305.113 E2F4(+) | 5349.029 SRF(+) | 5383.037 FOSL2(+) | 5503.508 MXI1(+) | 5511.595 CEBPB(+) | 5541.465 KLF6(+) | 5555.881 ETV4(+) | 5557.517 ETV6(+) | 5592.27 KLF15(+) | 5631.333 IRF1(+) | 5651.169 HIF1A(+) | 5698.256 KLF5(+) | 5754.738 NKX3-2(+) | 5798.433 MLX(+) | 5877.482 EHF(+) | 5986.126 CREM(+) | 6000.951 TGIF2(+) | 6105.157 CEBPA(+) | 6109.947 ZNF8(+) | 6114.394 NKX2-5(+) | 6140.55 TFDP1(+) | 6262.285 NFE2L2(+) | 6273.781 NFIL3(+) | 6379.031 MAFG(+) | 6420.883 FOXK2(+) | 6450.167 EBF1(+) | 6520.971 CUX1(+) | 6536.586 FOSL1(+) | 6546.942 TGIF1(+) | 6597.321 PRDM16(+) | 6619.049 XBP1(+) | 6759.762 KLF4(+) | 6801.321 MEF2C(+) | 6934.994 RELA(+) | 7008.795 SIX1(+) | 7085.599 BACH2(+) | 7380.565 BHLHA15(+) | 7699.024 BCL11A(+) | 7701.576 ATF3(+) | 8328.852 ATF4(+) | 8348.838 HSF1(+) | 8371.437 MAFB(+) | 8395.089 SPIB(+) | 8530.756 ZNF358(+) | 8602.69 ETV7(+) | 8855.706 JUNB(+) | 9077.182 AR(+) | 10003.77 GMEB2(+) | 10374.73 MAFF(+) | 10416.05 ELK1(+) | 11932.89 GTF2B(+) | 12644.69 FOSB(+) | 13132.04 CREB3(+) | 15508.6 NFIA(+) | 20265.22 KLF10(+) | 20506.35 STAT3(+) | 22590.98 NFIC(+) | 22977.15 CREB3L4(+) | 23218.69 ZNF222(+) | 23970.96 JUN(+) | 24631.8 ETV2(+) | 29242.83

wangjiawen2013 commented 4 years ago

np.sum(auc_mtx, axis=1).sort_values(): (it's too long to show the full dataset)

cflerin commented 4 years ago

Thanks for posting that output. Unfortunately it doesn't look like there's anything obviously wrong with your data (missing values, etc.), so I'm not sure why you are getting this index error.

akramdi commented 4 years ago

Hi,

I am here because I've encountered the exact same error with my data. In my case, setting num_workers to 1 did not work either (it results in the same error).. Is there any other solution to binarize the auc matrix ?

andrewni4313 commented 3 years ago

I have the same issue. It occurs on some datasets but not on others and size doesn't seem to be the issue as smaller datasets sometimes error out when larger ones don't.

andrewni4313 commented 3 years ago

Update so I ran binarization individually on every single regulon/column of auc_mtx and it errored out on a specific column:

for start in tqdm(range(0, len(auc_mtx.columns))):
    print(auc_mtx.iloc[:,start:start+1].columns)
    binary_mtx, auc_thresholds = binarize(auc_mtx.iloc[:,start:start+1], seed=42, num_workers=1)

56% 58/103 [09:32<08:00, 10.69s/it]
Index(['AI987944'], dtype='object', name='Regulon')
Index(['Arntl2'], dtype='object', name='Regulon')
Index(['Atf3'], dtype='object', name='Regulon')
Index(['Atf4'], dtype='object', name='Regulon')
Index(['Bcl11b'], dtype='object', name='Regulon')
Index(['Bcl3'], dtype='object', name='Regulon')
Index(['Bclaf1'], dtype='object', name='Regulon')
Index(['Bhlhe41'], dtype='object', name='Regulon')
Index(['Bmyc'], dtype='object', name='Regulon')
Index(['Brf1'], dtype='object', name='Regulon')
Index(['Cebpb'], dtype='object', name='Regulon')
Index(['Cebpz'], dtype='object', name='Regulon')
Index(['Chd2'], dtype='object', name='Regulon')
Index(['Creb1'], dtype='object', name='Regulon')
Index(['Creb5'], dtype='object', name='Regulon')
Index(['Crem'], dtype='object', name='Regulon')
Index(['Ddit3'], dtype='object', name='Regulon')
Index(['E2f1'], dtype='object', name='Regulon')
Index(['E2f2'], dtype='object', name='Regulon')
Index(['E2f4'], dtype='object', name='Regulon')
Index(['E2f7'], dtype='object', name='Regulon')
Index(['E2f8'], dtype='object', name='Regulon')
Index(['Egr2'], dtype='object', name='Regulon')
Index(['Ehf'], dtype='object', name='Regulon')
Index(['Elf1'], dtype='object', name='Regulon')
Index(['Elf2'], dtype='object', name='Regulon')
Index(['Elf4'], dtype='object', name='Regulon')
Index(['Ets1'], dtype='object', name='Regulon')
Index(['Etv1'], dtype='object', name='Regulon')
Index(['Etv6'], dtype='object', name='Regulon')
Index(['Ezh2'], dtype='object', name='Regulon')
Index(['Fli1'], dtype='object', name='Regulon')
Index(['Fos'], dtype='object', name='Regulon')
Index(['Fosb'], dtype='object', name='Regulon')
Index(['Fosl2'], dtype='object', name='Regulon')
Index(['Foxn3'], dtype='object', name='Regulon')
Index(['Foxo1'], dtype='object', name='Regulon')
Index(['Foxo3'], dtype='object', name='Regulon')
Index(['Foxp2'], dtype='object', name='Regulon')
Index(['Gabpa'], dtype='object', name='Regulon')
Index(['Gabpb1'], dtype='object', name='Regulon')
Index(['Gm14308'], dtype='object', name='Regulon')
Index(['Gtf2f1'], dtype='object', name='Regulon')
Index(['Gtf3c2'], dtype='object', name='Regulon')
Index(['Hdac2'], dtype='object', name='Regulon')
Index(['Hmga1'], dtype='object', name='Regulon')
Index(['Hoxb2'], dtype='object', name='Regulon')
Index(['Hsf1'], dtype='object', name='Regulon')
Index(['Irf1'], dtype='object', name='Regulon')
Index(['Irf2'], dtype='object', name='Regulon')
Index(['Irf7'], dtype='object', name='Regulon')
Index(['Irf8'], dtype='object', name='Regulon')
Index(['Irf9'], dtype='object', name='Regulon')
Index(['Jdp2'], dtype='object', name='Regulon')
Index(['Jun'], dtype='object', name='Regulon')
Index(['Junb'], dtype='object', name='Regulon')
Index(['Jund'], dtype='object', name='Regulon')
Index(['Kdm5a'], dtype='object', name='Regulon')
Index(['Klf2'], dtype='object', name='Regulon')

IndexError: index 0 is out of bounds for axis 0 with size 0

It errors out on Klf2 regulon. This doesn't seem to be specific to this regulon, however, as Klf2 works in other datasets I have tried.

When I removed this column/regulon the binarization proceeded as normal and did not error.

Here is the data for that one column: error_column.txt

Hopefully, a solution to the error can be found with this!

akramdi commented 3 years ago

Hi @andrewni4313 Here's my fix for this issue. I found that this error is caused by a failure of one of the two methods used to drive the threshold for binarization. Here's the guilty function:

https://github.com/aertslab/pySCENIC/blob/436561fd6179ececdf8ae6faa6da3dc4963a11e6/src/pyscenic/binarization.py#L18-L64

Two methods are implemented (HDT, BIC) but only one is used (method: str = 'hdt'). I edited this function so that the second method is used if the first one fails (done by catching the IndexError):

def derive_threshold(auc_mtx: pd.DataFrame, regulon_name: str, seed=None) -> float:
    '''
    Derive threshold on the AUC values of the given regulon to binarize the cells in two clusters: "on" versus "off"
    state of the regulator.

    :param auc_mtx: The dataframe with the AUC values for all cells and regulons (n_cells x n_regulons).
    :param regulon_name: the name of the regulon for which to predict the threshold.
    :param method: The method to use to decide if the distribution of AUC values for the given regulon is not unimodel.
        Can be either Hartigan's Dip Test (HDT) or Bayesian Information Content (BIC). The former method performs better
        but takes considerable more time to execute (40min for 350 regulons). The BIC compares the BIC for two Gaussian
        Mixture Models: single versus two components.
    :return: The threshold on the AUC values.
    '''
    assert auc_mtx is not None and not auc_mtx.empty
    assert regulon_name in auc_mtx.columns
    #assert method in {'hdt', 'bic'}

    data = auc_mtx[regulon_name].values
    if seed:
        np.random.seed(seed=seed)

    def isbimodal(data):
        try:
            # Use Hartigan's dip statistic to decide if distribution deviates from unimodality.
            _, pval, _ = diptst(np.msort(data))
            return (pval is not None) and (pval <= 0.05)
        except IndexError:
            # Compare Bayesian Information Content of two Gaussian Mixture Models.
            X = data.reshape(-1, 1)
            gmm2 = mixture.GaussianMixture(n_components=2, covariance_type='full', random_state=seed).fit(X)
            gmm1 = mixture.GaussianMixture(n_components=1, covariance_type='full', random_state=seed).fit(X)
            return gmm2.bic(X) <= gmm1.bic(X)

    if not isbimodal(data):
        # For a unimodal distribution the threshold is set as mean plus two standard deviations.
        return data.mean() + 2.0*data.std()
    else:
        # Fit a two component Gaussian Mixture model on the AUC distribution using an Expectation-Maximization algorithm
        # to identify the peaks in the distribution.
        gmm2 = mixture.GaussianMixture(n_components=2, covariance_type='full', random_state=seed).fit(data.reshape(-1, 1))
        # For a bimodal distribution the threshold is defined as the "trough" in between the two peaks.
        # This is solved as a minimization problem on the kernel smoothed density.
        return minimize_scalar(fun=stats.gaussian_kde(data),
                               bounds=sorted(gmm2.means_),
                               method='bounded').x[0]

This worked for me, I hope it helps! Maybe @cflerin can confirm if this solves the issue.

Best, Amira

wangjiawen2013 commented 1 year ago

Hi, I met an new error in this step (auc_mtx_binary = binarize(auc_mtx,num_workers=15) now, hope you kindly enough to check it :

multiprocessing.pool.RemoteTraceback: """ Traceback (most recent call last): File "/home/wangjw/programs/miniconda3/envs/scenic/lib/python3.7/multiprocessing/pool.py", line 121, in worker result = (True, func(*args, **kwds)) File "/home/wangjw/programs/miniconda3/envs/scenic/lib/python3.7/multiprocessing/pool.py", line 47, in starmapstar return list(itertools.starmap(args[0], args[1])) File "/home/wangjw/programs/miniconda3/envs/scenic/lib/python3.7/site-packages/pyscenic/binarization.py", line 52, in derivethreshold if not isbimodal(data, method): File "/home/wangjw/programs/miniconda3/envs/scenic/lib/python3.7/site-packages/pyscenic/binarization.py", line 43, in isbimodal , pval, = diptst(np.msort(data)) File "/home/wangjw/programs/miniconda3/envs/scenic/lib/python3.7/site-packages/pyscenic/diptest.py", line 53, in diptst d, (, idxs, left, , right, ) = dip_fn(dat, is_hist) File "/home/wangjw/programs/miniconda3/envs/scenic/lib/python3.7/site-packages/pyscenic/diptest.py", line 102, in dip_fn left_part, left_touchpoints = gcm(work_cdf-work_histogram, work_idxs) File "/home/wangjw/programs/miniconda3/envs/scenic/lib/python3.7/site-packages/pyscenic/diptest.py", line 29, in gcm distances = work_idxs[1:] - work_idxs[0] numpy.core._exceptions._UFuncNoLoopError: ufunc 'subtract' did not contain a loop with signature matching types (dtype('<U25'), dtype('<U25')) -> None """

The above exception was the direct cause of the following exception:

Traceback (most recent call last): File "binarize.py", line 6, in auc_mtx_binary = binarize(auc_mtx,num_workers=15) File "/home/wangjw/programs/miniconda3/envs/scenic/lib/python3.7/site-packages/pyscenic/binarization.py", line 79, in binarize thresholds = derive_thresholds(auc_mtx) File "/home/wangjw/programs/miniconda3/envs/scenic/lib/python3.7/site-packages/pyscenic/binarization.py", line 77, in derive_thresholds thrs = p.starmap( derive_threshold, [ (auc_mtx, c, seed) for c in auc_mtx.columns ] ) File "/home/wangjw/programs/miniconda3/envs/scenic/lib/python3.7/multiprocessing/pool.py", line 276, in starmap return self._map_async(func, iterable, starmapstar, chunksize).get() File "/home/wangjw/programs/miniconda3/envs/scenic/lib/python3.7/multiprocessing/pool.py", line 657, in get raise self._value numpy.core._exceptions.UFuncTypeError: ufunc 'subtract' did not contain a loop with signature matching types (dtype('<U25'), dtype('<U25')) -> None

Tripfantasy commented 8 months ago

@wangjiawen2013 Hello, have you found a fix for this error? Encountered it recently so I'm going to leave an additional note:

I think a likely issue is when converting the pd.dataframe generated from the aucell command to a .csv file it will default the 'cell' index to a culumn. Specifying the cell column as an index in pd.read_csv() allowed the binarization function to run successfully. Hope this helps!

bin_auc = pd.read_csv('bin_auc.csv',index_col='Cell')

wangjiawen2013 commented 8 months ago

@Tripfantasy I haven't used it for a long time. I'll try it later.

apal6 commented 2 days ago

Hi,

I am having the same error- Can you help ?

bin_auc = pd.read_csv('auc_mtx.csv',index_col='Cell')
auc_mtx_bi = binarize(bin_auc)

Error:

`---------------------------------------------------------------------------
RemoteTraceback                           Traceback (most recent call last)
RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/share/software/user/open/python/3.9.0/lib/python3.9/multiprocessing/pool.py", line 125, in worker
    result = (True, func(*args, **kwds))
  File "/share/software/user/open/python/3.9.0/lib/python3.9/multiprocessing/pool.py", line 51, in starmapstar
    return list(itertools.starmap(args[0], args[1]))
  File "/home/users/apal6/.local/lib/python3.9/site-packages/pyscenic/binarization.py", line 56, in derive_threshold
    if not isbimodal(data, method):
  File "/home/users/apal6/.local/lib/python3.9/site-packages/pyscenic/binarization.py", line 43, in isbimodal
    _, pval, _ = diptst(np.msort(data))
  File "/home/users/apal6/.local/lib/python3.9/site-packages/pyscenic/diptest.py", line 64, in diptst
    else (np.less(d, unif_dips).sum() + 1) / (np.float(numt) + 1)
  File "/share/software/user/open/py-numpy/1.24.2_py39/lib/python3.9/site-packages/numpy/__init__.py", line 305, in __getattr__
    raise AttributeError(__former_attrs__[attr])
AttributeError: module 'numpy' has no attribute 'float'.
`np.float` was a deprecated alias for the builtin `float`. To avoid this error in existing code, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:
    https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
"""

The above exception was the direct cause of the following exception:

AttributeError                            Traceback (most recent call last)
/tmp/ipykernel_5764/2118852679.py in <module>
----> 1 auc_mtx_bi = binarize(bin_auc)

~/.local/lib/python3.9/site-packages/pyscenic/binarization.py in binarize(auc_mtx, threshold_overides, seed, num_workers)
     92         return pd.Series(index=auc_mtx.columns, data=thrs)
     93 
---> 94     thresholds = derive_thresholds(auc_mtx)
     95     if threshold_overides is not None:
     96         thresholds[list(threshold_overides.keys())] = list(threshold_overides.values())

~/.local/lib/python3.9/site-packages/pyscenic/binarization.py in derive_thresholds(auc_mtx, seed)
     87     def derive_thresholds(auc_mtx, seed=seed):
     88         with Pool(processes=num_workers) as p:
---> 89             thrs = p.starmap(
     90                 derive_threshold, [(auc_mtx, c, seed) for c in auc_mtx.columns]
     91             )

/share/software/user/open/python/3.9.0/lib/python3.9/multiprocessing/pool.py in starmap(self, func, iterable, chunksize)
    370         `func` and (a, b) becomes func(a, b).
    371         '''
--> 372         return self._map_async(func, iterable, starmapstar, chunksize).get()
    373 
    374     def starmap_async(self, func, iterable, chunksize=None, callback=None,

/share/software/user/open/python/3.9.0/lib/python3.9/multiprocessing/pool.py in get(self, timeout)
    769             return self._value
    770         else:
--> 771             raise self._value
    772 
    773     def _set(self, i, obj):

/share/software/user/open/python/3.9.0/lib/python3.9/multiprocessing/pool.py in worker()
    123         job, i, func, args, kwds = task
    124         try:
--> 125             result = (True, func(*args, **kwds))
    126         except Exception as e:
    127             if wrap_exception and func is not _helper_reraises_exception:

/share/software/user/open/python/3.9.0/lib/python3.9/multiprocessing/pool.py in starmapstar()
     49 
     50 def starmapstar(args):
---> 51     return list(itertools.starmap(args[0], args[1]))
     52 
     53 #

~/.local/lib/python3.9/site-packages/pyscenic/binarization.py in derive_threshold()
     54             return gmm2.bic(X) <= gmm1.bic(X)
     55 
---> 56     if not isbimodal(data, method):
     57         # For a unimodal distribution the threshold is set as mean plus two standard deviations.
     58         return data.mean() + 2.0 * data.std()

~/.local/lib/python3.9/site-packages/pyscenic/binarization.py in isbimodal()
     41         if method == "hdt":
     42             # Use Hartigan's dip statistic to decide if distribution deviates from unimodality.
---> 43             _, pval, _ = diptst(np.msort(data))
     44             return (pval is not None) and (pval <= 0.05)
     45         else:

~/.local/lib/python3.9/site-packages/pyscenic/diptest.py in diptst()
     62         None
     63         if unif_dips.sum() == 0
---> 64         else (np.less(d, unif_dips).sum() + 1) / (np.float(numt) + 1)
     65     )
     66 

/share/software/user/open/py-numpy/1.24.2_py39/lib/python3.9/site-packages/numpy/__init__.py in __getattr__()
    303 
    304         if attr in __former_attrs__:
--> 305             raise AttributeError(__former_attrs__[attr])
    306 
    307         # Importing Tester requires importing all of UnitTest which is not a

AttributeError: module 'numpy' has no attribute 'float'.
`np.float` was a deprecated alias for the builtin `float`. To avoid this error in existing code, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:
    https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations`

Thank you, Aastha

Tripfantasy commented 2 days ago

AttributeError: module 'numpy' has no attribute 'float'. np.float was a deprecated alias for the builtin float. To avoid this error in existing code, use float by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use np.float64 here. The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations`

Hi @apal6 , this portion explains some dependency error between numpy and the version of scenic you're running. This doesn't seem to involve any issues with the input pandas dataframe.

You'll either want to downgrade the numpy version to something prior to 1.20.0, or modify the binarization function source code to use the builtin 'float()' rather than 'np.float()'. Depending on the scenic changelog, they may've fixed this in a newer version.

You could change source code in 'binarization.py' at this line and see if it fixes:

else (np.less(d, unif_dips).sum() + 1) / (np.float(numt) + 1)

to:

else (np.less(d, unif_dips).sum() + 1) / (float(numt) + 1)

Good luck!

apal6 commented 2 days ago

Thank you, I still don't see this particular code:

def binarize(
    auc_mtx: pd.DataFrame,
    threshold_overides: Optional[Mapping[str, float]] = None,
    seed=None,
    num_workers=1,
) -> (pd.DataFrame, pd.Series):
    """
    "Binarize" the supplied AUC matrix, i.e. decide if for each cells in the matrix a regulon is active or not based
    on the bimodal distribution of the AUC values for that regulon.

    :param auc_mtx: The dataframe with the AUC values for all cells and regulons (n_cells x n_regulons).
    :param threshold_overides: A dictionary that maps name of regulons to manually set thresholds.
    :return: A "binarized" dataframe and a series containing the AUC threshold used for each regulon.
    """

    def derive_thresholds(auc_mtx, seed=seed):
        with Pool(processes=num_workers) as p:
            thrs = p.starmap(
                derive_threshold, [(auc_mtx, c, seed) for c in auc_mtx.columns]
            )
        return pd.Series(index=auc_mtx.columns, data=thrs)

    thresholds = derive_thresholds(auc_mtx)
    if threshold_overides is not None:
        thresholds[list(threshold_overides.keys())] = list(threshold_overides.values())
    return (auc_mtx > thresholds).astype(int), thresholds

Tripfantasy commented 2 days ago

Thank you, I still don't see this particular code

My bad @apal6 , the code you want to change is actually in diptest.py:

def diptst(dat, is_hist=False, numt=1000):
    """diptest with pval"""
    # sample dip
    d, (_, idxs, left, _, right, _) = dip_fn(dat, is_hist)

    # simulate from null uniform
    unifs = np.random.uniform(size=numt * idxs.shape[0]).reshape([numt, idxs.shape[0]])
    unif_dips = np.apply_along_axis(dip_fn, 1, unifs, is_hist, True)

    # count dips greater or equal to d, add 1/1 to prevent a pvalue of 0
    pval = (
        None
        if unif_dips.sum() == 0
        else (np.less(d, unif_dips).sum() + 1) / (np.float(numt) + 1)
    )

aertslab / pySCENIC

Binarization error #174

Name Version Build Channel