Single-Cell Analysis in Python. Scales to >1M cells.

Hi everyone,

a while back, @giovp asked me to create a pull request that integrates analytical Pearson residuals in scanpy. We already discussed a bit with @LuckyMD and @ivirshup over at berenslab/umi-normalization#1 how to structure it, and now I made a version that should be ready for review.

As discussed earlier, this pull request implements two core methods:

• sc.pp.normalize_pearson_residuals(), which applies the method to adata.X. Overall, the function is very similar in structure to sc.pp.normalize_total() (support for layers, inplace operation etc).
• sc.pp.highly_variable_genes(flavor='pearson_residuals'), which selects genes based on Pearson residual variance. The "inner" function _highly_variable_pearson_residuals() is structured similarly to _highly_variable_seurat_v3() (support for multiple batches, median ranks for tie breaking). It includes the chunksize argument to allow for memory-efficient computation of the residual variance.

We discussed quite a lot how to implement a third function that would bundle gene selection, normalization by analytical residuals and PCA. This PR includes the two options that emerged at the end of that discussion, so now we have to choose ;)

• sc.pp.recipe_pearson_residuals() which does HVG selection and normalization both via Pearson residuals prior to PCA
• sc.pp.normalize_pearson_residuals_pca() which applies any HVG selection if the user previously added one to the adata object, and then normalizes via Pearson residuals and does PCA.

Both functions retain the raw input counts as adata.X and add fields for PCA/Normalization/HVG selection results (or return them) as applicable, most importantly the X_pca in adata.obsm['pearson_residuals_X_pca']. I hope this addresses some of the issues we discussed over at the other repo in a scanpy-y way.

Let me know what you think and where you think improvements are needed!

Cheers, Jan.

I know this (quite ancient) pull request has been open (#403), but I wasn't sure on its status. I think the consensus was to wait for sklearn to integrate the necessary changes? If that's still the case, then please feel free to remove this PR.

Here I make use of scipy's extremely nifty LinearOperator class to customize the dot product functions for an input sparse matrix. In this case, the 'custom' dot product performs implicit mean centering.

In my benchmarks, performing implicit mean centering in this way does not affect the runtime whatsoever. However, this approach has to use svds, for which randomized SVD is not implemented. So we have to use 'arpack', which can be significantly slower (but not intractably so.... in my hands, I could still do PCA on datasets of 200k+ cells in minutes, and it sure beats densifying the data, if you want more thorough benchmarks I am happy to generate them!).

The way I incorporated this functionality into scanpy/preprocessing/_simple.py might be questionable, and would love any suggestions or advice on how to better integrate this if there is interest in pushing this PR through. Let me know!

@falexwolf wrote:

Wow, this looks great! One remark for future PRs: We’re migrating to a new documentation style using type annotations

I'm still not convinced that we should use type annotations for Scanpy toplevel functions. People use Scanpy in Jupyter Lab and notebooks and not in Pycharm. Hence, there is no gain in the annotations, by contrast, the function annotations simply look super complicated and it's no longer feasible to grasp at first sight what's going on. This also regards the output of the help in Jupyter Lab and Notebooks.

So, while I think that for AnnData and everything in the background, type annotations may make sense for a few developers (not for me, as I'm doing everything on remote servers using emacs), it doesn't make sense for the Scanpy user.

Also, all the other big packages I work with all the time simply don't have it (numpy, seaborn, pandas, tensorflow) and it makes it harder and lengthier for contributors to contribute if they need to go through it.

Finally, I'm still not happy about how the automatically generated docs from the type annotation look: which is from Clearly, the automatically generated line with Union[...] is just way too complicated for a human to make sense of. The mix of auto-generated types in the docs and the manual annotations also looks inhomogeneous.

So, please let's stay away from having more type annotations and corresponding docstrings at this stage and let's simply continue imitating what all the major packages are doing.

Also: regarding your comment about the use of '``' vs. '' in the docs: again, I think it leads to an inhomogeneous appearance to have *two* types of markup for code-related things. I agree that the read-the-docs italicized default style for '' might be supoptimal, and I'll work on that if there is some time. But in general, I think there should be essentially one markup for code, as it's done in the tensorflow docs and a couple of other examples.

Happy to also discuss offline, @flying-sheep ;).

What stays the same:

• pip install scanpy
• pip install .
• pip install git+https://...
• you can install your deps with conda
• you can do a dev install

What changes:

• Please check the install docs, in short:
  • pip install -e .[every,single,extra] → flit install -s for dev installs
  • beni pyproject.toml > environment.yml for conda
• Extremely simple flit build and flit publish. Maybe install keyring to store your publish password, and you know everything you need to.
• flit build doesn’t clutter your dev directory with build/ and *.egg-info/ junk, it just creates dist/scanpy-*{.whl,.tar.gz}.
• No more obscure stuff nobody understands (MANIFEST.in, package_data, …)
• Centralized setup configuration in pyproject.toml instead of spread over multiple files

I tried to collect in one file the code used for plotting functions that use matplotlib scatter like sc.pl.tsne, sc.pl.pca and sc.pl.umap and others.

Also, I tried to annotate the code and improve the readability.

Currently, the code is on a separate file called scatter.py and not integrated into the API as this facilitates comparison with previous code.

Besides readability the proposed code can:

• Plot a large number of plots in multiple columms (instead of a long row of plots)
• Pass arguments directly to matplotlib.pyplot.scatter like vmax and vmin to adjust the color scale. When plotting multiple plots, this is useful to have a consistent range of values)

