Diagonal integration¶

Diagonal integration combines RNA and ATAC measured in different cells, with no pairing between them. If your RNA and ATAC come from the same cells, as in 10x Multiome, use the vertical tutorial. This tutorial runs iNMF and online_iNMF on D28, with 6,408 RNA cells and 4,606 ATAC cells, and then on data in your own format.

Most diagonal methods read ATAC as gene-activity scores, made beforehand with a tool such as Signac or ArchR. GLUE and Seurat_v5 read the peak matrix, and MultiMAP and Seurat_v3 need both. Seurat_v5 also needs RNA and ATAC from the same cells. mtb.describe_layout("diagonal") lists each method's ATAC files.

Methods run on Linux. On macOS or Windows, open this notebook in Colab.

1. Install¶

This cell installs multibench-sc. It keeps the numpy and pandas that are already installed, so Colab needs no restart.

Details

On Colab, choose a GPU runtime before you run the notebook: Runtime -> Change runtime type -> T4 GPU. On a CPU runtime, an environment that has a smaller CPU build gets that build, and training methods run slower.

In [1]:

import numpy, pandas
%pip install -q multibench-sc numpy=={numpy.__version__} pandas=={pandas.__version__}

import numpy, pandas %pip install -q multibench-sc numpy=={numpy.__version__} pandas=={pandas.__version__}

Note: you may need to restart the kernel to use updated packages.

2. Download the data and the environments¶

mtb.data.fetch downloads D28 (137 MB) once. Each method runs in its own environment. mtb.env.install downloads the environments for iNMF and online_iNMF, with no conda needed. The download is 0.9 GB.

Details

state is PACKED for an environment downloaded now and have for one that was already there. Without dry_run=False, mtb.env.install downloads nothing and returns the plan with its sizes.

Environments go to mtb.config.DEFAULT.envs_dir. To use another disk, set it before this cell.

From a terminal, multibench env install --methods iNMF,online_iNMF --packed --run does the same. With --category diagonal instead of --methods, it installs every diagonal environment: 9 envs, 21.6 GB to download on a CPU host, 28.8 GB on a GPU host.

In [2]:

import pandas as pd
import multibench as mtb

mtb.data.fetch("D28")

import pandas as pd import multibench as mtb mtb.data.fetch("D28")

downloading D28 (137 MB) from https://github.com/DSichang/scMultiBench/releases/download/data-v1/D28.tar.gz ...

Out[2]:

PosixPath('/tmp/mbtut3/home/.cache/multibench/data')

In [3]:

METHODS = ["iNMF", "online_iNMF"]
envs = mtb.env.install(METHODS, dry_run=False)
pd.DataFrame(envs)[["env", "methods", "state"]]

METHODS = ["iNMF", "online_iNMF"] envs = mtb.env.install(METHODS, dry_run=False) pd.DataFrame(envs)[["env", "methods", "state"]]

[env] downloading scmb_r (single build) 0.9 GB ...

[env]   0.1 GB of 0.9 GB (57 MB/s)

[env]   0.2 GB of 0.9 GB (66 MB/s)

[env]   0.3 GB of 0.9 GB (68 MB/s)

[env]   0.4 GB of 0.9 GB (72 MB/s)

[env]   0.5 GB of 0.9 GB (74 MB/s)

[env]   0.6 GB of 0.9 GB (76 MB/s)

[env]   0.6 GB of 0.9 GB (77 MB/s)

[env]   0.7 GB of 0.9 GB (77 MB/s)

[env]   0.8 GB of 0.9 GB (78 MB/s)

[env]   0.9 GB of 0.9 GB (78 MB/s)

[env] unpacking scmb_r -> /tmp/mbtut3/home/.cache/multibench/envs/scmb_r ...

Out[3]:

	env	methods	state
0	scmb_r	[iNMF, online_iNMF]	PACKED

3. Run the methods¶

run_all runs each method on D28 and scores its output with the scIB metrics. res.summary has one row per method with its status, run time and scores. Higher is better for every metric.

Details

Each method writes an embedding: a table of numbers with one row per cell. run_all scores it against the cell-type labels of D28.

The clustering metrics ARI, NMI, ASW, iASW, iF1 and cLISI measure how well the embedding separates the cell types. The batch metrics ASW_batch, GC and iLISI measure how well the batches mix. They appear only when the data has several batches.

