reproducing_ldsc_enrichment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "16750d51",
   "metadata": {},
   "source": [
    "# Reproducing the LDSC enricment analysis\n",
    "\n",
    "This workbook explains the process of reproducing the LDSC enrichment analysis as seen in our publication. The steps covered are:\n",
    "\n",
    "* Predicting genome-wide signal for the cell type of interest\n",
    "* Converting format of predicted signal for LDSC\n",
    "* Running LDSC enrichment analysis for a GWAS sumstats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91298a27",
   "metadata": {},
   "source": [
    "## 1. Predicting genome-wide signal for the cell type of interest\n",
    "\n",
    "First, predict genome-wide (chromosomes 1-22) for the 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 for three of Nott et al.'s cell types - Microglia, Oligodendrocytes and Neurons."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3008e597",
   "metadata": {},
   "source": [
    "## 2. Converting format of predicted signal for LDSC\n",
    "\n",
    "Second, we need to convert the bigwig containing the predicted h3k27ac signal to a format usable by LDSC - bed file.\n",
    "\n",
    "The first step is to make a `bedGraph` files of the `bigwig`.  - Use the `bigwigtobed` conda environment and run (update the path and histone marks/cell types used - currently set up for reproducing the results in the manuscript):\n",
    "\n",
    "```\n",
    "bash bin/mke_bed.sh\n",
    "```\n",
    "\n",
    "Next step is to make `bed-like` files from `bedGraphs`. Note, that this script uses a p-value cut-off of arc-sinh(1) to mark peaks for the bed-like file as we found this led to better enrichment than more stringent cut-offs. Use any environment with `python`, `pandas` and `numpy` installed to run this:\n",
    "\n",
    "```\n",
    "python bin/mke_bed_bedgraph.py\n",
    "```\n",
    "\n",
    "Now you have a bed file which LDSC can work with!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59e3b1a7",
   "metadata": {},
   "source": [
    "## 3. Running LDSC enrichment analysis for a GWAS sumstats\n",
    "\n",
    "The first step of running [LDSC](https://github.com/bulik/ldsc) is to create annotation files for the genome-wide signals (`bed` files) you want to check for GWAS enrichment in.\n",
    "\n",
    "Note that this step and all LDSC steps require reference files for LDSC to be downloaded and directory pointers updated in the scripts to get it to run. Please read each script carefully. Also note the code expects you to have cloned the LDSC repo to the working directory and created the conda environment from it.\n",
    "\n",
    "Use the `ldsc` conda environment (for all of these steps) and run:\n",
    "\n",
    "```\n",
    "bash bin/make_annot.sh\n",
    "```\n",
    "\n",
    "Next we need to generate the LDSC scores based on the created annotations. To do this, run:\n",
    "\n",
    "```\n",
    "bash bin/gen_ldsc_scores.sh\n",
    "```\n",
    "\n",
    "Then the final step is to looks for enrichment based on these generated LDSC scores for the GWAS of interest. Make sure to format your GWAS sumstats suitably for LDSC first, we did this with [MungeSumstats](https://github.com/neurogenomics/MungeSumstats) in R, with the following command:\n",
    "\n",
    "```\n",
    "gwas_paths <- MungeSumstats::format_sumstats(\n",
    "    ids='pth/to/gwas/sumstats/file.tsv.gz',\n",
    "    save_dir = 'save/path/processed_file.tsv',\n",
    "    convert_ref_genome = 'GRCh37',\n",
    "    compute_z = 'BETA',\n",
    "    save_format='LDSC')\n",
    "```\n",
    "\n",
    "Then just rename the tsv file extension to `.sumstats` to work with LDSC. \n",
    "\n",
    "Now you are finally ready to run the LDSC enrichment to get enrichment of the sumstats based on the annotations:\n",
    "\n",
    "```\n",
    "bash bin/run_ldsc_disease_enrich.sh \n",
    "```\n",
    "\n",
    "Note that this will return the enrichment of each baseline annotation and the other annotations you enter in the `ref-ld-chr` field. Also note these won't be labelled as you enter their names but instead just be at the end of the resulting tsv file as `L2_1`, `L2_2`, etc."
   ]
  }
 ],
 "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
}