{

 "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": [

    "#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 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": [

    "#Concatenate PBMC and BAL data for visualisations\n",

    "\n",

    "myeloid = anndata.concat([bal_myeloid, blood_myeloid], index_unique = None)"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": null,

   "metadata": {},

   "outputs": [],

   "source": [

    "##myeloid figure dotplot - fig 2A (left - RNA)\n",

    "\n",

    "sc.pl.DotPlot(myeloid, [\n",

    "'CLEC9A', 'CADM1', 'CLEC10A','CD1C', 'CD14', 'VCAN',\n",

    "    'CCR7', 'LAMP3',\n",

    "                            'AXL', 'SIGLEC6',\n",

    "                            'LILRA4', 'ITM2C', 'GZMB',\n",

    "    'IL1B', 'IER3', 'LDLR', 'CD83',\n",

    "                         'S100A12', 'CSF3R', \n",

    "                        'FCGR3A', 'MS4A7', 'LILRB1', 'CSF1R', 'CDKN1C',\n",

    "                         'C1QA', 'C1QB', 'C1QC', 'CCR1',\n",

    "    \n",

    "    'MARCO', 'MKI67', 'TOP2A'], \n",

    "              log = True, \n",

    "              groupby='full_clustering').style(cmap='Blues',dot_edge_color='black', dot_edge_lw=1).swap_axes(False).show(True)\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": null,

   "metadata": {},

   "outputs": [],

   "source": [

    "##myeloid figure dotplot - fig 2A (right - protein) - using only PBMC data (no CITEseq data for BAL)\n",

    "\n",

    "sc.pl.DotPlot(blood_myeloid, ['AB_CD141', 'AB_CLEC9A', 'AB_KIT','AB_BTLA',\n",

    "                            'AB_CD1C', 'AB_CD101', 'AB_FcERIa',\n",

    "                            'AB_CD5', \n",

    "                            'AB_CD123', 'AB_CD45RA', 'AB_CD304',\n",

    "                            'AB_CD14', 'AB_CD99', 'AB_CD64',  \n",

    "                            'AB_CR1', 'AB_ITGAM',\n",

    "                            'AB_CD16', 'AB_C5AR1',\n",

    "                              'AB_CX3CR1'], \n",

    "              log = True, \n",

    "              groupby='full_clustering',  expression_cutoff=0.15).style(cmap='YlOrRd',dot_edge_color='black', dot_edge_lw=1).swap_axes(False).show(True)\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": null,

   "metadata": {},

   "outputs": [],

   "source": [

    "#subset data again to blood monocyte and BAL macrophages\n",

    "\n",

    "monos = myeloid[myeloid.obs.full_clustering.isin(['CD83_CD14_mono', 'CD14_mono', 'CD16_mono', 'C1_CD16_mono', \n",

    "                                                 'bal_Mac']),:]"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": null,

   "metadata": {},

   "outputs": [],

   "source": [

    "#Figure 2C left/upper - analysis of healthy only\n",

    "healthy = monos[monos.obs.Status_on_day_collection_summary.isin(['bal_healthy', 'Healthy']),:]\n",

    "\n",

    "sc.pp.normalize_total(healthy, target_sum=1e4)\n",

    "\n",

    "sc.pp.log1p(healthy)\n",

    "\n",

    "sc.pp.highly_variable_genes(healthy, n_top_genes = 3000, flavor = 'seurat')\n",

    "\n",

    "healthy.raw = healthy\n",

    "\n",

    "healthy = healthy[:, healthy.var.highly_variable]\n",

    "\n",

    "sc.pp.scale(healthy, max_value=10)\n",

    "\n",

    "sc.tl.pca(healthy, svd_solver='arpack')\n",

    "\n",

    "sc.external.pp.harmony_integrate(healthy, 'sample_id', basis='X_pca', adjusted_basis='X_pca_harmony')\n",

    "\n",

    "sc.pp.neighbors(healthy, n_neighbors=10, n_pcs=50, use_rep = 'X_pca_harmony')\n",

    "\n",

    "sc.tl.paga(healthy, groups='full_clustering')\n",

    "\n",

    "#recolor clusters for consistency with bar graph\n",

    "\n",

    "healthy.uns['full_clustering_colors'][0] = '#76DDC9'\n",

    "healthy.uns['full_clustering_colors'][1] = '#E28686'\n",

    "healthy.uns['full_clustering_colors'][2] = '#C2A3E2'\n",

    "healthy.uns['full_clustering_colors'][3] = '#991111'\n",

    "healthy.uns['full_clustering_colors'][4] = '#ffff00'\n",

    "\n",

    "\n",

    "sc.pl.paga(test, color=['full_clustering'], threshold = 0.11, save = 'healthy_myeloid_PAGA.pdf',\n",

    "           fontsize = 0, node_size_scale = 10, min_edge_width = 5, frameon = False)"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": null,

   "metadata": {},

   "outputs": [],

   "source": [

    "#Figure 2C left/lower - limit to only covid samples\n",

    "covid_mono = covid_mono[covid_mono.obs.Status_on_day_collection_summary.isin(['Asymptomatic', 'Critical', \n",

    "                                                            'Mild', 'Moderate', 'Severe',\n",

    "                                                            'bal_mild', 'bal_severe']),:]\n",

    "\n",

    "sc.pp.normalize_total(covid_mono, target_sum=1e4)\n",

    "\n",

    "sc.pp.log1p(covid_mono)\n",

    "\n",

    "sc.pp.highly_variable_genes(covid_mono, n_top_genes = 3000, flavor = 'seurat')\n",

    "\n",

    "covid_mono.raw = covid_mono\n",

    "\n",

    "covid_mono = covid_mono[:, covid_mono.var.highly_variable]\n",

    "\n",

    "sc.pp.scale(covid_mono, max_value=10)\n",

    "\n",

    "sc.tl.pca(covid_mono, svd_solver='arpack')\n",

    "\n",

    "sc.external.pp.harmony_integrate(covid_mono, 'sample_id', basis='X_pca', adjusted_basis='X_pca_harmony')\n",

    "\n",

    "sc.pp.neighbors(covid_mono, n_neighbors=10, n_pcs=50, use_rep = 'X_pca_harmony')\n",

    "\n",

    "sc.tl.paga(covid_mono, groups='full_clustering')\n",

    "\n",

    "#Again, recolor for consistency with bar graph\n",

    "\n",

    "covid_mono.uns['full_clustering_colors'][0] = '#76DDC9'\n",

    "covid_mono.uns['full_clustering_colors'][1] = '#E28686'\n",

    "covid_mono.uns['full_clustering_colors'][2] = '#C2A3E2'\n",

    "covid_mono.uns['full_clustering_colors'][3] = '#991111'\n",

    "covid_mono.uns['full_clustering_colors'][4] = '#ffff00'\n",

    "\n",

    "\n",

    "sc.pl.paga(covid_mono, color=['full_clustering'], threshold = 0.2, save = 'covid_myeloid_PAGA.pdf',\n",

    "           fontsize = 0, node_size_scale = 10, min_edge_width = 5, frameon = False)"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": null,

   "metadata": {},

   "outputs": [],

   "source": [

    "#Pseudotime plot\n",

    "\n",

    "covid_mono.uns['iroot'] = np.flatnonzero(covid_mono.obs['full_clustering']  == 'CD83_CD14_mono')[0]\n",

    "\n",

    "\n",

    "sc.pp.log1p(covid_mono)\n",

    "sc.pp.scale(covid_mono)\n",

    "           \n",

    "           gene_names = ['NFKBIA',  'KLF6', 'VIM', 'CD14', 'S100A8',\n",

    "    'HLA-DPA1', 'HLA-DPB1', 'FCGR3A','CTSC', 'CTSL', \n",

    "    'CCR1','RSAD2', 'CCL2', 'CXCL10', 'CCL7', 'TNFSF10']                              # monocyte\n",

    "\n",

    "paths = [('erythrocytes', ['mono1', 'mono2', 'mono3', 'mono4', 'bal_Mac']),\n",

    "             ('neutrophils', ['mono1', 'mono2', 'mono3', 'mono4', 'bal_Mac']),\n",

    "             ('monocytes', ['mono1', 'mono2', 'mono3', 'mono4', 'bal_Mac'])]\n",

    "\n",

    "test.obs['distance'] = test.obs['dpt_pseudotime']\n",

    "\n",

    "_, axs = plt.subplots(ncols=3, figsize=(6, 5), gridspec_kw={'wspace': 0.05, 'left': 0.12})\n",

    "plt.subplots_adjust(left=0.05, right=0.98, top=0.82, bottom=0.2)\n",

    "for ipath, (descr, path) in enumerate(paths):\n",

    "    _, data = sc.pl.paga_path(\n",

    "            test, path, gene_names,\n",

    "            show_node_names=False,\n",

    "            ax=axs[ipath],\n",

    "            ytick_fontsize=12,\n",

    "            left_margin=0.15,\n",

    "            n_avg=50,\n",

    "            annotations=['distance'],\n",

    "            show_yticks=True if ipath==0 else False,\n",

    "            show_colorbar=False,\n",

    "            color_map='coolwarm',\n",

    "            color_maps_annotations={'distance': 'viridis'},\n",

    "            title='{} path'.format(descr),\n",

    "            return_data=True,\n",

    "            show=False, normalize_to_zero_one = True)"

   ]

  }

 ],

 "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

}