See cells 15 and 15 in this example: https://gist.github.com/fidelram/8b43f786e7519bcfb7ffc0d5ccdbb0fe If the admins would like to merge these changes I can replaced the previous functions.

An example on how to use the code:

import scanpy.plotting.tools.scatter as spl
spl.tsne(adata, color='louvain')

Further examples here

As discussed with @falexwolf this PR introduces a new Ingest class to process new small pieces of data.

sc.pp.neighbors(adata) # adata is huge with 1M observations

ingest = sc.Ingest(adata) # represents the existing data, learned annotations, structure and exposes it to functionality that allows to ingest new data very quickly

adata_small.obsm['X_model'] = model(adata_small.X)

ingest.neighbors(adata_small) # adata_small with just 1000 observations

now, we have the updated neighbors graph with 1,001,000 observations we want to do the same things as always

by leveraging the neighbors of the new data within the old data, map the new data into the embedding (umap), by just computing a correction to the existing embedding: a new data point gets the mean position of its k nearest neighbors

ingest.umap(adata_small)

update the clustering (mapping the 1000 observations into the existing clusters): a new data point maps into a cluster if the majority of its neighbors is a member of the cluster

ingest.louvain(adata_small)

Could you add Moran's I calculation to Scanpy? It could be used in scIB and to also find variable genes across embedding (could be an alternative to SEMITONES that takes a while to be computed).

Added restrict_to parameter to leiden by using louvain code as template. Tests are not yet provided.

A simple example of execution and checks:

# First split on cluster 4
sc.tl.leiden(adata, restrict_to=('leiden_res0.4', ['4']), resolution=0.6,
    key_added='leiden_res0.4_4_sub')

# Additional split
sc.tl.leiden(adata, restrict_to=('leiden_res0.4_4_sub', ['1', '2', '3', '4,4']),
    resolution=0.6, key_added='leiden_res0.4_4_add_sub')

# All partitions together
sc.pl.tsne(adata, color=['leiden_res0.4', 'leiden_res0.4_4_sub',
    'leiden_res0.4_4_add_sub'])

# Partition size check
## Original size of clusters
adata.obs['leiden_res0.4'].value_counts()
0    932
1    853
3    676
2    676
4    338
5     57
Name: leiden_res0.4, dtype: int64

