{
  "cells": [
    {
      "cell_type": "raw",
      "metadata": {},
      "source": [
        "---\n",
        "description: Assignment of cell identities based on gene expression\n",
        "  patterns using reference data.\n",
        "subtitle:  Scanpy Toolkit\n",
        "title:  Celltype prediction\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div>\n",
        "\n",
        "> **Note**\n",
        ">\n",
        "> Code chunks run Python commands unless it starts with `%%bash`, in\n",
        "> which case, those chunks run shell commands.\n",
        "\n",
        "</div>\n",
        "\n",
        "Celltype prediction can either be performed on indiviudal cells where\n",
        "each cell gets a predicted celltype label, or on the level of clusters.\n",
        "All methods are based on similarity to other datasets, single cell or\n",
        "sorted bulk RNAseq, or uses known marker genes for each cell type.\\\n",
        "Ideally celltype predictions should be run on each sample separately and\n",
        "not using the integrated data. In this case we will select one sample\n",
        "from the Covid data, `ctrl_13` and predict celltype by cell on that\n",
        "sample.\\\n",
        "Some methods will predict a celltype to each cell based on what it is\n",
        "most similar to, even if that celltype is not included in the reference.\n",
        "Other methods include an uncertainty so that cells with low similarity\n",
        "scores will be unclassified.\\\n",
        "There are multiple different methods to predict celltypes, here we will\n",
        "just cover a few of those.\n",
        "\n",
        "Here we will use a reference PBMC dataset that we get from scanpy\n",
        "datasets and classify celltypes based on two methods:\n",
        "\n",
        "-   Using scanorama for integration just as in the integration lab, and\n",
        "    then do label transfer based on closest neighbors.\n",
        "-   Using ingest to project the data onto the reference data and\n",
        "    transfer labels.\n",
        "-   Using Celltypist to predicted with a pretrained pbmc model or with\n",
        "    an own model based on the same reference data as the other methods.\n",
        "\n",
        "First, lets load required libraries"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: libraries\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import scanpy as sc\n",
        "import matplotlib.pyplot as plt\n",
        "import warnings\n",
        "import os\n",
        "import subprocess\n",
        "\n",
        "warnings.simplefilter(action=\"ignore\", category=Warning)\n",
        "\n",
        "# verbosity: errors (0), warnings (1), info (2), hints (3)\n",
        "sc.settings.verbosity = 2\n",
        "sc.settings.set_figure_params(dpi=80)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's read in the saved Covid-19 data object from the clustering step."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: fetch-data\n",
        "# download pre-computed data if missing or long compute\n",
        "fetch_data = True\n",
        "\n",
        "# url for source and intermediate data\n",
        "path_data = \"https://nextcloud.dc.scilifelab.se/public.php/webdav\"\n",
        "curl_upass = \"zbC5fr2LbEZ9rSE:scRNAseq2025\"\n",
        "\n",
        "path_results = \"data/covid/results\"\n",
        "if not os.path.exists(path_results):\n",
        "    os.makedirs(path_results, exist_ok=True)\n",
        "\n",
        "path_file = \"data/covid/results/scanpy_covid_qc_dr_int_cl.h5ad\"\n",
        "if fetch_data and not os.path.exists(path_file):\n",
        "    file_url = os.path.join(path_data, \"covid/results_scanpy/scanpy_covid_qc_dr_int_cl.h5ad\")\n",
        "    subprocess.call([\"curl\", \"-u\", curl_upass, \"-o\", path_file, file_url ])    \n",
        "\n",
        "adata = sc.read_h5ad(path_file)\n",
        "adata"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: check-data\n",
        "adata.uns['log1p']['base']=None\n",
        "print(adata.shape)\n",
        "# have only variable genes in X, use raw instead.\n",
        "adata = adata.raw.to_adata()\n",
        "print(adata.shape)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Subset one patient."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: subset\n",
        "adata = adata[adata.obs[\"sample\"] == \"ctrl_13\",:]\n",
        "print(adata.shape)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: count-cells\n",
        "adata.obs[\"leiden_0.6\"].value_counts()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Some clusters have very few cells from this individual, so any cluster\n",
        "comparisons may be biased by this."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: plot-umap\n",
        "sc.pl.umap(\n",
        "    adata, color=[\"leiden_0.6\"], palette=sc.pl.palettes.default_20\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Reference data\n",
        "\n",
        "Load the reference data from `scanpy.datasets`. It is the annotated and\n",
        "processed pbmc3k dataset from 10x."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: fetch-ref\n",
        "adata_ref = sc.datasets.pbmc3k_processed() \n",
        "\n",
        "adata_ref.obs['sample']='pbmc3k'\n",
        "\n",
        "print(adata_ref.shape)\n",
        "adata_ref.obs"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As you can see, the celltype annotation is in the metadata column\n",
        "`louvain`, so that is the column we will have to use for classification."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: plot-ref\n",
        "sc.pl.umap(adata_ref, color='louvain')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make sure we have the same genes in both datset by taking the\n",
        "intersection"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: intersect-genes\n",
        "# before filtering genes, store the full matrix in raw.\n",
        "adata.raw = adata\n",
        "# also store the umap in a new slot as it will get overwritten\n",
        "adata.obsm[\"X_umap_uncorr\"] = adata.obsm[\"X_umap\"]\n",
        "\n",
        "print(adata_ref.shape[1])\n",
        "print(adata.shape[1])\n",
        "var_names = adata_ref.var_names.intersection(adata.var_names)\n",
        "print(len(var_names))\n",
        "\n",
        "adata_ref = adata_ref[:, var_names]\n",
        "adata = adata[:, var_names]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First we need to rerun pca and umap with the same gene set for both\n",
        "datasets."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: process-ref\n",
        "sc.pp.pca(adata_ref)\n",
        "sc.pp.neighbors(adata_ref)\n",
        "sc.tl.umap(adata_ref)\n",
        "sc.pl.umap(adata_ref, color='louvain')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: process-data\n",
        "sc.pp.pca(adata)\n",
        "sc.pp.neighbors(adata)\n",
        "sc.tl.umap(adata)\n",
        "sc.pl.umap(adata, color='leiden_0.6')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Integrate with scanorama"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: integrate\n",
        "import scanorama\n",
        "\n",
        "#subset the individual dataset to the same variable genes as in MNN-correct.\n",
        "alldata = dict()\n",
        "alldata['ctrl']=adata\n",
        "alldata['ref']=adata_ref\n",
        "\n",
        "#convert to list of AnnData objects\n",
        "adatas = list(alldata.values())\n",
        "\n",
        "# run scanorama.integrate\n",
        "scanorama.integrate_scanpy(adatas, dimred = 50)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: scanorama-results\n",
        "# add in sample info\n",
        "adata_ref.obs['sample']='pbmc3k'\n",
        "\n",
        "# create a merged scanpy object and add in the scanorama \n",
        "adata_merged = alldata['ctrl'].concatenate(alldata['ref'], batch_key='sample', batch_categories=['ctrl','pbmc3k'])\n",
        "\n",
        "embedding = np.concatenate([ad.obsm['X_scanorama'] for ad in adatas], axis=0)\n",
        "adata_merged.obsm['Scanorama'] = embedding"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: process-scanorama\n",
        "#run  umap.\n",
        "sc.pp.neighbors(adata_merged, n_pcs =50, use_rep = \"Scanorama\")\n",
        "sc.tl.umap(adata_merged)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: plot-scanorama\n",
        "sc.pl.umap(adata_merged, color=[\"sample\",\"louvain\"])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Label transfer\n",
        "\n",
        "Using the functions from the [Spatial\n",
        "tutorial](https://scanpy.readthedocs.io/en/stable/tutorials/spatial/integration-scanorama.html)\n",
        "from Scanpy we will calculate normalized cosine distances between the\n",
        "two datasets and tranfer labels to the celltype with the highest scores."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: transfer\n",
        "from sklearn.metrics.pairwise import cosine_distances\n",
        "\n",
        "distances = 1 - cosine_distances(\n",
        "    adata_merged[adata_merged.obs['sample'] == \"pbmc3k\"].obsm[\"Scanorama\"],\n",
        "    adata_merged[adata_merged.obs['sample'] == \"ctrl\"].obsm[\"Scanorama\"],\n",
        ")\n",
        "\n",
        "def label_transfer(dist, labels, index):\n",
        "    lab = pd.get_dummies(labels)\n",
        "    class_prob = lab.to_numpy().T @ dist\n",
        "    norm = np.linalg.norm(class_prob, 2, axis=0)\n",
        "    class_prob = class_prob / norm\n",
        "    class_prob = (class_prob.T - class_prob.min(1)) / class_prob.ptp(1)\n",
        "    # convert to df\n",
        "    cp_df = pd.DataFrame(\n",
        "        class_prob, columns=lab.columns\n",
        "    )\n",
        "    cp_df.index = index\n",
        "    # classify as max score\n",
        "    m = cp_df.idxmax(axis=1)\n",
        "    \n",
        "    return m\n",
        "\n",
        "class_def = label_transfer(distances, adata_ref.obs.louvain, adata.obs.index)\n",
        "\n",
        "# add to obs section of the original object\n",
        "adata.obs['label_trans'] = class_def\n",
        "\n",
        "sc.pl.umap(adata, color=\"label_trans\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: plot-transfer\n",
        "# add to merged object.\n",
        "adata_merged.obs[\"label_trans\"] = pd.concat(\n",
        "    [class_def, adata_ref.obs[\"louvain\"]], axis=0\n",
        ").tolist()\n",
        "\n",
        "sc.pl.umap(adata_merged, color=[\"sample\",\"louvain\",'label_trans'])\n",
        "#plot only ctrl cells.\n",
        "sc.pl.umap(adata_merged[adata_merged.obs['sample']=='ctrl'], color='label_trans')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now plot how many cells of each celltypes can be found in each cluster."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: distributions\n",
        "tmp = pd.crosstab(adata.obs['leiden_0.6'],adata.obs['label_trans'], normalize='index')\n",
        "tmp.plot.bar(stacked=True).legend(bbox_to_anchor=(1.8, 1),loc='upper right')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Ingest\n",
        "\n",
        "Another method for celltype prediction is Ingest, for more information,\n",
        "please look at\n",
        "https://scanpy-tutorials.readthedocs.io/en/latest/integrating-data-using-ingest.html"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: ingest\n",
        "sc.tl.ingest(adata, adata_ref, obs='louvain')\n",
        "sc.pl.umap(adata, color=['louvain','leiden_0.6'], wspace=0.5)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As you can see, ingest has created a new umap for us, so to get\n",
        "consistent plotting, lets revert back to the old one for further\n",
        "plotting:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: old-umap\n",
        "adata.obsm[\"X_umap\"] = adata.obsm[\"X_umap_uncorr\"]\n",
        "\n",
        "sc.pl.umap(adata, color=['louvain','leiden_0.6'], wspace=0.5)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now plot how many cells of each celltypes can be found in each cluster."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: distributions-ingest\n",
        "tmp = pd.crosstab(adata.obs['leiden_0.6'],adata.obs['louvain'], normalize='index')\n",
        "tmp.plot.bar(stacked=True).legend(bbox_to_anchor=(1.8, 1),loc='upper right')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Celltypist\n",
        "\n",
        "[Celltypist](https://www.celltypist.org/) provides pretrained models for\n",
        "classification for many different human tissues and celltypes. Here, we\n",
        "are following the steps of this\n",
        "[tutorial](https://colab.research.google.com/github/Teichlab/celltypist/blob/main/docs/notebook/celltypist_tutorial.ipynb),\n",
        "with some adaptations for this dataset. So please check out the tutorial\n",
        "for more detail."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: fetch-celltypist\n",
        "import celltypist\n",
        "from celltypist import models\n",
        "\n",
        "# there are many different models, we will only download 2 of them for now.\n",
        "models.download_models(force_update = False, model = 'Immune_All_Low.pkl')\n",
        "models.download_models(force_update = False, model = 'Immune_All_High.pkl')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now select the model you want to use and show the info:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: celltypist-model\n",
        "model = models.Model.load(model = 'Immune_All_High.pkl')\n",
        "\n",
        "model"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To infer celltype labels to our cells, we first need to convert back to\n",
        "the full matrix. OBS! For celltypist we want to have log1p normalised\n",
        "expression to 10,000 counts per cell. Which we already have in\n",
        "`adata.raw.X`, check by summing up the data, it should sum to 10K."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: from-raw\n",
        "adata = adata.raw.to_adata() \n",
        "adata.X.expm1().sum(axis = 1)[:10]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: celltypist\n",
        "predictions = celltypist.annotate(adata, model = 'Immune_All_High.pkl', majority_voting = True)\n",
        "\n",
        "predictions.predicted_labels"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The first column `predicted_labels` is the predictions made for each\n",
        "individual cell, while `majority_voting` is done for local subclusters,\n",
        "the clustering identities are in column `over_clustering`.\n",
        "\n",
        "Now we convert the predictions to an anndata object."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: plot-celltypist\n",
        "adata = predictions.to_adata()\n",
        "\n",
        "sc.pl.umap(adata, color = ['leiden_0.6', 'predicted_labels', 'majority_voting'], legend_loc = 'on data')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div>\n",
        "\n",
        "> **Task**\n",
        ">\n",
        "> Rerun predictions with Celltypist, but use another model, for instance\n",
        "> `Immune_All_High.pkl`, or any other model you find relevant, you can\n",
        "> find a list of models [here](https://www.celltypist.org/models). How\n",
        "> do the results differ for you?\n",
        "\n",
        "</div>\n",
        "\n",
        "### Celltypist custom model\n",
        "\n",
        "We can also train our own model on any reference data that we want to\n",
        "use. In this case we will use the pbmc data in `adata_ref` to train a\n",
        "model.\n",
        "\n",
        "Celltypist requires the data to be in the format of log1p normalised\n",
        "expression to 10,000 counts per cell, we can check if that is the case\n",
        "for the object we have:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: prep-refdata\n",
        "adata_ref.raw.X.expm1().sum(axis = 1)[:10]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "These should all sum up to 10K, which is not the case, probably since\n",
        "some genes were removed after normalizing. Wo we will have to start from\n",
        "the raw counts of that dataset instead. Before we selected the data\n",
        "`pbmc3k_processed`, but now we will instead use `pbmc3k`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: prep-refdata2\n",
        "adata_ref2 = sc.datasets.pbmc3k() \n",
        "adata_ref2"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This data is not annotated, so we will have to match the indices from\n",
        "the filtered and processed object. And add in the metadata with\n",
        "annotations."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: prep-refdata3\n",
        "adata_ref2 = adata_ref2[adata_ref.obs_names,:]\n",
        "adata_ref2.obs = adata_ref.obs\n",
        "adata_ref2"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we can normalize the matrix:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: process-refdata\n",
        "sc.pp.normalize_total(adata_ref2, target_sum = 1e4)\n",
        "sc.pp.log1p(adata_ref2)\n",
        "\n",
        "# check the sums again\n",
        "adata_ref2.X.expm1().sum(axis = 1)[:10]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And finally train the model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: celltypist-train\n",
        "new_model = celltypist.train(adata_ref2, labels = 'louvain', n_jobs = 10, feature_selection = True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we can run predictions on our data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: celltypist-pred\n",
        "predictions2 = celltypist.annotate(adata, model = new_model, majority_voting = True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Instead of converting the predictions to anndata we will just add\n",
        "another column in the `adata.obs` with these new predictions since the\n",
        "column names from the previous celltypist runs with clash."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: celltypist-save\n",
        "adata.obs[\"predicted_labels_ref\"] = predictions2.predicted_labels[\"predicted_labels\"]\n",
        "adata.obs[\"majority_voting_ref\"] = predictions2.predicted_labels[\"majority_voting\"]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: celltypist-plot\n",
        "sc.pl.umap(adata, color = ['predicted_labels', 'majority_voting','predicted_labels_ref', 'majority_voting_ref'], legend_loc = 'on data', ncols=2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compare results\n",
        "\n",
        "The predictions from ingest is stored in the column 'louvain' while we\n",
        "named the label transfer with scanorama as 'predicted'"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: plot-all\n",
        "sc.pl.umap(adata, color=['louvain','label_trans','majority_voting', 'majority_voting_ref'], wspace=0.5, ncols=3)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As you can see, the main celltypes are generally the same, but there are\n",
        "clearly differences, especially with regards to the cells predicted as\n",
        "either ILC/NK/CD8 T-cells.\n",
        "\n",
        "The only way to make sure which method you trust is to look at what\n",
        "genes the different celltypes express and use your biological knowledge\n",
        "to make decisions.\n",
        "\n",
        "## Gene set analysis\n",
        "\n",
        "Another way of predicting celltypes is to use the differentially\n",
        "expressed genes per cluster and compare to lists of known cell marker\n",
        "genes. This requires a list of genes that you trust and that is relevant\n",
        "for the tissue you are working on.\n",
        "\n",
        "You can either run it with a marker list from the ontology or a list of\n",
        "your choice as in the example below."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: fetch-markers\n",
        "path_file = 'data/human_cell_markers.txt'\n",
        "if not os.path.exists(path_file):\n",
        "    file_url = os.path.join(path_data, \"misc/human_cell_markers.txt\")\n",
        "    subprocess.call([\"curl\", \"-u\", curl_upass, \"-o\", path_file, file_url ])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: read-markers\n",
        "df = pd.read_table(path_file)\n",
        "df\n",
        "\n",
        "print(df.shape)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: prep-markers\n",
        "# Filter for number of genes per celltype\n",
        "df['nG'] = df.geneSymbol.str.split(\",\").str.len()\n",
        "\n",
        "df = df[df['nG'] > 5]\n",
        "df = df[df['nG'] < 100]\n",
        "d = df[df['cancerType'] == \"Normal\"]\n",
        "print(df.shape)\n",
        "\n",
        "# convert to dict.\n",
        "df.index = df.cellName\n",
        "gene_dict = df.geneSymbol.str.split(\",\").to_dict()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: dge\n",
        "# run differential expression per cluster\n",
        "sc.tl.rank_genes_groups(adata, 'leiden_0.6', method='wilcoxon', key_added = \"wilcoxon\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: gsea\n",
        "# do gene set overlap to the groups in the gene list and top 300 DEGs.\n",
        "import gseapy\n",
        "\n",
        "gsea_res = dict()\n",
        "pred = dict()\n",
        "\n",
        "for cl in adata.obs['leiden_0.6'].cat.categories.tolist():\n",
        "    print(cl)\n",
        "    glist = sc.get.rank_genes_groups_df(adata, group=cl, key='wilcoxon')[\n",
        "        'names'].squeeze().str.strip().tolist()\n",
        "    enr_res = gseapy.enrichr(gene_list=glist[:300],\n",
        "                             organism='Human',\n",
        "                             gene_sets=gene_dict,\n",
        "                             background=adata.shape[1],\n",
        "                             cutoff=1)\n",
        "    if enr_res.results.shape[0] == 0:\n",
        "        pred[cl] = \"Unass\"\n",
        "    else:\n",
        "        enr_res.results.sort_values(\n",
        "            by=\"P-value\", axis=0, ascending=True, inplace=True)\n",
        "        print(enr_res.results.head(2))\n",
        "        gsea_res[cl] = enr_res\n",
        "        pred[cl] = enr_res.results[\"Term\"][0]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: gsea-res\n",
        "# prediction per cluster\n",
        "pred"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: plot-gsea\n",
        "prediction = [pred[x] for x in adata.obs['leiden_0.6']]\n",
        "adata.obs[\"GS_overlap_pred\"] = prediction\n",
        "\n",
        "sc.pl.umap(adata, color='GS_overlap_pred')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div>\n",
        "\n",
        "> **Discuss**\n",
        ">\n",
        "> As you can see, it agrees to some extent with the predictions from the\n",
        "> methods above, but there are clear differences, which do you think\n",
        "> looks better?\n",
        "\n",
        "</div>\n",
        "\n",
        "## Save data\n",
        "\n",
        "We can finally save the object for use in future steps."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: save\n",
        "adata.write_h5ad('data/covid/results/scanpy_covid_qc_dr_int_cl_ct-ctrl13.h5ad')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Session info\n",
        "\n",
        "```{=html}\n",
        "<details>\n",
        "```\n",
        "```{=html}\n",
        "<summary>\n",
        "```\n",
        "Click here\n",
        "```{=html}\n",
        "</summary>\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#| label: session\n",
        "sc.logging.print_versions()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "```{=html}\n",
        "</details>\n",
        "```"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}