Combat populates adata.X with NANs so sc.pp.highly_variable_genes function outputs error

[auesro]

After running sc.pp.combat(adata, key='sample'), adata.X is full of NaNs and sc.pp.highly_variagle_genes(adata) fails. No issues (and no NaNs) if NOT running the combat correction.

sc.pp.combat(adata, key='sample')
sc.pp.highly_variable_genes(adata)

In [1]: sc.pp.combat(adata, key='sample')
/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/anndata/_core/anndata.py:21: FutureWarning: pandas.core.index is deprecated and will be removed in a future version.  The public classes are available in the top-level namespace.
  from pandas.core.index import RangeIndex
/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/sklearn/externals/six.py:31: FutureWarning: The module is deprecated in version 0.21 and will be removed in version 0.23 since we've dropped support for Python 2.7. Please rely on the official version of six (https://pypi.org/project/six/).
  "(https://pypi.org/project/six/).", FutureWarning)
scanpy==1.4.6 anndata==0.7.1 umap==0.4.1 numpy==1.18.1 scipy==1.4.1 pandas==1.0.3 scikit-learn==0.22.2.post1 statsmodels==0.11.1 python-igraph==0.8.0
Standardizing Data across genes.

Found 11 batches

Found 0 numerical variables:
	

Found 3 genes with zero variance.
Fitting L/S model and finding priors

Finding parametric adjustments

Adjusting data

/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/scanpy/preprocessing/_combat.py:338: RuntimeWarning: invalid value encountered in true_divide
  change = max((abs(g_new - g_old) / g_old).max(), (abs(d_new - d_old) / d_old).max())
/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/scanpy/preprocessing/_combat.py:338: RuntimeWarning: divide by zero encountered in true_divide
  change = max((abs(g_new - g_old) / g_old).max(), (abs(d_new - d_old) / d_old).max())

In [2]: sc.pp.highly_variable_genes(adata)
extracting highly variable genes
Traceback (most recent call last):

  File "<ipython-input-2-7727f5f928cd>", line 1, in <module>
    sc.pp.highly_variable_genes(adata)

  File "/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/scanpy/preprocessing/_highly_variable_genes.py", line 235, in highly_variable_genes
    flavor=flavor,

  File "/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/scanpy/preprocessing/_highly_variable_genes.py", line 65, in _highly_variable_genes_single_batch
    df['mean_bin'] = pd.cut(df['means'], bins=n_bins)

  File "/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/pandas/core/reshape/tile.py", line 265, in cut
    duplicates=duplicates,

  File "/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/pandas/core/reshape/tile.py", line 381, in _bins_to_cuts
    f"Bin edges must be unique: {repr(bins)}.\n"

ValueError: Bin edges must be unique: array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan]).
You can drop duplicate edges by setting the 'duplicates' kwarg

Versions:

scanpy==1.4.6 anndata==0.7.1 umap==0.4.1 numpy==1.18.1 scipy==1.4.1 pandas==1.0.3 scikit-learn==0.22.2.post1 statsmodels==0.11.1 python-igraph==0.8.0

[ivirshup]

Did you have any missing values in your array before running combat?

[auesro]

None at all apparently:

In [2]: adata.shape
Out[2]: (538, 15036)
In [3]: 538*15036
Out[3]: 8089368
In [4]: np.sum(np.isnan(adata.X))
Out[4]: 0
In [5]: np.sum(~np.isnan(adata.X))
Out[5]: 8089368

[LuckyMD]

It might have sth to do with this line:

Found 3 genes with zero variance.

Maybe there are a couple of genes where the combat model can't be fit due to having 0 variance in a batch or 0 mean across a few batches? Could you check where the NaN values are in your adata.X post combat?

[auesro]

In theory Combat knows how to take care of zero variance genes according to the
code.
Well, post-Combat apparently NaNs are everywhere:

np.sum(np.isnan(adata_Combat.X))
Out[2]: 8089368

np.sum(~np.isnan(adata_Combat.X))
Out[3]: 0

This is really weird if only 3 genes have zero variance, right? Could it have anything to do with this warnings?:

Standardizing Data across genes.

Found 11 batches

Found 0 numerical variables:
	

Found 3 genes with zero variance.
Fitting L/S model and finding priors

Finding parametric adjustments

Adjusting data

/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/scanpy/preprocessing/_combat.py:338: RuntimeWarning: invalid value encountered in true_divide
  change = max((abs(g_new - g_old) / g_old).max(), (abs(d_new - d_old) / d_old).max())
/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/scanpy/preprocessing/_combat.py:338: RuntimeWarning: divide by zero encountered in true_divide
  change = max((abs(g_new - g_old) / g_old).max(), (abs(d_new - d_old) / d_old).max())

