sc.pp.highly_variable_genes with flavor='seurat_v3' and batch_key: conflict with documentation
Nobody has claimed this yet.
- Dominant language
- Python
- Stars
- 2.6k
- Forks
- 779
- Avg merge
- 1d 4h
- Merged PRs (30d)
- 27
Description
- I have checked that this issue has not already been reported.
- I have confirmed this bug exists on the latest version of scanpy.
- (optional) I have confirmed this bug exists on the master branch of scanpy.
Hey, I've noticed another potential problem within the seurat_v3 flavor of sc.pp.highly_variable_genes(). The documentation of the batch_key argument says on how the genes are ranked
For all flavors, genes are first sorted by how many batches they are a HVG. For dispersion-based flavors ties are broken by normalized dispersion. If
flavor = 'seurat_v3', ties are broken by the median (across batches) rank based on within-batch normalized variance.
However, when genes are sorted after computing everything, the seurat_v3 method sorts first by the median ranks and then by how many batches a gene is highly variable (contrary to what the docstring says):
For comparison, the other flavors sort the other way around (as the docstring says):
Not sure which sorting would be correct here to make seurat_v3 match the original behavior in Seurat, but I found these comments by @adamgayoso in #1204 that mention that when using batch_key in the current implementation, results don't match Seurat.. maybe the sorting order is the problem here?
https://github.com/theislab/scanpy/pull/1204#issuecomment-645700601
https://github.com/theislab/scanpy/pull/1204#issuecomment-663879113
Here is a piece of code that shows the top5 genes after seurat_v3 gene selection and sorting by highly_variable_nbatches (which should then all be HVG, but they are not)
import numpy as np
import scanpy as sc
import anndata
import sys
sys.path.append("scanpy/preprocessing")
from _utils import _get_mean_var
np.random.seed(42)
adata = anndata.AnnData(np.random.randint(0,5,(100,100)))
adata.obs['batch'] = np.random.randint(0,5,(100))
adata.obs['batch'] = adata.obs['batch'].astype('category')
n_top_genes = 50
adata = adata.copy()
sc.pp.highly_variable_genes(adata,n_top_genes=n_top_genes,flavor='seurat_v3',inplace=True,batch_key='batch')
adata.var.sort_values(['highly_variable_nbatches'], ascending=[False]).iloc[:5,:]
highly_variable highly_variable_rank means variances variances_norm highly_variable_nbatches
87 True 8.0 1.76 2.446869 1.232373 5
9 False 28.0 1.96 2.281212 1.159891 5
78 True 24.0 1.95 2.209596 1.124666 5
30 True 19.0 2.00 2.202020 1.134560 4
14 False 25.0 2.14 2.162020 1.088266 4
Versions
sc.logging.print_versions()
WARNING: If you miss a compact list, please tryprint_header!
anndata 0.7.5
scanpy 1.7.1.dev2+g8c469411
sinfo 0.3.1
PIL 8.1.0
anndata 0.7.5
cffi 1.14.4
constants NA
cycler 0.10.0
cython_runtime NA
dateutil 2.8.1
get_version 2.1
google NA
h5py 2.10.0
highs_wrapper NA
joblib 1.0.0
kiwisolver 1.3.1
legacy_api_wrap 1.2
llvmlite 0.35.0
matplotlib 3.3.3
mpl_toolkits NA
natsort 7.1.0
numba 0.52.0
numexpr 2.7.2
numpy 1.19.5
packaging 20.8
pandas 1.2.0
pkg_resources NA
pyparsing 2.4.7
pytz 2020.5
scanpy 1.7.1.dev2+g8c469411
scipy 1.6.0
sinfo 0.3.1
sitecustomize NA
six 1.15.0
sklearn 0.24.0
skmisc 0.1.3
tables 3.6.1
typing_extensions NA
yaml 5.3.1
Python 3.8.0 (default, Oct 28 2019, 16:14:01) [GCC 8.3.0]
Linux-3.10.0-957.el7.x86_64-x86_64-with-glibc2.27
40 logical CPU cores, x86_64
Session information updated at 2021-03-10 17:37
Contributor guide
First steps
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- Open a pull request that references the issue number.
Research direction
Start with the seurat_v3 batch_key sorting logic in scanpy/preprocessing/_highly_variable_genes.py around lines 139–144, then compare it with the documented ordering and the other flavor around lines 505–510. Run the reproducer from the issue and review the linked #1204 discussion; done means the intended ordering is decided, documented, and covered by a regression test against the reported behavior.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- python
- Domain
- bioinformatics
- Issue type
- Bug
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Activity status
- Stale
- Clarity
- Mostly clear
- Newbie friendliness
- 35/100