Exercise_Prostate_Cancer.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Exercise_Prostate_Cancer.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E8V64RlaB8Ju"
      },
      "source": [
        "# **Transcriptional mediators of treatment resistance in lethal prostate cancer**\n",
        "*Nature Medicine* volume 27, pages 426–433 (2021) [PMID: 33664492](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7960507/)\n",
        "\n",
        "\n",
        "\n",
        "**Abstract**\n",
        "\n",
        "Metastatic castration-resistant prostate cancer is typically lethal, exhibiting intrinsic or acquired resistance to second-generation androgen-targeting therapies and minimal response to immune checkpoint inhibitors1. Cellular programs driving resistance in both cancer and immune cells remain poorly understood. We present single-cell transcriptomes from 14 patients with advanced prostate cancer, spanning all common metastatic sites. Irrespective of treatment exposure, adenocarcinoma cells pervasively coexpressed multiple androgen receptor isoforms, including truncated isoforms hypothesized to mediate resistance to androgen-targeting therapies2,3. Resistance to enzalutamide was associated with cancer cell–intrinsic epithelial–mesenchymal transition and transforming growth factor-β signaling. Small cell carcinoma cells exhibited divergent expression programs driven by transcriptional regulators promoting lineage plasticity and HOXB5, HOXB6 and NR1D2 (refs. 4,5,6). Additionally, a subset of patients had high expression of dysfunction markers on cytotoxic CD8+ T cells undergoing clonal expansion following enzalutamide treatment. Collectively, the transcriptional characterization of cancer and immune cells from human metastatic castration-resistant prostate cancer provides a basis for the development of therapeutic approaches complementing androgen signaling inhibition.\n",
        "\n",
        "![picture](https://media.springernature.com/full/springer-static/image/art%3A10.1038%2Fs41591-021-01244-6/MediaObjects/41591_2021_1244_Fig4_HTML.png?as=webp)\n",
        "**Fig. 4: Cytotoxic lymphocyte populations and clonal expansions in metastatic niches.**\n",
        "\n",
        "**a**, Subclustering of NK and T cells. Each dot represents a single cell projected onto uniform manifold approximation and projection (UMAP) space colored corresponding to clustering via the Louvain algorithm. Clusters are manually labeled with dominant phenotype/cell type from patterns of marker gene expression. Cluster colors are used throughout subpanels. **b**, Expression of select marker, effector and co-inhibitory receptor genes within cytotoxic clusters, CD16+ NK (n = 30), CD8+GNLY+ (n = 84), CD8+CXCR4+ (n = 157) and CD8+PD-1+ (n = 106). Nominal P values from a two-sided Mann–Whitney U-test are shown. TPM, transcripts per million. **c**, T cell clonotypes from scRNA-seq TCR reconstruction. Each bar represents cells sharing a reconstructed productive TCR CDR3 sequence. Bars are grouped by patient. SCC, small cell carcinoma. **d**, Changes in clonal fractions of cytotoxic T cell clonotypes in patient 01115655 following enzalutamide treatment. Each subplot corresponds to a single clonotype with TCRα and β CDR3 amino acid sequences paired from scRNA-seq. Clonal fractions for the same CDR3 sequences (matching at both nucleotide and amino acid level) inferred from TCR reconstruction in bulk RNA-seq are plotted. All detected single cells of the displayed clonotypes came from the CD8+PD-1+ T cell cluster. **e**, PDCD1 to CD8A expression ratio in bulk RNA-seq of paired biopsies from the same patient before and after enzalutamide treatment. Paired biopsies did not always derive from the same site; color of dots indicate biopsy site types. For patient 01115462, the pre-enzalutamide biopsy was collected from the sacrum and the post-enzalutamide biopsy was collected from a left ischial lesion. **f**, Expanded T cell clonotypes following enzalutamide in patients 01115462 and 01115467. Each dot represents a single CDR3 sequence detected in bulk RNA-seq. CDR3 sequences identified both before and after enzalutamide are connected by a broken gray line.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "z1Nc7IazJDjs"
      },
      "source": [
        "#imports\n",
        "import os\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import seaborn as sns\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.stats import zscore\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ByPIzLHoFQYj"
      },
      "source": [
        "## <strong>Data Import</strong>\n",
        "The supplementary data of this publication is available at the Broad Institute’s Single Cell Portal (“[SCP](https://singlecell.broadinstitute.org/single_cell)”). For this bootcamp, I have uploaded the data in our GitHub repository so that you don't need to create an SCP account. However, I encourage you to explore this resource by your own because it contains lots of interesting data that you can use for your projects and the user interface is very intuitive and allows you to perform some exploratory visualizations.\n",
        "\n",
        "<br>\n",
        "\n",
        "Let's have a look at the entry for this project: [SCP1244](https://singlecell.broadinstitute.org/single_cell/study/SCP1244/transcriptional-mediators-of-treatment-resistance-in-lethal-prostate-cancer)\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "C-_OCw6TFRM4"
      },
      "source": [
        "# import workshop data directly from a URL:\n",
        "df_all = pd.read_csv('https://raw.githubusercontent.com/MMRES-PyBootcamp/MMRES-python-bootcamp2021/master/prostate_cancer_data/scp_clustering.tsv',sep='\\t',index_col=0,skiprows=[1])\n",
        "df_nkt = pd.read_csv('https://raw.githubusercontent.com/MMRES-PyBootcamp/MMRES-python-bootcamp2021/master/prostate_cancer_data/scp_nk_t_clustering.tsv',sep='\\t',index_col=0,skiprows=[1])\n",
        "df_metadata = pd.read_csv('https://raw.githubusercontent.com/MMRES-PyBootcamp/MMRES-python-bootcamp2021/master/prostate_cancer_data/scp_metadata.tsv',sep='\\t',index_col=0,skiprows=[1])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MYwhaTeCNhxu"
      },
      "source": [
        "# download, uncompress and import workshop data\n",
        "os.getcwd()\n",
        "os.mkdir('./assets/')\n",
        "!wget --quiet -O -  https://github.com/MMRES-PyBootcamp/MMRES-python-bootcamp2021/raw/master/prostate_cancer_data/scp_tpm.tsv.gz | gunzip > ./assets/scp_tpm.tsv\n",
        "df_exp=pd.read_csv('./assets/scp_tpm.tsv',sep='\\t',index_col=0).T\n",
        "df_exp.index=df_exp.index.astype(np.int)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IBy-DR97gN2p"
      },
      "source": [
        "## <strong>Whole Dataset</strong>\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-dO6ezE2BzaK"
      },
      "source": [
        "df=df_exp.merge(df_all,right_index=True,left_index=True,how='inner')\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IOleqdMQX3nm"
      },
      "source": [
        "\n",
        "f, ax = plt.subplots(figsize=(6.5, 6.5))\n",
        "sns.boxplotplot(data=df,x='cluster dominant cell type',y='CD3G')\n",
        "ax.set_xticklabels()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "y4gULinqYazL"
      },
      "source": [
        "# let's order the cell types according to the expression of this marker gene\n",
        "s_celltypes=sorted([(col,val['CD3G'].mean()) for col,val in df.groupby('cluster dominant cell type')], key=lambda x: x[1],reverse=True)\n",
        "s_celltypes\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EyS0TSY7W2nu"
      },
      "source": [
        "\n",
        "ordered_cells=np.array(s_celltypes).T[0]\n",
        "\n",
        "f, ax = plt.subplots(figsize=(6.5, 6.5))\n",
        "sns.boxplot(data=df,x='cluster dominant cell type',y='CD3G',order=ordered_cells)\n",
        "ax.set_xticklabels(ordered_cells,rotation=90)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xKrb6oTParnu"
      },
      "source": [
        "f, ax = plt.subplots(figsize=(6.5, 6.5))\n",
        "\n",
        "sns.scatterplot(x=\"X\", y=\"Y\",\n",
        "                hue=\"cluster dominant cell type\",\n",
        "                hue_order=ordered_cells,\n",
        "                linewidth=0,\n",
        "                data=df, ax=ax)\n",
        "\n",
        "ax.legend(bbox_to_anchor=(1.1, 1.05))\n",
        "\n",
        "f, ax = plt.subplots(figsize=(6.5, 6.5))\n",
        "\n",
        "sns.scatterplot(x=\"X\", y=\"Y\",\n",
        "                hue=\"CD3G\",\n",
        "                palette=plt.cm.seismic,hue_norm=(0,500),hue_order=\"CD3G\",\n",
        "                linewidth=0,\n",
        "                data=df, ax=ax)\n",
        "\n",
        "ax.legend(bbox_to_anchor=(1.5, 1.05))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "t0o1mm6UcAKX"
      },
      "source": [
        "sns.relplot(\n",
        "    data=df, x=\"X\", y=\"Y\",\n",
        "    col=\"cluster dominant cell type\", hue=\"CD3G\", \n",
        "    col_wrap=4, palette=plt.cm.seismic,hue_norm=(0,500), col_order=ordered_cells,\n",
        "    kind=\"scatter\"\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aCKfoDKFgVwU"
      },
      "source": [
        "## <strong>Cytotoxic lymphocyte populations</strong>\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6KPAkOu4hF0T"
      },
      "source": [
        "df=df_exp.merge(df_nkt,right_index=True,left_index=True,how='inner')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "K-Bw58OsjWPv"
      },
      "source": [
        "color_dict={\n",
        "    'CD8+ PD1+ T cell':'#2a9a79', \n",
        "    'CD8+ CXCR4+ T cell':'#7eb459', \n",
        "    'CD8+ GNLY+ T cell':'#486aa5',\n",
        "    'CD16+ NK cell':'#7e3689', \n",
        "    'CD4+ Treg cell':'#ec9034', \n",
        "    'CD4+ T cell':'#fd7971'\n",
        "}"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "13x0h82OhMaF"
      },
      "source": [
        "# let's order the cell types according to the expression of CD3 marker gene\n",
        "s_celltypes=sorted([(col,val['CD3G'].mean()) for col,val in df.groupby('cluster dominant cell type')], key=lambda x: x[1],reverse=True)\n",
        "ordered_cells=np.array(s_celltypes).T[0]\n",
        "f, ax = plt.subplots(figsize=(6.5, 6.5))\n",
        "sns.boxplot(data=df,x='cluster dominant cell type',y='CD3G',order=ordered_cells, showfliers = False,palette=color_dict)\n",
        "ax.set_xticklabels(ordered_cells,rotation=90)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mDYncWwHjG_m"
      },
      "source": [
        "# let's order the cell types according to the expression of CD4 marker gene\n",
        "s_celltypes=sorted([(col,val['CD4'].mean()) for col,val in df.groupby('cluster dominant cell type')], key=lambda x: x[1],reverse=True)\n",
        "ordered_cells=np.array(s_celltypes).T[0]\n",
        "f, ax = plt.subplots(figsize=(6.5, 6.5))\n",
        "sns.boxplot(data=df,x='cluster dominant cell type',y='CD4',order=ordered_cells, showfliers = False,palette=color_dict)\n",
        "ax.set_xticklabels(ordered_cells,rotation=90)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Rfost3ygjJ3g"
      },
      "source": [
        "# let's order the cell types according to the expression of CD16 (FCGR3B) marker gene\n",
        "s_celltypes=sorted([(col,val['FCGR3B'].mean()) for col,val in df.groupby('cluster dominant cell type')], key=lambda x: x[1],reverse=True)\n",
        "ordered_cells=np.array(s_celltypes).T[0]\n",
        "f, ax = plt.subplots(figsize=(6.5, 6.5))\n",
        "sns.boxplot(data=df,x='cluster dominant cell type',y='FCGR3B',order=ordered_cells, showfliers = False,palette=color_dict)\n",
        "ax.set_xticklabels(ordered_cells,rotation=90)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rT4sGeRPi9ud"
      },
      "source": [
        "# let's order the cell types according to the expression of CD8 marker gene\n",
        "s_celltypes=sorted([(col,val['CD8A'].mean()) for col,val in df.groupby('cluster dominant cell type')], key=lambda x: x[1],reverse=True)\n",
        "ordered_cells=np.array(s_celltypes).T[0]\n",
        "f, ax = plt.subplots(figsize=(6.5, 6.5))\n",
        "sns.boxplot(data=df,x='cluster dominant cell type',y='CD8A',order=ordered_cells, showfliers = False,palette=color_dict)\n",
        "ax.set_xticklabels(ordered_cells,rotation=90)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dE8LnB-slZbV"
      },
      "source": [
        "ordered_cells=['CD16+ NK cell', 'CD8+ GNLY+ T cell', 'CD8+ CXCR4+ T cell','CD8+ PD1+ T cell']\n",
        "f, ax = plt.subplots(figsize=(3, 6.5))\n",
        "sns.stripplot(data=df,x='cluster dominant cell type',y='PDCD1',order=ordered_cells, palette=color_dict)\n",
        "ax.set_xticklabels(ordered_cells,rotation=90)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JEVTYpvxmXOM"
      },
      "source": [
        "ordered_cells=['CD16+ NK cell', 'CD8+ GNLY+ T cell', 'CD8+ CXCR4+ T cell','CD8+ PD1+ T cell']\n",
        "f, ax = plt.subplots(figsize=(3, 6.5))\n",
        "sns.stripplot(data=df,x='cluster dominant cell type',y='CX3CR1',order=ordered_cells, palette=color_dict)\n",
        "ax.set_xticklabels(ordered_cells,rotation=90)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ytk9ObpDp-op"
      },
      "source": [
        "## <strong>Heatmaps</strong>\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "W-dsRVLIqIcK"
      },
      "source": [
        "markers=['CD4','CD8A','CD3G','FCGR3B','PDCD1','HAVCR2','TOX','CX3CR1','CXCR4','GZMB']"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "n3cwfjKwmyZo"
      },
      "source": [
        "sns.heatmap(df[markers],cmap=plt.cm.seismic,vmin=0,vmax=500)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KhoLLAmQqfBj"
      },
      "source": [
        "clusters=df['cluster dominant cell type']\n",
        "row_colors = clusters.map(color_dict)\n",
        "sns.clustermap(df[markers].apply(lambda x:zscore(x),axis=0),cmap=plt.cm.seismic,vmin=-2,vmax=2,row_colors=row_colors.to_numpy())\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "22bV3HogpuZ1"
      },
      "source": [
        "sns.clustermap(df[markers],z_score=1,cmap=plt.cm.seismic,vmin=-2,vmax=2,row_colors=row_colors.to_numpy())\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nLxcyoXlrUWx"
      },
      "source": [
        "markers=['CD4','CD8A','FCGR3B','PDCD1','CXCR4','GNLY']\n",
        "#markers=['CD4','CD8A','FCGR3B']\n",
        "sns.clustermap(df[markers],z_score=1,cmap=plt.cm.seismic,vmin=-2,vmax=2,row_colors=row_colors.to_numpy())\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dOqhqZlUvMf6"
      },
      "source": [
        "## <strong> Statistics </strong>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sS-2KpMGu7_H"
      },
      "source": [
        "https://medium.com/insights-school/learn-basic-statistics-with-python-cc0f45275929\n"
      ]
    }
  ]
}