# Check if first split is correct (can be iterated for subsequent splits)
## Assignment of samples in original clusters to subsplit clusters
adata.obs.loc[(adata.obs['leiden_res0.4'].isin(['4'])),
    'leiden_res0.4_4_sub'].value_counts()
4,0    103
4,1     68
4,2     66
4,3     57
4,4     44
5        0
3        0
2        0
1        0
0        0
Name: leiden_res0.4_4_sub, dtype: int64
## Assignment of samples not in original clusters to subsplit clusters
adata.obs.loc[~(adata.obs['leiden_res0.4'].isin(['4'])),
    'leiden_res0.4_4_sub'].value_counts()
0      932
1      853
3      676
2      676
5       57
4,4      0
4,3      0
4,2      0
4,1      0
4,0      0
Name: leiden_res0.4_4_sub, dtype: int64

...

Hello,

It would appear that louvain-igraph has been obsoleted in favour of leidenalg, and the author makes a persuasive case as to the superiority of the new approach. To my untrained eye, the algorithm is conceptually similar to the Louvain modification used by Seurat, but introduces an extra collapsed network refinement step.

it should be easy to support this in Scanpy - the syntax appears to be identical to the old louvain innards, and I was able to construct a very minimal dummy function for testing by taking the key bits of sc.tl.louvain() and replacing louvain. with leidenalg.:

import leidenalg
import numpy as np
import pandas as pd
from scanpy import utils
from natsort import natsorted

def leiden(adata, use_weights=False, resolution=1, iterations=-1):
	g = utils.get_igraph_from_adjacency(adata.uns['neighbors']['connectivities'], directed=True)
	weights = None
	if use_weights:
		weights = np.array(g.es["weight"]).astype(np.float64)
	part = leidenalg.find_partition(
		g, leidenalg.RBConfigurationVertexPartition, 
		resolution_parameter = resolution, weights = weights, 
		n_iterations = iterations,
	)
	groups = np.array(part.membership)
	adata.obs['louvain'] = pd.Categorical(
		values=groups.astype('U'),
		categories=natsorted(np.unique(groups).astype('U')),
	)

As such, replacing any louvain. with leidenalg. in sc.tl.louvain() would do most of the work. Probably the only new thing that would need support would the the n_iterations parameter in leidenalg.find_partition(). The default value is 2, positive values control how many passes of the algorithm are performed. -1 just makes it run until it fails to improve the clustering.

I'm adding that expression atlas downloader now (#489), and wondering where the files should go.

pbmc68k_reduced and toggleswitch put the datasets relative to where scanpy is installed (via __file__). All other functions place the data relative to where the python process was started.

While I like not storing the same files all over a filesystem, I'm not sure in the scanpy installation directory is the right place to be storing data.

Thoughts?

Hi,

I would like to extract the expression matrix (genes and counts) from a h5ad file. How can I do this? I have searched the documentation but I couldn't find anything about this (maybe I missed it).

• [ X] Additional function parameters / changed functionality / changed defaults?

Hi ScanPy team

I emailed @ivirshup but others should be involved I think.

This function would be useful if we could specify the number of threads to use: https://scanpy.readthedocs.io/en/stable/generated/scanpy.tl.rank_genes_groups.html

Based on the number of items in the "groupby" field, we could use a basic split-merge approach here: each thread would take several of these items, the calculations are entirely independent of one another, and then when each is completed we would join + concatenate the results.

I'm happy to help write up a PR help (or participate), but I'd like to hear if this is something you'd be willing to prioritize. (It's related to a project whereby Fabian is the PI.)

Best, Evan

• [x] I have checked that this issue has not already been reported.
• [x] 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.

Note: Please read this guide detailing how to provide the necessary information for us to reproduce your bug.

Hi, everyone: Many users probably do not rely on pp.normalize_total for downstream analysis, but I found a strange default behavior that I think is worth mentioning. pp.normalize_total() normalized my .layers['counts'] as well The documentation is a bit murky; not sure if that is the expected behavior when layer is unspecified, but such default behavior would undermine anyone who wishes to save the count information before RPKM normalization.