res.failures says why a method failed. params={"Method": {"key": value}} sets a method's parameters, and mtb.params_for lists them. Many methods take none.

mtb.load_batch("out/D28") reloads these results later without running anything.

In [4]:

res = mtb.run_all("D28", "diagonal", methods=METHODS, out_dir="out/D28")
res.summary

res = mtb.run_all("D28", "diagonal", methods=METHODS, out_dir="out/D28") res.summary

[run_all] iNMF (diagonal/D28) ...

Fetching the method scripts from PYangLab/scMultiBench into /tmp/mbtut3/home/.cache/multibench/scMultiBench_ref ...

[run_all]   -> CHAIN_OK (1224.6s) ARI 0.480

[run_all] online_iNMF (diagonal/D28) ...

[run_all]   -> CHAIN_OK (90.7s) ARI 0.402

Out[4]:

	method	status	run_sec	output_kind	emb_shape	n_tunable	label_order	label_order_confidence	batch_source	n_batches	...	ASW	iASW	iF1	cLISI	ASW_batch	GC	iLISI	label_order_note	caveat	reason
0	iNMF	CHAIN_OK	1224.6	embedding	[11014, 20]	0	rna_cty.csv+atac_cty.csv	0.9979	file_of_origin	2	...	0.4946	0.4832	0.4423	0.9649	0.7940	0.7758	0.6409	None
1	online_iNMF	CHAIN_OK	90.7	embedding	[11014, 20]	0	rna_cty.csv+atac_cty.csv	0.9965	file_of_origin	2	...	0.4666	0.4589	0.4189	0.9470	0.7777	0.7579	0.5614	None

2 rows × 22 columns

4. Plot¶

res.plot() draws the scores as a bubble table. Circle size shows the rank within a column, and bigger is better. The fill compares a value with the other rows in the same column.

Details

The lightest fill is the lowest value in this figure, not zero.

Metrics are grouped by family: blue for dimension reduction and clustering, green for batch correction. Each family starts with an Overall bar. Its length and colour both show the family score.

A column whose rows all hold the same value is drawn grey, and the note under the figure names it.

In [5]:

res.plot()

res.plot()

Out[5]:

5. Your own data¶

Your data needs raw RNA counts and ATAC as two AnnData objects, with a cell type for each cell. mtb.io.export_dataset writes them as a dataset folder, and run_all runs on that folder. iNMF and online_iNMF read ATAC as gene-activity scores.

Details: export

The command line writes the same folder from two .h5ad files:

multibench convert rna.h5ad mydata/MYDIAG --rna X --atac-from atac.h5ad \
    --atac-kind gene_activity --labels obs:celltype --category diagonal

For ATAC as a peak matrix, pass atac_kind="peak". GLUE and Seurat_v5 read peaks. mtb.describe_layout("diagonal") lists the files each method needs.

The folder name, MYDIAG, is the dataset name for run_all, and data_path is the folder that holds it.

mtb.scan("MYDIAG", "diagonal", data_path="mydata") checks the folder and the environments without running anything. Its reason column says what is missing.

overwrite=True replaces the files of an earlier run of this cell. Without it, export_dataset raises FileExistsError rather than replace a file.

Here 60% of D28's RNA cells and 60% of its ATAC cells take the place of your data:

In [6]:

import scanpy as sc

d = mtb.config.DEFAULT.data_path / "D28"
rna = mtb.io.read_canonical(d / "rna.h5")
rna.obs["celltype"] = pd.read_csv(d / "rna_cty.csv")["x"].values
atac = mtb.io.read_canonical(d / "atac_gas.h5")
atac.obs["celltype"] = pd.read_csv(d / "atac_cty.csv")["x"].values
rna = sc.pp.subsample(rna, fraction=0.6, random_state=0, copy=True)
atac = sc.pp.subsample(atac, fraction=0.6, random_state=0, copy=True)
rna, atac

import scanpy as sc d = mtb.config.DEFAULT.data_path / "D28" rna = mtb.io.read_canonical(d / "rna.h5") rna.obs["celltype"] = pd.read_csv(d / "rna_cty.csv")["x"].values atac = mtb.io.read_canonical(d / "atac_gas.h5") atac.obs["celltype"] = pd.read_csv(d / "atac_cty.csv")["x"].values rna = sc.pp.subsample(rna, fraction=0.6, random_state=0, copy=True) atac = sc.pp.subsample(atac, fraction=0.6, random_state=0, copy=True) rna, atac

