reproducing_cell_type_motif_enrichment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a876ec18",
   "metadata": {},
   "source": [
    "# Reproducing the cell type motif enrichment\n",
    "\n",
    "This workbook explains the process of reproducing the cell type motif enrichment\n",
    "analysis as seen in our publication. The steps covered are:\n",
    "\n",
    "  * Predicting genome-wide signal for the cell type of interest\n",
    "  * Running motif analysis (calculating graidents with respect to the input)\n",
    "  * Running cell type enrichment analysis on identified transcription factors"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6d6a3084",
   "metadata": {},
   "source": [
    "## 1. Predicting genome-wide signal for the cell type of interest\n",
    "\n",
    "First, predict genome-wide (chromosomes 1-22) for a cell type of interest and save the result as a bigWig file (a separate one for each histone mark). See below for the command to use with the EnformerCelltyping conda env to make a bigwig for each histone mark for the microglia from Nott et al., 2019:\n",
    "\n",
    "```\n",
    "python ./bin/predict_genome.py -c Nott19_Microglia -p data/demo/Nott19_Microglia_128.bigWig -o ./model_results/predictions/\n",
    "```\n",
    "\n",
    "Note this should only really be run with a GPU as it takes 33 hours with a GPU. The script can be adapted to run in parallel. Also make sure to process your ATAC-seq data as described in section 1.1 before running the above on it.\n",
    "\n",
    "Second Note - if you want to predict in more than one cell type/tissue, I advise precomputing ands saving the DNA embedding by passign it through enformer, this will massively speed up genome-wide predictions for your cell type of interest so pays off if you (eventually) want to look at more than 1 cell type. This does require a substantial amount of disk space (~240 GB). To use this approach, first precomute the DNA embeddings with:\n",
    "\n",
    "```\n",
    "python ./bin/precompute_dna_embeddings.py \n",
    "```\n",
    "\n",
    "Then use these embeddings when predicting in your cell type of interest with:\n",
    "\n",
    "```\n",
    "python ./bin/predict_genome_precomp.py -c Nott19_Microglia -p data/demo/Nott19_Microglia_128.bigWig -o ./model_results/predictions/\n",
    "```\n",
    "\n",
    "You can change the DNA Embedding directory with the -d parameter.\n",
    "\n",
    "In our paper we analysed H3K27ac, H3K4me1 and H3K9me3 for seven cell types - Keratinocytes, Epimap's Astrocyte, Nott. et al. 2019 Astrocytes, Monocytes, Nott. et al. 2019 Microglia, Nott. et al. 2019 Neurons and Heart."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd6b0bb4",
   "metadata": {},
   "source": [
    "## 2. Running motif analysis\n",
    "\n",
    "Next get the top 10% of these peaks by calculcating gradients with respect to the input. Note, this takes\n",
    "a while to run and really requires a large RAM GPU (compledted with an 80GB GPU, taking 1-2 days per cell type and mark). Run the following using the `EnformerCelltyping` conda environment:\n",
    "\n",
    "```\n",
    "python ./bin/motif_analysis.py -cell Nott19_Microglia -hist h3k27ac\n",
    "```\n",
    "\n",
    "And similarly for all marks/cell types.\n",
    "\n",
    "\n",
    "Next run the motif analysis using [Homer](http://homer.ucsd.edu/homer/motif/), run:\n",
    "\n",
    "```\n",
    "bash ./bin/motif_analysis.sh Nott19_Microglia\n",
    "```\n",
    "\n",
    "Note that this uses the `bigwigtobed` conda environment which you can get both from the `./environment` folder yml file.\n",
    "\n",
    "To identify the peaks to use for motif analysis, the script ranks peaks by those most affected by the global chromatin accessibility. The DNA motifs in the top decile are passed to Homer to identify transcription factor enrichment for peaks that are affected by the global chromatin accessibility signal for the cell.\n",
    "\n",
    "This will create two files:\n",
    "\n",
    "* `./model_results/motif_analysis/homer_background_tfs.csv` - This is the background list of all TF's tested for with Homer\n",
    "* `./model_results/motif_analysis/all_cells_sig_gbl_motifs.csv` - The TF's Homer found significantly enriched in the motfisd from the histone mark peaks relating to the global cell type chromatin accessibility signal.\n",
    "\n",
    "These two files will then be used to run the cell type enrichment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f0948be",
   "metadata": {},
   "source": [
    "## 3. Running cell type enrichment analysis on identified transcription factors\n",
    "\n",
    "Run [EWCE](https://www.frontiersin.org/articles/10.3389/fnins.2016.00016/full) to look for cell type enrichment of the transcription factors identified from the motif analysis. The R script used to run this and plot the results is here:\n",
    "\n",
    "```\n",
    "Rscript ./bin/cell_type_enrichment.R\n",
    "```\n",
    "\n",
    "This uses the [Descartes](doi:10.1126/science.aba7721) whole body dataset so should give an cell type enrichment applicable to all, high level cell types."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03de285a",
   "metadata": {},
   "source": [
    "We ran this for 7 cell types for the manuscript and the results looked like this:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4d93647",
   "metadata": {},
   "source": [
    "![cell type enrichment analysis](images/tf_cell_specificity.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EnformerCelltyping",
   "language": "python",
   "name": "enformercelltyping"
  },
  "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}