Merging_ICI_T_cells.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merged T cell from 3 ICI studies\n",
    "EGAS00001004809, GSE179994, and GSE169246\n",
    "``` \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import scanpy as sc\n",
    "import scipy.io as sio\n",
    "import anndata as ad\n",
    "import os as os\n",
    "import seaborn as sns\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in T cells of EGAS100010040809 \n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in count matrix first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_EGAS100010040809 = sc.read_h5ad('/data/ICI_exprs/ICI_T_cell_collection/EGAS00001004809-1863-counts_cells_cohort1_T_cells.h5ad')\n",
    "print(adata_EGAS100010040809.shape)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in Meta data of each dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EGAS100010040809_meta = pd.read_csv('/data/ICI_exprs/EGAS00001004809/1870-BIOKEY_metaData_tcells_cohort1_web.csv', index_col=0, header=0)\n",
    "print(EGAS100010040809_meta.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loaded matrix and meta data has discrepancy in cell number, \n",
    "so we need to filter out the cells that are not in the meta data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find common cells from adata and meta\n",
    "common_cells = np.intersect1d(adata_EGAS100010040809.obs.index, EGAS100010040809_meta.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " # Add sample_id to EGAS100010040809_meta by parsing the sample name\n",
    "EGAS100010040809_meta['sample_id'] = ['_'.join(x.split('_')[:3]) for x in EGAS100010040809_meta.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EGAS100010040809_to_keep = EGAS100010040809_meta.loc[common_cells,['patient_id', 'timepoint','sample_id']]\n",
    "EGAS100010040809_to_keep['batch'] = ['EGAS100010040809']*len(EGAS100010040809_to_keep)\n",
    "EGAS100010040809_to_keep.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EGAS100010040809_adata = adata_EGAS100010040809[common_cells, :]\n",
    "EGAS100010040809_adata.obs = EGAS100010040809_to_keep.loc[common_cells, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in GSE179994_all.Tcell.rawCounts.h5ad\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GSE179994_adata = sc.read_h5ad('/data/ICI_exprs/GSE179994/GSE179994_all.Tcell.rawCounts.h5ad')\n",
    "print(GSE179994_adata.shape)\n",
    "# the data is a matrix only without any obs or var"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Process GSE179994_meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GSE179994_meta = pd.read_csv('/data/ICI_exprs/GSE179994/GSE179994_Tcell.metadata.tsv.gz', sep = '\\t', index_col = 0, compression = 'gzip', header = 0)\n",
    "print(GSE179994_meta.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace 'sample' with 'sample_id'\n",
    "GSE179994_meta.columns = GSE179994_meta.columns.str.replace('sample', 'sample_id')\n",
    "GSE179994_meta.columns = GSE179994_meta.columns.str.replace('patient', 'patient_id')\n",
    "print(GSE179994_meta.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# split sample_id to extract time point\n",
    "GSE179994_meta['timepoint'] = GSE179994_meta['sample_id'].str.split('.').str[1]\n",
    "# replace 'post' with 'on'\n",
    "GSE179994_meta['timepoint'] = GSE179994_meta['timepoint'].str.replace('pre', 'pre')\n",
    "GSE179994_meta['timepoint'] = GSE179994_meta['timepoint'].str.replace('post', 'on')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GSE179994_meta_to_keep = GSE179994_meta.loc[: ,['patient_id', 'timepoint','sample_id']]\n",
    "GSE179994_meta_to_keep['batch'] = ['GSE179994'] * GSE179994_meta_to_keep.shape[0]\n",
    "#align the index of GSE179994_meta_to_keep with GSE179994_adata\n",
    "GSE179994_meta_to_keep = GSE179994_meta_to_keep.loc[GSE179994_adata.obs_names, :]\n",
    "print(GSE179994_meta_to_keep.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GSE179994_adata.obs = GSE179994_meta_to_keep"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find common genes in Xs of both datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "common_genes = np.intersect1d(EGAS100010040809_adata.var_names, GSE179994_adata.var_names)\n",
    "print(len(common_genes))\n",
    "\n",
    "# filter out genes that are not in common_genes\n",
    "EGAS100010040809_adata = EGAS100010040809_adata[:, common_genes]\n",
    "GSE179994_adata = GSE179994_adata[:, common_genes]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EGAS100010040809_adata.raw = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GSE179994_adata.raw = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merge adatas from GSE179994_meta and GSE179994_all.Tcell.rawCounts.h5ad "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adata = EGAS100010040809_adata.concatenate(GSE179994_adata)\n",
    "adata = ad.concat([EGAS100010040809_adata, GSE179994_adata])\n",
    "print(adata.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(adata.obs['batch'])\n",
    "adata.obs['batch'].value_counts()\n",
    "#print(EGAS100010040809_adata.obs['batch'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add combined_meta to adata.ob\n",
    "adata.write('/data/ICI_exprs/Merged_GSE179994_EGAS100010040809_T_cell_.h5ad')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Process combined data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.settings.verbosity = 3             # verbosity: errors (0), warnings (1), info (2), hints (3)\n",
    "sc.logging.print_header()\n",
    "sc.settings.set_figure_params(dpi=80, facecolor='white')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Removing cells expressing <500 || >5000 genes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# removing cells expressing <500 || >5000 genes\n",
    "sc.pp.filter_cells(adata, min_genes=500)\n",
    "sc.pp.filter_cells(adata, max_genes=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print ('Dimention of adata: ' + str(adata.shape))\n",
    "print ('columns for adata.obs: ' + str(adata.obs.columns))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Removing cells containing <400 || >25000 UMIs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# removing cells containing <400 || >25000 UMIs\n",
    "sc.pp.filter_cells(adata, min_counts = 400)\n",
    "sc.pp.filter_cells(adata, max_counts = 25000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print ('Dimention of adata: ' + str(adata.shape))\n",
    "print ('columns for adata.obs: ' + str(adata.obs.columns))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprecess with respect to gene (var)\n",
    "#### Removing genes covered by <3 cells per sample\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# removing genes covered by <3 cells\n",
    "sc.pp.filter_genes(adata, min_cells= 10)\n",
    "adata.var_names_make_unique()\n",
    "adata.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove cells with high percentage of mitochondrial genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# label genes as mt\n",
    "adata.var['mt'] = adata.var_names.str.startswith('MT-')  \n",
    "\n",
    "# annotate cells with the percent of genes assigned as mt\n",
    "sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)\n",
    "\n",
    "# Here we keep cells with < 20% mito ratio\n",
    "adata = adata[adata.obs['pct_counts_mt'] < 20, :]\n",
    "adata.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot statistics regarding cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n",
    "             jitter=0.0, multi_panel=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Integrating Samples"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Normalization & Logarithmization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Log normalization scaled up to 10000\n",
    "print('Before normalization, the sum of first row of X: ' + str(adata.X[0,:].sum()))\n",
    "sc.pp.normalize_total(adata, target_sum=1e4)\n",
    "print('After normalization, the sum of first row of X: ' + str(adata.X[0,:].sum()))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Before log, the sum of first row of X: ' + str(adata.X[0,:].sum()))\n",
    "# Logarithmize adata\n",
    "sc.pp.log1p(adata, base=2)\n",
    "print('After log, the sum of first row of X: ' + str(adata.X[0,:].sum()))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Keep high variance genes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.highly_variable_genes(adata, n_top_genes=5000)\n",
    "sc.pl.highly_variable_genes(adata)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Keep track of original adata and update adata.X to  high variance genes only"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.raw = adata\n",
    "adata = adata[:, adata.var.highly_variable]\n",
    "print('adata dimensions of high variance genes: ' + str(adata.shape))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.write('/data/ICI_exprs/ICI_NHDP/Merged_GSE179994_EGAS100010040809_T_cell_5K_hvg.h5ad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = sc.read('/data/ICI_exprs/ICI_NHDP/Merged_GSE179994_EGAS100010040809_T_cell_5K_hvg.h5ad')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Unsupervised cell clustering & identification of major cell types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform PCA   \n",
    "sc.tl.pca(adata, svd_solver='arpack', n_comps=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(adata.obsm['X_pca'].shape)\n",
    "print(adata.varm['PCs'].shape)\n",
    "print(adata.uns['pca']['variance_ratio'].shape)\n",
    "print(adata.obs.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.neighbors(adata, n_neighbors=80, n_pcs=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the Leiden algorithm to find clusters\n",
    "sc.tl.leiden(adata, resolution=0.5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load/find cell cycle markers: T-test/T-cells\n",
    "sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')\n",
    "sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Unsupervised cell clustering & identification of major cell types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform UMAP\n",
    "sc.tl.umap(adata)\n",
    "sc.pl.umap(adata, color='leiden',legend_loc='on data')\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check batch effect  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color='batch')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perform Harmony analysis\n",
    "The data above show that there is a significant batch effects. Apply Harmony analysis to re-project cells into new PCA spact and re-perfrom clustering analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform batch correction using harmony, which works in the PCA space instead of the original gene space\n",
    "from scanpy.external.pp import harmony_integrate\n",
    "sc.external.pp.harmony_integrate(adata, 'batch',  max_iter_harmony=20, random_state=0)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clustering again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cluster cells again after batch correction\n",
    "sc.pp.neighbors(adata, n_neighbors=80, n_pcs=30, use_rep='X_pca_harmony')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.leiden(adata, resolution=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot UMAP\n",
    "sc.tl.umap(adata)\n",
    "sc.pl.umap(adata, color=['leiden', 'batch', 'timepoint'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save a copy with cluster labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.write('/data/ICI_exprs/ICI_NHDP/Merged_GSE179994_EGAS100010040809_T_cell_5K_hvg_with_cls_meta.h5ad')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check the distribution of T markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# color by gene expression of T cell markers\n",
    "sc.pl.umap(adata, color= ['CD3D', 'CD3E'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Label TNK subtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_cell_makers = {\n",
    "    'CD4'\t: ['CD4', 'IL7R'],\n",
    "    'CD8'\t: [ 'CD8A', 'CD8B'],\n",
    "    'Naïve'\t: ['TCF7', 'SELL', 'LEF1', 'CCR7'],\n",
    "    'Exhausted' : ['LAG3', 'TIGIT', 'PDCD1', 'HAVCR2', 'CTLA4'],\n",
    "    'Cytotoxic' : ['IL2', 'GZMA', 'GNLY', 'PRF1', 'GZMB', 'GZMK', 'IFNG', 'NKG7'],\n",
    "    'Treg' : ['IL2RA', 'FOXP3', 'IKZF2', 'IKZF4',  'TNFRSF18'],\n",
    "    # 'Gamma-delta' : ['TRGC1', 'TRGC2', 'TRDC'],\n",
    "    'Th17' : ['IL17A',  'CCR6', 'KLRB1'],  #'IL22',\n",
    "    'MAIT' : ['SLC4A10', 'KLRB1', 'IL7R', 'DPP4'],  \n",
    "    'ILC' :\t['KIT', 'IL1R1'],\n",
    "    'Th1' :\t['STAT4', 'IL12RB2', 'IFNG'],\n",
    "    'Th2' :\t['GATA3', 'STAT6', 'IL4'],\n",
    "    'Tfh'\t: ['MAF', 'CXCL13', 'CXCR5', 'PDCD1'],\n",
    "    'NK' :  ['XCL1', 'FCGR3A', 'KLRD1', 'KLRF1', 'NCAM1'],\n",
    "    'Proliferation' : ['MKI67', 'PCNA', 'STMN1']\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check if the markers are in the var names\n",
    "for cell_type, markers in T_cell_makers.items():\n",
    "    print (cell_type, \":\", markers)\n",
    "    print (\"number of match in var: \", str(sum(adata.raw.var_names.isin(markers))))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.dendrogram(adata, groupby='leiden')\n",
    "sc.pl.dotplot(adata, T_cell_makers, 'leiden', dendrogram=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cell_type, markers in T_cell_makers.items():\n",
    "    print (cell_type, \":\", markers)\n",
    "    sc.pl.umap(adata, color=markers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obsm.keys"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot PD-1 and potential target genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color= ['PDCD1', 'CXCL13', 'HAVCR2','CTLA4', \"PRDM1\"])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The enriched genes for each cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load/find cell cycle markers: T-test/T-cells\n",
    "sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')\n",
    "sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examine the distribution of PDCD1 and potential target genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the umap with the timepoint\n",
    "sc.pl.umap(adata, color=['timepoint', 'PDCD1', 'CXCL13', 'HAVCR2','CTLA4', 'batch'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show the distribution of the timepoint\n",
    "sc.pl.umap(adata, color = ['timepoint', 'batch'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare PDCD1 and CXCL13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract cells wtih PDCD1 > 0.5\n",
    "adata_PDCD1 = adata[adata.raw[:, 'PDCD1'].X > 0.5, :].copy()\n",
    "# plot violin plot with values > 1\n",
    "sc.pl.violin(adata_PDCD1, ['PDCD1', 'CXCL13', 'TIGIT', 'HAVCR2'], split=True, groupby='timepoint', jitter=0.0, multi_panel=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.violin(adata_PDCD1, ['PDCD1', 'CXCL13', 'TIGIT', 'HAVCR2'], split=True, groupby='batch',  jitter=0.0, multi_panel=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It appears that cells expressing these genes exhibit bimodal distribution.  Whether different mode correspond to different cell type is unclear"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examine the relationship between PD-1 and potential target genes pre-treatment\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot scatter plot for pre \n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(12, 4))\n",
    "axes[0].scatter(CD274_pre, cxcl13_pre, s=2, color='black')\n",
    "axes[0].set_xlabel('CD274')\n",
    "axes[0].set_ylabel('CXCL13')\n",
    "axes[1].scatter(pdcd1_pre, cxcl13_pre, s=2, color='black')\n",
    "axes[1].set_xlabel('PDCD1')\n",
    "# plt.title('Pre')\n",
    "axes[2].scatter(pdcd1_pre * CD274_pre, cxcl13_pre, s=2, color='black' )\n",
    "axes[2].set_xlabel('PDCD1 * CD274')\n",
    "fig.tight_layout()\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot relationship on treatment\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract from adata_sample_tpm \n",
    "cxcl13_on = adata_sample_tpm.X[adata_sample_tpm.obs['timepoint'] == 'On', adata_sample_tpm.var_names == 'CXCL13'] \n",
    "pdcd1_on = adata_sample_tpm.X[adata_sample_tpm.obs['timepoint'] == 'On', adata_sample_tpm.var_names == 'PDCD1']\n",
    "CD274_on = adata_sample_tpm.X[adata_sample_tpm.obs['timepoint'] == 'On', adata_sample_tpm.var_names == 'CD274']\n",
    "CTLA4_on = adata_sample_tpm.X[adata_sample_tpm.obs['timepoint'] == 'On', adata_sample_tpm.var_names == 'CTLA4']\n",
    "GZMK_on = adata_sample_tpm.X[adata_sample_tpm.obs['timepoint'] == 'On', adata_sample_tpm.var_names == 'GZMK']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot scatter plot for pre \n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(12, 4))\n",
    "axes[0].scatter(CD274_on, cxcl13_on, s=2, color='black')\n",
    "axes[0].set_xlabel('CD274')\n",
    "axes[0].set_ylabel('CXCL13')\n",
    "axes[1].scatter(pdcd1_on, cxcl13_on, s=2, color='black')\n",
    "axes[1].set_xlabel('PDCD1')\n",
    "axes[2].scatter(pdcd1_on * CD274_on, cxcl13_on, s=2, color='black' )\n",
    "axes[2].set_xlabel('PDCD1 * CD274')\n",
    "fig.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot scatter plot for pre \n",
    "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 4))\n",
    "axes[0].scatter(pdcd1_pre * CD274_pre, cxcl13_pre, s=2, color='black' )\n",
    "axes[0].set_xlabel('PDCD1 * CD274')\n",
    "axes[1].scatter(pdcd1_on * CD274_on, cxcl13_on, s=2, color='black' )\n",
    "axes[1].set_xlabel('PDCD1 * CD274')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  },
  "vscode": {
   "interpreter": {
    "hash": "6099bf950c4e2d0da650f89b132257c4802c810b0835ed278af9d4dcc1b49a4c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}