Out[6]:

(AnnData object with n_obs × n_vars = 3844 × 10272
     obs: 'celltype',
 AnnData object with n_obs × n_vars = 2763 × 10272
     obs: 'celltype')

Write the folder, then run the same methods on it:

In [7]:

mtb.io.export_dataset(rna, "mydata/MYDIAG", atac=atac, atac_kind="gene_activity",
                      labels="obs:celltype", category="diagonal", overwrite=True)

mine = mtb.run_all("MYDIAG", "diagonal", methods=METHODS, data_path="mydata",
                   out_dir="out/MYDIAG")
mine.summary

mtb.io.export_dataset(rna, "mydata/MYDIAG", atac=atac, atac_kind="gene_activity", labels="obs:celltype", category="diagonal", overwrite=True) mine = mtb.run_all("MYDIAG", "diagonal", methods=METHODS, data_path="mydata", out_dir="out/MYDIAG") mine.summary

[run_all] iNMF (diagonal/MYDIAG) ...

[run_all]   -> CHAIN_OK (669.4s) ARI 0.370

[run_all] online_iNMF (diagonal/MYDIAG) ...

[run_all]   -> CHAIN_OK (74.0s) ARI 0.359

Out[7]:

	method	status	run_sec	output_kind	emb_shape	n_tunable	label_order	label_order_confidence	batch_source	n_batches	...	ASW	iASW	iF1	cLISI	ASW_batch	GC	iLISI	label_order_note	caveat	reason
0	iNMF	CHAIN_OK	669.4	embedding	[6607, 20]	0	rna_cty.csv+atac_cty.csv	1.0000	file_of_origin	2	...	0.4869	0.4838	0.4111	0.9475	0.8240	0.7507	0.7009	None
1	online_iNMF	CHAIN_OK	74.0	embedding	[6607, 20]	0	rna_cty.csv+atac_cty.csv	0.9944	file_of_origin	2	...	0.4689	0.4439	0.3845	0.9421	0.7406	0.7847	0.7007	None

2 rows × 22 columns

In [8]:

mine.plot()

mine.plot()

Out[8]:

6. Stored scores¶

The package ships stored scores for 12 methods on D28. load_results reads them and mtb.plot.bubble draws them, without running anything.

Details

source="rerun" reads the package's own runs of the methods. source="published", the default, reads the published scIB tables, which hold 14 methods for D28.

The stored scores used the leidenalg backend for Leiden clustering, and run_all uses igraph by default. The two backends can move ARI by up to about 0.1. To compare your runs with these scores, set mtb.config.DEFAULT.leiden_flavor = "leidenalg" before run_all.

In [9]:

long = mtb.load_results("diagonal", dataset="D28", source="rerun")
mtb.plot.bubble(long)

long = mtb.load_results("diagonal", dataset="D28", source="rerun") mtb.plot.bubble(long)

/tmp/ipykernel_2087109/4272824004.py:1: DegenerateRerunWarning: This re-run row has ARI below 0.01, while the published table has above 0.2 for the same method and dataset: Conos/D28 (re-run 0.2.1, ARI 0.0004). Such a row most likely comes from a failed re-run, for example a collapsed embedding or a wrong label order, not from the method itself. Drop it before ranking: df[df.method != 'Conos'].
  long = mtb.load_results("diagonal", dataset="D28", source="rerun")

Out[9]:

Troubleshooting¶

res.failures lists each method that failed or was skipped. Its error column ends with the method's error output.

Details

symptom	fix
env.install refuses on macOS or Windows	methods run only on Linux: use Colab or a Linux machine
a method needs an NVIDIA GPU	choose a GPU runtime on Colab, or a machine with a GPU
a warning that values are not whole numbers	export raw counts, for example with rna="layer:counts"
... matrix/data as cells x features	the matrix is transposed: export it again with mtb.io.export_dataset
a method times out	raise timeout= in run_all
low label_order_confidence	several label files fit the cell count: check label_order_candidates in res.results
batch metrics use the wrong batches	res.rescore(batch=my_vector) scores again without running the methods

Next steps¶

• mtb.cite(METHODS) returns the citations for the benchmark and the methods you ran.
• The other tutorials: vertical, mosaic, cross.
• The guides: run, evaluate, plot and discover methods.
• The interactive explorer has the full benchmark's rankings, with no install needed.