scvelo-CD4T.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scanpy as sc\n",
    "import besca as bc\n",
    "import os\n",
    "import bbknn\n",
    "import numpy as np\n",
    "import scvelo as scv\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as ss\n",
    "from scipy import stats\n",
    "import itertools\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_dict = {'R': 'coral', 'TF': 'firebrick', 'R_nadj': 'lightskyblue','R_adj': 'royalblue'}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.logging.print_versions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading the h5ad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "WD='/Fullanalysis/'\n",
    "analysis_name='sw_besca2_immune'\n",
    "dcdata=sc.read(WD+'analyzed/'+analysis_name+'/'+analysis_name+'.annotated.filtered.h5ad')\n",
    "figdir='/Fullanalysis/figures/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outpath=WD+'/analyzed/sw_besca2_immune/velocity/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velosubset='CD4TnoTregnoProlif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.set_figure_params()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sc.pl.umap(dcdata, color=['leiden', 'leiden_TNK', 'Sample type'],legend_loc='on data', legend_fontsize=0.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.set_figure_params()\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (12,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(dcdata, color=['celltype3_rc_cells'], size=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(dcdata, color=['celltype3_pub'], size=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.set_figure_params()\n",
    "sc.pl.umap(dcdata,color=['CD8A','CD8B','CD4','CD3D'], color_map='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd4=dcdata[(dcdata.obs['celltype2_rc_cells'].isin(['CD4-positive, alpha-beta T cell']))].copy()\n",
    "cd4=cd4[(cd4.obs['celltype2_pub'].isin(['CD4-positive T cell']))].copy()\n",
    "cd8=dcdata[(dcdata.obs['celltype2_rc_cells'].isin(['CD8-positive, alpha-beta T cell']))].copy()\n",
    "cd8=cd8[(cd8.obs['celltype2_pub'].isin(['CD8-positive T cell']))].copy()\n",
    "myelo=dcdata[(dcdata.obs['celltype1_rc_cells'].isin(['myeloid leukocyte']))].copy()\n",
    "myelo=myelo[(myelo.obs['celltype1'].isin(['myeloid leukocyte']))].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Breakdown of cell types per experiment (sample)\n",
    "bc.pl.celllabel_quant_stackedbar(\n",
    "    cd4, count_variable=\"celltype3_pub\", subset_variable=\"Lesion\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Breakdown of cell types per experiment (sample)\n",
    "bc.pl.celllabel_quant_stackedbar(\n",
    "    cd8, count_variable=\"celltype3_pub\", subset_variable=\"Lesion\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(dcdata.obs['celltype2_merged'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(dcdata.obs['celltype2_pub'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dcdata=dcdata[(dcdata.obs['celltype2_rc_cells'].isin(['CD4-positive, alpha-beta T cell']))].copy()\n",
    "dcdata=dcdata[(dcdata.obs['celltype2_pub'].isin(['CD4-positive T cell']))].copy()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dcdata.obs['leidenf']=dcdata.obs['leiden'].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(dcdata, color=['leiden', 'celltype3_merged','celltype3', 'Sample type'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sc.pl.umap(dcdata, color=['celltype3_merged'],legend_fontsize=8) #14, 34, 3, 31\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(dcdata, color=['celltype2_pub'],legend_fontsize=8) #14, 34, 3, 31"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(dcdata, color=['celltype3_pub'],legend_fontsize=8) #14, 34, 3, 31\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(dcdata, color=['CD8A','CD4','GZMB','CCR7','XCL1','KLRB1','CXCR6','CCR6','SLC4A10','LTK'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(dcdata, color=['celltype2_rc_cells'],legend_fontsize=8) #14, 34, 3, 31\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path1='/looms/'\n",
    "myc=['P11_digest_CD45pos.loom', 'P11_PBMC.loom', 'P12_digest_CD45pos.loom', 'P12_PBMC.loom', \n",
    "     'P13_digest_CD45pos.loom', 'P13_PBMC.loom', 'P29_digest_CD45pos.loom', \n",
    "     'P29_PBMC.loom', 'P33_digest_CD45pos.loom', 'P33_PBMC.loom', 'P34_digest_CD45pos.loom',\n",
    "     'P34_PBMC.loom' ,'P38_digest_CD45pos.loom', \n",
    "     'P38_PBMC.loom', 'P40_digest_CD45pos.loom', 'P40_PBMC.loom' ,'P43_digest_2_CD45pos.loom', \n",
    "     'P43_digest_CD45pos.loom', 'P43_PBMC.loom', 'P63_digest_CD45pos.loom', 'P63_PBMC.loom', \n",
    "     'P64_digest_2_CD45pos.loom', 'P64_digest_CD45pos.loom', 'P67_digest_2_CD45pos.loom', \n",
    "     'P67_digest_CD45pos.loom' ,'P67_PBMC.loom' ,'P68_digest_CD45pos.loom' ,'P68_PBMC.loom', \n",
    "     'P69_digest_CD45pos.loom', 'P72_digest_CD45pos.loom', 'P72_PBMC.loom', 'P77_digest_CD45pos.loom', \n",
    "     'P79_digest_2_CD45pos.loom', 'P79_digest_CD45pos.loom', 'P79_PBMC.loom', '\n",
    "     'P82_digest_CD45pos_and_CD45neg.loom', 'P82_PBMC.loom P83_digest_CD45pos.loom', \n",
    "     'P83_PBMC.loom', 'P86_digest_CD45pos.loom', 'P86_PBMC.loom', 'P87_digest_CD45pos.loom', \n",
    "     'P87_PBMC.loom', 'P91_digest_CD45pos.loom', 'P91_PBMC.loom']\n",
    "m1=myc.copy()\n",
    "myck=myc.copy() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loom=scv.read(WD+'/looms/loom_afterconcat.h5ad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sc.pp.neighbors(dcdata)\n",
    "bbknn.bbknn(dcdata, batch_key='PatientID', use_rep='X_pca')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.diffmap(dcdata)\n",
    "\n",
    "\n",
    "sc.pp.neighbors(dcdata, n_neighbors=15, use_rep='X_diffmap')\n",
    "sc.tl.leiden(dcdata, resolution=0.25)\n",
    "\n",
    "\n",
    "sc.tl.paga(dcdata, groups='leiden')\n",
    "\n",
    "sc.pl.umap(dcdata,color='leiden')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.paga(dcdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "sc.tl.umap(dcdata, init_pos='paga')\n",
    "\n",
    "sc.pl.umap(dcdata,color=['leiden','leidenf'],legend_loc='on data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(dcdata,color=['celltype3'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scv.utils.clean_obs_names(dcdata)\n",
    "#scv.utils.clean_obs_names(loom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata=scv.utils.merge(dcdata, loom).copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#adata.write(WD+\"/looms/CD4Tcells-\"+velosubset+\"loom.h5ad\", compression=\"gzip\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata=sc.read(WD+\"/looms/CD4Tcells-\"+velosubset+\"loom.h5ad\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.set_figure_params(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata,color=['XCL1','CCR7','SELL','MKI67', 'STMN1','ENTPD1','PRF1','GZMB','TIGIT','IL2RA','FOXP3',\n",
    "                        'LAG3','CD38','TOX','CD4','CD8A','CD8B', 'CD3D'], \n",
    "           color_map='viridis',save='-subsettrajectorymap-'+velosubset+'.Markers1.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['celltype3_pub'],legend_fontsize=8, \n",
    "          save='-subsettrajectorymap-'+velosubset+'.celltype3_pub.svg') #14, 34, 3, 31\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['celltype3_pub'],legend_fontsize=8, \n",
    "          save='-subsettrajectorymap-'+velosubset+'.celltype3_pub.pdf') #14, 34, 3, 31\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figdir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['PatientID'],legend_fontsize=8,\n",
    "          save='-subsettrajectorymap-'+velosubset+'.PatientID.svg') #14, 34, 3, 31\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['Sample type'],legend_fontsize=8,\n",
    "          save='-subsettrajectorymap-'+velosubset+'.SampleType.svg') #14, 34, 3, 31\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata,color=['leiden'],legend_loc='on data',\n",
    "          save='-subsettrajectorymap-'+velosubset+'.leiden.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata,color=['RCat'],ncols=2,save='-subsettrajectorymap-'+velosubset+'.RCat.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata,color=['Lesion'],ncols=2,save='-subsettrajectorymap-'+velosubset+'.Lesion.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#velosubset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#adata=scv.read(WD+\"analyzed/sw_besca2_immune/velocity/Velo-DCs-All-PBMCandTIL-dyn.h5ad\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#adata=scv.read(WD+\"CD8Tcells-loom.h5ad\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sc.set_figure_params(6)\n",
    "#sc.pl.umap(adata,color=['PatientID'],ncols=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata,color=['leiden'],ncols=2, legend_loc='on data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.proportions(adata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.proportions(adata, groupby='leiden')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata=adata.copy() #[adata.obs.sample(frac=0.5).index].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#subdata=subdata[subdata.obs.sample(frac=0.5).index].copy()#[adata.obs.sample(frac=0.5).index].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata,color=['leiden'],ncols=2, legend_loc='on data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#subdata=adata.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata,color=['Lesion'],ncols=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs['PatientID']=adata.obs['PatientID'].astype('str')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bc.pl.celllabel_quant_boxplot(\n",
    "    adata,\n",
    "    count_variable=\"leiden\",\n",
    "    subset_variable=\"PatientID\",\n",
    "    condition_identifier=\"RCat\",\n",
    "    plot_percentage=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myperc=bc.tl.count_occurrence_subset_conditions(\n",
    "    adata,\n",
    "    'PatientID',\n",
    "    'Lesion',\n",
    "    count_variable=\"leiden\",\n",
    "    return_percentage=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myperc=bc.tl.count_occurrence_subset_conditions(\n",
    "    adata,\n",
    "    'Lesion',\n",
    "    'Sample type',\n",
    "    count_variable=\"leiden\",\n",
    "    return_percentage=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myperc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figdir+'Percentage-Lesion-'+velosubset+'.all.PatientID.svg'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Breakdown of cell types per experiment (sample)\n",
    "bc.pl.celllabel_quant_stackedbar(\n",
    "    adata, count_variable=\"leiden\", subset_variable=\"Lesion\"\n",
    ")\n",
    "plt.savefig(figdir+'Percentage-Lesion-'+velosubset+'.all.PatientID.svg', format='svg', bbox_inches=\"tight\", dpi=300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Breakdown of cell types per experiment (sample)\n",
    "bc.pl.celllabel_quant_stackedbar(\n",
    "    adata, count_variable=\"leiden\", subset_variable=\"RCat\"\n",
    ")\n",
    "plt.savefig(figdir+'Percentage-RCat-'+velosubset+'.all.PatientID.svg', format='svg', bbox_inches=\"tight\", dpi=300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Breakdown of cell types per experiment (sample)\n",
    "bc.pl.celllabel_quant_stackedbar(\n",
    "    adata, count_variable=\"leiden\", subset_variable=\"PatientID\"\n",
    ")\n",
    "plt.savefig(figdir+'Percentage-Patient-'+velosubset+'.all.PatientID.svg', format='svg', bbox_inches=\"tight\", dpi=300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata,color=['CD8A','CD8B','CD4'],ncols=3, color_map='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#catorder=['8','12','10','6','4','5','3','1']\n",
    "#catorder=corder=['2','5','10','7','13','8']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#adata=adata[adata.obs['leiden'].isin(catorder)].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pp.filter_and_normalize(subdata,min_shared_counts=20)\n",
    "scv.pp.neighbors(subdata)\n",
    "scv.pp.moments(subdata, n_pcs=30, n_neighbors=30)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.recover_dynamics(subdata, n_jobs=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.velocity(subdata,mode='dynamical')\n",
    "scv.tl.velocity_graph(subdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scv.utils.randomized_velocity(subdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.velocity_confidence(subdata)\n",
    "keys = 'velocity_length', 'velocity_confidence'\n",
    "scv.pl.scatter(subdata, c=keys, cmap='coolwarm', perc=[5, 95], save='Velo-len-conf-allCD4T.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sc.set_figure_params(10)\n",
    "#import matplotlib.pyplot as plt\n",
    "#plt.rcParams[\"figure.figsize\"] = (10,8)\n",
    "scv.pl.velocity_embedding_stream(subdata, basis='umap',color=['Sample type'], \n",
    "                                 save='-velo-'+velosubset+'-dyn-Sample_type.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.velocity_embedding_stream(subdata, basis='umap',color=['leiden'],\n",
    "                                 save='-velo-'+velosubset+'-dyn-leiden.svg')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.recover_latent_time(subdata)\n",
    "\n",
    "scv.pl.scatter(subdata, color='latent_time', color_map='gnuplot', size=80)\n",
    "\n",
    "sc.pl.umap(subdata, color=['root_cells','end_points'])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.scatter(subdata, color='velocity_pseudotime', color_map='gnuplot', size=80,\n",
    "              save='-velo-'+velosubset+'-pseudotime.svg')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes = subdata.var['fit_likelihood'].sort_values(ascending=False).index\n",
    "scv.pl.scatter(subdata, basis=top_genes[:15], color='leiden',ncols=5, frameon=False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata.obs.loc[((subdata.obs['root_cells']>0.95) & (subdata.obs['Sample type']  == 'PBMC')),'leiden'].value_counts()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata.obs.loc[((subdata.obs['end_points']>0.9) & (subdata.obs['Sample type']  == 'TIL')),'leiden'].value_counts()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata.uns['iroot'] = np.flatnonzero(((subdata.obs['leiden']  == '0')&(subdata.obs['root_cells']>0.95)&(subdata.obs['Sample type']  == 'PBMC')))[0]\n",
    "subdata.uns['iend'] = np.flatnonzero(((subdata.obs['leiden']  == '5')&(subdata.obs['end_points']>0.9)&(subdata.obs['Sample type']  == 'TIL')))[0]\n",
    "\n",
    "scv.tl.recover_latent_time(subdata, root_key='iroot', end_key='iend')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.scatter(subdata, color=['latent_time','velocity_pseudotime'], color_map='gnuplot', size=80, save='-velo-'+velosubset+'-pseudotime-latent.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.paga(subdata, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime', threshold_root_end_prior=0.8)\n",
    "\n",
    "\n",
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(subdata,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1, )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#subcat='TILonly'\n",
    "subdata.write(WD+'analyzed/'+analysis_name+'/velocity/Velo-'+velosubset+'-dyn.calculatedVelo.h5ad')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Read DE genes ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdde=['CD4T-all-TILonly',\n",
    "     \"CD4T-AllnoTregMAIT-TILonly\",\"CD4T-effector-TILonly\",\"CD4T-naiveCM-TILonly\"]\n",
    "ffigdir='/Fullanalysis/analyzed/sw_besca2_immune/figures/publication/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alldesR=pd.DataFrame(None)\n",
    "alldesNR=pd.DataFrame(None)\n",
    "for mycd in cdde:\n",
    "\n",
    "    cdcNRgenes=pd.read_csv(ffigdir+'DEanalysis_wilcoxon.'+mycd+'.NRhigh.tsv',sep='\\t')\n",
    "    cdcRgenes=pd.read_csv(ffigdir+'DEanalysis_wilcoxon.'+mycd+'.Rhigh.tsv',sep='\\t')\n",
    "\n",
    "    cdcNRgenes['Comparison']=mycd\n",
    "    cdcRgenes['Comparison']=mycd\n",
    "    alldesR=pd.concat([alldesR, cdcRgenes])\n",
    "    alldesNR=pd.concat([alldesNR, cdcNRgenes])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alldesR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alldesNR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run again focusing on main transition clusters ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velosubset='subCD4TnoTregnoProlif_v2'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sc.pl.umap(subdata,color=['leiden'],ncols=2, legend_loc='on data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata=adata.copy() #[adata.obs.sample(frac=0.5).index].copy()\n",
    "\n",
    "###subdata=subdata[subdata.obs['leiden']!='5'].copy()\n",
    "subdata=subdata[subdata.obs['leiden']!='14'].copy()\n",
    "#subdata=subdata[subdata.obs['leiden']!='9'].copy()\n",
    "subdata=subdata[subdata.obs['leiden']!='7'].copy()\n",
    "subdata=subdata[subdata.obs['leiden']!='3'].copy()\n",
    "subdata=subdata[subdata.obs['leiden']!='4'].copy()\n",
    "subdata=subdata[subdata.obs['leiden']!='12'].copy()\n",
    "#subdata=subdata[subdata.obs['leiden']!='10'].copy()\n",
    "\n",
    "#subdata=subdata[subdata.obs.sample(frac=0.5).index].copy()#[adata.obs.sample(frac=0.5).index].copy()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata,color=['leiden'],ncols=2, legend_loc='on data',\n",
    "           save='-subsettrajectorymap-'+velosubset+'.Leiden.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata,color=['celltype3_pub'],ncols=2,save='-subsettrajectorymap-'+velosubset+'.celltype3_pub.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata,color=['CD8A','CD4','FOXP3','MKI67'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### get subsets of the data per response ###\n",
    "#rdata=subdata[subdata.obs['RCat'].isin(['R'])].copy()\n",
    "#tfdata=subdata[subdata.obs['RCat'].isin(['TF'])].copy()\n",
    "#rbdata=subdata[subdata.obs['RCat'].isin(['R', 'TF'])].copy()\n",
    "\n",
    "#nrdata=subdata[subdata.obs['RCat'].isin(['NR_adj'])].copy()\n",
    "#nrnadjdata=subdata[subdata.obs['RCat'].isin(['NR_nadj'])].copy()\n",
    "#nrbdata=subdata[subdata.obs['RCat'].isin(['NR_adj','NR_nadj'])].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pp.filter_and_normalize(subdata,min_shared_counts=20, n_top_genes=800)\n",
    "scv.pp.neighbors(subdata)\n",
    "scv.pp.moments(subdata, n_pcs=30, n_neighbors=30)\n",
    "\n",
    "scv.tl.recover_dynamics(subdata, n_jobs=20)\n",
    "\n",
    "scv.tl.velocity(subdata,mode='dynamical')\n",
    "scv.tl.velocity_graph(subdata)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scv.tl.velocity_graph(subdata, mode_neighbors='connectivities') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#scv.utils.randomized_velocity(subdata)\n",
    "\n",
    "scv.tl.velocity_confidence(subdata)\n",
    "keys = 'velocity_length', 'velocity_confidence'\n",
    "scv.pl.scatter(subdata, c=keys, cmap='coolwarm', perc=[5, 95])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.velocity_embedding_stream(subdata, basis='umap',color=['Sample type'], \n",
    "                                 save='-velo-'+velosubset+'-dyn-Sample_type.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sc.set_figure_params(10)\n",
    "#import matplotlib.pyplot as plt\n",
    "#plt.rcParams[\"figure.figsize\"] = (10,8)\n",
    "scv.pl.velocity_embedding_stream(subdata, basis='umap',color=['Sample type'], \n",
    "                                 save='-velo-'+velosubset+'-dyn-Sample_type.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.velocity_embedding_stream(subdata, basis='umap',color=['leiden'], \n",
    "                                 save='-velo-'+velosubset+'-dyn-leiden.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.velocity_embedding_stream(subdata, basis='umap',color=['leiden'], \n",
    "                                 save='-velo-'+velosubset+'-dyn-leiden.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### No initialisation latent time \n",
    "\n",
    "scv.tl.recover_latent_time(subdata)\n",
    "scv.pl.scatter(subdata, color='latent_time', color_map='gnuplot', size=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata, color=['root_cells','end_points'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.scatter(subdata, color='velocity_pseudotime', color_map='gnuplot', size=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes = subdata.var['fit_likelihood'].sort_values(ascending=False).index\n",
    "scv.pl.scatter(subdata, basis=top_genes[:15], color='leiden',ncols=5, frameon=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata.obs.loc[((subdata.obs['root_cells']>0.38) & (subdata.obs['Sample type']  == 'PBMC')),'leiden'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata.obs.loc[((subdata.obs['end_points']>0.9) & (subdata.obs['Sample type']  == 'TIL')),'leiden'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata.uns['iroot'] = np.flatnonzero(((subdata.obs['leiden']  == '0')&(subdata.obs['root_cells']>0.38)&(subdata.obs['Sample type']  == 'PBMC')))[0]\n",
    "subdata.uns['iend'] = np.flatnonzero(((subdata.obs['leiden']  == '5')&(subdata.obs['end_points']>0.9)&(subdata.obs['Sample type']  == 'TIL')))[0]\n",
    "\n",
    "scv.tl.recover_latent_time(subdata, root_key='iroot', end_key='iend')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "scv.pl.scatter(subdata, color=['latent_time','velocity_pseudotime'], color_map='gnuplot', size=80,\n",
    "               save='-velo-'+velosubset+'-latent_time.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Reanalyse per response subset ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbdata_new=subdata[subdata.obs.RCat.isin(['TF','R'])].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrbdata_new=subdata[subdata.obs.RCat.isin(['NR_adj','NR_nadj'])].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "scv.pp.neighbors(rbdata_new)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pp.moments(rbdata_new, n_pcs=30, n_neighbors=30)\n",
    "\n",
    "scv.tl.recover_dynamics(rbdata_new, n_jobs=20)\n",
    "\n",
    "scv.tl.velocity(rbdata_new,mode='dynamical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.velocity_graph(rbdata_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.paga(rbdata_new, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime', threshold_root_end_prior=0.8)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(rbdata_new,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1,\n",
    "           save='-paga-'+velosubset+'-rbdata-dyn-leiden.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.scatter(rbdata_new, color=['latent_time','velocity_pseudotime'], color_map='gnuplot', size=80,\n",
    "              save='-latent_time-'+velosubset+'-rbdata-dyn-leiden.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.velocity_embedding_stream(rbdata_new, basis='umap',color=['Sample type'], \n",
    "                                save='-velo-'+velosubset+'-rbdata-dyn-sample_type.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbdata_new.obs['clusters']=rbdata_new.obs['leiden'].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes = rbdata_new.var['fit_likelihood'].sort_values(ascending=False).index[:50]\n",
    "scv.pl.heatmap(rbdata_new, var_names=top_genes, sortby='latent_time', col_color='clusters',\n",
    "               n_convolve=100,font_scale=0.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbgenes=rbdata_new.var.loc[rbdata_new.var['velocity_genes']==True,:]['fit_likelihood'].sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrbdata_new.obs['clusters']=nrbdata_new.obs['leiden'].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pp.neighbors(nrbdata_new)\n",
    "\n",
    "scv.pp.moments(nrbdata_new, n_pcs=30, n_neighbors=30)\n",
    "\n",
    "scv.tl.recover_dynamics(nrbdata_new, n_jobs=20)\n",
    "\n",
    "scv.tl.velocity(nrbdata_new,mode='dynamical')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "scv.tl.velocity_graph(nrbdata_new)\n",
    "\n",
    "scv.tl.paga(nrbdata_new, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime', threshold_root_end_prior=0.8)\n",
    "\n",
    "\n",
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(nrbdata_new,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1, \n",
    "           save='-paga-'+velosubset+'-nrbdata-dyn-leiden.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.velocity_embedding_stream(nrbdata_new, basis='umap',color=['Sample type'],\n",
    "                                save='-velo-'+velosubset+'-nrbdata-dyn-sample_type.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrbgenes=nrbdata_new.var.loc[nrbdata_new.var['velocity_genes']==True,:]['fit_likelihood'].sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velogenes=pd.concat([rbgenes, nrbgenes], axis=1)\n",
    "velogenes.columns=['R','NR']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velogenes.sort_values(by='NR', ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbgenes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(rbdata_new,color='leiden')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velogenes.to_csv(WD+'analyzed/'+analysis_name+'/velocity/Velo-genes-per-response-CD4T.tsv', sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes=set(alldesR['Unnamed: 0']).intersection(set(rbgenes[rbgenes>0.1].index))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sc.pl.umap(rbdata_new,color=top_genes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(rbdata_new,color=['ITGB1','BIRC5','CTLA4','ICA1'], ncols=2,\n",
    "           save='-rbdata_new-'+velosubset+'-velogenes.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.heatmap(rbdata_new, var_names=top_genes, sortby='latent_time', col_color='leiden',\n",
    "               n_convolve=100,font_scale=0.6,  save='-rbdata_new-'+velosubset+'.latentimeVeloGenes.Leiden.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes=set(alldesNR['Unnamed: 0']).intersection(set(nrbgenes[nrbgenes>0.1].index))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(nrbdata_new,color=top_genes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(nrbdata_new,color=['NR4A2','NFKBIA','FOSB','ZNF331'], ncols=2,\n",
    "           save='-nrbdata_new-'+velosubset+'-velogenes.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.heatmap(nrbdata_new, var_names=top_genes, sortby='latent_time', col_color='celltype3_pub',\n",
    "               n_convolve=100,font_scale=0.6,save='-nrbdata_new-'+velosubset+'.latentimeVeloGenes.Leiden.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### For the top R patient check velocity "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bc.tl.count_occurrence(rbdata_new, count_variable=\"PatientID\", add_percentage=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bc.tl.count_occurrence(nrbdata_new, count_variable=\"PatientID\", add_percentage=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p34=subdata[subdata.obs.PatientID=='34'].copy()\n",
    "p34=p34[p34.obs['Sample type']=='TIL']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pp.neighbors(p34)\n",
    "\n",
    "\n",
    "scv.pp.moments(p34, n_pcs=30, n_neighbors=30)\n",
    "\n",
    "scv.tl.recover_dynamics(p34, n_jobs=20)\n",
    "\n",
    "scv.tl.velocity(p34,mode='dynamical')\n",
    "\n",
    "scv.tl.velocity_graph(p34)\n",
    "\n",
    "scv.tl.paga(p34, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime', threshold_root_end_prior=0.8)\n",
    "\n",
    "\n",
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(p34,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1,\n",
    "           save='-paga-'+velosubset+'-p34-dyn-leiden.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.velocity_embedding_stream(p34, basis='umap',color=['Sample type'],\n",
    "                                save='-velo-'+velosubset+'-p34-dyn-sample_type.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p67=subdata[subdata.obs.PatientID=='67'].copy()\n",
    "p67=p67[p67.obs['Sample type']=='TIL']\n",
    "\n",
    "scv.pp.neighbors(p67)\n",
    "\n",
    "\n",
    "scv.pp.moments(p67, n_pcs=30, n_neighbors=30)\n",
    "\n",
    "scv.tl.recover_dynamics(p67, n_jobs=20)\n",
    "\n",
    "scv.tl.velocity(p67,mode='dynamical')\n",
    "\n",
    "scv.tl.velocity_graph(p67)\n",
    "\n",
    "scv.tl.paga(p67, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime', threshold_root_end_prior=0.8)\n",
    "\n",
    "\n",
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(p67,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1,\n",
    "           save='-paga-'+velosubset+'-p67-dyn-leiden.svg')\n",
    "\n",
    "scv.pl.velocity_embedding_stream(p67, basis='umap',color=['Sample type'],\n",
    "                                save='-velo-'+velosubset+'-p67-dyn-sample_type.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p79=subdata[subdata.obs.PatientID=='79'].copy()\n",
    "p79=p79[p79.obs['Sample type']=='TIL']\n",
    "\n",
    "scv.pp.neighbors(p79)\n",
    "\n",
    "\n",
    "scv.pp.moments(p79, n_pcs=30, n_neighbors=30)\n",
    "\n",
    "scv.tl.recover_dynamics(p79, n_jobs=20)\n",
    "\n",
    "scv.tl.velocity(p79,mode='dynamical')\n",
    "\n",
    "scv.tl.velocity_graph(p79)\n",
    "\n",
    "scv.tl.paga(p79, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime', threshold_root_end_prior=0.8)\n",
    "\n",
    "\n",
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(p79,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1,\n",
    "           save='-paga-'+velosubset+'-p67-dyn-leiden.svg')\n",
    "\n",
    "scv.pl.velocity_embedding_stream(p79, basis='umap',color=['Sample type'],\n",
    "                                save='-velo-'+velosubset+'-p79-dyn-sample_type.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p64=subdata[subdata.obs.PatientID=='64'].copy()\n",
    "p64=p64[p64.obs['Sample type']=='TIL']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "scv.pp.neighbors(p64)\n",
    "\n",
    "\n",
    "scv.pp.moments(p64, n_pcs=30, n_neighbors=30)\n",
    "\n",
    "scv.tl.recover_dynamics(p64, n_jobs=20)\n",
    "\n",
    "scv.tl.velocity(p64,mode='dynamical')\n",
    "\n",
    "scv.tl.velocity_graph(p64)\n",
    "\n",
    "scv.tl.paga(p64, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime', threshold_root_end_prior=0.8)\n",
    "\n",
    "\n",
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(p64,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1,\n",
    "           save='-paga-'+velosubset+'-p64-dyn-leiden.svg')\n",
    "\n",
    "scv.pl.velocity_embedding_stream(p64, basis='umap',color=['Sample type'],\n",
    "                                save='-velo-'+velosubset+'-p64-dyn-sample_type.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p77=subdata[subdata.obs.PatientID=='77'].copy()\n",
    "p77\n",
    "\n",
    "scv.pp.neighbors(p77)\n",
    "\n",
    "\n",
    "scv.pp.moments(p77, n_pcs=30, n_neighbors=30)\n",
    "\n",
    "scv.tl.recover_dynamics(p77, n_jobs=20)\n",
    "\n",
    "scv.tl.velocity(p77,mode='dynamical')\n",
    "\n",
    "scv.tl.velocity_graph(p77)\n",
    "\n",
    "scv.tl.paga(p77, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime', threshold_root_end_prior=0.8)\n",
    "\n",
    "\n",
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(p77,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1,\n",
    "           save='-paga-'+velosubset+'-p77-dyn-leiden.svg')\n",
    "\n",
    "scv.pl.velocity_embedding_stream(p77, basis='umap',color=['Sample type'],\n",
    "                                save='-velo-'+velosubset+'-p77-dyn-sample_type.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p33=subdata[subdata.obs.PatientID=='33'].copy()\n",
    "p33=p33[p33.obs['Sample type']=='TIL']\n",
    "\n",
    "scv.pp.neighbors(p33)\n",
    "\n",
    "\n",
    "scv.pp.moments(p33, n_pcs=30, n_neighbors=30)\n",
    "\n",
    "scv.tl.recover_dynamics(p33, n_jobs=20)\n",
    "\n",
    "scv.tl.velocity(p33,mode='dynamical')\n",
    "\n",
    "scv.tl.velocity_graph(p33)\n",
    "\n",
    "scv.tl.paga(p33, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime', threshold_root_end_prior=0.8)\n",
    "\n",
    "\n",
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(p33,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1,\n",
    "           save='-paga-'+velosubset+'-p33-dyn-leiden.svg')\n",
    "\n",
    "scv.pl.velocity_embedding_stream(p33, basis='umap',color=['Sample type'],\n",
    "                                save='-velo-'+velosubset+'-p33-dyn-sample_type.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p91=subdata[subdata.obs.PatientID=='91'].copy()\n",
    "p91=p91[p91.obs['Sample type']=='TIL']\n",
    "\n",
    "scv.pp.neighbors(p91)\n",
    "\n",
    "\n",
    "scv.pp.moments(p91, n_pcs=30, n_neighbors=30)\n",
    "\n",
    "scv.tl.recover_dynamics(p91, n_jobs=20)\n",
    "\n",
    "scv.tl.velocity(p91,mode='dynamical')\n",
    "\n",
    "scv.tl.velocity_graph(p91)\n",
    "\n",
    "scv.tl.paga(p91, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime', threshold_root_end_prior=0.8)\n",
    "\n",
    "\n",
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(p91,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1,\n",
    "           save='-paga-'+velosubset+'-p91-dyn-leiden.svg')\n",
    "\n",
    "scv.pl.velocity_embedding_stream(p91, basis='umap',color=['Sample type'],\n",
    "                                save='-velo-'+velosubset+'-p91-dyn-sample_type.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "li={}\n",
    "li['NaiTcell']=['CCR7','IL7R','TCF7']\n",
    "li['MemTcell']=[ 'SELL', 'C1orf21', 'KLRB1', 'ARL4C']\n",
    "li['TExh']=['TNFRSF9', 'CSF1', 'TIGIT']\n",
    "\n",
    "sc.tl.score_genes(adata,li['NaiTcell'],score_name=\"LI_NAI\",use_raw=True)\n",
    "sc.tl.score_genes(adata,li['MemTcell'],score_name=\"LI_MEM\",use_raw=True)\n",
    "sc.tl.score_genes(adata,li['TExh'],score_name=\"LI_EXH\",use_raw=True)\n",
    "\n",
    "wu={}\n",
    "wu['CD4Tcell']=['IL6ST','CRIP1']\n",
    "wu['Teff']=['CX3CR1','GNLY', 'NKG7',  'GZMH', 'KLRD1', 'GZMB', 'PRF1', \n",
    "       'IFITM2', 'LITAF','ITGB2','GZMA','GPR56','KLRC2','GZMM','RAP1B'] ## KLRC2 is higher in Rs\n",
    "wu['TEM']=['GZMK','CCL4',   'DUSP2', 'CD74','DNAJB1','FOS','CCL3','IFNG'] ## DNAJ1, DUSP2, GZMK is higher in NRs\n",
    "wu['TRM']=['CCL4', 'XCL1',   'XCL2',   'ZNF683']\n",
    "wu['IL17']=['NCR3','KLRB1','LYAR','IL7R']\n",
    "\n",
    "sc.tl.score_genes(adata,wu['CD4Tcell'],score_name=\"WU_CD4\",use_raw=True)\n",
    "sc.tl.score_genes(adata,wu['Teff'],score_name=\"WU_TEFF\",use_raw=True)\n",
    "sc.tl.score_genes(adata,wu['TEM'],score_name=\"WU_TEM\",use_raw=True)\n",
    "sc.tl.score_genes(adata,wu['TRM'],score_name=\"WU_TRM\",use_raw=True)\n",
    "sc.tl.score_genes(adata,wu['IL17'],score_name=\"WU_IL17\",use_raw=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "li={}\n",
    "li['NaiTcell']=['CCR7','IL7R','TCF7']\n",
    "li['MemTcell']=[ 'SELL', 'C1orf21', 'KLRB1', 'ARL4C']\n",
    "li['TExh']=['TNFRSF9', 'CSF1', 'TIGIT']\n",
    "\n",
    "sc.tl.score_genes(subdata,li['NaiTcell'],score_name=\"LI_NAI\",use_raw=True)\n",
    "sc.tl.score_genes(subdata,li['MemTcell'],score_name=\"LI_MEM\",use_raw=True)\n",
    "sc.tl.score_genes(subdata,li['TExh'],score_name=\"LI_EXH\",use_raw=True)\n",
    "\n",
    "wu={}\n",
    "wu['CD4Tcell']=['IL6ST','CRIP1']\n",
    "wu['Teff']=['CX3CR1','GNLY', 'NKG7',  'GZMH', 'KLRD1', 'GZMB', 'PRF1', \n",
    "       'IFITM2', 'LITAF','ITGB2','GZMA','GPR56','KLRC2','GZMM','RAP1B'] ## KLRC2 is higher in Rs\n",
    "wu['TEM']=['GZMK','CCL4',   'DUSP2', 'CD74','DNAJB1','FOS','CCL3','IFNG'] ## DNAJ1, DUSP2, GZMK is higher in NRs\n",
    "wu['TRM']=['CCL4', 'XCL1',   'XCL2',   'ZNF683']\n",
    "wu['IL17']=['NCR3','KLRB1','LYAR','IL7R']\n",
    "\n",
    "sc.tl.score_genes(subdata,wu['CD4Tcell'],score_name=\"WU_CD4\",use_raw=True)\n",
    "sc.tl.score_genes(subdata,wu['Teff'],score_name=\"WU_TEFF\",use_raw=True)\n",
    "sc.tl.score_genes(subdata,wu['TEM'],score_name=\"WU_TEM\",use_raw=True)\n",
    "sc.tl.score_genes(subdata,wu['TRM'],score_name=\"WU_TRM\",use_raw=True)\n",
    "sc.tl.score_genes(subdata,wu['IL17'],score_name=\"WU_IL17\",use_raw=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['LI_NAI','LI_MEM','LI_EXH','WU_CD4','WU_CD4','WU_TRM','WU_TEM','WU_TEFF','WU_IL17'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['LI_NAI','LI_MEM','LI_EXH'], save='-markers-'+velosubset+'-Li.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata, color=['LI_NAI','LI_MEM','LI_EXH','WU_CD4','WU_CD4','WU_TRM','WU_TEM','WU_TEFF','WU_IL17'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata, color=['BTLA','TNFRSF1B','HAVCR2','ANXA1'], save='-velo-'+velosubset+'-interactiongenes.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#subcat='TILonly'\n",
    "subdata.write(WD+'analyzed/'+analysis_name+'/velocity/Velo-'+velosubset+'-dyn.calculatedVelo.final.h5ad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata=scv.read(WD+'analyzed/'+analysis_name+'/velocity/Velo-'+velosubset+'-dyn.calculatedVelo.final.h5ad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata.obs['celltype_new']=subdata.obs['leiden'].copy()\n",
    "\n",
    "subdata.obs['celltype_new']=subdata.obs['celltype_new'].replace('1', 'central memory')\n",
    "subdata.obs['celltype_new']=subdata.obs['celltype_new'].replace('0', 'central memory')\n",
    "subdata.obs['celltype_new']=subdata.obs['celltype_new'].replace('11', 'central memory')\n",
    "subdata.obs['celltype_new']=subdata.obs['celltype_new'].replace('2', 'effector memory')\n",
    "subdata.obs['celltype_new']=subdata.obs['celltype_new'].replace('5', 'effector memory, cellular stress')\n",
    "subdata.obs['celltype_new']=subdata.obs['celltype_new'].replace('13', 'activated, interferon')\n",
    "subdata.obs['celltype_new']=subdata.obs['celltype_new'].replace('6', 'effector')\n",
    "subdata.obs['celltype_new']=subdata.obs['celltype_new'].replace('8', 'follicular helper')\n",
    "subdata.obs['celltype_new']=subdata.obs['celltype_new'].replace('9', 'proliferating')\n",
    "subdata.obs['celltype_new']=subdata.obs['celltype_new'].replace('10', 'central memory/Treg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata, color='leiden',legend_loc='on data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata, color='celltype_new')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata, color='celltype3_pub')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "goi=['CDCA2','CDC20','DLGAP5','MKI67','BIRC5','AURKB','TYMS','STMN1',\n",
    "     'CD40LG','TMIGD2','ACTN1','KLF2','PLAC8','TCF7','SELL','CCR7','IL7R',  'PTGER2',\n",
    "      'PIK3R1','KLRG1','ZNF683','FCRL6','MATK','GNLY',\n",
    "     'XCL1','XCL2','CD200','CXCR3','SLAMF6','IL2','IL18R1','CXCR5','ID3','S1PR1','KLF3','IFNG','CCL4','CCL3','CXCL13',\n",
    "     'KLRF1','KLRB1','CRTAM','CD28','CAV1','IFITM1','CXCR6','KLRC2','ITGA1','CD101','CD69','FOS','GZMK','CD27',\n",
    "     'HLA-DQA1','GZMB','GZMA','PRF1','TNFRSF9','PRDM1','TOX','TIGIT','LAG3','HAVCR2','PDCD1',\n",
    "     'ENTPD1','CTLA4','ITGAE','CD38','DUSP4','DUSP1','ATF3','HSPA1A','HSPA1B']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labeloi=\"leiden\"\n",
    "\n",
    "### Perform DE cells of each celltype3 vs. all other cells\n",
    "DEgenes = bc.tl.dge.get_de(\n",
    "    adata, labeloi, demethod=\"wilcoxon\", topnr=5000, logfc=1, padj=0.05\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Export top marker genes per celltype to .tsv file to facilitate checking the annotation\n",
    "top5s=list()\n",
    "for myc in DEgenes.keys():\n",
    "    top5s=top5s+list(DEgenes[myc].sort_values(\"Log2FC\", ascending=False).iloc[0:8,:]['Name'])\n",
    "    DEgenes[myc].sort_values(\"Log2FC\", ascending=False).iloc[0:50,:].to_csv(WD + \"/figures/TopMarkergenes-\" + labeloi + \"_\"+ myc + \"_top50_wilcox.csv\", sep=\"\\t\",index=False)\n",
    "\n",
    "### Generate plots of the top marker genes \n",
    "sc.pl.dotplot(adata, var_names=top5s, groupby=labeloi, dot_max=0.8,vmax=3, dendrogram=True, save=\"TopMarkerplots-\"+labeloi+\".svg\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.matrixplot(adata, var_names=top5s, groupby=labeloi, standard_scale='var',  dendrogram=True, save=\"TopMarkerplots-\"+labeloi+\".svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp=adata[adata.obs.leiden!='14']\n",
    "#tmp=tmp[tmp.obs.leiden!='12']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myorder=['4', '3', '1', '0', '11', '13','2','10', '7', '12',  '6', '8', '5']\n",
    "myorders=[ '1', '0', '11', '13','2','10','9', '8', '6',  '5']\n",
    "\n",
    "goi=['CDCA2','MKI67','STMN1','MAL',\n",
    "     'TCF7','SELL','CCR7','IL7R','S1PR1','KLF3','IFIT1','RSAD2','ISG15','MX1','IFITM1',\n",
    "       'KLRG1','ZNF683','FCRL6',\n",
    "     'XCL1','IFNG','CCL4','CCL3',\n",
    "     'KLRC2','ITGA1','CD69','FOS','GZMK','CCL5','GZMB','GZMA','PRF1','DUSP1','ATF3','HSPA1A','HSPA1B','PRDM1','LAG3','HAVCR2',\n",
    "     'ENTPD1','CTLA4','CXCR6','ITGAE','CD38','PDCD1','CXCL13','TOX','DUSP4','TNFRSF9','TIGIT','LAYN','RTKN2', 'FOXP3', 'FANK1','IL2RA']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "goi=['CDCA2', 'MKI67', 'STMN1','LRRN3', 'CCR7','NPM1', 'SELL', 'IL7R', 'TCF7', 'KLRB1', 'FOXP3', 'RTKN2', 'IL2RA',\n",
    "'CXCR5', 'CXCL13', 'SLAMF6', 'TOX', 'RUNX3', 'PRDM1', 'IFNGR1','KLRG1','ZNF683','FCRL6','CXCR6',\n",
    "'IFIT1', 'IFITM1', 'ISG15', 'XCL1','XCL2', 'CD69', 'FOS','GZMK','CCL5', 'CCL4','GZMH', 'GZMB', 'GZMA','GNLY',\n",
    "     'NKG7', 'PRF1', 'PRDM1', 'IFNG','FANK1', 'RSAD2', 'MX1', 'GZMA', 'IL4I1', 'JUN',  'IGFL2', \n",
    "     'IL21', 'CAV1', 'PDCD1', 'LAG3', 'HAVCR1', 'ENTPD1', 'CTLA4', 'ITGAE','CD38','TNFRSF9','TIGIT','HSPA1A', 'HSPA1B', 'DUSP4', 'ATF3','DNAJB1']\n",
    "\n",
    "#LRRN3, CCR7, NPM1, SELL, PASK, IL7R, KLRB1, TNFSF13B, GZMK, CCL5, CCL4, GZMH, GZMA, GNLY, NGK7, CST7, PRDM1, IFNG, DNAJB1, HSPA1A, HSPA1B, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "goi=[ 'CCR7', 'SELL', 'IL7R', 'TCF7', 'KLRB1',\n",
    "     'FOXP3', 'RTKN2', 'IL2RA', 'PRDM1', 'RUNX3',  'IFNGR1',\n",
    "'IFIT1', 'IFITM1', 'ISG15','CDCA2', 'MKI67', 'STMN1','CXCR5', 'CXCL13', 'SLAMF6', 'TOX', 'PDCD1','HAVCR1','LAG3',  'ENTPD1', 'CTLA4', 'ITGAE','PRF1', \n",
    "     'XCL1', 'CD69', 'FOS', 'GZMB', 'GZMA',  'HSPA1A', 'HSPA1B', 'ATF3']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "root_path = os.getcwd()\n",
    "bescapath_full = os.path.dirname(bc.__file__)\n",
    "bescapath = os.path.split(bescapath_full)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gmt_file_anno= bescapath + '/besca/datasets/genesets/CellNames_scseqCMs6_sigs.gmt'\n",
    "\n",
    "mymarkers = bc.tl.sig.read_GMT_sign(gmt_file_anno,directed=False)\n",
    "mymarkers = bc.tl.sig.filter_siggenes(adata, mymarkers) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd4sigmarkers=mymarkers['CD4Tcell']+mymarkers['CD4Tcell_IL7Rmax']+mymarkers['NaiCD4Tcell']+mymarkers['CMCD4Tcell']+mymarkers['EMCD4Tcell']+mymarkers['RegTcell']+mymarkers['CytotoxCD4Tcell']+mymarkers['ProlifCD4Tcell']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#['CCR7', 'IL7R', 'S1PR1', 'KLF2', 'NKG7', 'IFNGR1', 'RUNX3', 'PRDM1', 'SLAMF6', 'TOX', 'BCL2', 'PSGL1', 'LY6C', 'CXCR6', 'CD27', 'HAVCR2', 'FR4', 'CXCR5', 'CCR7', 'IL7R', 'CD69', 'GZMB', 'GZMK']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd4m=list(set(cd4sigmarkers ).union(set(['CCR7', 'IL7R', 'S1PR1', 'KLF2', 'NKG7', 'IFNGR1', 'RUNX3', 'PRDM1', 'SLAMF6', 'TOX', 'BCL2', 'CXCR6', 'CD27', 'HAVCR2', 'CXCR5', 'CCR7', 'IL7R', 'CD69', 'GZMB', 'GZMK'])))\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#average_obs,fraction_obs=bc.get_means(tmp, 'leiden')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata, color='leiden')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "average_obs,fraction_obs=bc.get_means(tmp, 'leiden')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Plot parameters for publication \n",
    "def set_pub():    \n",
    "    small_size = 10\n",
    "    medium_size = 12\n",
    "    large_size = 14\n",
    "\n",
    "    resolution = 300 #in dpi\n",
    "    plt.rcParams['font.weight'] = 'normal'\n",
    "    #plt.rc('font', **{'family':'sans-serif','sans-serif':['Helvetica']})\n",
    "    plt.rc('axes', titlesize=large_size, titleweight = \"bold\")               # fontsize of the axes title\n",
    "    plt.rc('axes', labelsize=medium_size, labelweight = \"bold\")               # fontsize of the x and y labels\n",
    "    plt.rc('xtick', labelsize=small_size)               # fontsize of the tick labels\n",
    "    plt.rc('ytick', labelsize=small_size)               # fontsize of the tick labels\n",
    "    plt.rc('legend', fontsize=small_size, title_fontsize = medium_size)               # legend fontsize\n",
    "    plt.rc('figure', titlesize=large_size, titleweight = \"bold\")              # fontsize of the figure title\n",
    "    plt.rc('savefig', dpi=resolution)                   # higher res outputs\n",
    "\n",
    "    plt.rcParams['svg.fonttype'] = 'none'\n",
    "\n",
    "\n",
    "set_pub()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=sns.clustermap(average_obs.loc[myorder,cd4m],col_cluster=True, row_cluster=False,figsize=(20,5),\n",
    "                 cmap='viridis',metric='cosine', standard_scale=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "a=sns.clustermap(average_obs.loc[myorders,goi],col_cluster=True, row_cluster=False,figsize=(20,5),\n",
    "                 cmap='viridis',metric='cosine', standard_scale=1)\n",
    "a.savefig(figdir+'Heatmap-goi-'+velosubset+'.leiden.all.standard.pdf')\n",
    "a.savefig(figdir+'Heatmap-goi-'+velosubset+'.leiden.all.standard.svg', format='svg', bbox_inches=\"tight\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.set_figure_params(14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=sns.clustermap(average_obs.loc[myorder,goi],col_cluster=False, row_cluster=False,figsize=(20,5),\n",
    "                 cmap='viridis',metric='cosine', standard_scale=1)\n",
    "a.savefig(figdir+'Heatmap-goi-'+velosubset+'.leiden.all.standard.notclustered.pdf')\n",
    "a.savefig(figdir+'Heatmap-goi-'+velosubset+'.leiden.all.standard.notclustered.svg', format='svg', bbox_inches=\"tight\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=sns.clustermap(average_obs.loc[myorders,goi],col_cluster=False, row_cluster=False,figsize=(20,5),\n",
    "                 cmap='viridis',metric='cosine', standard_scale=1)\n",
    "a.savefig(figdir+'Heatmap-goi-'+velosubset+'.leiden.subdata.standard.notclustered.pdf')\n",
    "a.savefig(figdir+'Heatmap-goi-'+velosubset+'.leiden.subdata.standard.notclustered.svg', format='svg', bbox_inches=\"tight\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp.obs['celltype_new']=tmp.obs['leiden'].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('4', 'naive')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('3', 'naive')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('1', 'central memory')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('0', 'central memory')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('11', 'central memory')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('2', 'effector memory')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('5', 'effector memory, cellular stress')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('13', 'activated, IFN')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('6', 'effector')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('8', 'follicular helper')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('9', 'proliferating')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('10', 'central memory/Treg')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('7', 'Treg')\n",
    "tmp.obs['celltype_new']=tmp.obs['celltype_new'].replace('12', 'Treg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(tmp, color='celltype_new', save='-velo-'+velosubset+'-celltype_new.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=sns.clustermap(average_obs.loc[:,goi],col_cluster=True, row_cluster=True,figsize=(20,5),\n",
    "                 cmap='viridis',metric='cosine', standard_scale=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Generate plots of the top marker genes \n",
    "sc.pl.dotplot(tmp, var_names=goi, groupby=labeloi, dot_max=0.8,vmax=3, dendrogram=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=sns.clustermap(average_obs.loc[myorder,cd4m],col_cluster=True, row_cluster=False,figsize=(18,5),\n",
    "                 cmap='viridis',metric='cosine', standard_scale=1)\n",
    "#a.savefig(figdir+'Heatmap-goi-'+velosubset+'.celltype3_pub.tdata.standard.pdf')\n",
    "#a.savefig(figdir+'Heatmap-goi-'+velosubset+'.celltype3_pub.tdata.standard.svg', format='svg', bbox_inches=\"tight\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#['ZBTB16','IL7R','BCL2','MAL','LEF1','TCF7','S1PR1','KLF2','NELL2','AQP3','GNLY','HMGB2','STMN1','MKI67',\n",
    "#'SAMD3','MATK','IL18RAP','GZMK','GZMA','KLRD1','']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#myorder=['3', '4', '5', '0', '13', '1', '10', '9', '2', '8', '7', '6', '11','12']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.matrixplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#categories_order=myorder,\n",
    "sc.pl.matrixplot(tmp,var_names=goi,groupby='leiden',dendrogram=False, \n",
    "                 standard_scale='var', categories_order=myorder,\n",
    "                 save='Heatmap-goi-'+velosubset+'.leiden.tdata.nc.scanpy_zscore.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#categories_order=myorder,\n",
    "sc.pl.matrixplot(tmp,var_names=cd4m,groupby='leiden',dendrogram=False, \n",
    "                 standard_scale='var', categories_order=myorder,\n",
    "                 save='Heatmap-goi-'+velosubset+'.leiden.tdata.nc.scanpy_zscore.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#categories_order=myorder,\n",
    "#sc.pl.matrixplot(tmp,var_names=goi,groupby='leiden',dendrogram=False, \n",
    "#                 standard_scale='var', categories_order=myorder,\n",
    "#                 save='Heatmap-goi-'+velosubset+'.leiden.tdata.nc.scanpy_zscore.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#categories_order=myorder,\n",
    "sc.pl.matrixplot(tmp,var_names=goi,groupby='celltype3_pub',dendrogram=False, \n",
    "                 standard_scale='var', \n",
    "                 save='Heatmap-goi-'+velosubset+'.celltype3_pub.tdata.nc.scanpy_zscore.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###subdata.uns['iroot'] = np.flatnonzero(((subdata.obs['leiden']  == '3')&(subdata.obs['root_cells']>0.8)&(subdata.obs['Sample type']  == 'PBMC')))[0]\n",
    "#adata.uns['iroot'] =np.flatnonzero(((adata.obs['root_cells']>0.8)&((adata.obs['celltype3']  == 'naive thymus-derived CD8-positive, alpha-beta T cell')|(adata.obs['celltype3_TNK']  == 'NaiCD8Tcell'))))[0]\n",
    "###root_key='iroot'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###scv.tl.recover_latent_time(subdata, root_key='iroot', end_key='terminal_states_probs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###scv.pl.scatter(subdata, color='latent_time', color_map='gnuplot', size=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scv.tl.recover_latent_time(subdata, root_key='initial_states_probs', end_key='terminal_states_probs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.paga(subdata, groups='leiden', root_key='iroot', end_key='iend') #use_time_prior='velocity_pseudotime'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import cellrank as cr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(subdata,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scv.pl.velocity_graph(subdata, threshold=.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scv.tl.velocity_pseudotime(subdata)\n",
    "#scv.pl.scatter(subdata, color='velocity_pseudotime', cmap='gnuplot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata.obs['clusters']=subdata.obs.leiden.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes = subdata.var['fit_likelihood'].sort_values(ascending=False).index[:50]\n",
    "scv.pl.heatmap(subdata, var_names=top_genes, sortby='latent_time', col_color='clusters',\n",
    "               n_convolve=100,font_scale=0.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes = subdata.var['fit_likelihood'].sort_values(ascending=False).index[:50]\n",
    "scv.pl.heatmap(subdata, var_names=top_genes, sortby='latent_time', col_color='celltype3_pub',\n",
    "               n_convolve=100,font_scale=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.rank_dynamical_genes(subdata, groupby='clusters')\n",
    "df = scv.get_df(subdata, 'rank_dynamical_genes/names')\n",
    "df.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pp.filter_and_normalize(rbdata,min_shared_counts=20, n_top_genes=800)\n",
    "scv.pp.neighbors(rbdata)\n",
    "scv.pp.moments(rbdata)\n",
    "\n",
    "scv.tl.recover_dynamics(rbdata, n_jobs=20)\n",
    "scv.tl.velocity(rbdata,mode='dynamical')\n",
    "scv.tl.velocity_graph(rbdata)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.velocity_graph(rbdata) #mode_neighbors='connectivities' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.velocity_confidence(rbdata)\n",
    "keys = 'velocity_length', 'velocity_confidence'\n",
    "scv.pl.scatter(rbdata, c=keys, cmap='coolwarm', perc=[5, 95])\n",
    "\n",
    "scv.pl.velocity_embedding_stream(rbdata, basis='umap',color=['Sample type'], \n",
    "                                 save='-velo-'+velosubset+'-dyn-Sample_type.rbdata.svg')\n",
    "\n",
    "scv.pl.velocity_embedding_stream(rbdata, basis='umap',color=['leiden'], \n",
    "                                 save='-velo-'+velosubset+'-dyn-leiden.rbdata.svg')\n",
    "\n",
    "scv.tl.recover_latent_time(rbdata)\n",
    "\n",
    "scv.pl.scatter(rbdata, color='latent_time', color_map='gnuplot', size=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.scatter(rbdata, color='velocity_pseudotime', color_map='gnuplot', size=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbdata.obs.loc[((rbdata.obs['root_cells']>0.5) & (rbdata.obs['Sample type']  == 'PBMC')),'leiden'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbdata.obs.loc[((rbdata.obs['end_points']>0.9) & (rbdata.obs['Sample type']  == 'TIL')),'leiden'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbdata.uns['iroot'] = np.flatnonzero(((rbdata.obs['leiden']  == '13')&(rbdata.obs['root_cells']>0.5)&(rbdata.obs['Sample type']  == 'PBMC')))[0]\n",
    "rbdata.uns['iend'] = np.flatnonzero(((rbdata.obs['leiden']  == '1')&(rbdata.obs['end_points']>0.9)&(rbdata.obs['Sample type']  == 'TIL')))[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "scv.tl.recover_latent_time(rbdata, root_key='iroot', end_key='iend')\n",
    "\n",
    "scv.pl.scatter(rbdata, color='latent_time', color_map='gnuplot', size=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(rbdata, color='velocity_pseudotime')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.paga(rbdata, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime', threshold_root_end_prior=0.8)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(rbdata,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(rbdata, color=['root_cells','end_points'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "scv.pp.filter_and_normalize(nrbdata,min_shared_counts=20, n_top_genes=1000)\n",
    "scv.pp.neighbors(nrbdata)\n",
    "scv.pp.moments(nrbdata)\n",
    "\n",
    "scv.tl.recover_dynamics(nrbdata, n_jobs=20)\n",
    "scv.tl.velocity(nrbdata,mode='dynamical')\n",
    "#scv.tl.velocity_graph(nrbdata)\n",
    "scv.tl.velocity_graph(nrbdata, mode_neighbors='connectivities') \n",
    "\n",
    "\n",
    "scv.tl.velocity_confidence(nrbdata)\n",
    "keys = 'velocity_length', 'velocity_confidence'\n",
    "scv.pl.scatter(nrbdata, c=keys, cmap='coolwarm', perc=[5, 95])\n",
    "\n",
    "scv.pl.velocity_embedding_stream(nrbdata, basis='umap',color=['Sample type'], \n",
    "                                 save='-velo-'+velosubset+'-dyn-Sample_type.nrbdata.svg')\n",
    "\n",
    "scv.pl.velocity_embedding_stream(nrbdata, basis='umap',color=['leiden'], \n",
    "                                 save='-velo-'+velosubset+'-dyn-leiden.nrbdata.svg')\n",
    "\n",
    "scv.tl.recover_latent_time(nrbdata)\n",
    "\n",
    "scv.pl.scatter(nrbdata, color='latent_time', color_map='gnuplot', size=80)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrbdata.obs.loc[((nrbdata.obs['root_cells']>0.9) & (nrbdata.obs['Sample type']  == 'PBMC')),'leiden'].value_counts()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "nrbdata.obs.loc[((nrbdata.obs['end_points']>0.9) & (nrbdata.obs['Sample type']  == 'TIL')),'leiden'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrbdata.uns['iroot'] = np.flatnonzero(((nrbdata.obs['leiden']  == '0')&(nrbdata.obs['root_cells']>0.9)&(nrbdata.obs['Sample type']  == 'PBMC')))[0]\n",
    "nrbdata.uns['iend'] = np.flatnonzero(((nrbdata.obs['leiden']  == '5')&(nrbdata.obs['end_points']>0.9)&(nrbdata.obs['Sample type']  == 'TIL')))[0]\n",
    "\n",
    "scv.tl.recover_latent_time(nrbdata, root_key='iroot', end_key='iend')\n",
    "\n",
    "scv.pl.scatter(nrbdata, color='latent_time', color_map='gnuplot', size=80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "scv.tl.paga(nrbdata, groups='leiden', root_key='iroot', end_key='iend', \n",
    "            use_time_prior='velocity_pseudotime')\n",
    "\n",
    "\n",
    "sc.set_figure_params(10)\n",
    "scv.pl.paga(nrbdata,cluster_key=\"leiden\", basis='umap',node_size_scale=1, edge_width_scale=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(nrbdata, color=['root_cells','end_points'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(nrbdata, color=['latent_time','velocity_pseudotime'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#run cellrankanalysis\n",
    "cellrankanalysis=False\n",
    "if cellrankanalysis:\n",
    "    g = cr.estimators.GPCCA(vk)\n",
    "    g.fit(cluster_key=\"clusters\", n_states=[2, 6])\n",
    "    g.plot_macrostates(which=\"all\", discrete=True, legend_loc=\"right\", s=100)\n",
    "\n",
    "    g.predict_terminal_states()\n",
    "    g.plot_macrostates(which=\"terminal\", legend_loc=\"right\", s=100)\n",
    "\n",
    "    g.plot_macrostates(which=\"terminal\", discrete=False)\n",
    "\n",
    "    g.predict_initial_states()\n",
    "    g.plot_macrostates(which=\"initial\", legend_loc=\"right\", s=100)\n",
    "\n",
    "\n",
    "    g.fit(cluster_key=\"clusters\", n_states=11)\n",
    "\n",
    "    g.set_terminal_states(states=[\"Alpha\", \"Beta\", \"Epsilon\", \"Delta\"])\n",
    "    g.plot_macrostates(which=\"terminal\", legend_loc=\"right\", size=100)\n",
    "\n",
    "\n",
    "    g.compute_fate_probabilities()\n",
    "    g.plot_fate_probabilities(same_plot=False)\n",
    "    g.plot_fate_probabilities(same_plot=True)\n",
    "\n",
    "    driver_clusters = [\"Delta\", \"Fev+\", \"Ngn3 high EP\", \"Ngn3 low EP\"]\n",
    "\n",
    "    delta_df = g.compute_lineage_drivers(\n",
    "        lineages=[\"Delta\"], cluster_key=\"clusters\", clusters=driver_clusters\n",
    "    )\n",
    "    delta_df.head(10)\n",
    "\n",
    "    sc.pl.violin(adata, keys=[\"palantir_pseudotime\"], groupby=\"clusters\", rotation=90)\n",
    "\n",
    "    adata.obs[\"fate_probabilities_delta\"] = g.fate_probabilities[\"Delta\"].X.flatten()\n",
    "\n",
    "    sc.pl.embedding(\n",
    "        adata,\n",
    "        basis=\"umap\",\n",
    "        color=[\"fate_probabilities_delta\"] + list(delta_df.index[:8]),\n",
    "        color_map=\"viridis\",\n",
    "        s=50,\n",
    "        ncols=3,\n",
    "        vmax=\"p96\",\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata, color=['root_cells','end_points'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata.obs['clusters']=subdata.obs['leiden'].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes = subdata.var['fit_likelihood'].sort_values(ascending=False).index[:20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.set_figure_params()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "degenes=['PRF1', 'TMEM258', 'SPP1', 'MIF', 'HLA-DRB5', 'ATP5I', 'GZMB', 'PPDPF', 'PDCD1', 'MT2A', \n",
    "          'TUBA1A', 'PTGER4', 'JUNB', 'SERTAD1', 'TNFAIP3', 'PNRC1']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata, color=degenes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata, color=top_genes[0:12])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(subdata, color=top_genes[11:20])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp=subdata[subdata.obs.leiden.isin(['2','10','8','7'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(tmp, color='PatientID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(tmp, color='RCat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp.obs.PatientID='P'+tmp.obs.PatientID.astype(str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clustercomp=bc.tl.count_occurrence_subset(tmp, 'PatientID', count_variable ='leiden', return_percentage = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.set_figure_params()\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (12,5)\n",
    "import matplotlib \n",
    "\n",
    "#clustercomp=clustercomp[clustercomp.max().sort_values().index].copy()\n",
    "#clustercomp.transpose().plot(kind='bar', stacked=True)\n",
    "clustercomp[clustercomp.max().sort_values().index].transpose().plot.bar(stacked=True, figsize=(12, 8))\n",
    "plt.ylabel('Percentage donor per cluster')\n",
    "plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "#fig.savefig(figdir+'Clustercomp-'+velosubset+'.PatientID.pdf'),bbox_inches='tight')\n",
    "plt.savefig(figdir+'Clustercomp-'+velosubset+'.PatientID.svg', format='svg', bbox_inches=\"tight\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes = subdata.var['fit_likelihood'].sort_values(ascending=False).index[:30]\n",
    "scv.pl.heatmap(subdata, var_names=top_genes, sortby='latent_time', col_color='clusters',\n",
    "               font_scale=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scv.pl.velocity_embedding_stream(subdata, vkey=[ 'velocity_random'],color=['celltype3_pub'],legend_loc='right')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.tl.score_genes_cell_cycle(adata)\n",
    "scv.pl.scatter(adata, color_gradients=['S_score', 'G2M_score'], smooth=True, perc=[5, 95])\n",
    "s_genes, g2m_genes = scv.utils.get_phase_marker_genes(adata)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cellrank as cr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "warnings.simplefilter(\"ignore\", category=UserWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vk = cr.kernels.VelocityKernel(subdata)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vk.compute_transition_matrix()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ck = cr.kernels.ConnectivityKernel(subdata)\n",
    "ck.compute_transition_matrix()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ck.plot_projection()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vk.plot_projection()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vk.plot_random_walks(start_ixs={\"leiden\": \"3\"}, max_iter=200, seed=0,    dpi=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ck.plot_random_walks(start_ixs={\"leiden\": \"3\"}, max_iter=200, seed=0,    dpi=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vk.plot_projection(basis=\"umap\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes = subdata.var['fit_likelihood'].sort_values(ascending=False).index[:30]\n",
    "scv.pl.heatmap(subdata, var_names=top_genes, sortby='latent_time', col_color='clusters',\n",
    "               font_scale=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes = subdata.var[\"fit_likelihood\"].sort_values(ascending=False).index\n",
    "scv.pl.scatter(subdata, basis=top_genes[:10], ncols=5, frameon=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vk.write_to_adata()\n",
    "subdata.write(\n",
    "    \"Velo-\"+velosubset+\"_velocity_kernel.h5ad\", compression=\"gzip\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mysubset='All'\n",
    "subcat='PBMCandTIL'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#subcat='TILonly'\n",
    "subdata.write(WD+'analyzed/'+analysis_name+'/velocity/Velo-'+velosubset+'-'+mysubset+'-'+subcat+'-dyn.h5ad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdata=sc.read(WD+'analyzed/'+analysis_name+'/velocity/Velo-'+velosubset+'-'+mysubset+'-'+subcat+'-dyn.h5ad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddata=subdata.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddata.obs['RCat']=ddata.obs['Response3'].copy()\n",
    "ddata.obs['RCat']=list(ddata.obs['RCat'].replace('NR', 'NR_nadj'))\n",
    "ddata.obs['RCat']=list(ddata.obs['RCat'].replace('PD', 'NR_adj'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(ddata, color='leiden',legend_loc='on data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(ddata, color='celltype3_pub')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(ddata.obs['PatientID'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Remove two uncertain patients\n",
    "\n",
    "#ddata=ddata[ddata.obs['PatientID']!='P38'].copy()\n",
    "#ddata=ddata[ddata.obs['PatientID']!='P12'].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddata=ddata[ddata.obs['Sample type']=='TIL'].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ddata=ddata[ddata.obs['celltype3_pub']!='naive CD4-positive T cell'].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set(ddata.obs.celltype3_pub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddata.obs['PatientID']=ddata.obs['PatientID'].astype(str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddata.obs['R_clus']=ddata.obs['RCat'].astype(str)+'-'+[x.replace(\"-\", \"_\") for x in ddata.obs['leiden'].astype(str)]\n",
    "ddata.obs['Pat_clus']=ddata.obs['PatientID'].astype(str)+'-'+[x.replace(\"-\", \"_\") for x in ddata.obs['leiden'].astype(str)]\n",
    "dfPat = ddata.obs.groupby('Pat_clus')['velocity_length', 'velocity_confidence','latent_time'].mean().T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfPat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Resp=[ddata[ddata.obs['PatientID']==x.split('-')[0]].obs['RCat'][0] for x in list(dfPat.columns)]\n",
    "PatientID=[x.split('-')[0]  for x in list(dfPat.columns)]\n",
    "ClusterID=[x.split('-')[1]  for x in list(dfPat.columns)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mydf=pd.DataFrame([list(dfPat.loc['velocity_length',:]),Resp,PatientID,ClusterID]).transpose()\n",
    "mydf.index=dfPat.columns\n",
    "mydf.columns=['velocity_length','Response','PatientID','Cluster']\n",
    "\n",
    "mydfl=pd.DataFrame([list(dfPat.loc['latent_time',:]),Resp,PatientID,ClusterID]).transpose()\n",
    "mydfl.index=dfPat.columns\n",
    "mydfl.columns=['latent_time','Response','PatientID','Cluster']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mydf.groupby('PatientID', as_index=False)['velocity_length']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_dict = {'R': 'coral', 'TF': 'firebrick', 'NR_nadj': 'lightskyblue','NR_adj': 'royalblue'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Velocity test first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#corder=['8','6','1','21']\n",
    "mydf['velocity_length']=mydf['velocity_length'].astype('float')\n",
    "velMean=mydf.groupby('PatientID')['velocity_length'].mean()\n",
    "patid=[ddata[ddata.obs['PatientID']==x.split('-')[0]].obs['RCat'][0]  for x in list(velMean.index)]\n",
    "totest=pd.DataFrame([list(velMean),patid]).transpose()\n",
    "totest.columns=['Val','Cat']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ss.mannwhitneyu(list(totest.loc[totest['Cat']=='R',:]['Val']),list(totest.loc[totest['Cat']=='NR_nadj',:]['Val']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ss.mannwhitneyu(list(totest.loc[totest['Cat']=='TF',:]['Val']),list(totest.loc[totest['Cat']=='NR_adj',:]['Val']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(mydf['Response'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "totestk=totest.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mysubset='Final'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subcat='CD4T'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=sns.boxplot(y='Val',x='Cat',data=totestk,palette=color_dict,order=['R','TF','NR_nadj','NR_adj'])\n",
    "fig=sns.swarmplot(x='Cat',y='Val',data=totestk,color='black',order=['R','TF','NR_nadj','NR_adj'])\n",
    "fig.figure.savefig(figdir+'Velocitylen-'+velosubset+'-'+mysubset+'-'+subcat+'-per-response.pdf') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "totest=list(itertools.combinations(['R','TF','NR_nadj','NR_adj'], 2))\n",
    "pwilc={}\n",
    "pt={}\n",
    "for pairs in totest:\n",
    "    pwilc[pairs[0]+'-'+pairs[1]]=stats.mannwhitneyu(list(totestk.loc[totestk['Cat']==pairs[0],:]['Val']), \n",
    "                           list(totestk.loc[totestk['Cat']==pairs[1],:]['Val']))[1]\n",
    "    pt[pairs[0]+'-'+pairs[1]]=stats.ttest_ind(list(totestk.loc[totestk['Cat']==pairs[0],:]['Val']), \n",
    "                           list(totestk.loc[totestk['Cat']==pairs[1],:]['Val']))[1]\n",
    "\n",
    "myps=pd.DataFrame([pwilc,pt]).transpose()\n",
    "myps.columns=['MannWhitney','T-test']\n",
    "\n",
    "myps.to_csv(figdir+'Pval-Velocitylen-'+velosubset+'-'+mysubset+'-'+subcat+'-per-response.pdf',sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corder=list(set(mydf['Cluster']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in corder:\n",
    "    sns.boxplot(y='velocity_length',x='Response',data=mydf.loc[mydf['Cluster']==i,:],\n",
    "                palette=color_dict,order=['R','TF','NR_nadj','NR_adj'])\n",
    "    sns.swarmplot(x='Response',y='velocity_length',data=mydf.loc[mydf['Cluster']==i,:],\n",
    "                  color='black',order=['R','TF','NR_nadj','NR_adj'])\n",
    "    print('Cluster '+i)\n",
    "    plt.figure()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Latent time next"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mydfl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#corder=['8','6','1','21']\n",
    "mydfl['latent_time']=mydfl['latent_time'].astype('float')\n",
    "velMean=mydfl.groupby('PatientID')['latent_time'].mean()\n",
    "patid=[ddata[ddata.obs['PatientID']==x.split('-')[0]].obs['RCat'][0]  for x in list(velMean.index)]\n",
    "totest=pd.DataFrame([list(velMean),patid]).transpose()\n",
    "totest.columns=['Val','Cat']\n",
    "\n",
    "ss.mannwhitneyu(list(totest.loc[totest['Cat']=='R',:]['Val']),list(totest.loc[totest['Cat']=='NR_nadj',:]['Val']))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ss.mannwhitneyu(list(totest.loc[totest['Cat']=='TF',:]['Val']),list(totest.loc[totest['Cat']=='NR_adj',:]['Val']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "totestk=totest.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=sns.boxplot(y='Val',x='Cat',data=totestk,palette=color_dict,order=['R','TF','NR_nadj','NR_adj'])\n",
    "fig=sns.swarmplot(x='Cat',y='Val',data=totestk,color='black',order=['R','TF','NR_nadj','NR_adj'])\n",
    "fig.figure.savefig(figdir+'LatentTime-'+velosubset+'-'+mysubset+'-'+subcat+'-per-response.pdf') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corder=list(set(mydf['Cluster']))\n",
    "\n",
    "for i in corder:\n",
    "    sns.boxplot(y='latent_time',x='Response',data=mydfl.loc[mydf['Cluster']==i,:],\n",
    "                palette=color_dict,order=['R','TF','NR_nadj','NR_adj'])\n",
    "    sns.swarmplot(x='Response',y='latent_time',data=mydfl.loc[mydfl['Cluster']==i,:],\n",
    "                  color='black',order=['R','TF','NR_nadj','NR_adj'])\n",
    "    plt.title('Cluster ' +i)\n",
    "    print('Cluster '+i)\n",
    "    plt.figure()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velogenes=subdata.var.loc[(subdata.var['velocity_genes']==True)&(subdata.var['fit_likelihood']>.1),:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velogenes.sort_values('fit_likelihood', ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(list(velogenes.index))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_genes = subdata.var['fit_likelihood'].sort_values(ascending=False).index#[:300]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "invelo=[]\n",
    "for x in list(top_genes):\n",
    "    invelo=invelo+list(set([x]).intersection(set(velogenes.index)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.heatmap(subdata, var_names=invelo, sortby='latent_time', \n",
    "               col_color='leiden', n_convolve=30,\n",
    "               font_scale=0.5,  figsize=(8, 6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.heatmap(subdata, var_names=invelo[0:50], \n",
    "               sortby='latent_time', col_color='leiden', n_convolve=30,\n",
    "               font_scale=0.5,  figsize=(8, 6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scv.pl.heatmap(subdata, var_names=invelo[50:100], \n",
    "               sortby='latent_time', col_color='leiden', n_convolve=30,\n",
    "               font_scale=0.5,  figsize=(8, 6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (len(invelo)>100):\n",
    "    scv.pl.heatmap(subdata, var_names=invelo[100:len(invelo)], \n",
    "               sortby='latent_time', col_color='leiden', n_convolve=30,\n",
    "               font_scale=0.5,  figsize=(8, 6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velogenes.sort_values('fit_likelihood', ascending=False).to_csv(figdir+'/VelocityInfoPerGene_'+velosubset+'-'+mysubset+'_'+subcat+'_loom-dyn_allVeloGenesRanked.tsv',sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tkey='latent_time'\n",
    "var_names=invelo\n",
    "var_names = [name for name in var_names if name in adata.var_names]\n",
    "xkey='Ms'\n",
    "time = subdata.obs[tkey].values\n",
    "time = time[np.isfinite(time)]\n",
    "\n",
    "df = pd.DataFrame(subdata[:, var_names].layers[xkey][np.argsort(time)], columns=var_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_sort = np.argsort(np.argmax(df.values, axis=0))\n",
    "df = pd.DataFrame(df.values[:, max_sort], columns=df.columns[max_sort])\n",
    "df.index=pd.DataFrame([subdata.obs['celltype3_merged'][np.argsort(time)],\n",
    "                       subdata.obs['leiden'][np.argsort(time)],subdata.obs['latent_time'][np.argsort(time)]]).transpose().index\n",
    "df=pd.concat([df,pd.DataFrame([subdata.obs['celltype3_merged'][np.argsort(time)],\n",
    "                               subdata.obs['leiden'][np.argsort(time)],subdata.obs['latent_time'][np.argsort(time)]]).transpose()], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getMeans(adata,myvar):\n",
    "    gene_ids = adata.raw.var.index.values\n",
    "    clusters = adata.obs[myvar].cat.categories\n",
    "    obs = adata.raw[:,gene_ids].X.toarray()\n",
    "    obs = pd.DataFrame(obs,columns=gene_ids,index=adata.obs[myvar])\n",
    "    average_obs = obs.groupby(level=0).mean()\n",
    "    obs_bool = obs.astype(bool)\n",
    "    fraction_obs = obs_bool.groupby(level=0).sum()/obs_bool.groupby(level=0).count()\n",
    "    return(average_obs,fraction_obs)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avg,fct=getMeans(adata,'leiden')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! jupyter nbconvert --to html scvelo-CD4T.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "besca_dev_scvelo",
   "language": "python",
   "name": "besca_dev_scvelo"
  },
  "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.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}