Skip to content

Commit 7e413b1

Browse files
authored
Add ss_opm (#16)
* add ss_opm * fix ss_opm * add to wf * update script * address the review comments on ss_opm * Move `build_metadata()` and `to_sparse_csr()` into `ss_opm_common.py`, shared by both scripts instead of copied between them * Add `--cell_type_col`, so the `cell_ratio_*` features come from real cell types when a dataset has them. This task's file format does not carry cell types, so the default stays uniform * Add `--day_pattern` and document that a batch label which does not match yields day 0 rather than an error * Add `--n_epochs` and `--burnin_length_epoch`, capped during `viash test` with `info.test_default` * Report `ss_opm` as the method id rather than `ss_opm_predict` * Give `ss_opm_predict` a `test_setup` and build its model fixture in `test_resources.sh`
1 parent ce09e2f commit 7e413b1

11 files changed

Lines changed: 589 additions & 2 deletions

File tree

scripts/create_datasets/test_resources.sh

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -138,6 +138,19 @@ for name in bmmc_cite/normal bmmc_cite/swap bmmc_multiome/normal bmmc_multiome/s
138138
fi
139139
fi
140140

141+
echo "pre-train ss_opm on $name"
142+
if up_to_date $DATASET_DIR/$name/models/ss_opm/ $STATE; then
143+
echo " already up to date, skipping"
144+
else
145+
rm -rf $DATASET_DIR/$name/models/ss_opm/
146+
mkdir -p $DATASET_DIR/$name/models/ss_opm/
147+
viash run src/methods/ss_opm/ss_opm_train/config.vsh.yaml -- \
148+
--input_train_mod1 $DATASET_DIR/$name/train_mod1.h5ad \
149+
--input_train_mod2 $DATASET_DIR/$name/train_mod2.h5ad \
150+
--input_test_mod1 $DATASET_DIR/$name/test_mod1.h5ad \
151+
--output $DATASET_DIR/$name/models/ss_opm
152+
fi
153+
141154
# scbutterfly only does multiome (GEX<->ATAC)
142155
if [[ "$name" == bmmc_multiome/* ]]; then
143156
echo "pre-train scbutterfly on $name"
Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
__merge__: ../../../api/comp_method.yaml
2+
name: ss_opm
3+
label: SS-OPM
4+
summary: 1st place solution of the Kaggle Open Problems Multimodal Single-Cell Integration challenge.
5+
description: |
6+
Encoder-decoder MLP method using SVD-based dimensionality reduction for both inputs and
7+
targets, followed by batch-median correction. The encoder maps (optionally augmented)
8+
cell embeddings to a latent space; multiple decoder blocks predict target expression in
9+
the SVD-compressed space. The method was the winning solution of the NeurIPS 2021
10+
Open Problems Multimodal Single-Cell Integration Kaggle competition.
11+
references:
12+
doi:
13+
- 10.1101/2022.04.11.487796
14+
links:
15+
repository: https://github.com/shu65/open-problems-multimodal
16+
info:
17+
preferred_normalization: log_cp10k
18+
resources:
19+
- path: main.nf
20+
type: nextflow_script
21+
entrypoint: run_wf
22+
dependencies:
23+
- name: methods/ss_opm_train
24+
- name: methods/ss_opm_predict
25+
runners:
26+
- type: nextflow

src/methods/ss_opm/ss_opm/main.nf

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
workflow run_wf {
2+
take: input_ch
3+
main:
4+
output_ch = input_ch
5+
| ss_opm_train.run(
6+
fromState: ["input_train_mod1", "input_train_mod2", "input_test_mod1"],
7+
toState: ["input_model": "output"]
8+
)
9+
| ss_opm_predict.run(
10+
fromState: ["input_test_mod1", "input_model"],
11+
toState: ["output": "output"]
12+
)
13+
| map { tup ->
14+
[tup[0], [output: tup[1].output]]
15+
}
16+
17+
emit: output_ch
18+
}
Lines changed: 109 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,109 @@
1+
"""Helpers shared by ss_opm_train and ss_opm_predict."""
2+
3+
import numpy as np
4+
import pandas as pd
5+
import scipy.sparse
6+
7+
# the cell types the original ss_opm model was trained against. only used to name the
8+
# cell_ratio_* columns it expects; the ratios themselves are derived from the data when
9+
# cell type labels are available.
10+
CITE_CELL_TYPES = ["HSC", "EryP", "NeuP", "MasP", "MkP", "BP", "MoP"]
11+
12+
# number of batch singular-vector columns the cite model expects
13+
N_BATCH_SV = 8
14+
15+
16+
def to_sparse_csr(X):
17+
if scipy.sparse.issparse(X):
18+
return X.tocsr()
19+
return scipy.sparse.csr_matrix(X)
20+
21+
22+
def extract_day(batch, pattern=r"d(\d+)"):
23+
"""Pull the day out of a batch label.
24+
25+
The NeurIPS 2021 batches are named `s{site}d{day}`, e.g. `s1d2`. Datasets that
26+
label their batches differently yield NaN, which the caller fills with 0 -- the
27+
model then sees a single constant day rather than failing.
28+
"""
29+
return batch.astype(str).str.extract(pattern, expand=False).astype(float)
30+
31+
32+
def build_metadata(
33+
adata,
34+
task_type,
35+
cell_type_col=None,
36+
group_by_batch=True,
37+
day_pattern=r"d(\d+)",
38+
):
39+
"""Build the metadata frame ss_opm expects from an AnnData.
40+
41+
ss_opm was written against the Kaggle competition tables, which carry columns this
42+
task's API does not: `file_train_mod1.yaml` and `file_test_mod1.yaml` guarantee only
43+
`batch`. Everything else is either derived from `batch`, computed from the expression
44+
matrix, or filled with a neutral constant.
45+
46+
Parameters
47+
----------
48+
adata
49+
Input modality, with a `normalized` layer and `obs["batch"]`.
50+
task_type
51+
Either `"cite"` or `"multi"`; the cite model expects extra columns.
52+
cell_type_col
53+
Column in `adata.obs` holding cell type labels. When given, it drives both
54+
`cell_type` and the `cell_ratio_*` columns. When None -- the case for every
55+
dataset this task currently ships -- cell types are `"hidden"` and the ratios
56+
are uniform.
57+
group_by_batch
58+
Assign one group per batch. Set False to put every cell in group 0, which is
59+
what the predict path wants, since targets are absent and the group IDs are
60+
only used to look up target statistics.
61+
day_pattern
62+
Regex whose first capture group is the day within a batch label.
63+
"""
64+
obs = pd.DataFrame(index=adata.obs_names)
65+
66+
obs["batch"] = adata.obs["batch"].values
67+
obs["day"] = extract_day(adata.obs["batch"], day_pattern).fillna(0).values
68+
69+
# per-cell statistics from the normalized expression layer
70+
X = adata.layers["normalized"]
71+
X_dense = X.toarray() if scipy.sparse.issparse(X) else np.asarray(X, dtype=float)
72+
73+
obs["nonzero_ratio"] = (X_dense != 0).mean(axis=1)
74+
obs["nonzero_q25"] = np.percentile(X_dense, 25, axis=1)
75+
obs["nonzero_q50"] = np.percentile(X_dense, 50, axis=1)
76+
obs["nonzero_q75"] = np.percentile(X_dense, 75, axis=1)
77+
obs["mean"] = X_dense.mean(axis=1)
78+
obs["std"] = X_dense.std(axis=1)
79+
80+
if group_by_batch:
81+
batches = adata.obs["batch"].unique().tolist()
82+
obs["group"] = adata.obs["batch"].map({b: i for i, b in enumerate(batches)}).astype(int).values
83+
else:
84+
obs["group"] = 0
85+
86+
# cell type labels, when the caller can supply them
87+
if cell_type_col is not None and cell_type_col in adata.obs:
88+
obs["cell_type"] = adata.obs[cell_type_col].astype(str).values
89+
else:
90+
obs["cell_type"] = "hidden"
91+
92+
# donor and technology are not in this task's file format; gender_id defaults to 0
93+
obs["donor"] = 0
94+
obs["technology"] = "unknown"
95+
96+
if task_type == "cite":
97+
ratios = obs["cell_type"].value_counts(normalize=True)
98+
for cell_type in CITE_CELL_TYPES:
99+
obs[f"cell_ratio_{cell_type}"] = ratios.get(cell_type, 1.0 / len(CITE_CELL_TYPES))
100+
101+
batch_counts = adata.obs["batch"].value_counts()
102+
obs["cell_count"] = adata.obs["batch"].map(batch_counts).astype(float).values
103+
104+
# the originals are singular vectors of the full Kaggle batch matrix, which we
105+
# cannot reconstruct from a single dataset
106+
for i in range(N_BATCH_SV):
107+
obs[f"batch_sv{i}"] = 0.0
108+
109+
return obs
Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
__merge__: ../../../api/comp_method_predict.yaml
2+
name: ss_opm_predict
3+
4+
info:
5+
test_setup:
6+
with_model:
7+
input_model: resources_test/task_predict_modality/openproblems_neurips2021/bmmc_cite/swap/models/ss_opm
8+
arguments:
9+
- name: "--cell_type_col"
10+
type: string
11+
required: false
12+
description: |
13+
Column in `obs` holding cell type labels. ss_opm uses them for the `cell_ratio_*`
14+
features. This task's file format does not carry cell types, so it is unset by
15+
default and the ratios fall back to uniform; set it when plugging in a dataset
16+
that does have them.
17+
- name: "--day_pattern"
18+
type: string
19+
default: 'd(\d+)'
20+
description: |
21+
Regex whose first capture group is the collection day within a batch label. The
22+
default matches the NeurIPS 2021 `s{site}d{day}` naming. Batches that do not match
23+
get day 0, i.e. the model sees a single constant day rather than failing.
24+
resources:
25+
- type: python_script
26+
path: script.py
27+
- path: ../ss_opm_common.py
28+
engines:
29+
- type: docker
30+
image: openproblems/base_pytorch_nvidia:1
31+
setup:
32+
- type: docker
33+
run: pip install --no-cache-dir --no-deps git+https://github.com/shu65/open-problems-multimodal.git
34+
- type: python
35+
packages:
36+
- pyarrow
37+
- fastparquet
38+
runners:
39+
- type: executable
40+
- type: nextflow
41+
directives:
42+
label: [highmem, hightime, midcpu, highsharedmem, gpu]
Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
import sys
2+
import os
3+
import gc
4+
import pickle
5+
import numpy as np
6+
import pandas as pd
7+
import scipy.sparse
8+
import anndata as ad
9+
from ss_opm.model.encoder_decoder.encoder_decoder import EncoderDecoder
10+
11+
import torch
12+
device = 'cuda' if torch.cuda.is_available() else 'cpu'
13+
print(f'Using device: {device}', flush=True)
14+
15+
## VIASH START
16+
par = {
17+
'input_test_mod1': 'resources_test/task_predict_modality/openproblems_neurips2021/bmmc_cite/swap/test_mod1.h5ad',
18+
'input_model': 'output/models/ss_opm',
19+
'output': 'output/prediction.h5ad',
20+
'cell_type_col': None,
21+
'day_pattern': r'd(\d+)',
22+
}
23+
meta = {
24+
'name': 'ss_opm_predict',
25+
'resources_dir': 'src/methods/ss_opm',
26+
}
27+
## VIASH END
28+
29+
sys.path.append(meta['resources_dir'])
30+
from ss_opm_common import build_metadata, to_sparse_csr
31+
32+
# ---- Load task info ----
33+
with open(os.path.join(par['input_model'], 'task_info.pickle'), 'rb') as f:
34+
task_info = pickle.load(f)
35+
task_type = task_info['task_type']
36+
mod2 = task_info['mod2']
37+
dataset_id = task_info['dataset_id']
38+
print(f'Task type: {task_type}, mod2: {mod2}', flush=True)
39+
40+
# ---- Load test data ----
41+
print('Loading test data...', flush=True)
42+
input_test_mod1 = ad.read_h5ad(par['input_test_mod1'])
43+
test_inputs = to_sparse_csr(input_test_mod1.layers['normalized'])
44+
test_metadata = build_metadata(
45+
input_test_mod1,
46+
task_type,
47+
cell_type_col=par['cell_type_col'],
48+
group_by_batch=False,
49+
day_pattern=par['day_pattern'],
50+
)
51+
52+
# ---- Load model and preprocessing artifacts ----
53+
print('Loading model...', flush=True)
54+
with open(os.path.join(par['input_model'], 'pre_post_process.pickle'), 'rb') as f:
55+
pre_post_process = pickle.load(f)
56+
57+
model = EncoderDecoder(params=None)
58+
# PyTorch >=2.6 defaults weights_only=True, which blocks custom classes.
59+
# Patch torch.load to use weights_only=False for trusted local model files.
60+
import torch as _torch
61+
_orig_torch_load = _torch.load
62+
_torch.load = lambda *a, **kw: _orig_torch_load(*a, **{**kw, 'weights_only': False})
63+
model.load(os.path.join(par['input_model'], 'model'))
64+
_torch.load = _orig_torch_load
65+
model.params['device'] = device
66+
67+
mod2_var = pd.read_parquet(os.path.join(par['input_model'], 'mod2_var.parquet'))
68+
69+
# ---- Preprocess test inputs ----
70+
print('Preprocessing test data...', flush=True)
71+
preprocessed_test_inputs, _ = pre_post_process.preprocess(
72+
inputs_values=test_inputs,
73+
targets_values=None,
74+
metadata=test_metadata,
75+
)
76+
77+
# ---- Predict ----
78+
print('Predicting...', flush=True)
79+
y_pred = model.predict(
80+
x=test_inputs,
81+
preprocessed_x=preprocessed_test_inputs,
82+
metadata=test_metadata,
83+
)
84+
gc.collect()
85+
86+
# ---- Write output ----
87+
print('Writing output...', flush=True)
88+
# Prediction must be a sparse matrix to be compatible with all metrics.
89+
if not scipy.sparse.issparse(y_pred):
90+
y_pred = scipy.sparse.csr_matrix(y_pred)
91+
92+
output = ad.AnnData(
93+
layers={"normalized": y_pred},
94+
obs=input_test_mod1.obs,
95+
var=mod2_var,
96+
uns={
97+
"dataset_id": dataset_id,
98+
"method_id": "ss_opm",
99+
},
100+
)
101+
output.write_h5ad(par['output'], compression="gzip")
102+
print('Done!', flush=True)
Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
__merge__: ../../../api/comp_method_train.yaml
2+
name: ss_opm_train
3+
arguments:
4+
- name: "--n_epochs"
5+
type: integer
6+
default: 40
7+
description: Number of training epochs.
8+
info:
9+
test_default: 2
10+
- name: "--burnin_length_epoch"
11+
type: integer
12+
default: 10
13+
description: |
14+
Epochs before the training-length ratio starts ramping up. Must be below
15+
`--n_epochs`, otherwise every epoch stays at ratio 0.
16+
info:
17+
test_default: 0
18+
- name: "--cell_type_col"
19+
type: string
20+
required: false
21+
description: |
22+
Column in `obs` holding cell type labels. ss_opm uses them for the `cell_ratio_*`
23+
features. This task's file format does not carry cell types, so it is unset by
24+
default and the ratios fall back to uniform; set it when plugging in a dataset
25+
that does have them.
26+
- name: "--day_pattern"
27+
type: string
28+
default: 'd(\d+)'
29+
description: |
30+
Regex whose first capture group is the collection day within a batch label. The
31+
default matches the NeurIPS 2021 `s{site}d{day}` naming. Batches that do not match
32+
get day 0, i.e. the model sees a single constant day rather than failing.
33+
resources:
34+
- type: python_script
35+
path: script.py
36+
- path: ../ss_opm_common.py
37+
engines:
38+
- type: docker
39+
image: openproblems/base_pytorch_nvidia:1
40+
setup:
41+
- type: docker
42+
run: pip install --no-cache-dir --no-deps git+https://github.com/shu65/open-problems-multimodal.git
43+
- type: python
44+
packages:
45+
- pyarrow
46+
- fastparquet
47+
runners:
48+
- type: executable
49+
- type: nextflow
50+
directives:
51+
label: [highmem, hightime, midcpu, highsharedmem, gpu]

0 commit comments

Comments
 (0)