New metric: ksim

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
 

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]

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')

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

src/workflows/run_benchmark/main.nf

@@ -67,6 +67,7 @@ metrics = [
   kbet,
   kbet_pg,
   kbet_pg_label,
+  ksim,
   lisi,
   pcr
 ]