diff --git a/CHANGELOG.md b/CHANGELOG.md index 1d48dfc5..2e460fbd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -33,6 +33,7 @@ * Added ComBat-Seq method (PR #55). * Added `methods/scmerge2` component (PR #63). +* Added `metrics/ksim` component, reporting the `ksim_mean` and `ksim_accept_rate` metrics (PR #75). ## Minor changes diff --git a/src/metrics/ksim/config.vsh.yaml b/src/metrics/ksim/config.vsh.yaml new file mode 100644 index 00000000..ad0e9794 --- /dev/null +++ b/src/metrics/ksim/config.vsh.yaml @@ -0,0 +1,80 @@ +__merge__: /src/api/comp_metric.yaml +name: ksim +info: + metric_type: embedding + metrics: + - name: ksim_mean + label: kSIM mean + summary: Mean fraction of each cell's nearest neighbours that share its cell type + after batch correction. + description: | + The kSIM metric uses prior knowledge of cell type labels to assess whether local + neighbourhoods remain biologically consistent after batch correction (Li et al., + Nat Methods 2020). For each cell, its K nearest neighbours in the integrated + embedding are retrieved (the cell itself counts as its own first neighbour) and the + fraction of those neighbours that carry the same cell type label is recorded. + kSIM mean is the average of these per-cell fractions over the whole dataset. + + A high value means that cells of the same type stay locally grouped after + correction, whereas a low value indicates that the correction has disrupted true + biological structure, for example by overcorrecting batch effects. This variant is + continuous and therefore separates methods across the whole range of the scale. + + This implementation uses the `pegasus.calc_kSIM` function with its default + neighbourhood size of K=25. + references: + doi: + - 10.1038/s41592-020-0905-x + links: + homepage: https://pegasus.readthedocs.io/en/stable/ + documentation: https://pegasus.readthedocs.io/en/stable/api/pegasus.calc_kSIM.html + repository: https://github.com/lilab-bcb/pegasus + min: 0 + max: 1 + maximize: true + - name: ksim_accept_rate + label: kSIM acceptance rate + summary: Fraction of cells whose neighbourhood is almost entirely made up of their + own cell type after batch correction. + description: | + The kSIM acceptance rate is the thresholded companion of kSIM mean (Li et al., + Nat Methods 2020). The same per-cell fractions of same-cell-type nearest neighbours + are computed, but a cell is then counted as accepted only if that fraction reaches + at least `min_rate`. The acceptance rate is the proportion of accepted cells in the + dataset. + + Because it is a thresholded quantity, this variant is more stringent than kSIM mean: + it reports how many cells retain an essentially pure neighbourhood rather than how + pure neighbourhoods are on average. It can saturate near 0 for poorly separated + integrations, so it is best read alongside kSIM mean. + + This implementation uses the `pegasus.calc_kSIM` function with its defaults of + K=25 and min_rate=0.9, meaning that at least 23 of a cell's 25 neighbours must + share its cell type. + references: + doi: + - 10.1038/s41592-020-0905-x + links: + homepage: https://pegasus.readthedocs.io/en/stable/ + documentation: https://pegasus.readthedocs.io/en/stable/api/pegasus.calc_kSIM.html + repository: https://github.com/lilab-bcb/pegasus + min: 0 + max: 1 + maximize: true +resources: + - type: python_script + path: script.py + - path: /src/utils/read_anndata_partial.py +engines: + - type: docker + image: openproblems/base_python:1 + setup: + - type: python + pypi: + - pegasuspy + - setuptools<81 +runners: + - type: executable + - type: nextflow + directives: + label: [midtime, midmem, midcpu] diff --git a/src/metrics/ksim/script.py b/src/metrics/ksim/script.py new file mode 100644 index 00000000..7134fde7 --- /dev/null +++ b/src/metrics/ksim/script.py @@ -0,0 +1,61 @@ +import sys +import anndata as ad +import pegasus as pg +import pegasusio +from scipy.sparse import csr_matrix + +## VIASH START +par = { + 'input_integrated': 'resources_test/task_batch_integration/cxg_immune_cell_atlas/integrated_processed.h5ad', + 'input_solution': 'resources_test/task_batch_integration/cxg_immune_cell_atlas/solution.h5ad', + 'output': 'output.h5ad', +} + +meta = { + 'name': 'ksim', +} +## VIASH END + +sys.path.append(meta["resources_dir"]) +from read_anndata_partial import read_anndata + +n_threads = meta["cpus"] or -1 + +print('Read input...', flush=True) +adata = read_anndata(par['input_integrated'], obs='obs', obsm='obsm', uns='uns') +adata.obs = read_anndata(par['input_solution'], obs='obs').obs +adata.uns |= read_anndata(par['input_solution'], uns='uns').uns +print(adata, flush=True) + +print('Convert to pegasusio.MultimodalData...', flush=True) +adata.X = csr_matrix(adata.shape) +mmdata = pegasusio.MultimodalData(adata) + +print('Compute kSIM...', flush=True) +# K and min_rate are pegasus' own defaults, spelled out here so the metric stays +# stable across pegasus releases. Note that pegasus caps K at sqrt(n_obs), so the +# effective neighbourhood is smaller for datasets with fewer than 625 cells. +ksim_mean, ksim_accept_rate = pg.calc_kSIM( + mmdata, + attr='cell_type', + rep='emb', + K=25, + min_rate=0.9, + n_jobs=n_threads, +) +print('kSIM mean:', ksim_mean, flush=True) +print('kSIM acceptance rate:', ksim_accept_rate, flush=True) + +print('Create output AnnData object', flush=True) +output = ad.AnnData( + uns={ + 'dataset_id': adata.uns['dataset_id'], + 'normalization_id': adata.uns['normalization_id'], + 'method_id': adata.uns['method_id'], + 'metric_ids': [ 'ksim_mean', 'ksim_accept_rate' ], + 'metric_values': [ ksim_mean, ksim_accept_rate ] + } +) + +print('Write data to file', flush=True) +output.write_h5ad(par['output'], compression='gzip') diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml index ac1fe44c..2c0d1dd5 100644 --- a/src/workflows/run_benchmark/config.vsh.yaml +++ b/src/workflows/run_benchmark/config.vsh.yaml @@ -126,6 +126,7 @@ dependencies: - name: metrics/kbet - name: metrics/kbet_pg - name: metrics/kbet_pg_label + - name: metrics/ksim - name: metrics/lisi - name: metrics/pcr # data processors diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf index 970511ca..832f3c97 100644 --- a/src/workflows/run_benchmark/main.nf +++ b/src/workflows/run_benchmark/main.nf @@ -67,6 +67,7 @@ metrics = [ kbet, kbet_pg, kbet_pg_label, + ksim, lisi, pcr ]