{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scanpy as sc # v1.6\",\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "import os.path\n",
    "import anndata\n",
    "from matplotlib import rcParams\n",
    "import seaborn as sns\n",
    "import numba\n",
    "import mnnpy\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import full COVID-19 PBMC dataset\n",
    "os.chdir('/home/ngr18/covid/')\n",
    "covid_total = sc.read_h5ad('covid.h5ad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Subset PBMC data to myeloid populations and reorder categories (for dotplot visualisations)\n",
    "\n",
    "blood_myeloid = covid_total[covid_total.obs.full_clustering.isin(['CD83_CD14_mono', 'CD14_mono', \n",
    "                                                                  'CD16_mono', 'C1_CD16_mono',\n",
    "                                                                 'DC1', 'DC2', 'DC3', 'ASDC', 'pDC', 'DC_prolif',\n",
    "                                                                 'Mono_prolif']),:]\n",
    "\n",
    "blood_myeloid.obs.full_clustering = blood_myeloid.obs.full_clustering.cat.reorder_categories([\n",
    "'DC1', 'DC2', 'DC3', 'ASDC', 'pDC', 'DC_prolif', \n",
    "    'CD83_CD14_mono', 'CD14_mono', 'CD16_mono', 'C1_CD16_mono', 'Mono_prolif'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import BAL data GSE145926 (reannotated for DC subsets)\n",
    "\n",
    "os.chdir('/home/ngr18/covid/external_dataset/BAL')\n",
    "bal = sc.read_h5ad('full_bal.h5ad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Subset BAL data to myeloid populations\n",
    "\n",
    "bal_myeloid = bal[bal.obs.full_clustering.isin(['DC1', 'DC2', 'mat_DC', 'pDC', 'Mac', 'Prolif_mac']),:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Concatenate PBMC and BAL data\n",
    "\n",
    "myeloid = anndata.concat([bal_myeloid, blood_myeloid], index_unique = None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Normalise and scale complete object\n",
    "\n",
    "sc.pp.normalize_total(myeloid, target_sum=1e4)\n",
    "sc.pp.log1p(myeloid)\n",
    "sc.pp.scale(myeloid, max_value=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Add scores\n",
    "\n",
    "#https://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_TNFA_SIGNALING_VIA_NFKB.html\n",
    "TNF_genes = [\"ABCA1\", \t\"ACKR3\", \t\"AREG\", \t\"ATF3\", \t\"ATP2B1\", \t\"B4GALT1\", \t\"B4GALT5\", \t\"BCL2A1\", \t\"BCL3\", \t\"BCL6\", \t\"BHLHE40\", \t\"BIRC2\", \t\"BIRC3\", \t\"BMP2\", \t\"BTG1\", \t\"BTG2\", \t\"BTG3\", \t\"CCL2\", \t\"CCL20\", \t\"CCL4\", \t\"CCL5\", \t\"CCN1\", \t\"CCND1\", \t\"CCNL1\", \t\"CCRL2\", \t\"CD44\", \t\"CD69\", \t\"CD80\", \t\"CD83\", \t\"CDKN1A\", \t\"CEBPB\", \t\"CEBPD\", \t\"CFLAR\", \t\"CLCF1\", \t\"CSF1\", \t\"CSF2\", \t\"CXCL1\", \t\"CXCL10\", \t\"CXCL11\", \t\"CXCL2\", \t\"CXCL3\", \t\"CXCL6\", \t\"DDX58\", \t\"DENND5A\", \t\"DNAJB4\", \t\"DRAM1\", \t\"DUSP1\", \t\"DUSP2\", \t\"DUSP4\", \t\"DUSP5\", \t\"EDN1\", \t\"EFNA1\", \t\"EGR1\", \t\"EGR2\", \t\"EGR3\", \t\"EHD1\", \t\"EIF1\", \t\"ETS2\", \t\"F2RL1\", \t\"F3\", \t\"FJX1\", \t\"FOS\", \t\"FOSB\", \t\"FOSL1\", \t\"FOSL2\", \t\"FUT4\", \t\"G0S2\", \t\"GADD45A\", \t\"GADD45B\", \t\"GCH1\", \t\"GEM\", \t\"GFPT2\", \t\"GPR183\", \t\"HBEGF\", \t\"HES1\", \t\"ICAM1\", \t\"ICOSLG\", \t\"ID2\", \t\"IER2\", \t\"IER3\", \t\"IER5\", \t\"IFIH1\", \t\"IFIT2\", \t\"IFNGR2\", \t\"IL12B\", \t\"IL15RA\", \t\"IL18\", \t\"IL1A\", \t\"IL1B\", \t\"IL23A\", \t\"IL6\", \t\"IL6ST\", \t\"IL7R\", \t\"INHBA\", \t\"IRF1\", \t\"IRS2\", \t\"JAG1\", \t\"JUN\", \t\"JUNB\", \t\"KDM6B\", \t\"KLF10\", \t\"KLF2\", \t\"KLF4\", \t\"KLF6\", \t\"KLF9\", \t\"KYNU\", \t\"LAMB3\", \t\"LDLR\", \t\"LIF\", \t\"LITAF\", \t\"MAFF\", \t\"MAP2K3\", \t\"MAP3K8\", \t\"MARCKS\", \t\"MCL1\", \t\"MSC\", \t\"MXD1\", \t\"MYC\", \t\"NAMPT\", \t\"NFAT5\", \t\"NFE2L2\", \t\"NFIL3\", \t\"NFKB1\", \t\"NFKB2\", \t\"NFKBIA\", \t\"NFKBIE\", \t\"NINJ1\", \t\"NR4A1\", \t\"NR4A2\", \t\"NR4A3\", \t\"OLR1\", \t\"PANX1\", \t\"PDE4B\", \t\"PDLIM5\", \t\"PER1\", \t\"PFKFB3\", \t\"PHLDA1\", \t\"PHLDA2\", \t\"PLAU\", \t\"PLAUR\", \t\"PLEK\", \t\"PLK2\", \t\"PLPP3\", \t\"PMEPA1\", \t\"PNRC1\", \t\"PPP1R15A\", \t\"PTGER4\", \t\"PTGS2\", \t\"PTPRE\", \t\"PTX3\", \t\"RCAN1\", \t\"REL\", \t\"RELA\", \t\"RELB\", \t\"RHOB\", \t\"RIPK2\", \t\"RNF19B\", \t\"SAT1\", \t\"SDC4\", \t\"SERPINB2\", \t\"SERPINB8\", \t\"SERPINE1\", \t\"SGK1\", \t\"SIK1\", \t\"SLC16A6\", \t\"SLC2A3\", \t\"SLC2A6\", \t\"SMAD3\", \t\"SNN\", \t\"SOCS3\", \t\"SOD2\", \t\"SPHK1\", \t\"SPSB1\", \t\"SQSTM1\", \t\"STAT5A\", \t\"TANK\", \t\"TAP1\", \t\"TGIF1\", \t\"TIPARP\", \t\"TLR2\", \t\"TNC\", \t\"TNF\", \t\"TNFAIP2\", \t\"TNFAIP3\", \t\"TNFAIP6\", \t\"TNFAIP8\", \t\"TNFRSF9\", \t\"TNFSF9\", \t\"TNIP1\", \t\"TNIP2\", \t\"TRAF1\", \t\"TRIB1\", \t\"TRIP10\", \t\"TSC22D1\", \t\"TUBB2A\", \t\"VEGFA\", \t\"YRDC\", \t\"ZBTB10\", \t\"ZC3H12A\", \t\"ZFP36\"]\n",
    "\n",
    "#https://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_IL6_JAK_STAT3_SIGNALING\n",
    "JAK_genes = [\"A2M\", \t\"ACVR1B\", \t\"ACVRL1\", \t\"BAK1\", \t\"CBL\", \t\"CCL7\", \t\"CCR1\", \t\"CD14\", \t\"CD36\", \t\"CD38\", \t\"CD44\", \t\"CD9\", \t\"CNTFR\", \t\"CRLF2\", \t\"CSF1\", \t\"CSF2\", \t\"CSF2RA\", \t\"CSF2RB\", \t\"CSF3R\", \t\"CXCL1\", \t\"CXCL10\", \t\"CXCL11\", \t\"CXCL13\", \t\"CXCL3\", \t\"CXCL9\", \t\"DNTT\", \t\"EBI3\", \t\"FAS\", \t\"GRB2\", \t\"HAX1\", \t\"HMOX1\", \t\"IFNAR1\", \t\"IFNGR1\", \t\"IFNGR2\", \t\"IL10RB\", \t\"IL12RB1\", \t\"IL13RA1\", \t\"IL15RA\", \t\"IL17RA\", \t\"IL17RB\", \t\"IL18R1\", \t\"IL1B\", \t\"IL1R1\", \t\"IL1R2\", \t\"IL2RA\", \t\"IL2RG\", \t\"IL3RA\", \t\"IL4R\", \t\"IL6\", \t\"IL6ST\", \t\"IL7\", \t\"IL9R\", \t\"INHBE\", \t\"IRF1\", \t\"IRF9\", \t\"ITGA4\", \t\"ITGB3\", \t\"JUN\", \t\"LEPR\", \t\"LTB\", \t\"LTBR\", \t\"MAP3K8\", \t\"MYD88\", \t\"OSMR\", \t\"PDGFC\", \t\"PF4\", \t\"PIK3R5\", \t\"PIM1\", \t\"PLA2G2A\", \t\"PTPN1\", \t\"PTPN11\", \t\"PTPN2\", \t\"REG1A\", \t\"SOCS1\", \t\"SOCS3\", \t\"STAM2\", \t\"STAT1\", \t\"STAT2\", \t\"STAT3\", \t\"TGFB1\", \t\"TLR2\", \t\"TNF\", \t\"TNFRSF12A\", \t\"TNFRSF1A\", \t\"TNFRSF1B\", \t\"TNFRSF21\", \t\"TYK2\"]\n",
    "\n",
    "#https://www.gsea-msigdb.org/gsea/msigdb/cards/GO_RESPONSE_TO_TYPE_I_INTERFERON\n",
    "IFN_genes = [\"ABCE1\", \t\"ADAR\", \t\"BST2\", \t\"CACTIN\", \t\"CDC37\", \t\"CNOT7\", \t\"DCST1\", \t\"EGR1\", \t\"FADD\", \t\"GBP2\", \t\"HLA-A\", \t\"HLA-B\", \t\"HLA-C\", \t\"HLA-E\", \t\"HLA-F\", \t\"HLA-G\", \t\"HLA-H\", \t\"HSP90AB1\", \t\"IFI27\", \t\"IFI35\", \t\"IFI6\", \t\"IFIT1\", \t\"IFIT2\", \t\"IFIT3\", \t\"IFITM1\", \t\"IFITM2\", \t\"IFITM3\", \t\"IFNA1\", \t\"IFNA10\", \t\"IFNA13\", \t\"IFNA14\", \t\"IFNA16\", \t\"IFNA17\", \t\"IFNA2\", \t\"IFNA21\", \t\"IFNA4\", \t\"IFNA5\", \t\"IFNA6\", \t\"IFNA7\", \t\"IFNA8\", \t\"IFNAR1\", \t\"IFNAR2\", \t\"IFNB1\", \t\"IKBKE\", \t\"IP6K2\", \t\"IRAK1\", \t\"IRF1\", \t\"IRF2\", \t\"IRF3\", \t\"IRF4\", \t\"IRF5\", \t\"IRF6\", \t\"IRF7\", \t\"IRF8\", \t\"IRF9\", \t\"ISG15\", \t\"ISG20\", \t\"JAK1\", \t\"LSM14A\", \t\"MAVS\", \t\"METTL3\", \t\"MIR21\", \t\"MMP12\", \t\"MUL1\", \t\"MX1\", \t\"MX2\", \t\"MYD88\", \t\"NLRC5\", \t\"OAS1\", \t\"OAS2\", \t\"OAS3\", \t\"OASL\", \t\"PSMB8\", \t\"PTPN1\", \t\"PTPN11\", \t\"PTPN2\", \t\"PTPN6\", \t\"RNASEL\", \t\"RSAD2\", \t\"SAMHD1\", \t\"SETD2\", \t\"SHFL\", \t\"SHMT2\", \t\"SP100\", \t\"STAT1\", \t\"STAT2\", \t\"TBK1\", \t\"TREX1\", \t\"TRIM56\", \t\"TRIM6\", \t\"TTLL12\", \t\"TYK2\", \t\"UBE2K\", \t\"USP18\", \t\"WNT5A\", \t\"XAF1\", \t\"YTHDF2\", \t\"YTHDF3\", \t\"ZBP1\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Score genes\n",
    "\n",
    "sc.tl.score_genes(myeloid, TNF_genes, ctrl_size=50, gene_pool=None, n_bins=25, score_name='TNF_score', random_state=0, copy=False, use_raw=None)\n",
    "    \n",
    "sc.tl.score_genes(myeloid, JAK_genes, ctrl_size=50, gene_pool=None, n_bins=25, score_name='JAK_score', random_state=0, copy=False, use_raw=None)\n",
    "        \n",
    "sc.tl.score_genes(myeloid, IFN_genes, ctrl_size=50, gene_pool=None, n_bins=25, score_name='IFN_score', random_state=0, copy=False, use_raw=None)\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Export for visualisation in R\n",
    "\n",
    "myeloid.obs[['full_clustering', 'severity', 'TNF_score', 'JAK_score', 'IFN_score']].to_csv('TNF_IFN_JAK_scores_v2.csv')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "covid_py",
   "language": "python",
   "name": "covid_py"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