Minimal code sample (that we can copy&paste without having any data)

# Your code here
adata = sc.datasets.pbmc3k()
adata.layers['counts'] = adata.X
cell = adata.obs.index[1]
adata.var['mt'] = adata.var_names.str.startswith('MT-')  # annotate the group of mitochondrial genes as 'mt'
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)

print("Run 1: initial values after simple processing: ")
print('sum of count layer in designated cell: ', adata[cell,:].layers['counts'].sum())
print('obs[total_counts] value in cell: ', adata[cell,:].obs['total_counts'][0])
print('.X.sum() value in cell: ', adata[cell,:].X.sum())
print('sum of count layer of MALAT1 in cell: ', adata[cell,'MALAT1'].layers['counts'])
print('.X value of MALAT1 in cell: ', adata[cell,'MALAT1'].X)

print("\nRun 2: after sc.pp.normalize_total: ")
sc.pp.normalize_total(adata, target_sum=1e4)
print('sum of count layer in designated cell: ', adata[cell,:].layers['counts'].sum()) # Note that this changed too
print('obs[total_counts] value in cell: ', adata[cell,:].obs['total_counts'][0])
print('.X.sum() value in cell: ', adata[cell,:].X.sum())
print('sum of count layer of MALAT1 in cell: ', adata[cell,'MALAT1'].layers['counts'])
print('.X value of MALAT1 in cell: ', adata[cell,'MALAT1'].X)

adata = sc.datasets.pbmc3k()
adata.layers['counts'] = adata.X
cell = adata.obs.index[1]
adata.var['mt'] = adata.var_names.str.startswith('MT-')  # annotate the group of mitochondrial genes as 'mt'
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)

print("\nRun 3: normalization, specifing argument layer=None")
sc.pp.normalize_total(adata, target_sum=1e4, layer = None)
print('sum of count layer in designated cell: ', adata[cell,:].layers['counts'].sum())
print('obs[total_counts] value in cell: ', adata[cell,:].obs['total_counts'][0])
print('.X.sum() value in cell: ', adata[cell,:].X.sum())
print('sum of count layer of MALAT1 in cell: ', adata[cell,'MALAT1'].layers['counts'])
print('.X value of MALAT1 in cell: ', adata[cell,'MALAT1'].X)

#Output:
Run 1: initial values after simple processing: 
sum of count layer in designated cell:  4903.0
obs[total_counts] value in cell:  4903.0
.X.sum() value in cell:  4903.0
sum of count layer of MALAT1 in cell:    (0, 0)	142.0
.X value of MALAT1 in cell:    (0, 0)	142.0

Run 2: after sc.pp.normalize_total: 
normalizing counts per cell
    finished (0:00:00)
sum of count layer in designated cell:  10000.049
obs[total_counts] value in cell:  4903.0
.X.sum() value in cell:  10000.049
sum of count layer of MALAT1 in cell:    (0, 0)	289.61862
.X value of MALAT1 in cell:    (0, 0)	289.61862

Run 3: normalization, specifing argument layer=None
normalizing counts per cell
    finished (0:00:00)
sum of count layer in designated cell:  10000.049
obs[total_counts] value in cell:  4903.0
.X.sum() value in cell:  10000.049
sum of count layer of MALAT1 in cell:    (0, 0)	289.61862
.X value of MALAT1 in cell:    (0, 0)	289.61862

Versions

scanpy docs:

n_jobs : Optional[int] (default: None) Number of threads to use in training. All cores are used by default.

magic docs:

n_jobs (integer, optional, default: 1) – The number of jobs to use for the computation. If -1 all CPUs are used. If 1 is given, no parallel computing code is used at all, which is useful for debugging.

scanpy code:

