CrossSectional-Flare-Study.V3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scanpy as sc\n",
    "import anndata as ad\n",
    "from scipy.sparse import csr_matrix\n",
    "import logging\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import csv\n",
    "import scipy.stats as stats\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.api as sm\n",
    "import math\n",
    "import gc\n",
    "import random\n",
    "from sklearn.utils import shuffle\n",
    "import combat\n",
    "import patsy\n",
    "logging.basicConfig(level=logging.INFO)\n",
    "savepath = 'Lupus_study_adjusted_V3.h5ad'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "##################\n",
    "# Configure file #\n",
    "##################\n",
    "sc.settings.verbosity = 2\n",
    "sc.settings.autoshow = False\n",
    "sc.settings.set_figure_params(dpi=100, dpi_save=300, format='png', frameon=False, transparent=True, fontsize=16)\n",
    "plt.rcParams[\"image.aspect\"] = \"equal\"\n",
    "adata = sc.read(savepath, cache=True)\n",
    "print(adata)\n",
    "figdir = \"./figures.\"\n",
    "sc.settings.figdir = \"./figures.\"\n",
    "MASTERCOLORS = sc.pl.palettes.default_64\n",
    "MASTERCOLORS.remove(\"#FEFFE6\")\n",
    "MASTERCOLORS.remove(\"#FFFF00\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "disease_umap = sc.pl.umap(adata, color=['disease_cov', 'SLE status'], size=1, show=True, edgecolor=\"none\", palette=sc.pl.palettes.vega_20_scanpy, save='.disease.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.embedding_density(adata, basis='umap', groupby='disease_cov')\n",
    "sc.pl.embedding_density(adata, basis='umap', key='umap_density_disease_cov', show=True,  save='umap_density_disease_cov.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.embedding_density(adata, basis='umap', groupby='disease_cov')\n",
    "sc.pl.embedding_density(adata, basis='umap', key='umap_density_disease_cov', show=True,  save='umap_density_disease_cov.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "disease_umap = sc.pl.umap(adata, color='batch_cov', size=3,show=True, edgecolor=\"none\", palette=sc.pl.palettes.vega_20_scanpy, save='.batch.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot individual variability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Plot 16 plots for the individuals.\n",
    "colors = [\"#8ED1C6\",\"#FCF6B5\", \"#BEBAD9\", \"#F47F72\", \"#81B1D3\", \"#FBB463\", \"#B4D66C\", \"#F9CEE1\", \"#DAD9D9\", \"#BC80B7\", \"#CDE6C4\", \"#FEEC6E\", \"#E31F26\", \"#387EB9\", \"#4EAF49\", \"#984F9F\", \"#8ED1C6\",\"#FCF6B5\", \"#BEBAD9\"];\n",
    "fig,ax = plt.subplots(nrows=4,ncols=4, figsize=(7,7),sharex=True, sharey=True)\n",
    "fig.tight_layout()\n",
    "plt.subplots_adjust(wspace=-0.1, hspace=0)\n",
    "\n",
    "batch_name = adata.obs.batch_cov.unique()[0]\n",
    "batch = adata.obs.ind_cov_disease_cov[adata.obs.batch_cov==batch_name].unique();\n",
    "batch = np.sort(batch.categories.values)\n",
    "for ind_i in list(range(16)):\n",
    "    ind = batch[ind_i]\n",
    "    col = colors[ind_i]\n",
    "    sc.pl.umap(adata[adata.obs.ind_cov_disease_cov==ind,], color = \"ind_cov_disease_cov\", title=None, palette=[col,col], ax=ax[(ind_i-1)%4,math.floor(ind_i/4)], size=20, edgecolor=\"none\")\n",
    "    ax[(ind_i-1)%4,math.floor(ind_i/4)].get_xaxis().set_visible(False)\n",
    "    ax[(ind_i-1)%4,math.floor(ind_i/4)].get_yaxis().set_visible(False)\n",
    "    ax[(ind_i-1)%4,math.floor(ind_i/4)].get_legend().remove()\n",
    "    ax[(ind_i-1)%4,math.floor(ind_i/4)].set_title(\"\")\n",
    "    ax[(ind_i-1)%4,math.floor(ind_i/4)].set_aspect(\"equal\")\n",
    "\n",
    "fig.savefig(figdir+'/'+batch_name+'.ind.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rank genes and plot leiden groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')\n",
    "sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color='leiden', size=3, show=True, edgecolor=\"none\", save='.leiden.png', palette=MASTERCOLORS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Most expressed genes per leiden group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option('display.max_columns', None)\n",
    "unique_leiden = np.unique(adata.obs['leiden'].values)\n",
    "# Compile list of top genes\n",
    "GeneRanks = pd.DataFrame()\n",
    "for ii in range(len(unique_leiden)):\n",
    "    GeneRanks[str('leiden_' + str(ii))] = adata.var_names[np.flipud(np.argsort(np.mean(adata.X[adata.obs['leiden'] == str(ii)], axis=0)))]\n",
    "GeneRanks.to_csv('Flare_study_top_expression.csv')\n",
    "GeneRanks.head(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs['ct_cov'] = adata.obs['ct_cov'].astype('object')\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"0\"] = \"Classical Monocytes\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"1\"] = \"Naive T Cells\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"2\"] = \"Cytotoxic T Cells\" ## good                                                                                     \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"3\"] = \"B Cells\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"4\"] = \"Effector Memory T Cells\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"5\"] = \"Natural Killer Cells\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"6\"] = \"Naive T Cells\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"7\"] = \"Nonclassical Monocytes\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"8\"] = \"Classical Dendritic Cells\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"9\"] = \"Naive T Cells\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"10\"] = \"Megakaryocytes\" ## good                                                                                  \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"11\"] = \"Plasmacytoid Dendritic Cells\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"12\"] = \"Nonclassical Monocytes\" ## good                                                                                   \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"13\"] = \"Proliferating T Cells\" ## good  Macro                                                                                 \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"14\"] = \"Naive T Cells\" ## good                                                                                  \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"15\"] = \"RBCs\" ## good                                                                                   \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"16\"] = \"B Cells\" ## good                                                                                   \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"17\"] = \"Classical Monocytes\" ## good                                                                                  \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"18\"] = \"Progenitor Cells\" ## good                                                                              \n",
    "adata.obs['ct_cov'] = adata.obs.ct_cov.astype('category')\n",
    "sc.settings.set_figure_params(dpi=100, dpi_save=300, format='png', frameon=False, transparent=True, fontsize=16)\n",
    "plt.rcParams[\"image.aspect\"] = \"equal\"\n",
    "MasterORDER = ['Classical Monocytes','Nonclassical Monocytes', 'Classical Dendritic Cells', 'Plasmacytoid Dendritic Cells', 'Naive T Cells', 'Effector Memory T Cells', 'Cytotoxic T Cells', 'Proliferating T Cells', 'Natural Killer Cells', 'B Cells', 'Progenitor Cells', 'Megakaryocytes', 'RBCs']\n",
    "#colorrs = [\"#E58606\",\"#5D69B1\",\"#52BCA3\",\"#99C945\",\"#CC61B0\",\"#24796C\",\"#DAA51B\",\"#2F8AC4\",\"#764E9F\",\"#ED645A\",\"#CC3A8E\",'#BC23FF', '#D790FF']\n",
    "colorrs = [\"#4E79A7\",\"#A0CBE8\",\"#F28E2B\",\"#FFBE7D\",\"#8CD17D\",\"#B6992D\",\"#499894\",\"#E15759\",\"#FF9D9A\",\"#79706E\",\"#D37295\",\"#FABFD2\", '#000000',\"#B07AA1\",\"#D4A6C8\",\"#9D7660\",\n",
    "                 \"#E58606\", \"#5D69B1\", \"#24796C\", '#DAA51B', '#000000', '#99C945', '#ED645A']\n",
    "\n",
    "adata.obs['ct_cov'] = adata.obs['ct_cov'].cat.reorder_categories(MasterORDER)\n",
    "adata.uns['ct_cov_colors'] = colorrs\n",
    "celltype_umap = sc.pl.umap(adata, color='ct_cov', show=True, size=3, edgecolor=\"none\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example antibodies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['CD45RA|PTPRC|j95-28|pAbO', 'CD45RO|PTPRC|j95-19|pAbO', 'CD4|CD4|j95-14|pAbO', 'CD8|CD8A|j95-25|pAbO'], size=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Label cell populations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Proportions and Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Make proportion plots\n",
    "adata_obs_small = adata.obs\n",
    "ind_count = adata_obs_small.groupby(['ind_cov_disease_cov','ct_cov','disease_cov','pop_cov', 'ind_cov'])['ct_cov'].count()\n",
    "ind_count_sums = ind_count.groupby(level=[0]).sum()\n",
    "ind_count_sums = ind_count_sums.reset_index(name=\"counts\")\n",
    "ind_perc = ind_count/ind_count.groupby(level=[0]).transform(sum)*100\n",
    "ind_perc = ind_perc.reset_index(name=\"ct_perc\")\n",
    "# Add weights to WLS\n",
    "ind_perc['counts'] = ind_count.values.tolist()\n",
    "\n",
    "ind_perc['ind_count_sum'] = list(np.zeros(len(ind_count.values.tolist()),dtype=int))\n",
    "# Add total sums per individual to structure\n",
    "for ii in range(len(ind_count_sums)):\n",
    "    ind_perc['ind_count_sum'][ind_perc.ind_cov_disease_cov==ind_count_sums.ind_cov_disease_cov[ii]] = ind_count_sums.counts[ii]\n",
    "\n",
    "\n",
    "ind_perc.ind_cov_disease_cov = ind_perc.ind_cov_disease_cov.apply(lambda x: x.split('_')[0])\n",
    "ind_count_sums.ind_cov_disease_cov = ind_count_sums.ind_cov_disease_cov.apply(lambda x: x.split('_')[0])\n",
    "#ind_perc = ind_perc.set_index('ct_cov').join(cg_cov.set_index('ct_cov'))\n",
    "ind_perc.ind_cov_disease_cov = ind_perc.ind_cov_disease_cov.astype(\"str\")\n",
    "ind_perc.reset_index(inplace=True)\n",
    "ind_perc.ct_cov = ind_perc.ct_cov.astype('category')\n",
    "ind_perc.ct_cov = ind_perc.ct_cov.cat.reorder_categories(adata.obs.ct_cov.cat.categories.values)\n",
    "\n",
    "perc_plot = sns.catplot(x='disease_cov', y='ct_perc', order=[ \"Healthy\", \"Managed\", \"Treated\", \"Flare\"], hue='ct_cov', data=ind_perc, kind='violin', col_order=MasterORDER, col='ct_cov', col_wrap=3, cut=0, dodge=False, aspect=1, sharex=False, sharey=False, palette=colorrs)\n",
    "\n",
    "for ct_i in list(range(len(MasterORDER))):\n",
    "    ct = MasterORDER[ct_i]\n",
    "    sns.swarmplot(x=\"disease_cov\", y=\"ct_perc\", data=ind_perc[ind_perc.ct_cov == ct], order=[ \"Healthy\", \"Managed\", \"Treated\", \"Flare\"], color=\"0\", size=4, ax=perc_plot.axes[ct_i])\n",
    "    try:\n",
    "        sns.pointplot(x=\"disease_cov\", y=\"ct_perc\", hue=\"ind_cov\", data=ind_perc[ind_perc.ct_cov == ct], order=[ \"Healthy\", \"Managed\", \"Treated\", \"Flare\"], color=\"0\", scale=0.3, ax=perc_plot.axes[ct_i])\n",
    "    except:\n",
    "        continue\n",
    "    perc_plot.axes[ct_i].get_yaxis().label.set_visible(False)\n",
    "    perc_plot.axes[ct_i].get_xaxis().label.set_visible(False)\n",
    "    perc_plot.axes[ct_i].get_legend().remove()\n",
    "    perc_plot.set_xticklabels(rotation=90)\n",
    "    perc_plot.fig.subplots_adjust(wspace=2, hspace = 1)\n",
    "    \n",
    "perc_plot.savefig(figdir+\"/violin.ct_cov_figure6b.png\")\n",
    "perc_plot.savefig(figdir+\"/violin.ct_cov_figure6b.pdf\")\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Healthy vs. Managed WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Healthy', 'Managed'])]\n",
    "HEALTHYMANAGED_effect = {}\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease==\"Managed\"] = 0\n",
    "    disease.values[disease==\"Healthy\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    HEALTHYMANAGED_effect[ct] = effect\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Healthy vs. Treated WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Healthy', 'Treated'])]\n",
    "HEALTHYTREATED_effect = {}\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease==\"Treated\"] = 0\n",
    "    disease.values[disease==\"Healthy\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    HEALTHYTREATED_effect[ct] = effect\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Healthy vs. Flare WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Healthy', 'Flare'])]\n",
    "HEALTHYFLARE_effect = {}\n",
    "HEALTHYFLARE_FC = {}\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease==\"Flare\"] = 0\n",
    "    disease.values[disease==\"Healthy\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    HEALTHYFLARE_effect[ct] = effect\n",
    "    HEALTHYFLARE_FC[ct] = math.log2(sum(est.params)/est.params[0])\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Treated vs. Flare WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Treated', 'Flare'])]\n",
    "TREATEDFLARE_effect = {}\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease==\"Flare\"] = 0\n",
    "    disease.values[disease==\"Treated\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    TREATEDFLARE_effect[ct] = effect\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Treated vs. Managed WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Treated', 'Managed'])]\n",
    "TREATEDMANAGED_effect = {}\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease==\"Managed\"] = 0\n",
    "    disease.values[disease==\"Treated\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    TREATEDMANAGED_effect[ct] = effect\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Healthy vs. SLE WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Healthy', 'Treated', 'Untreated', 'Managed'])]\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease!=\"Healthy\"] = 0\n",
    "    disease.values[disease==\"Healthy\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Healthy vs. Treated&Flare WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Healthy', 'Treated', 'Untreated'])]\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease!=\"Healthy\"] = 0\n",
    "    disease.values[disease==\"Healthy\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getproportions(celltype, group1, group2):\n",
    "    print('{} proportion for celltype: {}: {}' .format(group1, celltype, ind_perc[ind_perc.disease_cov==group1][ind_perc.ct_cov==celltype].ct_perc.mean()))\n",
    "    print('{} proportion for celltype: {}: {}' .format(group2, celltype, ind_perc[ind_perc.disease_cov==group2][ind_perc.ct_cov==celltype].ct_perc.mean()))\n",
    "\n",
    "    \n",
    "getproportions(celltype='CD8+ T Cells',group1='Healthy', group2='Flare')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Feature plots for Untreated Treated pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "pairs = ['FLARE004', 'FLARE009', 'FLARE011', 'FLARE013', 'FLARE016', 'FLARE003', 'FLARE002', 'FLARE008']\n",
    "colors = [\"#000000\", \"#FF0000\"]\n",
    "for pair in pairs:\n",
    "    bdata = adata[adata.obs['ind_cov']==pair]\n",
    "    print(pair)\n",
    "    sc.pl.umap(bdata, size=10, show=True, color='disease_cov', save=str('TREATED_UNTREATED_PAIR'+ pair + '.disease.png'), palette=colors)\n",
    "    plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Map colors onto cell type labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colormapping = dict()\n",
    "unique_ct = adata.obs.ct_cov.cat.categories.tolist()\n",
    "for ii in range(len(unique_ct)):\n",
    "    colormapping[unique_ct[ii]] = adata.uns['ct_cov_colors'][ii]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Feature plots for marker genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Atypical memory B cells\n",
    "sc.pl.umap(adata, color=['CR2', 'CD27', 'FCRL5', 'ITGAX', 'NKG7'], show=False, save='.AtypicalBgenes.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## platelet effect\n",
    "sc.pl.umap(adata, color=[\"PF4\", \"SDPR\", \"GNG11\", \"PPBP\"], size=3, show=False, save='.platelet.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## cycling\n",
    "sc.pl.umap(adata, color=[\"KIAA0101\",\"STMN1\",\"TK1\",\"MKI67\"], size=3, show=False, save='.cycling.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## MACROPHAGES\n",
    "sc.pl.umap(adata, color=['CD163', 'HLA-DRB1', 'C1QA', 'IFITM3'], size=3, show=False, save=\".MACRO.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=[\"CD7\", 'PDCD1', 'CD28'],size=3, show=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## CD4, CD8 and NK axis\n",
    "sc.pl.umap(adata, color=[\"CD3D\",\"CD8A\",\"CD4\",\"NCAM1\", \"FCGR3A\", 'NKG7', 'IFNG', 'GZMB', 'PRF1'],size=3, show=True, save=\".TNK.png\")\n",
    "# Two NK sub populations: CD56high/lowCD16high/low, CD56lowCD16high\n",
    "# https://www.frontiersin.org/files/Articles/162361/fimmu-06-00567-HTML/image_m/fimmu-06-00567-g001.jpg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sc.pl.violin(adata, keys= ['GZMB', 'PRF1', 'SH2D1A', 'ICAM1','CD86', 'RAB27A', 'UNC13D', 'FHL3', 'FHL2'], groupby='disease_cov', show=False, use_raw=False, save='.T8emcyto1expression')\n",
    "#sc.pl.stacked_violin(adata, keys= ['GZMB', 'PRF1', 'SH2D1A', 'ICAM1','CD86', 'RAB27A', 'UNC13D', 'FHL3', 'FHL2'], groupby='disease_cov', show=False, use_raw=False, save='.T8emcyto1expression')\n",
    "sc.pl.umap(adata, color=['GZMB', 'PRF1', 'SH2D1A', 'ICAM1','CD86', 'RAB27A', 'UNC13D', 'FHL3', 'FHL2'],size=3, show=False, use_raw=True, save=\".gzmbgenes.png\")\n",
    "sc.pl.violin(adata, keys= ['GZMB', 'PRF1', 'SH2D1A', 'ICAM1','CD86', 'RAB27A', 'UNC13D', 'FHL3', 'FHL2'], groupby='disease_cov', show=False, use_raw=True, save='.T8emcyto1expression_raw')\n",
    "#sc.pl.stacked_violin(adata, var_names= ['GZMB', 'PRF1', 'SH2D1A', 'ICAM1','CD86', 'RAB27A', 'UNC13D', 'FHL3', 'FHL2'], groupby='disease_cov', show=False, use_raw=True, save='.T8emcyto1expression_raw')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Memory vs. helper vs. naive\n",
    "sc.pl.umap(adata, color=[\"CCR7\",\"IL7R\",\"S100A4\",\"CD58\", 'FAS', 'IL2RA'],show=False, size=3, save=\".TMemThTNaive.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## treg\n",
    "sc.pl.umap(adata, color=[\"FOXP3\",\"TNFRSF4\",\"ENTPD1\",\"CCR10\"],size=3, show=False, save=\".TREG.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## y chromosome gender effect\n",
    "sc.pl.umap(adata, color=[\"DDX3Y\", \"RPS4Y1\", \"FHIT\",\"TRAT1\"], size =3, show=False, save=\".Y.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## b cells, plasmablasts and pdcs\n",
    "sc.pl.umap(adata, color=['BTLA', 'P2RY8',\"MZB1\", \"CD19\", \"CD79A\", 'MS4A1', 'FCRL5', 'IL6', 'CR2'],size=3, show=False, save=\".B.png\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['ITGAM', 'CD33', 'FUT4', 'HLA-DRA', 'CD14', 'CD163'],size=3, show=False, save=\".B.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cytotoxic signature: Isolate cytotoxic CD3+ cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bdata = adata[adata.obs['ct_cov'].isin(['Cytotoxic T Cells'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(bdata, color=['CCR7', 'GZMB', 'PRF1', 'CD4', 'CD8A', 'leiden'], size=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(bdata, color=['CD4|CD4|j95-14|pAbO', 'CD8|CD8A|j95-25|pAbO'], size=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.asarray(bdata.obs['CD4|CD4|j95-14|pAbO'].tolist())\n",
    "Y = np.asarray(bdata.obs['CD8|CD8A|j95-25|pAbO'].tolist())\n",
    "f, ax = plt.subplots(figsize=(7, 7))\n",
    "ax = sns.scatterplot(X,Y)\n",
    "plt.xlabel('CD4')\n",
    "plt.ylabel('CD8A')\n",
    "joint_kws=dict(gridsize=50)\n",
    "sns.jointplot(X, Y, kind=\"hex\", color=\"b\", joint_kws= joint_kws);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.sum((X>np.nanmax(X)*.10) & (Y<np.nanmax(Y)*.10))) # CD4 only\n",
    "print(np.sum((X<np.nanmax(X)*.10) & (Y>np.nanmax(Y)*.10))) # CD8 only\n",
    "print(np.sum((X<np.nanmax(X)*.10) & (Y<np.nanmax(Y)*.10))) # Niether\n",
    "print(np.sum((X>np.nanmax(X)*.10) & (Y>np.nanmax(Y)*.10))) # CD4+CD8+"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reference populations were found by gating CD4 and CD8 antibody expression and then high resolution clustering was performed. A spearman r was computed for each leiden cluster to assign population identity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdata = sc.read('Lupus_study_cytoTCell.h5ad', cache=True)\n",
    "print(cdata)\n",
    "print('Number of CD4+ cytotoxic cells: {}'.format(np.sum(cdata.obs['ct_cov'].values=='CD4+ Cytotoxic T Cells')))\n",
    "print('Number of CD8+ cytotoxic cells: {}'.format(np.sum(cdata.obs['ct_cov'].values=='CD8+ Cytotoxic T Cells')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(cdata, color=['CCR7', 'GZMB', 'PRF1', 'CD4', 'CD8A', 'ct_cov'], size=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_gene_expression(cdata, celltypes, gene):\n",
    "    indlist  = cdata.obs['ind_cov_disease_cov'].unique().tolist()\n",
    "    ind_id = []; dz_id = []; ct_id = []; values = []; count = []\n",
    "    for celltype in celltypes:\n",
    "        for ii in range(len(indlist)):\n",
    "            values.append(cdata.raw.X[(cdata.obs['ct_cov'].values==celltype) & (cdata.obs['ind_cov_disease_cov'].values==indlist[ii]), cdata.raw.var_names==gene].mean())\n",
    "            count.append(cdata.raw.X[cdata.obs['ind_cov_disease_cov'].values==indlist[ii], cdata.raw.var_names==gene].shape[1])\n",
    "            ind_id.append(indlist[ii])\n",
    "            dz_id.append(cdata.obs['disease_cov'][cdata.obs['ind_cov_disease_cov']==indlist[ii]].values.unique()[0])\n",
    "            ct_id.append(celltype)\n",
    "\n",
    "    genexpression = pd.DataFrame(data={'Unique_ID': ind_id, 'disease_cov': dz_id, gene: values, 'Cell_Count':count, 'ct_cov':ct_id})    \n",
    "    genexpression.ct_cov = genexpression.ct_cov.astype('category')\n",
    "    perc_plot = sns.catplot(x='disease_cov', y=gene, order=[ \"Healthy\", \"Managed\", \"Treated\", \"Flare\"], hue='ct_cov', data=genexpression, kind='violin', col='ct_cov', col_wrap=2, cut=0, dodge=False, aspect=1, sharex=False, sharey=True)\n",
    "    for ct_i in range(len(celltypes)):\n",
    "        sns.swarmplot(x=\"disease_cov\", y=gene, data=genexpression[genexpression.ct_cov==celltypes[ct_i]], order=[ \"Healthy\", \"Managed\", \"Treated\", \"Flare\"], color=\"0\", size=6, ax=perc_plot.axes[ct_i])\n",
    "        try:\n",
    "            sns.pointplot(x=\"disease_cov\", y=gene, hue=\"Unique_ID\", data=genexpression[genexpression.Cell_Type==celltypes[ct_i]], order=[ \"Healthy\", \"Managed\", \"Treated\", \"Flare\"], color=\"0\", scale=0.3, ax=perc_plot.axes[ct_i])\n",
    "        except:\n",
    "            continue\n",
    "        perc_plot.axes[ct_i].get_xaxis().label.set_visible(False)\n",
    "        perc_plot.axes[ct_i].get_legend().remove()\n",
    "        perc_plot.set_xticklabels(rotation=90)\n",
    "        perc_plot.fig.subplots_adjust(wspace=0.5, hspace = 1)\n",
    "\n",
    "    return genexpression\n",
    "\n",
    "def get_stats(dataframe, groups, gene):\n",
    "    print('{} groups weighted by number of cells'.format(groups))\n",
    "    all_out = pd.DataFrame();\n",
    "    if 'SLE' in groups:\n",
    "        ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Treated', 'Managed', 'Healthy', 'Flare'])]\n",
    "    else:\n",
    "        ind_perc0 = ind_perc[ind_perc.disease_cov.isin(groups)]\n",
    "    Effect = {}\n",
    "    for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "        ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "        ct_perc = ind_perc0[gene][ind_perc0.ct_cov==ct]\n",
    "        weights = ind_perc0.Cell_Count[ind_perc0.ct_cov==ct]\n",
    "        disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "        disease = disease.astype(\"str\")\n",
    "        if 'SLE' in groups:\n",
    "            disease.values[disease=='Healthy'] = 0\n",
    "            disease.values[disease!='Healthy'] = 1\n",
    "        else:\n",
    "            disease.values[disease==groups[0]] = 0\n",
    "            disease.values[disease==groups[1]] = 1\n",
    "        disease = sm.add_constant(disease)\n",
    "        est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "        est=est.fit()\n",
    "        effect=est.params[1]\n",
    "        pval=est.pvalues[1]\n",
    "        Effect[ct] = effect\n",
    "        all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "        all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "    display(all_out)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "celltypes = ['CD4+ Cytotoxic T Cells', 'CD8+ Cytotoxic T Cells']\n",
    "genes     = ['IFITM1', 'CD69']\n",
    "#genes =['MT-ND1', 'MT-ND2', 'MT-CO1', 'MT-CO2', 'MT-ATP8', 'MT-ATP6', 'MT-CO3', 'MT-ND3', 'MT-ND4L', 'MT-ND4', 'MT-ND5', 'MT-ND6', 'MT-CYB']\n",
    "for ii in range(len(genes)):\n",
    "    ind_perc = plot_gene_expression(cdata, celltypes, genes[ii])\n",
    "    get_stats(ind_perc, ['Healthy', 'Managed'], genes[ii])\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Look at genes of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "celltypes = ['CD4+ Cytotoxic T Cells', 'CD8+ Cytotoxic T Cells']\n",
    "genes     = ['CD4', 'CD40LG', 'ZBTB7B', 'CD3G', 'TNFRSF8', 'TNFRSF4', 'TNFRSF9', 'ICOS', 'TBX21', 'CD8A', 'CD8B','CD27', 'PRF1', 'GZMH', 'GZMK', 'GZMM', 'RUNX3', 'PRDM1', 'CD28', 'KLRK1', 'CTLA4', 'PDCD1', 'BCL2', 'BCL2L11','BAX','EOMES']\n",
    "for ii in range(len(genes)):\n",
    "    plot_gene_expression(cdata, celltypes, genes[ii])\n",
    "    get_stats(ind_perc, ['Healthy', 'Managed'], genes[ii])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "genes = ['CD4', 'CD40LG', 'ZBTB7B', 'CD3G', 'TNFRSF8', 'TNFRSF4', 'TNFRSF9', 'ICOS', 'TBX21', 'CD8A', 'CD8B','CD27', 'PRF1', 'GZMH', 'GZMK', 'GZMM', 'RUNX3', 'PRDM1', 'CD28', 'KLRK1', 'CTLA4', 'PDCD1', 'BCL2', 'BCL2L11','BAX','EOMES']\n",
    "sc.pl.stacked_violin(cdata, groupby='ct_cov', use_raw=True,var_names=genes, show=True, save='Cytogene.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differentially expressed genes between CD4+ and CD8+ cytotoxic populations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.rank_genes_groups(cdata, groupby='ct_cov')\n",
    "sc.pl.rank_genes_groups(cdata, n_genes=50, fontsize=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differentially expressed genes between CD4+ and CD8+ cytotoxic populations (only SLE patients)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddata = cdata[cdata.obs['disease_cov'].isin(['Managed', 'Treated', 'Flare'])]\n",
    "sc.tl.rank_genes_groups(ddata, groupby='ct_cov',n_genes=25)\n",
    "sc.pl.rank_genes_groups(ddata, n_genes=25, fontsize=8)\n",
    "df = pd.DataFrame(ddata.uns['rank_genes_groups']['names'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "genes = df['CD4+ Cytotoxic T Cells'].tolist() + df['CD8+ Cytotoxic T Cells'].tolist()\n",
    "for ii in range(len(genes)):\n",
    "    plot_gene_expression(cdata, celltypes, genes[ii])\n",
    "    #get_stats(ind_perc, ['Healthy', 'Managed'], genes[ii])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differentially expressed genes between SLE and Healthy for CD4+ cytotoxic cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddata = cdata[cdata.obs['ct_cov']=='CD4+ Cytotoxic T Cells']\n",
    "sc.tl.rank_genes_groups(ddata, groupby='disease_cov',groups=['Healthy', 'Managed', 'Treated', 'Flare'],reference='Healthy', n_genes=25)\n",
    "sc.pl.rank_genes_groups(ddata, n_genes=25)\n",
    "sc.pl.rank_genes_groups_stacked_violin(ddata, groupby='disease_cov', n_genes=25,dendrogram=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame(ddata.uns['rank_genes_groups']['names'])\n",
    "genes = df['Managed'].tolist() + df['Treated'].tolist()+df['Flare'].tolist()\n",
    "for ii in range(len(genes)):\n",
    "    plot_gene_expression(cdata, celltypes, genes[ii])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differentially expressed genes between SLE and Healthy for CD8+ cytotoxic cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddata = cdata[cdata.obs['ct_cov']=='CD8+ Cytotoxic T Cells']\n",
    "sc.tl.rank_genes_groups(ddata, groupby='disease_cov',groups=['Healthy', 'Managed', 'Treated', 'Flare'],reference='Healthy', n_genes=50)\n",
    "sc.pl.rank_genes_groups(ddata, n_genes=25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gene rank scores for CD8+ cytotoxic vs. CD4+ cytotoxic and Healthy versus SLE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scores for differentially expressed genes in disease status\n",
    "sc.tl.rank_genes_groups(cdata, groupby='disease_cov', n_genes=len(mastergenes), groups=['Healthy', 'Managed', 'Treated', 'Flare'], reference='Healthy')\n",
    "genescores = cdata.uns['rank_genes_groups']['scores']\n",
    "genenames = cdata.uns['rank_genes_groups']['names']\n",
    "df = pd.DataFrame(genenames)\n",
    "dfscore = pd.DataFrame(genescores)\n",
    "dfscore = dfscore.rename(columns={'Managed':'Managed_score', 'Treated': 'Treated_score', 'Flare': 'Flare_score'})\n",
    "\n",
    "dfmanaged = pd.concat([df['Managed'], dfscore['Managed_score']], axis=1)\n",
    "dfmanaged = dfmanaged.set_index('Managed')\n",
    "dfmanaged.index.rename(name='Gene', inplace=True)\n",
    "\n",
    "dftreated = pd.concat([df['Treated'], dfscore['Treated_score']], axis=1)\n",
    "dftreated = dftreated.set_index('Treated')\n",
    "dftreated.index.rename(name='Gene', inplace=True)\n",
    "\n",
    "dfflare = pd.concat([df['Flare'], dfscore['Flare_score']], axis=1)\n",
    "dfflare = dfflare.set_index('Flare')\n",
    "dfflare.index.rename(name='Gene', inplace=True)\n",
    "\n",
    "# Scores for differentially expressed genes for CD8+ and CD4+ Cytotoxic T cells\n",
    "sc.tl.rank_genes_groups(cdata, groupby='ct_cov', n_genes=len(mastergenes),groups=['CD4+ Cytotoxic T Cells', 'CD8+ Cytotoxic T Cells'], reference='CD4+ Cytotoxic T Cells')\n",
    "genescores = cdata.uns['rank_genes_groups']['scores']\n",
    "genenames = cdata.uns['rank_genes_groups']['names']\n",
    "df2 = pd.DataFrame(genenames)\n",
    "df2score = pd.DataFrame(genescores)\n",
    "df2score = df2score.rename(columns={'CD8+ Cytotoxic T Cells':'CD4+<->CD8+_score'})\n",
    "df2 = pd.concat([df2, df2score], axis=1)\n",
    "df2 = df2.set_index('CD8+ Cytotoxic T Cells')\n",
    "df2.index.rename(name='Gene', inplace=True)\n",
    "\n",
    "\n",
    "dfmanaged = pd.concat([dfmanaged, df2], axis=1)\n",
    "dfmanaged = dfmanaged.sort_values(by=['Managed_score', 'CD4+<->CD8+_score'])\n",
    "dftreated = pd.concat([dftreated, df2], axis=1)\n",
    "dftreated = dftreated.sort_values(by=['Treated_score', 'CD4+<->CD8+_score'])\n",
    "dfflare = pd.concat([dfflare, df2], axis=1)\n",
    "dfflare = dfflare.sort_values(by=['Flare_score', 'CD4+<->CD8+_score'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.scatterplot(x='CD4+<->CD8+_score', y='Managed_score', data=dfmanaged, label='Managed')\n",
    "sns.scatterplot(x='CD4+<->CD8+_score', y='Treated_score', data=dftreated, label='Treated')\n",
    "sns.scatterplot(x='CD4+<->CD8+_score', y='Flare_score', data=dfflare, label='Flare')\n",
    "plt.ylabel('Healthy vs. SLE')\n",
    "plt.xlabel('CD4+ vs CD8+')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rename cytotoxic cells in main dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CD4 = cdata[cdata.obs['ct_cov']=='CD4+ Cytotoxic T Cells'].obs_names.tolist()\n",
    "CD8 = cdata[cdata.obs['ct_cov']=='CD8+ Cytotoxic T Cells'].obs_names.tolist()\n",
    "adata.obs['ct_cov'] = adata.obs['ct_cov'].astype('object')\n",
    "adata.obs['ct_cov'][adata.obs_names.isin(CD4)] = 'CD4+ Cytotoxic T Cells'\n",
    "adata.obs['ct_cov'][adata.obs_names.isin(CD8)] = 'CD8+ Cytotoxic T Cells'\n",
    "adata.obs['ct_cov'] = adata.obs.ct_cov.astype('category')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MasterORDER = ['Classical Monocytes','Nonclassical Monocytes', 'Classical Dendritic Cells', 'Plasmacytoid Dendritic Cells', 'Naive T Cells', 'Effector Memory T Cells', 'CD8+ Cytotoxic T Cells', 'CD4+ Cytotoxic T Cells', 'Proliferating T Cells', 'Natural Killer Cells', 'B Cells', 'Progenitor Cells', 'Megakaryocytes', 'RBCs']\n",
    "#colorrs = [\"#E58606\",\"#5D69B1\",\"#52BCA3\",\"#99C945\",\"#CC61B0\",\"#24796C\",\"#DAA51B\",\"#2F8AC4\",\"#764E9F\",\"#ED645A\",\"#CC3A8E\",'#BC23FF', '#D790FF']\n",
    "colorrs = [\"#4E79A7\",\"#A0CBE8\",\"#F28E2B\",\"#FFBE7D\",\"#8CD17D\",\"#B6992D\",\"#499894\",\"#E15759\",\"#FF9D9A\",\"#79706E\",\"#D37295\",\"#FABFD2\",\"#B07AA1\",\"#D4A6C8\",\"#9D7660\",\n",
    "                 \"#E58606\", \"#5D69B1\", \"#24796C\", '#DAA51B', '#000000', '#99C945', '#ED645A']\n",
    "\n",
    "adata.obs['ct_cov'] = adata.obs['ct_cov'].cat.reorder_categories(MasterORDER)\n",
    "adata.uns['ct_cov_colors'] = colorrs\n",
    "celltype_umap = sc.pl.umap(adata, color='ct_cov', show=True, size=3, edgecolor=\"none\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Make proportion plots\n",
    "adata_obs_small = adata.obs\n",
    "ind_count = adata_obs_small.groupby(['ind_cov_disease_cov','ct_cov','disease_cov','pop_cov', 'ind_cov'])['ct_cov'].count()\n",
    "ind_count_sums = ind_count.groupby(level=[0]).sum()\n",
    "ind_count_sums = ind_count_sums.reset_index(name=\"counts\")\n",
    "ind_perc = ind_count/ind_count.groupby(level=[0]).transform(sum)*100\n",
    "ind_perc = ind_perc.reset_index(name=\"ct_perc\")\n",
    "# Add weights to WLS\n",
    "ind_perc['counts'] = ind_count.values.tolist()\n",
    "\n",
    "ind_perc['ind_count_sum'] = list(np.zeros(len(ind_count.values.tolist()),dtype=int))\n",
    "# Add total sums per individual to structure\n",
    "for ii in range(len(ind_count_sums)):\n",
    "    ind_perc['ind_count_sum'][ind_perc.ind_cov_disease_cov==ind_count_sums.ind_cov_disease_cov[ii]] = ind_count_sums.counts[ii]\n",
    "\n",
    "# Filter samples that have < 100 cells\n",
    "ind_perc = ind_perc[ind_perc.ind_cov_disease_cov.isin(ind_count_sums.ind_cov_disease_cov[ind_count_sums.counts > 100])]\n",
    "\n",
    "\n",
    "ind_perc.ind_cov_disease_cov = ind_perc.ind_cov_disease_cov.apply(lambda x: x.split('_')[0])\n",
    "ind_count_sums.ind_cov_disease_cov = ind_count_sums.ind_cov_disease_cov.apply(lambda x: x.split('_')[0])\n",
    "#ind_perc = ind_perc.set_index('ct_cov').join(cg_cov.set_index('ct_cov'))\n",
    "ind_perc.ind_cov_disease_cov = ind_perc.ind_cov_disease_cov.astype(\"str\")\n",
    "ind_perc.reset_index(inplace=True)\n",
    "ind_perc.ct_cov = ind_perc.ct_cov.astype('category')\n",
    "ind_perc.ct_cov = ind_perc.ct_cov.cat.reorder_categories(adata.obs.ct_cov.cat.categories.values)\n",
    "\n",
    "perc_plot = sns.catplot(x='disease_cov', y='ct_perc', order=[ \"Healthy\", \"Managed\", \"Treated\", \"Flare\"], hue='ct_cov', data=ind_perc, kind='violin', col_order=MasterORDER, col='ct_cov', col_wrap=3, cut=0, dodge=False, aspect=1, sharex=False, sharey=False, palette=colorrs)\n",
    "\n",
    "for ct_i in list(range(len(MasterORDER))):\n",
    "    ct = MasterORDER[ct_i]\n",
    "    sns.swarmplot(x=\"disease_cov\", y=\"ct_perc\", data=ind_perc[ind_perc.ct_cov == ct], order=[ \"Healthy\", \"Managed\", \"Treated\", \"Flare\"], color=\"0\", size=4, ax=perc_plot.axes[ct_i])\n",
    "    try:\n",
    "        sns.pointplot(x=\"disease_cov\", y=\"ct_perc\", hue=\"ind_cov\", data=ind_perc[ind_perc.ct_cov == ct], order=[ \"Healthy\", \"Managed\", \"Treated\", \"Flare\"], color=\"0\", scale=0.3, ax=perc_plot.axes[ct_i])\n",
    "    except:\n",
    "        continue\n",
    "    perc_plot.axes[ct_i].get_yaxis().label.set_visible(False)\n",
    "    perc_plot.axes[ct_i].get_xaxis().label.set_visible(False)\n",
    "    perc_plot.axes[ct_i].get_legend().remove()\n",
    "    perc_plot.set_xticklabels(rotation=90)\n",
    "    perc_plot.fig.subplots_adjust(wspace=2, hspace = 1)\n",
    "    \n",
    "perc_plot.savefig(figdir+\"/violin.ct_cov_figure6b.png\")\n",
    "perc_plot.savefig(figdir+\"/violin.ct_cov_figure6b.pdf\")\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Healthy vs. Managed WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Healthy', 'Managed'])]\n",
    "HEALTHYMANAGED_effect = {}\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease==\"Managed\"] = 0\n",
    "    disease.values[disease==\"Healthy\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    HEALTHYMANAGED_effect[ct] = effect\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Healthy vs. Treated WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Healthy', 'Treated'])]\n",
    "HEALTHYTREATED_effect = {}\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease==\"Treated\"] = 0\n",
    "    disease.values[disease==\"Healthy\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    HEALTHYTREATED_effect[ct] = effect\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Healthy vs. Flare WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Healthy', 'Flare'])]\n",
    "HEALTHYFLARE_effect = {}\n",
    "HEALTHYFLARE_FC = {}\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease==\"Flare\"] = 0\n",
    "    disease.values[disease==\"Healthy\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    HEALTHYFLARE_effect[ct] = effect\n",
    "    HEALTHYFLARE_FC[ct] = math.log2(sum(est.params)/est.params[0])\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Treated vs. Flare WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Treated', 'Flare'])]\n",
    "TREATEDFLARE_effect = {}\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease==\"Flare\"] = 0\n",
    "    disease.values[disease==\"Treated\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    TREATEDFLARE_effect[ct] = effect\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Treated vs. Managed WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Treated', 'Managed'])]\n",
    "TREATEDMANAGED_effect = {}\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease==\"Managed\"] = 0\n",
    "    disease.values[disease==\"Treated\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    TREATEDMANAGED_effect[ct] = effect\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Healthy vs. SLE WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Healthy', 'Treated', 'Untreated', 'Managed'])]\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease!=\"Healthy\"] = 0\n",
    "    disease.values[disease==\"Healthy\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n",
    "\n",
    "# Statistical tests...\n",
    "print('as a proportion of total PBMC: Healthy vs. Treated&Flare WEIGHTED BY TOTAL PBMC COUNTS')\n",
    "all_out = pd.DataFrame();\n",
    "ind_perc0 = ind_perc[ind_perc.disease_cov.isin(['Healthy', 'Treated', 'Untreated'])]\n",
    "for ct_i in list(range(len(ind_perc0.ct_cov.cat.categories))):\n",
    "    ct = ind_perc0.ct_cov.cat.categories[ct_i]\n",
    "    ct_perc = ind_perc0.ct_perc[ind_perc0.ct_cov==ct]\n",
    "    weights = ind_perc0.ind_count_sum[ind_perc0.ct_cov==ct]\n",
    "    disease = ind_perc0.disease_cov[ind_perc0.ct_cov==ct]\n",
    "    disease = disease.astype(\"str\")\n",
    "    disease.values[disease!=\"Healthy\"] = 0\n",
    "    disease.values[disease==\"Healthy\"] = 1\n",
    "    disease = sm.add_constant(disease)\n",
    "    est=sm.WLS(ct_perc.astype(float), disease.astype(float), weights=weights)\n",
    "    #est=sm.OLS(ct_perc.astype(float), disease.astype(float))\n",
    "    est=est.fit()\n",
    "    effect=est.params[1]\n",
    "    pval=est.pvalues[1]\n",
    "    all_out = all_out.append(({\"Cell\":str(ct),\"Beta\":str(effect),\"Pval\":str(pval)}), ignore_index=True)\n",
    "    all_out['Pval'] = all_out['Pval'].astype(float)\n",
    "display(all_out)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gene Ontology analysis for each cell type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from Bio import Entrez\n",
    "from goatools.base import download_go_basic_obo\n",
    "from goatools.base import download_ncbi_associations\n",
    "from goatools.obo_parser import GODag\n",
    "from goatools.anno.genetogo_reader import Gene2GoReader\n",
    "from goatools.test_data.genes_NCBI_9606_ProteinCoding import GENEID2NT as GeneID2nt_hsa\n",
    "from goatools.goea.go_enrichment_ns import GOEnrichmentStudyNS\n",
    "from goatools.godag_plot import plot_gos, plot_results, plot_goid2goobj\n",
    "# Store gene summaries\n",
    "summary_dict = {}\n",
    "def getgeneids(gene_list):\n",
    "    Entrez.email = \"richard.perez@ucsf.edu\"\n",
    "    # Store symbol and id\n",
    "    geneid2symbol = {}\n",
    "    for genesymbol in gene_list:\n",
    "        try:\n",
    "            handle = Entrez.esearch(db=\"gene\", term=str(genesymbol + \"[GENE] AND Homo\"))\n",
    "            record = Entrez.read(handle)\n",
    "            geneid2symbol[int(record['IdList'][0])] = genesymbol\n",
    "            handle.close()\n",
    "        except:\n",
    "            continue\n",
    "    return geneid2symbol\n",
    "\n",
    "def id2symbol(GENES, celltypes, cdata):\n",
    "    for geneID in GENES:\n",
    "        Entrez.email = \"richard.perez@ucsf.edu\"\n",
    "        id_list= [str(geneID)]\n",
    "        request = Entrez.epost(\"gene\",id=\",\".join(id_list))\n",
    "        result = Entrez.read(request)\n",
    "        webEnv = result[\"WebEnv\"]\n",
    "        queryKey = result[\"QueryKey\"]\n",
    "        data = Entrez.esummary(db=\"gene\", webenv=webEnv, query_key=queryKey)\n",
    "        annotations = Entrez.read(data)\n",
    "        symbol  = annotations['DocumentSummarySet']['DocumentSummary'][0]['Name']\n",
    "        summary = annotations['DocumentSummarySet']['DocumentSummary'][0]['Summary']\n",
    "        print(symbol)\n",
    "        print(summary)\n",
    "        ind_perc = plot_gene_expression(cdata, celltypes, symbol)\n",
    "    return symbol, summary\n",
    "\n",
    "def id2symbolfast(gene_id):\n",
    "    Entrez.email = \"richard.perez@ucsf.edu\"\n",
    "    id_list= [str(gene_id)]\n",
    "    request = Entrez.epost(\"gene\",id=\",\".join(id_list))\n",
    "    result = Entrez.read(request)\n",
    "    webEnv = result[\"WebEnv\"]\n",
    "    queryKey = result[\"QueryKey\"]\n",
    "    data = Entrez.esummary(db=\"gene\", webenv=webEnv, query_key=queryKey)\n",
    "    annotations = Entrez.read(data)\n",
    "    symbol  = annotations['DocumentSummarySet']['DocumentSummary'][0]['Name']\n",
    "    summary = annotations['DocumentSummarySet']['DocumentSummary'][0]['Summary']\n",
    "    return symbol, summary\n",
    "\n",
    "def replace_id_with_symbol(results, summary_dict):\n",
    "    for ll in range(len(results['study_items'])):\n",
    "        study_item = np.asarray(results['study_items'][ll].split(','), dtype=int)\n",
    "        # Keep all gene summaries\n",
    "        for ii in range(len(study_item)):\n",
    "            symbol, summary = id2symbolfast(study_item[ii])\n",
    "            summary_dict[symbol] = summary\n",
    "        # Replace gene ID with gene symbol\n",
    "        study_item = [id2symbolfast(study_item[ii])[0] for ii in range(len(study_item))]\n",
    "        results['study_items'][ll]= study_item\n",
    "    return results, summary_dict\n",
    "\n",
    "def run_GO(gene_list, cutoff, title):\n",
    "    obo_fname = download_go_basic_obo()\n",
    "    fin_gene2go = download_ncbi_associations()\n",
    "    obodag = GODag(\"go-basic.obo\")\n",
    "\n",
    "    # Read NCBI's gene2go. Store annotations in a list of namedtuples\n",
    "    objanno = Gene2GoReader(fin_gene2go, taxids=[9606])\n",
    "\n",
    "    # Get namespace2association where:\n",
    "    #    namespace is:\n",
    "    #        BP: biological_process               \n",
    "    #        MF: molecular_function\n",
    "    #        CC: cellular_component\n",
    "    #    assocation is a dict:\n",
    "    #        key: NCBI GeneID\n",
    "    #        value: A set of GO IDs associated with that gene\n",
    "    ns2assoc = objanno.get_ns2assc()\n",
    "\n",
    "    for nspc, id2gos in ns2assoc.items():\n",
    "        print(\"{NS} {N:,} annotated human genes\".format(NS=nspc, N=len(id2gos)))\n",
    "\n",
    "    goeaobj = GOEnrichmentStudyNS(\n",
    "            GeneID2nt_hsa.keys(), # List of human protein-coding genes\n",
    "            ns2assoc, # geneid/GO associations\n",
    "            obodag, # Ontologies\n",
    "            propagate_counts = False,\n",
    "            alpha = cutoff, # default significance cut-off\n",
    "            methods = ['fdr_bh']) # defult multipletest correction method\n",
    "\n",
    "    geneid2symbol = getgeneids(gene_list)\n",
    "    # 'p_' means \"pvalue\". 'fdr_bh' is the multiple test method\n",
    "    geneids_study = geneid2symbol.keys()\n",
    "    goea_results_all = goeaobj.run_study(geneids_study)\n",
    "    goea_results_sig = [r for r in goea_results_all if r.p_fdr_bh < cutoff]\n",
    "    '''\n",
    "    This plot contains GOEA results:\n",
    "\n",
    "    GO terms colored by P-value:\n",
    "    pval < 0.005 (light red)\n",
    "    pval < 0.01 (light orange)\n",
    "    pval < 0.05 (yellow)\n",
    "    pval > 0.05 (grey) Study terms that are not statistically significant\n",
    "    GO terms with study gene counts printed. e.g., \"32 genes\"\n",
    "    '''\n",
    "    plot_results(title+\"{NS}.png\", goea_results_sig, id2symbol=geneid2symbol, study_items=20, items_p_line=5)\n",
    "    goeaobj.wr_xlsx(title+\".xlsx\", goea_results_sig, id2symbol=geneid2symbol)\n",
    "    return geneid2symbol, goea_results_sig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pd.set_option('display.max_columns', None)\n",
    "pd.set_option('display.max_colwidth', 100)\n",
    "ct_groups = np.unique(adata.obs['leiden'].tolist())\n",
    "for ii in range(len(ct_groups)):\n",
    "    try:\n",
    "        ddata = adata[adata.obs['leiden']==ct_groups[ii]]\n",
    "        sc.tl.rank_genes_groups(ddata, groupby='SLE status',groups=['Healthy', 'SLE'], reference='Healthy', n_genes=len(ddata.raw.var_names))\n",
    "        sc.pl.rank_genes_groups(ddata, n_genes=25, save='{}_SLEvsHealthy.png'.format(ct_groups[ii]))\n",
    "        genes = pd.DataFrame(ddata.uns['rank_genes_groups']['names'])\n",
    "        scores = pd.DataFrame(ddata.uns['rank_genes_groups']['scores'])\n",
    "        gene_list = list(np.ravel(genes[scores.values>(scores.mean()+(scores.std()*3))[0]].values.tolist()))\n",
    "        cutoff = 0.001\n",
    "        title = str(ct_groups[ii])+'_Lupus_Study'\n",
    "        gene_list = list(np.ravel(genes[scores.values>(scores.mean()+(scores.std()*3))[0]].values.tolist()))\n",
    "        run_GO(gene_list,cutoff,title)\n",
    "        results = pd.read_excel(title+\".xlsx\")\n",
    "        #results, summary_dict = replace_id_with_symbol(results, summary_dict)\n",
    "        print(ct_groups[ii])\n",
    "        print('Number of genes 3 standard deviations above mean: {}'.format(np.sum(scores.values>(scores.mean()+(scores.std()*3))[0])))\n",
    "        display(results)\n",
    "    except:\n",
    "        continue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EBV specific T cells: CD45RA-, HLA-DR+, CD38+ low SELL (CD62L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['CD69', 'HLA-DRA', 'HLA-DRB1', 'CD38', 'CD8A', 'IFNG', 'SELL', 'CD7', 'CD5', 'PAX5'], size=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['CD45RA|PTPRC|j95-28|pAbO', 'CD45RO|PTPRC|j95-19|pAbO', 'CD38|CD38|j95-01|pAbO', 'HLA-DR|CD74|j95-18|pAbO'], size=10)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bdata = adata[adata.obs['ct_cov'].isin(['CD8+ Cytotoxic T Cells'])]\n",
    "X = np.asarray(bdata.obs['CD45RA|PTPRC|j95-28|pAbO'].tolist())\n",
    "Y = np.asarray(bdata.obs['HLA-DR|CD74|j95-18|pAbO'].tolist())\n",
    "f, ax = plt.subplots(figsize=(7, 7))\n",
    "ax = sns.scatterplot(X,Y)\n",
    "plt.xlabel('CD45RA')\n",
    "plt.ylabel('HLA-DR')\n",
    "joint_kws=dict(gridsize=50)\n",
    "sns.jointplot(X, Y, kind=\"hex\", color=\"b\", joint_kws= joint_kws);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# helper T cell population distinct from TFH cells from Nature Medicine 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['IL10', 'PDCD1', 'CD4', 'CXCR3', 'ct_cov'], size=10)"
   ]
  }
 ],
 "metadata": {
  "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.7.4"
  },
  "toc-autonumbering": true,
  "toc-showcode": false,
  "toc-showmarkdowntxt": true,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 2
}