CrossSectional-Flare-Study.V4.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",
    "import scipy.io as io\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",
    "from sklearn.feature_selection import RFE\n",
    "from sklearn.feature_selection import SelectFromModel\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import roc_auc_score\n",
    "from sklearn.metrics import roc_curve\n",
    "from imblearn.over_sampling import SMOTE\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'\n",
    "nonormpath = 'Lupus_study_nonorm.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('Processed dataset: {}'.format(adata))\n",
    "figdir = \"./figures.\"\n",
    "sc.settings.figdir = \"./figures.\"\n",
    "MASTERCOLORS = sc.pl.palettes.default_64\n",
    "MASTERCOLORS.remove(\"#FEFFE6\")\n",
    "MASTERCOLORS.remove(\"#FFFF00\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs['ind_cov'][adata.obs['ind_cov']=='1771'] = '1771_1771'\n",
    "adata.obs['ind_cov'][adata.obs['ind_cov']=='1791'] = '1791_1791'\n",
    "adata.obs['ind_cov'][adata.obs['ind_cov']=='1240'] = '904425200_904425200'\n",
    "adata.obs['ind_cov'][adata.obs['ind_cov']=='FLARE003'] = '904194200_904194200'\n",
    "adata.obs['ind_cov'][adata.obs['ind_cov']=='FLARE008'] = '1763_1763'"
   ]
  },
  {
   "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": true
   },
   "outputs": [],
   "source": [
    "sc.tl.rank_genes_groups(adata, 'leiden', method='t-test', n_genes=25)\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(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Label cell populations"
   ]
  },
  {
   "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\"] = \"B Cells\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"1\"] = \"Cytotoxic T Cells\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"2\"] = \"Natural Killer Cells\" ## good                                                                                     \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"3\"] = \"Naive T Cells\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"4\"] = \"Classical Monocytes\" ## good    Effector Memory T Cells                                                                                 \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"5\"] = \"Classical Monocytes\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"6\"] = \"Classical Monocytes\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"7\"] = \"Cytotoxic T Cells\" ## good        Nonclassical Monocytes                                                                            \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"8\"] = \"Effector Memory T Cells\" ## good     Classical Dendritic Cells                                                                               \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"9\"] = \"Classical Monocytes\" ## good                                                                                    \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"10\"] = \"B Cells\" ## good     Megakaryocytes                                                                             \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"11\"] = \"Naive T Cells\" ## good Plasmacytoid Dendritic Cells                                                                                   \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"12\"] = \"Nonclassical Monocytes\" ## good                                                                                   \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"13\"] = \"Classical Monocytes\" ## good  Macro    Proliferating T Cells                                                                             \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"14\"] = \"Doublets\" ## good                                                                                  \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"15\"] = \"Naive T Cells\" ## good                                                                                   \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"16\"] = \"T-Reg T Cells\" ## good                                                                                   \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"17\"] = \"Doublets\" ## good                                                                                  \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"18\"] = \"Cytotoxic T Cells\" ## good  \n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"19\"] = \"Naive T Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"20\"] = \"Naive T Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"21\"] = \"Naive T Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"22\"] = \"Classical Dendritic Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"23\"] = \"Natural Killer Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"24\"] = \"Naive T Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"25\"] = \"Nonclassical Monocytes\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"26\"] = \"Classical Monocytes\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"27\"] = \"Progenitor Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"28\"] = \"Effector Memory T Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"29\"] = \"Proliferating T Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"30\"] = \"RBCs\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"31\"] = \"Classical Monocytes\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"32\"] = \"Plasmacytoid Dendritic Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"33\"] = \"Naive T Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"34\"] = \"B Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"35\"] = \"Naive T Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"36\"] = \"Naive T Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"37\"] = \"RBCs\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"38\"] = \"Naive T Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"39\"] = \"Megakaryocytes\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"40\"] = \"Naive T Cells\" ## good\n",
    "adata.obs['ct_cov'].loc[adata.obs.leiden == \"41\"] = \"B Cells\" ## good\n",
    "adata.obs['ct_cov'] = adata.obs.ct_cov.astype('category')"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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', 'T-Reg T Cells', 'Cytotoxic T Cells', 'Proliferating T Cells', 'Natural Killer Cells', 'B Cells', 'Progenitor Cells', 'Megakaryocytes', 'RBCs', 'Doublets']\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": [
    "# Remove RBCs and Doublets from further analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = adata[~adata.obs['ct_cov'].isin(['RBCs', 'Doublets'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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', 'T-Reg T Cells','Cytotoxic T Cells', 'Proliferating T Cells', 'Natural Killer Cells', 'B Cells', 'Progenitor Cells', 'Megakaryocytes']\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": [
    "# 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",
    "getproportions(celltype='Classical Monocytes',group1='Healthy', group2='Flare')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get differential gene expression for each cell population in psuedobulk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_pseudobulk(adata):\n",
    "    genes = adata.raw.var_names.tolist()\n",
    "    people = np.unique(adata.obs['ind_cov'].tolist())\n",
    "    pseudobulk = pd.DataFrame(index=genes, columns=people)\n",
    "    for ii in range(len(people)):\n",
    "        fdata = adata[adata.obs['ind_cov']==people[ii]]\n",
    "        pseudobulk[people[ii]] = np.ravel(np.sum(fdata.raw.X, axis=0)/len(fdata))\n",
    "    return pseudobulk\n",
    "\n",
    "def compute_LogFC(pseudobulk, SLE_people, healthy_people):\n",
    "    logFC = {}; pval  = {}\n",
    "    genes = pseudobulk.transpose().keys().tolist()\n",
    "    for ii in range(len(genes)):\n",
    "        h = pseudobulk[healthy_people].transpose()[genes[ii]].values\n",
    "        s = pseudobulk[SLE_people].transpose()[genes[ii]].values\n",
    "        _, p = stats.ranksums(s,h)\n",
    "        pval[genes[ii]] = p\n",
    "        s = np.nanmean(s)\n",
    "        h = np.nanmean(h)\n",
    "        logFC[genes[ii]] = np.log2(np.divide(s+10**-3,h+10**-3))\n",
    "    return logFC, pval\n",
    "\n",
    "def make_volcanoplot(logFC, pval, title):\n",
    "    num_std = 6 # Number of standard deviations above/below mean for gene log2 FC cutoff.\n",
    "    pval_cutoff = 9 # 10^-n (0.001/~39000) allows for 1/1000 false alarm rate corrected for # genes\n",
    "    s = np.std(np.asarray(list(logFC.values())))\n",
    "    upper = np.mean(np.asarray(list(logFC.values()))) + (np.std(np.asarray(list(logFC.values())))*num_std)\n",
    "    lower = np.mean(np.asarray(list(logFC.values()))) - (np.std(np.asarray(list(logFC.values())))*num_std)\n",
    "    genes = np.asarray(list(logFC.keys()))\n",
    "    x = np.asarray(list(logFC.values()))\n",
    "    y = np.log10(list(pval.values()))*-1\n",
    "    sc.settings.set_figure_params(dpi=150, dpi_save=300, format='png', frameon=False, transparent=True, fontsize=12)\n",
    "    # Plot non-significant values\n",
    "    sns.scatterplot(x[(x<=upper) & (x>=lower) & (y<=pval_cutoff)], y[(x<=upper) & (x>=lower) & (y<=pval_cutoff)], color=[0.5, 0.5, 0.5], alpha=0.2)\n",
    "    sns.scatterplot(x[(x<upper) & (y>pval_cutoff)], y[(x<upper) & (y>pval_cutoff)], color=[0.5, 0.5, 0.5], alpha=0.2)\n",
    "    sns.scatterplot(x[(x>lower) & (y>pval_cutoff)], y[(x>lower) & (y>pval_cutoff)], color=[0.5, 0.5, 0.5], alpha=0.2)\n",
    "    # Plot significant values\n",
    "    sns.scatterplot(x[(x>upper) & (y>pval_cutoff)], y[(x>upper) & (y>pval_cutoff)], color='#a1caf1')\n",
    "    sns.scatterplot(x[(x<lower) & (y>pval_cutoff)], y[(x<lower) & (y>pval_cutoff)], color='#a1caf1')\n",
    "    plt.xlabel('Log2 Fold Change')\n",
    "    plt.ylabel('-log10 p-value')\n",
    "    plt.title('Volcano plot: {}'.format(title))\n",
    "    for ii in range(len(genes[(x>upper) & (y>pval_cutoff)])):\n",
    "        plt.text(x[(x>upper) & (y>pval_cutoff)][ii], y[(x>upper) & (y>pval_cutoff)][ii], s=genes[(x>upper) & (y>pval_cutoff)][ii], fontsize=5)\n",
    "\n",
    "    for ii in range(len(genes[(x<lower) & (y>pval_cutoff)])):\n",
    "        plt.text(x[(x<lower) & (y>pval_cutoff)][ii], y[(x<lower) & (y>pval_cutoff)][ii], s=genes[(x<lower) & (y>pval_cutoff)][ii], fontsize=5)\n",
    "    \n",
    "    URG = genes[(x>upper) & (y>pval_cutoff)]\n",
    "    DRG = genes[(x<lower) & (y>pval_cutoff)]\n",
    "    return URG, DRG\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Log2 FC for Total PBMC psuedobulk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_totalPBMC_log2FC(adata):\n",
    "    healthy_people = np.unique(adata.obs['ind_cov'][adata.obs['SLE status']=='Healthy'].tolist())\n",
    "    SLE_people = np.unique(adata.obs['ind_cov'][adata.obs['SLE status']=='SLE'].tolist())\n",
    "    pseudobulk_totalPBMC = get_pseudobulk(adata)\n",
    "    logFC_totalPBMC, pval_totalPBMC = compute_LogFC(pseudobulk_totalPBMC, SLE_people, healthy_people)\n",
    "    data = pd.DataFrame(list(logFC_totalPBMC.values()), index=list(logFC_totalPBMC.keys()), columns=['Log2FC'])\n",
    "    data2 = pd.DataFrame(list(pval_totalPBMC.values()), index=list(pval_totalPBMC.keys()), columns=['p-value'])\n",
    "    data = pd.concat([data, data2], axis=1)\n",
    "    data = data.transpose()\n",
    "    data.to_csv('Log2FC-SLEvsHealthy-totalPBMC.csv')\n",
    "    return logFC_totalPBMC, pval_totalPBMC\n",
    "\n",
    "#logFC_totalPBMC, pval_totalPBMC = get_totalPBMC_log2FC(adata)\n",
    "\n",
    "tmp = pd.read_csv('Log2FC-SLEvsHealthy-totalPBMC.csv', index_col=0)\n",
    "logFC_totalPBMC = tmp.transpose()['Log2FC'].to_dict()\n",
    "pval_totalPBMC = tmp.transpose()['p-value'].to_dict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "URG_total, DRG_total = make_volcanoplot(logFC_totalPBMC, pval_totalPBMC, title='Total PBMCs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Log2 FC for each subpopulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_subpop_log2FC(adata):\n",
    "    ct_groups = np.unique(adata.obs['ct_cov'].tolist())\n",
    "    logFC_groups = {}; pval_groups = {};\n",
    "    for ii in range(len(ct_groups)):\n",
    "        print('Cell type: {}'.format(ct_groups[ii]))\n",
    "        fdata = adata[adata.obs['ct_cov']==ct_groups[ii]]\n",
    "        healthy_people = np.unique(fdata.obs['ind_cov'][fdata.obs['SLE status']=='Healthy'].tolist())\n",
    "        SLE_people = np.unique(fdata.obs['ind_cov'][fdata.obs['SLE status']=='SLE'].tolist())\n",
    "        pseudobulk = get_pseudobulk(fdata)\n",
    "        logFC, pval = compute_LogFC(pseudobulk, SLE_people, healthy_people)\n",
    "        logFC_groups[ct_groups[ii]] = logFC\n",
    "        pval_groups[ct_groups[ii]] = pval\n",
    "        data = pd.DataFrame(list(logFC_groups[ct_groups[ii]].values()), index=list(logFC_groups[ct_groups[ii]].keys()), columns=['Log2FC'])\n",
    "        data2 = pd.DataFrame(list(pval_groups[ct_groups[ii]].values()), index=list(pval_groups[ct_groups[ii]].keys()), columns=['p-value'])\n",
    "        data = pd.concat([data, data2], axis=1)\n",
    "        data = data.transpose()\n",
    "        data.to_csv('Log2FC-SLEvsHealthy-{}.csv'.format(ct_groups[ii]))\n",
    "        \n",
    "def load_subpop_log2FC(adata):\n",
    "    ct_groups = np.unique(adata.obs['ct_cov'].tolist())\n",
    "    logFC_groups = {}; pval_groups = {};\n",
    "    for ii in range(len(ct_groups)):\n",
    "        tmp = pd.read_csv('Log2FC-SLEvsHealthy-{}.csv'.format(ct_groups[ii]), index_col=0)\n",
    "        logFC_groups[ct_groups[ii]] = tmp.transpose()['Log2FC'].to_dict()\n",
    "        pval_groups[ct_groups[ii]] = tmp.transpose()['p-value'].to_dict()\n",
    "    return logFC_groups, pval_groups\n",
    "    \n",
    "#run_subpop_log2FC(adata)\n",
    "\n",
    "logFC_groups, pval_groups = load_subpop_log2FC(adata)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "URG = {}; DRG = {}; # Up and Down Regulated Genes for each population\n",
    "ct_groups = np.unique(adata.obs['ct_cov'].tolist())\n",
    "for ii in range(len(ct_groups)):\n",
    "    plt.figure()\n",
    "    URG[ct_groups[ii]], DRG[ct_groups[ii]] = make_volcanoplot(logFC_groups[ct_groups[ii]], pval_groups[ct_groups[ii]], ct_groups[ii])\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute Gene-Cell-Type specific expression per individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pseudobulk = get_pseudobulk(adata)\n",
    "corr_matrix = pd.DataFrame(index=list(pseudobulk.columns))\n",
    "for ii in range(len(URG_total)):\n",
    "    corr_matrix['{}-URG-totalPBMC'.format(URG_total[ii])] = pseudobulk.transpose()[URG_total[ii]].values\n",
    "for ii in range(len(DRG_total)):\n",
    "    corr_matrix['{}-DRG-totalPBMC'.format(DRG_total[ii])] = pseudobulk.transpose()[DRG_total[ii]].values\n",
    "    \n",
    "ct_groups = np.unique(adata.obs['ct_cov'].tolist())\n",
    "for ii in range(len(ct_groups)):\n",
    "    try:\n",
    "        print('Cell type: {}'.format(ct_groups[ii]))\n",
    "        fdata = adata[adata.obs['ct_cov']==ct_groups[ii]]\n",
    "        pseudobulk = get_pseudobulk(fdata)\n",
    "        for gg in range(len(URG[ct_groups[ii]])):   \n",
    "            corr_matrix['{0}-URG-{1}'.format(URG[ct_groups[ii]][gg], ct_groups[ii])] = pseudobulk.transpose()[URG[ct_groups[ii]][gg]].values\n",
    "        for gg in range(len(DRG[ct_groups[ii]])):   \n",
    "            corr_matrix['{0}-DRG-{1}'.format(DRG[ct_groups[ii]][gg], ct_groups[ii])] = pseudobulk.transpose()[DRG[ct_groups[ii]][gg]].values\n",
    "    except:\n",
    "        continue\n",
    "corr_matrix.to_csv('Expression-Composition-corrmatrix.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make signatures for cell type specific expression without pan signature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "ct_groups = np.unique(adata.obs['ct_cov'].tolist())\n",
    "URG_exclude = {}; DRG_exclude = {}\n",
    "for ii in range(len(ct_groups)):\n",
    "    listURG = list(URG[ct_groups[ii]])\n",
    "    listDRG = list(DRG[ct_groups[ii]])\n",
    "    removallist = list(set(URG[ct_groups[ii]]).intersection(URG_total))\n",
    "    removallistD = list(set(DRG[ct_groups[ii]]).intersection(DRG_total))\n",
    "    for item in removallist:\n",
    "        listURG.remove(item)\n",
    "    for item in removallistD:\n",
    "        listDRG.remove(item)\n",
    "    URG_exclude[ct_groups[ii]] = listURG\n",
    "    DRG_exclude[ct_groups[ii]] = listDRG\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute averaged gene expression of significant genes per person "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "pseudobulk = get_pseudobulk(adata)\n",
    "corr_matrix = pd.DataFrame(index=list(pseudobulk.columns))\n",
    "corr_matrix['URG-totalPBMC'] = pseudobulk.transpose()[URG_total].mean(axis=1).values\n",
    "corr_matrix['DRG-totalPBMC'] = pseudobulk.transpose()[DRG_total].mean(axis=1).values\n",
    "\n",
    "\n",
    "ct_groups = np.unique(adata.obs['ct_cov'].tolist())\n",
    "for ii in range(len(ct_groups)):\n",
    "    try:\n",
    "        print('Cell type: {}'.format(ct_groups[ii]))\n",
    "        fdata = adata[adata.obs['ct_cov']==ct_groups[ii]]\n",
    "        pseudobulk = get_pseudobulk(fdata)\n",
    "        #corr_matrix['URG-{}'.format(ct_groups[ii])] = pseudobulk.transpose()[URG[ct_groups[ii]]].mean(axis=1).values\n",
    "        #corr_matrix['DRG-{}'.format(ct_groups[ii])] = pseudobulk.transpose()[DRG[ct_groups[ii]]].mean(axis=1).values\n",
    "        corr_matrix['URG-{}'.format(ct_groups[ii])] = pseudobulk.transpose()[URG_exclude[ct_groups[ii]]].mean(axis=1).values\n",
    "        corr_matrix['DRG-{}'.format(ct_groups[ii])] = pseudobulk.transpose()[DRG_exclude[ct_groups[ii]]].mean(axis=1).values\n",
    "    except:\n",
    "        continue\n",
    "corr_matrix.to_csv('Expression-Composition-corrmatrix.csv')\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_matrix = pd.read_csv('Expression-Composition-corrmatrix.csv', index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option('display.max_columns', None)\n",
    "display(corr_matrix)\n",
    "# Clustermap this gene expression\n",
    "# Average pan expression genes across cell types and within cell types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# genes by cell type + total PBMC corr matrix. Be sure to average the people after pseudobulk. Split by healthy and SLE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Corrlation plot of corrmatrix expression values. Combine over pan gene PBMC and per cell type\n",
    "sc.settings.set_figure_params(dpi=250, dpi_save=300, format='png', frameon=False, transparent=True, fontsize=12)\n",
    "sns.set(font_scale=0.2)\n",
    "sns.clustermap(corr_matrix.dropna().transpose(), standard_scale=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute composition of each cell type per person"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Make proportion plots\n",
    "adata_obs_small = adata.obs\n",
    "ind_count = adata_obs_small.groupby(['ind_cov', 'ct_cov','SLE status'])['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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add composition changes\n",
    "people = corr_matrix.index.tolist()\n",
    "ct_groups = np.unique(adata.obs['ct_cov'].tolist())\n",
    "for jj in range(len(ct_groups)):\n",
    "    tmp = ind_perc[ind_perc['ct_cov']==ct_groups[jj]]\n",
    "    tmplist = np.empty((len(people),))\n",
    "    tmplist[:] = np.nan\n",
    "    for ii in range(len(people)):\n",
    "        try:\n",
    "            tmplist[ii] = tmp[tmp['ind_cov']==people[ii]]['ct_perc'].values[0]\n",
    "        except:\n",
    "            continue\n",
    "    corr_matrix['Proportion-{}'.format(ct_groups[jj])] = tmplist\n",
    "    \n",
    "# Remove empty categories    \n",
    "labels = list(corr_matrix.columns)\n",
    "for lbl in labels:\n",
    "    if np.sum(np.isnan(corr_matrix[lbl]))==len(corr_matrix):\n",
    "        print('Removing category: {} due to lack of data.'.format(lbl))\n",
    "        corr_matrix = corr_matrix.drop(lbl, axis='columns')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute correlations between URG, DRGs and cell type composition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlations = corr_matrix.corr(method='spearman')\n",
    "correlations[np.isnan(correlations)] = 0\n",
    "sns.set(font_scale=0.2)\n",
    "sns.clustermap(np.abs(correlations))\n",
    "plt.savefig(\"Correlation_plot.png\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_barplot(correlations, label):\n",
    "    sc.settings.set_figure_params(dpi=150, dpi_save=200, format='png', frameon=False, transparent=True, fontsize=5)\n",
    "    x = np.asarray(list((correlations[label].values)))\n",
    "    indices = np.flipud(np.argsort(x)); x = x[indices];\n",
    "    y = np.asarray(list(correlations[label].index)); y = y[indices]\n",
    "    x = np.concatenate((x[:20], x[-20:]))\n",
    "    y = np.concatenate((y[:20], y[-20:]))\n",
    "    g = sns.barplot(x, y)\n",
    "    plt.title('Correlations to {}'.format(label))\n",
    "\n",
    "def plot_gene_expression(cdata, celltypes, gene, normalize):\n",
    "    sns.set(font_scale=1)\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",
    "            normvalues = np.divide(np.asarray(values), np.asarray(count)) # Normalize by cell number\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",
    "    if normalize == True:\n",
    "        genexpression = pd.DataFrame(data={'Unique_ID': ind_id, 'disease_cov': dz_id, gene: normvalues, 'Cell_Count':count, 'ct_cov':ct_id})    \n",
    "    else:\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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='ISG15-URG-totalPBMC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='Proportion-Naive T Cells')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='NOG-DRG-totalPBMC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='Proportion-Plasmacytoid Dendritic Cells')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='Proportion-Classical Dendritic Cells')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='Proportion-Classical Monocytes')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='Proportion-Cytotoxic T Cells')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='Proportion-Proliferating T Cells')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ct_groups = np.unique(adata.obs['ct_cov'].tolist())\n",
    "genes = ['ISG15', 'IFI6', 'NOG']\n",
    "for gene in genes:\n",
    "    plot_gene_expression(adata, celltypes=['Classical Monocytes', 'Naive T Cells'], gene=gene, normalize=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "ct_groups = np.unique(adata.obs['ct_cov'].tolist())\n",
    "genes = ['ISG15', 'IFI6', 'NOG']\n",
    "for gene in genes:\n",
    "    plot_gene_expression(adata, celltypes=['Classical Monocytes', 'Naive T Cells'], gene=gene, normalize=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clinical covariate analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clincal = pd.read_csv('lupus.sledai.txt', sep='\\t')\n",
    "# Add study disease status\n",
    "disease_status = {}\n",
    "for ii in range(len(clincal['genotypeid'])):\n",
    "    try:\n",
    "        disease_status[clincal['genotypeid'][ii]] = 1 #np.unique(adata.obs['SLE status'][adata.obs['ind_cov']==clincal['genotypeid'][ii]])[0]\n",
    "    except:\n",
    "        disease_status[clincal['genotypeid'][ii]] = 1 #'SLE' # This person by definition of being in sheet has SLE\n",
    "        #print('Individual: {} not in study data'.format(clincal['genotypeid'][ii]))\n",
    "\n",
    "clincal['Disease_status'] = disease_status.values()\n",
    "display(clincal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove SLE patients that have no clinical data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_matrix_mod = corr_matrix[~corr_matrix.transpose().columns.isin(['FLARE002', 'FLARE004', 'FLARE005', 'FLARE006', \n",
    "                                                                     'FLARE007', 'FLARE009', 'FLARE011', 'FLARE013', \n",
    "                                                                     'FLARE014', 'FLARE015', 'FLARE016', 'FLARE017', \n",
    "                                                                     'FLARE018', 'FLARE019', 'FLARE020'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_covar(corr_matrix_mod, keys):\n",
    "    clincal = pd.read_csv('lupus.sledai.txt', sep='\\t')\n",
    "    inds = list(corr_matrix_mod.transpose().columns)\n",
    "    # Add study disease status\n",
    "    disease_status = {}\n",
    "    for ii in range(len(clincal['genotypeid'])):\n",
    "        try:\n",
    "            disease_status[clincal['genotypeid'][ii]] = 1 #np.unique(adata.obs['SLE status'][adata.obs['ind_cov']==clincal['genotypeid'][ii]])[0]\n",
    "        except:\n",
    "            disease_status[clincal['genotypeid'][ii]] = 1 #'SLE' # This person by definition of being in sheet has SLE\n",
    "            #print('Individual: {} not in study data'.format(clincal['genotypeid'][ii]))\n",
    "\n",
    "    clincal['Disease_status'] = disease_status.values()\n",
    "    for key in keys:\n",
    "        covar = {}\n",
    "        for ii in range(len(inds)):\n",
    "            if np.sum(clincal['genotypeid'].isin([inds[ii]]).values):\n",
    "                covar[inds[ii]] = clincal[clincal['genotypeid'].isin([inds[ii]]).values][key].values[0]\n",
    "            else:\n",
    "                covar[inds[ii]] = 0\n",
    "        corr_matrix_mod[key] = list(covar.values())\n",
    "    return corr_matrix_mod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_matrix_mod = build_covar(corr_matrix_mod, keys = ['SLAQ score', 'sleactivity', 'SLIC score', 'acrcsum', 'sledaiscore'])\n",
    "display(corr_matrix_mod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlations = corr_matrix_mod.corr(method='spearman')\n",
    "correlations[np.isnan(correlations)] = 0\n",
    "sns.set(font_scale=0.2)\n",
    "sns.clustermap(np.abs(correlations))\n",
    "plt.savefig(\"Correlation+clinical_plot.png\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='SLAQ score')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='sleactivity')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='SLIC score')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='sledaiscore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplot(correlations, label='acrcsum')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(font_scale=1)\n",
    "sns.regplot(corr_matrix_mod['NOG-DRG-totalPBMC'].values, corr_matrix_mod['SLAQ score'].values)\n",
    "#sns.regplot(corr_matrix_mod['NOG-DRG-Naive T Cells'].values, corr_matrix_mod['SLAQ score'].values)\n",
    "plt.xlabel('NOG expression')\n",
    "plt.ylabel('SLAQ score')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(font_scale=1)\n",
    "sns.regplot(corr_matrix_mod['NOG-DRG-totalPBMC'].values, corr_matrix_mod['acrcsum'].values)\n",
    "#sns.regplot(corr_matrix_mod['NOG-DRG-Naive T Cells'].values, corr_matrix_mod['SLAQ score'].values)\n",
    "plt.xlabel('NOG expression')\n",
    "plt.ylabel('SLAQ score')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic regression on clinical variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logistic_regression(corr_matrix, clincal_var, n_boot):\n",
    "    # Remove people with no clinical variables\n",
    "    corr_matrix_logreg = corr_matrix[~corr_matrix.transpose().columns.isin(['FLARE002', 'FLARE004', 'FLARE005', 'FLARE006', \n",
    "                                                                         'FLARE007', 'FLARE009', 'FLARE011', 'FLARE013', \n",
    "                                                                         'FLARE014', 'FLARE015', 'FLARE016', 'FLARE017', \n",
    "                                                                         'FLARE018', 'FLARE019', 'FLARE020'])]\n",
    "    # Add clinical outcome\n",
    "    corr_matrix_logreg = build_covar(corr_matrix_logreg, keys = [clincal_var])\n",
    "    # Drop NaN from dataframe\n",
    "    corr_matrix_logreg = corr_matrix_logreg.dropna()\n",
    "    # Dependent variable\n",
    "    Y = corr_matrix_logreg[clincal_var].values.ravel()\n",
    "    # Remove dependent variable from variables\n",
    "    corr_matrix_logreg = corr_matrix_logreg.drop(clincal_var, axis=1)\n",
    "    # observations x features\n",
    "    X = corr_matrix_logreg.as_matrix()\n",
    "    # Variables of interest\n",
    "    VOIs = {}; ROC_AUC = {}; FPR = {}; TPR = {}\n",
    "    for ii in range(n_boot):\n",
    "        # Split data\n",
    "        X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.2)\n",
    "\n",
    "        # Over sample cases and undersample controls to improve classification\n",
    "        os = SMOTE(random_state=0)\n",
    "        os_data_X, os_data_y = os.fit_sample(X_train, y_train)\n",
    "\n",
    "        # Setup logistic regression\n",
    "        logreg = LogisticRegression(max_iter=5000, penalty='elasticnet', solver='saga', l1_ratio=0.5)\n",
    "        # Find optimal parameters\n",
    "        #rfe = RFE(logreg, 3)\n",
    "        rfe = SelectFromModel(logreg, threshold='10*mean')\n",
    "        rfe = rfe.fit(os_data_X, os_data_y)\n",
    "        ind_var = np.asarray(list(corr_matrix_logreg.columns))\n",
    "\n",
    "        # Preserve optimate features only from train and test set\n",
    "        #os_data_X = os_data_X[:, corr_matrix_logreg.columns.isin(ind_var[rfe.ranking_==1])]\n",
    "        #X_test = X_test[:, corr_matrix_logreg.columns.isin(ind_var[rfe.ranking_==1])]\n",
    "        os_data_X = os_data_X[:, rfe.get_support()]\n",
    "        X_test = X_test[:, rfe.get_support()]\n",
    "        # Fit Logistic regression\n",
    "        logreg = LogisticRegression(max_iter=5000, penalty='elasticnet', solver='saga', l1_ratio=0.5)\n",
    "        logreg.fit(os_data_X, os_data_y)\n",
    "        #VOIs[ii] = list(ind_var[rfe.ranking_==1])\n",
    "        VOIs[ii] = list(ind_var[rfe.get_support()])\n",
    "        # Perform prediction\n",
    "        y_pred = logreg.predict(X_test)\n",
    "\n",
    "        # ROC curve\n",
    "        logit_roc_auc = roc_auc_score(y_test, logreg.predict(X_test))\n",
    "        ROC_AUC[ii] = logit_roc_auc \n",
    "        # Collect True and False Positive rates\n",
    "        fpr, tpr, thresholds = roc_curve(y_test, logreg.predict_proba(X_test)[:,1])\n",
    "        FPR[ii] = fpr\n",
    "        TPR[ii] = tpr\n",
    "    return VOIs, ROC_AUC, FPR, TPR\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_ROC(ROC_AUC, FPR, TPR):\n",
    "    plt.figure()\n",
    "    for ii in range(len(ROC_AUC)):\n",
    "        plt.plot(np.asarray(list(FPR.values()))[ii], np.asarray(list(TPR.values()))[ii], 'k-', alpha=0.5, linewidth=.25)\n",
    "    plt.plot([0, 1], [0, 1],'r--')\n",
    "    plt.xlim([0.0, 1.0])\n",
    "    plt.ylim([0.0, 1.05])\n",
    "    plt.xlabel('False Positive Rate')\n",
    "    plt.ylabel('True Positive Rate')\n",
    "    plt.title('ROC AUC average: {0}, median:{1}'.format(round(np.mean(list(ROC_AUC.values())),3), round(np.median(list(ROC_AUC.values())),3)))\n",
    "    plt.show()\n",
    "    \n",
    "def show_VOI(VOIs):\n",
    "    from palettable.colorbrewer.qualitative import Pastel1_7\n",
    "    flat_list = [item for sublist in list(VOIs.values()) for item in sublist]\n",
    "    unique_VOIs = np.unique(flat_list)\n",
    "    #unique_VOIs = np.unique(np.ravel(list(VOIs.values())))\n",
    "    VOI_count = {}\n",
    "    for ii in range(len(unique_VOIs)):\n",
    "        VOI_count[unique_VOIs[ii]] = np.sum(np.asarray(flat_list) == unique_VOIs[ii])\n",
    "    def my_autopct(pct):\n",
    "        return ('%1.1f%%' % pct) if pct > 5 else ''\n",
    "    labels = np.asarray(list(VOI_count.keys()))\n",
    "    labels[list(VOI_count.values()) < (np.sum(list(VOI_count.values()))*0.05)] = ''\n",
    "    my_circle=plt.Circle( (0,0), 0.7, color='white')\n",
    "    plt.pie(list(VOI_count.values()), labels=labels, colors=Pastel1_7.hex_colors, textprops={'fontsize': 7}, autopct=my_autopct)\n",
    "    p=plt.gcf()\n",
    "    p.gca().add_artist(my_circle)\n",
    "    plt.show()\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cases only\n",
    "corr_matrix_cases = corr_matrix[corr_matrix.transpose().columns.isin(np.unique(adata.obs['ind_cov'][adata.obs['SLE status']=='SLE'].tolist()))]\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_matrix\n",
    "# Cluster genes in expression/cluster this"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kidney involvment\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix_cases, clincal_var='kidney', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Malar rash"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix_cases, clincal_var='acrmalar', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discoid rash"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix_cases, clincal_var='acrdiscoid', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Photosensitivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix_cases, clincal_var='acrphotosensitiv', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mucosal ulcers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix_cases, clincal_var='acrmuculcers', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Arthritis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix_cases, clincal_var='acrarthritisc', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Serositis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix_cases, clincal_var='acrserositisc', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)\n",
    "# Why does RP4. Do cases analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lupus Nephritis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix_cases, clincal_var='acrlupneph', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Renal involvement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix_cases, clincal_var='acrrenalc', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# anti-ANA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix, clincal_var='acranac', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# anti-Smith"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix, clincal_var='acrantismith', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# anti-dsDNA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix_cases, clincal_var='acrantidsdna', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# anti-phospholipid antigen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix_cases, clincal_var='acrapla', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lymphopenia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR, = logistic_regression(corr_matrix_cases, clincal_var='acrlymphopenia', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Leukopenia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix, clincal_var='acrleukopenia', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thrombocytopenia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "VOIs, ROC_AUC, FPR, TPR = logistic_regression(corr_matrix, clincal_var='acrthrombocyto', n_boot=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear regression on sum of logistic regression predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def run_logreg_bank(corr_matrix, clincal, lupus_var, final_y, split):\n",
    "    acrdata = clincal.transpose()[clincal.columns.isin(lupus_var)].transpose()\n",
    "\n",
    "    # Remove people with no clinical variables\n",
    "    corr_matrix_logreg = corr_matrix[~corr_matrix.transpose().columns.isin(['FLARE002', 'FLARE004', 'FLARE005', 'FLARE006', \n",
    "                                                                         'FLARE007', 'FLARE009', 'FLARE011', 'FLARE013', \n",
    "                                                                         'FLARE014', 'FLARE015', 'FLARE016', 'FLARE017', \n",
    "                                                                         'FLARE018', 'FLARE019', 'FLARE020'])]\n",
    "    # List of ACR clinical covariates\n",
    "    clin_var = list(acrdata.columns) \n",
    "    # Add clinical covariates to dataframe\n",
    "    for clin in clin_var:\n",
    "        corr_matrix_logreg = build_covar(corr_matrix_logreg, keys = [clin])\n",
    "    # Drop people with incomplete data from dataframe.\n",
    "    corr_matrix_logreg = corr_matrix_logreg.dropna()\n",
    "    # ACR sum - dependent variable for linear regression\n",
    "    Y = corr_matrix_logreg[final_y].values.ravel()\n",
    "    # Remove dependent variable for linear regression from dataframe and clinical covariates\n",
    "    corr_matrix_logreg = corr_matrix_logreg.drop(final_y, axis=1)\n",
    "    clin_var.remove(final_y)\n",
    "\n",
    "    # Split data\n",
    "    X_train, X_test, y_linear_train, y_linear_test = train_test_split(corr_matrix_logreg, Y, test_size=split)\n",
    "    # Generate dictionaries to store outcome of each logistic regression\n",
    "    VOIs = {}; ROC_AUC = {}; FPR = {}; TPR = {}; Y_logreg_pred ={}\n",
    "    for lr in range(len(clin_var)):\n",
    "        try:\n",
    "            clin_tmp = clin_var.copy()\n",
    "            print('Clinical Variable: {}'.format(clin_tmp[lr]))\n",
    "            # Remove all dependent variables except one of interest.\n",
    "            clin_tmp.remove(clin_tmp[lr])\n",
    "            X_train_tmp = X_train.drop(clin_tmp, axis=1)\n",
    "            X_test_tmp  = X_test.drop(clin_tmp, axis=1)\n",
    "            # Dependent variable for logistic regression\n",
    "            Y_logreg_train = X_train_tmp[clin_var[lr]].values.ravel()\n",
    "            Y_logreg_test  = X_test_tmp[clin_var[lr]].values.ravel()\n",
    "            # Remove dependent varaible for logistic regression from dataframe\n",
    "            X_train_tmp = X_train_tmp.drop(clin_var[lr], axis=1)\n",
    "            X_test_tmp = X_test_tmp.drop(clin_var[lr], axis=1)\n",
    "            # Over sample cases and undersample controls to improve classification\n",
    "            os = SMOTE(random_state=0)\n",
    "            os_data_X, os_data_y = os.fit_sample(X_train_tmp, Y_logreg_train)\n",
    "            # Setup logistic regression\n",
    "            logreg = LogisticRegression(max_iter=2000)\n",
    "            # Find optimal parameters\n",
    "            rfe = RFE(logreg, 3)\n",
    "            rfe = rfe.fit(os_data_X, os_data_y)\n",
    "            ind_var = np.asarray(list(os_data_X.columns))\n",
    "\n",
    "            # Preserve optimate features only from train and test set\n",
    "            os_data_X = os_data_X[os_data_X.columns[os_data_X.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "            X_test_tmp = X_test_tmp[os_data_X.columns[os_data_X.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "\n",
    "            # Fit Logistic regression\n",
    "            logreg = LogisticRegression(max_iter=2000)\n",
    "            logreg.fit(os_data_X, os_data_y)\n",
    "            VOIs[clin_var[lr]] = list(ind_var[rfe.ranking_==1])\n",
    "            # Perform prediction\n",
    "            y_pred = logreg.predict(X_test_tmp)\n",
    "            Y_logreg_pred[clin_var[lr]] = y_pred\n",
    "            # ROC curve\n",
    "            logit_roc_auc = roc_auc_score(Y_logreg_test, y_pred)\n",
    "            ROC_AUC[clin_var[lr]] = logit_roc_auc \n",
    "            # Collect True and False Positive rates\n",
    "            fpr, tpr, thresholds = roc_curve(Y_logreg_test, logreg.predict_proba(X_test_tmp)[:,1])\n",
    "            FPR[clin_var[lr]] = fpr\n",
    "            TPR[clin_var[lr]] = tpr\n",
    "            print('ROC_AOC: {}'.format(ROC_AUC[clin_var[lr]]))\n",
    "            print('Variables: {}'.format(VOIs[clin_var[lr]]))\n",
    "        except:\n",
    "            print('Unable to model variable: {}'.format(clin_tmp[lr]))\n",
    "            continue\n",
    "    return VOIs, Y_logreg_pred, y_linear_test, ROC_AUC, FPR, TPR\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ACR Sum prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "final_y = 'acrcsum'\n",
    "# Keep only ACR data\n",
    "lupus_var = ['acrarthritisc','acrantidsdna', 'acrantismith', 'acrlupneph', 'acranac', final_y]\n",
    "VOIs, Y_logreg_pred, y_linear_test, ROC_AUC, FPR, TPR = run_logreg_bank(corr_matrix, clincal, lupus_var, final_y, split=0.2)\n",
    "\n",
    "model = sm.OLS(pd.DataFrame.from_dict(Y_logreg_pred).transpose().sum(axis=0).values, y_linear_test)\n",
    "results = model.fit()\n",
    "print(results.summary())\n",
    "print('Parameters: ', results.params)\n",
    "sns.regplot(pd.DataFrame.from_dict(Y_logreg_pred).transpose().sum(axis=0).values, y_linear_test)\n",
    "plt.xlabel('Prediction')\n",
    "plt.ylabel(final_y)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SLEDAI prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "final_y = 'sledaiscore'\n",
    "# Keep only ACR data\n",
    "lupus_var = ['acrarthritisc','acrantidsdna', 'acrantismith', 'acrlupneph', 'acranac', final_y]\n",
    "VOIs, Y_logreg_pred, y_linear_test, ROC_AUC, FPR, TPR = run_logreg_bank(corr_matrix, clincal, lupus_var, final_y, split=0.2)\n",
    "\n",
    "model = sm.OLS(pd.DataFrame.from_dict(Y_logreg_pred).transpose().sum(axis=0).values, y_linear_test)\n",
    "results = model.fit()\n",
    "print(results.summary())\n",
    "print('Parameters: ', results.params)\n",
    "sns.regplot(pd.DataFrame.from_dict(Y_logreg_pred).transpose().sum(axis=0).values, y_linear_test)\n",
    "plt.xlabel('Prediction')\n",
    "plt.ylabel(final_y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = sm.OLS(clincal['acrcsum'].values, clincal['sledaiscore'].values)\n",
    "results = model.fit()\n",
    "print(results.summary())\n",
    "print('Parameters: ', results.params)\n",
    "sns.regplot(clincal['acrcsum'].values, clincal['sledaiscore'].values)\n",
    "plt.xlabel('ACR sum')\n",
    "plt.ylabel('SLEDAI')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SLIC prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "final_y = 'SLIC score'\n",
    "# Keep only ACR data\n",
    "lupus_var = ['acrarthritisc','acrantidsdna', 'acrantismith', 'acrlupneph', 'acranac', final_y]\n",
    "VOIs, Y_logreg_pred, y_linear_test, ROC_AUC, FPR, TPR = run_logreg_bank(corr_matrix, clincal, lupus_var, final_y, split=0.2)\n",
    "\n",
    "model = sm.OLS(pd.DataFrame.from_dict(Y_logreg_pred).transpose().sum(axis=0).values, y_linear_test)\n",
    "results = model.fit()\n",
    "print(results.summary())\n",
    "print('Parameters: ', results.params)\n",
    "sns.regplot(pd.DataFrame.from_dict(Y_logreg_pred).transpose().sum(axis=0).values, y_linear_test)\n",
    "plt.xlabel('Prediction')\n",
    "plt.ylabel(final_y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = sm.OLS(clincal['acrcsum'].values, clincal['SLIC score'].values)\n",
    "results = model.fit()\n",
    "print(results.summary())\n",
    "print('Parameters: ', results.params)\n",
    "sns.regplot(clincal['acrcsum'].values, clincal['SLIC score'].values)\n",
    "plt.xlabel('ACR sum')\n",
    "plt.ylabel('SLIC score')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SLAQ prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_y = 'SLAQ score'\n",
    "# Keep only ACR data\n",
    "lupus_var = ['acrarthritisc','acrantidsdna', 'acrantismith', 'acrlupneph', 'acranac', final_y]\n",
    "VOIs, Y_logreg_pred, y_linear_test, ROC_AUC, FPR, TPR = run_logreg_bank(corr_matrix, clincal, lupus_var, final_y, split=0.2)\n",
    "\n",
    "model = sm.OLS(pd.DataFrame.from_dict(Y_logreg_pred).transpose().sum(axis=0).values, y_linear_test)\n",
    "results = model.fit()\n",
    "print(results.summary())\n",
    "print('Parameters: ', results.params)\n",
    "sns.regplot(pd.DataFrame.from_dict(Y_logreg_pred).transpose().sum(axis=0).values, y_linear_test)\n",
    "plt.xlabel('Prediction')\n",
    "plt.ylabel(final_y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = sm.OLS(clincal['acrcsum'].values, clincal['SLAQ score'].values)\n",
    "results = model.fit()\n",
    "print(results.summary())\n",
    "print('Parameters: ', results.params)\n",
    "sns.regplot(clincal['acrcsum'].values, clincal['SLAQ score'].values)\n",
    "plt.xlabel('ACR sum')\n",
    "plt.ylabel('SLAQ score')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ROC curves comparing disease status and associated scale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Patients may have SLE but exhibit no SLE symptoms at the time of blood collection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acrdata = clincal.transpose()[clincal.columns.isin(lupus_var)].transpose()\n",
    "corr_matrix_ROC = pd.DataFrame(index=list(corr_matrix.index))\n",
    "# Remove people with no clinical variables\n",
    "corr_matrix_ROC = corr_matrix_ROC[~corr_matrix_ROC.transpose().columns.isin(['FLARE002', 'FLARE004', 'FLARE005', 'FLARE006', \n",
    "                                                                     'FLARE007', 'FLARE009', 'FLARE011', 'FLARE013', \n",
    "                                                                     'FLARE014', 'FLARE015', 'FLARE016', 'FLARE017', \n",
    "                                                                     'FLARE018', 'FLARE019', 'FLARE020'])]\n",
    "# List of ACR clinical covariates\n",
    "clin_var = list(acrdata.columns)\n",
    "# Add clinical covariates to dataframe\n",
    "for clin in clin_var:\n",
    "    corr_matrix_ROC = build_covar(corr_matrix_ROC, keys = ['Disease_status', 'acrcsum', 'sledaiscore', 'SLIC score', 'SLAQ score', 'sleactivity'])\n",
    "# Drop people with incomplete data from dataframe.\n",
    "corr_matrix_ROC = corr_matrix_ROC.dropna()\n",
    "# ACR sum - dependent variable for linear regression\n",
    "Y = corr_matrix_ROC['Disease_status'].values.ravel()\n",
    "# Remove dependent variable for linear regression from dataframe and clinical covariates\n",
    "corr_matrix_ROC = corr_matrix_ROC.drop('Disease_status', axis=1)\n",
    "clin_var.remove(final_y)\n",
    "\n",
    "# Turn scales into diagnosis outcomes based study cutoffs\n",
    "corr_matrix_ROC['acrcsum'] = ((corr_matrix_ROC['acrcsum'] >=4)*1).values\n",
    "corr_matrix_ROC['sledaiscore'] = ((corr_matrix_ROC['sledaiscore'] >=3)*1).values\n",
    "corr_matrix_ROC['SLIC score'] = ((corr_matrix_ROC['SLIC score'] >=4)*1).values\n",
    "corr_matrix_ROC['SLAQ score'] = ((corr_matrix_ROC['SLAQ score'] >=4)*1).values\n",
    "corr_matrix_ROC['sleactivity'] = ((corr_matrix_ROC['sleactivity'] >=3)*1).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "scales = list(corr_matrix_ROC.columns)\n",
    "for ii in range(corr_matrix_ROC.shape[1]):\n",
    "    # ROC curve\n",
    "    logit_roc_auc = roc_auc_score(Y, corr_matrix_ROC[scales[ii]].values)\n",
    "    # Collect True and False Positive rates\n",
    "    fpr, tpr, thresholds = roc_curve(Y, corr_matrix_ROC[scales[ii]].values)\n",
    "\n",
    "    plt.figure()\n",
    "    plt.plot(fpr, tpr, 'k-', alpha=1, linewidth=2)\n",
    "    plt.plot([0, 1], [0, 1],'r--')\n",
    "    plt.xlim([0.0, 1.0])\n",
    "    plt.ylim([0.0, 1.05])\n",
    "    plt.xlabel('False Positive Rate')\n",
    "    plt.ylabel('True Positive Rate')\n",
    "    plt.title('ROC AUC: {0}, for {1} scale'.format(round(logit_roc_auc,3), scales[ii]))\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two Stage Logistic-Logistic model for prediction of clinical scales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_important_clinical_features(coef, labels):\n",
    "    from palettable.colorbrewer.qualitative import Pastel1_7\n",
    "    vals = np.asarray(list(coef.values())).mean(axis=0)/np.sum(np.asarray(list(coef.values())).mean(axis=0))\n",
    "    def my_autopct(pct):\n",
    "        return ('%1.1f%%' % pct) if pct > 0 else ''\n",
    "    my_circle=plt.Circle( (0,0), 0.7, color='white')\n",
    "    plt.pie(vals, labels=labels, colors=Pastel1_7.hex_colors, textprops={'fontsize': 7}, autopct=my_autopct)\n",
    "    p=plt.gcf()\n",
    "    p.gca().add_artist(my_circle)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ACR diagnosis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.linear_model import LinearRegression\n",
    "final_y = 'acrcsum'\n",
    "# Keep only ACR data\n",
    "lupus_var = ['acrarthritisc','acrantidsdna', 'acrantismith', 'acrlupneph', 'acranac', final_y]\n",
    "# Variables of interest\n",
    "VOIs_gene = {}; VOIs = {}; ROC_AUC = {}; FPR = {}; TPR = {}; R2 = {}; coef = {}\n",
    "model='LogReg'\n",
    "for ii in range(50):\n",
    "    #try:\n",
    "    VOIs_gene[ii], Y_logreg_pred, y_linear_test, ROC_AUC_, FPR_, TPR_ = run_logreg_bank(corr_matrix, clincal, lupus_var, final_y, split=0.6)\n",
    "    if model == 'LogReg':\n",
    "        # Convert ordinal scale to binary SLE Diagnosis or not\n",
    "        y_linear_test = (y_linear_test>=4)*1\n",
    "    logoutput = pd.DataFrame.from_dict(Y_logreg_pred)\n",
    "    X_train, X_test, y_linear_train, y_linear_test_2 = train_test_split(logoutput, y_linear_test, test_size=0.2)\n",
    "    if model == 'LogReg':\n",
    "        os = SMOTE(random_state=0)\n",
    "        X_train, y_linear_train = os.fit_sample(X_train, y_linear_train)\n",
    "        # Setup logistic regression\n",
    "        logreg = LogisticRegression(max_iter=1000)\n",
    "    if model == 'Linear':\n",
    "        logreg = LinearRegression()\n",
    "    # Find optimal parameters\n",
    "    rfe = RFE(logreg, 5)\n",
    "    rfe = rfe.fit(X_train, y_linear_train)\n",
    "    ind_var = np.asarray(list(X_train.columns))\n",
    "    # Preserve optimate features only from train and test set\n",
    "    X_train = X_train[X_train.columns[X_train.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "    X_test = X_test[X_train.columns[X_train.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "    if model == 'LogReg':\n",
    "        clf = LogisticRegression(max_iter=1000)\n",
    "    if model == 'Linear':\n",
    "        clf = LinearRegression()\n",
    "    clf.fit(X_train, y_linear_train)\n",
    "    print('VOIs: {}'.format(list(ind_var[rfe.ranking_==1])))\n",
    "    VOIs[ii] = list(ind_var[rfe.ranking_==1])\n",
    "    model_pred = clf.predict(X_test)\n",
    "    coef[ii] = clf.coef_.tolist()[0]\n",
    "    labels = list(X_test.columns)\n",
    "    if model == 'LogReg':\n",
    "        # ROC curve\n",
    "        logit_roc_auc = roc_auc_score(y_linear_test_2, clf.predict(X_test))\n",
    "        ROC_AUC[ii] = logit_roc_auc \n",
    "        # Collect True and False Positive rates\n",
    "        fpr, tpr, thresholds = roc_curve(y_linear_test_2, clf.predict_proba(X_test)[:,1])\n",
    "        FPR[ii] = fpr\n",
    "        TPR[ii] = tpr\n",
    "        print('2nd stage ROC AUC: {}'.format(logit_roc_auc))\n",
    "    if model == 'Linear':\n",
    "        R2[ii] = clf.score(X_test, y_linear_test_2)\n",
    "        print('R squared: {}'.format(clf.score(X_test, y_linear_test_2)))\n",
    "#except:\n",
    "    #    continue\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gene signatures underlying diagnosis of SLE with ACR scale\n",
    "VOI_genes = [list(np.ravel(list(VOIs_gene[ii].values()))) for ii in VOIs_gene]\n",
    "VOIs = [item for sublist in VOI_genes for item in sublist]\n",
    "VOIs = {'VOIs':VOIs  }\n",
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR) # split 0.6\n",
    "show_important_clinical_features(coef, labels)\n",
    "# Reach out to Yun Song"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cohort diagnosis (based on SLE or Healthy label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.linear_model import LinearRegression\n",
    "final_y = 'Disease_status'\n",
    "# Keep only ACR data\n",
    "lupus_var = ['acrarthritisc','acrantidsdna', 'acrantismith', 'acrlupneph', 'acranac', final_y]\n",
    "# Variables of interest\n",
    "VOIs_gene = {}; VOIs = {}; ROC_AUC = {}; FPR = {}; TPR = {}; R2 = {}; \n",
    "model='LogReg'\n",
    "for ii in range(50):\n",
    "    #try:\n",
    "    VOIs_gene[ii], Y_logreg_pred, y_linear_test, ROC_AUC_, FPR_, TPR_ = run_logreg_bank(corr_matrix, clincal, lupus_var, final_y, split=0.6)\n",
    "    logoutput = pd.DataFrame.from_dict(Y_logreg_pred)\n",
    "    X_train, X_test, y_linear_train, y_linear_test_2 = train_test_split(logoutput, y_linear_test, test_size=0.2)\n",
    "    if model == 'LogReg':\n",
    "        os = SMOTE(random_state=0)\n",
    "        X_train, y_linear_train = os.fit_sample(X_train, y_linear_train)\n",
    "        # Setup logistic regression\n",
    "        logreg = LogisticRegression(max_iter=1000)\n",
    "    if model == 'Linear':\n",
    "        logreg = LinearRegression()\n",
    "    # Find optimal parameters\n",
    "    rfe = RFE(logreg, 5)\n",
    "    rfe = rfe.fit(X_train, y_linear_train)\n",
    "    ind_var = np.asarray(list(X_train.columns))\n",
    "    # Preserve optimate features only from train and test set\n",
    "    X_train = X_train[X_train.columns[X_train.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "    X_test = X_test[X_train.columns[X_train.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "    if model == 'LogReg':\n",
    "        clf = LogisticRegression(max_iter=1000)\n",
    "    if model == 'Linear':\n",
    "        clf = LinearRegression()\n",
    "    clf.fit(X_train, y_linear_train)\n",
    "    print('VOIs: {}'.format(list(ind_var[rfe.ranking_==1])))\n",
    "    VOIs[ii] = list(ind_var[rfe.ranking_==1])\n",
    "    coef[ii] = clf.coef_.tolist()[0]\n",
    "    labels = list(X_test.columns)\n",
    "    model_pred = clf.predict(X_test)\n",
    "    if model == 'LogReg':\n",
    "        # ROC curve\n",
    "        logit_roc_auc = roc_auc_score(y_linear_test_2, clf.predict(X_test))\n",
    "        ROC_AUC[ii] = logit_roc_auc \n",
    "        # Collect True and False Positive rates\n",
    "        fpr, tpr, thresholds = roc_curve(y_linear_test_2, clf.predict_proba(X_test)[:,1])\n",
    "        FPR[ii] = fpr\n",
    "        TPR[ii] = tpr\n",
    "        print('2nd stage ROC AUC: {}'.format(logit_roc_auc))\n",
    "    if model == 'Linear':\n",
    "        R2[ii] = clf.score(X_test, y_linear_test_2)\n",
    "        print('R squared: {}'.format(clf.score(X_test, y_linear_test_2)))\n",
    "#except:\n",
    "    #    continue\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gene signatures underlying diagnosis of SLE with ACR scale\n",
    "VOI_genes = [list(np.ravel(list(VOIs_gene[ii].values()))) for ii in VOIs_gene]\n",
    "VOIs = [item for sublist in VOI_genes for item in sublist]\n",
    "VOIs = {'VOIs':VOIs}\n",
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR) # split 0.6\n",
    "show_important_clinical_features(coef, labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SLEDAI Diagnosis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.linear_model import LinearRegression\n",
    "final_y = 'acrcsum'\n",
    "# Keep only ACR data\n",
    "lupus_var = ['acrarthritisc','acrantidsdna', 'acrantismith', 'acrlupneph', 'acranac', final_y]\n",
    "# Variables of interest\n",
    "VOIs_gene = {}; VOIs = {}; ROC_AUC = {}; FPR = {}; TPR = {}; R2 = {}; \n",
    "model='LogReg'\n",
    "for ii in range(50):\n",
    "    #try:\n",
    "    VOIs_gene[ii], Y_logreg_pred, y_linear_test, ROC_AUC_, FPR_, TPR_ = run_logreg_bank(corr_matrix, clincal, lupus_var, final_y, split=0.6)\n",
    "    if model == 'LogReg':\n",
    "        # Convert ordinal scale to binary SLE Diagnosis or not\n",
    "        y_linear_test = (y_linear_test>=4)*1\n",
    "    logoutput = pd.DataFrame.from_dict(Y_logreg_pred)\n",
    "    X_train, X_test, y_linear_train, y_linear_test_2 = train_test_split(logoutput, y_linear_test, test_size=0.2)\n",
    "    if model == 'LogReg':\n",
    "        os = SMOTE(random_state=0)\n",
    "        X_train, y_linear_train = os.fit_sample(X_train, y_linear_train)\n",
    "        # Setup logistic regression\n",
    "        logreg = LogisticRegression(max_iter=1000)\n",
    "    if model == 'Linear':\n",
    "        logreg = LinearRegression()\n",
    "    # Find optimal parameters\n",
    "    rfe = RFE(logreg, 5)\n",
    "    rfe = rfe.fit(X_train, y_linear_train)\n",
    "    ind_var = np.asarray(list(X_train.columns))\n",
    "    # Preserve optimate features only from train and test set\n",
    "    X_train = X_train[X_train.columns[X_train.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "    X_test = X_test[X_train.columns[X_train.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "    if model == 'LogReg':\n",
    "        clf = LogisticRegression(max_iter=1000)\n",
    "    if model == 'Linear':\n",
    "        clf = LinearRegression()\n",
    "    clf.fit(X_train, y_linear_train)\n",
    "    print('VOIs: {}'.format(list(ind_var[rfe.ranking_==1])))\n",
    "    VOIs[ii] = list(ind_var[rfe.ranking_==1])\n",
    "    coef[ii] = clf.coef_.tolist()[0]\n",
    "    labels = list(X_test.columns)\n",
    "    model_pred = clf.predict(X_test)\n",
    "    if model == 'LogReg':\n",
    "        # ROC curve\n",
    "        logit_roc_auc = roc_auc_score(y_linear_test_2, clf.predict(X_test))\n",
    "        ROC_AUC[ii] = logit_roc_auc \n",
    "        # Collect True and False Positive rates\n",
    "        fpr, tpr, thresholds = roc_curve(y_linear_test_2, clf.predict_proba(X_test)[:,1])\n",
    "        FPR[ii] = fpr\n",
    "        TPR[ii] = tpr\n",
    "        print('2nd stage ROC AUC: {}'.format(logit_roc_auc))\n",
    "    if model == 'Linear':\n",
    "        R2[ii] = clf.score(X_test, y_linear_test_2)\n",
    "        print('R squared: {}'.format(clf.score(X_test, y_linear_test_2)))\n",
    "#except:\n",
    "    #    continue\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gene signatures underlying diagnosis of SLE with ACR scale\n",
    "VOI_genes = [list(np.ravel(list(VOIs_gene[ii].values()))) for ii in VOIs_gene]\n",
    "VOIs = [item for sublist in VOI_genes for item in sublist]\n",
    "VOIs = {'VOIs':VOIs}\n",
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR) # split 0.6\n",
    "show_important_clinical_features(coef, labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SLAQ Diagnosis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.linear_model import LinearRegression\n",
    "final_y = 'SLAQ score'\n",
    "# Keep only ACR data\n",
    "lupus_var = ['acrarthritisc','acrantidsdna', 'acrantismith', 'acrlupneph', 'acranac', final_y]\n",
    "# Variables of interest\n",
    "VOIs_gene = {}; VOIs = {}; ROC_AUC = {}; FPR = {}; TPR = {}; R2 = {}; \n",
    "model='LogReg'\n",
    "for ii in range(50):\n",
    "    #try:\n",
    "    VOIs_gene[ii], Y_logreg_pred, y_linear_test, ROC_AUC_, FPR_, TPR_ = run_logreg_bank(corr_matrix, clincal, lupus_var, final_y, split=0.6)\n",
    "    if model == 'LogReg':\n",
    "        # Convert ordinal scale to binary SLE Diagnosis or not\n",
    "        y_linear_test = (y_linear_test>=4)*1\n",
    "    logoutput = pd.DataFrame.from_dict(Y_logreg_pred)\n",
    "    X_train, X_test, y_linear_train, y_linear_test_2 = train_test_split(logoutput, y_linear_test, test_size=0.2)\n",
    "    if model == 'LogReg':\n",
    "        os = SMOTE(random_state=0)\n",
    "        X_train, y_linear_train = os.fit_sample(X_train, y_linear_train)\n",
    "        # Setup logistic regression\n",
    "        logreg = LogisticRegression(max_iter=1000)\n",
    "    if model == 'Linear':\n",
    "        logreg = LinearRegression()\n",
    "    # Find optimal parameters\n",
    "    rfe = RFE(logreg, 5)\n",
    "    rfe = rfe.fit(X_train, y_linear_train)\n",
    "    ind_var = np.asarray(list(X_train.columns))\n",
    "    # Preserve optimate features only from train and test set\n",
    "    X_train = X_train[X_train.columns[X_train.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "    X_test = X_test[X_train.columns[X_train.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "    if model == 'LogReg':\n",
    "        clf = LogisticRegression(max_iter=1000)\n",
    "    if model == 'Linear':\n",
    "        clf = LinearRegression()\n",
    "    clf.fit(X_train, y_linear_train)\n",
    "    print('VOIs: {}'.format(list(ind_var[rfe.ranking_==1])))\n",
    "    VOIs[ii] = list(ind_var[rfe.ranking_==1])\n",
    "    coef[ii] = clf.coef_.tolist()[0]\n",
    "    labels = list(X_test.columns)\n",
    "    model_pred = clf.predict(X_test)\n",
    "    if model == 'LogReg':\n",
    "        # ROC curve\n",
    "        logit_roc_auc = roc_auc_score(y_linear_test_2, clf.predict(X_test))\n",
    "        ROC_AUC[ii] = logit_roc_auc \n",
    "        # Collect True and False Positive rates\n",
    "        fpr, tpr, thresholds = roc_curve(y_linear_test_2, clf.predict_proba(X_test)[:,1])\n",
    "        FPR[ii] = fpr\n",
    "        TPR[ii] = tpr\n",
    "        print('2nd stage ROC AUC: {}'.format(logit_roc_auc))\n",
    "    if model == 'Linear':\n",
    "        R2[ii] = clf.score(X_test, y_linear_test_2)\n",
    "        print('R squared: {}'.format(clf.score(X_test, y_linear_test_2)))\n",
    "#except:\n",
    "    #    continue\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gene signatures underlying diagnosis of SLE with ACR scale\n",
    "VOI_genes = [list(np.ravel(list(VOIs_gene[ii].values()))) for ii in VOIs_gene]\n",
    "VOIs = [item for sublist in VOI_genes for item in sublist]\n",
    "VOIs = {'VOIs':VOIs}\n",
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR) # split 0.6\n",
    "show_important_clinical_features(coef, labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SLIC Diagnosis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.linear_model import LinearRegression\n",
    "final_y = 'SLIC score'\n",
    "# Keep only ACR data\n",
    "lupus_var = ['acrarthritisc','acrantidsdna', 'acrantismith', 'acrlupneph', 'acranac', final_y]\n",
    "# Variables of interest\n",
    "VOIs_gene = {}; VOIs = {}; ROC_AUC = {}; FPR = {}; TPR = {}; R2 = {}; \n",
    "model='LogReg'\n",
    "for ii in range(50):\n",
    "    #try:\n",
    "    VOIs_gene[ii], Y_logreg_pred, y_linear_test, ROC_AUC_, FPR_, TPR_ = run_logreg_bank(corr_matrix, clincal, lupus_var, final_y, split=0.6)\n",
    "    if model == 'LogReg':\n",
    "        # Convert ordinal scale to binary SLE Diagnosis or not\n",
    "        y_linear_test = (y_linear_test>=4)*1\n",
    "    logoutput = pd.DataFrame.from_dict(Y_logreg_pred)\n",
    "    X_train, X_test, y_linear_train, y_linear_test_2 = train_test_split(logoutput, y_linear_test, test_size=0.2)\n",
    "    if model == 'LogReg':\n",
    "        os = SMOTE(random_state=0)\n",
    "        X_train, y_linear_train = os.fit_sample(X_train, y_linear_train)\n",
    "        # Setup logistic regression\n",
    "        logreg = LogisticRegression(max_iter=1000)\n",
    "    if model == 'Linear':\n",
    "        logreg = LinearRegression()\n",
    "    # Find optimal parameters\n",
    "    rfe = RFE(logreg, 5)\n",
    "    rfe = rfe.fit(X_train, y_linear_train)\n",
    "    ind_var = np.asarray(list(X_train.columns))\n",
    "    # Preserve optimate features only from train and test set\n",
    "    X_train = X_train[X_train.columns[X_train.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "    X_test = X_test[X_train.columns[X_train.columns.isin(ind_var[rfe.ranking_==1])]]\n",
    "    if model == 'LogReg':\n",
    "        clf = LogisticRegression(max_iter=1000)\n",
    "    if model == 'Linear':\n",
    "        clf = LinearRegression()\n",
    "    clf.fit(X_train, y_linear_train)\n",
    "    print('VOIs: {}'.format(list(ind_var[rfe.ranking_==1])))\n",
    "    VOIs[ii] = list(ind_var[rfe.ranking_==1])\n",
    "    coef[ii] = clf.coef_.tolist()[0]\n",
    "    labels = list(X_test.columns)\n",
    "    model_pred = clf.predict(X_test)\n",
    "    if model == 'LogReg':\n",
    "        # ROC curve\n",
    "        logit_roc_auc = roc_auc_score(y_linear_test_2, clf.predict(X_test))\n",
    "        ROC_AUC[ii] = logit_roc_auc \n",
    "        # Collect True and False Positive rates\n",
    "        fpr, tpr, thresholds = roc_curve(y_linear_test_2, clf.predict_proba(X_test)[:,1])\n",
    "        FPR[ii] = fpr\n",
    "        TPR[ii] = tpr\n",
    "        print('2nd stage ROC AUC: {}'.format(logit_roc_auc))\n",
    "    if model == 'Linear':\n",
    "        R2[ii] = clf.score(X_test, y_linear_test_2)\n",
    "        print('R squared: {}'.format(clf.score(X_test, y_linear_test_2)))\n",
    "#except:\n",
    "    #    continue\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gene signatures underlying diagnosis of SLE with ACR scale\n",
    "VOI_genes = [list(np.ravel(list(VOIs_gene[ii].values()))) for ii in VOIs_gene]\n",
    "VOIs = [item for sublist in VOI_genes for item in sublist]\n",
    "VOIs = {'VOIs':VOIs}\n",
    "show_VOI(VOIs)\n",
    "plot_ROC(ROC_AUC, FPR, TPR) # split 0.6\n",
    "show_important_clinical_features(coef, labels)"
   ]
  },
  {
   "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": [
    "# 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": [
    "## 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": [
    "## Memory vs. Helper vs. Naive\n",
    "sc.pl.umap(adata, color=[\"CCR7\",\"IL7R\",\"S100A4\",\"CD58\", 'FAS', 'IL2RA'],show=True, size=3, save=\".TMemThTNaive.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## T-Reg\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": "markdown",
   "metadata": {},
   "source": [
    "# Protein Log2 FC change analysis "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subset data to contain only lupus cohort 3 +- flare and compute log2 FC change for total PBMC \n",
    "# and each cell type population. Correlate with clinical variables and scales in addition to proportion cell types.\n",
    "# Possibly use log2FC genes and add them into correlation plot without recomputing them.\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); 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['ct_cov'].tolist())\n",
    "for ii in range(len(ct_groups)):\n",
    "    try:\n",
    "        ddata = adata[adata.obs['ct_cov']==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"
   ]
  }
 ],
 "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
}