do_DEG_call_celltypes.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a802d7e-4549-41b6-b815-4212dca83995",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scanpy as sc\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30f79f47-a697-4786-a125-f288eba9b67e",
   "metadata": {},
   "source": [
    "# Final STAG AML data after all cleaning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15b4f064-6caa-4c54-a037-619a74830f98",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = sc.read(\"STAG_AML_final.h5ad\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8211b035-bd70-471c-93e3-d60d05b598c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sc.tl.leiden(adata, resolution=1.2) more granularity to tease out LMPP probably"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ebfe98d-2c55-4070-b975-f7dac9caff28",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (10, 10)\n",
    "sc.pl.umap(\n",
    "    adata,\n",
    "    color=[\"leiden\"],\n",
    "    s=50,\n",
    "    frameon=False,\n",
    "    show=False,\n",
    "    legend_loc=\"on data\",\n",
    "    legend_fontsize=20,\n",
    "    legend_fontoutline=4,\n",
    ")\n",
    "plt.savefig(f\"STAG_AML_final_leiden.jpg\", bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30f25094-a6c4-4ad8-a1dc-f1441d727b89",
   "metadata": {},
   "source": [
    "# do DE on leiden clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f485f71b-7ba4-4380-a472-ae19e05cfa88",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.uns[\"log1p\"][\"base\"] = None\n",
    "sc.tl.rank_genes_groups(adata, \"leiden\", method=\"wilcoxon\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c65861a6-4fc0-47ad-90b0-dba7eba3bad0",
   "metadata": {},
   "outputs": [],
   "source": [
    "result = adata.uns[\"rank_genes_groups\"]\n",
    "groups = result[\"names\"].dtype.names\n",
    "deg = pd.DataFrame(\n",
    "    {\n",
    "        group + \"_\" + key[:4]: result[key][group]\n",
    "        for group in groups\n",
    "        for key in [\"names\", \"pvals_adj\", \"logfoldchanges\"]\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd3e1d1a-ac9f-4845-a507-54968d14b893",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 6\n",
    "deg_sub = deg[[f\"{i}_name\", f\"{i}_pval\", f\"{i}_logf\"]].copy()\n",
    "deg_sub[\"logpval_mult_logf\"] = -np.log10(deg_sub[f\"{i}_pval\"] + 10 ** (-300)) * np.abs(\n",
    "    deg_sub[f\"{i}_logf\"]\n",
    ")\n",
    "\n",
    "deg_sub = (\n",
    "    deg_sub[(deg_sub[f\"{i}_pval\"] < 0.001) & (deg_sub[f\"{i}_logf\"] > 1.5)]\n",
    "    .sort_values(by=f\"{i}_logf\", ascending=False)\n",
    "    .head(20)\n",
    ")\n",
    "deg_sub.head(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17da2513-578f-47cc-8387-1a888c03ff0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "genes = deg_sub[deg_sub.columns[0]].tolist()[:20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f597ebcb-b50a-4532-bbba-b5c2f0ed6d11",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (5, 5)\n",
    "sc.pl.umap(\n",
    "    adata,\n",
    "    color=genes,\n",
    "    s=20,\n",
    "    frameon=False,\n",
    "    show=False,\n",
    "    legend_loc=\"on data\",\n",
    "    legend_fontsize=10,\n",
    "    ncols=5,\n",
    ")\n",
    "plt.savefig(f\"STAG_AML_final_DEGs_leiden_c{i}.jpg\", bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e579a5e-d798-44e5-81c0-b3d17e7a0da4",
   "metadata": {},
   "source": [
    "# name leiden clusters based on DEGs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c80e64b-b5ce-4a7e-9d71-adc1bf771c44",
   "metadata": {},
   "outputs": [],
   "source": [
    "leiden_cell = {\n",
    "    \"0\": \"HSC\",\n",
    "    \"1\": \"Cycling\",\n",
    "    \"2\": \"GMP\",\n",
    "    \"3\": \"GMP\",\n",
    "    \"4\": \"HSC\",\n",
    "    \"5\": \"CD14+CD16+ Mono.\",\n",
    "    \"6\": \"Early Mono.\",\n",
    "    \"7\": \"Naive CD4_CD8 T\",\n",
    "    \"8\": \"Neutrophil\",\n",
    "    \"9\": \"HSC\",\n",
    "    \"10\": \"Neutrophil\",\n",
    "    \"11\": \"NK_Eff_Mem T\",\n",
    "    \"12\": \"CD16 Mono.\",\n",
    "    \"13\": \"MEP\",\n",
    "    \"14\": \"Erythrocytes\",\n",
    "    \"15\": \"Cycling\",\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c07ad27-8db7-42d7-8913-0d2a2c3d098c",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs[\"cell_type\"] = adata.obs.leiden.apply(lambda x: leiden_cell[x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c81e2d4-ab3f-424e-8054-d48a6b571acf",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (10, 10)\n",
    "sc.pl.umap(\n",
    "    adata,\n",
    "    color=[\"cell_type\"],\n",
    "    s=50,\n",
    "    frameon=False,\n",
    "    show=False,\n",
    "    legend_loc=\"on data\",\n",
    "    legend_fontsize=20,\n",
    "    legend_fontoutline=4,\n",
    ")\n",
    "plt.savefig(f\"STAG_AML_final_cell_types.jpg\", bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5116d071-0df1-4956-8745-4df47f132037",
   "metadata": {},
   "source": [
    "# do DE on cell_type clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8a0a43a-3eb4-49c5-b861-3e4166932cc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.uns[\"log1p\"][\"base\"] = None\n",
    "sc.tl.rank_genes_groups(\n",
    "    adata, \"cell_type\", method=\"wilcoxon\", key_added=\"cell_type_DEGs\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c23976c7-25f9-4284-9ffd-b3ad8bcc7b11",
   "metadata": {},
   "outputs": [],
   "source": [
    "result = adata.uns[\"cell_type_DEGs\"]\n",
    "groups = result[\"names\"].dtype.names\n",
    "deg = pd.DataFrame(\n",
    "    {\n",
    "        group + \"_\" + key[:4]: result[key][group]\n",
    "        for group in groups\n",
    "        for key in [\"names\", \"pvals_adj\", \"logfoldchanges\"]\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4357ccfd-4390-4bd3-b65e-ff85cfe8a9b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_types = adata.obs.cell_type.cat.categories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c67970e3-bf65-4b1a-ae46-9ac9c2f07929",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(cell_types), cell_types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02c645de-0ac7-4786-a297-f24c434f4ce7",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.write_h5ad(\"STAG_AML_annotated_DEGs.h5ad\", compression=\"gzip\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b569e8d2-ab7f-40bd-bd79-a1e78fa75287",
   "metadata": {},
   "source": [
    "# aggregate all cell_type DEG in one dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2412e14-74d4-4075-ba34-77374e00e305",
   "metadata": {},
   "outputs": [],
   "source": [
    "ALL_DEGs = {}\n",
    "for i in range(len(cell_types)):\n",
    "    i = cell_types[i]\n",
    "    deg_sub = deg[[f\"{i}_name\", f\"{i}_pval\", f\"{i}_logf\"]].copy()\n",
    "    deg_sub[\"logpval_mult_logf\"] = -np.log10(\n",
    "        deg_sub[f\"{i}_pval\"] + 10 ** (-300)\n",
    "    ) * np.abs(deg_sub[f\"{i}_logf\"])\n",
    "\n",
    "    deg_sub = deg_sub[\n",
    "        (deg_sub[f\"{i}_pval\"] < 0.000001) & (deg_sub[f\"{i}_logf\"] > 1.5)\n",
    "    ].sort_values(by=f\"{i}_logf\", ascending=False)\n",
    "    ALL_DEGs[i] = deg_sub.head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adf3f46f-5a87-4c07-a27c-4341c713fa15",
   "metadata": {},
   "source": [
    "# plot tops DEGs of each cell_type cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "105662d9-c450-49a1-af43-42c9af4e5029",
   "metadata": {},
   "outputs": [],
   "source": [
    "for c in range(len(cell_types)):\n",
    "    col = 5\n",
    "    row = 4\n",
    "    wid = 5\n",
    "    deg_sub = ALL_DEGs[cell_types[c]]\n",
    "    genes = deg_sub[deg_sub.columns[0]].tolist()\n",
    "    fig, ax = plt.subplots(\n",
    "        row,\n",
    "        col,\n",
    "        figsize=(col * wid, row * wid),\n",
    "        gridspec_kw={\"wspace\": 0.01, \"hspace\": 0.1},\n",
    "    )\n",
    "    axr = ax.ravel()\n",
    "    for i, g in enumerate(genes):\n",
    "        sc.pl.umap(adata, color=g, frameon=False, s=10, ax=axr[i], show=False)\n",
    "        axr[i].set_title(f\"{g}\", fontsize=15)\n",
    "        cbar = axr[i].collections[0].colorbar\n",
    "        if cbar != None:\n",
    "            cbar.remove()\n",
    "    for j in np.arange(i + 1, len(axr)):\n",
    "        axr[j].remove()\n",
    "    plt.savefig(f\"STAG_AML_final_DEGs_{cell_types[c]}.jpg\", bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "514beeae-a0b2-4a00-86db-2e79d58e19d7",
   "metadata": {},
   "source": [
    "# plot tops DEGs from the entire AML dataset on top of our dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "923af4be-a13d-4db8-8291-9bcddf5ce343",
   "metadata": {},
   "outputs": [],
   "source": [
    "deg_full_data = pd.read_csv(\"../all_csvs/LMPP_DEG.csv\", index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7b176f2-149e-4ddf-be70-66cf36754cc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "deg_full_data = pd.read_csv(\"../all_csvs/GMP_DEG.csv\" ,index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "624f1cdc-915b-4af6-9bde-0530faaf670a",
   "metadata": {},
   "outputs": [],
   "source": [
    "deg_full_data = pd.read_csv(\"../all_csvs/Prog_Mk_DEG.csv\", index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5865ffe4-d1ca-4555-bc59-ac4df5e3cb9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "genes = adata.var[adata.var.index.isin(deg_full_data.index[:5])].index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f002c3fe-10d5-44e6-8768-714956d55c51",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (4, 4)\n",
    "sc.pl.umap(adata, color=genes, s=20, frameon=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b818b875-a2fc-437e-91fb-2790dc0b16b9",
   "metadata": {},
   "source": [
    "# find which clusters have genes of interest as DEG in the entire AML dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b31956b3-7d9d-4270-9bba-31ebae02d724",
   "metadata": {},
   "outputs": [],
   "source": [
    "deg_full_data = pd.read_csv(\"../all_csvs/AML_predicted_celltype_DEG.csv\", index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01a7460b-7d1e-48f2-a8a7-c71af4398a60",
   "metadata": {},
   "outputs": [],
   "source": [
    "genes = deg_sub[deg_sub.columns[0]].tolist()[:40]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01ebd8b5-dd29-44b9-8d65-255d77ffb7ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "for gene in genes:\n",
    "    result = deg_full_data.eq(gene).stack()\n",
    "    print(22 * \"-\", gene, 20 * \"-\")\n",
    "\n",
    "    if result[result].shape[0] > 0:\n",
    "        for i in range(7):  # print up to 7 clusters\n",
    "            row_index, column_index = result[result].index[i]\n",
    "            values = deg_full_data.loc[\n",
    "                row_index,\n",
    "                [\n",
    "                    column_index.split(\"_name\")[0] + \"_pval\",\n",
    "                    column_index.split(\"_name\")[0] + \"_logf\",\n",
    "                ],\n",
    "            ].values\n",
    "            logfold = values[1]\n",
    "            pval = values[0]\n",
    "\n",
    "            if pval < 0.0001 and logfold > 1.2:\n",
    "                print(np.round(logfold, 2), column_index, np.round(pval, 5))"
   ]
  }
 ],
 "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}