[auesro]

OK, I got rid of a few genes that were not expressed in the dataset (its a subset of cells from the full dataset) with sc.pp.filter_genes(adata, min_counts=1) and now the zero variance genes are gone but still same warning, same NaNs and same error when trying to run sc.pp.highly_variable_genes():

In [1]: sc.pp.combat(adata_Combat, key='sample')
/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/anndata/_core/anndata.py:21: FutureWarning: pandas.core.index is deprecated and will be removed in a future version.  The public classes are available in the top-level namespace.
  from pandas.core.index import RangeIndex
/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/sklearn/externals/six.py:31: FutureWarning: The module is deprecated in version 0.21 and will be removed in version 0.23 since we've dropped support for Python 2.7. Please rely on the official version of six (https://pypi.org/project/six/).
  "(https://pypi.org/project/six/).", FutureWarning)
scanpy==1.4.6 anndata==0.7.1 umap==0.4.1 numpy==1.18.1 scipy==1.4.1 pandas==1.0.3 scikit-learn==0.22.2.post1 statsmodels==0.11.1 python-igraph==0.8.0
Standardizing Data across genes.

Found 11 batches

Found 0 numerical variables:
	

Fitting L/S model and finding priors

Finding parametric adjustments

/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/scanpy/preprocessing/_combat.py:338: RuntimeWarning: divide by zero encountered in true_divide
  change = max((abs(g_new - g_old) / g_old).max(), (abs(d_new - d_old) / d_old).max())
Adjusting data

In [2]: np.sum(~np.isnan(adata_Combat.X))
Out[2]: 0

In [3]: np.sum(np.isnan(adata_Combat.X))
Out[3]: 7644442

In [4]: sc.pp.highly_variable_genes(adata_Combat)
extracting highly variable genes
Traceback (most recent call last):

  File "<ipython-input-4-a706aaf6f1f8>", line 1, in <module>
    sc.pp.highly_variable_genes(adata_Combat)

  File "/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/scanpy/preprocessing/_highly_variable_genes.py", line 235, in highly_variable_genes
    flavor=flavor,

  File "/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/scanpy/preprocessing/_highly_variable_genes.py", line 65, in _highly_variable_genes_single_batch
    df['mean_bin'] = pd.cut(df['means'], bins=n_bins)

  File "/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/pandas/core/reshape/tile.py", line 265, in cut
    duplicates=duplicates,

  File "/home/auesro/anaconda3/envs/Scanpy/lib/python3.7/site-packages/pandas/core/reshape/tile.py", line 381, in _bins_to_cuts
    f"Bin edges must be unique: {repr(bins)}.\n"

ValueError: Bin edges must be unique: array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan]).
You can drop duplicate edges by setting the 'duplicates' kwarg

[LuckyMD]

It's hard to say what's going on here. There seems to be sth happening in sc.pp.combat() so I would evaluate this separately from sc.pp.highly_variable_genes(). It would be important to have a minimal reproducible example it seems.

[auesro]

It's only 538 cells, so it's a small file, I can send it to you, if you want to have a look.

…

________________________________ From: MalteDLuecken <notifications@github.com> Sent: Monday, April 20, 2020 2:42:07 PM To: theislab/scanpy <scanpy@noreply.github.com> Cc: Augusto Escalante <ae_rodriguez_@hotmail.com>; Author <author@noreply.github.com> Subject: Re: [theislab/scanpy] Combat populates adata.X with NANs so sc.pp.highly_variable_genes function outputs error (#1172) It's hard to say what's going on here. There seems to be sth happening in sc.pp.combat() so I would evaluate this separately from sc.pp.highly_variable_genes(). It would be important to have a minimal reproducible example it seems. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#1172 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AFZKH4VCL573RRAU55AHS23RNQ7J7ANCNFSM4ML4AVXA>.

[LuckyMD]

I'm afraid i don't have the bandwidth right now...

@Marius1311 do you have any idea what's going on here? I believe you wrote the function, right?

[auesro]

Ok, I think I fixed it. Maybe it is mentioned somewhere or implicitly understood (sorry I am too new to this) but Combat will not do well if the number of cells in one of the batches is 1.
One of my batches mantained only 1 cell after filtering and subsetting, removing this one sample from the analysis solved the Combat problem: no NaNs, so that highly_variable_genes could worked as expected.

[LuckyMD]

You can't fit the variance on 1 cell, so that makes sense ;).

[gokceneraslan]

I think it's a nice corner case we should handle. Can you file another bug about having 1cell batches in combat or highly_variable_genes with batch_key option?

[auesro]

@LuckyMD of course your are right...
@gokceneraslan sure, I can, it would be nice just to get a warning or have Combat halt the processing on first checking that there is a batch with just 1 cell.