https://github.com/scverse/scanpy/blob/536ed15bc73ab5d1131c0d530dd9d4f2dc9aee36/scanpy/external/pp/_magic.py#L164

https://github.com/scverse/scanpy/blob/536ed15bc73ab5d1131c0d530dd9d4f2dc9aee36/scanpy/_settings.py#L82

I'm guessing the scanpy docs are wrong when they say "All cores are used by default." ?

Bug: var_names -parameter for sc.tl.dendrogram -function is not used properly.

Currently:

• Hierarchical clustering is calculated on all of the var_names (genes) when var_names is not None

Fix:

• Subset of genes defined by var_names is now used

In addition:

• When all of the values of some row of rep_df (or mean_df) are equal, df.T.corr() is not defined for that row resulting in NaNs in correlation matrix.
• This is quite common with a subset of genes var_names, e.g. all 0 in all cells (these cells have already passed quality control/filtering at this point of downstream analysis).
• This throws an error in distance.squareform(1-corr_matrix): ValueError: Distance matrix 'X' must be symmetric.
• Fix: In this case add 'dummy' feature rep_df["dummy"] = -1 to make sure that at least one feature in a row is distinct.
• Notice, this addition affects (increases) the correlation between the rows. However, it should affect all rows equally and hence the hierarchy stays as is.

• [x] I have checked that this issue has not already been reported.
• [x] 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.

Note: Please read this guide detailing how to provide the necessary information for us to reproduce your bug.

Minimal code sample (that we can copy&paste without having any data)

Write any anndata with pearson residuals in uns

ad_all.write(filename='output/10x_h5/ad_all_2cello.h5ad')

The pearson_residual_df looks like this, with 38291 rows (obs) and 5000 columns (features) :

{'theta': 100,
 'clip': None,
 'computed_on': 'adata.X',
 'pearson_residuals_df': gene_name                             A2M  AADACL2-AS1      AAK1     ABCA1  \
 barcode                                                                      
 GAACGTTCACACCGAC-1-placenta_81  -1.125285    -1.159130 -3.921314 -2.533474   
 TATACCTGTTAGCTAC-1-placenta_81  -1.091364     3.267127 -1.806667 -2.109586   
 CTCAAGAGTGACTGTT-1-placenta_81  -1.074943    12.272920 -1.948798 -2.735791   
 TTCATTGTCACGAACT-1-placenta_81  -1.098699    -1.131765  3.481171  4.472371   
 TATCAGGCAGCTCATA-1-placenta_81  -1.107734    -1.141064 -0.571775 -2.813671   
 ...                                   ...          ...       ...       ...   
 CACAACATCGGCGATC-1-placenta_314 -0.115585    -0.119107 -0.434686 -0.303945   
 AGCCAGCGTGCCCAGT-1-placenta_314 -0.097424    -0.100394 -0.366482 -0.256219   
 CCGGTGAGTGTTCGAT-1-placenta_314 -0.110334    -0.113696 -0.414971 -0.290148   
 AGGTCATAGCCTGACC-1-placenta_314 -0.115585    -0.119107 -0.434686 -0.303945   
 TTTATGCCAAAGGGTC-1-placenta_314 -0.112876    -0.116316 -0.424515 -0.296827 

Unable to create attribute (object header message is too large)

Above error raised while writing key 'pearson_residuals_df' of <class 'h5py._hl.group.Group'> to /

Versions

• [x] I have checked that this issue has not already been reported.
• [x] I have confirmed this bug exists on the latest version of scanpy.
• [x] (optional) I have confirmed this bug exists on the master branch of scanpy.

When sorting the adata object using adata.var.sort_index(inplace=True) the adata.X is not sorting accordingly (see example below).

import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib.pyplot as plt
sc.settings.verbosity = 3
sc.logging.print_header()

adata = sc.datasets.pbmc3k()
adata.to_df()['MPO'].sum()
adata.var.sort_index(inplace=True)
adata.to_df()['MPO'].sum()

Versions