scvi-tools.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering and differential expression analysis using scVI-tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import scanpy as sc\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import json\n",
    "import anndata\n",
    "import scvi\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_name = 'PBMC1' # modify this\n",
    "data_in_path = './data/{}/filtered/10X/'.format(dataset_name)\n",
    "data_out_path = './data/{}/scvi-tools/'.format(dataset_name)\n",
    "default_path = '{}default/'.format(data_out_path)\n",
    "celltypist_path = '{}celltypist/'.format(data_out_path)\n",
    "antibody_path = '{}antibody/'.format(data_out_path)\n",
    "nclusters_celltypist_path = './data/{}/celltypist/nclusters.json'.format(dataset_name)\n",
    "nclusters_antibody_path = './data/{}/antibody_annotation/nclusters_postproc.json'.format(dataset_name)\n",
    "\n",
    "top_number_of_markers = 500\n",
    "min_cluster_size = 0\n",
    "\n",
    "with open(nclusters_celltypist_path) as f:\n",
    "    nclusters_celltypist = json.load(f)['nclusters']\n",
    "with open(nclusters_antibody_path) as f:\n",
    "    nclusters_antibody = json.load(f)['nclusters']\n",
    "\n",
    "nclusters_threshold = 1/10\n",
    "min_ncluster_celltypist = round(nclusters_celltypist - nclusters_threshold*nclusters_celltypist)\n",
    "max_ncluster_celltypist = round(nclusters_celltypist + nclusters_threshold*nclusters_celltypist)\n",
    "min_ncluster_antibody = round(nclusters_antibody - nclusters_threshold*nclusters_antibody)\n",
    "max_ncluster_antibody = round(nclusters_antibody + nclusters_threshold*nclusters_antibody)\n",
    "\n",
    "if not os.path.exists(default_path):\n",
    "    os.makedirs(default_path)\n",
    "\n",
    "if not os.path.exists(celltypist_path):\n",
    "    os.makedirs(celltypist_path)\n",
    "\n",
    "if not os.path.exists(antibody_path):\n",
    "    os.makedirs(antibody_path)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data loading and preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dataset loading "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = sc.read_10x_mtx(\n",
    "    data_in_path,\n",
    "    var_names='gene_symbols',\n",
    "    cache=False\n",
    ")\n",
    "adata.var_names_make_unique()\n",
    "adata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Studying feature variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix = adata.X\n",
    "matrix = matrix.todense()\n",
    "neg_variances = np.sort(-np.var(matrix, axis=0))\n",
    "sorted_log_variances = [np.log(-i) for i in neg_variances.T][:2000]\n",
    "plt.scatter([i for i in range(len(sorted_log_variances))], sorted_log_variances)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data normalization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.layers[\"counts\"] = adata.X.copy()\n",
    "sc.pp.normalize_total(adata, target_sum=1e4)\n",
    "sc.pp.log1p(adata)\n",
    "adata.raw = adata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Feature selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "j8gBNU_7txZU"
   },
   "outputs": [],
   "source": [
    "sc.pp.highly_variable_genes(\n",
    "    adata,\n",
    "    n_top_genes=500, # chosen based on elbow plot above\n",
    "    subset=True,\n",
    "    layer=\"counts\",\n",
    "    flavor=\"seurat_v3\",\n",
    ")\n",
    "sc.pl.highly_variable_genes(adata)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "RAWmYvjkZ0bW"
   },
   "source": [
    "## Model training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "xPnlTctkZ0bX"
   },
   "outputs": [],
   "source": [
    "scvi.model.SCVI.setup_anndata(\n",
    "    adata,\n",
    "    layer=\"counts\",\n",
    ")\n",
    "\n",
    "model = scvi.model.SCVI(adata)\n",
    "model.train(use_gpu=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "iVy8Eo1-Z0be"
   },
   "source": [
    "## Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save latent representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "R4FwiBpUtQJa"
   },
   "outputs": [],
   "source": [
    "latent = model.get_latent_representation()\n",
    "adata.obsm[\"X_scVI\"] = latent\n",
    "adata.layers[\"scvi_normalized\"] = model.get_normalized_expression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot pca explained variance ratio to choose number of components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "S_IJVCZknAfC"
   },
   "outputs": [],
   "source": [
    "sc.tl.pca(adata, svd_solver=\"arpack\", n_comps=20, use_highly_variable=True)\n",
    "sc.pl.pca_variance_ratio(adata, log=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clustering with default parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform the clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ME5yLdHtxHsb"
   },
   "outputs": [],
   "source": [
    "sc.pp.neighbors(\n",
    "    adata,\n",
    "    n_pcs=10, # chosen based on elbow plot above\n",
    "    use_rep=\"X_scVI\"\n",
    ")\n",
    "sc.tl.leiden(adata, key_added=\"leiden_scVI\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the clustering in the PCA space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.pca(\n",
    "    adata,\n",
    "    color=[\"leiden_scVI\"],\n",
    "    size=20,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the clustering results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "default_clusters_df = pd.DataFrame(adata.obs['leiden_scVI'])\n",
    "default_clusters_df = default_clusters_df.rename(columns={'leiden_scVI': 'cluster'})\n",
    "default_clusters_df.index.names = ['cell']\n",
    "default_clusters_df['cluster'] = default_clusters_df['cluster'].astype(int) + 1\n",
    "default_clusters_df.to_csv(default_path+'clustering_labels.csv')\n",
    "default_num_clusters = len(default_clusters_df['cluster'].unique())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differential expression on default clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform the differential expression analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_cells = np.array([])\n",
    "adata_train = adata.copy()\n",
    "for cluster in range(default_num_clusters):\n",
    "    cluster_cells = default_clusters_df[default_clusters_df['cluster'] == cluster]\n",
    "    cluster_cells = cluster_cells.sample(frac=0.2)\n",
    "    test_cells = np.concatenate([test_cells, cluster_cells.index.values])\n",
    "    adata_train = adata_train[~adata_train.obs.index.isin(cluster_cells.index)].copy()\n",
    "pd.DataFrame(test_cells, columns=['cell']).to_csv(f'{default_path}/test_cells.csv', index=False)\n",
    "model.adata = adata_train\n",
    "default_de_df = model.differential_expression(\n",
    "    groupby=\"leiden_scVI\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "default_markers = {}\n",
    "default_cats = adata_train.obs.leiden_scVI.cat.categories\n",
    "for i, c in enumerate(default_cats):\n",
    "    cell_type_df = default_de_df.loc[default_de_df.group1 == c]\n",
    "    # if uncommented there are not enough genes\n",
    "    #cell_type_df = cell_type_df[cell_type_df[\"bayes_factor\"] > 3]\n",
    "    #cell_type_df = cell_type_df[cell_type_df[\"non_zeros_proportion1\"] > 0.1]\n",
    "    default_markers[c] = cell_type_df.sort_values('lfc_mean', ascending=False).index.tolist()[:top_number_of_markers]\n",
    "\n",
    "default_markers_df = pd.DataFrame(columns=['gene', 'cluster', 'rank'])\n",
    "for i in range(default_num_clusters):\n",
    "    tmp_list = [ [default_markers[str(i)][j], i, j] for j in range(top_number_of_markers)]\n",
    "    tmp_df = pd.DataFrame(tmp_list, columns=['gene', 'cluster', 'rank'])\n",
    "    default_markers_df = pd.concat([default_markers_df, tmp_df])\n",
    "default_markers_df['cluster'] += 1\n",
    "default_markers_df['rank'] += 1\n",
    "default_markers_df.to_csv(default_path+'markers.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clustering tuning resolution according to celltypist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Nummber of clusters to find: {}\".format(nclusters_celltypist))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ids of clusters bigger than 40 cells\n",
    "mapping = pd.read_csv('./data/{}/celltypist/celltypist_mapping.csv'.format(dataset_name).format(dataset_name))\n",
    "counts = pd.read_csv('./data/{}/celltypist/celltypist_annotation_counts.csv'.format(dataset_name))\n",
    "mapping_counts = mapping.merge(counts, left_on='go', right_on='cluster.ids')\n",
    "mapping_counts = mapping_counts[mapping_counts['count'] > min_cluster_size]\n",
    "clusters_ids_to_keep = mapping_counts['id']\n",
    "\n",
    "# get barcodes of cells in clusters bigger than 40 cells\n",
    "celltypist_labels_df = pd.read_csv('./data/{}/celltypist/celltypist_labels.csv'.format(dataset_name))\n",
    "celltypist_labels_df = celltypist_labels_df[celltypist_labels_df['cluster.ids'].isin(clusters_ids_to_keep)]\n",
    "barcodes_to_keep = celltypist_labels_df['cell']\n",
    "barcodes_to_keep = [barcode[:-2] for barcode in barcodes_to_keep]\n",
    "subset_cells = adata.obs_names.isin(barcodes_to_keep)\n",
    "adata_celltypist = adata[subset_cells, :].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scvi.model.SCVI.setup_anndata(\n",
    "    adata_celltypist,\n",
    "    layer=\"counts\",\n",
    ")\n",
    "\n",
    "model = scvi.model.SCVI(adata_celltypist)\n",
    "model.train(use_gpu=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save latent representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latent = model.get_latent_representation()\n",
    "adata_celltypist.obsm[\"X_scVI\"] = latent\n",
    "adata_celltypist.layers[\"scvi_normalized\"] = model.get_normalized_expression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot pca explained variance ratio to choose number of components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.pca(adata_celltypist, svd_solver=\"arpack\", n_comps=20, use_highly_variable=True)\n",
    "sc.pl.pca_variance_ratio(adata_celltypist, log=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform the clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.neighbors(\n",
    "    adata_celltypist,\n",
    "    n_pcs=10, # chosen based on elbow plot above\n",
    "    use_rep=\"X_scVI\"\n",
    ")\n",
    "sc.tl.leiden(adata_celltypist, key_added=\"leiden_scVI\")\n",
    "\n",
    "max_resolution = 3\n",
    "min_resolution = 0\n",
    "\n",
    "while True:\n",
    "    resolution = (max_resolution + min_resolution)/2\n",
    "    print(\"Trying resolution: {}\".format(resolution))\n",
    "    sc.tl.leiden(adata_celltypist, key_added=\"leiden_scVI\", resolution=resolution)\n",
    "    num_clusters = adata_celltypist.obs.leiden_scVI.values.categories.nunique()\n",
    "    print(\"Number of clusters found: {}\".format(num_clusters))\n",
    "\n",
    "    if num_clusters >= min_ncluster_celltypist and num_clusters <= max_ncluster_celltypist:\n",
    "        break\n",
    "    elif num_clusters <  min_ncluster_celltypist:\n",
    "        min_resolution = resolution\n",
    "    else:\n",
    "        max_resolution = resolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the clustering in the PCA space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.pca(\n",
    "    adata_celltypist,\n",
    "    color=[\"leiden_scVI\"],\n",
    "    size=20,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the clustering results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "celltypist_clusters_df = pd.DataFrame(adata_celltypist.obs['leiden_scVI'])\n",
    "celltypist_clusters_df = celltypist_clusters_df.rename(columns={'leiden_scVI': 'cluster'})\n",
    "celltypist_clusters_df.index.names = ['cell']\n",
    "celltypist_clusters_df['cluster'] = celltypist_clusters_df['cluster'].astype(int) + 1\n",
    "celltypist_clusters_df.to_csv(celltypist_path+'clustering_labels.csv')\n",
    "celltypist_num_clusters = len(celltypist_clusters_df['cluster'].unique())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differential expression on clusters tuned according to celltypist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform the differential expression analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "celltypist_de_df = model.differential_expression(\n",
    "    groupby=\"leiden_scVI\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "celltypist_markers = {}\n",
    "cats = adata_celltypist.obs.leiden_scVI.cat.categories\n",
    "for i, c in enumerate(cats):\n",
    "    cell_type_df = celltypist_de_df.loc[celltypist_de_df.group1 == c]\n",
    "    # if uncommented there are not enough genes\n",
    "    #cell_type_df = cell_type_df[cell_type_df[\"bayes_factor\"] > 3]\n",
    "    #cell_type_df = cell_type_df[cell_type_df[\"non_zeros_proportion1\"] > 0.1]\n",
    "    celltypist_markers[c] = cell_type_df.sort_values('lfc_mean', ascending=False).index.tolist()[:top_number_of_markers]\n",
    "\n",
    "celltypist_markers_df = pd.DataFrame(columns=['gene', 'cluster', 'rank'])\n",
    "for i in range(celltypist_num_clusters):\n",
    "    tmp_list = [ [celltypist_markers[str(i)][j], i, j] for j in range(top_number_of_markers)]\n",
    "    tmp_df = pd.DataFrame(tmp_list, columns=['gene', 'cluster', 'rank'])\n",
    "    celltypist_markers_df = pd.concat([celltypist_markers_df, tmp_df],  ignore_index=True)\n",
    "\n",
    "celltypist_markers_df['cluster'] += 1\n",
    "celltypist_markers_df['rank'] += 1\n",
    "celltypist_markers_df.to_csv(celltypist_path+'markers.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clustering tuning resolution according to protein surface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Nummber of clusters to find: {}\".format(nclusters_antibody))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get barcodes of cells labelled using protein surface\n",
    "antibody_labels_df = pd.read_csv('./data/{}/antibody_annotation/antibody_labels_postproc.csv'.format(dataset_name))\n",
    "barcodes_to_keep = antibody_labels_df['cell']\n",
    "\n",
    "subset_cells = adata.obs_names.isin(barcodes_to_keep)\n",
    "adata_antibody = adata[subset_cells, :].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scvi.model.SCVI.setup_anndata(\n",
    "    adata_antibody,\n",
    "    layer=\"counts\",\n",
    ")\n",
    "\n",
    "model = scvi.model.SCVI(adata_antibody)\n",
    "model.train(use_gpu=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save latent representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latent = model.get_latent_representation()\n",
    "adata_antibody.obsm[\"X_scVI\"] = latent\n",
    "adata_antibody.layers[\"scvi_normalized\"] = model.get_normalized_expression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot pca explained variance ratio to choose number of components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.pca(adata_antibody, svd_solver=\"arpack\", n_comps=20, use_highly_variable=True)\n",
    "sc.pl.pca_variance_ratio(adata_antibody, log=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform the clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.neighbors(\n",
    "    adata_antibody,\n",
    "    n_pcs=10, # chosen based on elbow plot above\n",
    "    use_rep=\"X_scVI\"\n",
    ")\n",
    "sc.tl.leiden(adata_antibody, key_added=\"leiden_scVI\")\n",
    "\n",
    "max_resolution = 3\n",
    "min_resolution = 0\n",
    "\n",
    "while True:\n",
    "    resolution = (max_resolution + min_resolution)/2\n",
    "    print(\"Trying resolution: {}\".format(resolution))\n",
    "    sc.tl.leiden(adata_antibody, key_added=\"leiden_scVI\", resolution=resolution)\n",
    "    num_clusters = adata_antibody.obs.leiden_scVI.values.categories.nunique()\n",
    "    print(\"Number of clusters found: {}\".format(num_clusters))\n",
    "\n",
    "    if num_clusters >= min_ncluster_antibody and num_clusters <= max_ncluster_antibody:\n",
    "        break\n",
    "    elif num_clusters < min_ncluster_antibody:\n",
    "        min_resolution = resolution\n",
    "    else:\n",
    "        max_resolution = resolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the clustering in the PCA space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.pca(\n",
    "    adata_antibody,\n",
    "    color=[\"leiden_scVI\"],\n",
    "    size=20,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the clustering results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "antibody_clusters_df = pd.DataFrame(adata_antibody.obs['leiden_scVI'])\n",
    "antibody_clusters_df = antibody_clusters_df.rename(columns={'leiden_scVI': 'cluster'})\n",
    "antibody_clusters_df.index.names = ['cell']\n",
    "antibody_clusters_df['cluster'] = antibody_clusters_df['cluster'].astype(int) + 1\n",
    "antibody_clusters_df.to_csv(antibody_path+'clustering_labels.csv')\n",
    "antibody_num_clusters = len(antibody_clusters_df['cluster'].unique())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differential expression on clusters tuned according to protein surface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform the differential expression analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "antibody_de_df = model.differential_expression(\n",
    "    groupby=\"leiden_scVI\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "antibody_markers = {}\n",
    "cats = adata_antibody.obs.leiden_scVI.cat.categories\n",
    "for i, c in enumerate(cats):\n",
    "    cell_type_df = antibody_de_df.loc[antibody_de_df.group1 == c]\n",
    "    # if uncommented there are not enough genes\n",
    "    #cell_type_df = cell_type_df[cell_type_df[\"bayes_factor\"] > 3]\n",
    "    #cell_type_df = cell_type_df[cell_type_df[\"non_zeros_proportion1\"] > 0.1]\n",
    "    antibody_markers[c] = cell_type_df.sort_values('lfc_mean', ascending=False).index.tolist()[:top_number_of_markers]\n",
    "\n",
    "antibody_markers_df = pd.DataFrame(columns=['gene', 'cluster', 'rank'])\n",
    "for i in range(antibody_num_clusters):\n",
    "    tmp_list = [ [antibody_markers[str(i)][j], i, j] for j in range(top_number_of_markers)]\n",
    "    tmp_df = pd.DataFrame(tmp_list, columns=['gene', 'cluster', 'rank'])\n",
    "    antibody_markers_df = pd.concat([antibody_markers_df, tmp_df],  ignore_index=True)\n",
    "\n",
    "antibody_markers_df['cluster'] += 1\n",
    "antibody_markers_df['rank'] += 1\n",
    "antibody_markers_df.to_csv(antibody_path+'markers.csv', index=False)"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "name": "api_overview.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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.11.5"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}