diff --git a/CHANGELOG.md b/CHANGELOG.md index b11ecf4d5..5dba76e76 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -193,6 +193,7 @@ * `run_benchmark`: write the commit the workflow ran from and the launch time into `task_info.yaml`, instead of publishing `_viash.yaml` verbatim (PR #136). +* Remove `metrics/bras` , `metrics/n_inconsitent_peaks` and `control_methods/shuffle_by_cell_type` and `control_methods/batch` ## BUG FIXES @@ -260,3 +261,5 @@ * Fix check cell_type_in_both to > 0 (instead of 1) and some spelling mistakes in print. + +* Updated anndata version to 0.12 in lisi metric to support NULL writing (PR #139) \ No newline at end of file diff --git a/scripts/run_benchmark/wehi_hpc/images.txt b/scripts/run_benchmark/wehi_hpc/images.txt index eb63d04e5..4eb946518 100644 --- a/scripts/run_benchmark/wehi_hpc/images.txt +++ b/scripts/run_benchmark/wehi_hpc/images.txt @@ -2,8 +2,6 @@ docker://ghcr.io/openproblems-bio/openproblems/utils/extract_uns_metadata:build_ docker://ghcr.io/openproblems-bio/task_cyto_batch_integration/control_methods/no_integration:build_setup_run_hpc docker://ghcr.io/openproblems-bio/task_cyto_batch_integration/control_methods/perfect_integration:build_setup_run_hpc docker://ghcr.io/openproblems-bio/task_cyto_batch_integration/control_methods/shuffle_integration_globally:build_setup_run_hpc -docker://ghcr.io/openproblems-bio/task_cyto_batch_integration/control_methods/shuffle_integration_within_batch:build_setup_run_hpc -docker://ghcr.io/openproblems-bio/task_cyto_batch_integration/control_methods/shuffle_integration_within_cell_type:build_setup_run_hpc docker://ghcr.io/openproblems-bio/task_cyto_batch_integration/methods/batchadjust_all_controls:build_setup_run_hpc docker://ghcr.io/openproblems-bio/task_cyto_batch_integration/methods/batchadjust_one_control:build_setup_run_hpc docker://ghcr.io/openproblems-bio/task_cyto_batch_integration/methods/combat:build_setup_run_hpc diff --git a/src/control_methods/shuffle_integration_within_batch/config.vsh.yaml b/src/control_methods/shuffle_integration_within_batch/config.vsh.yaml deleted file mode 100644 index bc9103896..000000000 --- a/src/control_methods/shuffle_integration_within_batch/config.vsh.yaml +++ /dev/null @@ -1,23 +0,0 @@ -__merge__: /src/api/comp_control_method.yaml -name: shuffle_integration_within_batch -label: Shuffle Integration Within Batches -summary: Randomly reassign cells to any samples within the same batch. -description: | - This negative-control method randomly shuffles cells within each batch. - Cells lose their sample and cell type identity but remain in their original batch. - - Example: A cell from mouse 3 batch 1 may be reassigned to mouse 4 batch 1, - but not to another batch (e.g., mouse 6 batch 2). -# status: disabled -resources: - - type: python_script - path: script.py - - path: /src/control_methods/utils.py -engines: - - type: docker - image: openproblems/base_python:1 -runners: - - type: executable - - type: nextflow - directives: - label: [lowtime,lowmem,lowcpu] diff --git a/src/control_methods/shuffle_integration_within_batch/script.py b/src/control_methods/shuffle_integration_within_batch/script.py deleted file mode 100644 index 928425f9b..000000000 --- a/src/control_methods/shuffle_integration_within_batch/script.py +++ /dev/null @@ -1,66 +0,0 @@ -import sys - -import anndata as ad - -## VIASH START -par = { - "input_unintegrated": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad", - "output_integrated_split1": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/control_integrated_split1.h5ad", - "output_integrated_split2": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/control_integrated_split2.h5ad", -} -meta = { - "name": "shuffle_integration_within_batch", - "resources_dir": "src/control_methods", -} -## VIASH END - -print("Importing helper functions", flush=True) -sys.path.append(meta["resources_dir"]) -from utils import _randomize_features - -print("Reading and preparing input files", flush=True) -adata = ad.read_h5ad(par["input_unintegrated"]) -adata_split1 = adata[(adata.obs.is_control > 0) | (adata.obs.split == 1)].copy() -adata_split2 = adata[(adata.obs.is_control > 0) | (adata.obs.split == 2)].copy() - -print("Randomise features - split 1", flush=True) -adata_split1.obs["batch_str"] = adata_split1.obs["batch"].astype(str) -integrated = _randomize_features( - adata_split1.layers["preprocessed"], - partition=adata_split1.obs["batch"], -) - -# create new anndata -output_split1 = ad.AnnData( - obs=adata_split1.obs[[]], - var=adata_split1.var[[]], - layers={"integrated": integrated}, - uns={ - "dataset_id": adata_split1.uns["dataset_id"], - "method_id": meta["name"], - "parameters": {}, - }, -) - -print("Randomise features - split 2", flush=True) -adata_split2.obs["batch_str"] = adata_split2.obs["batch"].astype(str) -integrated = _randomize_features( - adata_split2.layers["preprocessed"], - partition=adata_split2.obs["batch"], -) - -# create new anndata -output_split2 = ad.AnnData( - obs=adata_split2.obs[[]], - var=adata_split2.var[[]], - layers={"integrated": integrated}, - uns={ - "dataset_id": adata_split2.uns["dataset_id"], - "method_id": meta["name"], - "parameters": {}, - }, -) - -print("Write output AnnData to file", flush=True) -output_split1.write_h5ad(par["output_integrated_split1"], compression="gzip") -output_split2.write_h5ad(par["output_integrated_split2"], compression="gzip") diff --git a/src/control_methods/shuffle_integration_within_cell_type/config.vsh.yaml b/src/control_methods/shuffle_integration_within_cell_type/config.vsh.yaml deleted file mode 100644 index a7825e41d..000000000 --- a/src/control_methods/shuffle_integration_within_cell_type/config.vsh.yaml +++ /dev/null @@ -1,24 +0,0 @@ -__merge__: /src/api/comp_control_method.yaml -name: shuffle_integration_within_cell_type -label: Shuffle Integration Within Cell Type -summary: Randomly reassign cells to any cell types -description: | - This negative-control method randomly shuffles cells within each cell type. - Cells retain their cell type identity but lose their sample and batch identity. - - Example: A B cell remains a B cell but if may be reassigned from mouse 3 batch 1 to - mouse 4 batch 1 or mouse 6 batch 2. - -# status: disabled -resources: - - type: python_script - path: script.py - - path: /src/control_methods/utils.py -engines: - - type: docker - image: openproblems/base_python:1 -runners: - - type: executable - - type: nextflow - directives: - label: [lowtime,lowmem,lowcpu] diff --git a/src/control_methods/shuffle_integration_within_cell_type/script.py b/src/control_methods/shuffle_integration_within_cell_type/script.py deleted file mode 100644 index 842595788..000000000 --- a/src/control_methods/shuffle_integration_within_cell_type/script.py +++ /dev/null @@ -1,66 +0,0 @@ -import sys - -import anndata as ad - -## VIASH START -par = { - "input_unintegrated": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad", - "output_integrated_split1": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/control_integrated_split1.h5ad", - "output_integrated_split2": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/control_integrated_split2.h5ad", -} -meta = { - "name": "shuffle_integration_within_cell_type", - "resources_dir": "src/control_methods", -} -## VIASH END - -print("Importing helper functions", flush=True) -sys.path.append(meta["resources_dir"]) -from utils import _randomize_features - -print("Reading and preparing input files", flush=True) -adata = ad.read_h5ad(par["input_unintegrated"]) -adata_split1 = adata[(adata.obs.is_control > 0) | (adata.obs.split == 1)].copy() -adata_split2 = adata[(adata.obs.is_control > 0) | (adata.obs.split == 2)].copy() - -print("Randomise features - split 1", flush=True) -adata_split1.obs["batch_str"] = adata_split1.obs["batch"].astype(str) -integrated = _randomize_features( - adata_split1.layers["preprocessed"], - partition=adata_split1.obs["cell_type"], -) - -# create new anndata -output_split1 = ad.AnnData( - obs=adata_split1.obs[[]], - var=adata_split1.var[[]], - layers={"integrated": integrated}, - uns={ - "dataset_id": adata_split1.uns["dataset_id"], - "method_id": meta["name"], - "parameters": {}, - }, -) - -print("Randomise features - split 2", flush=True) -adata_split2.obs["batch_str"] = adata_split2.obs["batch"].astype(str) -integrated = _randomize_features( - adata_split2.layers["preprocessed"], - partition=adata_split2.obs["cell_type"], -) - -# create new anndata -output_split2 = ad.AnnData( - obs=adata_split2.obs[[]], - var=adata_split2.var[[]], - layers={"integrated": integrated}, - uns={ - "dataset_id": adata_split2.uns["dataset_id"], - "method_id": meta["name"], - "parameters": {}, - }, -) - -print("Write output AnnData to file", flush=True) -output_split1.write_h5ad(par["output_integrated_split1"], compression="gzip") -output_split2.write_h5ad(par["output_integrated_split2"], compression="gzip") diff --git a/src/metrics/bras/config.vsh.yaml b/src/metrics/bras/config.vsh.yaml deleted file mode 100644 index 5e0a0bd4e..000000000 --- a/src/metrics/bras/config.vsh.yaml +++ /dev/null @@ -1,95 +0,0 @@ -__merge__: ../../api/comp_metric.yaml - -# A unique identifier for your component (required). -# Can contain only lowercase letters or underscores. -name: bras -status: disabled -# Metadata for your component -info: - metrics: - # A unique identifier for your metric (required). - # Can contain only lowercase letters or underscores. - - name: bras - # A relatively short label, used when rendering visualisarions (required) - label: Bras - # A one sentence summary of how this metric works (required). Used when - # rendering summary tables. - summary: "Batch removal adapted silhouette (BRAS)." - # A multi-line description of how this component works (required). Used - # when rendering reference documentation. - description: | - BRAS evaluates batch effect removal with respect to batch ids within each label (cell type annotations derived from manual gating), - using a modified silhouette score that accounts for nested batch effects. - Unlike standard silhouette, BRAS computes between-cluster distances using the between_cluster_distances method rather than nearest-cluster approach. - The implementation is based on the scib-metrics package. The final bras score is the average of the bras scores for both data splits. - Cosine similarity is used as a distance measure. - A higher scores indicates better batch mixing. - references: - doi: - - 10.1101/2025.01.21.634098 - bibtex: - - | - @article{rautenstrauch2025metrics, - title={Metrics Matter: Why We Need to Stop Using Silhouette in Single-Cell Benchmarking}, - author={Rautenstrauch, Pia and Ohler, Uwe}, - journal={bioRxiv}, - pages={2025--01}, - year={2025}, - publisher={Cold Spring Harbor Laboratory} - } - links: - # URL to the documentation for this metric (required). - documentation: https://scib-metrics.readthedocs.io/en/stable/generated/scib_metrics.bras.html - # URL to the code repository for this metric (required). - repository: https://github.com/YosefLab/scib-metrics/blob/main/src/scib_metrics/metrics/_silhouette.py - # The minimum possible value for this metric (required) - min: 0 - # The maximum possible value for this metric (required) - max: 1 - # Whether a higher value represents a 'better' solution (required) - maximize: true - -# Resources required to run the component -resources: - # The script of your component (required) - - type: python_script - path: script.py - - path: /src/utils/helper_functions.py - - -engines: - # Specifications for the Docker image for this component. - # testing gpu jax version - # - type: docker - # image: openproblems/base_pytorch_nvidia:1.1 - # setup: - # - type: python - # packages: - # - jax[cuda_12_pip] - # - scib-metrics~=0.5.6 - - type: docker - image: python:3.11 - setup: - - type: apt - packages: - - procps - - type: python - packages: - - jax~=0.6.2 - - jaxlib~=0.6.2 - - anndata~=0.11.0 - - scanpy~=1.11.0 - - scib-metrics~=0.5.6 - - pyyaml - - requests - - jsonschema - github: - - "openproblems-bio/core#subdirectory=packages/python/openproblems" - -runners: - # This platform allows running the component natively - - type: executable - # Allows turning the component into a Nextflow module / pipeline. - - type: nextflow - directives: - label: [midtime,midmem,highcpu] diff --git a/src/metrics/bras/script.py b/src/metrics/bras/script.py deleted file mode 100644 index db3428a13..000000000 --- a/src/metrics/bras/script.py +++ /dev/null @@ -1,87 +0,0 @@ -import sys - -import anndata as ad -import numpy as np -from scib_metrics import bras - -## VIASH START -# The following code has been auto-generated by Viash. -par = { - "input_unintegrated": r"resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad", - "input_integrated_split1": r"resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split1.h5ad", - "input_integrated_split2": r"resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split2.h5ad", - "output": r"resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/harmonypy_bras_score.h5ad", -} -meta = { - "name": r"bras", - "functionality_name": r"bras", - "resources_dir": r"src/utils", -} - -## VIASH END - -sys.path.append(meta["resources_dir"]) -from helper_functions import ( - get_obs_var_for_integrated, - remove_unlabelled, - subset_markers_tocorrect, - subset_nocontrols, -) - -print("Reading input files", flush=True) -integrated_s1 = ad.read_h5ad(par["input_integrated_split1"]) -integrated_s2 = ad.read_h5ad(par["input_integrated_split2"]) -unintegrated = ad.read_h5ad(par["input_unintegrated"]) - -print("Formatting input files", flush=True) -integrated_s1, integrated_s2 = get_obs_var_for_integrated( - integrated_s1, integrated_s2, unintegrated -) - -integrated_s1 = subset_nocontrols(integrated_s1) -integrated_s1 = subset_markers_tocorrect(integrated_s1) -integrated_s1 = subset_nocontrols(integrated_s1) -integrated_s1 = remove_unlabelled(integrated_s1) - -integrated_s2 = subset_nocontrols(integrated_s2) -integrated_s2 = subset_markers_tocorrect(integrated_s2) -integrated_s2 = subset_nocontrols(integrated_s2) -integrated_s2 = remove_unlabelled(integrated_s2) - -print("Compute metrics", flush=True) -batch_labels_s1 = integrated_s1.obs["batch"].values -ct_labels_s1 = integrated_s1.obs["cell_type"].values - -bras_s1 = bras( - integrated_s1.layers["integrated"], - labels=ct_labels_s1, - batch=batch_labels_s1, - metric="euclidean", - chunk_size=512, -) - -batch_labels_s2 = integrated_s2.obs["batch"].values -ct_labels_s2 = integrated_s2.obs["cell_type"].values - -bras_s2 = bras( - integrated_s2.layers["integrated"], - labels=ct_labels_s2, - batch=batch_labels_s2, - metric="euclidean", - chunk_size=512, -) - -bras_score = np.mean([bras_s1, bras_s2]) - -print("Write output AnnData to file", flush=True) -output = ad.AnnData( - uns={ - "dataset_id": integrated_s1.uns["dataset_id"], - "method_id": integrated_s1.uns["method_id"], - "metric_ids": ["bras"], - "metric_values": [bras_score], - "bras_s1": bras_s1, - "bras_s2": bras_s2, - } -) -output.write_h5ad(par["output"], compression="gzip") diff --git a/src/metrics/lisi/config.vsh.yaml b/src/metrics/lisi/config.vsh.yaml index 253f0cf5c..592984858 100644 --- a/src/metrics/lisi/config.vsh.yaml +++ b/src/metrics/lisi/config.vsh.yaml @@ -110,7 +110,7 @@ engines: packages: - jax~=0.6.2 - jaxlib~=0.6.2 - - anndata~=0.11.0 + - anndata~=0.12.0 - scanpy~=1.11.0 - scib-metrics~=0.5.7 - pyyaml diff --git a/src/metrics/n_inconsistent_peaks/config.vsh.yaml b/src/metrics/n_inconsistent_peaks/config.vsh.yaml deleted file mode 100644 index 9c8bae0ec..000000000 --- a/src/metrics/n_inconsistent_peaks/config.vsh.yaml +++ /dev/null @@ -1,93 +0,0 @@ - -__merge__: ../../api/comp_metric.yaml - -name: n_inconsistent_peaks -status: disabled - -info: - metrics: - - name: n_inconsistent_peaks - label: Number of inconsistent peaks Global - # A one sentence summary of how this metric works (required). Used when - # rendering summary tables. - summary: "Comparison of the number of marker‑expression peaks between validation and batch‑normalized data." - # A multi-line description of how this component works (required). Used - # when rendering reference documentation. - description: | - The metric compares the number of marker expression peaks between batch integrated technical replicates (split 1 and split 2). - The metric is calculated as the absolute difference between the number of peaks in the technical replicates. - The marker expression profiles are first z-score scaled (with common mu and sigma) and then smoothed using kernel density estimation (KDE) (`scipy.stats.gaussian_kde`). - Finally, peaks are identified using the `scipy.signal.find_peaks` function. - For peak calling, the `prominence` parameter is set to 0.01 and the `height` parameter is set to 0.1. - references: - doi: - - 10.1038/s41592-019-0686-2 - links: - # URL to the documentation for this metric (required). - documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html#scipy.signal.find_peaks - # URL to the code repository for this metric (required). - repository: https://github.com/scipy/scipy/blob/main/scipy/signal/_peak_finding.py - # The minimum possible value for this metric (required) - min: 0 - # The maximum possible value for this metric (required) - max: +.inf - # Whether a higher value represents a 'better' solution (required) - maximize: false - - - name: n_inconsistent_peaks_ct - label: Number of inconsistent peaks (Cell Type) - summary: "Comparison of the number of cell‑type marker‑expression peaks between validation and batch‑normalized data." - description: | - The metric compares the number of cell type specific marker expression peaks between the validation and batch-normalized data. - The number of peaks is calculated using the `scipy.signal.find_peaks` function. - The metric is calculated as the absolute difference between the number of peaks in the validation and batch-normalized data. - The (cell type) marker expression profiles are first smoothed using kernel density estimation (KDE) (`scipy.stats.gaussian_kde`), - and then peaks are then identified using the `scipy.signal.find_peaks` function. - For peak calling, the `prominence` parameter is set to 0.1 and the `height` parameter is set to 0.05*max_density. - references: - doi: - - 10.1038/s41592-019-0686-2 - links: - # URL to the documentation for this metric (required). - documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html#scipy.signal.find_peaks - # URL to the code repository for this metric (required). - repository: https://github.com/scipy/scipy/blob/v1.15.2/scipy/signal/_peak_finding.py#L0-L1 - # The minimum possible value for this metric (required) - min: 0 - # The maximum possible value for this metric (required) - max: +.inf - # Whether a higher value represents a 'better' solution (required) - maximize: false - -# Component-specific parameters (optional) -# arguments: -# - name: "--n_neighbors" -# type: "integer" -# default: 5 -# description: Number of neighbors to use. - -# Resources required to run the component -resources: - # The script of your component (required) - - type: python_script - path: script.py - - path: helper.py - - path: /src/utils/helper_functions.py - -engines: - # Specifications for the Docker image for this component. - - type: docker - image: openproblems/base_python:1 - # Add custom dependencies here (optional). For more information, see - # https://viash.io/reference/config/engines/docker/#setup . - # setup: - # - type: python - # packages: numpy<2 - -runners: - # This platform allows running the component natively - - type: executable - # Allows turning the component into a Nextflow module / pipeline. - - type: nextflow - directives: - label: [midtime,midmem,midcpu] diff --git a/src/metrics/n_inconsistent_peaks/helper.py b/src/metrics/n_inconsistent_peaks/helper.py deleted file mode 100644 index 9ffb3a579..000000000 --- a/src/metrics/n_inconsistent_peaks/helper.py +++ /dev/null @@ -1,47 +0,0 @@ -import numpy as np -from scipy.stats import gaussian_kde -from scipy.signal import find_peaks - -def get_kde_density(expression_array): - ''' - Returns the density of the array using a gaussian kernel density estimation. - - Inputs: - expression_array: array of values (1D) representing the marker expression - - Outputs: - density: array of values representing the density of marker expression - ''' - - min_val = expression_array.min() - max_val = expression_array.max() - marker_values = np.reshape(expression_array, (1,-1)) # Reshape array for KDE - kde = gaussian_kde(marker_values, bw_method='scott') - x_grid = np.linspace(min_val, max_val, 100) - density = kde(x_grid) - #Plot, for debugging - # sns.scatterplot(x=x_grid, y=density) - # plt.show() - return density - -def call_peaks(density): - ''' - Returns the peaks of the density using scipy.signal.find_peaks. - - Inputs: - density: array of values representing the density of marker expression - - Outputs: - peaks: array of values representing the peaks of the density - ''' - - height_trsh = 0.1 - prom_trsh = 0.01 - - peaks, _ = find_peaks(density, - prominence=prom_trsh, - height= height_trsh - ) - num_peaks = len(peaks) - - return num_peaks \ No newline at end of file diff --git a/src/metrics/n_inconsistent_peaks/script.py b/src/metrics/n_inconsistent_peaks/script.py deleted file mode 100644 index 6e1c858c2..000000000 --- a/src/metrics/n_inconsistent_peaks/script.py +++ /dev/null @@ -1,124 +0,0 @@ -import sys - -import anndata as ad -import numpy as np - -## VIASH START -# The following code has been auto-generated by Viash. -par = { - 'input_unintegrated': 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad', - 'input_integrated_split1': 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split1.h5ad', - 'input_integrated_split2': 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split2.h5ad', - 'output': 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/score.h5ad' -} -meta = { - 'name': 'n_inconsistent_peaks', -} - - -## VIASH END - -sys.path.append(meta["resources_dir"]) -from helper import call_peaks, get_kde_density -from helper_functions import ( - get_obs_var_for_integrated, - remove_unlabelled, - subset_markers_tocorrect, - subset_nocontrols, -) - -print("Reading input files", flush=True) -integrated_s1 = ad.read_h5ad(par["input_integrated_split1"]) -integrated_s2 = ad.read_h5ad(par["input_integrated_split2"]) -unintegrated = ad.read_h5ad(par["input_unintegrated"]) - -print("Formatting input files", flush=True) -integrated_s1, integrated_s2 = get_obs_var_for_integrated( - integrated_s1, integrated_s2, unintegrated -) - -integrated_s1 = subset_nocontrols(integrated_s1) -integrated_s1 = subset_markers_tocorrect(integrated_s1) -integrated_s1 = subset_nocontrols(integrated_s1) -integrated_s1 = remove_unlabelled(integrated_s1) - -integrated_s2 = subset_nocontrols(integrated_s2) -integrated_s2 = subset_markers_tocorrect(integrated_s2) -integrated_s2 = subset_nocontrols(integrated_s2) -integrated_s2 = remove_unlabelled(integrated_s2) - -print('Compute metric (All cells)', flush=True) -donor_list = integrated_s1.obs['donor'].unique() -n_inconsistent_peaks = 0 - -for donor in donor_list: - s1_view = integrated_s1[integrated_s1.obs['donor'] == donor] - s2_view = integrated_s2[integrated_s2.obs['donor'] == donor] - - for marker in s1_view.var.index: - marker_expression_s1_unscaled = np.array(s1_view[:,marker].layers["integrated"]) - marker_expression_s2_unscaled = np.array(s2_view[:,marker].layers["integrated"]) - - pooled = np.concatenate([marker_expression_s1_unscaled, marker_expression_s2_unscaled]) - mu, sd = pooled.mean(), pooled.std() - marker_expression_s1 = (marker_expression_s1_unscaled - mu) / (sd) - marker_expression_s2 = (marker_expression_s2_unscaled - mu) / (sd) - - density_s1 = get_kde_density(marker_expression_s1) - peaks_s1 = call_peaks(density_s1) - density_s2 = get_kde_density(marker_expression_s2) - peaks_s2 = call_peaks(density_s2) - - if peaks_s1 != peaks_s2: - n_inconsistent_peaks += abs(peaks_s1 - peaks_s2) - -print('Compute metric (per cell type)', flush=True) -n_inconsistent_peaks_ct = 0 - -for donor in donor_list: - s1_view = integrated_s1[integrated_s1.obs['donor'] == donor] - s2_view = integrated_s2[integrated_s2.obs['donor'] == donor] - celltype_list = s1_view.obs['cell_type'].unique() - - for celltype in celltype_list: - s1_view_ct = s1_view[s1_view.obs['cell_type'] == celltype] - s2_view_ct = s2_view[s2_view.obs['cell_type'] == celltype] - - if s1_view_ct.shape[0] < 100 or s2_view_ct.shape[0] < 100: - print(donor,celltype,"skipped cause less than 100 cells are present in either split 1 or split 2 dataset") - continue - - for marker in s1_view_ct.var.index: - marker_expression_s1_unscaled = np.array(s1_view_ct[:, marker].layers["integrated"]) - marker_expression_s2_unscaled = np.array(s2_view_ct[:, marker].layers["integrated"]) - - pooled = np.concatenate([marker_expression_s1_unscaled, marker_expression_s2_unscaled]) - mu, sd = pooled.mean(), pooled.std() - marker_expression_s1 = (marker_expression_s1_unscaled - mu) / (sd) - marker_expression_s2 = (marker_expression_s2_unscaled - mu) / (sd) - - density_s1 = get_kde_density(marker_expression_s1) - peaks_s1 = call_peaks(density_s1) - density_s2 = get_kde_density(marker_expression_s2) - peaks_s2 = call_peaks(density_s2) - - if peaks_s1 != peaks_s2: - n_inconsistent_peaks_ct += abs(peaks_s1 - peaks_s2) - - - -uns_metric_ids = ['n_inconsistent_peaks', 'n_inconsistent_peaks_ct'] -uns_metric_values = [ float(n_inconsistent_peaks), float(n_inconsistent_peaks_ct) ] - -print("Write output AnnData to file", flush=True) -output = ad.AnnData( - uns={ - 'dataset_id': integrated_s1.uns['dataset_id'], - 'method_id': integrated_s1.uns['method_id'], - 'metric_ids': uns_metric_ids, - 'metric_values': uns_metric_values - } -) -output.write_h5ad(par['output'], compression='gzip') - -print(uns_metric_ids, uns_metric_values) \ No newline at end of file diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml index b1f99d7d9..82b4ac030 100644 --- a/src/workflows/run_benchmark/config.vsh.yaml +++ b/src/workflows/run_benchmark/config.vsh.yaml @@ -83,8 +83,6 @@ dependencies: - name: utils/extract_uns_metadata repository: op - name: control_methods/shuffle_integration_globally - - name: control_methods/shuffle_integration_within_batch - - name: control_methods/shuffle_integration_within_cell_type - name: control_methods/no_integration - name: control_methods/perfect_integration - name: methods/harmonypy diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf index dee4ca18f..c30f82863 100644 --- a/src/workflows/run_benchmark/main.nf +++ b/src/workflows/run_benchmark/main.nf @@ -10,8 +10,6 @@ workflow auto { // construct list of methods and control methods methods = [ shuffle_integration_globally, - shuffle_integration_within_batch, - shuffle_integration_within_cell_type, no_integration, perfect_integration, combat,