tutorial/scc_analysis_colocalization_conditions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cell state dynamics of skin squamous cell carcinoma dependent on colocalized cell population\n",
    "Exdyn can estimate cell state dynamics dependent on multivariate extrinsic factors such as colocalization with neighboring cells and chromatin accessibility. Exdyn provide several downstream analysis for exploring the contributing fetures of the extrinsic factors. Here, we estimate cell state dynamics of human skin squamous cell carcinoma dependen on the colocalization scores with cells other than cancer cells, which were estimated by DeepCOLOR (). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model setup and optimization\n",
    "Exdyn accepts AnnData object for dynamics estimation. The model parameters such as the number of neural network layers can be specified through a dictionary of the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! pip install exdyn torchtext==0.6.0\n",
    "! wget -O tumor_adata.h5ad https://www.dropbox.com/scl/fi/nsxr4y5t5seb5i8zs5dfh/tumor_adata.h5ad?rlkey=zizkuf2kbdmepy0fwhptfg883&dl=0\n",
    "! wget -O coloc_adata.h5ad https://www.dropbox.com/scl/fi/ymkgwx4fc14cqgkr8ctyk/coloc_adata.h5ad?rlkey=c6k06dzrz51ym9dremx8hcisb&dl=0\n",
    "! wget -O lt_df.csv https://www.dropbox.com/scl/fi/863c4mgfkgxe4qk2ndh6j/lt_df.csv?rlkey=5f8lw3rbhrb7bnnzt1k7a7ox0&dl=0\n",
    "import scanpy as sc\n",
    "from exdyn import workflow, modules, visualization\n",
    "import numpy as np\n",
    "import torch\n",
    "import random\n",
    "seed = 42\n",
    "\n",
    "random.seed(seed)\n",
    "np.random.seed(seed)\n",
    "torch.manual_seed(seed)\n",
    "torch.backends.cudnn.benchmark = False\n",
    "torch.backends.cudnn.deterministic = True\n",
    "\n",
    "\n",
    "adata_file_path = 'tumor_adata.h5ad'\n",
    "coloc_adata_file_path = 'coloc_adata.h5ad'\n",
    "lt_df_file_path = 'lt_df.csv'\n",
    "adata = sc.read_h5ad(adata_file_path)\n",
    "top_corr_genes, gene_names = workflow.extract_top_corr_genes(adata)\n",
    "adata.var['dynamics_genes'] = False\n",
    "adata.var['dynamics_genes'][top_corr_genes] = True\n",
    "\n",
    "model_params = {\n",
    "            'x_dim': 100,\n",
    "            'z_dim': 10,\n",
    "            'enc_z_h_dim': 128, 'enc_d_h_dim': 128, 'dec_z_h_dim': 128,\n",
    "            'num_enc_z_layers': 2, 'num_enc_d_layers': 2,\n",
    "            'num_dec_z_layers': 2, 'use_vamp': False, 'loss_mode': 'nb'\n",
    "}\n",
    "checkpoint_dirname = 'checkpoint'\n",
    "adata.obs['sample'] = 'Sample'\n",
    "est_adata, model = workflow.conduct_cvicdyf_inference(adata, model_params, checkpoint_dirname, batch_size=128, two_step=False, dyn_mode=False, epoch=1000, patience=50, module=modules.Cvicdyf, use_highly_variable=False, batch_key='sample', condition_key='coloc', cond_in_obsm=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize purturbation of cell state dynamics by extrinsic factors\n",
    "Exdyn can estimate the purturbation effect on cell state dynamics by extrinsic factors through singular value decompotion of Jacobian matrix.\n",
    "Firstly, you can visualize the cell state dynamics given original extrinsic factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est_adata.obs['Annotation'] = est_adata.obs['level3_celltype']\n",
    "visualization.plot_mean_flow(est_adata, cluster_key='Annotation', legend_loc='on data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, you can conduct SVD of Jacobian matrix and identify the pricipal purtubations and the effects on cell state dynamics. Here, you can visualize the purtubed dynamics for top 3 modes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_modes = 4\n",
    "est_adata = workflow.conduct_jacobian_analysis(est_adata, model, top_modes=top_modes)\n",
    "from matplotlib import pyplot as plt\n",
    "fig, axes = plt.subplots(1, top_modes, figsize=(6 * top_modes, 5 * 1))\n",
    "for mode in range(top_modes):\n",
    "   ax = axes[mode] \n",
    "   visualization.diff_visualize(est_adata.obsm[f'dd_dc_embed{mode}'], est_adata.obsm['X_umap'], ax)\n",
    "   ax.set_title(f\"Dynamics difference of mode {mode + 1}\")\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note that the contribution of each feature of extrinsic factors to each mode of principal purturbation can be found in `adata.uns['dd_dc_v']`. This can be useful for identify the fetures which can cause specific cell state transitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est_adata.uns['dd_dc_v']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also quantify the induction of each cell cluster by each purturbation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mode = 3\n",
    "sign = 1\n",
    "cluster = 'TSK'\n",
    "est_adata = workflow.analyze_flux_into_clusters(est_adata, cluster, mode, sign, 'Annotation')\n",
    "sc.pl.umap(est_adata, color=['Annotation', f'mode{mode * sign}_flux_into_{cluster}'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyze cell cell interacraction\n",
    "In the case where the extrinsic factors are colozalizaiton scores with neighboring cells, you can analyze cell cell interaction which mediate specific mode of principal purturbation patterns of cell state dynamics. Firstly, you can annotate and visualize the contribution of neighboring cells to each mode of the purturbations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coloc_adata = sc.read_h5ad(coloc_adata_file_path)\n",
    "cov_labels = []\n",
    "for i in range(top_modes):\n",
    "    cov_label = f'coloc_cov_mode{i}'\n",
    "    coloc_adata.obs[cov_label] = est_adata.uns['dd_dc_v'][:, i]\n",
    "    cov_labels.append(cov_label)\n",
    "\n",
    "sc.pl.umap(coloc_adata, color=cov_labels + ['level2_celltype'], legend_loc='on data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we extract cancer cells directed TSK by mode 1 perturbation and quantify the activity of each ligand for the perturbation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "top_tsk_adata = est_adata[est_adata.obs[f'mode{mode * sign}_flux_into_{cluster}'] > np.quantile(est_adata.obs[f'mode{mode * sign}_flux_into_{cluster}'], 0.95)]\n",
    "top_tsk_adata = workflow.calc_covdiff(top_tsk_adata, model, mode, sign, 'coloc', 'sample', sep_num=5)\n",
    "condiff_vec = pd.Series(top_tsk_adata.layers[f'cov_diff_{sign * mode}'].mean(axis=0), top_tsk_adata.var_names)\n",
    "lt_df = pd.read_csv(lt_df_file_path, index_col=0)\n",
    "ligacs = workflow.analyze_ligand_activity(condiff_vec, lt_df)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to derive plausible candidates for cell cell interaction, you can intersects genes with high ligand activity with positive DEGs for contributing neighboring cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coloc_adata.obs['target_coloc'] = (coloc_adata.obs['coloc_cov_mode1'] > np.quantile(coloc_adata.obs['coloc_cov_mode1'], 0.9)).astype(str)\n",
    "sc.tl.rank_genes_groups(coloc_adata, 'target_coloc', groups=['True'], reference='False', method='wilcoxon')\n",
    "deg_df = sc.get.rank_genes_groups_df(coloc_adata, 'True')\n",
    "pos_degs = deg_df.query('pvals_adj < 0.05 and logfoldchanges > 1').names\n",
    "\n",
    "ligacs[np.intersect1d(pos_degs, ligacs.index)].sort_values(ascending=False)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "functorch",
   "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.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}