diff --git a/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/diff_expr.py b/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/diff_expr.py
index f6d3c5a16..05234762c 100644
--- a/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/diff_expr.py
+++ b/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/diff_expr.py
@@ -35,9 +35,11 @@
]
+# flake8: noqa: D103
+
+
def compute_memento_estimators_from_precomputed_stats(estimators_df: pl.DataFrame) -> pl.DataFrame:
- """
- Computes the mean and standard error of the mean (SEM) for each feature in the estimators DataFrame.
+ """Computes the mean and standard error of the mean (SEM) for each feature in the estimators DataFrame.
This function takes a DataFrame containing precomputed statistics for each feature, including the number of observations,
sum, sum of squares, and size factor. It calculates the mean and SEM for each feature based on these statistics.
@@ -112,16 +114,20 @@ def compute_all(
n_processes: int,
covariates_str: Optional[str] = None,
) -> Tuple[pd.DataFrame, pstats.Stats]:
- default_covariates = CUBE_LOGICAL_DIMS_OBS
+ default_covariates: List[str] = []
if covariates_str is None:
covariates = default_covariates
else:
covariates = covariates_str.split(",")
+
+ # make treatment variable be in the first column because that is needed for the design matrix
+ # NOTE: if covariates == [] then variables for the design matrix will only contain the treatment column
+ variables = [treatment] + [covariate for covariate in covariates if covariate != treatment]
+
with tiledb.open(os.path.join(cube_path, OBS_GROUPS_ARRAY), "r") as obs_groups_array:
obs_groups_df = obs_groups_array.query(cond=query_filter or None).df[:]
- if covariates != default_covariates:
- obs_groups_df = obs_groups_df[covariates + [treatment, "obs_group_joinid", "n_obs"]]
+ obs_groups_df = obs_groups_df[variables + ["obs_group_joinid", "n_obs"]]
distinct_treatment_values = obs_groups_df[treatment].nunique()
assert distinct_treatment_values == 2, "treatment must have exactly 2 distinct values"
@@ -135,8 +141,6 @@ def compute_all(
f"computing for {len(obs_groups_df)} obs groups ({obs_groups_df.n_obs.sum()} cells) and {len(features)} features using {n_feature_groups} processes, {len(features) // n_feature_groups} features/process"
)
- # make treatment variable be in the first column of the design matrix
- variables = [treatment] + [covariate for covariate in covariates if covariate != treatment]
selected_vars_groups_groupby = obs_groups_df.groupby(variables, observed=True)
agg_dict = {i: "first" for i in variables}
@@ -308,9 +312,7 @@ def de_wls(
n: npt.NDArray[np.float32],
v: npt.NDArray[np.float32],
) -> Tuple[np.float32, np.float32, np.float32]:
- """
- Perform DE for each gene using Weighted Least Squares (i.e., a weighted Linear Regression model)
- """
+ """Perform DE for each gene using Weighted Least Squares (i.e., a weighted Linear Regression model)."""
coef = de_wls_fit(X, y, n)
z, pv = de_wls_stats(X, v, coef)
diff --git a/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/scientific_validation_notebooks/memento_vanilla_testcases.ipynb b/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/scientific_validation_notebooks/memento_vanilla_testcases.ipynb
new file mode 100644
index 000000000..dd55f5a67
--- /dev/null
+++ b/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/scientific_validation_notebooks/memento_vanilla_testcases.ipynb
@@ -0,0 +1,883 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "9670cc2f-903a-4753-9a12-123c80859f77",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "import memento\n",
+ "import pandas as pd\n",
+ "\n",
+ "import cellxgene_census"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "73efecd1-faaa-48f0-b230-d460f5b745a5",
+ "metadata": {},
+ "outputs": [
+ {
+ "ename": "ModuleNotFoundError",
+ "evalue": "No module named 'itables'",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)",
+ "Cell \u001b[0;32mIn[3], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mitables\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m init_notebook_mode\n\u001b[1;32m 2\u001b[0m init_notebook_mode(all_interactive\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n",
+ "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'itables'"
+ ]
+ }
+ ],
+ "source": [
+ "from itables import init_notebook_mode\n",
+ "\n",
+ "init_notebook_mode(all_interactive=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "ef602f2a-7c34-4729-bde9-9588b3bee9d3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "census = cellxgene_census.open_soma(census_version=\"2023-12-15\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "c8506c34-2e45-4f06-b19c-f4f4cbae12c7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# guo, disease == \"COVID-19\"\n",
+ "dataset_guo = \"ae5341b8-60fb-4fac-86db-86e49ee66287\"\n",
+ "# aruna\n",
+ "dataset_aruna = \"59b69042-47c2-47fd-ad03-d21beb99818f\"\n",
+ "# wilk, slide-seq\n",
+ "dataset_wilk = \"055ca631-6ffb-40de-815e-b931e10718c0\"\n",
+ "\n",
+ "cell_type_monocyte = \"classical monocyte\"\n",
+ "cell_type_t_cell = \"central memory CD4-positive, alpha-beta T cell\"\n",
+ "cell_type_nk = \"natural killer cell\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fffe0e77-112e-419f-b592-4a6c0b825e22",
+ "metadata": {},
+ "source": [
+ "## Perform 5-10 more concordance comparisons with a controlled environment \n",
+ "\n",
+ "Steps:\n",
+ "- Look for example cases\n",
+ " - 1-2 known comparisons, single-dataset single-donor, no batch effects.\n",
+ " - 2-4 known comparisons, simple batch-effects with multi-donors but no multi-assay nor multi-dataset.\n",
+ " - 2-4 known comparisons, complex batch-effects with multi-datasets and multi-assay.\n",
+ "- For each example\n",
+ " - Run vanilla memento\n",
+ " - Run pre-comuputed memento\n",
+ " - Perform comparison of p-values between the two versions as in plot below\n",
+ "- Validation passes if there is linear correlation, with spearman > 0.9\n",
+ "- Brownie points if there is a high overlap between between significant genes at p < 0.001, do hypergeometric test\n",
+ "\n",
+ "### Example cases\n",
+ "\n",
+ "#### 1-2 known comparisons, single-dataset single-donor, no batch effects.\n",
+ "\n",
+ "**a) Single-dataset, single-donor, no batch effects.**\n",
+ "\n",
+ "- Collection: A Web Portal and Workbench for Biological Dissection of Single Cell COVID-19 Host Responses\n",
+ "- Dataset: Individual Single-Cell RNA-seq PBMC Data from Arunachalam et al.\n",
+ " - Assay: 10X\n",
+ "- Comparison:\n",
+ " - classical monocytes vs T-cells in one donor\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "0286a638-9b78-41c6-8e10-b6dba69558d2",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "datasets = [dataset_aruna]\n",
+ "donor_id = \"cov17\"\n",
+ "disease = \"normal\"\n",
+ "cell_types = [cell_type_monocyte, cell_type_t_cell]\n",
+ "\n",
+ "adata = cellxgene_census.get_anndata(\n",
+ " census=census,\n",
+ " organism=\"homo_sapiens\",\n",
+ " obsm_layers=[\"scvi\"],\n",
+ " obs_value_filter=f\"dataset_id in {datasets} and disease == '{disease}' and cell_type in {cell_types} and donor_id == '{donor_id}'\",\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "561ca819-07d9-4b72-8e82-5940f6e27c0c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# adata.obs[\"donor_id\"].value_counts()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "5b0727bf-fb32-403b-b674-9f8398a784fb",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "# scanpy.pp.neighbors(adata, use_rep=\"scvi\")\n",
+ "# scanpy.tl.umap(adata)\n",
+ "# scanpy.pl.umap(adata, color=[\"cell_type\"])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5bf9391a-2320-4a89-b0f3-4942ca18b7d0",
+ "metadata": {},
+ "source": [
+ "Running memento vanilla"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "9d61902a-2d06-4c6d-97c9-c541d62d577b",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/psridharan/miniconda3/envs/vanilla-memento/lib/python3.10/site-packages/memento/main.py:181: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_numeric without passing `errors` and catch exceptions explicitly instead\n",
+ " df[col] = pd.to_numeric(df[col], errors='ignore')\n",
+ "/Users/psridharan/miniconda3/envs/vanilla-memento/lib/python3.10/site-packages/tiledb/ctx.py:560: UserWarning: TileDB is a multithreading library and deadlocks are likely if fork() is called after a TileDB context has been created (such as for array access). To safely use TileDB with multiprocessing or concurrent.futures, choose 'spawn' as the start method for child processes. For example: multiprocessing.set_start_method('spawn').\n",
+ " warnings.warn(\n",
+ "[Parallel(n_jobs=12)]: Using backend LokyBackend with 12 concurrent workers.\n",
+ "[Parallel(n_jobs=12)]: Done 26 tasks | elapsed: 3.0s\n",
+ "[Parallel(n_jobs=12)]: Done 176 tasks | elapsed: 4.7s\n",
+ "[Parallel(n_jobs=12)]: Done 576 tasks | elapsed: 8.7s\n",
+ "/Users/psridharan/miniconda3/envs/vanilla-memento/lib/python3.10/site-packages/joblib/externals/loky/process_executor.py:752: UserWarning: A worker stopped while some jobs were given to the executor. This can be caused by a too short worker timeout or by a memory leak.\n",
+ " warnings.warn(\n",
+ "[Parallel(n_jobs=12)]: Done 1276 tasks | elapsed: 17.8s\n",
+ "[Parallel(n_jobs=12)]: Done 2176 tasks | elapsed: 29.1s\n",
+ "[Parallel(n_jobs=12)]: Done 3276 tasks | elapsed: 41.6s\n",
+ "[Parallel(n_jobs=12)]: Done 4264 tasks | elapsed: 56.4s\n",
+ "[Parallel(n_jobs=12)]: Done 5312 out of 5335 | elapsed: 1.1min remaining: 0.3s\n",
+ "[Parallel(n_jobs=12)]: Done 5335 out of 5335 | elapsed: 1.2min finished\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Setup\n",
+ "\n",
+ "# Assuming sequenced to 50%, dataset specific number\n",
+ "adata.var.index = adata.var[\"feature_name\"].values\n",
+ "adata.obs[\"q\"] = 0.15\n",
+ "\n",
+ "\n",
+ "# Classical monocyte encoded as 1\n",
+ "adata.obs[\"treatment\"] = (adata.obs[\"cell_type\"] == cell_types[0]).astype(int)\n",
+ "\n",
+ "# Setup memento\n",
+ "memento.setup_memento(adata, q_column=\"q\", trim_percent=0.1) # trim_percent tunes cell size calculation\n",
+ "memento.create_groups(adata, label_columns=[\"treatment\"])\n",
+ "memento.compute_1d_moments(adata, min_perc_group=0.9)\n",
+ "group_metadata = memento.get_groups(adata)\n",
+ "\n",
+ "treatment_df = group_metadata[[\"treatment\"]]\n",
+ "\n",
+ "memento.ht_1d_moments(\n",
+ " adata,\n",
+ " # covariate=covariate_df,\n",
+ " treatment=treatment_df,\n",
+ " num_boot=5000,\n",
+ " verbose=1,\n",
+ " num_cpus=12,\n",
+ " resample_rep=False,\n",
+ " approx=False,\n",
+ ")\n",
+ "\n",
+ "result = memento.get_1d_ht_result(adata)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "bd72f75b-a699-4510-9b9f-9fe0f9706052",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# result"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "bf6027c5-e344-43f2-a18c-f22a0e731f85",
+ "metadata": {},
+ "source": [
+ "**b) Single-dataset, single-donor, no batch effects.**\n",
+ "\n",
+ "- Collection: A Web Portal and Workbench for Biological Dissection of Single Cell COVID-19 Host Responses\n",
+ "- Dataset: Individual Single-Cell RNA-seq PBMC Data from Arunachalam et al.\n",
+ " - Assay: 10X\n",
+ "- Comparison:\n",
+ " - classical monocytes vs natural killer cells in one donor\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "7190dd91-f4ea-40eb-ab85-42d974c3eda1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "datasets = [dataset_aruna]\n",
+ "donor_id = \"cov17\"\n",
+ "disease = \"normal\"\n",
+ "cell_types = [cell_type_monocyte, cell_type_nk]\n",
+ "\n",
+ "adata = cellxgene_census.get_anndata(\n",
+ " census=census,\n",
+ " organism=\"homo_sapiens\",\n",
+ " obsm_layers=[\"scvi\"],\n",
+ " obs_value_filter=f\"dataset_id in {datasets} and disease == '{disease}' and cell_type in {cell_types} and donor_id == '{donor_id}'\",\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "id": "68146cff-495f-4788-b731-ad12c3b024f9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# adata.obs[\"donor_id\"].value_counts()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "id": "aaf6cd74-79b9-47a3-8bc8-aeab0b40011f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# scanpy.pp.neighbors(adata, use_rep=\"scvi\")\n",
+ "# scanpy.tl.umap(adata)\n",
+ "# scanpy.pl.umap(adata, color=[\"cell_type\"])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "d85e27d8-e85c-4ea7-89c4-4ff6211207f6",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/psridharan/miniconda3/envs/vanilla-memento/lib/python3.10/site-packages/memento/main.py:181: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_numeric without passing `errors` and catch exceptions explicitly instead\n",
+ " df[col] = pd.to_numeric(df[col], errors='ignore')\n",
+ "[Parallel(n_jobs=12)]: Using backend LokyBackend with 12 concurrent workers.\n",
+ "/Users/psridharan/miniconda3/envs/vanilla-memento/lib/python3.10/site-packages/tiledb/ctx.py:560: UserWarning: TileDB is a multithreading library and deadlocks are likely if fork() is called after a TileDB context has been created (such as for array access). To safely use TileDB with multiprocessing or concurrent.futures, choose 'spawn' as the start method for child processes. For example: multiprocessing.set_start_method('spawn').\n",
+ " warnings.warn(\n",
+ "[Parallel(n_jobs=12)]: Done 26 tasks | elapsed: 2.9s\n",
+ "[Parallel(n_jobs=12)]: Done 176 tasks | elapsed: 4.8s\n",
+ "/Users/psridharan/miniconda3/envs/vanilla-memento/lib/python3.10/site-packages/joblib/externals/loky/process_executor.py:752: UserWarning: A worker stopped while some jobs were given to the executor. This can be caused by a too short worker timeout or by a memory leak.\n",
+ " warnings.warn(\n",
+ "[Parallel(n_jobs=12)]: Done 426 tasks | elapsed: 8.0s\n",
+ "[Parallel(n_jobs=12)]: Done 1108 tasks | elapsed: 18.2s\n",
+ "[Parallel(n_jobs=12)]: Done 2008 tasks | elapsed: 31.4s\n",
+ "[Parallel(n_jobs=12)]: Done 3108 tasks | elapsed: 46.8s\n",
+ "[Parallel(n_jobs=12)]: Done 4300 tasks | elapsed: 1.1min\n",
+ "[Parallel(n_jobs=12)]: Done 5224 out of 5247 | elapsed: 1.3min remaining: 0.4s\n",
+ "[Parallel(n_jobs=12)]: Done 5247 out of 5247 | elapsed: 1.3min finished\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Setup\n",
+ "\n",
+ "# Assuming sequenced to 50%, dataset specific number\n",
+ "adata.var.index = adata.var[\"feature_name\"].values\n",
+ "adata.obs[\"q\"] = 0.15\n",
+ "\n",
+ "\n",
+ "# Classical monocyte encoded as 1\n",
+ "adata.obs[\"treatment\"] = (adata.obs[\"cell_type\"] == cell_types[0]).astype(int)\n",
+ "\n",
+ "# Setup memento\n",
+ "memento.setup_memento(adata, q_column=\"q\", trim_percent=0.1) # trim_percent tunes cell size calculation\n",
+ "memento.create_groups(adata, label_columns=[\"treatment\"])\n",
+ "memento.compute_1d_moments(adata, min_perc_group=0.9)\n",
+ "group_metadata = memento.get_groups(adata)\n",
+ "\n",
+ "treatment_df = group_metadata[[\"treatment\"]]\n",
+ "\n",
+ "memento.ht_1d_moments(\n",
+ " adata,\n",
+ " # covariate=covariate_df,\n",
+ " treatment=treatment_df,\n",
+ " num_boot=5000,\n",
+ " verbose=1,\n",
+ " num_cpus=12,\n",
+ " resample_rep=False,\n",
+ " approx=False,\n",
+ ")\n",
+ "\n",
+ "result = memento.get_1d_ht_result(adata)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "id": "065daced-0e91-428a-802c-06d2cdcb518b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# result"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f9b27000-b471-4aa1-adb1-9039e770b918",
+ "metadata": {},
+ "source": [
+ "#### 2-4 known comparisons, simple batch-effects with multi-datasets but no multi-assay.\n",
+ "\n",
+ "**a) Single-dataset, two donors, controlling for donors**\n",
+ "\n",
+ "- Collection: A Web Portal and Workbench for Biological Dissection of Single Cell COVID-19 Host Responses\n",
+ "- Dataset: Individual Single-Cell RNA-seq PBMC Data from Arunachalam et al.\n",
+ " - Assay: 10X\n",
+ "- Comparison:\n",
+ " - classical monocytes vs T-cells in two donors, including donors as covariates.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "id": "2b315f72-c8f5-48f7-b906-5d733d10bf92",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "datasets = [dataset_aruna]\n",
+ "donor_id = [\"cov17\", \"cov18\"]\n",
+ "disease = \"normal\"\n",
+ "cell_types = [cell_type_monocyte, cell_type_nk]\n",
+ "\n",
+ "adata = cellxgene_census.get_anndata(\n",
+ " census=census,\n",
+ " organism=\"homo_sapiens\",\n",
+ " obsm_layers=[\"scvi\"],\n",
+ " obs_value_filter=f\"dataset_id in {datasets} and disease == '{disease}' and cell_type in {cell_types} and donor_id in {donor_id}\",\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "id": "10e1eafd-866c-467f-9ebc-b7184875044d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# adata.obs[\"donor_id\"].value_counts()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "id": "471dc42e-6d0e-4659-ba99-ce3b8eb0767e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# scanpy.pp.neighbors(adata, use_rep=\"scvi\")\n",
+ "# scanpy.tl.umap(adata)\n",
+ "# scanpy.pl.umap(adata, color=[\"cell_type\"])\n",
+ "# scanpy.pl.umap(adata, color=[\"donor_id\"])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "id": "54d6c680-4288-418c-8d7c-be4384bdee61",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/pgarcia-nieto/scripts/repos/scrna-parameter-estimation/memento/main.py:181: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_numeric without passing `errors` and catch exceptions explicitly instead\n",
+ " df[col] = pd.to_numeric(df[col], errors='ignore')\n",
+ "[Parallel(n_jobs=12)]: Using backend LokyBackend with 12 concurrent workers.\n",
+ "[Parallel(n_jobs=12)]: Done 29 tasks | elapsed: 0.7s\n",
+ "[Parallel(n_jobs=12)]: Done 328 tasks | elapsed: 5.4s\n",
+ "/Users/pgarcia-nieto/venv_uce/lib/python3.11/site-packages/joblib/externals/loky/process_executor.py:752: UserWarning: A worker stopped while some jobs were given to the executor. This can be caused by a too short worker timeout or by a memory leak.\n",
+ " warnings.warn(\n",
+ "/Users/pgarcia-nieto/venv_uce/lib/python3.11/site-packages/tiledb/ctx.py:560: UserWarning: TileDB is a multithreading library and deadlocks are likely if fork() is called after a TileDB context has been created (such as for array access). To safely use TileDB with multiprocessing or concurrent.futures, choose 'spawn' as the start method for child processes. For example: multiprocessing.set_start_method('spawn').\n",
+ " warnings.warn(\n",
+ "[Parallel(n_jobs=12)]: Done 828 tasks | elapsed: 14.1s\n",
+ "[Parallel(n_jobs=12)]: Done 1528 tasks | elapsed: 26.8s\n",
+ "[Parallel(n_jobs=12)]: Done 2428 tasks | elapsed: 42.1s\n",
+ "[Parallel(n_jobs=12)]: Done 3384 tasks | elapsed: 1.0min\n",
+ "[Parallel(n_jobs=12)]: Done 4034 tasks | elapsed: 1.2min\n",
+ "[Parallel(n_jobs=12)]: Done 4784 tasks | elapsed: 1.4min\n",
+ "[Parallel(n_jobs=12)]: Done 4834 out of 4857 | elapsed: 1.4min remaining: 0.4s\n",
+ "[Parallel(n_jobs=12)]: Done 4857 out of 4857 | elapsed: 1.4min finished\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Setup\n",
+ "\n",
+ "# Assuming sequenced to 50%, dataset specific number\n",
+ "adata.var.index = adata.var[\"feature_name\"].values\n",
+ "adata.obs[\"q\"] = 0.15\n",
+ "\n",
+ "\n",
+ "# Classical monocyte encoded as 1\n",
+ "adata.obs[\"treatment\"] = (adata.obs[\"cell_type\"] == cell_types[0]).astype(int)\n",
+ "\n",
+ "# Setup memento\n",
+ "memento.setup_memento(adata, q_column=\"q\", trim_percent=0.1) # trim_percent tunes cell size calculation\n",
+ "memento.create_groups(adata, label_columns=[\"treatment\", \"donor_id\"])\n",
+ "memento.compute_1d_moments(adata, min_perc_group=0.9)\n",
+ "group_metadata = memento.get_groups(adata)\n",
+ "\n",
+ "treatment_df = group_metadata[[\"treatment\"]]\n",
+ "covariate_df = pd.get_dummies(group_metadata[[\"donor_id\"]], drop_first=True).astype(float)\n",
+ "covariate_df -= covariate_df.mean() # covariates for Lin's estimator\n",
+ "\n",
+ "# Include interactions\n",
+ "covariate_df[\"interaction\"] = treatment_df.iloc[:, 0] * covariate_df.iloc[:, 0]\n",
+ "\n",
+ "memento.ht_1d_moments(\n",
+ " adata,\n",
+ " covariate=covariate_df,\n",
+ " treatment=treatment_df,\n",
+ " num_boot=5000,\n",
+ " verbose=1,\n",
+ " num_cpus=12,\n",
+ " resample_rep=False,\n",
+ " approx=False,\n",
+ ")\n",
+ "\n",
+ "result = memento.get_1d_ht_result(adata)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "id": "952d9f79-a060-46b1-aaad-719b94557696",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ "\n",
+ " \n",
+ " \n",
+ " | gene | \n",
+ " tx | \n",
+ " de_coef | \n",
+ " de_se | \n",
+ " de_pval | \n",
+ " dv_coef | \n",
+ " dv_se | \n",
+ " dv_pval | \n",
+ "
\n",
+ " | Loading... (need help?) |
\n",
+ "\n",
+ "
\n",
+ "\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " treatment | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | sg^0^cov07 | \n",
+ " 0 | \n",
+ "
\n",
+ " \n",
+ " | sg^1^cov07 | \n",
+ " 1 | \n",
+ "
\n",
+ " \n",
+ " | sg^0^cov08 | \n",
+ " 0 | \n",
+ "
\n",
+ " \n",
+ " | sg^1^cov08 | \n",
+ " 1 | \n",
+ "
\n",
+ " \n",
+ " | sg^0^cov09 | \n",
+ " 0 | \n",
+ "
\n",
+ " \n",
+ " | sg^1^cov09 | \n",
+ " 1 | \n",
+ "
\n",
+ " \n",
+ " | sg^1^cov17 | \n",
+ " 1 | \n",
+ "
\n",
+ " \n",
+ " | sg^0^cov17 | \n",
+ " 0 | \n",
+ "
\n",
+ " \n",
+ " | sg^0^cov18 | \n",
+ " 0 | \n",
+ "
\n",
+ " \n",
+ " | sg^1^cov18 | \n",
+ " 1 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " donor_id_cov08 | \n",
+ " donor_id_cov09 | \n",
+ " donor_id_cov17 | \n",
+ " donor_id_cov18 | \n",
+ " interaction | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | sg^0^cov07 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.0 | \n",
+ "
\n",
+ " \n",
+ " | sg^1^cov07 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ "
\n",
+ " \n",
+ " | sg^0^cov08 | \n",
+ " 0.8 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | sg^1^cov08 | \n",
+ " 0.8 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " 0.8 | \n",
+ "
\n",
+ " \n",
+ " | sg^0^cov09 | \n",
+ " -0.2 | \n",
+ " 0.8 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.0 | \n",
+ "
\n",
+ " \n",
+ " | sg^1^cov09 | \n",
+ " -0.2 | \n",
+ " 0.8 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ "
\n",
+ " \n",
+ " | sg^1^cov17 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " 0.8 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ "
\n",
+ " \n",
+ " | sg^0^cov17 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " 0.8 | \n",
+ " -0.2 | \n",
+ " -0.0 | \n",
+ "
\n",
+ " \n",
+ " | sg^0^cov18 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " 0.8 | \n",
+ " -0.0 | \n",
+ "
\n",
+ " \n",
+ " | sg^1^cov18 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " -0.2 | \n",
+ " 0.8 | \n",
+ " -0.2 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " gene | \n",
+ " tx | \n",
+ " de_coef | \n",
+ " de_se | \n",
+ " de_pval | \n",
+ " dv_coef | \n",
+ " dv_se | \n",
+ " dv_pval | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 937 | \n",
+ " ENSG00000000419 | \n",
+ " treatment | \n",
+ " 0.364687 | \n",
+ " 0.072641 | \n",
+ " 0.000086 | \n",
+ " -0.078413 | \n",
+ " 0.299268 | \n",
+ " 0.792242 | \n",
+ "
\n",
+ " \n",
+ " | 1668 | \n",
+ " ENSG00000001497 | \n",
+ " treatment | \n",
+ " 0.039923 | \n",
+ " 0.077712 | \n",
+ " 0.607279 | \n",
+ " 0.199968 | \n",
+ " 0.302095 | \n",
+ " 0.498100 | \n",
+ "
\n",
+ " \n",
+ " | 104 | \n",
+ " ENSG00000001629 | \n",
+ " treatment | \n",
+ " 0.566675 | \n",
+ " 0.073385 | \n",
+ " 0.000075 | \n",
+ " -0.052604 | \n",
+ " 0.231213 | \n",
+ " 0.819436 | \n",
+ "
\n",
+ " \n",
+ " | 1650 | \n",
+ " ENSG00000001631 | \n",
+ " treatment | \n",
+ " 0.174073 | \n",
+ " 0.070125 | \n",
+ " 0.013597 | \n",
+ " -0.138273 | \n",
+ " 0.276711 | \n",
+ " 0.593681 | \n",
+ "
\n",
+ " \n",
+ " | 1662 | \n",
+ " ENSG00000002549 | \n",
+ " treatment | \n",
+ " 2.020447 | \n",
+ " 0.051671 | \n",
+ " 0.000064 | \n",
+ " 1.234160 | \n",
+ " 0.233148 | \n",
+ " 0.001589 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " coef | \n",
+ " z | \n",
+ " pval | \n",
+ "
\n",
+ " \n",
+ " | feature_id | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | (ENSG00000000003,) | \n",
+ " 1.891206 | \n",
+ " 15.524457 | \n",
+ " 2.370035e-54 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000419,) | \n",
+ " -0.550701 | \n",
+ " -9.138716 | \n",
+ " 6.319813e-20 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000457,) | \n",
+ " -0.029194 | \n",
+ " -0.281349 | \n",
+ " 7.784424e-01 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000460,) | \n",
+ " -0.235270 | \n",
+ " -1.346930 | \n",
+ " 1.780027e-01 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000938,) | \n",
+ " -4.475637 | \n",
+ " -40.313418 | \n",
+ " 0.000000e+00 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " coef | \n",
+ " z | \n",
+ " pval | \n",
+ "
\n",
+ " \n",
+ " | feature_id | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | (ENSG00000000419,) | \n",
+ " -0.550701 | \n",
+ " -9.138716 | \n",
+ " 6.319813e-20 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000001497,) | \n",
+ " -0.266663 | \n",
+ " -4.249166 | \n",
+ " 2.145684e-05 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000001629,) | \n",
+ " -0.769980 | \n",
+ " -12.426287 | \n",
+ " 1.881732e-35 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000001631,) | \n",
+ " -0.377861 | \n",
+ " -6.240487 | \n",
+ " 4.362094e-10 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000002549,) | \n",
+ " -2.272766 | \n",
+ " -50.409757 | \n",
+ " 0.000000e+00 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x = -1 * np.log10(vanilla_memento_result_df[\"de_pval\"])\n",
+ "y = -1 * np.log10(filtered_precomputed_memento_result_df[\"pval\"])\n",
+ "\n",
+ "plt.scatter(x, y, label=\"p values\")\n",
+ "\n",
+ "plt.plot(x, x, color=\"black\", label=\"Line with Slope 1\")\n",
+ "\n",
+ "plt.xlabel(\"vanilla\")\n",
+ "plt.ylabel(\"precomputed\")\n",
+ "\n",
+ "plt.title(\"Scatter Plot of p values\")\n",
+ "plt.legend()\n",
+ "\n",
+ "plt.grid(True)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "937e5d6c-4c88-4e47-a247-5037792ecf56",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.13"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/scientific_validation_notebooks/precomputed_memento_validation_testcase2.ipynb b/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/scientific_validation_notebooks/precomputed_memento_validation_testcase2.ipynb
new file mode 100644
index 000000000..09e549169
--- /dev/null
+++ b/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/scientific_validation_notebooks/precomputed_memento_validation_testcase2.ipynb
@@ -0,0 +1,714 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "6caacabf-1002-47fb-b38d-96fcb4f10133",
+ "metadata": {},
+ "source": [
+ "## Scientific validation of memento implementations\n",
+ "\n",
+ "This notebook runs through test cases for differential expression where the p-value outputs\n",
+ "of vanilla memento implementation are compared to the p-value outputs of precomputed memento\n",
+ "implementation\n",
+ "\n",
+ "**NOTE**:\n",
+ "- Run this notebook from the `cellxgene-census` repo in the `psridharan/memento-sci-val-1` branch\n",
+ "- Download estimator cube from `s3://psridharan-tmp/memento/memento-cube-census-data-2023-12-15/`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 31,
+ "id": "ae923ea5-7c72-4743-aea1-c854593fc1ec",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "import memento\n",
+ "import numpy as np\n",
+ "\n",
+ "import cellxgene_census\n",
+ "import cellxgene_census.experimental.diffexp.memento.diff_expr as precomputed_memento\n",
+ "\n",
+ "census = cellxgene_census.open_soma(census_version=\"2023-12-15\")\n",
+ "precomputed_memento_estimator_cube_path = (\n",
+ " \"/Users/psridharan/code/cellxgene-census/memento-cubes/memento-cube-census-data-2023-12-15\"\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "dea834b2-0db9-49c2-a206-2251863dc891",
+ "metadata": {},
+ "source": [
+ "## TEST CASE 2: Single Dataset, Single Donor, No Batch Effects\n",
+ "- Collection: A Web Portal and Workbench for Biological Dissection of Single Cell COVID-19 Host Responses\n",
+ "- Dataset: Individual Single-Cell RNA-seq PBMC Data from Arunachalam et al.\n",
+ " - Assay: 10X\n",
+ "- Comparison:\n",
+ " - classical monocytes vs T-cells in one donor"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 32,
+ "id": "59ccd78c-0161-4f5b-a3ad-485943465cda",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "{'census_obs_query': \"is_primary_data == True and dataset_id in ['59b69042-47c2-47fd-ad03-d21beb99818f'] and disease in ['normal'] and cell_type in ['classical monocyte', 'central memory CD4-positive, alpha-beta T cell'] and donor_id == 'cov17'\",\n",
+ " 'precomputed_memento_query': \"dataset_id in ['59b69042-47c2-47fd-ad03-d21beb99818f'] and disease_ontology_term_id in ['PATO:0000461'] and cell_type_ontology_term_id in ['CL:0000860', 'CL:0000904'] and donor_id == 'cov17'\"}"
+ ]
+ },
+ "execution_count": 32,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# data\n",
+ "dataset_aruna = \"59b69042-47c2-47fd-ad03-d21beb99818f\"\n",
+ "disease_info_map = {\"normal\": \"PATO:0000461\"}\n",
+ "\n",
+ "cell_type_monocyte = \"classical monocyte\"\n",
+ "cell_type_t_cell = \"central memory CD4-positive, alpha-beta T cell\"\n",
+ "cell_type_info_map = {cell_type_monocyte: \"CL:0000860\", cell_type_t_cell: \"CL:0000904\"}\n",
+ "\n",
+ "\n",
+ "# query params\n",
+ "datasets = [dataset_aruna]\n",
+ "\n",
+ "diseases = list(disease_info_map.keys())\n",
+ "disease_ontology_ids = list(disease_info_map.values())\n",
+ "\n",
+ "donor_id = \"cov17\"\n",
+ "\n",
+ "cell_types = list(cell_type_info_map.keys())\n",
+ "cell_type_ontology_ids = list(cell_type_info_map.values())\n",
+ "\n",
+ "# a test case is encapsulated as a census query and as a query to precomputed memento cube.\n",
+ "# The precomputed memento cube is the same as the census query except that the precomputed memento cube\n",
+ "# stores ontology term IDs and such the query should be formulated with ontology term IDs\n",
+ "test_case = {\n",
+ " \"census_obs_query\": f\"is_primary_data == True and dataset_id in {datasets} and disease in {diseases} and cell_type in {cell_types} and donor_id == '{donor_id}'\",\n",
+ " \"precomputed_memento_query\": f\"dataset_id in {datasets} and disease_ontology_term_id in {disease_ontology_ids} and cell_type_ontology_term_id in {cell_type_ontology_ids} and donor_id == '{donor_id}'\",\n",
+ "}\n",
+ "\n",
+ "test_case"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3bd2c1f8-52eb-49f6-b5fd-937323618db8",
+ "metadata": {},
+ "source": [
+ "#### Get anndata object for test case"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 33,
+ "id": "3fdd1931-32e4-4874-9f6c-c496099cc38d",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "cell_type dataset_id donor_id\n",
+ "classical monocyte 59b69042-47c2-47fd-ad03-d21beb99818f cov17 2050\n",
+ "central memory CD4-positive, alpha-beta T cell 59b69042-47c2-47fd-ad03-d21beb99818f cov17 498\n",
+ "Name: count, dtype: int64"
+ ]
+ },
+ "execution_count": 33,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "adata = cellxgene_census.get_anndata(\n",
+ " census=census,\n",
+ " organism=\"homo_sapiens\",\n",
+ " obsm_layers=[\"scvi\"],\n",
+ " obs_value_filter=test_case[\"census_obs_query\"],\n",
+ ")\n",
+ "\n",
+ "adata.obs[[\"cell_type\", \"dataset_id\", \"donor_id\"]].value_counts()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7dbdfafa-2b6d-438d-baf5-f485e162a8ab",
+ "metadata": {},
+ "source": [
+ "#### Run vanilla memento on test case"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "id": "3bd6853a-1896-4b57-bd2b-f313add62be0",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/psridharan/miniconda3/envs/vanilla-memento/lib/python3.10/site-packages/memento/main.py:181: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_numeric without passing `errors` and catch exceptions explicitly instead\n",
+ " df[col] = pd.to_numeric(df[col], errors='ignore')\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Assuming sequenced to 50%, dataset specific number\n",
+ "adata.var.index = adata.var[\"feature_id\"].values # use feature_id to refer to the gene\n",
+ "adata.obs[\"q\"] = 0.15\n",
+ "\n",
+ "\n",
+ "# Classical monocyte encoded as 1\n",
+ "adata.obs[\"treatment\"] = (adata.obs[\"cell_type\"] == \"classical monocyte\").astype(int)\n",
+ "\n",
+ "# Setup memento\n",
+ "memento.setup_memento(adata, q_column=\"q\", trim_percent=0.1) # trim_percent tunes cell size calculation\n",
+ "memento.create_groups(adata, label_columns=[\"treatment\"])\n",
+ "memento.compute_1d_moments(adata, min_perc_group=0.9)\n",
+ "group_metadata = memento.get_groups(adata)\n",
+ "\n",
+ "treatment_df = group_metadata[[\"treatment\"]]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 35,
+ "id": "5f441f7d-7f38-4644-b4ad-e7e2c0cb2383",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " treatment | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | sg^1 | \n",
+ " 1 | \n",
+ "
\n",
+ " \n",
+ " | sg^0 | \n",
+ " 0 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " gene | \n",
+ " tx | \n",
+ " de_coef | \n",
+ " de_se | \n",
+ " de_pval | \n",
+ " dv_coef | \n",
+ " dv_se | \n",
+ " dv_pval | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 1122 | \n",
+ " ENSG00000000419 | \n",
+ " treatment | \n",
+ " 0.567164 | \n",
+ " 0.140124 | \n",
+ " 0.001202 | \n",
+ " -0.251849 | \n",
+ " 0.562401 | \n",
+ " 0.629074 | \n",
+ "
\n",
+ " \n",
+ " | 2045 | \n",
+ " ENSG00000001497 | \n",
+ " treatment | \n",
+ " 0.178305 | \n",
+ " 0.181198 | \n",
+ " 0.312737 | \n",
+ " -0.397483 | \n",
+ " 0.568722 | \n",
+ " 0.471106 | \n",
+ "
\n",
+ " \n",
+ " | 136 | \n",
+ " ENSG00000001629 | \n",
+ " treatment | \n",
+ " 0.573830 | \n",
+ " 0.147534 | \n",
+ " 0.001664 | \n",
+ " 0.150754 | \n",
+ " 0.466089 | \n",
+ " 0.741252 | \n",
+ "
\n",
+ " \n",
+ " | 2020 | \n",
+ " ENSG00000001631 | \n",
+ " treatment | \n",
+ " 0.356685 | \n",
+ " 0.137996 | \n",
+ " 0.010998 | \n",
+ " -0.062119 | \n",
+ " 0.595587 | \n",
+ " 0.898020 | \n",
+ "
\n",
+ " \n",
+ " | 369 | \n",
+ " ENSG00000002330 | \n",
+ " treatment | \n",
+ " 0.293128 | \n",
+ " 0.131607 | \n",
+ " 0.027994 | \n",
+ " -0.000774 | \n",
+ " 0.581164 | \n",
+ " 0.998200 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " coef | \n",
+ " z | \n",
+ " pval | \n",
+ "
\n",
+ " \n",
+ " | feature_id | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | (ENSG00000000003,) | \n",
+ " 2.141824 | \n",
+ " 4.180162 | \n",
+ " 2.913019e-05 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000419,) | \n",
+ " -0.782643 | \n",
+ " -6.067012 | \n",
+ " 1.303122e-09 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000457,) | \n",
+ " -0.087643 | \n",
+ " -0.472080 | \n",
+ " 6.368700e-01 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000460,) | \n",
+ " -0.330125 | \n",
+ " -0.906739 | \n",
+ " 3.645451e-01 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000938,) | \n",
+ " -4.277515 | \n",
+ " -15.819741 | \n",
+ " 2.274272e-56 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " coef | \n",
+ " z | \n",
+ " pval | \n",
+ "
\n",
+ " \n",
+ " | feature_id | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | (ENSG00000000419,) | \n",
+ " -0.782643 | \n",
+ " -6.067012 | \n",
+ " 1.303122e-09 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000001497,) | \n",
+ " -0.425436 | \n",
+ " -2.739698 | \n",
+ " 6.149557e-03 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000001629,) | \n",
+ " -0.833652 | \n",
+ " -6.176022 | \n",
+ " 6.573672e-10 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000001631,) | \n",
+ " -0.596159 | \n",
+ " -4.674299 | \n",
+ " 2.949596e-06 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000002330,) | \n",
+ " -0.509686 | \n",
+ " -4.300392 | \n",
+ " 1.704964e-05 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x = -1 * np.log10(vanilla_memento_result_df[\"de_pval\"])\n",
+ "y = -1 * np.log10(filtered_precomputed_memento_result_df[\"pval\"])\n",
+ "\n",
+ "plt.scatter(x, y, label=\"p values\")\n",
+ "\n",
+ "plt.plot(x, x, color=\"black\", label=\"Line with Slope 1\")\n",
+ "\n",
+ "plt.xlabel(\"vanilla\")\n",
+ "plt.ylabel(\"precomputed\")\n",
+ "\n",
+ "plt.title(\"Scatter Plot of p values\")\n",
+ "plt.legend()\n",
+ "\n",
+ "plt.grid(True)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "937e5d6c-4c88-4e47-a247-5037792ecf56",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.13"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/scientific_validation_notebooks/precomputed_memento_validation_testcase3.ipynb b/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/scientific_validation_notebooks/precomputed_memento_validation_testcase3.ipynb
new file mode 100644
index 000000000..200a0c95b
--- /dev/null
+++ b/api/python/cellxgene_census/src/cellxgene_census/experimental/diffexp/memento/scientific_validation_notebooks/precomputed_memento_validation_testcase3.ipynb
@@ -0,0 +1,718 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "6caacabf-1002-47fb-b38d-96fcb4f10133",
+ "metadata": {},
+ "source": [
+ "## Scientific validation of memento implementations\n",
+ "\n",
+ "This notebook runs through test cases for differential expression where the p-value outputs\n",
+ "of vanilla memento implementation are compared to the p-value outputs of precomputed memento\n",
+ "implementation\n",
+ "\n",
+ "**NOTE**:\n",
+ "- Run this notebook from the `cellxgene-census` repo in the `psridharan/memento-sci-val-1` branch\n",
+ "- Download estimator cube from `s3://psridharan-tmp/memento/memento-cube-census-data-2023-12-15/`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "id": "ae923ea5-7c72-4743-aea1-c854593fc1ec",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "import memento\n",
+ "import numpy as np\n",
+ "\n",
+ "import cellxgene_census\n",
+ "import cellxgene_census.experimental.diffexp.memento.diff_expr as precomputed_memento\n",
+ "\n",
+ "census = cellxgene_census.open_soma(census_version=\"2023-12-15\")\n",
+ "precomputed_memento_estimator_cube_path = (\n",
+ " \"/Users/psridharan/code/cellxgene-census/memento-cubes/memento-cube-census-data-2023-12-15\"\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "dea834b2-0db9-49c2-a206-2251863dc891",
+ "metadata": {},
+ "source": [
+ "## TEST CASE 3: Single Dataset, Single Donor, No Batch Effects\n",
+ "- Collection: A Web Portal and Workbench for Biological Dissection of Single Cell COVID-19 Host Responses\n",
+ "- Dataset: Individual Single-Cell RNA-seq PBMC Data from Arunachalam et al.\n",
+ " - Assay: 10X\n",
+ "- Comparison:\n",
+ " - classical monocytes vs natural killer cells in one donor"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "id": "59ccd78c-0161-4f5b-a3ad-485943465cda",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "{'census_obs_query': \"is_primary_data == True and dataset_id in ['59b69042-47c2-47fd-ad03-d21beb99818f'] and disease in ['normal'] and cell_type in ['classical monocyte', 'natural killer cell'] and donor_id == 'cov17'\",\n",
+ " 'precomputed_memento_query': \"dataset_id in ['59b69042-47c2-47fd-ad03-d21beb99818f'] and disease_ontology_term_id in ['PATO:0000461'] and cell_type_ontology_term_id in ['CL:0000860', 'CL:0000623'] and donor_id == 'cov17'\"}"
+ ]
+ },
+ "execution_count": 19,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# data\n",
+ "dataset_aruna = \"59b69042-47c2-47fd-ad03-d21beb99818f\"\n",
+ "disease_info_map = {\"normal\": \"PATO:0000461\"}\n",
+ "\n",
+ "cell_type_monocyte = \"classical monocyte\"\n",
+ "cell_type_nk = \"natural killer cell\"\n",
+ "cell_type_info_map = {\n",
+ " cell_type_monocyte: \"CL:0000860\",\n",
+ " cell_type_nk: \"CL:0000623\",\n",
+ "}\n",
+ "\n",
+ "\n",
+ "# query params\n",
+ "datasets = [dataset_aruna]\n",
+ "\n",
+ "diseases = list(disease_info_map.keys())\n",
+ "disease_ontology_ids = list(disease_info_map.values())\n",
+ "\n",
+ "donor_id = \"cov17\"\n",
+ "\n",
+ "cell_types = list(cell_type_info_map.keys())\n",
+ "cell_type_ontology_ids = list(cell_type_info_map.values())\n",
+ "\n",
+ "# a test case is encapsulated as a census query and as a query to precomputed memento cube.\n",
+ "# The precomputed memento cube is the same as the census query except that the precomputed memento cube\n",
+ "# stores ontology term IDs and such the query should be formulated with ontology term IDs\n",
+ "test_case = {\n",
+ " \"census_obs_query\": f\"is_primary_data == True and dataset_id in {datasets} and disease in {diseases} and cell_type in {cell_types} and donor_id == '{donor_id}'\",\n",
+ " \"precomputed_memento_query\": f\"dataset_id in {datasets} and disease_ontology_term_id in {disease_ontology_ids} and cell_type_ontology_term_id in {cell_type_ontology_ids} and donor_id == '{donor_id}'\",\n",
+ "}\n",
+ "\n",
+ "test_case"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3bd2c1f8-52eb-49f6-b5fd-937323618db8",
+ "metadata": {},
+ "source": [
+ "#### Get anndata object for test case"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "id": "3fdd1931-32e4-4874-9f6c-c496099cc38d",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "cell_type dataset_id donor_id\n",
+ "natural killer cell 59b69042-47c2-47fd-ad03-d21beb99818f cov17 2312\n",
+ "classical monocyte 59b69042-47c2-47fd-ad03-d21beb99818f cov17 2050\n",
+ "Name: count, dtype: int64"
+ ]
+ },
+ "execution_count": 20,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "adata = cellxgene_census.get_anndata(\n",
+ " census=census,\n",
+ " organism=\"homo_sapiens\",\n",
+ " obsm_layers=[\"scvi\"],\n",
+ " obs_value_filter=test_case[\"census_obs_query\"],\n",
+ ")\n",
+ "\n",
+ "adata.obs[[\"cell_type\", \"dataset_id\", \"donor_id\"]].value_counts()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7dbdfafa-2b6d-438d-baf5-f485e162a8ab",
+ "metadata": {},
+ "source": [
+ "#### Run vanilla memento on test case"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "id": "3bd6853a-1896-4b57-bd2b-f313add62be0",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/psridharan/miniconda3/envs/vanilla-memento/lib/python3.10/site-packages/memento/main.py:181: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_numeric without passing `errors` and catch exceptions explicitly instead\n",
+ " df[col] = pd.to_numeric(df[col], errors='ignore')\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Assuming sequenced to 50%, dataset specific number\n",
+ "adata.var.index = adata.var[\"feature_id\"].values # use feature_id to refer to the gene\n",
+ "adata.obs[\"q\"] = 0.15\n",
+ "\n",
+ "\n",
+ "# Classical monocyte encoded as 1\n",
+ "adata.obs[\"treatment\"] = (adata.obs[\"cell_type\"] == \"classical monocyte\").astype(int)\n",
+ "\n",
+ "# Setup memento\n",
+ "memento.setup_memento(adata, q_column=\"q\", trim_percent=0.1) # trim_percent tunes cell size calculation\n",
+ "memento.create_groups(adata, label_columns=[\"treatment\"])\n",
+ "memento.compute_1d_moments(adata, min_perc_group=0.9)\n",
+ "group_metadata = memento.get_groups(adata)\n",
+ "\n",
+ "treatment_df = group_metadata[[\"treatment\"]]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "id": "5f441f7d-7f38-4644-b4ad-e7e2c0cb2383",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " treatment | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | sg^1 | \n",
+ " 1 | \n",
+ "
\n",
+ " \n",
+ " | sg^0 | \n",
+ " 0 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " gene | \n",
+ " tx | \n",
+ " de_coef | \n",
+ " de_se | \n",
+ " de_pval | \n",
+ " dv_coef | \n",
+ " dv_se | \n",
+ " dv_pval | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 1121 | \n",
+ " ENSG00000000419 | \n",
+ " treatment | \n",
+ " 0.338290 | \n",
+ " 0.076319 | \n",
+ " 0.000742 | \n",
+ " -0.098810 | \n",
+ " 0.275596 | \n",
+ " 0.720456 | \n",
+ "
\n",
+ " \n",
+ " | 1429 | \n",
+ " ENSG00000000938 | \n",
+ " treatment | \n",
+ " 0.729090 | \n",
+ " 0.034047 | \n",
+ " 0.000004 | \n",
+ " 0.025821 | \n",
+ " 0.227044 | \n",
+ " 0.913817 | \n",
+ "
\n",
+ " \n",
+ " | 1526 | \n",
+ " ENSG00000001084 | \n",
+ " treatment | \n",
+ " -0.014172 | \n",
+ " 0.095271 | \n",
+ " 0.883423 | \n",
+ " 0.113890 | \n",
+ " 0.282276 | \n",
+ " 0.674865 | \n",
+ "
\n",
+ " \n",
+ " | 2039 | \n",
+ " ENSG00000001497 | \n",
+ " treatment | \n",
+ " -0.116158 | \n",
+ " 0.097860 | \n",
+ " 0.236753 | \n",
+ " -0.202672 | \n",
+ " 0.321127 | \n",
+ " 0.525895 | \n",
+ "
\n",
+ " \n",
+ " | 132 | \n",
+ " ENSG00000001629 | \n",
+ " treatment | \n",
+ " 0.087122 | \n",
+ " 0.075774 | \n",
+ " 0.245351 | \n",
+ " -0.279639 | \n",
+ " 0.217676 | \n",
+ " 0.193561 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " coef | \n",
+ " z | \n",
+ " pval | \n",
+ "
\n",
+ " \n",
+ " | feature_id | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | (ENSG00000000419,) | \n",
+ " 0.865093 | \n",
+ " 20.125887 | \n",
+ " 4.378676e-90 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000457,) | \n",
+ " 0.391230 | \n",
+ " 4.385940 | \n",
+ " 1.154858e-05 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000460,) | \n",
+ " 0.343338 | \n",
+ " 1.998222 | \n",
+ " 4.569254e-02 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000938,) | \n",
+ " 1.246206 | \n",
+ " 67.162439 | \n",
+ " 0.000000e+00 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000971,) | \n",
+ " -1.589614 | \n",
+ " -15.843008 | \n",
+ " 1.571209e-56 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " coef | \n",
+ " z | \n",
+ " pval | \n",
+ "
\n",
+ " \n",
+ " | feature_id | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | (ENSG00000000419,) | \n",
+ " 0.865093 | \n",
+ " 20.125887 | \n",
+ " 4.378676e-90 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000000938,) | \n",
+ " 1.246206 | \n",
+ " 67.162439 | \n",
+ " 0.000000e+00 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000001084,) | \n",
+ " 0.531871 | \n",
+ " 8.897293 | \n",
+ " 5.722510e-19 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000001497,) | \n",
+ " 0.392839 | \n",
+ " 6.355219 | \n",
+ " 2.081294e-10 | \n",
+ "
\n",
+ " \n",
+ " | (ENSG00000001629,) | \n",
+ " 0.602583 | \n",
+ " 12.930142 | \n",
+ " 3.042548e-38 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "